Steihaug-Toint truncated conjugate gradient method

truncated_conjugate_gradient_descent(M, f, grad_f, Hess_f, p=rand(M), X=rand(M; vector_at=p);
    kwargs...
)
truncated_conjugate_gradient_descent(M, mho::ManifoldHessianObjective, p=rand(M), X=rand(M; vector_at=p);
    kwargs...
)
truncated_conjugate_gradient_descent(TpM::TangentSpace, trmo::TrustRegionModelObjective, p, X;
    kwargs...
)

solve the trust-region subproblem

\[\begin{align*} \operatorname{arg\,min}_{Y ∈ T_p\mathcal{M}}&\ m_p(Y) = f(p) + ⟨\operatorname{grad}f(p), Y⟩_p + \frac{1}{2} ⟨\mathcal{H}_p[Y], Y⟩_p\\ \text{such that} & \ \lVert Y \rVert_p ≤ Δ \end{align*}\]

on a manifold $\mathcal{M}$ by using the Steihaug-Toint truncated conjugate-gradient (tCG) method. This can be done in-place of X.

For a description of the algorithm and theorems offering convergence guarantees, see [ABG06, CGT00].

Input

• M::AbstractManifold: a Riemannian manifold $\mathcal{M}$
• f: a cost function $f: \mathcal{M}→ ℝ$ implemented as (M, p) -> v
• grad_f: the (Riemannian) gradient $\operatorname{grad}f: \mathcal{M} → T\mathcal{M}$ of f as a function (M, p) -> X or a function (M, X, p) -> X computing X in-place
• Hess_f: the (Riemannian) Hessian $\operatorname{Hess}f: T_{p}\mathcal{M} → T_{p}\mathcal{M}$ of f as a function (M, p, X) -> Y or a function (M, Y, p, X) -> Y computing Y in-place
• p::P: a point on the manifold $\mathcal{M}$
• X::T: a tangent vector at the point $p$ on the manifold $\mathcal{M}$

Instead of the three functions, you either provide a ManifoldHessianObjective mho which is then used to build the trust region model, or a TrustRegionModelObjective trmo directly.

Keyword arguments

• callbacks::D = Dict{Symbol,Function}(): provided callback functions given either as a single function (symbol, problem, state, k) called in every hook or as a (vector of) pairs :hook => function, which are processed by process_callbacks_arg. As key you can either pass single symbol or an array of symbols to indicate a callback should be added in multiple places
• evaluation::AbstractEvaluationType=AllocatingEvaluation(): specify whether the functions that return an array, for example a point or a tangent vector, work by allocating their 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.
• preconditioner: a preconditioner for the Hessian H. This is either an allocating function (M, p, X) -> Y or an in-place function (M, Y, p, X) -> Y, see evaluation, and by default set to the identity.
• θ=1.0: the superlinear convergence target rate of $1+θ$
• κ=0.1: the linear convergence target rate.
• project!=copyto!: for numerical stability it is possible to project onto the tangent space after every iteration. the function has to work inplace of Y, that is (M, Y, p, X) -> Y, where X and Y can be the same memory.
• randomize=false: whether to start the iteration from the given X (true, this disables preconditioning) or from the zero vector (false, X then only provides the memory for the result)
• stopping_criterion::StoppingCriterion=StopAfterIteration(manifold_dimension(base_manifold(TpM)))|StopWhenResidualIsReducedByFactorOrPower(; κ=κ, θ=θ)|StopWhenTrustRegionIsExceeded()|StopWhenCurvatureIsNegative()|StopWhenModelIncreased(): a functor indicating that the stopping criterion is fulfilled
• trust_region_radius=injectivity_radius(M) / 4 (or 1.0 if the injectivity radius is infinite): the initial trust-region radius

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.

See also

trust_regions

truncated_conjugate_gradient_descent(M, f, grad_f, Hess_f, p=rand(M), X=rand(M; vector_at=p);
    kwargs...
)
truncated_conjugate_gradient_descent(M, mho::ManifoldHessianObjective, p=rand(M), X=rand(M; vector_at=p);
    kwargs...
)
truncated_conjugate_gradient_descent(TpM::TangentSpace, trmo::TrustRegionModelObjective, p, X;
    kwargs...
)

solve the trust-region subproblem

\[\begin{align*} \operatorname{arg\,min}_{Y ∈ T_p\mathcal{M}}&\ m_p(Y) = f(p) + ⟨\operatorname{grad}f(p), Y⟩_p + \frac{1}{2} ⟨\mathcal{H}_p[Y], Y⟩_p\\ \text{such that} & \ \lVert Y \rVert_p ≤ Δ \end{align*}\]

on a manifold $\mathcal{M}$ by using the Steihaug-Toint truncated conjugate-gradient (tCG) method. This can be done in-place of X.

For a description of the algorithm and theorems offering convergence guarantees, see [ABG06, CGT00].

Input

• M::AbstractManifold: a Riemannian manifold $\mathcal{M}$
• f: a cost function $f: \mathcal{M}→ ℝ$ implemented as (M, p) -> v
• grad_f: the (Riemannian) gradient $\operatorname{grad}f: \mathcal{M} → T\mathcal{M}$ of f as a function (M, p) -> X or a function (M, X, p) -> X computing X in-place
• Hess_f: the (Riemannian) Hessian $\operatorname{Hess}f: T_{p}\mathcal{M} → T_{p}\mathcal{M}$ of f as a function (M, p, X) -> Y or a function (M, Y, p, X) -> Y computing Y in-place
• p::P: a point on the manifold $\mathcal{M}$
• X::T: a tangent vector at the point $p$ on the manifold $\mathcal{M}$

Instead of the three functions, you either provide a ManifoldHessianObjective mho which is then used to build the trust region model, or a TrustRegionModelObjective trmo directly.

Keyword arguments

• callbacks::D = Dict{Symbol,Function}(): provided callback functions given either as a single function (symbol, problem, state, k) called in every hook or as a (vector of) pairs :hook => function, which are processed by process_callbacks_arg. As key you can either pass single symbol or an array of symbols to indicate a callback should be added in multiple places
• evaluation::AbstractEvaluationType=AllocatingEvaluation(): specify whether the functions that return an array, for example a point or a tangent vector, work by allocating their 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.
• preconditioner: a preconditioner for the Hessian H. This is either an allocating function (M, p, X) -> Y or an in-place function (M, Y, p, X) -> Y, see evaluation, and by default set to the identity.
• θ=1.0: the superlinear convergence target rate of $1+θ$
• κ=0.1: the linear convergence target rate.
• project!=copyto!: for numerical stability it is possible to project onto the tangent space after every iteration. the function has to work inplace of Y, that is (M, Y, p, X) -> Y, where X and Y can be the same memory.
• randomize=false: whether to start the iteration from the given X (true, this disables preconditioning) or from the zero vector (false, X then only provides the memory for the result)
• stopping_criterion::StoppingCriterion=StopAfterIteration(manifold_dimension(base_manifold(TpM)))|StopWhenResidualIsReducedByFactorOrPower(; κ=κ, θ=θ)|StopWhenTrustRegionIsExceeded()|StopWhenCurvatureIsNegative()|StopWhenModelIncreased(): a functor indicating that the stopping criterion is fulfilled
• trust_region_radius=injectivity_radius(M) / 4 (or 1.0 if the injectivity radius is infinite): the initial trust-region radius

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.

See also

trust_regions

TruncatedConjugateGradientState <: AbstractHessianSolverState

describe the Steihaug-Toint truncated conjugate-gradient method, with

Fields

Let T denote the type of a tangent vector and R <: Real.

• callbacks::D: provided callback functions given as a dictionary with symbols as keys
• δ::T: the conjugate gradient search direction
• δHδ, YPδ, δPδ, YPY: temporary inner products with Hδ and preconditioned inner products.
• Hδ, HY: temporary results of the Hessian applied to δ and Y, respectively.
• project!: for numerical stability it is possible to project onto the tangent space after every iteration. the function has to work inplace of Y, that is (M, Y, p, X) -> Y, where X and Y can be the same memory.
• randomize: whether to start from the given tangent vector (true, without preconditioning) or from the zero vector (false)
• residual::T: the gradient of the model $m(Y)$
• stop::StoppingCriterion: a functor indicating that the stopping criterion is fulfilled
• trust_region_radius::R: the trust-region radius
• X::T: the gradient $\operatorname{grad}f(p)$
• Y::T: current iterate tangent vector
• z::T: the preconditioned residual
• z_r::R: inner product of the residual and z

Constructor

TruncatedConjugateGradientState(TpM::TangentSpace; X=rand(TpM), kwargs...)

Initialize the TCG state.

Input

• TpM: a TangentSpace

Keyword arguments

• callbacks::D = Dict{Symbol,Function}(): provided callback functions given as a dictionary with symbols as keys
• κ=0.1
• project!::F=copyto!: initialize the numerical stabilization to just copy the result
• randomize=false
• θ=1.0
• trust_region_radius=injectivity_radius(base_manifold(TpM)) / 4 (or 1.0 if the injectivity radius is infinite)
• stopping_criterion::StoppingCriterion=StopAfterIteration(manifold_dimension(base_manifold(TpM)))|StopWhenResidualIsReducedByFactorOrPower(; κ=κ, θ=θ)|StopWhenTrustRegionIsExceeded()|StopWhenCurvatureIsNegative()|StopWhenModelIncreased(): a functor indicating that the stopping criterion is fulfilled
• X::T =rand(TpM): a tangent vector at the point $p$ on the manifold $\mathcal{M}$

See also

truncated_conjugate_gradient_descent, trust_regions

StopWhenResidualIsReducedByFactorOrPower <: StoppingCriterion

A functor for testing if the norm of residual at the current iterate is reduced either by a power of 1+θ or by a factor κ compared to the norm of the initial residual. The criterion hence reads

$\lVert r_k \rVert_{p} ≦ \lVert r_0 \rVert_{p^{(0)}} \min \bigl( κ, \lVert r_0 \rVert_{p^{(0)}}^{θ} \bigr)$.

Fields

• κ: the reduction factor
• θ: part of the reduction power
• at_iteration::Int: an integer indicating at which iteration the stopping criterion last indicated to stop, which might also be before the solver started (0). Any negative value indicates that this was not yet the case.

Constructor

StopWhenResidualIsReducedByFactorOrPower(; κ=0.1, θ=1.0)

Initialize the StopWhenResidualIsReducedByFactorOrPower functor to indicate to stop after the norm of the current residual is less than or equal to the smaller of the norm of the initial residual to the power of 1+θ and the norm of the initial residual times κ.

See also

truncated_conjugate_gradient_descent, trust_regions

StopWhenTrustRegionIsExceeded <: StoppingCriterion

A stopping criterion to stop when next iterate is larger than the trust-region radius $Δ ≤ \lVert Y^{(k)} \rVert_{p^{(k)}}$.

This can for example be used in the Steihaug-Toint truncated conjugate gradient method as a subsolver for trust_regions.

Fields

• at_iteration::Int: an integer indicating at which iteration the stopping criterion last indicated to stop, which might also be before the solver started (0). Any negative value indicates that this was not yet the case.
• trr the trust region radius
• YPY the computed squared norm of $Y$.

Constructor

StopWhenTrustRegionIsExceeded()

initialize the StopWhenTrustRegionIsExceeded functor to indicate to stop after the norm of the next iterate is greater than the trust-region radius.

See also

truncated_conjugate_gradient_descent, trust_regions

StopWhenCurvatureIsNegative <: StoppingCriterion

A functor for testing if the curvature of the model is negative, $⟨δ_k, \operatorname{Hess} F(p)[δ_k]⟩_p ≦ 0$. In this case, the model is not strictly convex, and the stepsize as computed does not yield a reduction of the model.

Fields

• at_iteration::Int: an integer indicating at which iteration the stopping criterion last indicated to stop, which might also be before the solver started (0). Any negative value indicates that this was not yet the case.
• value store the value of the inner product.

Constructor

StopWhenCurvatureIsNegative()

See also

truncated_conjugate_gradient_descent, trust_regions

StopWhenModelIncreased <: StoppingCriterion

A functor for testing if the model value increased.

Fields

• at_iteration::Int: an integer indicating at which iteration the stopping criterion last indicated to stop, which might also be before the solver started (0). Any negative value indicates that this was not yet the case.
• model_value store the last model value
• inc_model_value store the model value that increased

Constructor

StopWhenModelIncreased()

See also

truncated_conjugate_gradient_descent, trust_regions

set_parameter!(c::StopWhenResidualIsReducedByFactorOrPower, :ResidualPower, v)

Update the residual Power θ to v.

set_parameter!(c::StopWhenResidualIsReducedByFactorOrPower, :ResidualFactor, v)

Update the residual Factor κ to v.