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 initialised 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 initialised 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.
• 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
• minimum_acceptable_model_improvement::Real = eps(number_eltype(p))

See also

gradient_descent, LevenbergMarquardt

AbstractLinearSurrogateObjective{E <: AbstractEvaluationType, O <: AbstractManifoldObjective{E}} <: AbstractManifoldObjective{E}

Provide a linear surrogate model for the given AbstractManifoldObjective O of the form

\[μ_p(X) = \frac{1}{2}\lVert \mathcal{L}(X) + y \rVert_{2}^2 + \frac{λ}{2}\lVert X \rVert_{p}^2, \qquad\text{ for }X ∈ T_{p}\mathcal{M}, λ ≥ 0,\]

where $\mathcal{L}$ is a linear operator on the tangent space at a point $p ∈ M$ that maps into some vector space $V$ and $y ∈ V$ is a fixed vector in that space and $\lVert ⋅ \rVert$ is a norm on $V$.

Both $\mathcal{L}$ and $y$ are derived from the objective O and usually depend on the base point $p ∈ M$.

Besides the usual methods defined for AbstractManifoldObjective that may be implemented like get_cost and get_gradient, the following methods should be implemented for a concrete subtype of AbstractLinearSurrogateObjective

• get_linear_operator to compute/evaluate the linear operator $\mathcal{L}$
• get_vector_field to compute/evaluate the vector $y$
• get_objective to provide access to the underlying objective O

See also the NormalEquationsObjective for the corresponding normal equations.

AbstractLevenbergMarquardtLinearSurrogateObjective{E<:AbstractEvaluationType}

Abstract supertype for Levenberg-Marquardt surrogates like LevenbergMarquardtLinearSurrogateObjective and LevenbergMarquardtLinearSurrogateCoordinatesObjective.

LevenbergMarquardtLinearSurrogateObjective{E<:AbstractEvaluationType, VF<:AbstractManifoldFirstOrderObjective{E}, R} <: AbstractLevenbergMarquardtLinearSurrogateObjective{E}

Given an ManifoldNonlinearLeastSquaresObjective objective and a penalty $λ$, this objective represents the penalized objective for the sub-problem to solve within every step of the Levenberg-Marquardt algorithm following the ideas of [TMHF00] given by

\[μ_p(X) = \frac{1}{2}\lVert \mathcal{L}(X) + y \rVert_{2}^2 + \frac{λ}{2}\lVert X \rVert_{p}^2, \qquad\text{ for }X ∈ T_{p}\mathcal{M}, λ ≥ 0,\]

where $X ∈ T_{p}\mathcal{M}$, $λ ≥ 0$ is the damping or penalty term, $\mathcal{L}: T_{p}\mathcal{M} → ℝ^n$ is a linear operator, and $y = y(p) ∈ ℝ^n$ is a vector field. For the derivation of the Riemannian case, see [BB26].

In order to build a surrogate also for the robustified Levenberg-Marquardt, introduce $α = 1 - \sqrt{1 + 2 \frac{ρ''(p)}{ρ'(p)}\lVert F(p) \rVert_{2}^2}$ and set $y = \frac{\sqrt{ρ'(p)}}{1-α}F(p)$ and $\mathcal{L}(X) = CJ_F(p)[X]$ with

\[C = \sqrt{ρ'(p)}(I-αP), \qquad P = \frac{F(p)F(p)^\mathrm{T}}{\lVert F(p) \rVert_{2}^2},\]

where $F(p) ∈ ℝ^n$ is the vector of residuals at point $p ∈ M$ and $J_F(p): T_{p}\mathcal{M} → ℝ^n$ is the Jacobian. These two can be accessed with get_vector_field for $y$ and get_linear_operator for $\mathcal{L}$, respectively. For technical details on the scaling using $α$, especially how the threshold and mode act as safeguards, see get_LevenbergMarquardt_scaling

Fields

• objective: the ManifoldNonlinearLeastSquaresObjective to penalize
• penalty::Real: the damping term $λ$
• threshold::Real: threshold $ε$ for stabilization of $α$ as $α ≤ 1-ε$, see get_LevenbergMarquardt_scaling
• mode::Symbol: which ode to use to stabilize α, see the internal helper get_LevenbergMarquardt_scaling
• value_cache: a vector to store the residuals $F(p)$ at the current point p internally to avoid re-computations

Constructor

LevenbergMarquardtLinearSurrogateObjective(objective; penalty::Real = 1e-6, threshold::Real = 1e-4, mode::Symbol = :Strict)

LevenbergMarquardtLinearSurrogateCoordinatesObjective{E<:AbstractEvaluationType, VF<:AbstractManifoldFirstOrderObjective{E}, R} <: AbstractLevenbergMarquardtLinearSurrogateObjective{E}

A subobjective similar to LevenbergMarquardtLinearSurrogateObjective but which uses coordinate-based Jacobians in a single, selected basis instead of being centered around linear operators.

Fields

• objective: the ManifoldNonlinearLeastSquaresObjective to penalize
• penalty::Real: the damping term $λ$
• threshold::Real: stabilization $ε$ for $α ≤ 1-ε$ in the rescaling of the residual and jacobian, see get_LevenbergMarquardt_scaling
• mode::Symbol: which mode to use to stabilize α, see the internal helper get_LevenbergMarquardt_scaling
• value_cache: a vector to store the residuals $F(p)$ at the current point p internally to avoid recomputations
• jacobian_cache: a vector to store the coordinate-based Jacobian of the residuals at the current point p internally to avoid recomputations. If the Jacobian is used as a linear operator, this is just a vector of nothings.

Constructor

LevenbergMarquardtLinearSurrogateCoordinatesObjective(objective; penalty::Real = 1e-6, threshold::Real = 1e-4, mode::Symbol = :Strict)

NormalEquationsObjective{E <: AbstractEvaluationType, O<: AbstractLinearSurrogateObjective{E}} <: AbstractSymmetricLinearSystemObjective{E}

A AbstractLinearSurrogateObjective might be overdetermined, and it usually is overdetermined, e.g. for the case of the LevenbergMarquardt algorithm. For this case, one considers the normal equations.

This wrapper provides the same three functions as the wrapped surrogate

• get_linear_operator to compute/evaluate the linear operator $\mathcal{L}$
• get_vector_field to compute/evaluate the vector $y$
• get_objective to provide access to the underlying surrogate

so that we obtain a symmetric linear system of equations, that can be

• solved with an iterative method like conjugate_gradient_descent or conjugate_residual
• solved as a linear system in a basis of the corresponding tangent space.

CoordinatesNormalSystemState <: AbstractManoptSolverState

A solver state indicating that we solve the LevenbergMarquardtLinearSurrogateObjective using a linear system in coordinates of the tangent space at the current iterate

Fields

• A an $n×n$ matrix to store the normal equations linear from get_linear_operator in coordinates, where n is the number of coordinates
• b a $n$ vector storing the right hand side of the normal equations in coordinates
• basis::AbstractBasis
• linsolve a functor (A,b) -> c to solve the linear system or (c, A, b) -> c depending on the evaluation type specified in solve!

Constructor

CoordinatesNormalSystemState(
    M::AbstractManifold, p = rand(M);
    evaluation = InplaceEvaluation(),
    linsolve = default_lm_lin_solve!,
    basis = DefaultOrthonormalBasis(),
    A = nothing
)

Construct the state, where not providing a memory for A uses the eltype of p to determine the element type of the matrix to store.

residual_scaling, operator_scaling = get_LevenbergMarquardt_scaling(ρ_prime::Real, ρ_double_prime::Real, FSq::Real, threshold::Real=1.0e-5, mode::Symbol=:Strict)

Compute the scalings for the residual $y$ and within the operator $C$ that are required for the robust rescaling within LevenbergMarquardts get_vector_field and get_linear_operator, respectively. Here FSq denotes $s = \lVert F(p) \rVert_{2}^2$ of the residual vector function $F$ evaluated at some point $p$, and ρ_prime=ρ'(s) and ρ_double_prime=ρ''(s) denote the current AbstractRobustifierFunctions first and second derivative evaluated at $s$.

The value for $α$ is given by

\[ α = 1 - \sqrt{1 + 2\frac{ρ''(s)}{ρ'(s)}s}\]

and hence the scaling of the residual and the within the projection of the operator are

\[\frac{\sqrt{ρ'(s)}}{1-α} \qquad\text{ and }\qquad \begin{cases} \frac{α}{s} & \text{ if } s ≠ 0\\\\ 0 & \text{ else,}\end{cases}\]

respectively.

Numerical stability

For a unique solution that is a minimizer in a Levenberg-Marquardt step, we require α < 1 and [TMHF00] recommends to bound this even by $1-ε$ for some threshold $ε > 0$.

Furthermore if $ρ'(s) + 2ρ''(s)⋅s ≤ 0$ the Hessian is also indefinite. This can be caught by making sure the argument of the $√$ is ensured to be non-negative.

The Ceres solver even omits the second term in the square root already if $ρ(s)'' < 0$ for stability reason, which means setting $α = 0$. In the case $s = 0$ we also set the operator scaling $α / s = 0$.

This function offers two modes

• :Normal keeps negative $ρ''(s) < 0$ but makes sure the square root is well-defined.
• :Strict (default) set $α = 0$ when $ρ''(s) < 0$ or when $s = 0$

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.

ZeroTangentVector

A small internal helper type to represent the zero tangent vector with two advantages

• we avoid to allocate it
• we can dispatch on it