Conjugate residual solver in a tangent space

conjugate_residual(TpM::TangentSpace, A, b, X=zero_vector(TpM))
conjugate_residual(TpM::TangentSpace, slso::AbstractSymmetricLinearSystemObjective, X=zero_vector(TpM))
conjugate_residual!(TpM::TangentSpace, A, b, X)
conjugate_residual!(TpM::TangentSpace, slso::AbstractSymmetricLinearSystemObjective, X)

Compute the solution of $\mathcal{A}(p)[X] + b(p) = 0_p$, where

• $\mathcal{A}$ is a linear, symmetric operator on $T_{p}\mathcal{M}$
• $b$ is a vector field on the manifold
• $X ∈ T_{p}\mathcal{M}$ is a tangent vector
• $0_p$ is the zero vector $T_{p}\mathcal{M}$.

This implementation follows Algorithm 3 in [LY24] and is initialized with $X^{(0)}$ as the given tangent vector X (or as the zero vector if warm_start=false) and

• the initial residual $r^{(0)} = -b(p) - \mathcal{A}(p)[X^{(0)}]$
• the initial conjugate direction $d^{(0)} = r^{(0)}$
• initialize $Y^{(0)} = \mathcal{A}(p)[r^{(0)}]$

It then performs the following steps at iteration $k=0,…$ until the stopping_criterion is fulfilled.

1. compute a step size $α_k = \displaystyle\frac{⟨ r^{(k)}, \mathcal{A}(p)[r^{(k)}] ⟩_p}{⟨ \mathcal{A}(p)[d^{(k)}], \mathcal{A}(p)[d^{(k)}] ⟩_p}$
2. do a step $X^{(k+1)} = X^{(k)} + α_kd^{(k)}$
3. update the residual $r^{(k+1)} = r^{(k)} - α_k Y^{(k)}$
4. compute $Z = \mathcal{A}(p)[r^{(k+1)}]$
5. Update the conjugate coefficient $β_k = \displaystyle\frac{⟨ r^{(k+1)}, \mathcal{A}(p)[r^{(k+1)}] ⟩_p}{⟨ r^{(k)}, \mathcal{A}(p)[r^{(k)}] ⟩_p}$
6. Update the conjugate direction $d^{(k+1)} = r^{(k+1)} + β_kd^{(k)}$
7. Update $Y^{(k+1)} = Z + β_k Y^{(k)}$

Note that the right hand side of Step 7 is the same as evaluating $\mathcal{A}[d^{(k+1)}]$, but avoids the actual evaluation

Input

• TpM the TangentSpace as the domain
• A a symmetric linear operator on the tangent space (M, p, X) -> Y
• b a vector field on the manifold (M, p) -> X
• X the initial tangent vector

Keyword arguments

• 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.
• 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
• stopping_criterion::StoppingCriterion=StopAfterIteration(manifold_dimension(M))|StopWhenRelativeResidualLess(c,1e-8): a functor indicating that the stopping criterion is fulfilled where $c = \lVert b \rVert$ is the norm of the vector field b at p.
• warm_start=true: whether to reuse the initial X to warm start the solver (true) or to start from the zero vector (false), see ConjugateResidualState.

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.

conjugate_residual(TpM::TangentSpace, A, b, X=zero_vector(TpM))
conjugate_residual(TpM::TangentSpace, slso::AbstractSymmetricLinearSystemObjective, X=zero_vector(TpM))
conjugate_residual!(TpM::TangentSpace, A, b, X)
conjugate_residual!(TpM::TangentSpace, slso::AbstractSymmetricLinearSystemObjective, X)

Compute the solution of $\mathcal{A}(p)[X] + b(p) = 0_p$, where

• $\mathcal{A}$ is a linear, symmetric operator on $T_{p}\mathcal{M}$
• $b$ is a vector field on the manifold
• $X ∈ T_{p}\mathcal{M}$ is a tangent vector
• $0_p$ is the zero vector $T_{p}\mathcal{M}$.

This implementation follows Algorithm 3 in [LY24] and is initialized with $X^{(0)}$ as the given tangent vector X (or as the zero vector if warm_start=false) and

• the initial residual $r^{(0)} = -b(p) - \mathcal{A}(p)[X^{(0)}]$
• the initial conjugate direction $d^{(0)} = r^{(0)}$
• initialize $Y^{(0)} = \mathcal{A}(p)[r^{(0)}]$

It then performs the following steps at iteration $k=0,…$ until the stopping_criterion is fulfilled.

1. compute a step size $α_k = \displaystyle\frac{⟨ r^{(k)}, \mathcal{A}(p)[r^{(k)}] ⟩_p}{⟨ \mathcal{A}(p)[d^{(k)}], \mathcal{A}(p)[d^{(k)}] ⟩_p}$
2. do a step $X^{(k+1)} = X^{(k)} + α_kd^{(k)}$
3. update the residual $r^{(k+1)} = r^{(k)} - α_k Y^{(k)}$
4. compute $Z = \mathcal{A}(p)[r^{(k+1)}]$
5. Update the conjugate coefficient $β_k = \displaystyle\frac{⟨ r^{(k+1)}, \mathcal{A}(p)[r^{(k+1)}] ⟩_p}{⟨ r^{(k)}, \mathcal{A}(p)[r^{(k)}] ⟩_p}$
6. Update the conjugate direction $d^{(k+1)} = r^{(k+1)} + β_kd^{(k)}$
7. Update $Y^{(k+1)} = Z + β_k Y^{(k)}$

Note that the right hand side of Step 7 is the same as evaluating $\mathcal{A}[d^{(k+1)}]$, but avoids the actual evaluation

Input

• TpM the TangentSpace as the domain
• A a symmetric linear operator on the tangent space (M, p, X) -> Y
• b a vector field on the manifold (M, p) -> X
• X the initial tangent vector

Keyword arguments

• 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.
• 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
• stopping_criterion::StoppingCriterion=StopAfterIteration(manifold_dimension(M))|StopWhenRelativeResidualLess(c,1e-8): a functor indicating that the stopping criterion is fulfilled where $c = \lVert b \rVert$ is the norm of the vector field b at p.
• warm_start=true: whether to reuse the initial X to warm start the solver (true) or to start from the zero vector (false), see ConjugateResidualState.

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.

ConjugateResidualState{T,R,TStop<:StoppingCriterion,C<:AbstractDict{Symbol}} <: AbstractManoptSolverState

A state for the conjugate_residual solver.

Fields

• X::T: the iterate
• callbacks::C: the callbacks dictionary
• r::T: the residual $r = -b(p) - \mathcal{A}(p)[X]$
• d::T: the conjugate direction
• Ar::T, Ad::T: storages for $\mathcal{A}(p)[r]$, $\mathcal{A}(p)[d]$
• rAr::R: internal field for storing $⟨ r, \mathcal{A}(p)[r] ⟩$
• α::R: a step length
• β::R: the conjugate coefficient
• stop::StoppingCriterion: a functor indicating that the stopping criterion is fulfilled
• warm_start: whether to warm start or not when reusing this state, i.e.
  • true (default): means we reuse the values in X on initialization and set the remaining terms accordingly. This involves one call to the objective's linear system and right hand side.
  • false: Initialize X to the zero vector and hence d=r=-b(p), which saves one evaluation of the linear operator.

Constructor

ConjugateResidualState(TpM::TangentSpace, slso::AbstractSymmetricLinearSystemObjective; kwargs...)

Initialize the state with default values.

Keyword arguments

• r=-get_gradient(TpM, slso, X)
• d=copy(TpM, r)
• Ar=get_hessian(TpM, slso, X, r)
• Ad=copy(TpM, Ar)
• α::R=0.0
• β::R=0.0
• callbacks::D = Dict{Symbol,Function}(): provided callback functions given as a dictionary with symbols as keys
• stopping_criterion::StoppingCriterion=StopAfterIteration(manifold_dimension(M))|StopWhenGradientNormLess(1e-8): 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}$
• warm_start=true: whether to reuse the values in X for the initialization (true) or to start from the zero vector (false), see the field description above.

See also

conjugate_residual

StopWhenRelativeResidualLess <: StoppingCriterion

Stop when the relative residual in the conjugate_residual is below a certain threshold, i.e.

\[\displaystyle\frac{\lVert r^{(k)} \rVert}{c} < ε,\]

where $c = \lVert b \rVert$ of the initial vector from the vector field in $\mathcal{A}(p)[X] + b(p) = 0_p$, from the conjugate_residual

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.
• c: the initial norm
• ε: the threshold
• norm_r: the last computed norm of the residual

Constructor

StopWhenRelativeResidualLess(c, ε; norm_r = 2*c*ε)

Initialize the stopping criterion.