Update syntax

README.md

@@ -2,62 +2,53 @@
 [![Coverage Status](https://coveralls.io/repos/matthieugomez/LeastSquaresOptim.jl/badge.svg?branch=master&service=github)](https://coveralls.io/github/matthieugomez/LeastSquaresOptim.jl?branch=master)
 ## Motivation
 
-This package solves non linear least squares optimization problems. The package is inspired by the [Ceres library](http://ceres-solver.org/nnls_solving.html). 
+This package solves non linear least squares optimization problems. The package is inspired by the [Ceres library](http://ceres-solver.org/nnls_solving.html).  This package is written with large scale problems in mind (in particular for sparse Jacobians). 
 
 
 ## Simple Syntax
-
 The symple syntax mirrors the `Optim.jl` syntax
 
 ```julia
 using LeastSquaresOptim
 function rosenbrock(x)
 	[1 - x[1], 100 * (x[2]-x[1]^2)]
 end
-optimize(rosenbrock, zeros(2), Dogleg())
-optimize(rosenbrock, zeros(2), LevenbergMarquardt())
+x0 = zeros(2)
+optimize(rosenbrock, x0, Dogleg())
+optimize(rosenbrock, x0, LevenbergMarquardt())
 ```
-You can also add the options : `x_tol`, `f_tol`, `g_tol`, `iterations` and `Δ` (initial radius).
-
-The gradient of `f` is computed using automatic differenciation, through the package `ForwardDiff.jl`. This means that the function `f` should work when `x` is a `DualNumber`.
+You can also add the options : `x_tol`, `f_tol`, `g_tol`, `iterations`, `Δ` (initial radius), and `autodiff`.
 
 
-## Non Allocating Syntax
-This package is written with large scale problems in mind. In particular, memory is allocated once and for all at the start of the function call ; objects are updated in place at each method iteration.  The advanced syntax allows to take full advantage of this.
+## In-place Syntax
 
-1. To find `x` that minimizes `f'(x)f(x)`, construct a `LeastSquaresProblem` object with:
+The alternative syntax allows to use an in place functions or to specify the jacobian. Construct a `LeastSquaresProblem` object with:
  - `x` an initial set of parameters.
  - `f!(out, x)` that writes `f(x)` in `out`.
  - the option `output_length` to specify the length of the output vector. 
- - `g!` a function such that `g!(out, x)` writes the jacobian at x in `out`. Otherwise, the jacobian will be computed with the `ForwardDiff.jl` package
- - `y` a preallocation for `f`
- - `J` a preallocation for the jacobian
-
-
-	A simple example:
-	```julia
-	using LeastSquaresOptim
-	function rosenbrock_f!(out, x)
-	    out[1] = 1 - x[1]
-	    out[2] = 100 * (x[2]-x[1]^2)
-	end
-	optimize!(LeastSquaresProblem(x = zeros(2), f! = rosenbrock_f!, output_length = 2))
-
-	# if you want to use gradient
-	function rosenbrock_g!(J, x)
-	    J[1, 1] = -1
-	    J[1, 2] = 0
-	    J[2, 1] = -200 * x[1]
-	    J[2, 2] = 100
-	end
-	optimize!(LeastSquaresProblem(x = zeros(2), f! = rosenbrock_f!, g! = rosenbrock_g!, output_length = 2))
-	```
-
-2. You can also specify a particular least square solver (a least square optimization method proceeds by solving successively linear least squares problems `min||Ax - b||^2`). 
-	```julia
-	optimize!(LeastSquaresProblem(x = x, f! = rosenbrock_f!, output_length = 2), LeastSquaresOptim.LSMR())
-	```
+ - Optionally, `g!` a function such that `g!(out, x)` writes the jacobian at x in `out`. Otherwise, the jacobian will be computed following the `:autodiff` argument.
+
+```julia
+using LeastSquaresOptim
+function rosenbrock_f!(out, x)
+ out[1] = 1 - x[1]
+ out[2] = 100 * (x[2]-x[1]^2)
+end
+optimize!(LeastSquaresProblem(x = zeros(2), f! = rosenbrock_f!, output_length = 2, autodiff = :central), Dogleg())
+
+# if you want to use gradient
+function rosenbrock_g!(J, x)
+    J[1, 1] = -1
+    J[1, 2] = 0
+    J[2, 1] = -200 * x[1]
+    J[2, 2] = 100
+end
+optimize!(LeastSquaresProblem(x = zeros(2), f! = rosenbrock_f!, g! = rosenbrock_g!, output_length = 2), Dogleg())
+```
 
+## Choice of Optimizer / Least Square Solver
+- You can specify two least squares optimizers, `Dogleg()` and `LevenbergMarquardt()`
+- You can specify three least squares solvers that will be used within the optimizer
 	- `LeastSquaresOptim.QR()`. Available for dense jacobians
 	- `LeastSquaresOptim.Cholesky()`. Available for dense jacobians
 	- `LeastSquaresOptim.LSMR()`. A conjugate gradient method ([LSMR]([http://web.stanford.edu/group/SOL/software/lsmr/) with diagonal preconditioner). The jacobian can be of any type that defines the following interface is defined:
@@ -66,27 +57,22 @@ This package is written with large scale problems in mind. In particular, memory
 		- `colsumabs2!(x, A)` updates x to the sum of squared elements of each column
 		- `size(A, d)` returns the nominal dimensions along the dth axis in the equivalent matrix representation of A.
 		- `eltype(A)` returns the element type implicit in the equivalent matrix representation of A.
-
-		Similarly, `x` or `f(x)` may be custom types. An example of the interface to define can be found in the package [SparseFactorModels.jl](https://github.com/matthieugomez/SparseFactorModels.jl).
-
+	Similarly, `x` or `f(x)` may be custom types. An example of the interface to define can be found in the package [SparseFactorModels.jl](https://github.com/matthieugomez/SparseFactorModels.jl).
 	For the `LSMR` solver, you can optionally specifying a function `preconditioner!` and a matrix `P` such that `preconditioner(x, J, P)` updates `P` as a preconditioner for `J'J` in the case of a Dogleg optimization method, and such that `preconditioner(x, J, λ, P)` updates `P` as a preconditioner for `J'J + λ` in the case of LevenbergMarquardt optimization method. By default, the preconditioner is chosen as the diagonal of of the matrix `J'J`. The preconditioner can be any type that supports `A_ldiv_B!(x, P, y)`
 
-	The `optimizers` and `solvers` are presented in more depth in the [Ceres documentation](http://ceres-solver.org/nnls_solving.html). For dense jacobians, the default options are `Dogle()` and `QR()`. For sparse jacobians, the default options are  `LevenbergMarquardt()` and `LSMR()`. 
+The `optimizers` and `solvers` are presented in more depth in the [Ceres documentation](http://ceres-solver.org/nnls_solving.html). For dense jacobians, the default options are `Dogle()` and `QR()`. For sparse jacobians, the default options are  `LevenbergMarquardt()` and `LSMR()`. 
 
-3. You can even avoid initial allocations by directly passing a `LeastSquaresProblemAllocated` to the `optimize!` function. Such an object bundles a `LeastSquaresProblem` object with a few storage objects. This allows to save memory when repeatedly solving non linear least square problems.
-	```julia
-	rosenbrock = LeastSquaresProblemAllocated(x, fcur, rosenbrock_f!, J, rosenbrock_g!; 
-	                                          LeastSquaresOptim.Dogleg(), LeastSquaresOptim.QR())
-	optimize!(rosenbrock)
-	```
+```julia
+optimize!(LeastSquaresProblem(x = x, f! = rosenbrock_f!, output_length = 2), Dogleg())
+optimize!(LeastSquaresProblem(x = x, f! = rosenbrock_f!, output_length = 2), Dogleg(LeastSquaresOptim.QR()))
+```
 
 
 
 ## Related packages
 Related:
 - [LSqfit.jl](https://github.com/JuliaOpt/LsqFit.jl) is a higher level package to fit curves (i.e. models of the form y = f(x, β))
 - [Optim.jl](https://github.com/JuliaOpt/Optim.jl) solves general optimization problems.
-- [IterativeSolvers.jl](https://github.com/JuliaLang/IterativeSolvers.jl) includes several iterative solvers for linear least squares.
 - [NLSolve.jl](https://github.com/EconForge/NLsolve.jl) solves non linear equations by least squares minimization.
 
 

src/solver/dense_cholesky.jl

@@ -16,7 +16,7 @@ function DenseCholeskyAllocatedSolver(cholm::Tc) where {Tc <: StridedMatrix}
 end
 
 
-function AbstractAllocatedSolver(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer, ::Cholesky) where {Tx, Ty, Tf, TJ <: StridedVecOrMat, Tg}
+function AbstractAllocatedSolver(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::AbstractOptimizer{Cholesky}) where {Tx, Ty, Tf, TJ <: StridedVecOrMat, Tg}
     return DenseCholeskyAllocatedSolver(Array{eltype(nls.J)}(undef, length(nls.x), length(nls.x)))
 end
 

src/solver/dense_qr.jl

@@ -22,7 +22,7 @@ end
 ##
 ##############################################################################
 
-function AbstractAllocatedSolver(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::Dogleg, solver::QR) where {Tx, Ty, Tf, TJ <: StridedVecOrMat, Tg}
+function AbstractAllocatedSolver(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::Dogleg{QR}) where {Tx, Ty, Tf, TJ <: StridedVecOrMat, Tg}
     return DenseQRAllocatedSolver(similar(nls.J), _zeros(nls.y))
 end
 
@@ -43,7 +43,7 @@ end
 ##
 ##############################################################################
 
-function AbstractAllocatedSolver(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::LevenbergMarquardt, solver::QR) where {Tx, Ty, Tf, TJ <: StridedVecOrMat, Tg}
+function AbstractAllocatedSolver(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::LevenbergMarquardt{QR}) where {Tx, Ty, Tf, TJ <: StridedVecOrMat, Tg}
     qrm = zeros(eltype(nls.J), length(nls.y) + length(nls.x), length(nls.x))
     u = zeros(length(nls.y) + length(nls.x))
     return DenseQRAllocatedSolver(qrm, u)

src/solver/iterative_lsmr.jl

@@ -125,7 +125,7 @@ end
 ## eltype, size, mul!
 ##
 ##############################################################################
-function getpreconditioner(nls::LeastSquaresProblem, optimizer::Dogleg, solver::LSMR{Nothing, Nothing})
+function getpreconditioner(nls::LeastSquaresProblem, optimizer::Dogleg{LSMR{Nothing, Nothing}})
     preconditioner! = function(x, J, out)
         colsumabs2!(out._, J)
         Tout = eltype(out._)
@@ -135,8 +135,8 @@ function getpreconditioner(nls::LeastSquaresProblem, optimizer::Dogleg, solver::
     preconditioner = InverseDiagonal(_zeros(nls.x))
     return preconditioner!, preconditioner
 end
-function getpreconditioner(nls::LeastSquaresProblem, optimizer::Dogleg, solver::LSMR)
-    return solver.preconditioner!, solver.preconditioner
+function getpreconditioner(nls::LeastSquaresProblem, optimizer::Dogleg{LSMR{T1, T2}}) where {T1, T2}
+    return optimizer.solver.preconditioner!, optimizer.solver.preconditioner
 end
 
 struct LSMRAllocatedSolver{Tx0, Tx1, Tx2, Tx22, Tx3, Tx4, Tx5, Ty} <: AbstractAllocatedSolver
@@ -151,8 +151,8 @@ struct LSMRAllocatedSolver{Tx0, Tx1, Tx2, Tx22, Tx3, Tx4, Tx5, Ty} <: AbstractAl
 end
 
 
-function AbstractAllocatedSolver(nls::LeastSquaresProblem, optimizer::Dogleg, solver::LSMR)
-    preconditioner!, preconditioner = getpreconditioner(nls, optimizer, solver)
+function AbstractAllocatedSolver(nls::LeastSquaresProblem, optimizer::Dogleg{LSMR{T1, T2}}) where {T1, T2}
+    preconditioner!, preconditioner = getpreconditioner(nls, optimizer)
     LSMRAllocatedSolver(preconditioner!, preconditioner, _zeros(nls.x), _zeros(nls.x), 
         _zeros(nls.x), _zeros(nls.x), _zeros(nls.x), _zeros(nls.y))
 end
@@ -193,7 +193,7 @@ end
 
 ##
 ##############################################################################
-function getpreconditioner(nls::LeastSquaresProblem, optimizer::LevenbergMarquardt, ::LSMR{Nothing, Nothing})
+function getpreconditioner(nls::LeastSquaresProblem, optimizer::LevenbergMarquardt{LSMR{Nothing, Nothing}})
     preconditioner! = function(x, J, damp, out)
         colsumabs2!(out._, J)
         Tout = eltype(out._)
@@ -205,10 +205,11 @@ function getpreconditioner(nls::LeastSquaresProblem, optimizer::LevenbergMarquar
     return preconditioner!, preconditioner
 end
 
-function getpreconditioner(nls::LeastSquaresProblem, optimizer::LevenbergMarquardt, solver::LSMR)
-    return solver.preconditioner!, solver.preconditioner
+function getpreconditioner(nls::LeastSquaresProblem, optimizer::LevenbergMarquardt{LSMR{T1, T2}}) where {T1, T2}
+    return optimizer.solver.preconditioner!, optimizer.solver.preconditioner
 end
 
+
 struct LSMRDampenedAllocatedSolver{Tx0, Tx1, Tx2, Tx22, Tx3, Tx4, Tx5, Tx6, Ty} <: AbstractAllocatedSolver
     preconditioner!::Tx0
     preconditioner::Tx1
@@ -221,8 +222,8 @@ struct LSMRDampenedAllocatedSolver{Tx0, Tx1, Tx2, Tx22, Tx3, Tx4, Tx5, Tx6, Ty}
     u::Ty
 end
 
-function AbstractAllocatedSolver(nls::LeastSquaresProblem, optimizer::LevenbergMarquardt, solver::LSMR)
-    preconditioner!, preconditioner = getpreconditioner(nls, optimizer, solver)
+function AbstractAllocatedSolver(nls::LeastSquaresProblem,  optimizer::LevenbergMarquardt{LSMR{T1, T2}}) where {T1, T2}
+    preconditioner!, preconditioner = getpreconditioner(nls, optimizer)
     LSMRDampenedAllocatedSolver(preconditioner!, preconditioner, _zeros(nls.x), _zeros(nls.x), _zeros(nls.x), _zeros(nls.x), _zeros(nls.x),  _zeros(nls.x), _zeros(nls.y))
 end
 

src/types.jl

@@ -1,6 +1,6 @@
- ##############################################################################
+##############################################################################
 ##
-## Non Linear Least Squares
+## Non Linear Least Squares Problem
 ##
 ##############################################################################
 
@@ -33,7 +33,6 @@ function LeastSquaresProblem(x::Tx, y::Ty, f!::Tf, J::TJ, g!::Tg) where {Tx, Ty,
 end
 
 
-
 function LeastSquaresProblem(;x = error("initial x required"), y = nothing, f! = error("initial f! required"), g! = nothing, J = nothing, output_length = 0, autodiff = :central)
     if typeof(y) == Nothing
         if output_length == 0
@@ -64,18 +63,11 @@ function LeastSquaresProblem(;x = error("initial x required"), y = nothing, f! =
 end
 
 
-
 ###############################################################################
 ##
-## Non Linear Least Squares Allocated
-## groups a LeastSquaresProblem with allocations
+## Optimizer and SOlver
 ##
 ##############################################################################
-# optimizer
-abstract type AbstractOptimizer end
-struct Dogleg <: AbstractOptimizer end
-struct LevenbergMarquardt <: AbstractOptimizer end
-
 
 
 # solver
@@ -88,6 +80,24 @@ struct LSMR{T1, T2} <: AbstractSolver
 end
 LSMR() = LSMR(nothing, nothing)
 
+
+# optimizer
+abstract type AbstractOptimizer{T} end
+struct Dogleg{T} <: AbstractOptimizer{T}
+    solver::T
+end
+Dogleg() = Dogleg(nothing)
+struct LevenbergMarquardt{T} <: AbstractOptimizer{T}
+    solver::T
+end
+LevenbergMarquardt() = LevenbergMarquardt(nothing)
+
+
+_solver(x::AbstractOptimizer) = x.solver
+_solver(x::Nothing) = nothing
+
+
+
 ## for dense matrices, default to cholesky ; otherwise LSMR
 function default_solver(x::AbstractSolver, J)
     if (typeof(x) <: QR) && (typeof(J) <: SparseMatrixCSC)
@@ -99,11 +109,19 @@ default_solver(::Nothing, J::StridedVecOrMat) = QR()
 default_solver(::Nothing, J) = LSMR()
 
 ## for LSMR, default to levenberg_marquardt ; otherwise dogleg
-default_optimizer(x::AbstractOptimizer, y) = x
-default_optimizer(::Nothing, ::LSMR) = LevenbergMarquardt()
-default_optimizer(::Nothing, ::AbstractSolver) = Dogleg()
+default_optimizer(x::Dogleg, y::AbstractSolver) = Dogleg(y)
+default_optimizer(x::LevenbergMarquardt, y::AbstractSolver) = LevenbergMarquardt(y)
+default_optimizer(::Nothing, ::LSMR) = LevenbergMarquardt(LSMR())
+default_optimizer(::Nothing, x) = Dogleg(x)
+
+
 
 
+##############################################################################
+##
+## Non Linear Least Squares Problem Allocated
+##
+##############################################################################
 
 abstract type AbstractAllocatedOptimizer end
 abstract type AbstractAllocatedSolver end
@@ -119,22 +137,15 @@ mutable struct LeastSquaresProblemAllocated{Tx, Ty, Tf, TJ, Tg, Toptimizer <: Ab
 end
 
 # Constructor
-function LeastSquaresProblemAllocated(nls::LeastSquaresProblem, optimizer::Union{Nothing, AbstractOptimizer}, solver::Union{Nothing, AbstractSolver})
-    solver = default_solver(solver, nls.J)
+function LeastSquaresProblemAllocated(nls::LeastSquaresProblem, optimizer::Union{Nothing, AbstractOptimizer})
+    solver = default_solver(_solver(optimizer), nls.J)
     optimizer = default_optimizer(optimizer, solver)
     LeastSquaresProblemAllocated(
-        nls.x, nls.y, nls.f!, nls.J, nls.g!, AbstractAllocatedOptimizer(nls, optimizer), AbstractAllocatedSolver(nls, optimizer, solver))
+        nls.x, nls.y, nls.f!, nls.J, nls.g!, AbstractAllocatedOptimizer(nls, optimizer), AbstractAllocatedSolver(nls, optimizer))
 end
 function LeastSquaresProblemAllocated(args...; kwargs...)
     LeastSquaresProblemAllocated(LeastSquaresProblem(args...); kwargs...)
 end
-
-# optimize
-function optimize!(nls::LeastSquaresProblem, optimizer::Union{Nothing, AbstractOptimizer} = nothing, solver::Union{Nothing, AbstractSolver} = nothing; kwargs...)
-    nlsp = LeastSquaresProblemAllocated(nls, optimizer, solver)
-    optimize!(nlsp; kwargs...)
-end
-
 ###############################################################################
 ##
 ## Optim-like syntax
@@ -145,6 +156,20 @@ function optimize(f, x, t::AbstractOptimizer; autodiff = :central, kwargs...)
     optimize!(LeastSquaresProblem(x = deepcopy(x), f! = (out, x) -> copyto!(out, f(x)), output_length = length(f(x)), autodiff = autodiff), t; kwargs...)
 end
 
+function optimize!(nls::LeastSquaresProblem, optimizer::Union{Nothing, AbstractOptimizer} = nothing; kwargs...)
+    optimize!(LeastSquaresProblemAllocated(nls, optimizer); kwargs...)
+end
+
+
+
+
+Base.@deprecate optimize!(nls::LeastSquaresProblem, optimizer::AbstractOptimizer, solver::Union{Nothing, AbstractSolver}; kwargs...) optimize!(nls, typeof{optimizer}(solver); kwargs...)
+Base.@deprecate optimize!(nls::LeastSquaresProblem, solver::AbstractSolver; kwargs...) optimize!(nls, Dogleg{solver}; kwargs...)
+Base.@deprecate LeastSquaresProblemAllocated(nls::LeastSquaresProblem, optimizer::AbstractOptimizer, solver::AbstractSolver) LeastSquaresProblemAllocated(nls, default_optimizer(optimizer, solver))
+
+
+
+
 ###############################################################################
 ##
 ## Result of Non Linear Least Squares

test/nonlinearfitting.jl

@@ -1454,15 +1454,15 @@ end
 
 tests =  [misra1a, Chwirut2, Chwitrut1, Lanczos3, Gauss1, Gauss2, DanWood, Misra1b, MGH09, Thurber, BoxBOD, Rat42, MGH10, Eckerle4, Rat43, Bennett5]
 
-for (optimizer, optimizer_abbr) in ((LeastSquaresOptim.Dogleg(), :dl), (LeastSquaresOptim.LevenbergMarquardt(), :lm))
+for (optimizer, optimizer_abbr) in ((LeastSquaresOptim.Dogleg, :dl), (LeastSquaresOptim.LevenbergMarquardt, :lm))
     n = 0
     N = 0
     for test in tests
         global name, data, parameters, f, solution = test()
       #  @show solution
         for j in 1:size(parameters, 2)
             global nls = LeastSquaresProblem(x = parameters[:, j], f! = (fcur, x) -> ff!(fcur, x, f, data), output_length = size(data, 1))
-            global r = optimize!(nls, optimizer, LeastSquaresOptim.QR(), x_tol = 1e-50, f_tol = 1e-36, g_tol = 1e-50)
+            global r = optimize!(nls, optimizer(LeastSquaresOptim.QR()), x_tol = 1e-50, f_tol = 1e-36, g_tol = 1e-50)
             n += norm(r.minimizer - solution) <= 1e-3
             N += 1
             @test !isnan(mean(r.minimizer) )

test/nonlinearleastsquares.jl

@@ -88,7 +88,7 @@ function factor()
     return name, f!, g!, x
 end
 iter = 0
-for (optimizer, optimizer_abbr) in ((LeastSquaresOptim.Dogleg(), :dl), (LeastSquaresOptim.LevenbergMarquardt(), :lm))
+for (optimizer, optimizer_abbr) in ((LeastSquaresOptim.Dogleg, :dl), (LeastSquaresOptim.LevenbergMarquardt, :lm))
     for (solver, solver_abbr) in ((LeastSquaresOptim.QR(), :qr), (LeastSquaresOptim.LSMR(), :iter))
         global iter += 1
         global name, f!, g!, x = factor()
@@ -99,7 +99,7 @@ for (optimizer, optimizer_abbr) in ((LeastSquaresOptim.Dogleg(), :dl), (LeastSqu
             J = sparse(J)
         end
         global nls = LeastSquaresProblem(x = x, y = fcur, f! = f!, J = J, g! = g!)
-        global r = optimize!(nls, optimizer, solver)
+        global r = optimize!(nls, optimizer(solver))
         if iter == 1
             show(r)
         end