Conjugate Residual Solver in a Tangent space

conjugate_residual(TpM::TangentSpace, A, b, p=rand(TpM))
conjugate_residual(TpM::TangentSpace, slso::SymmetricLinearSystemObjective, p=rand(TpM))
conjugate_residual!(TpM::TangentSpace, A, b, p)
conjugate_residual!(TpM::TangentSpace, slso::SymmetricLinearSystemObjective, p)

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 initalised with $X^{(0)}$ as the zero vector 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)[X^{(0)}]$

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

1. compute a step size $α_k = \displaystyle\frac{\langle r^{(k)}, \mathcal A(p)[r^{(k)}] \rangle_p}{\langle \mathcal A(p)[d^{(k)}], \mathcal A(p)[d^{(k)}] \rangle_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{\langle r^{(k+1)}, \mathcal A(p)[r^{(k+1)}] \rangle_p}{\langle r^{(k)}, \mathcal A(p)[r^{(k)}] \rangle_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 tangent space (M, p) -> X

Keyword arguments

• evaluation=AllocatingEvaluation specify whether A and b are implemented allocating or in-place
• stopping_criterion::StoppingCriterion=StopAfterIteration(manifold_dimension(TpM))|StopWhenRelativeResidualLess(c,1e-8), where c is the norm of $\lVert b \rVert$.

Output

the obtained (approximate) minimizer $X^*$. To obtain the whole final state of the solver, see get_solver_return for details.

ConjugateResidualState{T,R,TStop<:StoppingCriterion} <: AbstractManoptSolverState

A state for the conjugate_residual solver.

Fields

• X::T: the iterate
• r::T: the residual $r = -b(p) - \mathcal A(p)[X]$
• d::T: the conjugate direction
• Ar::T, Ad::T: storages for $\mathcal A$
• rAr::R: internal field for storing $⟨ r, \mathcal A(p)[r] ⟩$
• α::R: a step length
• β::R: the conjugate coefficient
• stop::TStop: a StoppingCriterion for the solver

Constructor

    function ConjugateResidualState(
        TpM::TangentSpace,
        slso::SymmetricLinearSystemObjective;
        X=rand(TpM),
        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,
        stopping_criterion=StopAfterIteration(manifold_dimension(TpM)) |
                           StopWhenGradientNormLess(1e-8),
        kwargs...,
)

Initialise the state with default values.

SymmetricLinearSystemObjective{E<:AbstractEvaluationType,TA,T} <: AbstractManifoldObjective{E}

Model the objective

\[f(X) = \frac{1}{2} \lVert \mathcal A[X] + b \rVert_p^2,\qquad X ∈ T_p\mathcal M,\]

defined on the tangent space $T_p\mathcal M$ at $p$ on the manifold $\mathcal M$.

In other words this is an objective to solve $\mathcal A(p)[X] = -b(p)$ for some linear symmetric operator and a vector function. Note the minus on the right hand side, which makes this objective especially tailored for (iteratively) solving Newton-like equations.

Fields

• A!!: a symmetric, linear operator on the tangent space
• b!!: a gradient function

where A!! can work as an allocating operator (M, p, X) -> Y or an in-place one (M, Y, p, X) -> Y, and similarly b!! can either be a function (M, p) -> X or (M, X, p) -> X. The first variants allocate for the result, the second variants work in-place.

Constructor

SymmetricLinearSystemObjective(A, b; evaluation=AllocatingEvaluation())

Generate the objective specifying whether the two parts work allocating or in-place.

StopWhenRelativeResidualLess <: StoppingCriterion

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

\[ \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 indicates at which iteration (including i=0) the stopping criterion was fulfilled and is -1 while it is not fulfilled.
• c: the initial norm
• ε: the threshold
• norm_rk: the last computed norm of the residual

Constructor

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

Initialise the stopping criterion.