CMA-ES (#373)

* initial work on CMA-ES

* new reference

* transport p vectors, better printing and documentation

* cond stopping criterion and some docs

* add test file to runtests

* even more stopping criteria

* Add TolFunCondition

* poorly conditioned test case and some docs

* renaming

* a bit of docs

* tests for stopping criteria

* Performance improvements

* bump tolerance

* Renaming, new test

* some docs

* Apply suggestions from code review

Co-authored-by: Ronny Bergmann <[email protected]>

* address some review comments

* addressing review

* coverage

* add defaults, more asserts

* Update docs

* minor docs improvements

---------

Co-authored-by: Ronny Bergmann <[email protected]>

Author: mateuszbaran

Changelog.md

@@ -5,6 +5,12 @@ All notable Changes to the Julia package `Manopt.jl` will be documented in this
 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.4.59] - unreleased
+
+### Added
+
+* A Riemannian variant of the CMA-ES (Covariance Matrix Adaptation Evolutionary Strategy) algorithm, `cma_es`.
+
 ## [0.4.58] - March 18, 2024
 
 ### Added

docs/make.jl

@@ -170,6 +170,7 @@ makedocs(;
             "Alternating Gradient Descent" => "solvers/alternating_gradient_descent.md",
             "Augmented Lagrangian Method" => "solvers/augmented_Lagrangian_method.md",
             "Chambolle-Pock" => "solvers/ChambollePock.md",
+            "CMA-ES" => "solvers/cma_es.md",
             "Conjugate gradient descent" => "solvers/conjugate_gradient_descent.md",
             "Convex bundle method" => "solvers/convex_bundle_method.md",
             "Cyclic Proximal Point" => "solvers/cyclic_proximal_point.md",

docs/src/references.bib

@@ -264,6 +264,21 @@ @article{ChambollePock:2011
     VOLUME       = {40}
 }
 
+
+@article{ColuttoFruhaufFuchsScherzer:2010,
+	title = {The {CMA}-{ES} on {Riemannian} {Manifolds} to {Reconstruct} {Shapes} in 3-{D} {Voxel} {Images}},
+	volume = {14},
+	issn = {1941-0026, 1089-778X},
+	url = {http://ieeexplore.ieee.org/document/5299260/},
+	doi = {10.1109/TEVC.2009.2029567},
+	number = {2},
+	journal = {IEEE Transactions on Evolutionary Computation},
+	author = {Colutto, S. and Fruhauf, F. and Fuchs, M. and Scherzer, O.},
+	month = apr,
+	year = {2010},
+	pages = {227--245},
+}
+
 @book{ConnGouldToint:2000,
     DOI       = {10.1137/1.9780898719857},
     YEAR      = {2000},
@@ -412,6 +427,15 @@ @article{HagerZhang:2006
     VOLUME       = {2},
     YEAR         = {2006}
 }
+
+@article{Hansen:2023,
+	AUTHOR = {Hansen, Nikolaus},
+	TITLE = {The {CMA} {Evolution} {Strategy}: {A} {Tutorial}},
+    JOURNAL = {ArXiv Preprint},
+	URL = {http://arxiv.org/abs/1604.00772},
+    NUMBER = {1604.00772},
+	YEAR = {2023},
+}
 @article{HestenesStiefel:1952,
     AUTHOR    = {M.R. Hestenes and E. Stiefel},
     DOI       = {10.6028/jres.049.044},

docs/src/solvers/cma_es.md

@@ -0,0 +1,57 @@
+# Covariance matrix adaptation evolutionary strategy
+
+```@meta
+CurrentModule = Manopt
+```
+
+The CMA-ES algorithm has been implemented based on [Hansen:2023](@cite) with basic Riemannian adaptations, related to transport of covariance matrix and its update vectors. Other attempts at adapting CMA-ES to Riemannian optimzation include [ColuttoFruhaufFuchsScherzer:2010](@cite).
+The algorithm is suitable for global optimization.
+
+Covariance matrix transport between consecutive mean points is handled by `eigenvector_transport!` function which is based on the idea of transport of matrix eigenvectors.
+
+```@docs
+cma_es
+```
+
+## State
+
+```@docs
+CMAESState
+```
+
+## Stopping Criteria
+
+```@docs
+StopWhenBestCostInGenerationConstant
+StopWhenCovarianceIllConditioned
+StopWhenEvolutionStagnates
+StopWhenPopulationCostConcentrated
+StopWhenPopulationDiverges
+StopWhenPopulationStronglyConcentrated
+```
+
+## [Technical details](@id sec-cma-es-technical-details)
+
+The [`cma_es`](@ref) solver requires the following functions of a manifold to be available
+
+* A [`retract!`](@extref ManifoldsBase :doc:`retractions`)`(M, q, p, X)`; it is recommended to set the [`default_retraction_method`](@extref `ManifoldsBase.default_retraction_method-Tuple{AbstractManifold}`) to a favourite retraction. If this default is set, a `retraction_method=` does not have to be specified.
+* A [`vector_transport_to!`](@extref ManifoldsBase :doc:`vector_transports`)`M, Y, p, X, q)`; it is recommended to set the [`default_vector_transport_method`](@extref `ManifoldsBase.default_vector_transport_method-Tuple{AbstractManifold}`) to a favourite retraction. If this default is set, a `vector_transport_method=` does not have to be specified.
+* A [`copyto!`](@extref `Base.copyto!-Tuple{AbstractManifold, Any, Any}`)`(M, q, p)` and [`copy`](@extref `Base.copy-Tuple{AbstractManifold, Any}`)`(M,p)` for points and similarly `copy(M, p, X)` for tangent vectors.
+* [`get_coordinates!`](@extref `ManifoldsBase.get_coordinates-Tuple{AbstractManifold, Any, Any, ManifoldsBase.AbstractBasis}`)`(M, Y, p, X, b)` and [`get_vector!`](@extref `ManifoldsBase.get_vector-Tuple{AbstractManifold, Any, Any, ManifoldsBase.AbstractBasis}`)`(M, X, p, c, b)` with respect to the [`AbstractBasis`](@extref `ManifoldsBase.AbstractBasis`) `b` provided, which is [`DefaultOrthonormalBasis`](@extref `ManifoldsBase.DefaultOrthonormalBasis`) by default from the `basis=` keyword.
+* An [`is_flat`](@extref `ManifoldsBase.is_flat-Tuple{AbstractManifold}`)`(M)`.
+
+## Internal helpers
+
+You may add new methods to `eigenvector_transport!` if you know a more optimized implementation
+for your manifold.
+
+```@docs
+Manopt.eigenvector_transport!
+```
+
+## Literature
+
+```@bibliography
+Pages = ["cma_es.md"]
+Canonical=false
+```

docs/src/solvers/quasi_Newton.md

@@ -125,7 +125,7 @@ The [`quasi_Newton`](@ref) solver requires the following functions of a manifold
 * A [`vector_transport_to!`](@extref ManifoldsBase :doc:`vector_transports`)`M, Y, p, X, q)`; it is recommended to set the [`default_vector_transport_method`](@extref `ManifoldsBase.default_vector_transport_method-Tuple{AbstractManifold}`) to a favourite retraction. If this default is set, a `vector_transport_method=` or `vector_transport_method_dual=` (for ``\mathcal N``) does not have to be specified.
 * By default quasi Newton uses [`ArmijoLinesearch`](@ref) which requires [`max_stepsize`](@ref)`(M)` to be set and an implementation of [`inner`](@extref `ManifoldsBase.inner-Tuple{AbstractManifold, Any, Any, Any}`)`(M, p, X)`.
 * the [`norm`](@extref `LinearAlgebra.norm-Tuple{AbstractManifold, Any, Any}`) as well, to stop when the norm of the gradient is small, but if you implemented `inner`, the norm is provided already.
-* A [`copyto!](@extref `Base.copyto!-Tuple{AbstractManifold, Any, Any}`)`(M, q, p)` and [`copy`](@extref `Base.copy-Tuple{AbstractManifold, Any}`)`(M,p)` for points and similarly `copy(M, p, X)` for tangent vectors.
+* A [`copyto!`](@extref `Base.copyto!-Tuple{AbstractManifold, Any, Any}`)`(M, q, p)` and [`copy`](@extref `Base.copy-Tuple{AbstractManifold, Any}`)`(M,p)` for points and similarly `copy(M, p, X)` for tangent vectors.
 * By default the tangent vector storing the gradient is initialized calling [`zero_vector`](@extref `ManifoldsBase.zero_vector-Tuple{AbstractManifold, Any}`)`(M,p)`.
 
 Most Hessian approximations further require [`get_coordinates`](@extref `ManifoldsBase.get_coordinates-Tuple{AbstractManifold, Any, Any, ManifoldsBase.AbstractBasis}`)`(M, p, X, b)` with respect to the [`AbstractBasis`](@extref `ManifoldsBase.AbstractBasis`) `b` provided, which is [`DefaultOrthonormalBasis`](@extref `ManifoldsBase.DefaultOrthonormalBasis`) by default from the `basis=` keyword.

src/Manopt.jl

@@ -18,7 +18,20 @@ using Colors
 using DataStructures: CircularBuffer, capacity, length, push!, size, isfull
 using Dates: Millisecond, Nanosecond, Period, canonicalize, value
 using LinearAlgebra:
-    Diagonal, I, eigen, eigvals, tril, Symmetric, dot, cholesky, eigmin, opnorm
+    cond,
+    Diagonal,
+    I,
+    Eigen,
+    eigen,
+    eigen!,
+    eigvals,
+    tril,
+    Symmetric,
+    dot,
+    cholesky,
+    eigmin,
+    opnorm,
+    mul!
 using ManifoldDiff:
     adjoint_differential_log_argument,
     adjoint_differential_log_argument!,
@@ -93,6 +106,7 @@ using ManifoldsBase:
     inner,
     inverse_retract,
     inverse_retract!,
+    is_flat,
     is_point,
     is_vector,
     log,
@@ -102,6 +116,7 @@ using ManifoldsBase:
     mid_point!,
     norm,
     number_eltype,
+    number_of_coordinates,
     power_dimensions,
     project,
     project!,
@@ -128,10 +143,10 @@ using Markdown
 using Preferences:
     @load_preference, @set_preferences!, @has_preference, @delete_preferences!
 using Printf
-using Random: shuffle!, rand, randperm
+using Random: AbstractRNG, default_rng, shuffle!, rand, randn!, randperm
 using Requires
 using SparseArrays
-using Statistics: cor, cov, mean, std
+using Statistics: cor, cov, mean, median, std
 
 include("plans/plan.jl")
 # solvers general framework
@@ -142,6 +157,7 @@ include("solvers/alternating_gradient_descent.jl")
 include("solvers/augmented_Lagrangian_method.jl")
 include("solvers/convex_bundle_method.jl")
 include("solvers/ChambollePock.jl")
+include("solvers/cma_es.jl")
 include("solvers/conjugate_gradient_descent.jl")
 include("solvers/cyclic_proximal_point.jl")
 include("solvers/difference_of_convex_algorithm.jl")
@@ -416,6 +432,8 @@ export adaptive_regularization_with_cubics,
     convex_bundle_method!,
     ChambollePock,
     ChambollePock!,
+    cma_es,
+    cma_es!,
     conjugate_gradient_descent,
     conjugate_gradient_descent!,
     cyclic_proximal_point,
@@ -477,18 +495,24 @@ export StopAfter,
     StopWhenAll,
     StopWhenAllLanczosVectorsUsed,
     StopWhenAny,
+    StopWhenBestCostInGenerationConstant,
     StopWhenChangeLess,
     StopWhenCostLess,
     StopWhenCostNaN,
+    StopWhenCovarianceIllConditioned,
     StopWhenCurvatureIsNegative,
     StopWhenEntryChangeLess,
+    StopWhenEvolutionStagnates,
     StopWhenGradientChangeLess,
     StopWhenGradientNormLess,
     StopWhenFirstOrderProgress,
     StopWhenIterateNaN,
     StopWhenLagrangeMultiplierLess,
     StopWhenModelIncreased,
+    StopWhenPopulationCostConcentrated,
     StopWhenPopulationConcentrated,
+    StopWhenPopulationDiverges,
+    StopWhenPopulationStronglyConcentrated,
     StopWhenSmallerOrEqual,
     StopWhenStepsizeLess,
     StopWhenSubgradientNormLess,

src/plans/stopping_criterion.jl

@@ -901,9 +901,22 @@ mutable struct StopWhenAny{TCriteria<:Tuple} <: StoppingCriterionSet
     StopWhenAny(c::Vector{StoppingCriterion}) = new{typeof(tuple(c...))}(tuple(c...), "")
     StopWhenAny(c::StoppingCriterion...) = new{typeof(c)}(c, "")
 end
+
+# _fast_any(f, tup::Tuple) is functionally equivalent to any(f, tup) but on Julia 1.10
+# this implementation is faster on heterogeneous tuples
+@inline _fast_any(f, tup::Tuple{}) = true
+@inline _fast_any(f, tup::Tuple{T}) where {T} = f(tup[1])
+@inline function _fast_any(f, tup::Tuple)
+    if f(tup[1])
+        return true
+    else
+        return _fast_any(f, tup[2:end])
+    end
+end
+
 function (c::StopWhenAny)(p::AbstractManoptProblem, s::AbstractManoptSolverState, i::Int)
     (i == 0) && (c.reason = "") # reset on init
-    if any(subC -> subC(p, s, i), c.criteria)
+    if _fast_any(subC -> subC(p, s, i), c.criteria)
         c.reason = string((get_reason(subC) for subC in c.criteria)...)
         return true
     end
@@ -957,6 +970,8 @@ function Base.:|(s1::StopWhenAny, s2::T) where {T<:StoppingCriterion}
     return StopWhenAny(s1.criteria..., s2)
 end
 
+is_active_stopping_criterion(c::StoppingCriterion) = !isempty(c.reason)
+
 @doc raw"""
     get_active_stopping_criteria(c)
 
@@ -974,7 +989,7 @@ end
 # for non-array containing stopping criteria, the recursion ends in either
 # returning nothing or an 1-element array containing itself
 function get_active_stopping_criteria(c::sC) where {sC<:StoppingCriterion}
-    if c.reason != ""
+    if is_active_stopping_criterion(c)
         return [c] # recursion top
     else
         return []