Fix a few bugs in the existing code.

Author: kellertuer

docs/src/optimization_packages/manopt.md

@@ -1,7 +1,9 @@
 # Manopt.jl
 
-[Manopt.jl](https://github.com/JuliaManifolds/Manopt.jl) is a package with implementations of a variety of optimization solvers on manifolds supported by
-[Manifolds](https://github.com/JuliaManifolds/Manifolds.jl).
+[Manopt.jl](https://github.com/JuliaManifolds/Manopt.jl) is a package providing solvers
+for optimization problems defined on Riemannian manifolds.
+The implementation is based on [ManifoldsBase.jl](https://github.com/JuliaManifolds/ManifoldsBase.jl) interface and can hence be used for all maniolds defined in
+[Manifolds](https://github.com/JuliaManifolds/Manifolds.jl) or any other manifold implemented using the interface.
 
 ## Installation: OptimizationManopt.jl
 
@@ -29,7 +31,7 @@ The common kwargs `maxiters`, `maxtime` and `abstol` are supported by all the op
 function or `OptimizationProblem`.
 
 !!! note
-    
+
     The `OptimizationProblem` has to be passed the manifold as the `manifold` keyword argument.
 
 ## Examples

lib/OptimizationManopt/src/OptimizationManopt.jl

@@ -65,6 +65,7 @@ function call_manopt_optimizer(
         loss,
         gradF,
         x0;
+        hessF=nothing, # ignore that keyword for this solver
         kwargs...)
     opts = Manopt.gradient_descent(M,
         loss,
@@ -84,6 +85,7 @@ function call_manopt_optimizer(M::ManifoldsBase.AbstractManifold, opt::NelderMea
         loss,
         gradF,
         x0;
+        hessF=nothing, # ignore that keyword for this solver
         kwargs...)
     opts = NelderMead(M, loss; return_state = true, kwargs...)
     minimizer = Manopt.get_solver_result(opts)
@@ -98,6 +100,7 @@ function call_manopt_optimizer(M::ManifoldsBase.AbstractManifold,
         loss,
         gradF,
         x0;
+        hessF=nothing, # ignore that keyword for this solver
         kwargs...)
     opts = Manopt.conjugate_gradient_descent(M,
         loss,
@@ -106,7 +109,6 @@ function call_manopt_optimizer(M::ManifoldsBase.AbstractManifold,
         return_state = true,
         kwargs...
     )
-    # we unwrap DebugOptions here
     minimizer = Manopt.get_solver_result(opts)
     return (; minimizer = minimizer, minimum = loss(M, minimizer), options = opts)
 end
@@ -119,6 +121,7 @@ function call_manopt_optimizer(M::ManifoldsBase.AbstractManifold,
         loss,
         gradF,
         x0;
+        hessF=nothing, # ignore that keyword for this solver
         population_size::Int = 100,
         kwargs...)
     swarm = [x0, [rand(M) for _ in 1:(population_size - 1)]...]
@@ -136,10 +139,10 @@ function call_manopt_optimizer(M::Manopt.AbstractManifold,
         loss,
         gradF,
         x0;
+        hessF=nothing, # ignore that keyword for this solver
         kwargs...
 )
     opts = quasi_Newton(M, loss, gradF, x0; return_state = true, kwargs...)
-    # we unwrap DebugOptions here
     minimizer = Manopt.get_solver_result(opts)
     return (; minimizer = minimizer, minimum = loss(M, minimizer), options = opts)
 end
@@ -151,9 +154,9 @@ function call_manopt_optimizer(M::ManifoldsBase.AbstractManifold,
         loss,
         gradF,
         x0;
+        hessF=nothing, # ignore that keyword for this solver
         kwargs...)
     opt = cma_es(M, loss, x0; return_state = true, kwargs...)
-    # we unwrap DebugOptions here
     minimizer = Manopt.get_solver_result(opt)
     return (; minimizer = minimizer, minimum = loss(M, minimizer), options = opt)
 end
@@ -165,9 +168,9 @@ function call_manopt_optimizer(M::ManifoldsBase.AbstractManifold,
         loss,
         gradF,
         x0;
+        hessF=nothing, # ignore that keyword for this solver
         kwargs...)
     opt = convex_bundle_method(M, loss, gradF, x0; return_state = true, kwargs...)
-    # we unwrap DebugOptions here
     minimizer = Manopt.get_solver_result(opt)
     return (; minimizer = minimizer, minimum = loss(M, minimizer), options = opt)
 end
@@ -205,7 +208,6 @@ function call_manopt_optimizer(M::ManifoldsBase.AbstractManifold,
     else
         trust_regions(M, loss, gradF, hessF, x0; return_state = true, kwargs...)
     end
-    # we unwrap DebugOptions here
     minimizer = Manopt.get_solver_result(opt)
     return (; minimizer = minimizer, minimum = loss(M, minimizer), options = opt)
 end
@@ -217,9 +219,9 @@ function call_manopt_optimizer(M::ManifoldsBase.AbstractManifold,
         loss,
         gradF,
         x0;
+        hessF=nothing, # ignore that keyword for this solver
         kwargs...)
     opt = Frank_Wolfe_method(M, loss, gradF, x0; return_state = true, kwargs...)
-    # we unwrap DebugOptions here
     minimizer = Manopt.get_solver_result(opt)
     return (; minimizer = minimizer, minimum = loss(M, minimizer), options = opt)
 end
@@ -238,6 +240,9 @@ function SciMLBase.requireshessian(opt::Union{
 end
 
 function build_loss(f::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`,
+    # or it is not, then we need M.
     return function (::AbstractManifold, θ)
         x = f.f(θ, prob.p)
         cb(x, θ)
@@ -351,13 +356,11 @@ function SciMLBase.__solve(cache::OptimizationCache{
         stopping_criterion = Manopt.StopAfterIteration(500)
     end
 
-    # TODO: With the new keyword warnings we can not just always pass down hessF!
     opt_res = call_manopt_optimizer(manifold, cache.opt, _loss, gradF, cache.u0;
         solver_kwarg..., stopping_criterion = stopping_criterion, hessF)
 
     asc = get_stopping_criterion(opt_res.options)
-    # TODO: Switch to `has_converged` once that was released.
-    opt_ret = Manopt.indicates_convergence(asc) ? ReturnCode.Success : ReturnCode.Failure
+    opt_ret = Manopt.has_converged(asc) ? ReturnCode.Success : ReturnCode.Failure
 
     return SciMLBase.build_solution(cache,
         cache.opt,

lib/OptimizationManopt/test/runtests.jl

@@ -151,7 +151,7 @@ R2 = Euclidean(2)
 
         opt = OptimizationManopt.AdaptiveRegularizationCubicOptimizer()
 
-        #TODO: This autodiff currently provides a Hessian that seem to not procide a Hessian
+        #TODO: This autodiff currently provides a Hessian that seem to not provide a Hessian
         # ARC Fails but also AD before that warns. So it passes _some_ hessian but a wrong one, even in format
         optprob = OptimizationFunction(rosenbrock, AutoForwardDiff())
         prob = OptimizationProblem(optprob, x0, p; manifold = R2)
@@ -175,23 +175,6 @@ R2 = Euclidean(2)
         @test sol.minimum < 0.1
     end
 
-    # @testset "Circle example from Manopt" begin
-    #     Mc = Circle()
-    #     pc = 0.0
-    #     data = [-π / 4, 0.0, π / 4]
-    #     fc(y, _) = 1 / 2 * sum([distance(M, y, x)^2 for x in data])
-    #     sgrad_fc(G, y, _) = G .= -log(Mc, y, rand(data))
-
-    #     opt = OptimizationManopt.StochasticGradientDescentOptimizer()
-
-    #     optprob = OptimizationFunction(fc, grad = sgrad_fc)
-    #     prob = OptimizationProblem(optprob, pc; manifold = Mc)
-
-    #     sol = Optimization.solve(prob, opt)
-
-    #     @test all([is_point(Mc, q, true) for q in [q1, q2, q3, q4, q5]])
-    # end
-
     @testset "Custom constraints" begin
         cons(res, x, p) = (res .= [x[1]^2 + x[2]^2, x[1] * x[2]])
 
@@ -204,49 +187,4 @@ R2 = Euclidean(2)
         #TODO: What is this?
         @test_throws SciMLBase.IncompatibleOptimizerError Optimization.solve(prob_cons, opt)
     end
-
-    @testset "SPD Manifold" begin
-        M = SymmetricPositiveDefinite(5)
-        m = 100
-        σ = 0.005
-        q = Matrix{Float64}(I, 5, 5) .+ 2.0
-        data2 = [exp(M, q, σ * rand(M; vector_at = q)) for i in 1:m]
-
-        f(x, p = nothing) = sum(distance(M, x, data2[i])^2 for i in 1:m)
-
-        optf = OptimizationFunction(f, Optimization.AutoFiniteDiff())
-        prob = OptimizationProblem(optf, data2[1]; manifold = M, maxiters = 1000)
-
-        opt = OptimizationManopt.GradientDescentOptimizer()
-        @time sol = Optimization.solve(prob, opt)
-
-        @test sol.u≈q rtol=1e-2
-
-        function closed_form_solution(M::SymmetricPositiveDefinite, L, U, p, X)
-            # extract p^1/2 and p^{-1/2}
-            (p_sqrt_inv, p_sqrt) = Manifolds.spd_sqrt_and_sqrt_inv(p)
-            # Compute D & Q
-            e2 = eigen(p_sqrt_inv * X * p_sqrt_inv) # decompose Sk  = QDQ'
-            D = Diagonal(1.0 .* (e2.values .< 0))
-            Q = e2.vectors
-            Uprime = Q' * p_sqrt_inv * U * p_sqrt_inv * Q
-            Lprime = Q' * p_sqrt_inv * L * p_sqrt_inv * Q
-            P = cholesky(Hermitian(Uprime - Lprime))
-
-            z = P.U' * D * P.U + Lprime
-            return p_sqrt * Q * z * Q' * p_sqrt
-        end
-        N = m
-        U = mean(data2)
-        L = inv(sum(1 / N * inv(matrix) for matrix in data2))
-
-        opt = OptimizationManopt.FrankWolfeOptimizer()
-        optf = OptimizationFunction(f, Optimization.AutoFiniteDiff())
-        prob = OptimizationProblem(optf, data2[1]; manifold = M)
-
-        @time sol = Optimization.solve(
-            prob, opt, sub_problem = (M, p, X) -> closed_form_solution(M, p, L, U, X),
-            maxiters = 1000)
-        @test sol.u≈q rtol=1e-2
-    end
 end