GitHub

Generalized Cauchy direction subsolver

The generalized Cauchy direction (GCD) subsolver 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 subsolver 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 subsolver. 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 $X$.

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

Simple stepsize limiting

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

Manopt.MaxStepsizeInDirectionSubsolver — Type
struct MaxStepsizeInDirectionSubsolver end

Helper container for finding the maximum stepsize in a direction. Stores the manifold M, container for the list of bounds F_list, and the bound indices.

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.

source
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 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.
source

Internal types and method

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.

source
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.

source
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.
source
Manopt.GeneralizedCauchyDirectionSubsolver — Type
GeneralizedCauchyDirectionSubsolver{TM <: AbstractManifold, TP, T_HA <: AbstractQuasiNewtonDirectionUpdate, TFU <: AbstractSegmentHessianUpdater}

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

source

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^{ op} 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.

source
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.

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

Evaluate the quadratic form associated with the stored Hessian approximation.

source
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} B 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 $B$ is ha.A.

source
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^{ op} B c$ where $c$ are the coordinates of the tangent vector X at p (in the basis d.basis) and $B$ is d.matrix.

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

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

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

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

source
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.

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

Evaluate the quadratic form associated with the stored Hessian approximation.

source
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.

source
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.

source

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.

source
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.

source
Manopt.GenericSegmentHessianUpdater — Type
struct GenericSegmentHessianUpdater <: AbstractSegmentHessianUpdater end

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

source
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.

source
get_bounds_index(::AbstractManifold)

Get the bound indices of manifold M. Standard manifolds don't have bounds, so Base.OneTo(1) is returned.

source
Manopt.get_stepsize_bound — Function
get_stepsize_bound(M::Hyperrectangle, x, 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:

  1. If d[i] > 0, the formula reads (M.ub[i] - p[i]) / d[i].
  2. If d[i] < 0, the formula reads (M.lb[i] - p[i]) / d[i].
  3. If d[i] == 0, the result is Inf.
source
get_stepsize_bound(M::AbstractManifold, x, d, i)

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

source
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.

source
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.

source
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.

source
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.

source
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₂ - 2·cs₁ᵀ M₂₁ cy₂ - cs₁ᵀ M₂₂ cs₂$ using the blocks $M₁₁$, $M₂₁$, $M₂₂$ stored in gh and the current scale $θ$. Returns the scalar value.

source
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.

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