Levenberg-Marquardt

LevenbergMarquardt(M, f, jacobian_f, p, num_components=-1; kwargs...)
LevenbergMarquardt(M, vgf, p; kwargs...)
LevenbergMarquardt(M, nlso, p; kwargs...)
LevenbergMarquardt!(M, f, jacobian_f, p, num_components=-1; kwargs...)
LevenbergMarquardt!(M, vgf, p; kwargs...)
LevenbergMarquardt!(M, nlso, p; kwargs...)

compute the the Riemannian Levenberg-Marquardt algorithm [Pee93, AOT22, BB26] to solve

\[\operatorname*{arg\,min}_{p ∈ \mathcal{M}} f(p), \qquad f(p) = \frac{1}{2} \sum_{i=1}^{m} ρ_i \bigl( \lVert F_i(p) \rVert^2 \bigr)\]

where $F_i: \mathcal{M} → ℝ^{n_i}$ is the $i$th block component of length $n_i > 0$ and each $ρ_i: ℝ → ℝ$ is a robustifier function, cf. AbstractRobustifierFunction, for each such a block component. The overall residual function is denoted by $F: \mathcal{M} → ℝ^{n}$ with $n = \sum_{i=1}^{m} n_i$ and concatenates all block components.

The second block of signatures perform the optimization in-place of p.

The regularization parameter is updated using a generalized scheme proposed in [Fan06], Eq. (2.2). See also [Yua15] for other schemes. The generalized scheme offers separate thresholds for the acceptance of new points (candidate_acceptance_threshold), decreasing the regularization parameter (damping_reduction_threshold) and increasing the regularization parameter (damping_increase_threshold).

Input

• M::AbstractManifold: a Riemannian manifold $\mathcal{M}$
• f: a residual function $f: \mathcal{M}→ℝ^m$. The residual function can be provided in two different ways
  • as a single function returning a vector $f(p) ∈ ℝ^m$
  • as a vector of functions, where each single function returns a scalar $f_i(p) ∈ ℝ$
  The type is determined by the function_type= keyword argument.
• jacobian_f: the Jacobian of $f$. The Jacobian can be provided in three different ways
  • as a single function returning a vector of gradient vectors $\bigl(\operatorname{grad} f_i(p)\bigr)_{i=1}^m$
  • as a vector of functions, where each single function returns a gradient vector $\operatorname{grad} f_i(p)$, $i=1,…,m$
  • as a single function returning a (coefficient) matrix $J ∈ ℝ^{m×d}$, where $d$ is the dimension of the manifold.
  These coefficients are given with respect to an AbstractBasis of the tangent space at p. The type is determined by the jacobian_type= keyword argument.
• p::P: a point on the manifold $\mathcal{M}$
• num_components: length $m$ of the vector returned by the cost function. By default its value is -1 which means that it is determined automatically by calling f one additional time. This is only possible when evaluation is AllocatingEvaluation, for mutating evaluation this value must be explicitly specified.

You can also provide the cost and its Jacobian already as a VectorGradientFunction vgf to indicate you are working on a single block, Alternatively, passing a ManifoldNonlinearLeastSquaresObjective nlso also works.

Keyword arguments

If you provide f and its jacobian

• evaluation::AbstractEvaluationType=AllocatingEvaluation(): specify whether the functions that return an array, for example a point or a tangent vector, work by allocating its result (AllocatingEvaluation) or whether they modify their input argument to return the result therein (InplaceEvaluation). Since usually the first argument is the manifold, the modified argument is the second.
• function_type=FunctionVectorialType: an AbstractVectorialType specifying the type of cost function provided.
• jacobian_type=FunctionVectorialType: an AbstractVectorialType specifying the type of Jacobian provided.

as well as then these are already combined in a single VectorGradientFunction vgf

• robustifier::AbstractRobustifierFunction=IdentityRobustifier(): for the robust variant, specify how the robustification is meant to take place.
  • if you provide a single vectorial function and its Jacobian, a single robustifer is applied to every component function of this vectorial function (each component is a block in the sum)
  • if you provide a vector of VectorGradientFunctions, each needs a robustifier.
• evaluation::AbstractEvaluationType=AllocatingEvaluation(): specify whether the functions that return an array, for example a point or a tangent vector, work by allocating its result (AllocatingEvaluation) or whether they modify their input argument to return the result therein (InplaceEvaluation). Since usually the first argument is the manifold, the modified argument is the second.

as well as in general using the model improvement parameter $m_k$ in several places, cf [BB26]

• candidate_acceptance_threshold=0.2: sufficient model improvement $η ∈ (0,1)$, i.e. $m_k > η$ to accept a candidate point

• damping_increase_factor=5.0: factor $β_{\text{i}}$ to increase damping, when the model is inaccurate

• damping_increase_threshold=candidate_acceptance_threshold: threshold $η_{\text{l}}$ the value $m_k$has to be below to increase damping. The default yields, that we increase damping when we reject a candidate.

• damping_reduction_factor= 1 / damping_increase_factor: factor $β_{\text{d}}$ to reduce damping, when the model is accurate

• damping_reduction_threshold=Inf: threshold $β_{\text{d}}$ to reduce damping, when the model is accurate The default means, that we never reduce damping.

• damping_term_min = 0.1: lower bound $μ_{\text{l}}$ for the damping $μ_k$ throughout the iterations

• damping_term_max = Inf: upper bound $μ_{\text{u}}$ for the damping $μ_k$ throughout the iterations

• initial_damping_term=damping_term_min: initial damping $μ_0$

• initial_residual_values = zeros(m): a cache for the vector of residuals, m is the number of residual blocks

• initial_jacobian_matrices: a cache for the evaluated Jacobians (currently only used if use_unified_basis = true, then initialized to a vector of jacobian matrices, otherwise ignored)

• retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)): a retraction $\operatorname{retr}$ to use, see the section on retractions

• scaling_threshold = 1.0e-6: a threshold ε to bound the scaling parameter α in the robust case away from 1, see get_LevenbergMarquardt_scaling

• scaling_mode = :Strict: specify the scaling stabilization mode, see get_LevenbergMarquardt_scaling

• stopping_criterion::StoppingCriterion=StopAfterIteration(500)|StopWhenGradientNormLess(1.0e-12)|StopWhenStepsizeLess(1.0e-12): a functor indicating that the stopping criterion is fulfilled

• sub_objective : specify the objective for the surrogate sub problem to solver in every iteration. This is set depending on the use_unified_basis

  • if true to the LevenbergMarquardtLinearSurrogateCoordinatesObjective which always works in coordinates of one single basis per tangent space and allows to cache Jacobian evaluations
  • if false to the LevenbergMarquardtLinearSurrogateObjective that can work either with linear operators or in (even different) coordinates.

  This keyword is ignored if you set the sub_problem and/or sub_state keyword directly

• sub_problem =DefaultManoptProblem(TangentSpace(M, p), sub_objective): specify the sub problem to be solved. This should usually be phrased on the tangent space at the current iterate

• sub_state =ConjugateResidualState(TangentSpace(M, p), sub_objective): specify the solver for the surrogate, see also conjugate_residual

• use_unified_basis = false: specify to use a single basis for all Jacobian evaluations at a certain iterate, see sub_objective this requires that all Jacobians involved are of type CoefficientVectorialType, since only then a jacobian can be represented as a matrix, and then here unified in the sense that all use the same basis.

All other keyword arguments are passed to decorate_state! for state decorators or decorate_objective! for objective decorators, respectively.

Output

The obtained approximate minimizer $p^*$. To obtain the whole final state of the solver, see get_solver_return for details, especially the return_state= keyword.

LevenbergMarquardt(M, f, jacobian_f, p, num_components=-1; kwargs...)
LevenbergMarquardt(M, vgf, p; kwargs...)
LevenbergMarquardt(M, nlso, p; kwargs...)
LevenbergMarquardt!(M, f, jacobian_f, p, num_components=-1; kwargs...)
LevenbergMarquardt!(M, vgf, p; kwargs...)
LevenbergMarquardt!(M, nlso, p; kwargs...)

compute the the Riemannian Levenberg-Marquardt algorithm [Pee93, AOT22, BB26] to solve

\[\operatorname*{arg\,min}_{p ∈ \mathcal{M}} f(p), \qquad f(p) = \frac{1}{2} \sum_{i=1}^{m} ρ_i \bigl( \lVert F_i(p) \rVert^2 \bigr)\]

where $F_i: \mathcal{M} → ℝ^{n_i}$ is the $i$th block component of length $n_i > 0$ and each $ρ_i: ℝ → ℝ$ is a robustifier function, cf. AbstractRobustifierFunction, for each such a block component. The overall residual function is denoted by $F: \mathcal{M} → ℝ^{n}$ with $n = \sum_{i=1}^{m} n_i$ and concatenates all block components.

The second block of signatures perform the optimization in-place of p.

The regularization parameter is updated using a generalized scheme proposed in [Fan06], Eq. (2.2). See also [Yua15] for other schemes. The generalized scheme offers separate thresholds for the acceptance of new points (candidate_acceptance_threshold), decreasing the regularization parameter (damping_reduction_threshold) and increasing the regularization parameter (damping_increase_threshold).

Input

• M::AbstractManifold: a Riemannian manifold $\mathcal{M}$
• f: a residual function $f: \mathcal{M}→ℝ^m$. The residual function can be provided in two different ways
  • as a single function returning a vector $f(p) ∈ ℝ^m$
  • as a vector of functions, where each single function returns a scalar $f_i(p) ∈ ℝ$
  The type is determined by the function_type= keyword argument.
• jacobian_f: the Jacobian of $f$. The Jacobian can be provided in three different ways
  • as a single function returning a vector of gradient vectors $\bigl(\operatorname{grad} f_i(p)\bigr)_{i=1}^m$
  • as a vector of functions, where each single function returns a gradient vector $\operatorname{grad} f_i(p)$, $i=1,…,m$
  • as a single function returning a (coefficient) matrix $J ∈ ℝ^{m×d}$, where $d$ is the dimension of the manifold.
  These coefficients are given with respect to an AbstractBasis of the tangent space at p. The type is determined by the jacobian_type= keyword argument.
• p::P: a point on the manifold $\mathcal{M}$
• num_components: length $m$ of the vector returned by the cost function. By default its value is -1 which means that it is determined automatically by calling f one additional time. This is only possible when evaluation is AllocatingEvaluation, for mutating evaluation this value must be explicitly specified.

You can also provide the cost and its Jacobian already as a VectorGradientFunction vgf to indicate you are working on a single block, Alternatively, passing a ManifoldNonlinearLeastSquaresObjective nlso also works.

Keyword arguments

If you provide f and its jacobian

• evaluation::AbstractEvaluationType=AllocatingEvaluation(): specify whether the functions that return an array, for example a point or a tangent vector, work by allocating its result (AllocatingEvaluation) or whether they modify their input argument to return the result therein (InplaceEvaluation). Since usually the first argument is the manifold, the modified argument is the second.
• function_type=FunctionVectorialType: an AbstractVectorialType specifying the type of cost function provided.
• jacobian_type=FunctionVectorialType: an AbstractVectorialType specifying the type of Jacobian provided.

as well as then these are already combined in a single VectorGradientFunction vgf

• robustifier::AbstractRobustifierFunction=IdentityRobustifier(): for the robust variant, specify how the robustification is meant to take place.
  • if you provide a single vectorial function and its Jacobian, a single robustifer is applied to every component function of this vectorial function (each component is a block in the sum)
  • if you provide a vector of VectorGradientFunctions, each needs a robustifier.
• evaluation::AbstractEvaluationType=AllocatingEvaluation(): specify whether the functions that return an array, for example a point or a tangent vector, work by allocating its result (AllocatingEvaluation) or whether they modify their input argument to return the result therein (InplaceEvaluation). Since usually the first argument is the manifold, the modified argument is the second.

as well as in general using the model improvement parameter $m_k$ in several places, cf [BB26]

• candidate_acceptance_threshold=0.2: sufficient model improvement $η ∈ (0,1)$, i.e. $m_k > η$ to accept a candidate point

• damping_increase_factor=5.0: factor $β_{\text{i}}$ to increase damping, when the model is inaccurate

• damping_increase_threshold=candidate_acceptance_threshold: threshold $η_{\text{l}}$ the value $m_k$has to be below to increase damping. The default yields, that we increase damping when we reject a candidate.

• damping_reduction_factor= 1 / damping_increase_factor: factor $β_{\text{d}}$ to reduce damping, when the model is accurate

• damping_reduction_threshold=Inf: threshold $β_{\text{d}}$ to reduce damping, when the model is accurate The default means, that we never reduce damping.

• damping_term_min = 0.1: lower bound $μ_{\text{l}}$ for the damping $μ_k$ throughout the iterations

• damping_term_max = Inf: upper bound $μ_{\text{u}}$ for the damping $μ_k$ throughout the iterations

• initial_damping_term=damping_term_min: initial damping $μ_0$

• initial_residual_values = zeros(m): a cache for the vector of residuals, m is the number of residual blocks

• initial_jacobian_matrices: a cache for the evaluated Jacobians (currently only used if use_unified_basis = true, then initialized to a vector of jacobian matrices, otherwise ignored)

• retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)): a retraction $\operatorname{retr}$ to use, see the section on retractions

• scaling_threshold = 1.0e-6: a threshold ε to bound the scaling parameter α in the robust case away from 1, see get_LevenbergMarquardt_scaling

• scaling_mode = :Strict: specify the scaling stabilization mode, see get_LevenbergMarquardt_scaling

• stopping_criterion::StoppingCriterion=StopAfterIteration(500)|StopWhenGradientNormLess(1.0e-12)|StopWhenStepsizeLess(1.0e-12): a functor indicating that the stopping criterion is fulfilled

• sub_objective : specify the objective for the surrogate sub problem to solver in every iteration. This is set depending on the use_unified_basis

  • if true to the LevenbergMarquardtLinearSurrogateCoordinatesObjective which always works in coordinates of one single basis per tangent space and allows to cache Jacobian evaluations
  • if false to the LevenbergMarquardtLinearSurrogateObjective that can work either with linear operators or in (even different) coordinates.

  This keyword is ignored if you set the sub_problem and/or sub_state keyword directly

• sub_problem =DefaultManoptProblem(TangentSpace(M, p), sub_objective): specify the sub problem to be solved. This should usually be phrased on the tangent space at the current iterate

• sub_state =ConjugateResidualState(TangentSpace(M, p), sub_objective): specify the solver for the surrogate, see also conjugate_residual

• use_unified_basis = false: specify to use a single basis for all Jacobian evaluations at a certain iterate, see sub_objective this requires that all Jacobians involved are of type CoefficientVectorialType, since only then a jacobian can be represented as a matrix, and then here unified in the sense that all use the same basis.

All other keyword arguments are passed to decorate_state! for state decorators or decorate_objective! for objective decorators, respectively.

Output

The obtained approximate minimizer $p^*$. To obtain the whole final state of the solver, see get_solver_return for details, especially the return_state= keyword.

LevenbergMarquardtState{P,T} <: AbstractGradientSolverState

Describes a Gradient based descent algorithm, with

Fields

• damping_term: current value of the damping term
• damping_term_min: lower bound for the damping term
• damping_term_max: upper bound for the damping term
• damping_increase_factor: improvement quotient exceeds damping_reduction_threshold.
• damping_reduction_threshold: threshold for the improvement quotient above which the damping term is reduced by multiplying it with β_reduction.
• damping_increase_threshold : threshold for the improvement quotient below which the damping term is increased by multiplying it with β.
• direction: the current search direction, which is the solution of the linearized subproblem in each iteration.
• candidate_acceptance_threshold: Scaling factor for the sufficient cost decrease threshold required to accept new proposal points. Allowed range: 0 < η < 1.
• callbacks: the callbacks dictionary
• jacobian_matrices: the current Jacobian of $F$ in matrix form per block, hence a vector of matrices. This is (by default) set to nothing if another representation is used.
• minimum_acceptable_model_improvement: the minimum improvement in the model function that is required to accept a new point; if this is not met, the new point is rejected and the damping term is increased.
• p::P: a point on the manifold $\mathcal{M}$ storing the current iterate
• retraction_method::AbstractRetractionMethod: a retraction $\operatorname{retr}$ to use, see the section on retractions
• residual_values: values of the residuals calculated in the solver setup or the previous iteration
• stop::StoppingCriterion: a functor indicating that the stopping criterion is fulfilled
• sub_problem::Union{AbstractManoptProblem, F}: specify a problem for a solver or a closed form solution function, which can be allocating or in-place.
• sub_state::Union{AbstractManoptProblem, F}: a state to specify the sub solver to use. For a closed form solution, this indicates the type of function.
• X::T: a tangent vector at the point $p$ on the manifold $\mathcal{M}$

Constructor

LevenbergMarquardtState(M, sub_problem, sub_state, initial_residual_values, initial_jacobian; kwargs...)

Generate the Levenberg-Marquardt solver state.

Keyword arguments

The following fields are keyword arguments

• candidate_acceptance_threshold = 0.2,
• damping_increase_factor = 5.0
• damping_reduction_factor = 0.5
• damping_term_min = 0.1
• damping_term_max = Inf
• damping_term = damping_term_min
• damping_reduction_threshold = Inf
• damping_increase_threshold = candidate_acceptance_threshold
• direction = copy(M, p, X)
• p =rand(M)
• X =zero_vector(M, p)
• retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)): a retraction $\operatorname{retr}$ to use, see the section on retractions
• stopping_criterion::StoppingCriterion=StopAfterIteration(200)|StopWhenGradientNormLess(1e-12)|StopWhenStepsizeLess(1e-12): a functor indicating that the stopping criterion is fulfilled
• callbacks::D = Dict{Symbol,Function}(): provided callback functions given as a dictionary with symbols as keys
• minimum_acceptable_model_improvement::Real = eps(number_eltype(p))

See also

gradient_descent, LevenbergMarquardt

default_lm_lin_solve!(sk, JJ::AbstractMatrix, grad_f_c)

Solve the linear system of equations of the normal equations JJ \ grad_f_c where JJ is a symmetric positive definite matrix and save the result to sk. In case of numerical errors the PosDefException is caught and the default symmetric solver (Symmetric(JJ) \ grad_f_c) is used.

The function is intended to be used with LevenbergMarquardt.

add_normal_vector_field!(M::AbstractManifold, X, o::AbstractFirstOrderVectorFunction, r::AbstractRobustifierFunction, p)
add_normal_vector_field!(M::AbstractManifold, c, o::AbstractFirstOrderVectorFunction, r::AbstractRobustifierFunction, p, B::AbstractBasis)

Add the contribution of o / r to the normal linear operator tangent vector in X or c. See get_normal_vector_field for the mathematical details. Note that this is done per every block (vectorial function with its robustifier) of the underlying ManifoldNonlinearLeastSquaresObjective and summed up.

See also get_normal_linear_operator for evaluating the corresponding linear operator of the (normal) linear system, and get_LevenbergMarquardt_scaling for details on the scaling and computation of $C$.

add_normal_linear_operator!(
    M::AbstractManifold, A::AbstractMatrix, o::AbstractFirstOrderVectorFunction,
    r::AbstractRobustifierFunction, p, basis::AbstractBasis;
    value_cache = get_value(M, o, p), ε::Real, mode::Symbol
)

Add the contribution of a single block (vectorial function with its robustifier) to the linear normal operator, i.e. compute $A += J_F^*(p)[C^T C J_F(p)[X]]$ in-place of A for the given block. See get_normal_linear_operator for details

add_linear_operator_coord!(
    M::AbstractManifold, y::AbstractVector, lmsco::LevenbergMarquardtLinearSurrogateCoordinatesObjective, p, cX::AbstractVector
)

Add the (Triggs correction, residual-like) linear operator corresponding to the lmsco surrogate to vector y. It is assumed that lmsco.value_cache has been filled in step_solver! of LevenbergMarquardt, so we can just use it here.