Generalized Cauchy direction sub solver

The generalized Cauchy direction (GCD) sub solver is a component in optimization algorithms that handle problems with bound constraints [BBP26]. It solves the following problem

\[\begin{align*} \operatorname*{arg\,min}_{Y ∈ T_p (D \times \mathcal{M})}&\ m_p(Y), \qquad m_p(Y) = ⟨X_g, Y⟩_p + \frac{1}{2} ⟨\mathcal{H}_p[Y], Y⟩_p\\ \text{such that}& \ \exp_p(Y) = \exp_p(\alpha X) \in D \times \mathcal{M} \text{ for some } \alpha \in [0, A] \end{align*}\]

where $X=(X_{\mathrm{D}}, X_{\mathcal{M}})$ is a given direction, the exponential map handles projection of the tangent vector when reaching the boundary, $D$ is a box domain (Hyperrectangle), $\mathcal{M}$ is a Riemannian manifold, $X_g$ is the gradient of a scalar function $f$ at point $p=(p_{\mathrm{D}}, p_{\mathcal{M}})$, $A$ is the maximum allowed step size on $\mathcal{M}$ at point $p=(p_{\mathrm{D}}, p_{\mathcal{M}})$ in direction $X_{\mathcal{M}}$ (infinity is supported) and $\mathcal{H}_p$ is a linear operator that approximates the Hessian of $f$ at $p$.

Additionally, the sub solver indicates whether the selected direction $Y$ reaches the boundary of $D$ at some point, in which case the subsequent step size selection in direction $Y$ needs to be limited to the interval $[0, s_{\max}]$, where the number $1 ≤ s_{\max} ≤ ∞$ is also returned by the sub solver. Note that the value $s_{\max}=1$ is obtained when the minimum lies at the boundary of $D$, while larger values indicate that we are further away from the boundary along the selected direction $Y$.

The solver is currently primarily intended for internal use by optimization algorithms that require bound-constrained sub problem solutions.

Simple step size limiting

In case there is no Hessian approximation available, a simple step size limiting procedure can be used to limit the step size in direction $X$ to the maximum allowed by the boundary of $D$ and the maximum allowed step size on $\mathcal{M}$. This procedure is available using the following:

Manopt.MaxStepsizeInDirectionSubsolver — Type

MaxStepsizeInDirectionSubsolver{TFT, TBI}

Helper container for finding the maximum stepsize in a direction.

Fields

F_list::Vector{TFT}: the bounds met along the direction, as (stepsize, index) pairs
bounds_indices::TBI: the indices of the components that are bounded

Constructor

MaxStepsizeInDirectionSubsolver(M::AbstractManifold, p)

Initialize the MaxStepsizeInDirectionSubsolver for manifold M and point p. The F_list is initialized to be empty and will be populated during the search for the maximum stepsize in a direction. Floating point type of the elements bounds in F_list is determined by the number type of p.

The MaxStepsizeInDirectionSubsolver can be reused for multiple different points and directions on the same manifold, but it is not thread-safe.

Manopt.find_max_stepsize_in_direction — Function

find_max_stepsize_in_direction(M::AbstractManifold, gcd::MaxStepsizeInDirectionSubsolver, p, d)

Find the maximum stepsize that can be performed from point p in direction d.

The function returns a pair (status, maxstepsize) where status is a symbol describing the result of the search, and `maxstepsizeis the maximum stepsize that can be taken frompin the directiond`.

The status can be one of the following:

:found_limited if a finite bound exists; max_stepsize is then the smallest positive stepsize bound,
:found_unlimited if no finite bound exists; max_stepsize is then Inf,
:not_found if no positive stepsize bound exists at all; max_stepsize is then NaN.

Internal types and methods

These symbols are directly used by solvers to compute the descent direction corresponding to the Generalized Cauchy direction.

Manopt.has_anisotropic_max_stepsize — Function

Manopt.has_anisotropic_max_stepsize(::Hyperrectangle)

Returns true, as Hyperrectangle manifold requires generalized Cauchy point computation in solvers.

has_anisotropic_max_stepsize(M::AbstractManifold)

Return true if M has a max_stepsize that depends on the direction, for example if M is a Hyperrectangle-like manifold with corners, or a product of it with a standard manifold. Otherwise return false.

Manopt.find_generalized_cauchy_direction! — Function

find_generalized_cauchy_direction!(
    M::AbstractManifold,
    gcd::GeneralizedCauchyDirectionSubsolver, d_out, p, d, X
)

Find generalized Cauchy direction looking from point p on manifold M in direction d and save it to d_out. Gradient of the objective at p is X.

The function returns a pair (status, maxstepsize) where status is a symbol describing the result of the search, and `maxstepsizeis the maximum stepsize that can be taken in the directiond_out`.

The status can be one of the following:

:found_limited if the point was found and we can perform a step of length at most 1 in direction d_out afterwards,
:found_unlimited if the point was found and we can perform a step of length at most max_stepsize(M, p) in direction d_out afterwards,
:not_found if the search cannot be performed in direction d.

Manopt.GeneralizedCauchyDirectionSubsolver — Type

GeneralizedCauchyDirectionSubsolver{TX, T_HA, TFU <: AbstractSegmentHessianUpdater, TFT, TBI, TO}

Helper container for generalized Cauchy direction search. Stores the cached original descent direction (d_original), the quasi-Newton direction update ha, the hessian_segment_updater, which computes certain values of the Hessian while advancing segments, the list of bounds F_list, the bound indices bounds_indices, and the heap ordering. Instances are reused across segments during find_generalized_cauchy_direction! to avoid allocations.

These symbols are used to evaluate the Hessian approximation at specific tangent vectors during the generalized Cauchy direction computation.

Manopt.hessian_value — Function

hessian_value(d::QuasiNewtonMatrixDirectionUpdate, M, p, X::UnitVector, Y)

Evaluate the quadratic form associated with the stored quasi-Newton matrix. Returns the scalar $c_b^{\mathrm{T}} B c$ where $c_b$ are the coordinates of the UnitVector X at p (assumed to correspond to the basis d.basis), $c$ are the coordinates of the tangent vector Y at p (in the basis d.basis) and $B$ is d.matrix.

hessian_value(gh::QuasiNewtonLimitedMemoryBoxDirectionUpdate, M::AbstractManifold, p, X::UnitVector, Y)

Compute $⟨X, B Y⟩$, where $B$ is the (1, 1)-Hessian represented by gh, where X is the UnitVector.

hessian_value(ha::LevenbergMarquardtBoxSubsolver, M, p, X::UnitVector, Y)

Evaluate the quadratic form associated with the stored Hessian approximation.

hessian_value(ha::CoordinatesNormalSystemState, M, p, X::UnitVector, Y)

Evaluate the quadratic form associated with the stored Hessian approximation. Returns the scalar $c_b^{\top} A c$ where $c_b$ are the coordinates of the UnitVector X at p (assumed to correspond to the basis ha.basis), $c$ are the coordinates of the tangent vector Y at p (in the basis ha.basis) and $A$ is ha.A.

Manopt.hessian_value_diag — Function

hessian_value_diag(d::QuasiNewtonMatrixDirectionUpdate, M, p, X)

Evaluate the quadratic form associated with the stored quasi-Newton matrix. Returns the scalar $c^{\mathrm{T}} B c$ where $c$ are the coordinates of the tangent vector X at p (in the basis d.basis) and $B$ is d.matrix.

hessian_value_diag(d::QuasiNewtonMatrixDirectionUpdate, M, p, X::UnitVector)

Evaluate the quadratic form associated with the stored quasi-Newton matrix. Returns the scalar $c^{\mathrm{T}} B c$ where $c$ are the coordinates of the UnitVector X at p (in the basis d.basis) and $B$ is d.matrix.

hessian_value_diag(gh::QuasiNewtonLimitedMemoryBoxDirectionUpdate, M::AbstractManifold, p, X)

Compute $⟨X, B X⟩$, where $B$ is the (1, 1)-Hessian represented by gh.

hessian_value_diag(gh::QuasiNewtonLimitedMemoryBoxDirectionUpdate, M::AbstractManifold, p, X::UnitVector)

Compute $⟨X, B X⟩$, where $B$ is the (1, 1)-Hessian represented by gh, and X is the UnitVector.

hessian_value_diag(ha::LevenbergMarquardtBoxSubsolver, M, p, X)

Evaluate the quadratic form associated with the stored Hessian approximation.

hessian_value_diag(ha::CoordinatesNormalSystemState, M::AbstractManifold, p, X::UnitVector)

Evaluate the quadratic form associated with the Hessian approximation of a CoordinatesNormalSystemState. Returns the scalar $c^{\top} A c$ where $c$ are the coordinates of the UnitVector X at p (in the basis ha.basis) and $A$ is ha.A.

hessian_value_diag(ha::CoordinatesNormalSystemState, M::AbstractManifold, p, X)

Evaluate the quadratic form associated with the Hessian approximation of a CoordinatesNormalSystemState. Returns the scalar $c^{\top} A c$ where $c$ are the coordinates of X at p (in the basis ha.basis) and $A$ is ha.A.

These are internal symbols used to manage and manipulate bound constraints during the GCD computation.

Manopt.init_updater! — Function

init_updater!(::AbstractManifold, hessian_segment_updater::AbstractSegmentHessianUpdater, p, d, ha)

Method for initialization of AbstractSegmentHessianUpdater hessian_segment_updater just before the loop that examines subsequent intervals for GCD.

Manopt.UnitVector — Type

UnitVector{TB}

A type representing a unit tangent vector on a Hyperrectangle-like manifold with corners, or a product of it with a standard manifold. The field index stores the index of the element equal to 1. All other elements are equal to 0.

Manopt.to_coordinate_index — Function

to_coordinate_index(M::AbstractManifold, b::UnitVector{Int}, B::AbstractBasis)

Get the index of coordinate equal to 1 of UnitVector b with respect to AbstractBasis B.

to_coordinate_index(M::ProductManifold, b::UnitVector{Tuple{Int, Int}}, B::AbstractBasis)

Get the index of coordinate equal to 1 of UnitVector b with respect to AbstractBasis B.

Manopt.AbstractSegmentHessianUpdater — Type

abstract type AbstractSegmentHessianUpdater end

Abstract type for methods that calculate f' and f'' in the GCD calculation in subsequent line segments in GeneralizedCauchyDirectionSubsolver.

Manopt.GenericSegmentHessianUpdater — Type

GenericSegmentHessianUpdater{TX} <: AbstractSegmentHessianUpdater

Generic f' and f'' calculation that only relies on hessian_value but is relatively slow for high-dimensional domains.

Fields

d_z::TX: the accumulated direction of the segments walked so far
d_tmp::TX: the current direction, with the components that reached a bound set to zero

Manopt.get_bounds_index — Function

get_bounds_index(::Hyperrectangle)

Get the bound indices of Hyperrectangle M. They are the same as the indices of the lower (or upper) bounds.

get_bounds_index(::AbstractManifold)

Get the bound indices of manifold M. Standard manifolds don't have bounds, so Base.OneTo(0), that is an empty range, is returned.

Manopt.get_stepsize_bound — Function

get_stepsize_bound(M::Hyperrectangle, p, d, i)

Get the upper bound on moving in direction d from point p on Hyperrectangle M, for the bound index i. There are three cases:

If d[i] > 0, the formula reads (M.ub[i] - p[i]) / d[i].
If d[i] < 0, the formula reads (M.lb[i] - p[i]) / d[i].
If d[i] == 0, the result is Inf.

get_stepsize_bound(M::AbstractManifold, p, d, i)

Get the upper bound on moving in direction d from point p on manifold M, for the bound index i.

Manopt.get_at_bound_index — Function

Manopt.get_at_bound_index(::Hyperrectangle, X, b)

Extract the element of tangent vector X to a point on Hyperrectangle at index b.

Manopt.set_stepsize_bound! — Function

Manopt.set_stepsize_bound!(M::Hyperrectangle, d_out, p, d, t_current::Real)

Limit the per-component stepsize in d_out to the bound imposed by the box constraints.

For each element i in the tangent vector d_out, if the stepsize bound in direction d for that element is less than t_current, set the element of d_out to the distance from p[i] to the bound in the direction of d[i]. If the stepsize bound is non-positive, set the element to 0.

set_stepsize_bound!(M::AbstractManifold, d_out, p, d, t_current::Real)

Limit the per-component stepsize in d_out to the bound imposed by the box constraints.

For each component at index i in the tangent vector d_out, if the stepsize bound in direction d for that component is less than t_current, set the element of d_out to the distance from p[i] to the bound in the direction of d[i]. If the stepsize bound is non-positive, set the element to 0.

By default it does not modify d_out because most manifolds don't have direction-specific stepsize bounds, and general anisotropic bounds are handled differently.

Manopt.set_zero_at_index! — Function

Manopt.set_zero_at_index!(M::Hyperrectangle, d, i)

Set element of tangent vector d on Hyperrectangle at index i to 0.

set_zero_at_index!(M::ProductManifold, d, i::Tuple{Int,Any})

Set the element of the i[1]th component of d at bound index i[2] to zero.

Manopt.LimitedMemorySegmentHessianUpdater — Type

LimitedMemorySegmentHessianUpdater{TV <: AbstractVector} <: AbstractSegmentHessianUpdater

Hessian value calculation for generalized Cauchy direction line segments that is optimized for QuasiNewtonLimitedMemoryBoxDirectionUpdate. It relies on a specific Hessian structure.

Fields

p_s::TV: the inner products of the stored $s_k$ with the current direction
p_y::TV: the inner products of the stored $y_k$ with the current direction
c_s::TV: the inner products of the stored $s_k$ with the accumulated direction $d_z$
c_y::TV: the inner products of the stored $y_k$ with the accumulated direction $d_z$

Manopt.hessian_value_from_inner_products — Function

hessian_value_from_inner_products(gh::QuasiNewtonLimitedMemoryBoxDirectionUpdate, iss::Real, cy1, cs1, cy2, cs2)

Evaluate the quadratic form defined by the current blockwise Hessian approximation stored in gh, given precomputed coordinate vectors.

Arguments:

iss: inner product of original vectors.
cy1, cy2: coordinates of $y$-like vectors in the $Y_k$ basis.
cs1, cs2: coordinates of $s$-like vectors in the scaled $S_k$ basis.

The result is $θ·iss - cy₁ᵀ M₁₁ cy₂ - cs₁ᵀ M₂₁ cy₂ - cs₂ᵀ M₂₁ cy₁ - cs₁ᵀ M₂₂ cs₂$ using the blocks $M₁₁$, $M₂₁$, $M₂₂$ stored in gh and the current scale $θ$. Returns the scalar value.

Manopt.update_current_scale! — Function

update_current_scale!(M::AbstractManifold, p, gh::QuasiNewtonLimitedMemoryBoxDirectionUpdate)

Refresh the scaling factor and blockwise Hessian approximation stored in gh using the nonzero curvature pairs currently in memory.

Identifies the most recent index with nonzero $ρ_i$ to scale the initial Hessian guess by $ρ_i‖y_i‖^2 / θ$.
Builds $L_k$ and $S_k^\top S_k$ from the stored $(s_i, y_i)$ pairs and updates the block matrices $M₁₁$, $M₂₁$, and $M₂₂$ via the blockwise inverse formula.
If all $ρ_i$ vanish, resets current_scale to the inverse of initial_scale and clears the block matrices.

Returns the mutated gh.

Literature

[BBP26]: M. Baran, R. Bergmann and P. Przybysz. A Riemannian quasi-Newton algorithm for optimization with Euclidean bounds (2026), arXiv:2605.10573 [math.OC].