Generalized Cauchy direction and box constraints

.github/workflows/ci.yml

@@ -18,7 +18,6 @@ jobs:
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.julia-version }}
-          arch: x64
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1

.github/workflows/nightly.yml

@@ -15,7 +15,6 @@ jobs:
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.julia-version }}
-          arch: x64
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         env:

Changelog.md

@@ -6,6 +6,19 @@ The file was started with Version `0.4`.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## Unreleased
+
+### Added
+
+* `nonpositive_curvature_behavior` for `QuasiNewtonLimitedMemoryDirectionUpdate` that determines how transported (y, s) vector pairs are treated after transport; if their inner product gets too low, it may lead to non-positive-definite Hessians which needs to be avoided. This resolves issue #549. (#554)
+* `GeneralizedCauchyDirectionFinder` for handling direction selection in the presence of box (`Hyperrectangle`) constraints in quasi-Newton methods. This allows for L-BFGS-B-style box constraint handling. (#554)
+* New stopping criteria: `StopWhenRelativeAPosterioriCostChangeLessOrEqual` and `StopWhenProjectedNegativeGradientNormLess`. (#554)
+
+### Fixed
+
+* Line searches consistently respect `stop_when_stepsize_exceeds` keyword argument as a hard limit. (#554)
+* `StopWhenChangeLess` falsely claimed to indicate convergence. This is now fixed.
+
 ## [0.5.32] January 15, 2026
 
 ### Fixed
@@ -46,7 +59,7 @@ Moved the documentation glossaries to using the new [Glossaries.jl](https://gith
 * Removed `atol` from `DebugFeasibility` and instead use the one newly added `atol` from the `ConstrainedManifoldObjective`. (#546)
 * Move from CompatHelper to dependabot to keep track of dependency updates in Julia packages. (#547)
 * moved the `ManoptTestSuite` module to a sub module `Manopt.Test` within `Manopt.jl`,
-so it can be easier reused by others as well (#550)
+  so it can be easier reused by others as well (#550)
 * moved to using a `Project.toml` for tests and an overall `[Workspace]`.
   This also allows finally to run single test files without installing all packages manually, but instead just switching to and instantiating the test environment. (#550)
 * for compatibility, state also `[source]` entries consistently in the sub `Project.toml` files. (#550)

docs/Project.toml

@@ -23,8 +23,8 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 RipQP = "1e40b3f8-35eb-4cd8-8edd-3e515bb9de08"
 
-[sources.Manopt]
-path = ".."
+[sources]
+Manopt = {path = ".."}
 
 [compat]
 BenchmarkTools = "1.3"

docs/make.jl

@@ -45,6 +45,7 @@ tutorials_menu =
     "Implement a solver" => "tutorials/ImplementASolver.md",
     "Optimize on your own manifold" => "tutorials/ImplementOwnManifold.md",
     "Do constrained optimization" => "tutorials/ConstrainedOptimization.md",
+    "Do optimization with bounds" => "tutorials/BoxDomain.md",
 ]
 # Check whether all tutorials are rendered, issue a warning if not (and quarto if not set)
 all_tutorials_exist = true
@@ -189,6 +190,7 @@ makedocs(;
             "Douglas—Rachford" => "solvers/DouglasRachford.md",
             "Exact Penalty Method" => "solvers/exact_penalty_method.md",
             "Frank-Wolfe" => "solvers/FrankWolfe.md",
+            "Generalized Cauchy direction subsolver" => "solvers/generalized_cauchy_direction_subsolver.md",
             "Gradient Descent" => "solvers/gradient_descent.md",
             "Interior Point Newton" => "solvers/interior_point_Newton.md",
             "Levenberg–Marquardt" => "solvers/LevenbergMarquardt.md",

docs/src/extensions.md

@@ -45,10 +45,12 @@ x_opt = quasi_Newton(
 )
 ```
 
-In general this defines the following new [stepsize](@ref Stepsize)
+In general this defines the following new [stepsize](@ref Stepsize) with helper functions for setting and getting the maximum step size:
 
 ```@docs
 Manopt.LineSearchesStepsize
+Manopt.linesearches_get_max_alpha
+Manopt.linesearches_set_max_alpha
 ```
 
 ## Manifolds.jl

docs/src/references.bib

@@ -260,6 +260,33 @@ @book{Boumal:2023
     ISBN      = {978-1-00-916616-4}
 }
 
+@article{ByrdNocedalSchnabel:1994,
+	title = {Representations of quasi-{Newton} matrices and their use in limited memory methods},
+	volume = {63},
+	issn = {1436-4646},
+	doi = {10.1007/BF01582063},
+	number = {1},
+	journal = {Mathematical Programming},
+	author = {Byrd, Richard H. and Nocedal, Jorge and Schnabel, Robert B.},
+	month = jan,
+	year = {1994},
+	pages = {129--156},
+}
+
+@article{ByrdLuNocedalZhu:1995,
+	title = {A {Limited} {Memory} {Algorithm} for {Bound} {Constrained} {Optimization}},
+	volume = {16},
+	issn = {1064-8275},
+	doi = {10.1137/0916069},
+	number = {5},
+	journal = {SIAM Journal on Scientific Computing},
+	author = {Byrd, Richard H. and Lu, Peihuang and Nocedal, Jorge and Zhu, Ciyou},
+	month = sep,
+	year = {1995},
+	note = {Publisher: Society for Industrial and Applied Mathematics},
+	pages = {1190--1208},
+}
+
 % --- C
 %
 %
@@ -871,4 +898,17 @@ @article{ZhangSra:2018
     TITLE      = {Towards Riemannian accelerated gradient methods},
     URL        = {https://arxiv.org/abs/1806.02812},
     YEAR       = {2018},
-}
+}
+
+@article{ZhuByrdLuNocedal:1997,
+	title = {Algorithm 778: {L}-{BFGS}-{B}: {Fortran} subroutines for large-scale bound-constrained optimization},
+	volume = {23},
+	issn = {0098-3500},
+	doi = {10.1145/279232.279236},
+	number = {4},
+	journal = {ACM Trans. Math. Softw.},
+	author = {Zhu, Ciyou and Byrd, Richard H. and Lu, Peihuang and Nocedal, Jorge},
+	month = dec,
+	year = {1997},
+	pages = {550--560},
+}

docs/src/solvers/generalized_cauchy_direction_subsolver.md

@@ -0,0 +1,63 @@
+# Generalized Cauchy Direction subsolver
+
+The Generalized Cauchy Direction (GCD) subsolver is a component in optimization algorithms that handle problems with bound constraints. It solves the following problem
+
+```math
+\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`](@extref Manifolds.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 \leq s_{\max} \leq \infty$ 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.
+
+## Internal types and method
+
+### Symbols related to the GCD computation
+
+These symbols are directly used by solvers to compute the descent direction corresponding to the Generalized Cauchy direction.
+
+```@docs
+Manopt.has_anisotropic_max_stepsize
+Manopt.find_generalized_cauchy_direction!
+Manopt.GeneralizedCauchyDirectionFinder
+```
+
+### Symbols related to the Hessian approximation
+
+These symbols are used to evaluate the Hessian approximation at specific tangent vectors during the generalized Cauchy direction computation.
+
+```@docs
+Manopt.hessian_value
+Manopt.hessian_value_diag
+```
+
+### Symbols related to bound handling
+
+These are internal symbols used to manage and manipulate bound constraints during the GCP computation.
+
+```@docs
+Manopt.init_updater!
+Manopt.UnitVector
+Manopt.to_coordinate_index
+Manopt.AbstractSegmentHessianUpdater
+Manopt.GenericSegmentHessianUpdater
+Manopt.get_bounds_index
+Manopt.get_stepsize_bound
+Manopt.get_at_bound_index
+Manopt.set_stepsize_bound!
+Manopt.set_zero_at_index!
+```
+
+### Symbols related to specific Hessian approximations
+
+```@docs
+Manopt.LimitedMemorySegmentHessianUpdater
+Manopt.hessian_value_from_inner_products
+Manopt.set_M_current_scale!
+```

docs/src/solvers/quasi_Newton.md

@@ -88,6 +88,7 @@ QuasiNewtonLimitedMemoryDirectionUpdate
 QuasiNewtonCautiousDirectionUpdate
 Manopt.initialize_update!
 QuasiNewtonPreconditioner
+QuasiNewtonLimitedMemoryBoxDirectionUpdate
 ```
 
 ## Hessian update rules

ext/ManoptLineSearchesExt.jl

@@ -5,6 +5,37 @@ import Manopt: LineSearchesStepsize
 using ManifoldsBase
 using LineSearches
 
+Manopt.linesearches_get_max_alpha(ls::LineSearches.HagerZhang) = ls.alphamax
+Manopt.linesearches_get_max_alpha(ls::LineSearches.MoreThuente) = ls.alphamax
+
+function Manopt.linesearches_set_max_alpha(ls::LineSearches.HagerZhang{T, Tm}, max_alpha::T) where {T, Tm}
+    return HagerZhang{T, Tm}(
+        delta = ls.delta,
+        sigma = ls.sigma,
+        alphamax = max_alpha,
+        rho = ls.rho,
+        epsilon = ls.epsilon,
+        gamma = ls.gamma,
+        linesearchmax = ls.linesearchmax,
+        psi3 = ls.psi3,
+        display = ls.display,
+        mayterminate = ls.mayterminate,
+        cache = ls.cache,
+        check_flatness = ls.check_flatness,
+    )
+end
+function Manopt.linesearches_set_max_alpha(ls::LineSearches.MoreThuente{T}, max_alpha::T) where {T}
+    return MoreThuente{T}(
+        f_tol = ls.f_tol,
+        gtol = ls.gtol,
+        x_tol = ls.x_tol,
+        alphamin = ls.alphamin,
+        alphamax = max_alpha,
+        maxfev = ls.maxfev,
+        cache = ls.cache
+    )
+end
+
 function (cs::Manopt.LineSearchesStepsize)(
         mp::AbstractManoptProblem,
         s::AbstractManoptSolverState,
@@ -24,6 +55,19 @@ function (cs::Manopt.LineSearchesStepsize)(
     # guess initial alpha
     α0 = cs.initial_guess(mp, s, k, cs.last_stepsize, η; lf0 = fp, Dlf0 = dphi_0)
 
+    # handle stepsize limit
+    local ls # COV_EXCL_LINE
+    if :stop_when_stepsize_exceeds in keys(kwargs)
+        new_max_alpha = min(
+            kwargs[:stop_when_stepsize_exceeds],
+            Manopt.linesearches_get_max_alpha(cs.linesearch),
+        )
+        ls = Manopt.linesearches_set_max_alpha(cs.linesearch, new_max_alpha)
+        α0 = min(α0, new_max_alpha)
+    else
+        ls = cs.linesearch
+    end
+
     # perform actual line-search
 
     function ϕ(α)
@@ -41,7 +85,7 @@ function (cs::Manopt.LineSearchesStepsize)(
         return Manopt.get_cost_and_differential(mp, p_tmp, Y_tmp; Y = X_tmp)
     end
 
-    α, fp = cs.linesearch(ϕ, dϕ, ϕdϕ, α0, fp, dphi_0)
+    α, fp = ls(ϕ, dϕ, ϕdϕ, α0, fp, dphi_0)
     cs.last_stepsize = α
     return α
 end

ext/ManoptManifoldsExt/ManoptManifoldsExt.jl

@@ -2,7 +2,7 @@ module ManoptManifoldsExt
 
 using ManifoldsBase: exp, log, ParallelTransport, vector_transport_to
 using Manopt
-using Manopt: _math, _tex, ManifoldDefaultsFactory, _produce_type
+using Manopt: _math, _tex, ManifoldDefaultsFactory, _produce_type, get_stepsize_bound
 import Manopt:
     max_stepsize,
     get_gradient,