Fix some issues with CG and certain coefficient rules when using restart (#604)

* refactor CGRule
* Fix some issues with CG + HZ
* provide a proper test + formatting
* fix CG test
* Refactor CG rules.
* Fix 2 cases where we implicitly relied on storage being updated in case it was not filled and now fill it explicitly.
* bump version.

---------

Co-authored-by: Ronny Bergmann <git@ronnybergmann.net>

Author: mateuszbaran

Changelog.md

@@ -6,6 +6,14 @@ 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).
 
+## [0.5.37] May 5, 2026
+
+### Changed
+
+* The default restart rule for `conjugate_gradient_descent` is now `RestartOnNonDescent` instead of `NeverRestart`, which makes the algorithm more robust to non-convexity and numerical issues. The old default can still be used by explicitly passing `restart_condition=NeverRestart()`. (#604)
+* `HagerZhangCoefficientRule` now has a safeguard against the denominator being too close to zero (the `denom_threshold` field). By default it is set to 1.0e-10. You can set it to a lower positive value (or even zero) to weaken the safeguard, but it is recommended to keep it to avoid numerical issues. (#604)
+* introduce for all `Rule`s also a variant without being encapsulated in a memory, where the old values have to be passed as keywords. This is now used by the `ConjugateGradientBealeRestartRule` when evaluating its inner rule. (#604)
+
 ## [0.5.36] April 24, 2026
 
 ### Added

Project.toml

@@ -1,6 +1,6 @@
 name = "Manopt"
 uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
-version = "0.5.36"
+version = "0.5.37"
 authors = [{family-names = "Bergmann", given-names = "Ronny", alias = "kellertuer", city = "Trondheim", affiliation = "Norwegian University of Science and Technology", country = "NO", email = "manopt@ronnybergmann.net", orcid = "https://orcid.org/0000-0001-8342-7218", website = "https://ronnybergmann.net"}]
 
 [workspace]

src/plans/conjugate_gradient_plan.jl

@@ -500,10 +500,11 @@ end
     get_storage(a::AbstractStateAction, key::Symbol)
 
 Return the internal value of the [`AbstractStateAction`](@ref) `a` at the
-`Symbol` `key`.
+`Symbol` `key`. Returns `nothing` if the key does not exist
 """
-get_storage(a::AbstractStateAction, key::Symbol) = a.values[key]
-
+function get_storage(a::AbstractStateAction, key::Symbol)
+    return get(a.values, key, nothing)
+end
 """
     get_storage(a::AbstractStateAction, ::PointStorageKey{key}) where {key}
 

src/plans/solver_state.jl

@@ -142,7 +142,7 @@ function conjugate_gradient_descent!(
         mgo::O,
         p;
         coefficient::Union{DirectionUpdateRule, ManifoldDefaultsFactory} = ConjugateDescentCoefficient(),
-        restart_condition::AbstractRestartCondition = NeverRestart(),
+        restart_condition::AbstractRestartCondition = RestartOnNonDescent(),
         retraction_method::AbstractRetractionMethod = default_retraction_method(M, typeof(p)),
         stepsize::Union{Stepsize, ManifoldDefaultsFactory} = default_stepsize(
             M, ConjugateGradientDescentState; retraction_method = retraction_method
@@ -195,7 +195,8 @@ function step_solver!(amp::AbstractManoptProblem, cgs::ConjugateGradientDescentS
     cgs.δ .-= cgs.X
     if (cgs.restart_condition(amp, cgs, k))
         # restart solver; set dir to -grad
-        cgs.δ = -copy(get_manifold(amp), cgs.p, cgs.X)
+        copyto!(M, cgs.δ, cgs.X)
+        cgs.δ .*= -1
         update_storage!(cgs.coefficient.storage, amp, cgs)
         cgs.β = 0.0
     end

src/solvers/conjugate_gradient_descent.jl

@@ -1,7 +1,8 @@
 using Manopt, Manifolds, Test
 
 struct DummyCGCoeff <: DirectionUpdateRule end
-(u::DummyCGCoeff)(p, s, k) = 0.2
+(::DummyCGCoeff)(pr, st, k; kwagrs...) = 0.2
+(::Manopt.DirectionUpdateRuleStorage{DummyCGCoeff})(pr, st, k) = 0.2
 Manopt.update_rule_storage_points(::DummyCGCoeff) = Tuple{}
 Manopt.update_rule_storage_vectors(::DummyCGCoeff) = Tuple{}
 
@@ -16,22 +17,26 @@ Manopt.update_rule_storage_vectors(::DummyCGCoeff) = Tuple{}
         p0 = [1.0, 0.0]
         pr = DefaultManoptProblem(M, ManifoldGradientObjective(f, grad_f))
         cgs2 = ConjugateGradientDescentState(
-            M;
-            p = p0,
+            M; p = p0,
             stopping_criterion = StopAfterIteration(2),
             stepsize = Manopt.ConstantStepsize(M, 1.0),
             coefficient = dur2,
         )
         cgs2.X = [0.0, 0.2]
+        # Fake update history to get a certain old X and old p
+        cgs2.coefficient(pr, cgs2, 0)
+        # the inner check is 0.2 which is still less than 0.3
         @test cgs2.coefficient(pr, cgs2, 1) != 0
         cgs3 = ConjugateGradientDescentState(
-            M;
-            p = p0,
+            M; p = p0,
             stopping_criterion = StopAfterIteration(2),
             stepsize = Manopt.ConstantStepsize(M, 1.0),
             coefficient = dur3,
         )
         cgs3.X = [0.0, 0.2]
+        # Fake update history to get a certain old X and old p
+        cgs3.coefficient(pr, cgs3, 0)
+        # then we are above the threshold 0.1 (namely at 0.2) and we get a descent step
         @test cgs3.coefficient(pr, cgs3, 1) == 0
     end
     @testset "representation and summary of Coefficients" begin

test/plans/test_conjugate_gradient_plan.jl

@@ -404,4 +404,24 @@ using ManifoldDiff: grad_distance
         grad_f(M, p) = sum(1 / n * grad_distance.(Ref(M), data, Ref(p)))
         @test conjugate_gradient_descent(M, f, grad_f, data[1]) isa PoincareBallPoint
     end
+
+    @testset "Issue #603: CG with HZ rule on a numerically challenging problem" begin
+        M = Sphere(2)
+        p0 = [1.0, 0.0, 0.0]
+
+        a = [0.0, 1.0, 0.0]
+
+        f(M, p) = 0.5 * norm(p - a)^2
+        grad_f(M, p) = project(M, p, p - a)
+
+        cgs = conjugate_gradient_descent(
+            M,
+            f,
+            grad_f,
+            p0;
+            coefficient = ConjugateGradientBealeRestart(HagerZhangCoefficient()),
+            return_state = true,
+        )
+        @test norm(M, cgs.p, grad_f(M, cgs.p)) < 1.0e-8
+    end
 end