GitHub

Manopt.jl

Manopt.jl is a package providing solvers for optimization problems defined on Riemannian manifolds. The implementation is based on ManifoldsBase.jl interface and can hence be used for all maniolds defined in Manifolds or any other manifold implemented using the interface.

Installation: OptimizationManopt.jl

To use the Optimization.jl interface to Manopt, install the OptimizationManopt package:

import Pkg;
Pkg.add("OptimizationManopt");

Methods

The following methods are available for the OptimizationManopt package:

The common kwargs maxiters, maxtime and abstol are supported by all the optimizers. Solver specific kwargs from Manopt can be passed to the solve function or OptimizationProblem.

Passing maxiters/maxtime bounds the run while keeping the convergence part of the solver's own Manopt default stopping criterion active, so a run that converges before hitting the bound reports retcode = Success rather than MaxIters/MaxTime. For the derivative-free solvers (NelderMeadOptimizer, ParticleSwarmOptimizer, CMAESOptimizer) and ConvexBundleOptimizer, Manopt's own stopping criteria never indicate convergence, so a run bounded only by maxiters/maxtime reports MaxIters/MaxTime. A Manopt stopping_criterion passed as a solver kwarg is combined with the requested maxiters/maxtime bounds and replaces the default convergence criterion.

Note

The OptimizationProblem has to be passed the manifold as the manifold keyword argument.

Reexported Manopt.jl API

The optimizers listed above are defined by OptimizationManopt itself; they are thin wrappers whose names this package owns. From Manopt.jl only the Manopt module binding is re-exported, so that using OptimizationManopt is enough to write Manopt.ArmijoLinesearch(M) as the examples below do.

Manopt's solver options are passed to solve and keep that qualified spelling:

  • Stepsizes, e.g. Manopt.ArmijoLinesearch, Manopt.WolfePowellLinesearch
  • Stopping criteria, e.g. Manopt.StopAfterIteration, Manopt.StopWhenGradientNormLess
  • Quasi-Newton update rules, e.g. Manopt.InverseBFGS, Manopt.SR1
  • Debug and record actions, e.g. Manopt.DebugCost, Manopt.RecordIterate

Manopt's own surface is roughly 480 names and includes solve!, BFGS, NelderMead and getindex, so it is deliberately not blanket-reexported: doing so shadowed the SciML solve! and collided with Optim.jl's optimizer names.

Anything else from Manopt.jl must be reached through the Manopt module or imported from Manopt directly.

Examples

The Rosenbrock function on the Euclidean manifold can be optimized using the GradientDescentOptimizer as follows:

using OptimizationBase, OptimizationManopt, Manifolds, Manopt, LinearAlgebra, ADTypes, Zygote
rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
x0 = zeros(2)
p = [1.0, 100.0]

R2 = Euclidean(2)

stepsize = Manopt.ArmijoLinesearch(R2)
opt = OptimizationManopt.GradientDescentOptimizer()

optf = OptimizationFunction(rosenbrock, ADTypes.AutoZygote())

prob = OptimizationProblem(
    optf, x0, p; manifold = R2, stepsize = stepsize, maxiters = 25000)

sol = OptimizationBase.solve(prob, opt)
retcode: Success
u: 2-element Vector{Float64}:
 0.9999999960724049
 0.9999999921503789
Note

Plain gradient descent zig-zags slowly down Rosenbrock's narrow curved valley, so it needs tens of thousands of iterations (a few seconds 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:

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)

function closed_form_solution!(M::SymmetricPositiveDefinite, q, 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
    copyto!(M, q, p_sqrt * Q * z * Q' * p_sqrt)
    return q
end
N = m
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 = 5000)

opt = OptimizationManopt.FrankWolfeOptimizer()
sol = OptimizationBase.solve(
    prob, opt, sub_problem = (M, q, p, X) -> closed_form_solution!(M, q, L, U, p, X),
    evaluation = Manopt.InplaceEvaluation())
retcode: Success
u: 5×5 Matrix{Float64}:
 2.99987  1.99996  2.00009  2.00015  2.00023
 1.99996  3.00011  2.00025  2.0003   2.00042
 2.00009  2.00025  3.00023  2.00046  2.00043
 2.00015  2.0003   2.00046  3.00055  2.00055
 2.00023  2.00042  2.00043  2.00055  3.00068

This example is based on the example in the Manopt and Weber and Sra'22.

The following example is adapted from the Rayleigh Quotient example in ManoptExamples.jl. We solve the Rayleigh quotient problem on the Sphere manifold:

using OptimizationBase, OptimizationManopt
using Manifolds, LinearAlgebra
using Manopt

n = 1000
A = Symmetric(randn(n, n) / n)
manifold = Sphere(n - 1)

cost(x, p = nothing) = -x' * A * x
egrad(G, x, p = nothing) = (G .= -2 * A * x)

optf = OptimizationFunction(cost, grad = egrad)
x0 = rand(manifold)
prob = OptimizationProblem(optf, x0, manifold = manifold, maxiters = 5000)

sol = solve(prob, GradientDescentOptimizer())
retcode: Success
u: 1000-element Vector{Float64}:
 -0.011766617078942507
  0.0011301736460171375
  0.037212052134910495
 -0.0010854042476820477
 -0.02910637958470677
 -0.010027411279128299
  0.018528347162582327
 -0.04658913149259999
 -0.00602503999436661
  0.004507448669318531
  ⋮
 -0.008969868061594862
 -0.031182428554254973
  0.046921107433974985
  0.015196010094263558
 -0.06296720176033874
  0.051309014021110474
  0.09275834349850441
  0.014485601708695119
 -0.06077110571906837

Let's check that this indeed corresponds to the minimum eigenvalue of the matrix A.

@show eigmin(A)
@show sol.objective
-0.06346378612505356