Bound constrained nonlinear least squares

In this part we see how to solve problems of the form:

\[\begin{align*} \min\limits_{\theta} & \frac{1}{2}\|r(\theta)\|^2 \\ & \theta_l\le\theta\le\theta_u \end{align*}\]

Compared to the unconstrained case, there are essentially two differences:

  • the solve() function has an extra parameter, bc that store bound constraints.
  • The solver category is different and Abstract_Solver_Conf, it replaced by Abstract_BC_Solver_Conf, where BC stands for bound constrained.

Furthers details can be found there: solve(nls::AbstractNLS, θ_init::AbstractVector, bc::BoundConstraints, conf::Abstract_BC_Solver_Conf)

You can reproduce the results below using sandbox/example_Rosenbrock.jl

Problem definition

Identical to the unconstrained case. We have:

nls = create_NLS_problem_using_ForwardDiff(2 => 2) do θ
  sqrt(2)* [ 1-θ[1], 10*(θ[2]-θ[1]^2) ]
end

Choose a solver

Algorithm parameters are defined by sub-typing Abstract_BC_Solver_Conf. This structure is then used to identify the selected algorithm. For the moment there is only one implementation, the a bound constrained version of the classical Levenberg-Marquardt method:

conf = LevenbergMarquardt_BC_Conf()
LevenbergMarquardt_BC_Conf(LevenbergMarquardt_Conf(1000, 1.0e-8, 1.0e-8, 0.001), 10, NLS_Solver.Kunisch_Rendl_Conf(50))

We also need a starting point for the $\theta$. Here we start at point $(3,3)$:

θ_init = Float64[3,3]
2-element Vector{Float64}:
 3.0
 3.0

The solve() function

To solve the problem you simply have to call the solve() function. For bound constrained problems, this function has the following prototype

function solve(nls::AbstractNLS,
               θ_init::AbstractVector,
               bc::BoundConstraints,
               conf::Abstract_Solver_Conf)::Abstract_Solver_Result

In our case, if we want $2 \le \theta_i \le 4$ this gives

θl = Float64[2,2]
θu = Float64[4,4]

bc = BoundConstraints(θl,θu)

result = solve(nls, θ_init, bc, conf)
LevenbergMarquardt_Result{Float64}(true, 5, 1.0000000000000002, [2.0, 3.999999999974335])

Using solver result

The solve() function returns a Abstract_BC_Solver_Result sub-typed structure that contains algorithm result.

In peculiar you can check if the method has converged

@assert converged(result)

and get the optimal θ

θ_solution = solution(result)
2-element ReadOnlyArrays.ReadOnlyArray{Float64, 1, Vector{Float64}}:
 2.0
 3.999999999974335