Fix OptimizationManopt reporting MaxIters even when converged

docs/src/optimization_packages/manopt.md

@@ -52,11 +52,17 @@ opt = OptimizationManopt.GradientDescentOptimizer()
 optf = OptimizationFunction(rosenbrock, ADTypes.AutoZygote())
 
 prob = OptimizationProblem(
-    optf, x0, p; manifold = R2, stepsize = stepsize)
+    optf, x0, p; manifold = R2, stepsize = stepsize, maxiters = 25000)
 
 sol = Optimization.solve(prob, opt)
 ```
 
+!!! note
+
+    Plain gradient descent zig-zags slowly down Rosenbrock's narrow curved valley, so it
+    needs tens of thousands of iterations (still well under a second here) to bring the
+    gradient norm below the default tolerance and report `retcode = Success`.
+
 The box-constrained Karcher mean problem on the SPD manifold with the Frank-Wolfe algorithm can be solved as follows:
 
 ```@example Manopt
@@ -90,10 +96,12 @@ U = mean(data2)
 L = inv(sum(1 / N * inv(matrix) for matrix in data2))
 
 optf = OptimizationFunction(f, ADTypes.AutoZygote())
-prob = OptimizationProblem(optf, U; manifold = M, maxiters = 1000)
+prob = OptimizationProblem(optf, U; manifold = M, maxiters = 5000)
 
+opt = OptimizationManopt.FrankWolfeOptimizer()
 sol = Optimization.solve(
-    prob, opt, sub_problem = (M, q, p, X) -> closed_form_solution!(M, q, L, U, p, X))
+    prob, opt, sub_problem = (M, q, p, X) -> closed_form_solution!(M, q, L, U, p, X),
+    evaluation = Manopt.InplaceEvaluation())
 ```
 
 This example is based on the [example](https://juliamanifolds.github.io/ManoptExamples.jl/stable/examples/Riemannian-mean/) in the Manopt and [Weber and Sra'22](https://doi.org/10.1007/s10107-022-01840-5).
@@ -115,7 +123,7 @@ egrad(G, x, p = nothing) = (G .= -2 * A * x)
 
 optf = OptimizationFunction(cost, grad = egrad)
 x0 = rand(manifold)
-prob = OptimizationProblem(optf, x0, manifold = manifold)
+prob = OptimizationProblem(optf, x0, manifold = manifold, maxiters = 5000)
 
 sol = solve(prob, GradientDescentOptimizer())
 ```

lib/OptimizationManopt/src/OptimizationManopt.jl

@@ -44,6 +44,19 @@ function __map_optimizer_args!(
 
     if !isnothing(abstol)
         push!(criteria, _default_convergence_criterion(opt, manifold, abstol))
+    elseif !isnothing(maxiters) || !isnothing(maxtime)
+        # Without this, `criteria` above would contain *only* `StopAfterIteration`/
+        # `StopAfter`, which never `indicates_convergence` (Manopt.jl semantics), so
+        # `Manopt.has_converged` could never be true and the run would always report
+        # `ReturnCode.MaxIters`/`MaxTime` below, no matter how well it actually converged.
+        # Falling back to each solver's own default tolerance (rather than one constant)
+        # keeps this close to what Manopt would have done with no `stopping_criterion`
+        # override at all; optimizers whose real default isn't a single tolerance-based
+        # criterion opt out via `_default_fallback_abstol` returning `nothing`.
+        fallback_abstol = _default_fallback_abstol(opt)
+        if !isnothing(fallback_abstol)
+            push!(criteria, _default_convergence_criterion(opt, manifold, fallback_abstol))
+        end
     end
 
     if !isnothing(reltol)
@@ -284,6 +297,24 @@ function _default_convergence_criterion(::AbstractManoptOptimizer, M, abstol)
     return Manopt.StopWhenChangeLess(M, abstol)
 end
 
+# Fallback tolerance used, via `_default_convergence_criterion` above, only when the user
+# supplies `maxiters`/`maxtime` without an explicit `abstol` (see `__map_optimizer_args!`).
+# Each value mirrors that solver's own top-level default tolerance in Manopt.jl, so this
+# stays close to "what Manopt would have done with no override at all" instead of imposing
+# one generic number. Optimizers whose real Manopt default is a composite criterion that
+# isn't representable as a single `_default_convergence_criterion` call (`NelderMead`'s
+# `StopWhenPopulationConcentrated`, `CMAES`'s multi-part `default_cma_es_stopping_criterion`,
+# `ConvexBundle`'s `StopWhenLagrangeMultiplierLess`, not gradient-norm) opt out with
+# `nothing`, leaving their pre-existing behavior unchanged.
+_default_fallback_abstol(::GradientDescentOptimizer) = 1.0e-8
+_default_fallback_abstol(::ConjugateGradientDescentOptimizer) = 1.0e-8
+_default_fallback_abstol(::QuasiNewtonOptimizer) = 1.0e-6
+_default_fallback_abstol(::TrustRegionsOptimizer) = 1.0e-6
+_default_fallback_abstol(::AdaptiveRegularizationCubicOptimizer) = 1.0e-9
+_default_fallback_abstol(::FrankWolfeOptimizer) = 1.0e-6
+_default_fallback_abstol(::ParticleSwarmOptimizer) = 1.0e-4
+_default_fallback_abstol(::AbstractManoptOptimizer) = nothing
+
 function build_loss(f::OptimizationBase.OptimizationFunction, prob, cb)
     # TODO: I do not understand this. Why is the manifold not used?
     # Either this is an Euclidean cost, then we should probably still call `embed`,

lib/OptimizationManopt/test/core_tests.jl

@@ -51,6 +51,28 @@ R2 = Euclidean(2)
         @test sol.objective < 0.2
     end
 
+    @testset "Gradient descent maxiters without abstol still reports Success" begin
+        # Regression test: passing only `maxiters` (no `abstol`) used to fully replace
+        # Manopt's own default `stopping_criterion`, leaving `StopAfterIteration` as the
+        # *only* criterion. Since that criterion never `indicates_convergence`, the run
+        # was structurally guaranteed to report `ReturnCode.MaxIters`, even once it had
+        # actually converged. `maxiters` here is set well above what's needed so the run
+        # converges before hitting the cap.
+        x0 = zeros(2)
+        p = [1.0, 100.0]
+
+        stepsize = Manopt.ArmijoLinesearch(R2)
+        opt = OptimizationManopt.GradientDescentOptimizer()
+
+        optprob = OptimizationFunction(rosenbrock, OptimizationBase.AutoForwardDiff())
+        prob = OptimizationProblem(
+            optprob, x0, p; manifold = R2, stepsize = stepsize, maxiters = 25_000
+        )
+        sol = OptimizationBase.solve(prob, opt)
+        @test sol.objective < 1.0e-10
+        @test SciMLBase.successful_retcode(sol)
+    end
+
     @testset "Nelder-Mead" begin
         x0 = zeros(2)
         p = [1.0, 100.0]