diff --git a/docs/src/optimization_packages/manopt.md b/docs/src/optimization_packages/manopt.md index f80c3c0f7..1b42191a8 100644 --- a/docs/src/optimization_packages/manopt.md +++ b/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()) ``` diff --git a/lib/OptimizationManopt/src/OptimizationManopt.jl b/lib/OptimizationManopt/src/OptimizationManopt.jl index 5b7ce8f95..1d8f94848 100644 --- a/lib/OptimizationManopt/src/OptimizationManopt.jl +++ b/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`, diff --git a/lib/OptimizationManopt/test/core_tests.jl b/lib/OptimizationManopt/test/core_tests.jl index e74a0749f..2b008b4fb 100644 --- a/lib/OptimizationManopt/test/core_tests.jl +++ b/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]