From 8baa2f0a6f071389170df755dae7cf28715df746 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Fri, 27 Sep 2024 18:01:15 +0200 Subject: [PATCH 01/29] Initial sketch of (LT)MADS --- src/Manopt.jl | 1 + src/plans/mesh_adaptive_plan.jl | 37 ++++++++++++++++++ src/plans/plan.jl | 1 + src/solvers/mesh_adaptive_direct_seach.jl | 47 +++++++++++++++++++++++ 4 files changed, 86 insertions(+) create mode 100644 src/plans/mesh_adaptive_plan.jl create mode 100644 src/solvers/mesh_adaptive_direct_seach.jl diff --git a/src/Manopt.jl b/src/Manopt.jl index 4973e04b76..e2a930e7c0 100644 --- a/src/Manopt.jl +++ b/src/Manopt.jl @@ -198,6 +198,7 @@ include("solvers/FrankWolfe.jl") include("solvers/gradient_descent.jl") include("solvers/interior_point_Newton.jl") include("solvers/LevenbergMarquardt.jl") +include("solvers/mesh_adaptive_direct_seach.jl") include("solvers/particle_swarm.jl") include("solvers/primal_dual_semismooth_Newton.jl") include("solvers/proximal_bundle_method.jl") diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl new file mode 100644 index 0000000000..4f82db5222 --- /dev/null +++ b/src/plans/mesh_adaptive_plan.jl @@ -0,0 +1,37 @@ +""" + AbstractMeshPollFunctiom + +An abstract type for common “poll” strategies in the [`mesh_adaptive_direct_search`](@ref) +solver. +A subtype of this The functor has to fulfill + +* be callable as `poll!(problem, state, mesh_size)` and modify the state + +as well as + +* provide a `get_poll_success(poll!)` function that indicates whether the last poll was successful in finding a new candidate. +""" +abstract type AbstractMeshPollFunctiom end + +""" + MeshAdaptiveState <: AbstractManoptSolverState + +For a search step as a functor that is a subtype of this type, +the following is expected + +* be callable as `search!(problem, state, mesh_size)` and modify the states iterate +* return a (maybe) new `meshsize` value. +* provide a `get_search_success(search!)` that returns whether the last search was sucressul. +""" +abstract type AbstractMeshSearchFunction end + +""" + MeshAdaptiveState <: AbstractManoptSolverState +""" +mutable struct MeshAdaptiveSearchState{P,PT,ST,TStop<:StoppingCriterion} <: + AbstractManoptSolverState + p::P + stop::TStop + poll::PT + search::ST +end diff --git a/src/plans/plan.jl b/src/plans/plan.jl index 9a5465654c..438891c02b 100644 --- a/src/plans/plan.jl +++ b/src/plans/plan.jl @@ -132,6 +132,7 @@ include("exact_penalty_method_plan.jl") include("frank_wolfe_plan.jl") include("interior_point_Newton_plan.jl") include("quasi_newton_plan.jl") +include("mesh_adaptive_plan.jl") include("nonlinear_least_squares_plan.jl") include("difference_of_convex_plan.jl") include("Douglas_Rachford_plan.jl") diff --git a/src/solvers/mesh_adaptive_direct_seach.jl b/src/solvers/mesh_adaptive_direct_seach.jl new file mode 100644 index 0000000000..586e0bc130 --- /dev/null +++ b/src/solvers/mesh_adaptive_direct_seach.jl @@ -0,0 +1,47 @@ +""" + LowerTriangularPoll <: AbstractMeshPollFunctiom + +# Fields +* `last_p` – the last point in case we have to transport +* `q` – a temporary memor for a point on the manifold +* `basis` – a basis of the tangent space at `last_p` +* `directions` – a set of mesh directions. +* `L` a lower triangular matrix in coefficients of the basis +* `retraction_method` – a method to perform the retractiom +* `vector_transport_method` – a method to perform the vector transort +""" +mutable struct LowerTriangularPoll{ + P, + T, + R, + VT<:AbstractVector{T}, + B<:AbstractBasis, + A, + RM<:AbstractRetractionMethod, + VTM<:AbstractVectorTransportMethod, +} <: AbstractMeshPollFunctiom + last_p::P + q::P + basis::B + directiions::VT + L::A + mesh_size::R + retraction_method::RM + vector_transiort_method::VTM +end + +""" + DefaultSearch <: AbstractMeshPollFunctiom + +# Fields +* `q` – a temporary memor for a point on the manifold +* `X`: the search direction +* `retraction_method` – a method to perform the retractiom +* `vector_transport_method` – a method to perform the vector transort +""" +mutable struct DefaultSearch{P,T,R} <: AbstractMeshSearchFunction + q::P + X::T + redcution_factor::R + search_direction::T +end From 4c6d7407c081cab307455b2f2640bd3e5d76a30e Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Thu, 2 Jan 2025 21:43:34 +0100 Subject: [PATCH 02/29] Design concrete search and poll structs further and document them. --- docs/src/references.bib | 7 ++ src/plans/mesh_adaptive_plan.jl | 24 +++--- src/solvers/mesh_adaptive_direct_seach.jl | 98 +++++++++++++++++------ 3 files changed, 91 insertions(+), 38 deletions(-) diff --git a/docs/src/references.bib b/docs/src/references.bib index c7294640ea..2348e23e0e 100644 --- a/docs/src/references.bib +++ b/docs/src/references.bib @@ -325,6 +325,13 @@ @article{DiepeveenLellmann:2021 VOLUME = {14}, YEAR = {2021}, } +@online{Dreisigmeyer:2007, + AUTHOR = {Dreisigmayer, David W.}, + PUBLISHER = {Optimization Online}, + TITLE = {Direct Search Alogirthms over Riemannian Manifolds}, + URL = {https://optimization-online.org/?p=9134}, + YEAR = {2007} +} @article{DuranMoelleSbertCremers:2016, AUTHOR = {Duran, J. and Moeller, M. and Sbert, C. and Cremers, D.}, TITLE = {Collaborative Total Variation: A General Framework for Vectorial TV Models}, diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index 4f82db5222..03160751c8 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -1,36 +1,36 @@ """ - AbstractMeshPollFunctiom + AbstractMeshPollFunction An abstract type for common “poll” strategies in the [`mesh_adaptive_direct_search`](@ref) solver. -A subtype of this The functor has to fulfill +A subtype of this The functor has to fullfil -* be callable as `poll!(problem, state, mesh_size)` and modify the state +* be callable as `poll!(problem, mesh_size)` and modify the state as well as -* provide a `get_poll_success(poll!)` function that indicates whether the last poll was successful in finding a new candidate. +* provide a `get_poll_success(poll!)` function that indicates whether the last poll was successful in finding a new candidate, +this returns the last sucessful mesh vector used. """ abstract type AbstractMeshPollFunctiom end """ - MeshAdaptiveState <: AbstractManoptSolverState - -For a search step as a functor that is a subtype of this type, -the following is expected + AbstractMeshSearchFunction -* be callable as `search!(problem, state, mesh_size)` and modify the states iterate -* return a (maybe) new `meshsize` value. -* provide a `get_search_success(search!)` that returns whether the last search was sucressul. +Should be callable as search!(problem, mesh_size) """ abstract type AbstractMeshSearchFunction end """ MeshAdaptiveState <: AbstractManoptSolverState + + """ -mutable struct MeshAdaptiveSearchState{P,PT,ST,TStop<:StoppingCriterion} <: +mutable struct MeshAdaptiveSearchState{P,F<:Real,M,PT,ST,TStop<:StoppingCriterion} <: AbstractManoptSolverState p::P + mesh_size::F + poll_size::F stop::TStop poll::PT search::ST diff --git a/src/solvers/mesh_adaptive_direct_seach.jl b/src/solvers/mesh_adaptive_direct_seach.jl index 586e0bc130..e56a636b5c 100644 --- a/src/solvers/mesh_adaptive_direct_seach.jl +++ b/src/solvers/mesh_adaptive_direct_seach.jl @@ -1,47 +1,93 @@ """ - LowerTriangularPoll <: AbstractMeshPollFunctiom + LowerTriangularAdaptivePoll <: AbstractMesh + +Generate a mesh based on Section 6 and 7 of [Deisigmeyer:2007](@ref) # Fields -* `last_p` – the last point in case we have to transport -* `q` – a temporary memor for a point on the manifold -* `basis` – a basis of the tangent space at `last_p` -* `directions` – a set of mesh directions. -* `L` a lower triangular matrix in coefficients of the basis -* `retraction_method` – a method to perform the retractiom -* `vector_transport_method` – a method to perform the vector transort +* `p::P`: a point on the manifold, where the mesh is build in the tangent space +* `q::P`: a memory for a new point/candidate +* `already_updated::Int`: a poll counter ``l_c``, to check whether the random_vector ``b_l`` was already created +* `random_vector`: a ``d``-dimensional random vector ``b_l``` +* `random_index`: a random index ``ι`` +* `mesh` +* `basis`: a basis of the current tangent space with respect to which the mesh is stored +* `last_poll::T` the last successful poll direction stored as a tangent vector. + initiliased to the zero vector and reset to the zero vector after moving to a new tangent space. +* `vector_transport_method`: """ -mutable struct LowerTriangularPoll{ +mutable struct LowerTriangularAdaptivePoll{ P, T, - R, - VT<:AbstractVector{T}, - B<:AbstractBasis, - A, - RM<:AbstractRetractionMethod, + F<:Real, + V<:AbstractVector{F}, + M<:AbstractMatrix{F}, + I<:Int, + B, VTM<:AbstractVectorTransportMethod, -} <: AbstractMeshPollFunctiom - last_p::P + RM<:AbstractRetractionMethod, +} <: AbstractMeshPollFunction + p::P q::P + already_updated::I + random_vector::V + random_index::I + mesh::M basis::B - directiions::VT - L::A - mesh_size::R + last_poll::T + last_poll_improved::Bool retraction_method::RM - vector_transiort_method::VTM + vector_transport_method::VTM end +function LowerTriangularAdaptivePoll( + M, + p=rand(M); + basis=DefaultOrthonormalBasis(), + retraction_method=default_retraction_method(M), + vector_transport_method=default_vector_transport_method(M), +) + d = manifold_dimension(M) + b_l = zeros(d) + D_k = zeros(d, d) + last_poll = zero_vector(M, p) + return LowerTriangularAdaptiveMesh( + p, + copy(M, p), + 0, + b_l, + 0, + D_k, + basis, + last_poll, + false, + retraction_method, + vector_transport_method, + ) +end +function get_poll_success(poll!::LowerTriangularAdaptivePoll) + return nothing +end +function (poll!::LowerTriangularAdaptivePoll)(amp::AbstractManoptProblem, stepsize) + return M = get_manifold(amp) + # Implement the code from Dreisigmeyer p. 16/17 about mesh generation +end """ - DefaultSearch <: AbstractMeshPollFunctiom + DefaultSearch <: AbstractMeshSearchFunction # Fields -* `q` – a temporary memor for a point on the manifold + +* `q`: a temporary memory for a point on the manifold * `X`: the search direction +* `last_seach_improved::Bool` indicate whether the last search was succesfull, i.e. improved the cost. * `retraction_method` – a method to perform the retractiom -* `vector_transport_method` – a method to perform the vector transort """ -mutable struct DefaultSearch{P,T,R} <: AbstractMeshSearchFunction +mutable struct DefaultSearch{P,T} <: AbstractMeshSearchFunction q::P X::T - redcution_factor::R - search_direction::T + last_seach_improved::Bool +end + +function (seach!::DefaultSearch)(amp::AbstractManoptProblem, mesh_size) + return M = get_manifold(amp) + # Implement the code from Dreisigmeyer p. 17 about search generation end From fb945c2e05a9448afb77152dd5e4a295c0708e9f Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Thu, 2 Jan 2025 21:48:45 +0100 Subject: [PATCH 03/29] Add remaining todos. --- src/plans/mesh_adaptive_plan.jl | 3 +++ src/solvers/mesh_adaptive_direct_seach.jl | 15 +++++++++++---- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index 03160751c8..822f4187ab 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -35,3 +35,6 @@ mutable struct MeshAdaptiveSearchState{P,F<:Real,M,PT,ST,TStop<:StoppingCriterio poll::PT search::ST end + +# TODO: Stopping critertion based on poll_size +# TODO: Show for state \ No newline at end of file diff --git a/src/solvers/mesh_adaptive_direct_seach.jl b/src/solvers/mesh_adaptive_direct_seach.jl index e56a636b5c..3d597dfc16 100644 --- a/src/solvers/mesh_adaptive_direct_seach.jl +++ b/src/solvers/mesh_adaptive_direct_seach.jl @@ -65,7 +65,7 @@ function LowerTriangularAdaptivePoll( ) end function get_poll_success(poll!::LowerTriangularAdaptivePoll) - return nothing + return poll!.last_poll_improved end function (poll!::LowerTriangularAdaptivePoll)(amp::AbstractManoptProblem, stepsize) return M = get_manifold(amp) @@ -81,13 +81,20 @@ end * `last_seach_improved::Bool` indicate whether the last search was succesfull, i.e. improved the cost. * `retraction_method` – a method to perform the retractiom """ -mutable struct DefaultSearch{P,T} <: AbstractMeshSearchFunction +mutable struct DefaultMeshAdaptiveDirectSearch{P,T} <: AbstractMeshSearchFunction q::P X::T last_seach_improved::Bool end - -function (seach!::DefaultSearch)(amp::AbstractManoptProblem, mesh_size) +function get_search_success(search!::DefaultMeshAdaptiveDirectSearch) + return search!.last_seach_improved +end +function (search!::DefaultMeshAdaptiveDirectSearch)(amp::AbstractManoptProblem, mesh_size) return M = get_manifold(amp) # Implement the code from Dreisigmeyer p. 17 about search generation end + +# TODO: lower_triangular_mesh_adaptive_direct_search highlevel interface + +# TODO: Init solver: Do a first poll already? Probably good idea. +# TODO: step_solver to call search, poll and update both sizes. \ No newline at end of file From fa6004d046acf5c5514fa893c71e052b99118a56 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Fri, 3 Jan 2025 09:33:03 +0100 Subject: [PATCH 04/29] continue docs. --- src/plans/mesh_adaptive_plan.jl | 4 ++-- src/solvers/mesh_adaptive_direct_seach.jl | 6 ++++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index 822f4187ab..47897f681e 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -3,7 +3,7 @@ An abstract type for common “poll” strategies in the [`mesh_adaptive_direct_search`](@ref) solver. -A subtype of this The functor has to fullfil +A subtype of this The functor has to fulllfil * be callable as `poll!(problem, mesh_size)` and modify the state @@ -12,7 +12,7 @@ as well as * provide a `get_poll_success(poll!)` function that indicates whether the last poll was successful in finding a new candidate, this returns the last sucessful mesh vector used. """ -abstract type AbstractMeshPollFunctiom end +abstract type AbstractMeshPollFunction end """ AbstractMeshSearchFunction diff --git a/src/solvers/mesh_adaptive_direct_seach.jl b/src/solvers/mesh_adaptive_direct_seach.jl index 3d597dfc16..2e79edc404 100644 --- a/src/solvers/mesh_adaptive_direct_seach.jl +++ b/src/solvers/mesh_adaptive_direct_seach.jl @@ -1,7 +1,9 @@ """ - LowerTriangularAdaptivePoll <: AbstractMesh + LowerTriangularAdaptivePoll <: AbstractMeshPollFunction + +Generate a mesh (poll step) based on Section 6 and 7 of [Deisigmeyer:2007](@ref) + -Generate a mesh based on Section 6 and 7 of [Deisigmeyer:2007](@ref) # Fields * `p::P`: a point on the manifold, where the mesh is build in the tangent space From 435aea11150d592c682a18970114cde70e381964 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Fri, 3 Jan 2025 11:53:06 +0100 Subject: [PATCH 05/29] Implement most of the logic, just not yet the updates(vector transports when we move --- src/plans/mesh_adaptive_plan.jl | 25 ++- src/solvers/mesh_adaptive_direct_seach.jl | 196 +++++++++++++++++++--- 2 files changed, 190 insertions(+), 31 deletions(-) diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index 47897f681e..226fbbd430 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -5,36 +5,47 @@ An abstract type for common “poll” strategies in the [`mesh_adaptive_direct_ solver. A subtype of this The functor has to fulllfil -* be callable as `poll!(problem, mesh_size)` and modify the state +* be callable as `poll!(problem, mesh_size; kwargs...)` and modify the state as well as * provide a `get_poll_success(poll!)` function that indicates whether the last poll was successful in finding a new candidate, -this returns the last sucessful mesh vector used. +this returns the last successful mesh vector used. + +The `kwargs...` could include +* `scale_mesh=1.0`: to rescale the mesh globally +* `max_stepsize=Inf`: avoid exceeding a step size beyon this, e.g. injectivity radius. + any vector longer than this should be shortened to the provided max stepsize. """ abstract type AbstractMeshPollFunction end """ AbstractMeshSearchFunction -Should be callable as search!(problem, mesh_size) +Should be callable as search!(problem, mesh_size, X; kwargs...) + +where `X` is the last succesful poll direction, if that exists and the zero vector otherwise. """ abstract type AbstractMeshSearchFunction end """ - MeshAdaptiveState <: AbstractManoptSolverState + MeshAdaptiveDirectSearchState <: AbstractManoptSolverState +* `p`: current iterate +* `q`: temp (old) iterate """ -mutable struct MeshAdaptiveSearchState{P,F<:Real,M,PT,ST,TStop<:StoppingCriterion} <: +mutable struct MeshAdaptiveDirectSearchState{P,F<:Real,M,PT,ST,TStop<:StoppingCriterion} <: AbstractManoptSolverState p::P + q::P mesh_size::F - poll_size::F + scale_mesh::F + max_stepsize::F stop::TStop poll::PT search::ST end # TODO: Stopping critertion based on poll_size -# TODO: Show for state \ No newline at end of file +# TODO: Show for state diff --git a/src/solvers/mesh_adaptive_direct_seach.jl b/src/solvers/mesh_adaptive_direct_seach.jl index 2e79edc404..aae9ae4a73 100644 --- a/src/solvers/mesh_adaptive_direct_seach.jl +++ b/src/solvers/mesh_adaptive_direct_seach.jl @@ -1,7 +1,14 @@ """ LowerTriangularAdaptivePoll <: AbstractMeshPollFunction -Generate a mesh (poll step) based on Section 6 and 7 of [Deisigmeyer:2007](@ref) +Generate a mesh (poll step) based on Section 6 and 7 of [Deisigmeyer:2007](@ref), +with two small modifications: +* The mesh can be scaled globally so instead of ``Δ_0^m=1`` a certain different scale is used +* Any poll direction can be rescaled if it is too long. This is to not exceed the inhectivity radius for example. + +# Functor + + (p::LowerTriangularAdaptivePoll)(problem, mesh_size; scale_mesh=1.0, max_stepsize=inf) @@ -11,11 +18,12 @@ Generate a mesh (poll step) based on Section 6 and 7 of [Deisigmeyer:2007](@ref) * `already_updated::Int`: a poll counter ``l_c``, to check whether the random_vector ``b_l`` was already created * `random_vector`: a ``d``-dimensional random vector ``b_l``` * `random_index`: a random index ``ι`` -* `mesh` +* `mesh`: a vector of tangent vectors storing the mesh. * `basis`: a basis of the current tangent space with respect to which the mesh is stored -* `last_poll::T` the last successful poll direction stored as a tangent vector. - initiliased to the zero vector and reset to the zero vector after moving to a new tangent space. -* `vector_transport_method`: +* `X::T` the last successful poll direction stored as a tangent vector. + initialised to the zero vector and reset to the zero vector after moving to a new tangent space. +$(_var(:Field, :retraction_method)) +$(_var(:Field, :vector_transport_method)) """ mutable struct LowerTriangularAdaptivePoll{ P, @@ -30,12 +38,12 @@ mutable struct LowerTriangularAdaptivePoll{ } <: AbstractMeshPollFunction p::P q::P - already_updated::I + poll_counter::I random_vector::V random_index::I mesh::M basis::B - last_poll::T + X::T last_poll_improved::Bool retraction_method::RM vector_transport_method::VTM @@ -50,8 +58,8 @@ function LowerTriangularAdaptivePoll( ) d = manifold_dimension(M) b_l = zeros(d) - D_k = zeros(d, d) - last_poll = zero_vector(M, p) + D_k = zeros(d, d + 1) + X = zero_vector(M, p) return LowerTriangularAdaptiveMesh( p, copy(M, p), @@ -60,43 +68,183 @@ function LowerTriangularAdaptivePoll( 0, D_k, basis, - last_poll, + X, false, retraction_method, vector_transport_method, ) end -function get_poll_success(poll!::LowerTriangularAdaptivePoll) - return poll!.last_poll_improved +function get_poll_success(ltap::LowerTriangularAdaptivePoll) + return ltap.last_poll_improved end -function (poll!::LowerTriangularAdaptivePoll)(amp::AbstractManoptProblem, stepsize) - return M = get_manifold(amp) - # Implement the code from Dreisigmeyer p. 16/17 about mesh generation +function get_poll_direction(ltap::LowerTriangularAdaptivePoll) + return ltap.X +end +function get_poll_point(ltap::LowerTriangularAdaptivePoll) + return ltap.p +end +function (ltap::LowerTriangularAdaptivePoll)( + amp::AbstractManoptProblem, mesh_size; scale_mesh=1.0, max_stepsize=inf +) + M = get_manifold(amp) + n = manifold_dimension(M) + l = -log(4, mesh_size) + S = (-2^l + 1):(2^l - 1) + if ltap.poll_counter <= l # we did not yet generate a b_l on this scale + ltap.poll_counter += 1 + ltap.random_index = rand(1:n) + ltap.random_vector + for i in 1:n + if i == r + ltap.random_vector[i] = rand([-2^l, 2^l]) + else + ltap.random_vector[i] = rand(S) + end + end + end #otherwise we already created ltap.randomvector for this mesh size + # Generate L lower triangular, (n-1)x(n-1) in M + for i in 1:(n - 1) + for j in n - 1 + poll.mesh[i, j] = (j > i) ? 0.0 : ((i == j) ? rand([-2^l, 2^l]) : rand(S)) + end + end + # Shift to construct n × n matrix B + # (bottom left) + ltap.mesh[(ltap.random_index + 1):n, (1:n)] = poll.mesh[ + (ltap.random_index):(n - 1), (1:n) + ] + # zero row above + ltap.mesh[ltap.random_index, (1:n)] .= 0 + # left column: random vector + ltap.mesh[:, n] .= ltap.random_vector + # set last column to minus the sum. + ltap.mesh[:, n + 1] .= -1 .* sum(ltap.mesh[:, 1:n]; dims=2) + # Permute B (first n columns) + ltap.mesh[:, 1:n] .= ltap.mesh[:, randperm(n)] + # look for best + ltap.last_poll_improved = false + i = 0 + c = get_cost(amp, ltap.p) + while (i < (n + 1)) && !(ltap.last_poll_improved) + i = i + 1 # runs for the last time for i=n+1 and hence the sum. + # get vector – scale mesh + get_vector!(M, ltap.X, p, scale_mesh .* ltap.mesh[:, i], ltap.basis) + # shorten if necessary + if norm(M, ltap, p, ltap.X) > max_stepsize + ltap.X = max_stepsize & norm(M, ltap, p, ltap.X) * ltap.X + end + retract!(M, ltap.q, ltap.p, ltap.X, ltap.retraction_method) + if get_cost(amp, ltap.q) < c + copyto!(M, ltap.p, ltap, q) + ltap.last_poll_improved = true + # this also breaks while + end + end + # clear temp vector if we did not improve. + !(ltap.last_poll_improved) && (zero_vector!(M, ltap.X, p)) + return ltap end """ - DefaultSearch <: AbstractMeshSearchFunction + DefaultMeshAdaptiveDirectSearch <: AbstractMeshSearchFunction + +# Functor + + (s::DefaultMeshAdaptiveDirectSearch)(problem, mesh_size, X; scale_mesh=1.0, max_stepsize=inf) + # Fields * `q`: a temporary memory for a point on the manifold * `X`: the search direction * `last_seach_improved::Bool` indicate whether the last search was succesfull, i.e. improved the cost. -* `retraction_method` – a method to perform the retractiom +$(_var(:Field, :retraction_method)) """ mutable struct DefaultMeshAdaptiveDirectSearch{P,T} <: AbstractMeshSearchFunction - q::P + p::P X::T - last_seach_improved::Bool + last_search_improved::Bool + retraction_method::RM +end +function DefaultMeshAdaptiveDirectSearch( + M, p=rand(M); X=zero_vector(M, p), retraction_method=default_retaction_method(M) +) + return DefaultMeshAdaptiveDirectSearch(p, X, false, retracttion_method) end function get_search_success(search!::DefaultMeshAdaptiveDirectSearch) - return search!.last_seach_improved + return search!.last_search_improved end -function (search!::DefaultMeshAdaptiveDirectSearch)(amp::AbstractManoptProblem, mesh_size) - return M = get_manifold(amp) +function get_search_point(search!::DefaultMeshAdaptiveDirectSearch) + return search!.last_search_improved +end +function (dmads::DefaultMeshAdaptiveDirectSearch)( + amp::AbstractManoptProblem, mesh_size, p, X; scale_mesh=1.0, max_stepsize=inf +) + M = get_manifold(amp) + dmads.X = 4 * mesh_size * scale_mesh * X + if norm(M, p, dmads.X) > max_stepsize + dmads.X = max_stepsize / norm(M, dmads.p, dmads.X) * dmads.X + end + retract!(M, dmads.p, p, dmads.X, dmads.retraction_method) + dmads.last_search_improved = get_cost(amp, dmads.q) < get_cost(amp, p) # Implement the code from Dreisigmeyer p. 17 about search generation + return dmads end # TODO: lower_triangular_mesh_adaptive_direct_search highlevel interface -# TODO: Init solver: Do a first poll already? Probably good idea. -# TODO: step_solver to call search, poll and update both sizes. \ No newline at end of file +# Init already do a poll, since the first search requires a poll +function initialize_solver!( + amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearchState +) + # do one poll step + return madss.poll( + amp, madss.mesh_size; scale_mesh=madss.scale_mesh, max_stepsize=madss.max_stepsize + ) +end +# TODO: step_solver to call search, poll and update both sizes. +function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearchState, k) + M = get_manifolds(amp) + n = manifold_dimension(M) + copyto!(M, madss.q, madss.p) # copy p to store previous q + # search + # TODO: update search if we moved -> PT X + # update_search!(poll, madss.p) + madss.search( + amp, + madss.mesh_size, + madss.mesh_size, + get_poll_direction(madss.poll); + scale_mesh=madss.scale_mesh, + max_stepsize=madss.max_stepsize, + ) + # For sucesful search, copy over iterate - skip poll + if get_search_success(madss.seachr) + copyto!(M, madss.p, get_search_point(madss.search)) + else #search was not sucessful: poll + #TODO: update poll basis -> vector transport from poll.p to madss.p + # * at least poll.X + # * better also matrix + # * probably also basis if cached + # + # update_poll!(poll, madss.p) + # + madss.poll( + amp, + madss.mesh_size; + scale_mesh=madss.scale_mesh, + max_stepsize=madss.max_stepsize, + ) + # For succesfull poll, copy over iterate + if get_poll_success(madss.poll) + copyto!(M, madss.p, get_poll_point(madss.search)) + end + end + # If neither found a better candidate -> reduce step size, we might be close already! + if !(get_poll_success(madss.poll)) && !(get_search_success(madss.search)) + madss.mesh_size /= 4 + elseif madss.mesh_size < 0.25 # else + madss.mesh_size *= 4 # Coarsen the mesh but not beyond 1 + end + # Update poll size parameter + return madss.poll_size = n * sqrt(madss.mesh_size) +end From 88b96834034fcfa9656e47b90bcc7f94cd7e1092 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Fri, 3 Jan 2025 11:53:52 +0100 Subject: [PATCH 06/29] forgot to store poll_size. --- src/plans/mesh_adaptive_plan.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index 226fbbd430..766a13831e 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -42,6 +42,7 @@ mutable struct MeshAdaptiveDirectSearchState{P,F<:Real,M,PT,ST,TStop<:StoppingCr mesh_size::F scale_mesh::F max_stepsize::F + poll_size::F stop::TStop poll::PT search::ST From 3c2b53755b89cc9455118eb5d93b1f14730a14bf Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 5 Jan 2025 16:11:36 +0100 Subject: [PATCH 07/29] first MADS variant that includes all necessary functions. --- src/Manopt.jl | 5 + src/plans/cache.jl | 2 - src/plans/mesh_adaptive_plan.jl | 330 +++++++++++++++++++++- src/solvers/mesh_adaptive_direct_seach.jl | 266 +++++------------ 4 files changed, 402 insertions(+), 201 deletions(-) diff --git a/src/Manopt.jl b/src/Manopt.jl index 2a46588a22..92232b9b63 100644 --- a/src/Manopt.jl +++ b/src/Manopt.jl @@ -431,6 +431,8 @@ export WolfePowellLinesearch, WolfePowellBinaryLinesearch export AbstractStateAction, StoreStateAction export has_storage, get_storage, update_storage! export objective_cache_factory +export AbstractMeshPollFunction, LowerTriangularAdaptivePoll +export AbstractMeshSearchFunction, DefaultMeshAdaptiveDirectSearch # # Direction Update Rules export DirectionUpdateRule @@ -480,6 +482,8 @@ export adaptive_regularization_with_cubics, interior_point_Newton!, LevenbergMarquardt, LevenbergMarquardt!, + mesh_adaptive_direct_search, + mesh_adaptive_direct_search!, NelderMead, NelderMead!, particle_swarm, @@ -541,6 +545,7 @@ export StopAfter, StopWhenKKTResidualLess, StopWhenLagrangeMultiplierLess, StopWhenModelIncreased, + StopWhenPollSizeLess, StopWhenPopulationCostConcentrated, StopWhenPopulationConcentrated, StopWhenPopulationDiverges, diff --git a/src/plans/cache.jl b/src/plans/cache.jl index 098c0c2e8f..bde956a889 100644 --- a/src/plans/cache.jl +++ b/src/plans/cache.jl @@ -250,8 +250,6 @@ which function evaluations to cache. number of (least recently used) calls to cache * `cache_sizes=Dict{Symbol,Int}()`: a named tuple or dictionary specifying the sizes individually for each cache. - - """ struct ManifoldCachedObjective{E,P,O<:AbstractManifoldObjective{<:E},C<:NamedTuple{}} <: AbstractDecoratedManifoldObjective{E,P} diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index 766a13831e..7099b1bda2 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -22,12 +22,246 @@ abstract type AbstractMeshPollFunction end """ AbstractMeshSearchFunction -Should be callable as search!(problem, mesh_size, X; kwargs...) +Should be callable as search!(problem, mesh_size, p, X; kwargs...) -where `X` is the last succesful poll direction, if that exists and the zero vector otherwise. +where `X` is the last succesful poll direction from the tangent space at `p`` +if that exists and the zero vector otherwise. """ abstract type AbstractMeshSearchFunction end +""" + StopWhenPollSizeLess <: StoppingCriterion + +stores a threshold when to stop looking at the poll mesh size of an [`MeshAdaptiveDirectSearchState`](@ref). + +# Constructor + + StopWhenPollSizeLess(ε) + +initialize the stopping criterion to a threshold `ε`. +""" +mutable struct StopWhenPollSizeLess{F} <: StoppingCriterion + threshold::F + last_poll_size::F + at_iteration::Int + function StopWhenPollSizeLess(ε::F) where {F<:Real} + return new{F}(ε, zero(ε), -1) + end +end + +# +# Specific Polls +# +""" + LowerTriangularAdaptivePoll <: AbstractMeshPollFunction + +Generate a mesh (poll step) based on Section 6 and 7 of [Deisigmeyer:2007](@ref), +with two small modifications: +* The mesh can be scaled globally so instead of ``Δ_0^m=1`` a certain different scale is used +* Any poll direction can be rescaled if it is too long. This is to not exceed the inhectivity radius for example. + +# Functor + + (p::LowerTriangularAdaptivePoll)(problem, mesh_size; scale_mesh=1.0, max_stepsize=inf) + + + +# Fields + +* `p::P`: a point on the manifold, where the mesh is build in the tangent space +* `q::P`: a memory for a new point/candidate +* `random_vector`: a ``d``-dimensional random vector ``b_l``` +* `random_index`: a random index ``ι`` +* `mesh`: a vector of tangent vectors storing the mesh. +* `basis`: a basis of the current tangent space with respect to which the mesh is stored +* `X::T` the last successful poll direction stored as a tangent vector. + initialised to the zero vector and reset to the zero vector after moving to a new tangent space. +$(_var(:Field, :retraction_method)) +$(_var(:Field, :vector_transport_method)) +""" +mutable struct LowerTriangularAdaptivePoll{ + P, + T, + F<:Real, + V<:AbstractVector{F}, + M<:AbstractMatrix{F}, + I<:Int, + B, + VTM<:AbstractVectorTransportMethod, + RM<:AbstractRetractionMethod, +} <: AbstractMeshPollFunction + p::P + q::P + random_vector::V + random_index::I + mesh::M + basis::B + X::T + last_poll_improved::Bool + retraction_method::RM + vector_transport_method::VTM +end + +function LowerTriangularAdaptivePoll( + M, + p=rand(M); + basis=DefaultOrthonormalBasis(), + retraction_method=default_retraction_method(M), + vector_transport_method=default_vector_transport_method(M), +) + d = manifold_dimension(M) + b_l = zeros(d) + D_k = zeros(d, d + 1) + X = zero_vector(M, p) + return LowerTriangularAdaptiveMesh( + p, + copy(M, p), + 0, + b_l, + 0, + D_k, + basis, + X, + false, + retraction_method, + vector_transport_method, + ) +end +function get_poll_success(ltap::LowerTriangularAdaptivePoll) + return ltap.last_poll_improved +end +function get_poll_direction(ltap::LowerTriangularAdaptivePoll) + return ltap.X +end +function get_poll_basepoint(ltap::LowerTriangularAdaptivePoll) + return ltap.p +end +function update_poll_basepoint!(M, ltap::LowerTriangularAdaptivePoll{P}, p::P) where {P} + vector_transport_to!(M, ltap.X, ltap.p, ltap.X, p, ltap.vector_transport_method) + return copyto!(M, ltap.p, p) +end +function show(io::IO, ltap::LowerTriangularAdaptivePoll) + s = "LowerTriangularAdaptivePoll using `basis=`$(ltap.basis), `retraction_method=`$(ltap.retraction_method), and `vector_transport_method=`$(ltap.vector_transport_method)" + return print(io, s) +end +function (ltap::LowerTriangularAdaptivePoll)( + amp::AbstractManoptProblem, mesh_size; scale_mesh=1.0, max_stepsize=inf +) + M = get_manifold(amp) + n = manifold_dimension(M) + l = -log(4, mesh_size) + S = (-2^l + 1):(2^l - 1) + # A random index ι + ltap.random_index = rand(1:n) + # generate a random b_l vector + for i in 1:n + ltap.random_vector[i] = rand(i == ltap.random_index ? [-2^l, 2^l] : S) + end + # Generate L lower triangular, (n-1)x(n-1) in M + for i in 1:(n - 1) + for j in (n - 1) + ltap.mesh[i, j] = (j > i) ? 0.0 : rand((i == j) ? [-2^l, 2^l] : S) + end + end + # Shift to construct n × n matrix B + # (bottom left) + ltap.mesh[(ltap.random_index + 1):n, (1:n)] = poll.mesh[ + (ltap.random_index):(n - 1), (1:n) + ] + # zero row above + ltap.mesh[ltap.random_index, (1:n)] .= 0 + # left column: random vector + ltap.mesh[:, n] .= ltap.random_vector + # set last column to minus the sum. + ltap.mesh[:, n + 1] .= -1 .* sum(ltap.mesh[:, 1:n]; dims=2) + # Permute B (first n columns) + ltap.mesh[:, 1:n] .= ltap.mesh[:, randperm(n)] + # look for best + ltap.last_poll_improved = false + i = 0 + c = get_cost(amp, ltap.p) + while (i < (n + 1)) && !(ltap.last_poll_improved) + i = i + 1 # runs for the last time for i=n+1 and hence the sum. + # get vector – scale mesh + get_vector!(M, ltap.X, ltap.p, scale_mesh .* ltap.mesh[:, i], ltap.basis) + # shorten if necessary + if norm(M, ltap, ltap.p, ltap.X) > max_stepsize + ltap.X = max_stepsize & norm(M, ltap, ltap.p, ltap.X) * ltap.X + end + retract!(M, ltap.q, ltap.p, ltap.X, ltap.retraction_method) + if get_cost(amp, ltap.q) < c + copyto!(M, ltap.p, ltap.q) + ltap.last_poll_improved = true + # this also breaks while + end + end + # clear temp vector if we did not improve. + !(ltap.last_poll_improved) && (zero_vector!(M, ltap.X, p)) + return ltap +end + +# +# Specific Searches +# + +""" + DefaultMeshAdaptiveDirectSearch <: AbstractMeshSearchFunction + +# Functor + + (s::DefaultMeshAdaptiveDirectSearch)(problem, mesh_size, X; scale_mesh=1.0, max_stepsize=inf) + +# Fields + +* `q`: a temporary memory for a point on the manifold +* `X`: information to perform the search, e.g. the last direction found by poll. +* `last_seach_improved::Bool` indicate whether the last search was succesfull, i.e. improved the cost. +$(_var(:Field, :retraction_method)) + +# Constructor + + DefaultMeshAdaptiveDirectSearch(M::AbstractManifold, p=rand(M); kwargs...) + +## Keyword awrguments + +* `X::T=zero_vector(M,p) +$(_var(:Keyword, :retraction_method)) +""" +mutable struct DefaultMeshAdaptiveDirectSearch{P,T,RM} <: AbstractMeshSearchFunction + p::P + X::T + last_search_improved::Bool + retraction_method::RM +end +function DefaultMeshAdaptiveDirectSearch( + M, p=rand(M); X=zero_vector(M, p), retraction_method=default_retaction_method(M) +) + return DefaultMeshAdaptiveDirectSearch(p, X, false, retraction_method) +end +function get_search_success(search!::DefaultMeshAdaptiveDirectSearch) + return search!.last_search_improved +end +function get_search_point(search!::DefaultMeshAdaptiveDirectSearch) + return search!.last_search_improved +end +function show(io::IO, dmads::DefaultMeshAdaptiveDirectSearch) + s = "DefaultMeshAdaptiveDirectSearch using `retraction_method=`$(dmads.retraction_method)" + return print(io, s) +end +function (dmads::DefaultMeshAdaptiveDirectSearch)( + amp::AbstractManoptProblem, mesh_size, p, X; scale_mesh=1.0, max_stepsize=inf +) + M = get_manifold(amp) + dmads.X = 4 * mesh_size * scale_mesh * X + if norm(M, p, dmads.X) > max_stepsize + dmads.X = max_stepsize / norm(M, dmads.p, dmads.X) * dmads.X + end + retract!(M, dmads.p, p, dmads.X, dmads.retraction_method) + dmads.last_search_improved = get_cost(amp, dmads.q) < get_cost(amp, p) + # Implement the code from Dreisigmeyer p. 17 about search generation + return dmads +end + """ MeshAdaptiveDirectSearchState <: AbstractManoptSolverState @@ -47,6 +281,94 @@ mutable struct MeshAdaptiveDirectSearchState{P,F<:Real,M,PT,ST,TStop<:StoppingCr poll::PT search::ST end +function MeshAdaptiveDirectSearchState( + M::AbstractManifold, + p=rand(M); + mesh_basis::B=DefaultOrthonormalBasis(), + mesh_size=injectivity_radius(M) / 4, + scale_mesh=1.0, + max_stepsize=injectivity_radius(M), + poll_size=manifold_dimension(M) * sqrt(mesh_size), + stopping_criterion::SC=StopAfterIteration(500) | StopWhenPollSizeLess(1e-7), + retraction_method=default_retraction_method(M, typeof(p)), + vector_transport_method=default_vector_transport_method(M, typeof(p)), + poll::PT=LowerTriangularAdaptivePoll( + M, + copy(M, p); + basis=mesh_basis, + retraction_method=retraction_method, + vector_transport_method=vector_transport_method, + ), + search::ST=DefaultMeshAdaptiveDirectSearch( + M, copy(M, p); retraction_method=retraction_method + ), +) where { + PT<:AbstractMeshPollFunction, + ST<:AbstractMeshSearchFunction, + SC<:StoppingCriterion, + B<:AbstractBasis, +} + return MeshAdaptiveDirectSearchState( + p, + copy(p), + mesh_size, + scale_mesh, + max_stepsize, + poll_size, + stopping_criterion, + poll, + search, + ) +end + +function show(io::IO, mads::MeshAdaptiveDirectSearchState) + i = get_count(mads, :Iterations) + Iter = (i > 0) ? "After $i iterations\n" : "" + s = """ + # Solver state for `Manopt.jl`s mesh adaptive direct search + $Iter + + ## Parameters + * `mesh_size = ` $(mads.mesh_size) + * `scale_mesh = ` $(mads.scale_mesh) + * `max_stepsize = ` $(mads.max_stepsize) + * `poll_size = ` $(mads.poll_size) + * `poll:` $(repr(mads.poll))` + * `search:` $(repr(mads.poll))` + + ## Stopping criterion + $(status_summary(mads.stop)) + + This indicates convergence: $Conv""" + return print(io, s) +end -# TODO: Stopping critertion based on poll_size -# TODO: Show for state +get_solver_result(ips::MeshAdaptiveDirectSearchState) = ips.p + +function (c::StopWhenPollSizeLess)( + p::AbstractManoptProblem, s::MeshAdaptiveDirectSearchState, k::Int +) + if k == 0 # reset on init + c.at_iteration = -1 + end + c.last_poll_size = s.poll_size + if c.last_poll_size < c.threshold && k > 0 + c.at_iteration = k + return true + end + return false +end +function get_reason(c::StopWhenPollSizeLess) + if (c.last_poll_size < c.threshold) && (c.at_iteration >= 0) + return "The algorithm computed a poll step size ($(c.last_poll_size)) less than $(c.threshold).\n" + end + return "" +end +function status_summary(c::StopWhenPollSizeLess) + has_stopped = (c.at_iteration >= 0) + s = has_stopped ? "reached" : "not reached" + return "Poll step size s < $(c.threshold):\t$s" +end +function show(io::IO, c::StopWhenPollSizeLess) + return print(io, "StopWhenPollSizeLess($(c.threshold))\n $(status_summary(c))") +end diff --git a/src/solvers/mesh_adaptive_direct_seach.jl b/src/solvers/mesh_adaptive_direct_seach.jl index aae9ae4a73..688cb0f35d 100644 --- a/src/solvers/mesh_adaptive_direct_seach.jl +++ b/src/solvers/mesh_adaptive_direct_seach.jl @@ -1,197 +1,81 @@ -""" - LowerTriangularAdaptivePoll <: AbstractMeshPollFunction - -Generate a mesh (poll step) based on Section 6 and 7 of [Deisigmeyer:2007](@ref), -with two small modifications: -* The mesh can be scaled globally so instead of ``Δ_0^m=1`` a certain different scale is used -* Any poll direction can be rescaled if it is too long. This is to not exceed the inhectivity radius for example. - -# Functor - - (p::LowerTriangularAdaptivePoll)(problem, mesh_size; scale_mesh=1.0, max_stepsize=inf) +_doc_mads = """ + mesh_adaptive_direct_search(M, f, p=rand(M); kwargs...) + mesh_adaptive_direct_search!(M, f, p; kwargs...) -# Fields -* `p::P`: a point on the manifold, where the mesh is build in the tangent space -* `q::P`: a memory for a new point/candidate -* `already_updated::Int`: a poll counter ``l_c``, to check whether the random_vector ``b_l`` was already created -* `random_vector`: a ``d``-dimensional random vector ``b_l``` -* `random_index`: a random index ``ι`` -* `mesh`: a vector of tangent vectors storing the mesh. -* `basis`: a basis of the current tangent space with respect to which the mesh is stored -* `X::T` the last successful poll direction stored as a tangent vector. - initialised to the zero vector and reset to the zero vector after moving to a new tangent space. -$(_var(:Field, :retraction_method)) -$(_var(:Field, :vector_transport_method)) """ -mutable struct LowerTriangularAdaptivePoll{ - P, - T, - F<:Real, - V<:AbstractVector{F}, - M<:AbstractMatrix{F}, - I<:Int, - B, - VTM<:AbstractVectorTransportMethod, - RM<:AbstractRetractionMethod, -} <: AbstractMeshPollFunction - p::P - q::P - poll_counter::I - random_vector::V - random_index::I - mesh::M - basis::B - X::T - last_poll_improved::Bool - retraction_method::RM - vector_transport_method::VTM -end -function LowerTriangularAdaptivePoll( - M, - p=rand(M); - basis=DefaultOrthonormalBasis(), - retraction_method=default_retraction_method(M), - vector_transport_method=default_vector_transport_method(M), -) - d = manifold_dimension(M) - b_l = zeros(d) - D_k = zeros(d, d + 1) - X = zero_vector(M, p) - return LowerTriangularAdaptiveMesh( - p, - copy(M, p), - 0, - b_l, - 0, - D_k, - basis, - X, - false, - retraction_method, - vector_transport_method, - ) -end -function get_poll_success(ltap::LowerTriangularAdaptivePoll) - return ltap.last_poll_improved -end -function get_poll_direction(ltap::LowerTriangularAdaptivePoll) - return ltap.X +@doc "$(_doc_mads)" +mesh_adaptive_direct_search(M::AbstractManifold, args...; kwargs...) + +function mesh_adaptive_direct_search(M::AbstractManifold, f; kwargs...) + return mesh_adaptive_direct_search(M, f, rand(M); kwargs...) end -function get_poll_point(ltap::LowerTriangularAdaptivePoll) - return ltap.p +function mesh_adaptive_direct_search(M::AbstractManifold, f, p; kwargs...) + mco = ManifoldCostObjective(f) + return mesh_adaptive_direct_search(M, mco; kwargs...) end -function (ltap::LowerTriangularAdaptivePoll)( - amp::AbstractManoptProblem, mesh_size; scale_mesh=1.0, max_stepsize=inf +function mesh_adaptive_direct_search( + M::AbstractManifold, mco::AbstractManifoldCostObjective, p; kwargs... ) - M = get_manifold(amp) - n = manifold_dimension(M) - l = -log(4, mesh_size) - S = (-2^l + 1):(2^l - 1) - if ltap.poll_counter <= l # we did not yet generate a b_l on this scale - ltap.poll_counter += 1 - ltap.random_index = rand(1:n) - ltap.random_vector - for i in 1:n - if i == r - ltap.random_vector[i] = rand([-2^l, 2^l]) - else - ltap.random_vector[i] = rand(S) - end - end - end #otherwise we already created ltap.randomvector for this mesh size - # Generate L lower triangular, (n-1)x(n-1) in M - for i in 1:(n - 1) - for j in n - 1 - poll.mesh[i, j] = (j > i) ? 0.0 : ((i == j) ? rand([-2^l, 2^l]) : rand(S)) - end - end - # Shift to construct n × n matrix B - # (bottom left) - ltap.mesh[(ltap.random_index + 1):n, (1:n)] = poll.mesh[ - (ltap.random_index):(n - 1), (1:n) - ] - # zero row above - ltap.mesh[ltap.random_index, (1:n)] .= 0 - # left column: random vector - ltap.mesh[:, n] .= ltap.random_vector - # set last column to minus the sum. - ltap.mesh[:, n + 1] .= -1 .* sum(ltap.mesh[:, 1:n]; dims=2) - # Permute B (first n columns) - ltap.mesh[:, 1:n] .= ltap.mesh[:, randperm(n)] - # look for best - ltap.last_poll_improved = false - i = 0 - c = get_cost(amp, ltap.p) - while (i < (n + 1)) && !(ltap.last_poll_improved) - i = i + 1 # runs for the last time for i=n+1 and hence the sum. - # get vector – scale mesh - get_vector!(M, ltap.X, p, scale_mesh .* ltap.mesh[:, i], ltap.basis) - # shorten if necessary - if norm(M, ltap, p, ltap.X) > max_stepsize - ltap.X = max_stepsize & norm(M, ltap, p, ltap.X) * ltap.X - end - retract!(M, ltap.q, ltap.p, ltap.X, ltap.retraction_method) - if get_cost(amp, ltap.q) < c - copyto!(M, ltap.p, ltap, q) - ltap.last_poll_improved = true - # this also breaks while - end - end - # clear temp vector if we did not improve. - !(ltap.last_poll_improved) && (zero_vector!(M, ltap.X, p)) - return ltap + q = copy(M, p) + return mesh_adaptive_direct_search!(M, mco, q; kwargs...) end -""" - DefaultMeshAdaptiveDirectSearch <: AbstractMeshSearchFunction - -# Functor - - (s::DefaultMeshAdaptiveDirectSearch)(problem, mesh_size, X; scale_mesh=1.0, max_stepsize=inf) - -# Fields - -* `q`: a temporary memory for a point on the manifold -* `X`: the search direction -* `last_seach_improved::Bool` indicate whether the last search was succesfull, i.e. improved the cost. -$(_var(:Field, :retraction_method)) -""" -mutable struct DefaultMeshAdaptiveDirectSearch{P,T} <: AbstractMeshSearchFunction - p::P - X::T - last_search_improved::Bool - retraction_method::RM -end -function DefaultMeshAdaptiveDirectSearch( - M, p=rand(M); X=zero_vector(M, p), retraction_method=default_retaction_method(M) -) - return DefaultMeshAdaptiveDirectSearch(p, X, false, retracttion_method) -end -function get_search_success(search!::DefaultMeshAdaptiveDirectSearch) - return search!.last_search_improved +@doc "$(_doc_mads)" +mesh_adaptive_direct_search!(M::AbstractManifold, args...; kwargs...) +function mesh_adaptive_direct_search!(M::AbstractManifold, f, p; kwargs...) + mco = ManifoldCostObjective(f) + return mesh_adaptive_direct_search!(M, mco, p; kwargs...) end -function get_search_point(search!::DefaultMeshAdaptiveDirectSearch) - return search!.last_search_improved -end -function (dmads::DefaultMeshAdaptiveDirectSearch)( - amp::AbstractManoptProblem, mesh_size, p, X; scale_mesh=1.0, max_stepsize=inf -) - M = get_manifold(amp) - dmads.X = 4 * mesh_size * scale_mesh * X - if norm(M, p, dmads.X) > max_stepsize - dmads.X = max_stepsize / norm(M, dmads.p, dmads.X) * dmads.X - end - retract!(M, dmads.p, p, dmads.X, dmads.retraction_method) - dmads.last_search_improved = get_cost(amp, dmads.q) < get_cost(amp, p) - # Implement the code from Dreisigmeyer p. 17 about search generation - return dmads +function mesh_adaptive_direct_search!( + M::AbstractManifold, + mco::AbstractManifoldCostObjective, + p; + mesh_basis::B=DefaultOrthonormalBasis(), + mesh_size=injectivity_radius(M) / 4, + scale_mesh=1.0, + max_stepsize=injectivity_radius(M), + poll_size=manifold_dimension(M) * sqrt(mesh_size), + stopping_criterion::StoppingCriterion=StopAfterIteration(500) | + StopWhenPollSizeLess(1e-7), + retraction_method::AbstractRetractionMethod=default_retraction_method(M, eltype(p)), + vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method( + M, eltype(p) + ), + poll::PT=LowerTriangularAdaptivePoll( + M, + copy(M, p); + basis=mesh_basis, + retraction_method=retraction_method, + vector_transport_method=vector_transport_method, + ), + search::ST=DefaultMeshAdaptiveDirectSearch( + M, copy(M, p); retraction_method=retraction_method + ), + kwargs..., #collect rest +) where {B<:AbstractBasis,PT<:AbstractMeshPollFunction,ST<:AbstractMeshSearchFunction} + dmco = decorate_objective!(M, mco; kwargs...) + dmp = DefaultManoptProblem(M, dmco) + madss = MeshAdaptiveDirectSearchState( + M; + mesh_basis=mesh_basis, + mesh_size=mesh_size, + scale_mesh=scale_mesh, + max_stepsize=max_stepsize, + poll_size=poll_size, + stopping_criterion=stopping_criterion, + retraction_method=retraction_method, + vector_transport_method=vector_transport_method, + poll=poll, + search=search, + ) + dmadss = decorate_state!(madss; kwargs...) + solve!(dmp, dmadss) + return get_solver_return(get_objective(dmp), dmadss) end -# TODO: lower_triangular_mesh_adaptive_direct_search highlevel interface - # Init already do a poll, since the first search requires a poll function initialize_solver!( amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearchState @@ -201,33 +85,24 @@ function initialize_solver!( amp, madss.mesh_size; scale_mesh=madss.scale_mesh, max_stepsize=madss.max_stepsize ) end -# TODO: step_solver to call search, poll and update both sizes. function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearchState, k) M = get_manifolds(amp) n = manifold_dimension(M) copyto!(M, madss.q, madss.p) # copy p to store previous q # search - # TODO: update search if we moved -> PT X - # update_search!(poll, madss.p) madss.search( amp, madss.mesh_size, - madss.mesh_size, + get_poll_basepoint(madss.poll), get_poll_direction(madss.poll); scale_mesh=madss.scale_mesh, max_stepsize=madss.max_stepsize, ) - # For sucesful search, copy over iterate - skip poll + # For succesful search, copy over iterate - skip poll if get_search_success(madss.seachr) copyto!(M, madss.p, get_search_point(madss.search)) else #search was not sucessful: poll - #TODO: update poll basis -> vector transport from poll.p to madss.p - # * at least poll.X - # * better also matrix - # * probably also basis if cached - # - # update_poll!(poll, madss.p) - # + update_poll_basepoint!(M, madss.poll, madss.p) madss.poll( amp, madss.mesh_size; @@ -236,7 +111,7 @@ function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearc ) # For succesfull poll, copy over iterate if get_poll_success(madss.poll) - copyto!(M, madss.p, get_poll_point(madss.search)) + copyto!(M, madss.p, get_poll_basepoint(madss.search)) end end # If neither found a better candidate -> reduce step size, we might be close already! @@ -246,5 +121,6 @@ function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearc madss.mesh_size *= 4 # Coarsen the mesh but not beyond 1 end # Update poll size parameter - return madss.poll_size = n * sqrt(madss.mesh_size) + madss.poll_size = n * sqrt(madss.mesh_size) + return madss end From 354c8ffd2ff34639d0dda80800eb1bef4f9dec2a Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 5 Jan 2025 16:29:52 +0100 Subject: [PATCH 08/29] extend docs. --- docs/make.jl | 1 + docs/src/references.bib | 4 +- docs/src/solvers/index.md | 1 + .../src/solvers/mesh_adaptive_direct_seach.md | 45 +++++++++++++++++++ src/plans/mesh_adaptive_plan.jl | 2 +- 5 files changed, 50 insertions(+), 3 deletions(-) create mode 100644 docs/src/solvers/mesh_adaptive_direct_seach.md diff --git a/docs/make.jl b/docs/make.jl index fb8fb85f27..0baa04ed41 100755 --- a/docs/make.jl +++ b/docs/make.jl @@ -200,6 +200,7 @@ makedocs(; "Gradient Descent" => "solvers/gradient_descent.md", "Interior Point Newton" => "solvers/interior_point_Newton.md", "Levenberg–Marquardt" => "solvers/LevenbergMarquardt.md", + "MADS" => "solvers/mesh_adaptive_direct_seach.md", "Nelder–Mead" => "solvers/NelderMead.md", "Particle Swarm Optimization" => "solvers/particle_swarm.md", "Primal-dual Riemannian semismooth Newton" => "solvers/primal_dual_semismooth_Newton.md", diff --git a/docs/src/references.bib b/docs/src/references.bib index 2348e23e0e..d1c9745755 100644 --- a/docs/src/references.bib +++ b/docs/src/references.bib @@ -325,9 +325,9 @@ @article{DiepeveenLellmann:2021 VOLUME = {14}, YEAR = {2021}, } -@online{Dreisigmeyer:2007, +@techreport{Dreisigmeyer:2007, AUTHOR = {Dreisigmayer, David W.}, - PUBLISHER = {Optimization Online}, + INSTITUTION = {Optimization Online}, TITLE = {Direct Search Alogirthms over Riemannian Manifolds}, URL = {https://optimization-online.org/?p=9134}, YEAR = {2007} diff --git a/docs/src/solvers/index.md b/docs/src/solvers/index.md index 3227a8644f..d9ebcc3c32 100644 --- a/docs/src/solvers/index.md +++ b/docs/src/solvers/index.md @@ -20,6 +20,7 @@ In that list a 🏅 is used to indicate state-of-the-art solvers, that usually p For derivative free only function evaluations of ``f`` are used. +* [Mesh adapive direct search](mesh_adaptive_direct_seach.md) performs a mesh based exploration (poll) and search. * [Nelder-Mead](NelderMead.md) a simplex based variant, that is using ``d+1`` points, where ``d`` is the dimension of the manifold. * [Particle Swarm](particle_swarm.md) 🫏 use the evolution of a set of points, called swarm, to explore the domain of the cost and find a minimizer. * [CMA-ES](cma_es.md) uses a stochastic evolutionary strategy to perform minimization robust to local minima of the objective. diff --git a/docs/src/solvers/mesh_adaptive_direct_seach.md b/docs/src/solvers/mesh_adaptive_direct_seach.md new file mode 100644 index 0000000000..c712e066e5 --- /dev/null +++ b/docs/src/solvers/mesh_adaptive_direct_seach.md @@ -0,0 +1,45 @@ +# Mesh adaptive direct search (MADS) + + +```@meta +CurrentModule = Manopt +``` + +```@docs + mesh_adaptive_direct_search + mesh_adaptive_direct_search! +``` + +## State + +```@docs + MeshAdaptiveDirectSearchState +``` + +## Poll + +```@docs + AbstractMeshPollFunction + LowerTriangularAdaptivePoll +``` + +## Search + +```@docs + AbstractMeshSearchFunction + DefaultMeshAdaptiveDirectSearch +``` + +## Additional stopping criteria + +```@docs + StopWhenPollSizeLess +``` + +## Technical details + +The [`mesh_adaptive_direct_search`](@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. +* Within the default initialization [`rand`](@extref Base.rand-Tuple{AbstractManifold})`(M)` is used to generate the initial population +* 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. diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index 7099b1bda2..ffd51ad16a 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -55,7 +55,7 @@ end """ LowerTriangularAdaptivePoll <: AbstractMeshPollFunction -Generate a mesh (poll step) based on Section 6 and 7 of [Deisigmeyer:2007](@ref), +Generate a mesh (poll step) based on Section 6 and 7 of [Dreisigmeyer:2007](@cite), with two small modifications: * The mesh can be scaled globally so instead of ``Δ_0^m=1`` a certain different scale is used * Any poll direction can be rescaled if it is too long. This is to not exceed the inhectivity radius for example. From 7b8ad8d79cab646a46abeea7a3e9219c3fb08878 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Mon, 6 Jan 2025 09:09:50 +0100 Subject: [PATCH 09/29] Fix a few typos. --- Changelog.md | 8 +++++++- docs/make.jl | 2 +- docs/src/about.md | 1 + docs/src/solvers/index.md | 2 +- ...ive_direct_seach.md => mesh_adaptive_direct_search.md} | 0 5 files changed, 10 insertions(+), 3 deletions(-) rename docs/src/solvers/{mesh_adaptive_direct_seach.md => mesh_adaptive_direct_search.md} (100%) diff --git a/Changelog.md b/Changelog.md index 2aab0343fd..41465f2ec4 100644 --- a/Changelog.md +++ b/Changelog.md @@ -5,7 +5,13 @@ 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.5.5] Januaey 4, 2025 +## [0.5.6] unreleased + +### Added + +* A mesh adaptive direct search algorithm (MADS), for now with the LTMADS variant using a lower triangular random matrix in the poll step. + +## [0.5.5] January 4, 2025 ### Added diff --git a/docs/make.jl b/docs/make.jl index 0baa04ed41..20eb70d596 100755 --- a/docs/make.jl +++ b/docs/make.jl @@ -200,7 +200,7 @@ makedocs(; "Gradient Descent" => "solvers/gradient_descent.md", "Interior Point Newton" => "solvers/interior_point_Newton.md", "Levenberg–Marquardt" => "solvers/LevenbergMarquardt.md", - "MADS" => "solvers/mesh_adaptive_direct_seach.md", + "MADS" => "solvers/mesh_adaptive_direct_search.md", "Nelder–Mead" => "solvers/NelderMead.md", "Particle Swarm Optimization" => "solvers/particle_swarm.md", "Primal-dual Riemannian semismooth Newton" => "solvers/primal_dual_semismooth_Newton.md", diff --git a/docs/src/about.md b/docs/src/about.md index 1ea018a27a..c82bddb48a 100644 --- a/docs/src/about.md +++ b/docs/src/about.md @@ -10,6 +10,7 @@ Thanks to the following contributors to `Manopt.jl`: * [Constantin Ahlmann-Eltze](https://const-ae.name) implemented the [gradient and differential `check` functions](helpers/checks.md) * [Renée Dornig](https://github.com/r-dornig) implemented the [particle swarm](solvers/particle_swarm.md), the [Riemannian Augmented Lagrangian Method](solvers/augmented_Lagrangian_method.md), the [Exact Penalty Method](solvers/exact_penalty_method.md), as well as the [`NonmonotoneLinesearch`](@ref). These solvers are also the first one with modular/exchangable sub solvers. * [Willem Diepeveen](https://www.maths.cam.ac.uk/person/wd292) implemented the [primal-dual Riemannian semismooth Newton](solvers/primal_dual_semismooth_Newton.md) solver. +* [Sander Engen Oddsen](https://github.com/oddsen) contributed to the implementation of the [LTMADS](solvers/mesh_adaptive_direct_search.md) solver. * [Hajg Jasa](https://www.ntnu.edu/employees/hajg.jasa) implemented the [convex bundle method](solvers/convex_bundle_method.md) and the [proximal bundle method](solvers/proximal_bundle_method.md) and a default subsolver each of them. * Even Stephansen Kjemsås contributed to the implementation of the [Frank Wolfe Method](solvers/FrankWolfe.md) solver. * Mathias Ravn Munkvold contributed most of the implementation of the [Adaptive Regularization with Cubics](solvers/adaptive-regularization-with-cubics.md) solver as well as its [Lanczos](@ref arc-Lanczos) subsolver diff --git a/docs/src/solvers/index.md b/docs/src/solvers/index.md index d9ebcc3c32..bd6900e6c7 100644 --- a/docs/src/solvers/index.md +++ b/docs/src/solvers/index.md @@ -20,9 +20,9 @@ In that list a 🏅 is used to indicate state-of-the-art solvers, that usually p For derivative free only function evaluations of ``f`` are used. -* [Mesh adapive direct search](mesh_adaptive_direct_seach.md) performs a mesh based exploration (poll) and search. * [Nelder-Mead](NelderMead.md) a simplex based variant, that is using ``d+1`` points, where ``d`` is the dimension of the manifold. * [Particle Swarm](particle_swarm.md) 🫏 use the evolution of a set of points, called swarm, to explore the domain of the cost and find a minimizer. +* [Mesh adapive direct search](mesh_adaptive_direct_search.md) performs a mesh based exploration (poll) and search. * [CMA-ES](cma_es.md) uses a stochastic evolutionary strategy to perform minimization robust to local minima of the objective. ## First order diff --git a/docs/src/solvers/mesh_adaptive_direct_seach.md b/docs/src/solvers/mesh_adaptive_direct_search.md similarity index 100% rename from docs/src/solvers/mesh_adaptive_direct_seach.md rename to docs/src/solvers/mesh_adaptive_direct_search.md From 939d7b0de3abfd1a9b9636e752599d76eb31a6b4 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Tue, 7 Jan 2025 09:59:01 +0100 Subject: [PATCH 10/29] Fix two typos. --- src/plans/constrained_plan.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/plans/constrained_plan.jl b/src/plans/constrained_plan.jl index 875d12efbb..2bf002fb57 100644 --- a/src/plans/constrained_plan.jl +++ b/src/plans/constrained_plan.jl @@ -27,9 +27,9 @@ A common supertype for fucntors that model constraint functions with slack. This supertype additionally provides access for the fields * `μ::T` the dual for the inequality constraints -* `s::T` the slack parametyer, and +* `s::T` the slack parameter, and * `β::R` the the barrier parameter -which is also of typee `T`. +which is also of type `T`. """ abstract type AbstractConstrainedSlackFunctor{T,R} end From 3a4e2acd8da9714ea422c12ac04d2d143afbe7a4 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Wed, 8 Jan 2025 17:43:04 +0100 Subject: [PATCH 11/29] Fix typos add a first running, but failing test. --- src/plans/mesh_adaptive_plan.jl | 39 ++++++++++--------- src/solvers/mesh_adaptive_direct_seach.jl | 8 ++-- test/runtests.jl | 1 + .../test_mesh_adaptive_direct_search.jl | 17 ++++++++ 4 files changed, 43 insertions(+), 22 deletions(-) create mode 100644 test/solvers/test_mesh_adaptive_direct_search.jl diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index ffd51ad16a..8da894045b 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -113,10 +113,9 @@ function LowerTriangularAdaptivePoll( b_l = zeros(d) D_k = zeros(d, d + 1) X = zero_vector(M, p) - return LowerTriangularAdaptiveMesh( + return LowerTriangularAdaptivePoll( p, copy(M, p), - 0, b_l, 0, D_k, @@ -165,7 +164,7 @@ function (ltap::LowerTriangularAdaptivePoll)( end # Shift to construct n × n matrix B # (bottom left) - ltap.mesh[(ltap.random_index + 1):n, (1:n)] = poll.mesh[ + ltap.mesh[(ltap.random_index + 1):n, (1:n)] = ltap.mesh[ (ltap.random_index):(n - 1), (1:n) ] # zero row above @@ -185,8 +184,8 @@ function (ltap::LowerTriangularAdaptivePoll)( # get vector – scale mesh get_vector!(M, ltap.X, ltap.p, scale_mesh .* ltap.mesh[:, i], ltap.basis) # shorten if necessary - if norm(M, ltap, ltap.p, ltap.X) > max_stepsize - ltap.X = max_stepsize & norm(M, ltap, ltap.p, ltap.X) * ltap.X + if norm(M, ltap.p, ltap.X) > max_stepsize + ltap.X = max_stepsize / norm(M, ltap.p, ltap.X) * ltap.X end retract!(M, ltap.q, ltap.p, ltap.X, ltap.retraction_method) if get_cost(amp, ltap.q) < c @@ -196,7 +195,7 @@ function (ltap::LowerTriangularAdaptivePoll)( end end # clear temp vector if we did not improve. - !(ltap.last_poll_improved) && (zero_vector!(M, ltap.X, p)) + !(ltap.last_poll_improved) && (zero_vector!(M, ltap.X, ltap.p)) return ltap end @@ -257,8 +256,10 @@ function (dmads::DefaultMeshAdaptiveDirectSearch)( dmads.X = max_stepsize / norm(M, dmads.p, dmads.X) * dmads.X end retract!(M, dmads.p, p, dmads.X, dmads.retraction_method) - dmads.last_search_improved = get_cost(amp, dmads.q) < get_cost(amp, p) - # Implement the code from Dreisigmeyer p. 17 about search generation + dmads.last_search_improved = get_cost(amp, dmads.p) < get_cost(amp, p) + if dmads.last_search_improved + copyto!(M, dmads.p, p) + end return dmads end @@ -269,7 +270,7 @@ end * `q`: temp (old) iterate """ -mutable struct MeshAdaptiveDirectSearchState{P,F<:Real,M,PT,ST,TStop<:StoppingCriterion} <: +mutable struct MeshAdaptiveDirectSearchState{P,F<:Real,PT,ST,SC<:StoppingCriterion} <: AbstractManoptSolverState p::P q::P @@ -277,18 +278,18 @@ mutable struct MeshAdaptiveDirectSearchState{P,F<:Real,M,PT,ST,TStop<:StoppingCr scale_mesh::F max_stepsize::F poll_size::F - stop::TStop + stop::SC poll::PT search::ST end function MeshAdaptiveDirectSearchState( M::AbstractManifold, - p=rand(M); + p::P=rand(M); mesh_basis::B=DefaultOrthonormalBasis(), - mesh_size=injectivity_radius(M) / 4, - scale_mesh=1.0, - max_stepsize=injectivity_radius(M), - poll_size=manifold_dimension(M) * sqrt(mesh_size), + mesh_size::F=injectivity_radius(M) / 4, + scale_mesh::F=1.0, + max_stepsize::F=injectivity_radius(M), + poll_size::F=manifold_dimension(M) * sqrt(mesh_size), stopping_criterion::SC=StopAfterIteration(500) | StopWhenPollSizeLess(1e-7), retraction_method=default_retraction_method(M, typeof(p)), vector_transport_method=default_vector_transport_method(M, typeof(p)), @@ -303,12 +304,14 @@ function MeshAdaptiveDirectSearchState( M, copy(M, p); retraction_method=retraction_method ), ) where { + P, + F, PT<:AbstractMeshPollFunction, ST<:AbstractMeshSearchFunction, SC<:StoppingCriterion, B<:AbstractBasis, } - return MeshAdaptiveDirectSearchState( + return MeshAdaptiveDirectSearchState{P,F,PT,ST,SC}( p, copy(p), mesh_size, @@ -320,6 +323,7 @@ function MeshAdaptiveDirectSearchState( search, ) end +get_iterate(mads::MeshAdaptiveDirectSearchState) = mads.p function show(io::IO, mads::MeshAdaptiveDirectSearchState) i = get_count(mads, :Iterations) @@ -338,8 +342,7 @@ function show(io::IO, mads::MeshAdaptiveDirectSearchState) ## Stopping criterion $(status_summary(mads.stop)) - - This indicates convergence: $Conv""" + """ return print(io, s) end diff --git a/src/solvers/mesh_adaptive_direct_seach.jl b/src/solvers/mesh_adaptive_direct_seach.jl index 688cb0f35d..078053bc1c 100644 --- a/src/solvers/mesh_adaptive_direct_seach.jl +++ b/src/solvers/mesh_adaptive_direct_seach.jl @@ -39,7 +39,7 @@ function mesh_adaptive_direct_search!( max_stepsize=injectivity_radius(M), poll_size=manifold_dimension(M) * sqrt(mesh_size), stopping_criterion::StoppingCriterion=StopAfterIteration(500) | - StopWhenPollSizeLess(1e-7), + StopWhenPollSizeLess(1e-12), retraction_method::AbstractRetractionMethod=default_retraction_method(M, eltype(p)), vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method( M, eltype(p) @@ -86,7 +86,7 @@ function initialize_solver!( ) end function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearchState, k) - M = get_manifolds(amp) + M = get_manifold(amp) n = manifold_dimension(M) copyto!(M, madss.q, madss.p) # copy p to store previous q # search @@ -99,7 +99,7 @@ function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearc max_stepsize=madss.max_stepsize, ) # For succesful search, copy over iterate - skip poll - if get_search_success(madss.seachr) + if get_search_success(madss.search) copyto!(M, madss.p, get_search_point(madss.search)) else #search was not sucessful: poll update_poll_basepoint!(M, madss.poll, madss.p) @@ -111,7 +111,7 @@ function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearc ) # For succesfull poll, copy over iterate if get_poll_success(madss.poll) - copyto!(M, madss.p, get_poll_basepoint(madss.search)) + copyto!(M, madss.p, get_poll_basepoint(madss.poll)) end end # If neither found a better candidate -> reduce step size, we might be close already! diff --git a/test/runtests.jl b/test/runtests.jl index 80947d9fa1..bab88aca10 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -55,6 +55,7 @@ include("utils/example_tasks.jl") include("solvers/test_gradient_descent.jl") include("solvers/test_interior_point_Newton.jl") include("solvers/test_Levenberg_Marquardt.jl") + include("solvers/test_mesh_adaptive_direct_search.jl") include("solvers/test_Nelder_Mead.jl") include("solvers/test_proximal_bundle_method.jl") include("solvers/test_proximal_point.jl") diff --git a/test/solvers/test_mesh_adaptive_direct_search.jl b/test/solvers/test_mesh_adaptive_direct_search.jl new file mode 100644 index 0000000000..8849e45509 --- /dev/null +++ b/test/solvers/test_mesh_adaptive_direct_search.jl @@ -0,0 +1,17 @@ +using Manifolds, Manopt, Test, LinearAlgebra + +@testset "Mesh Adaptive Direct Search" begin + # A small spectral procrustes example + A = [1.0 0.0; 1.0 1.0; 0.0 1.0] + W = 1 / sqrt(2) .* [1.0 -1.0; 1.0 1.0] + B = A * W + M = Rotations(2) + p0 = [1.0 0.0; 0.0 1.0] + f(M, p) = opnorm(B - A * p) + p_s = mesh_adaptive_direct_search( + M, + f, + p0; + debug=[:Iteration, :Cost, " ", :poll_size, " ", :mesh_size, " ", :Stop, "\n"], + ) +end From 9167d93a9cbd373fd9f8bc2231f96e0e41760ddf Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Thu, 9 Jan 2025 07:45:02 +0100 Subject: [PATCH 12/29] Stabilize I still seems to reduce the size too fast --- src/plans/mesh_adaptive_plan.jl | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index 8da894045b..0b4203d29c 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -92,6 +92,7 @@ mutable struct LowerTriangularAdaptivePoll{ } <: AbstractMeshPollFunction p::P q::P + poll_counter::I random_vector::V random_index::I mesh::M @@ -116,6 +117,7 @@ function LowerTriangularAdaptivePoll( return LowerTriangularAdaptivePoll( p, copy(M, p), + 0, b_l, 0, D_k, @@ -150,25 +152,30 @@ function (ltap::LowerTriangularAdaptivePoll)( n = manifold_dimension(M) l = -log(4, mesh_size) S = (-2^l + 1):(2^l - 1) - # A random index ι - ltap.random_index = rand(1:n) - # generate a random b_l vector - for i in 1:n - ltap.random_vector[i] = rand(i == ltap.random_index ? [-2^l, 2^l] : S) - end + if ltap.poll_counter <= l # we did not yet generate a b_l on this scale + ltap.poll_counter += 1 + # A random index ι + ltap.random_index = rand(1:n) + # generate a random b_l vector + for i in 1:n + ltap.random_vector[i] = rand(i == ltap.random_index ? [-2^l, 2^l] : S) + end + end #otherwise we already created ltap.randomvector for this mesh size # Generate L lower triangular, (n-1)x(n-1) in M for i in 1:(n - 1) - for j in (n - 1) + for j in 1:(n - 1) ltap.mesh[i, j] = (j > i) ? 0.0 : rand((i == j) ? [-2^l, 2^l] : S) end end # Shift to construct n × n matrix B # (bottom left) - ltap.mesh[(ltap.random_index + 1):n, (1:n)] = ltap.mesh[ - (ltap.random_index):(n - 1), (1:n) - ] - # zero row above - ltap.mesh[ltap.random_index, (1:n)] .= 0 + if n > 1 + ltap.mesh[(ltap.random_index + 1):n, (1:n)] = ltap.mesh[ + (ltap.random_index):(n - 1), (1:n) + ] + # zero row above + ltap.mesh[ltap.random_index, (1:n)] .= 0 + end # left column: random vector ltap.mesh[:, n] .= ltap.random_vector # set last column to minus the sum. From 2b298246459e028e90aac623d375ea74898d315d Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Thu, 9 Jan 2025 17:47:09 +0100 Subject: [PATCH 13/29] Finally found the bug in scaling the mesh to be the culprit --- src/plans/mesh_adaptive_plan.jl | 88 +++++++++++-------- src/solvers/mesh_adaptive_direct_seach.jl | 39 ++++---- .../test_mesh_adaptive_direct_search.jl | 3 +- 3 files changed, 74 insertions(+), 56 deletions(-) diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index 0b4203d29c..ac159c6027 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -3,7 +3,7 @@ An abstract type for common “poll” strategies in the [`mesh_adaptive_direct_search`](@ref) solver. -A subtype of this The functor has to fulllfil +A subtype of this The functor has to fulfil * be callable as `poll!(problem, mesh_size; kwargs...)` and modify the state @@ -24,7 +24,7 @@ abstract type AbstractMeshPollFunction end Should be callable as search!(problem, mesh_size, p, X; kwargs...) -where `X` is the last succesful poll direction from the tangent space at `p`` +where `X` is the last successful poll direction from the tangent space at `p`` if that exists and the zero vector otherwise. """ abstract type AbstractMeshSearchFunction end @@ -58,7 +58,7 @@ end Generate a mesh (poll step) based on Section 6 and 7 of [Dreisigmeyer:2007](@cite), with two small modifications: * The mesh can be scaled globally so instead of ``Δ_0^m=1`` a certain different scale is used -* Any poll direction can be rescaled if it is too long. This is to not exceed the inhectivity radius for example. +* Any poll direction can be rescaled if it is too long. This is to not exceed the injectivity radius for example. # Functor @@ -68,8 +68,8 @@ with two small modifications: # Fields -* `p::P`: a point on the manifold, where the mesh is build in the tangent space -* `q::P`: a memory for a new point/candidate +* `base_point::P`: a point on the manifold, where the mesh is build in the tangent space +* `candidate::P`: a memory for a new point/candidate * `random_vector`: a ``d``-dimensional random vector ``b_l``` * `random_index`: a random index ``ι`` * `mesh`: a vector of tangent vectors storing the mesh. @@ -90,8 +90,8 @@ mutable struct LowerTriangularAdaptivePoll{ VTM<:AbstractVectorTransportMethod, RM<:AbstractRetractionMethod, } <: AbstractMeshPollFunction - p::P - q::P + base_point::P + candidate::P poll_counter::I random_vector::V random_index::I @@ -135,11 +135,19 @@ function get_poll_direction(ltap::LowerTriangularAdaptivePoll) return ltap.X end function get_poll_basepoint(ltap::LowerTriangularAdaptivePoll) - return ltap.p + return ltap.base_point +end +function get_poll_best_candidate(ltap::LowerTriangularAdaptivePoll) + return ltap.candidate end function update_poll_basepoint!(M, ltap::LowerTriangularAdaptivePoll{P}, p::P) where {P} - vector_transport_to!(M, ltap.X, ltap.p, ltap.X, p, ltap.vector_transport_method) - return copyto!(M, ltap.p, p) + vector_transport_to!( + M, ltap.X, ltap.base_point, ltap.X, p, ltap.vector_transport_method + ) + copyto!(M, ltap.base_point, p) + # reset candidate as well + copyto!(M, ltap.candidate, ltap.base_point) + return ltap end function show(io::IO, ltap::LowerTriangularAdaptivePoll) s = "LowerTriangularAdaptivePoll using `basis=`$(ltap.basis), `retraction_method=`$(ltap.retraction_method), and `vector_transport_method=`$(ltap.vector_transport_method)" @@ -176,7 +184,7 @@ function (ltap::LowerTriangularAdaptivePoll)( # zero row above ltap.mesh[ltap.random_index, (1:n)] .= 0 end - # left column: random vector + # second to last column: random vector ltap.mesh[:, n] .= ltap.random_vector # set last column to minus the sum. ltap.mesh[:, n + 1] .= -1 .* sum(ltap.mesh[:, 1:n]; dims=2) @@ -185,24 +193,34 @@ function (ltap::LowerTriangularAdaptivePoll)( # look for best ltap.last_poll_improved = false i = 0 - c = get_cost(amp, ltap.p) + c = get_cost(amp, ltap.base_point) while (i < (n + 1)) && !(ltap.last_poll_improved) i = i + 1 # runs for the last time for i=n+1 and hence the sum. # get vector – scale mesh - get_vector!(M, ltap.X, ltap.p, scale_mesh .* ltap.mesh[:, i], ltap.basis) + get_vector!( + M, + ltap.X, + ltap.base_point, + mesh_size * scale_mesh .* ltap.mesh[:, i], + ltap.basis, + ) # shorten if necessary - if norm(M, ltap.p, ltap.X) > max_stepsize - ltap.X = max_stepsize / norm(M, ltap.p, ltap.X) * ltap.X + if norm(M, ltap.base_point, ltap.X) > max_stepsize + ltap.X = max_stepsize / norm(M, ltap.base_point, ltap.X) * ltap.X end - retract!(M, ltap.q, ltap.p, ltap.X, ltap.retraction_method) - if get_cost(amp, ltap.q) < c - copyto!(M, ltap.p, ltap.q) + retract!(M, ltap.candidate, ltap.base_point, ltap.X, ltap.retraction_method) + if get_cost(amp, ltap.candidate) < c + # Keep old point and search direction, since the update will come later only + # copyto!(M, ltap.base_point, ltap.candidate) ltap.last_poll_improved = true # this also breaks while end end - # clear temp vector if we did not improve. - !(ltap.last_poll_improved) && (zero_vector!(M, ltap.X, ltap.p)) + # clear temp vector if we did not improve – set to zero vector and “clear” candidate. + if !(ltap.last_poll_improved) + zero_vector!(M, ltap.X, ltap.base_point) + copyto!(M, ltap.candidate, ltap.base_point) + end return ltap end @@ -228,13 +246,14 @@ $(_var(:Field, :retraction_method)) DefaultMeshAdaptiveDirectSearch(M::AbstractManifold, p=rand(M); kwargs...) -## Keyword awrguments +## Keyword arguments * `X::T=zero_vector(M,p) $(_var(:Keyword, :retraction_method)) """ mutable struct DefaultMeshAdaptiveDirectSearch{P,T,RM} <: AbstractMeshSearchFunction p::P + q::P X::T last_search_improved::Bool retraction_method::RM @@ -242,13 +261,13 @@ end function DefaultMeshAdaptiveDirectSearch( M, p=rand(M); X=zero_vector(M, p), retraction_method=default_retaction_method(M) ) - return DefaultMeshAdaptiveDirectSearch(p, X, false, retraction_method) + return DefaultMeshAdaptiveDirectSearch(p, copy(M, p), X, false, retraction_method) end -function get_search_success(search!::DefaultMeshAdaptiveDirectSearch) - return search!.last_search_improved +function get_search_success(dmads::DefaultMeshAdaptiveDirectSearch) + return dmads.last_search_improved end -function get_search_point(search!::DefaultMeshAdaptiveDirectSearch) - return search!.last_search_improved +function get_search_point(dmads::DefaultMeshAdaptiveDirectSearch) + return dmads.p end function show(io::IO, dmads::DefaultMeshAdaptiveDirectSearch) s = "DefaultMeshAdaptiveDirectSearch using `retraction_method=`$(dmads.retraction_method)" @@ -262,10 +281,10 @@ function (dmads::DefaultMeshAdaptiveDirectSearch)( if norm(M, p, dmads.X) > max_stepsize dmads.X = max_stepsize / norm(M, dmads.p, dmads.X) * dmads.X end - retract!(M, dmads.p, p, dmads.X, dmads.retraction_method) - dmads.last_search_improved = get_cost(amp, dmads.p) < get_cost(amp, p) + retract!(M, dmads.q, p, dmads.X, dmads.retraction_method) + dmads.last_search_improved = get_cost(amp, dmads.q) < get_cost(amp, p) if dmads.last_search_improved - copyto!(M, dmads.p, p) + copyto!(M, dmads.p, dmads.q) end return dmads end @@ -280,7 +299,6 @@ end mutable struct MeshAdaptiveDirectSearchState{P,F<:Real,PT,ST,SC<:StoppingCriterion} <: AbstractManoptSolverState p::P - q::P mesh_size::F scale_mesh::F max_stepsize::F @@ -319,15 +337,7 @@ function MeshAdaptiveDirectSearchState( B<:AbstractBasis, } return MeshAdaptiveDirectSearchState{P,F,PT,ST,SC}( - p, - copy(p), - mesh_size, - scale_mesh, - max_stepsize, - poll_size, - stopping_criterion, - poll, - search, + p, mesh_size, scale_mesh, max_stepsize, poll_size, stopping_criterion, poll, search ) end get_iterate(mads::MeshAdaptiveDirectSearchState) = mads.p diff --git a/src/solvers/mesh_adaptive_direct_seach.jl b/src/solvers/mesh_adaptive_direct_seach.jl index 078053bc1c..f7f1b7e688 100644 --- a/src/solvers/mesh_adaptive_direct_seach.jl +++ b/src/solvers/mesh_adaptive_direct_seach.jl @@ -34,12 +34,12 @@ function mesh_adaptive_direct_search!( mco::AbstractManifoldCostObjective, p; mesh_basis::B=DefaultOrthonormalBasis(), - mesh_size=injectivity_radius(M) / 4, + mesh_size=1.0, scale_mesh=1.0, max_stepsize=injectivity_radius(M), - poll_size=manifold_dimension(M) * sqrt(mesh_size), + poll_size=manifold_dimension(M) * mesh_size, stopping_criterion::StoppingCriterion=StopAfterIteration(500) | - StopWhenPollSizeLess(1e-12), + StopWhenPollSizeLess(1e-10), retraction_method::AbstractRetractionMethod=default_retraction_method(M, eltype(p)), vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method( M, eltype(p) @@ -80,27 +80,34 @@ end function initialize_solver!( amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearchState ) + M = get_manifold(amp) # do one poll step - return madss.poll( + madss.poll( amp, madss.mesh_size; scale_mesh=madss.scale_mesh, max_stepsize=madss.max_stepsize ) + if get_poll_success(madss.poll) + copyto!(M, madss.p, get_poll_best_candidate(madss.poll)) + end + return madss end function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearchState, k) M = get_manifold(amp) n = manifold_dimension(M) - copyto!(M, madss.q, madss.p) # copy p to store previous q - # search - madss.search( - amp, - madss.mesh_size, - get_poll_basepoint(madss.poll), - get_poll_direction(madss.poll); - scale_mesh=madss.scale_mesh, - max_stepsize=madss.max_stepsize, - ) - # For succesful search, copy over iterate - skip poll + # search if the last poll or last search was sucessful + if get_search_success(madss.search) || get_poll_success(madss.poll) + madss.search( + amp, + madss.mesh_size, + get_poll_best_candidate(madss.poll), + get_poll_direction(madss.poll); + scale_mesh=madss.scale_mesh, + max_stepsize=madss.max_stepsize, + ) + end + # For succesful search, copy over iterate - skip poll, but update base if get_search_success(madss.search) copyto!(M, madss.p, get_search_point(madss.search)) + update_poll_basepoint!(M, madss.poll, madss.p) else #search was not sucessful: poll update_poll_basepoint!(M, madss.poll, madss.p) madss.poll( @@ -111,7 +118,7 @@ function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearc ) # For succesfull poll, copy over iterate if get_poll_success(madss.poll) - copyto!(M, madss.p, get_poll_basepoint(madss.poll)) + copyto!(M, madss.p, get_poll_best_candidate(madss.poll)) end end # If neither found a better candidate -> reduce step size, we might be close already! diff --git a/test/solvers/test_mesh_adaptive_direct_search.jl b/test/solvers/test_mesh_adaptive_direct_search.jl index 8849e45509..83ff66b30c 100644 --- a/test/solvers/test_mesh_adaptive_direct_search.jl +++ b/test/solvers/test_mesh_adaptive_direct_search.jl @@ -12,6 +12,7 @@ using Manifolds, Manopt, Test, LinearAlgebra M, f, p0; - debug=[:Iteration, :Cost, " ", :poll_size, " ", :mesh_size, " ", :Stop, "\n"], + # debug=[:Iteration, :Cost, " ", :poll_size, " ", :mesh_size, " ", :Stop, "\n"], ) + @test distance(M, p_s, W) < 1e-9 end From 57ee145088a24b52e0ff2569218048c580e60e98 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Thu, 9 Jan 2025 17:53:32 +0100 Subject: [PATCH 14/29] Fix state print a bit. --- src/plans/mesh_adaptive_plan.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index ac159c6027..fed7765d72 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -350,12 +350,12 @@ function show(io::IO, mads::MeshAdaptiveDirectSearchState) $Iter ## Parameters - * `mesh_size = ` $(mads.mesh_size) - * `scale_mesh = ` $(mads.scale_mesh) - * `max_stepsize = ` $(mads.max_stepsize) - * `poll_size = ` $(mads.poll_size) - * `poll:` $(repr(mads.poll))` - * `search:` $(repr(mads.poll))` + * mesh_size: $(mads.mesh_size) + * scale_mesh: $(mads.scale_mesh) + * max_stepsize: $(mads.max_stepsize) + * poll_size: $(mads.poll_size) + * poll: $(repr(mads.poll)) + * search: $(repr(mads.search)) ## Stopping criterion $(status_summary(mads.stop)) From a7e9f8c81139038945ccdd34fc5b77de4f9f21ac Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Thu, 9 Jan 2025 20:17:53 +0100 Subject: [PATCH 15/29] change poll and mesh size to be internal parameters. --- src/Manopt.jl | 2 +- src/plans/mesh_adaptive_plan.jl | 22 ++++++++++++------- ...each.jl => mesh_adaptive_direct_search.jl} | 13 ++++++----- 3 files changed, 23 insertions(+), 14 deletions(-) rename src/solvers/{mesh_adaptive_direct_seach.jl => mesh_adaptive_direct_search.jl} (94%) diff --git a/src/Manopt.jl b/src/Manopt.jl index 92232b9b63..c2174340aa 100644 --- a/src/Manopt.jl +++ b/src/Manopt.jl @@ -198,7 +198,7 @@ include("solvers/FrankWolfe.jl") include("solvers/gradient_descent.jl") include("solvers/interior_point_Newton.jl") include("solvers/LevenbergMarquardt.jl") -include("solvers/mesh_adaptive_direct_seach.jl") +include("solvers/mesh_adaptive_direct_search.jl") include("solvers/particle_swarm.jl") include("solvers/primal_dual_semismooth_Newton.jl") include("solvers/proximal_bundle_method.jl") diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index fed7765d72..f6aef13816 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -64,8 +64,6 @@ with two small modifications: (p::LowerTriangularAdaptivePoll)(problem, mesh_size; scale_mesh=1.0, max_stepsize=inf) - - # Fields * `base_point::P`: a point on the manifold, where the mesh is build in the tangent space @@ -78,6 +76,15 @@ with two small modifications: initialised to the zero vector and reset to the zero vector after moving to a new tangent space. $(_var(:Field, :retraction_method)) $(_var(:Field, :vector_transport_method)) + +# Constructor + + LowerTriangularAdaptivePoll(M, p=rand(M); kwargs...) + +* `basis=`[`DefaultOrthonormalBasis`](@extref `ManifoldsBase.DefaultOrthonormalBasis`) +$(_var(:Keyword, :retraction_method)) +$(_var(:Keyword, :vector_transport_method)) +$(_var(:Keyword, :X)) """ mutable struct LowerTriangularAdaptivePoll{ P, @@ -109,11 +116,11 @@ function LowerTriangularAdaptivePoll( basis=DefaultOrthonormalBasis(), retraction_method=default_retraction_method(M), vector_transport_method=default_vector_transport_method(M), + X=zero_vector(M, p), ) d = manifold_dimension(M) b_l = zeros(d) D_k = zeros(d, d + 1) - X = zero_vector(M, p) return LowerTriangularAdaptivePoll( p, copy(M, p), @@ -311,10 +318,8 @@ function MeshAdaptiveDirectSearchState( M::AbstractManifold, p::P=rand(M); mesh_basis::B=DefaultOrthonormalBasis(), - mesh_size::F=injectivity_radius(M) / 4, - scale_mesh::F=1.0, - max_stepsize::F=injectivity_radius(M), - poll_size::F=manifold_dimension(M) * sqrt(mesh_size), + scale_mesh::F=injectivity_radius(M) / 2, + max_stepsize=injectivity_radius(M), stopping_criterion::SC=StopAfterIteration(500) | StopWhenPollSizeLess(1e-7), retraction_method=default_retraction_method(M, typeof(p)), vector_transport_method=default_vector_transport_method(M, typeof(p)), @@ -336,8 +341,9 @@ function MeshAdaptiveDirectSearchState( SC<:StoppingCriterion, B<:AbstractBasis, } + poll_s = manifold_dimension(M) * 1.0 return MeshAdaptiveDirectSearchState{P,F,PT,ST,SC}( - p, mesh_size, scale_mesh, max_stepsize, poll_size, stopping_criterion, poll, search + p, 1.0, scale_mesh, max_stepsize, poll_s, stopping_criterion, poll, search ) end get_iterate(mads::MeshAdaptiveDirectSearchState) = mads.p diff --git a/src/solvers/mesh_adaptive_direct_seach.jl b/src/solvers/mesh_adaptive_direct_search.jl similarity index 94% rename from src/solvers/mesh_adaptive_direct_seach.jl rename to src/solvers/mesh_adaptive_direct_search.jl index f7f1b7e688..2a995d12f4 100644 --- a/src/solvers/mesh_adaptive_direct_seach.jl +++ b/src/solvers/mesh_adaptive_direct_search.jl @@ -4,6 +4,13 @@ _doc_mads = """ mesh_adaptive_direct_search(M, f, p=rand(M); kwargs...) mesh_adaptive_direct_search!(M, f, p; kwargs...) + +# Keyword arguments + +* `basis=`[`DefaultOrthonormalBasis`](@extref `ManifoldsBase.DefaultOrthonormalBasis`) +$(_var(:Keyword, :retraction_method)) +$(_var(:Keyword, :vector_transport_method)) +$(_var(:Keyword, :X)) """ @doc "$(_doc_mads)" @@ -34,10 +41,8 @@ function mesh_adaptive_direct_search!( mco::AbstractManifoldCostObjective, p; mesh_basis::B=DefaultOrthonormalBasis(), - mesh_size=1.0, - scale_mesh=1.0, + scale_mesh=injectivity_radius(M) / 2, max_stepsize=injectivity_radius(M), - poll_size=manifold_dimension(M) * mesh_size, stopping_criterion::StoppingCriterion=StopAfterIteration(500) | StopWhenPollSizeLess(1e-10), retraction_method::AbstractRetractionMethod=default_retraction_method(M, eltype(p)), @@ -61,10 +66,8 @@ function mesh_adaptive_direct_search!( madss = MeshAdaptiveDirectSearchState( M; mesh_basis=mesh_basis, - mesh_size=mesh_size, scale_mesh=scale_mesh, max_stepsize=max_stepsize, - poll_size=poll_size, stopping_criterion=stopping_criterion, retraction_method=retraction_method, vector_transport_method=vector_transport_method, From e430d739d2372b394727b913c5cd8871416b6321 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Thu, 9 Jan 2025 20:41:23 +0100 Subject: [PATCH 16/29] unify naming and add docstrings to all new (small) functions --- src/plans/mesh_adaptive_plan.jl | 68 ++++++++++++++++++---- src/solvers/mesh_adaptive_direct_search.jl | 24 ++++---- 2 files changed, 69 insertions(+), 23 deletions(-) diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index f6aef13816..c59ee7a8dd 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -9,7 +9,7 @@ A subtype of this The functor has to fulfil as well as -* provide a `get_poll_success(poll!)` function that indicates whether the last poll was successful in finding a new candidate, +* provide a `is_successful(poll!)` function that indicates whether the last poll was successful in finding a new candidate, this returns the last successful mesh vector used. The `kwargs...` could include @@ -135,19 +135,47 @@ function LowerTriangularAdaptivePoll( vector_transport_method, ) end -function get_poll_success(ltap::LowerTriangularAdaptivePoll) +""" + is_successful(ltap::LowerTriangularAdaptivePoll) + +Return whether the last [`LowerTriangularAdaptivePoll`](@ref) step was successful +""" +function is_successful(ltap::LowerTriangularAdaptivePoll) return ltap.last_poll_improved end -function get_poll_direction(ltap::LowerTriangularAdaptivePoll) +""" + get_descent_direction(ltap::LowerTriangularAdaptivePoll) + +Return the direction of the last [`LowerTriangularAdaptivePoll`](@ref) that yields a descent of the cost. +If the poll was not successful, the zero vector is returned +""" +function get_descent_direction(ltap::LowerTriangularAdaptivePoll) return ltap.X end -function get_poll_basepoint(ltap::LowerTriangularAdaptivePoll) +""" + get_basepoint(ltap::LowerTriangularAdaptivePoll) + +Return the base point of the tangent space, where the mash for the [`LowerTriangularAdaptivePoll`](@ref) is build in. +""" +function get_basepoint(ltap::LowerTriangularAdaptivePoll) return ltap.base_point end -function get_poll_best_candidate(ltap::LowerTriangularAdaptivePoll) +""" + get_candidate(ltap::LowerTriangularAdaptivePoll) + +Return the candidate of the last successful [`LowerTriangularAdaptivePoll`](@ref). +If the poll was unsuccessful, the base point is returned. +""" +function get_candidate(ltap::LowerTriangularAdaptivePoll) return ltap.candidate end -function update_poll_basepoint!(M, ltap::LowerTriangularAdaptivePoll{P}, p::P) where {P} +""" + update_basepoint!(M, ltap::LowerTriangularAdaptivePoll, p) + +Update the base point of the [`LowerTriangularAdaptivePoll`](@ref). +This especially also updates the basis, that is used to build a (new) mesh. +""" +function update_basepoint!(M, ltap::LowerTriangularAdaptivePoll{P}, p::P) where {P} vector_transport_to!( M, ltap.X, ltap.base_point, ltap.X, p, ltap.vector_transport_method ) @@ -157,7 +185,7 @@ function update_poll_basepoint!(M, ltap::LowerTriangularAdaptivePoll{P}, p::P) w return ltap end function show(io::IO, ltap::LowerTriangularAdaptivePoll) - s = "LowerTriangularAdaptivePoll using `basis=`$(ltap.basis), `retraction_method=`$(ltap.retraction_method), and `vector_transport_method=`$(ltap.vector_transport_method)" + s = "LowerTriangularAdaptivePoll on a basis $(ltap.basis), the retraction_method $(ltap.retraction_method), and the vector_transport_method $(ltap.vector_transport_method)" return print(io, s) end function (ltap::LowerTriangularAdaptivePoll)( @@ -270,10 +298,20 @@ function DefaultMeshAdaptiveDirectSearch( ) return DefaultMeshAdaptiveDirectSearch(p, copy(M, p), X, false, retraction_method) end -function get_search_success(dmads::DefaultMeshAdaptiveDirectSearch) +""" + is_successful(dmads::DefaultMeshAdaptiveDirectSearch) + +Return whether the last [`DefaultMeshAdaptiveDirectSearch`](@ref) was succesful. +""" +function is_successful(dmads::DefaultMeshAdaptiveDirectSearch) return dmads.last_search_improved end -function get_search_point(dmads::DefaultMeshAdaptiveDirectSearch) +""" + get_candidate(dmads::DefaultMeshAdaptiveDirectSearch) + +Return the last candidate a [`DefaultMeshAdaptiveDirectSearch`](@ref) found +""" +function get_candidate(dmads::DefaultMeshAdaptiveDirectSearch) return dmads.p end function show(io::IO, dmads::DefaultMeshAdaptiveDirectSearch) @@ -299,8 +337,16 @@ end """ MeshAdaptiveDirectSearchState <: AbstractManoptSolverState -* `p`: current iterate -* `q`: temp (old) iterate +# Fields + +$(_var(:Field, :p; add=[:as_Iterate])) +* `mesh_size`: the current (internal) mesh size +* `scale_mesh`: the current scaling of the internal mesh size, yields the actual mesh size used +* `max_stepsize`: an upper bound for the longest step taken in looking for a candidate in either poll or search +* `poll_size` +$(_var(:Field, :stopping_criterion, "stop")) +* `poll::`[`AbstractMeshPollFunction`]: a poll step (functor) to perform +* `search::`[`AbstractMeshSearchFunction`}(@ref) a search step (functor) to perform """ mutable struct MeshAdaptiveDirectSearchState{P,F<:Real,PT,ST,SC<:StoppingCriterion} <: diff --git a/src/solvers/mesh_adaptive_direct_search.jl b/src/solvers/mesh_adaptive_direct_search.jl index 2a995d12f4..70377088b1 100644 --- a/src/solvers/mesh_adaptive_direct_search.jl +++ b/src/solvers/mesh_adaptive_direct_search.jl @@ -88,8 +88,8 @@ function initialize_solver!( madss.poll( amp, madss.mesh_size; scale_mesh=madss.scale_mesh, max_stepsize=madss.max_stepsize ) - if get_poll_success(madss.poll) - copyto!(M, madss.p, get_poll_best_candidate(madss.poll)) + if is_successful(madss.poll) + copyto!(M, madss.p, get_candidate(madss.poll)) end return madss end @@ -97,22 +97,22 @@ function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearc M = get_manifold(amp) n = manifold_dimension(M) # search if the last poll or last search was sucessful - if get_search_success(madss.search) || get_poll_success(madss.poll) + if is_successful(madss.search) || is_successful(madss.poll) madss.search( amp, madss.mesh_size, - get_poll_best_candidate(madss.poll), - get_poll_direction(madss.poll); + get_candidate(madss.poll), + get_descent_direction(madss.poll); scale_mesh=madss.scale_mesh, max_stepsize=madss.max_stepsize, ) end # For succesful search, copy over iterate - skip poll, but update base - if get_search_success(madss.search) - copyto!(M, madss.p, get_search_point(madss.search)) - update_poll_basepoint!(M, madss.poll, madss.p) + if is_successful(madss.search) + copyto!(M, madss.p, get_candidate(madss.search)) + update_basepoint!(M, madss.poll, madss.p) else #search was not sucessful: poll - update_poll_basepoint!(M, madss.poll, madss.p) + update_basepoint!(M, madss.poll, madss.p) madss.poll( amp, madss.mesh_size; @@ -120,12 +120,12 @@ function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearc max_stepsize=madss.max_stepsize, ) # For succesfull poll, copy over iterate - if get_poll_success(madss.poll) - copyto!(M, madss.p, get_poll_best_candidate(madss.poll)) + if is_successful(madss.poll) + copyto!(M, madss.p, get_candidate(madss.poll)) end end # If neither found a better candidate -> reduce step size, we might be close already! - if !(get_poll_success(madss.poll)) && !(get_search_success(madss.search)) + if !(is_successful(madss.poll)) && !(is_successful(madss.search)) madss.mesh_size /= 4 elseif madss.mesh_size < 0.25 # else madss.mesh_size *= 4 # Coarsen the mesh but not beyond 1 From 5a59142959ce01261daf206ffe0b286d1a412a92 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Thu, 9 Jan 2025 21:38:31 +0100 Subject: [PATCH 17/29] Fix docs. --- docs/src/solvers/mesh_adaptive_direct_search.md | 17 +++++++++++++++++ src/plans/mesh_adaptive_plan.jl | 15 ++++++++++++--- 2 files changed, 29 insertions(+), 3 deletions(-) diff --git a/docs/src/solvers/mesh_adaptive_direct_search.md b/docs/src/solvers/mesh_adaptive_direct_search.md index c712e066e5..deb7da8650 100644 --- a/docs/src/solvers/mesh_adaptive_direct_search.md +++ b/docs/src/solvers/mesh_adaptive_direct_search.md @@ -23,6 +23,16 @@ CurrentModule = Manopt LowerTriangularAdaptivePoll ``` +as well as the internal functions + +```@docs +Manopt.get_descent_direction(::LowerTriangularAdaptivePoll) +Manopt.is_successful(::LowerTriangularAdaptivePoll) +Manopt.get_candidate(::LowerTriangularAdaptivePoll) +Manopt.get_basepoint(::LowerTriangularAdaptivePoll) +Manopt.update_basepoint!(M, ltap::LowerTriangularAdaptivePoll{P}, p::P) where {P} +``` + ## Search ```@docs @@ -30,6 +40,13 @@ CurrentModule = Manopt DefaultMeshAdaptiveDirectSearch ``` +as well as the internal functions + +```@docs +Manopt.is_successful(::DefaultMeshAdaptiveDirectSearch) +Manopt.get_candidate(::DefaultMeshAdaptiveDirectSearch) +``` + ## Additional stopping criteria ```@docs diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index c59ee7a8dd..55966fe87f 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -7,10 +7,13 @@ A subtype of this The functor has to fulfil * be callable as `poll!(problem, mesh_size; kwargs...)` and modify the state -as well as +as well as to provide functions -* provide a `is_successful(poll!)` function that indicates whether the last poll was successful in finding a new candidate, -this returns the last successful mesh vector used. +* `is_successful(poll!)` that indicates whether the last poll was successful in finding a new candidate +* `get_basepoint(poll!)` that returns the base point at which the mesh is build +* `get_candidate(poll!)` that returns the last found candidate if the poll was successful. Otherwise the base point is returned +* `get_descent_direction(poll!)` the the vector that points from the base point to the candidate. If the last poll was not successful, the zero vector is returned +* `update_basepoint!(M, poll!, p)` that updates the base point to `p` and all necessary internal data to a new point to build a mesh at The `kwargs...` could include * `scale_mesh=1.0`: to rescale the mesh globally @@ -26,6 +29,12 @@ Should be callable as search!(problem, mesh_size, p, X; kwargs...) where `X` is the last successful poll direction from the tangent space at `p`` if that exists and the zero vector otherwise. + + +Besides that the following functions should be implemented + +* `is_successful(search!)` that indicates whether the last search was successful in finding a new candidate +* `get_candidate(search!)` that returns the last found candidate """ abstract type AbstractMeshSearchFunction end From aff790094d409f8ce6313d6f7bffa14fc735a7ba Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 26 Jan 2025 15:50:41 +0100 Subject: [PATCH 18/29] work on code coverage. --- src/Manopt.jl | 1 + src/plans/mesh_adaptive_plan.jl | 33 +++++++---- src/solvers/mesh_adaptive_direct_search.jl | 3 +- test/plans/test_mesh_adaptive_plan.jl | 59 +++++++++++++++++++ test/runtests.jl | 1 + .../test_mesh_adaptive_direct_search.jl | 24 +++++++- test/test_deprecated.jl | 4 +- 7 files changed, 109 insertions(+), 16 deletions(-) create mode 100644 test/plans/test_mesh_adaptive_plan.jl diff --git a/src/Manopt.jl b/src/Manopt.jl index c2174340aa..711e8d6117 100644 --- a/src/Manopt.jl +++ b/src/Manopt.jl @@ -341,6 +341,7 @@ export AbstractGradientSolverState, InteriorPointNewtonState, LanczosState, LevenbergMarquardtState, + MeshAdaptiveDirectSearchState, NelderMeadState, ParticleSwarmState, PrimalDualSemismoothNewtonState, diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index 55966fe87f..58bf62ca55 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -16,6 +16,7 @@ as well as to provide functions * `update_basepoint!(M, poll!, p)` that updates the base point to `p` and all necessary internal data to a new point to build a mesh at The `kwargs...` could include + * `scale_mesh=1.0`: to rescale the mesh globally * `max_stepsize=Inf`: avoid exceeding a step size beyon this, e.g. injectivity radius. any vector longer than this should be shortened to the provided max stepsize. @@ -76,20 +77,22 @@ with two small modifications: # Fields * `base_point::P`: a point on the manifold, where the mesh is build in the tangent space +* `basis`: a basis of the current tangent space with respect to which the mesh is stored * `candidate::P`: a memory for a new point/candidate +* `mesh`: a vector of tangent vectors storing the mesh. * `random_vector`: a ``d``-dimensional random vector ``b_l``` * `random_index`: a random index ``ι`` -* `mesh`: a vector of tangent vectors storing the mesh. -* `basis`: a basis of the current tangent space with respect to which the mesh is stored -* `X::T` the last successful poll direction stored as a tangent vector. - initialised to the zero vector and reset to the zero vector after moving to a new tangent space. $(_var(:Field, :retraction_method)) $(_var(:Field, :vector_transport_method)) +* `X::T` the last successful poll direction stored as a tangent vector. + initialised to the zero vector and reset to the zero vector after moving to a new tangent space. # Constructor LowerTriangularAdaptivePoll(M, p=rand(M); kwargs...) +## Keyword arguments + * `basis=`[`DefaultOrthonormalBasis`](@extref `ManifoldsBase.DefaultOrthonormalBasis`) $(_var(:Keyword, :retraction_method)) $(_var(:Keyword, :vector_transport_method)) @@ -194,11 +197,16 @@ function update_basepoint!(M, ltap::LowerTriangularAdaptivePoll{P}, p::P) where return ltap end function show(io::IO, ltap::LowerTriangularAdaptivePoll) - s = "LowerTriangularAdaptivePoll on a basis $(ltap.basis), the retraction_method $(ltap.retraction_method), and the vector_transport_method $(ltap.vector_transport_method)" + s = """LowerTriangularAdaptivePoll + with + * basis on the tangent space: $(ltap.basis) + * retraction_method: $(ltap.retraction_method) + * vector_transport_method: $(ltap.vector_transport_method) + """ return print(io, s) end function (ltap::LowerTriangularAdaptivePoll)( - amp::AbstractManoptProblem, mesh_size; scale_mesh=1.0, max_stepsize=inf + amp::AbstractManoptProblem, mesh_size; scale_mesh=1.0, max_stepsize=Inf ) M = get_manifold(amp) n = manifold_dimension(M) @@ -303,7 +311,7 @@ mutable struct DefaultMeshAdaptiveDirectSearch{P,T,RM} <: AbstractMeshSearchFunc retraction_method::RM end function DefaultMeshAdaptiveDirectSearch( - M, p=rand(M); X=zero_vector(M, p), retraction_method=default_retaction_method(M) + M, p=rand(M); X=zero_vector(M, p), retraction_method=default_retraction_method(M) ) return DefaultMeshAdaptiveDirectSearch(p, copy(M, p), X, false, retraction_method) end @@ -324,11 +332,14 @@ function get_candidate(dmads::DefaultMeshAdaptiveDirectSearch) return dmads.p end function show(io::IO, dmads::DefaultMeshAdaptiveDirectSearch) - s = "DefaultMeshAdaptiveDirectSearch using `retraction_method=`$(dmads.retraction_method)" + s = """DefaultMeshAdaptiveDirectSearch + with + * retraction_method: $(dmads.retraction_method) + """ return print(io, s) end function (dmads::DefaultMeshAdaptiveDirectSearch)( - amp::AbstractManoptProblem, mesh_size, p, X; scale_mesh=1.0, max_stepsize=inf + amp::AbstractManoptProblem, mesh_size, p, X; scale_mesh=1.0, max_stepsize=Inf ) M = get_manifold(amp) dmads.X = 4 * mesh_size * scale_mesh * X @@ -415,8 +426,8 @@ function show(io::IO, mads::MeshAdaptiveDirectSearchState) * scale_mesh: $(mads.scale_mesh) * max_stepsize: $(mads.max_stepsize) * poll_size: $(mads.poll_size) - * poll: $(repr(mads.poll)) - * search: $(repr(mads.search)) + * poll:\n $(replace(repr(mads.poll), "\n" => "\n ")[1:end-3]) + * search:\n $(replace(repr(mads.search), "\n" => "\n ")[1:end-3]) ## Stopping criterion $(status_summary(mads.stop)) diff --git a/src/solvers/mesh_adaptive_direct_search.jl b/src/solvers/mesh_adaptive_direct_search.jl index 70377088b1..e4d676b214 100644 --- a/src/solvers/mesh_adaptive_direct_search.jl +++ b/src/solvers/mesh_adaptive_direct_search.jl @@ -64,7 +64,8 @@ function mesh_adaptive_direct_search!( dmco = decorate_objective!(M, mco; kwargs...) dmp = DefaultManoptProblem(M, dmco) madss = MeshAdaptiveDirectSearchState( - M; + M, + p; mesh_basis=mesh_basis, scale_mesh=scale_mesh, max_stepsize=max_stepsize, diff --git a/test/plans/test_mesh_adaptive_plan.jl b/test/plans/test_mesh_adaptive_plan.jl new file mode 100644 index 0000000000..d5bd0de18e --- /dev/null +++ b/test/plans/test_mesh_adaptive_plan.jl @@ -0,0 +1,59 @@ +using ManifoldsBase, Manifolds, Manopt, Test + +@testset "Test Mesh Adaptive Plan" begin + M = ManifoldsBase.DefaultManifold(3) + f(M, p) = norm(p) + mesh_size = 1.0 + cmp = DefaultManoptProblem(M, ManifoldCostObjective(f)) + @testset "Test Poll Accessors" begin + ltap = LowerTriangularAdaptivePoll(M, zeros(3)) + # On init there was not yet a check + @test !Manopt.is_successful(ltap) + @test Manopt.get_descent_direction(ltap) == ltap.X + @test Manopt.get_candidate(ltap) == ltap.candidate + p2 = [2.0, 0.0, 0.0] + @test Manopt.update_basepoint!(M, ltap, p2) === ltap + @test Manopt.get_basepoint(ltap) == p2 + @test startswith(repr(ltap), "LowerTriangularAdaptivePoll\n") + # test call + Random.seed!(42) + ltap(cmp, mesh_size) + # check that this was successful + @test Manopt.is_successful(ltap) + # test call2 scale down! + Random.seed!(42) + ltap(cmp, mesh_size; max_stepsize=1.0) + # check that this was successful as well + @test Manopt.is_successful(ltap) + #... and short enough + @test norm(M, p2, Manopt.get_descent_direction(ltap)) <= 1.0 + end + + @testset "Test Search Accessors" begin + p = ones(3) + dmads = DefaultMeshAdaptiveDirectSearch(M, p) + @test !Manopt.is_successful(dmads) + @test Manopt.get_candidate(dmads) == dmads.p + @test startswith(repr(dmads), "DefaultMeshAdaptiveDirectSearch\n") + X = -ones(3) + # This step would bring us to zero, but we only allow a max step 1.0 + dmads(cmp, 1.0, p, X; max_stepsize=1.0) + # and that should still improve + @test Manopt.is_successful(dmads) + end + + @testset "State Accessors" begin + p = ones(3) + mads = MeshAdaptiveDirectSearchState(M, p) + @test startswith( + repr(mads), "# Solver state for `Manopt.jl`s mesh adaptive direct search\n" + ) + @test get_iterate(mads) == p + @test get_solver_result(mads) == p + end + + @testset "Stopping Criteria" begin + sps = StopWhenPollSizeLess(1.0) + @test startswith(repr(sps), "StopWhenPollSizeLess(1.0)") + end +end diff --git a/test/runtests.jl b/test/runtests.jl index bab88aca10..3dbe0443d3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,6 +15,7 @@ include("utils/example_tasks.jl") include("plans/test_cache.jl") include("plans/test_conjugate_residual_plan.jl") include("plans/test_interior_point_newton_plan.jl") + include("plans/test_mesh_adaptive_plan.jl") include("plans/test_nelder_mead_plan.jl") include("plans/test_nonlinear_least_squares_plan.jl") include("plans/test_gradient_plan.jl") diff --git a/test/solvers/test_mesh_adaptive_direct_search.jl b/test/solvers/test_mesh_adaptive_direct_search.jl index 83ff66b30c..7b1a3c5f64 100644 --- a/test/solvers/test_mesh_adaptive_direct_search.jl +++ b/test/solvers/test_mesh_adaptive_direct_search.jl @@ -1,4 +1,4 @@ -using Manifolds, Manopt, Test, LinearAlgebra +using Manifolds, Manopt, Test, LinearAlgebra, Random @testset "Mesh Adaptive Direct Search" begin # A small spectral procrustes example @@ -8,11 +8,29 @@ using Manifolds, Manopt, Test, LinearAlgebra M = Rotations(2) p0 = [1.0 0.0; 0.0 1.0] f(M, p) = opnorm(B - A * p) - p_s = mesh_adaptive_direct_search( + s = mesh_adaptive_direct_search( M, f, p0; # debug=[:Iteration, :Cost, " ", :poll_size, " ", :mesh_size, " ", :Stop, "\n"], + return_state=true, ) - @test distance(M, p_s, W) < 1e-9 + @test distance(M, get_solver_result(s), W) < 1e-9 + @test startswith(get_reason(s), "The algorithm computed a poll step size") + # + # + # A bit larger example inplace + # A small spectral Procrustes example + A2 = [1.0 0.0 0.0; 1.0 1.0 1.0; 0.0 1.0 2.0; 1.0 1.0 1.0] + α = π / 8 + W2 = [cos(α) -sin(α) 0.0; sin(α) cos(α) 0.0; 0.0 0.0 1.0] + B2 = A2 * W2 + M2 = Rotations(3) + p1 = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] + f2(M, p) = opnorm(B2 - A2 * p) + Random.seed!(42) + # start with a very small mesh size - yields a more exact result + p_s2 = mesh_adaptive_direct_search!(M2, f2, p1; scale_mesh=0.1) + @test isapprox(M, p_s2, p1) + @test distance(M2, p_s2, W2) < 1e-8 end diff --git a/test/test_deprecated.jl b/test/test_deprecated.jl index 3b0b1585b4..b523a88aaf 100644 --- a/test/test_deprecated.jl +++ b/test/test_deprecated.jl @@ -1,3 +1,5 @@ using Manopt, Manifolds, ManifoldsBase, Test -@testset "test deprecated definitions still work" begin end +@testset "test deprecated definitions still work" begin + # after breaking releases this is usually empty. +end From df0f0420e7bca00a56e0155bb41fb7b9739b66e9 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 26 Jan 2025 16:10:29 +0100 Subject: [PATCH 19/29] Cover a final line. --- test/plans/test_mesh_adaptive_plan.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/plans/test_mesh_adaptive_plan.jl b/test/plans/test_mesh_adaptive_plan.jl index d5bd0de18e..39dee947c1 100644 --- a/test/plans/test_mesh_adaptive_plan.jl +++ b/test/plans/test_mesh_adaptive_plan.jl @@ -54,6 +54,7 @@ using ManifoldsBase, Manifolds, Manopt, Test @testset "Stopping Criteria" begin sps = StopWhenPollSizeLess(1.0) + @test get_reason(sps) === "" @test startswith(repr(sps), "StopWhenPollSizeLess(1.0)") end end From a8d47e4c188f334d212cd89ddc895a7d73170de8 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sun, 26 Jan 2025 21:14:24 +0100 Subject: [PATCH 20/29] improve typing and performance a little --- src/plans/mesh_adaptive_plan.jl | 65 ++++++++++++++++----------- test/plans/test_mesh_adaptive_plan.jl | 2 +- 2 files changed, 39 insertions(+), 28 deletions(-) diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index 58bf62ca55..c8a30a444c 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -11,24 +11,27 @@ as well as to provide functions * `is_successful(poll!)` that indicates whether the last poll was successful in finding a new candidate * `get_basepoint(poll!)` that returns the base point at which the mesh is build -* `get_candidate(poll!)` that returns the last found candidate if the poll was successful. Otherwise the base point is returned -* `get_descent_direction(poll!)` the the vector that points from the base point to the candidate. If the last poll was not successful, the zero vector is returned -* `update_basepoint!(M, poll!, p)` that updates the base point to `p` and all necessary internal data to a new point to build a mesh at +* `get_candidate(poll!)` that returns the last found candidate if the poll was successful. + Otherwise the base point is returned +* `get_descent_direction(poll!)` the the vector that points from the base point to the candidate. + If the last poll was not successful, the zero vector is returned +* `update_basepoint!(M, poll!, p)` that updates the base point to `p` and all necessary + internal data to a new point to build a mesh at The `kwargs...` could include * `scale_mesh=1.0`: to rescale the mesh globally -* `max_stepsize=Inf`: avoid exceeding a step size beyon this, e.g. injectivity radius. - any vector longer than this should be shortened to the provided max stepsize. +* `max_stepsize=Inf`: avoid exceeding a step size beyond this value, e.g. injectivity radius. + any vector longer than this should be shortened to the provided maximum step size. """ abstract type AbstractMeshPollFunction end """ AbstractMeshSearchFunction -Should be callable as search!(problem, mesh_size, p, X; kwargs...) +Should be callable as `search!(problem, mesh_size, p, X; kwargs...)` -where `X` is the last successful poll direction from the tangent space at `p`` +where `X` is the last successful poll direction from the tangent space at `p` if that exists and the zero vector otherwise. @@ -123,11 +126,11 @@ mutable struct LowerTriangularAdaptivePoll{ end function LowerTriangularAdaptivePoll( - M, + M::AbstractManifold, p=rand(M); - basis=DefaultOrthonormalBasis(), - retraction_method=default_retraction_method(M), - vector_transport_method=default_vector_transport_method(M), + basis::AbstractBasis=DefaultOrthonormalBasis(), + retraction_method::AbstractRetractionMethod=default_retraction_method(M), + vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method(M), X=zero_vector(M, p), ) d = manifold_dimension(M) @@ -206,7 +209,7 @@ function show(io::IO, ltap::LowerTriangularAdaptivePoll) return print(io, s) end function (ltap::LowerTriangularAdaptivePoll)( - amp::AbstractManoptProblem, mesh_size; scale_mesh=1.0, max_stepsize=Inf + amp::AbstractManoptProblem, mesh_size::Real; scale_mesh::Real=1.0, max_stepsize::Real=Inf ) M = get_manifold(amp) n = manifold_dimension(M) @@ -217,14 +220,19 @@ function (ltap::LowerTriangularAdaptivePoll)( # A random index ι ltap.random_index = rand(1:n) # generate a random b_l vector - for i in 1:n - ltap.random_vector[i] = rand(i == ltap.random_index ? [-2^l, 2^l] : S) - end + rand!(ltap.random_vector, S) + ltap.random_vector[ltap.random_index] = rand((-2^l, 2^l)) end #otherwise we already created ltap.randomvector for this mesh size # Generate L lower triangular, (n-1)x(n-1) in M for i in 1:(n - 1) for j in 1:(n - 1) - ltap.mesh[i, j] = (j > i) ? 0.0 : rand((i == j) ? [-2^l, 2^l] : S) + if i < j + ltap.mesh[i, j] = rand(S) + elseif i == j + ltap.mesh[i, j] = rand((-2^l, 2^l)) + else + ltap.mesh[i, j] = 0 + end end end # Shift to construct n × n matrix B @@ -257,8 +265,9 @@ function (ltap::LowerTriangularAdaptivePoll)( ltap.basis, ) # shorten if necessary - if norm(M, ltap.base_point, ltap.X) > max_stepsize - ltap.X = max_stepsize / norm(M, ltap.base_point, ltap.X) * ltap.X + ltap_X_norm = norm(M, ltap.base_point, ltap.X) + if ltap_X_norm > max_stepsize + ltap.X .*= max_stepsize / ltap_X_norm end retract!(M, ltap.candidate, ltap.base_point, ltap.X, ltap.retraction_method) if get_cost(amp, ltap.candidate) < c @@ -285,13 +294,13 @@ end # Functor - (s::DefaultMeshAdaptiveDirectSearch)(problem, mesh_size, X; scale_mesh=1.0, max_stepsize=inf) + (s::DefaultMeshAdaptiveDirectSearch)(problem, mesh_size::Real, X; scale_mesh::Real=1.0, max_stepsize::Real=inf) # Fields * `q`: a temporary memory for a point on the manifold * `X`: information to perform the search, e.g. the last direction found by poll. -* `last_seach_improved::Bool` indicate whether the last search was succesfull, i.e. improved the cost. +* `last_seach_improved::Bool` indicate whether the last search was successful, i.e. improved the cost. $(_var(:Field, :retraction_method)) # Constructor @@ -342,9 +351,11 @@ function (dmads::DefaultMeshAdaptiveDirectSearch)( amp::AbstractManoptProblem, mesh_size, p, X; scale_mesh=1.0, max_stepsize=Inf ) M = get_manifold(amp) - dmads.X = 4 * mesh_size * scale_mesh * X - if norm(M, p, dmads.X) > max_stepsize - dmads.X = max_stepsize / norm(M, dmads.p, dmads.X) * dmads.X + dmads.X .= (4 * mesh_size * scale_mesh) .* X + + dmads_X_norm = norm(M, p, dmads.X) + if dmads_X_norm > max_stepsize + dmads.X .*= max_stepsize / dmads_X_norm end retract!(M, dmads.q, p, dmads.X, dmads.retraction_method) dmads.last_search_improved = get_cost(amp, dmads.q) < get_cost(amp, p) @@ -369,7 +380,7 @@ $(_var(:Field, :stopping_criterion, "stop")) * `search::`[`AbstractMeshSearchFunction`}(@ref) a search step (functor) to perform """ -mutable struct MeshAdaptiveDirectSearchState{P,F<:Real,PT,ST,SC<:StoppingCriterion} <: +mutable struct MeshAdaptiveDirectSearchState{P,F<:Real,PT<:AbstractMeshPollFunction,ST<:AbstractMeshSearchFunction,SC<:StoppingCriterion} <: AbstractManoptSolverState p::P mesh_size::F @@ -385,10 +396,10 @@ function MeshAdaptiveDirectSearchState( p::P=rand(M); mesh_basis::B=DefaultOrthonormalBasis(), scale_mesh::F=injectivity_radius(M) / 2, - max_stepsize=injectivity_radius(M), + max_stepsize::F=injectivity_radius(M), stopping_criterion::SC=StopAfterIteration(500) | StopWhenPollSizeLess(1e-7), - retraction_method=default_retraction_method(M, typeof(p)), - vector_transport_method=default_vector_transport_method(M, typeof(p)), + retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), + vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method(M, typeof(p)), poll::PT=LowerTriangularAdaptivePoll( M, copy(M, p); diff --git a/test/plans/test_mesh_adaptive_plan.jl b/test/plans/test_mesh_adaptive_plan.jl index 39dee947c1..1ad7e00e7a 100644 --- a/test/plans/test_mesh_adaptive_plan.jl +++ b/test/plans/test_mesh_adaptive_plan.jl @@ -1,4 +1,4 @@ -using ManifoldsBase, Manifolds, Manopt, Test +using ManifoldsBase, Manifolds, Manopt, Test, Random @testset "Test Mesh Adaptive Plan" begin M = ManifoldsBase.DefaultManifold(3) From 1d6454eae1cb09c321de567017b92bcf8aa84f8f Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sun, 26 Jan 2025 21:20:39 +0100 Subject: [PATCH 21/29] formatting --- src/plans/mesh_adaptive_plan.jl | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index c8a30a444c..a012b7e818 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -130,7 +130,9 @@ function LowerTriangularAdaptivePoll( p=rand(M); basis::AbstractBasis=DefaultOrthonormalBasis(), retraction_method::AbstractRetractionMethod=default_retraction_method(M), - vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method(M), + vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method( + M + ), X=zero_vector(M, p), ) d = manifold_dimension(M) @@ -209,7 +211,10 @@ function show(io::IO, ltap::LowerTriangularAdaptivePoll) return print(io, s) end function (ltap::LowerTriangularAdaptivePoll)( - amp::AbstractManoptProblem, mesh_size::Real; scale_mesh::Real=1.0, max_stepsize::Real=Inf + amp::AbstractManoptProblem, + mesh_size::Real; + scale_mesh::Real=1.0, + max_stepsize::Real=Inf, ) M = get_manifold(amp) n = manifold_dimension(M) @@ -380,8 +385,13 @@ $(_var(:Field, :stopping_criterion, "stop")) * `search::`[`AbstractMeshSearchFunction`}(@ref) a search step (functor) to perform """ -mutable struct MeshAdaptiveDirectSearchState{P,F<:Real,PT<:AbstractMeshPollFunction,ST<:AbstractMeshSearchFunction,SC<:StoppingCriterion} <: - AbstractManoptSolverState +mutable struct MeshAdaptiveDirectSearchState{ + P, + F<:Real, + PT<:AbstractMeshPollFunction, + ST<:AbstractMeshSearchFunction, + SC<:StoppingCriterion, +} <: AbstractManoptSolverState p::P mesh_size::F scale_mesh::F @@ -399,7 +409,9 @@ function MeshAdaptiveDirectSearchState( max_stepsize::F=injectivity_radius(M), stopping_criterion::SC=StopAfterIteration(500) | StopWhenPollSizeLess(1e-7), retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), - vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method(M, typeof(p)), + vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method( + M, typeof(p) + ), poll::PT=LowerTriangularAdaptivePoll( M, copy(M, p); From 83da62efbf993f75e62d2d43f294da9821e69dff Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Mon, 27 Jan 2025 19:57:31 +0100 Subject: [PATCH 22/29] fix some typos, add some types --- src/plans/mesh_adaptive_plan.jl | 23 +++++++++++++++------- src/solvers/mesh_adaptive_direct_search.jl | 12 +++++------ 2 files changed, 22 insertions(+), 13 deletions(-) diff --git a/src/plans/mesh_adaptive_plan.jl b/src/plans/mesh_adaptive_plan.jl index a012b7e818..0d43d7e74a 100644 --- a/src/plans/mesh_adaptive_plan.jl +++ b/src/plans/mesh_adaptive_plan.jl @@ -227,7 +227,7 @@ function (ltap::LowerTriangularAdaptivePoll)( # generate a random b_l vector rand!(ltap.random_vector, S) ltap.random_vector[ltap.random_index] = rand((-2^l, 2^l)) - end #otherwise we already created ltap.randomvector for this mesh size + end #otherwise we already created ltap.random_vector for this mesh size # Generate L lower triangular, (n-1)x(n-1) in M for i in 1:(n - 1) for j in 1:(n - 1) @@ -305,7 +305,7 @@ end * `q`: a temporary memory for a point on the manifold * `X`: information to perform the search, e.g. the last direction found by poll. -* `last_seach_improved::Bool` indicate whether the last search was successful, i.e. improved the cost. +* `last_search_improved::Bool` indicate whether the last search was successful, i.e. improved the cost. $(_var(:Field, :retraction_method)) # Constructor @@ -314,10 +314,11 @@ $(_var(:Field, :retraction_method)) ## Keyword arguments -* `X::T=zero_vector(M,p) +* `X::T=zero_vector(M, p) $(_var(:Keyword, :retraction_method)) """ -mutable struct DefaultMeshAdaptiveDirectSearch{P,T,RM} <: AbstractMeshSearchFunction +mutable struct DefaultMeshAdaptiveDirectSearch{P,T,RM<:AbstractRetractionMethod} <: + AbstractMeshSearchFunction p::P q::P X::T @@ -325,14 +326,17 @@ mutable struct DefaultMeshAdaptiveDirectSearch{P,T,RM} <: AbstractMeshSearchFunc retraction_method::RM end function DefaultMeshAdaptiveDirectSearch( - M, p=rand(M); X=zero_vector(M, p), retraction_method=default_retraction_method(M) + M::AbstractManifold, + p=rand(M); + X=zero_vector(M, p), + retraction_method::AbstractRetractionMethod=default_retraction_method(M), ) return DefaultMeshAdaptiveDirectSearch(p, copy(M, p), X, false, retraction_method) end """ is_successful(dmads::DefaultMeshAdaptiveDirectSearch) -Return whether the last [`DefaultMeshAdaptiveDirectSearch`](@ref) was succesful. +Return whether the last [`DefaultMeshAdaptiveDirectSearch`](@ref) was successful. """ function is_successful(dmads::DefaultMeshAdaptiveDirectSearch) return dmads.last_search_improved @@ -353,7 +357,12 @@ function show(io::IO, dmads::DefaultMeshAdaptiveDirectSearch) return print(io, s) end function (dmads::DefaultMeshAdaptiveDirectSearch)( - amp::AbstractManoptProblem, mesh_size, p, X; scale_mesh=1.0, max_stepsize=Inf + amp::AbstractManoptProblem, + mesh_size::Real, + p, + X; + scale_mesh::Real=1.0, + max_stepsize::Real=Inf, ) M = get_manifold(amp) dmads.X .= (4 * mesh_size * scale_mesh) .* X diff --git a/src/solvers/mesh_adaptive_direct_search.jl b/src/solvers/mesh_adaptive_direct_search.jl index e4d676b214..1d831f7d03 100644 --- a/src/solvers/mesh_adaptive_direct_search.jl +++ b/src/solvers/mesh_adaptive_direct_search.jl @@ -41,8 +41,8 @@ function mesh_adaptive_direct_search!( mco::AbstractManifoldCostObjective, p; mesh_basis::B=DefaultOrthonormalBasis(), - scale_mesh=injectivity_radius(M) / 2, - max_stepsize=injectivity_radius(M), + scale_mesh::Real=injectivity_radius(M) / 2, + max_stepsize::Real=injectivity_radius(M), stopping_criterion::StoppingCriterion=StopAfterIteration(500) | StopWhenPollSizeLess(1e-10), retraction_method::AbstractRetractionMethod=default_retraction_method(M, eltype(p)), @@ -97,7 +97,7 @@ end function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearchState, k) M = get_manifold(amp) n = manifold_dimension(M) - # search if the last poll or last search was sucessful + # search if the last poll or last search was successful if is_successful(madss.search) || is_successful(madss.poll) madss.search( amp, @@ -108,11 +108,11 @@ function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearc max_stepsize=madss.max_stepsize, ) end - # For succesful search, copy over iterate - skip poll, but update base + # For successful search, copy over iterate - skip poll, but update base if is_successful(madss.search) copyto!(M, madss.p, get_candidate(madss.search)) update_basepoint!(M, madss.poll, madss.p) - else #search was not sucessful: poll + else #search was not successful: poll update_basepoint!(M, madss.poll, madss.p) madss.poll( amp, @@ -120,7 +120,7 @@ function step_solver!(amp::AbstractManoptProblem, madss::MeshAdaptiveDirectSearc scale_mesh=madss.scale_mesh, max_stepsize=madss.max_stepsize, ) - # For succesfull poll, copy over iterate + # For successful poll, copy over iterate if is_successful(madss.poll) copyto!(M, madss.p, get_candidate(madss.poll)) end From 0c6322b8780af12a4347a03b833c6a04dcea555f Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Tue, 4 Feb 2025 12:15:55 +0100 Subject: [PATCH 23/29] A bit of work on typos. --- .vale.ini | 3 ++- Changelog.md | 2 +- Readme.md | 2 +- docs/make.jl | 2 +- docs/src/about.md | 4 ++-- docs/src/index.md | 4 ++-- docs/src/solvers/index.md | 2 +- docs/styles/config/vocabularies/Manopt/accept.txt | 1 + src/plans/debug.jl | 4 ++-- src/plans/record.jl | 6 +++--- src/solvers/augmented_Lagrangian_method.jl | 2 +- src/solvers/convex_bundle_method.jl | 2 +- src/solvers/proximal_bundle_method.jl | 2 +- test/solvers/test_difference_of_convex.jl | 4 ++-- tutorials/CountAndCache.qmd | 4 ++-- tutorials/EmbeddingObjectives.qmd | 2 +- tutorials/HowToDebug.qmd | 8 ++++---- tutorials/HowToRecord.qmd | 8 ++++---- tutorials/ImplementASolver.qmd | 4 ++-- tutorials/ImplementOwnManifold.qmd | 4 ++-- tutorials/InplaceGradient.qmd | 2 +- tutorials/Optimize.qmd | 4 ++-- 22 files changed, 39 insertions(+), 37 deletions(-) diff --git a/.vale.ini b/.vale.ini index ff2bb0bce6..3f4c7cfd64 100644 --- a/.vale.ini +++ b/.vale.ini @@ -7,6 +7,7 @@ Packages = Google [formats] # code blocks with Julia in Markdown do not yet work well qmd = md +jl = md [docs/src/*.md] BasedOnStyles = Vale, Google @@ -39,7 +40,7 @@ TokenIgnores = \$(.+)\$,\[.+?\]\(@(ref|id|cite).+?\),`.+`,``.*``,\s{4}.+\n Google.Units = false #wto ignore formats= for now. TokenIgnores = \$(.+)\$,\[.+?\]\(@(ref|id|cite).+?\),`.+`,``.*``,\s{4}.+\n -[tutorials/*.md] ; actually .qmd for the first, second autogenerated +[tutorials/*.qmd] ; actually .qmd for the first, second autogenerated BasedOnStyles = Vale, Google ; ignore (1) math (2) ref and cite keys (3) code in docs (4) math in docs (5,6) indented blocks TokenIgnores = (\$+[^\n$]+\$+) diff --git a/Changelog.md b/Changelog.md index 41465f2ec4..7ba5eba06c 100644 --- a/Changelog.md +++ b/Changelog.md @@ -128,7 +128,7 @@ In general we introduce a few factories, that avoid having to pass the manifold * the previous `stabilize=true` is now set with `(project!)=embed_project!` in general, and if the manifold is represented by points in the embedding, like the sphere, `(project!)=project!` suffices * the new default is `(project!)=copyto!`, so by default no projection/stabilization is performed. -* the positional argument `p` (usually the last or the third to last if subsolvers existed) has been moved to a keyword argument `p=` in all State constructors +* the positional argument `p` (usually the last or the third to last if sub solvers existed) has been moved to a keyword argument `p=` in all State constructors * in `NelderMeadState` the `population` moved from positional to keyword argument as well, * the way to initialise sub solvers in the solver states has been unified In the new variant * the `sub_problem` is always a positional argument; namely the last one diff --git a/Readme.md b/Readme.md index 36c3733ca4..fd3684e62a 100644 --- a/Readme.md +++ b/Readme.md @@ -36,7 +36,7 @@ In Julia you can get started by just typing using Pkg; Pkg.add("Manopt"); ``` -and then checkout the [Get started: optimize!](https://manoptjl.org/stable/tutorials/Optimize/) tutorial. +and then checkout the [🏔️ Get started with Manopt.jl](https://manoptjl.org/stable/tutorials/Optimize/) tutorial. ## Related packages diff --git a/docs/make.jl b/docs/make.jl index 20eb70d596..689fce6ecc 100755 --- a/docs/make.jl +++ b/docs/make.jl @@ -35,7 +35,7 @@ tutorials_in_menu = !("--exclude-tutorials" ∈ ARGS) # (a) setup the tutorials menu – check whether all files exist tutorials_menu = "How to..." => [ - "🏔️ Get started: optimize." => "tutorials/Optimize.md", + "🏔️ Get started with Manopt.jl." => "tutorials/Optimize.md", "Speedup using in-place computations" => "tutorials/InplaceGradient.md", "Use automatic differentiation" => "tutorials/AutomaticDifferentiation.md", "Define objectives in the embedding" => "tutorials/EmbeddingObjectives.md", diff --git a/docs/src/about.md b/docs/src/about.md index c82bddb48a..79e2b25339 100644 --- a/docs/src/about.md +++ b/docs/src/about.md @@ -29,8 +29,8 @@ to clone/fork the repository or open an issue. * [ExponentialFamilyProjection.jl](https://github.com/ReactiveBayes/ExponentialFamilyProjection.jl) package uses `Manopt.jl` to project arbitrary functions onto the closest exponential family distributions. The package also integrates with [`RxInfer.jl`](https://github.com/ReactiveBayes/RxInfer.jl) to enable Bayesian inference in a larger set of probabilistic models. * [Caesar.jl](https://github.com/JuliaRobotics/Caesar.jl) within non-Gaussian factor graph inference algorithms -Is a package missing? [Open an issue](https://github.com/JuliaManifolds/Manopt.jl/issues/new)! -It would be great to collect anything and anyone using Manopt.jl +If you are missing a package, that uses `Manopt.jl`, please [open an issue](https://github.com/JuliaManifolds/Manopt.jl/issues/new). +It would be great to collect anything and anyone using Manopt.jl in this list. ## Further packages diff --git a/docs/src/index.md b/docs/src/index.md index 1eb56a1430..93e9496344 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -20,7 +20,7 @@ or in other words: find the point ``p`` on the manifold, where ``f`` reaches its It belongs to the “Manopt family”, which includes [Manopt](https://manopt.org) (Matlab) and [pymanopt.org](https://www.pymanopt.org/) (Python). If you want to delve right into `Manopt.jl` read the -[🏔️ Get started: optimize.](tutorials/Optimize.md) tutorial. +[🏔️ Get started with Manopt.jl.](tutorials/Optimize.md) tutorial. `Manopt.jl` makes it easy to use an algorithm for your favourite manifold as well as a manifold for your favourite algorithm. It already provides @@ -94,7 +94,7 @@ The notation in the documentation aims to follow the same [notation](https://jul ### Visualization To visualize and interpret results, `Manopt.jl` aims to provide both easy plot functions as well as [exports](helpers/exports.md). Furthermore a system to get [debug](plans/debug.md) during the iterations of an algorithms as well as [record](plans/record.md) capabilities, for example to record a specified tuple of values per iteration, most prominently [`RecordCost`](@ref) and -[`RecordIterate`](@ref). Take a look at the [🏔️ Get started: optimize.](tutorials/Optimize.md) tutorial on how to easily activate this. +[`RecordIterate`](@ref). Take a look at the [🏔️ Get started with Manopt.jl.](tutorials/Optimize.md) tutorial on how to easily activate this. ## Literature diff --git a/docs/src/solvers/index.md b/docs/src/solvers/index.md index bd6900e6c7..f59b3faeb6 100644 --- a/docs/src/solvers/index.md +++ b/docs/src/solvers/index.md @@ -203,7 +203,7 @@ also use the third (lowest level) and just call solve!(problem, state) ``` -### Closed-form subsolvers +### Closed-form sub solvers If a subsolver solution is available in closed form, `ClosedFormSubSolverState` is used to indicate that. diff --git a/docs/styles/config/vocabularies/Manopt/accept.txt b/docs/styles/config/vocabularies/Manopt/accept.txt index 75246fb123..afcf45d687 100644 --- a/docs/styles/config/vocabularies/Manopt/accept.txt +++ b/docs/styles/config/vocabularies/Manopt/accept.txt @@ -23,6 +23,7 @@ Cartis canonicalization canonicalized Constantin +[Cc]ubics Dai deactivatable Diepeveen diff --git a/src/plans/debug.jl b/src/plans/debug.jl index 298eb1d2b2..1d205580d3 100644 --- a/src/plans/debug.jl +++ b/src/plans/debug.jl @@ -158,7 +158,7 @@ Whether internal variables are updates is determined by `always_update`. This method does not perform any print itself but relies on it's children's print. -It also sets the subsolvers active parameter, see |`DebugWhenActive`}(#ref). +It also sets the sub solvers active parameter, see |`DebugWhenActive`}(#ref). Here, the `activattion_offset` can be used to specify whether it refers to _this_ iteration, the `i`th, when this call is _before_ the iteration, then the offset should be 0, for the _next_ iteration, that is if this is called _after_ an iteration, it has to be set to 1. @@ -185,7 +185,7 @@ function (d::DebugEvery)(p::AbstractManoptProblem, st::AbstractManoptSolverState elseif d.always_update d.debug(p, st, -1) end - # set activity for this iterate in subsolvers + # set activity for this iterate in sub solvers set_parameter!( st, :SubState, diff --git a/src/plans/record.jl b/src/plans/record.jl index a8ccabdd04..585b0a234b 100644 --- a/src/plans/record.jl +++ b/src/plans/record.jl @@ -243,7 +243,7 @@ function (re::RecordEvery)( elseif re.always_update re.record(amp, ams, 0) end - # Set activity to activate or deactivate subsolvers + # Set activity to activate or deactivate sub solvers # note that since recording is happening at the end # sets activity for the _next_ iteration set_parameter!( @@ -390,7 +390,7 @@ getindex(r::RecordGroup, i) = get_record(r, i) @doc raw""" RecordSubsolver <: RecordAction -Record the current subsolvers recording, by calling [`get_record`](@ref) +Record the current sub solvers recording, by calling [`get_record`](@ref) on the sub state with # Fields @@ -428,7 +428,7 @@ status_summary(::RecordSubsolver) = ":Subsolver" record action that only records if the `active` boolean is set to true. This can be set from outside and is for example triggered by |`RecordEvery`](@ref) on recordings of the subsolver. -While this is for subsolvers maybe not completely necessary, recording values that +While this is for sub solvers maybe not completely necessary, recording values that are never accessible, is not that useful. # Fields diff --git a/src/solvers/augmented_Lagrangian_method.jl b/src/solvers/augmented_Lagrangian_method.jl index d26006f916..76e64a3596 100644 --- a/src/solvers/augmented_Lagrangian_method.jl +++ b/src/solvers/augmented_Lagrangian_method.jl @@ -473,7 +473,7 @@ function augmented_Lagrangian_method!( ), sub_problem::AbstractManoptProblem=DefaultManoptProblem( M, - # pass down objective type to subsolvers + # pass down objective type to sub solvers decorate_objective!( M, ManifoldGradientObjective(sub_cost, sub_grad; evaluation=evaluation); diff --git a/src/solvers/convex_bundle_method.jl b/src/solvers/convex_bundle_method.jl index 4a19e036f2..99ebf64373 100644 --- a/src/solvers/convex_bundle_method.jl +++ b/src/solvers/convex_bundle_method.jl @@ -474,7 +474,7 @@ end # # -# Dispatching on different types of subsolvers +# Dispatching on different types of sub solvers # (a) closed form allocating function _convex_bundle_subsolver!( M, bms::ConvexBundleMethodState{P,T,F,ClosedFormSubSolverState{AllocatingEvaluation}} diff --git a/src/solvers/proximal_bundle_method.jl b/src/solvers/proximal_bundle_method.jl index 583705ce6c..10773c916c 100644 --- a/src/solvers/proximal_bundle_method.jl +++ b/src/solvers/proximal_bundle_method.jl @@ -432,7 +432,7 @@ get_solver_result(pbms::ProximalBundleMethodState) = pbms.p_last_serious # # -# Dispatching on different types of subsolvers +# Dispatching on different types of sub solvers # (a) closed form allocating function _proximal_bundle_subsolver!( M, pbms::ProximalBundleMethodState{P,T,F,ClosedFormSubSolverState{AllocatingEvaluation}} diff --git a/test/solvers/test_difference_of_convex.jl b/test/solvers/test_difference_of_convex.jl index 99e6231516..f405552d36 100644 --- a/test/solvers/test_difference_of_convex.jl +++ b/test/solvers/test_difference_of_convex.jl @@ -192,7 +192,7 @@ import Manifolds: inner ) end @testset "Running the closed form solution solvers" begin - # make them a bit by providing subsolvers as functions + # make them a bit by providing sub solvers as functions function dca_sub(M, p, X) q = copy(M, p) lin_s = LinearizedDCCost(g, copy(M, p), copy(M, p, X)) @@ -216,7 +216,7 @@ import Manifolds: inner @test isapprox(M, p11, p12) @test f(M, p11) ≈ 0.0 atol = 1e-15 - # fake them a bit by providing subsolvers as functions + # fake them a bit by providing sub solvers as functions function prox_g(M, λ, p) q = copy(M, p) prox = ProximalDCCost(g, copy(M, p), λ) diff --git a/tutorials/CountAndCache.qmd b/tutorials/CountAndCache.qmd index 4e80cf5cbf..ea73d76468 100644 --- a/tutorials/CountAndCache.qmd +++ b/tutorials/CountAndCache.qmd @@ -5,7 +5,7 @@ author: Ronny Bergmann In this tutorial, we want to investigate the caching and counting (statistics) features of [Manopt.jl](https://manoptjl.org). We reuse the optimization tasks from the -introductory tutorial [Get started: optimize!](https://manoptjl.org/stable/tutorials/Optimize!.html). +introductory tutorial [🏔️ Get started with Manopt.jl](https://manoptjl.org/stable/tutorials/Optimize!.html). ## Introduction @@ -67,7 +67,7 @@ using ManifoldDiff: grad_distance ## Counting -We first define our task, the Riemannian Center of Mass from the [Get started: optimize!](https://manoptjl.org/stable/tutorials/Optimize!.html) tutorial. +We first define our task, the Riemannian Center of Mass from the [🏔️ Get started with Manopt.jl](https://manoptjl.org/stable/tutorials/Optimize!.html) tutorial. ```{julia} n = 100 diff --git a/tutorials/EmbeddingObjectives.qmd b/tutorials/EmbeddingObjectives.qmd index a3b7112a2e..cfdcdaee47 100644 --- a/tutorials/EmbeddingObjectives.qmd +++ b/tutorials/EmbeddingObjectives.qmd @@ -181,7 +181,7 @@ distance(M, q1, q2) This conversion also works for the gradients of constraints, and is passed down to -subsolvers by default when these are created using the Euclidean objective $f$, $\nabla f$ and $\nabla^2 f$. +sub solvers by default when these are created using the Euclidean objective $f$, $\nabla f$ and $\nabla^2 f$. ## Summary diff --git a/tutorials/HowToDebug.qmd b/tutorials/HowToDebug.qmd index 9dcc6b1297..cfc22609b1 100644 --- a/tutorials/HowToDebug.qmd +++ b/tutorials/HowToDebug.qmd @@ -106,7 +106,7 @@ Note that the number (`25`) yields that all but `:Start` and `:Stop` are only di ## Subsolver debug -Subsolvers have a `sub_kwargs` keyword, such that you can pass keywords to the sub solver as well. This works well if you do not plan to change the subsolver. If you do you can wrap your own `solver_state=` argument in a [`decorate_state!`](@ref) and pass a `debug=` password to this function call. +Sub solvers have a `sub_kwargs` keyword, such that you can pass keywords to the sub solver as well. This works well if you do not plan to change the subsolver. If you do you can wrap your own `solver_state=` argument in a [`decorate_state!`](@ref) and pass a `debug=` password to this function call. Keywords in a keyword have to be passed as pairs (`:debug => [...]`). For most debugs, there further exists a longer form to specify the format to print. @@ -123,10 +123,10 @@ p3 = exact_penalty_method( ); ``` -The different lengths of the dotted lines come from the fact that ---at least in the beginning--- the subsolver performs a few steps and each subsolvers step prints a dot. +The different lengths of the dotted lines come from the fact that ---at least in the beginning--- the subsolver performs a few steps and each sub solvers step prints a dot. -For this issue, there is the next symbol (similar to the `:Stop`) to indicate that a debug set is a subsolver set `:WhenActive`, which introduces a [`DebugWhenActive`](@ref) that is only activated when the outer debug is actually active, or inother words [`DebugEvery`](@ref) is active itself. -Furthermore, we want to print the iteration number _before_ printing the subsolvers steps, so we put this into a `Pair`, but we can leave the remaining ones as single +For this issue, there is the next symbol (similar to the `:Stop`) to indicate that a debug set is a subsolver set `:WhenActive`, which introduces a [`DebugWhenActive`](@ref) that is only activated when the outer debug is actually active, or another words [`DebugEvery`](@ref) is active itself. +Furthermore, we want to print the iteration number _before_ printing the sub solvers steps, so we put this into a `Pair`, but we can leave the remaining ones as single entries. Finally we also prefix `:Stop` with `" | "` and print the iteration number at the time we stop. We get diff --git a/tutorials/HowToRecord.qmd b/tutorials/HowToRecord.qmd index c396dddf34..774af3de2f 100644 --- a/tutorials/HowToRecord.qmd +++ b/tutorials/HowToRecord.qmd @@ -15,7 +15,7 @@ Several predefined recordings exist, for example [`RecordCost`](https://manoptjl For fields of the `State` the recording can also be done [`RecordEntry`](https://manoptjl.org/stable/plans/record/#Manopt.RecordEvery). For other recordings, for example more advanced computations before storing a value, an own `RecordAction` can be defined. -We illustrate these using the gradient descent from the [Get started: optimize!](https://manoptjl.org/stable/tutorials/Optimize!.html) tutorial. +We illustrate these using the gradient descent from the [Get started: optimize](https://manoptjl.org/stable/tutorials/Optimize!.html) tutorial. Here the focus is put on ways to investigate the behaviour during iterations by using Recording techniques. @@ -204,9 +204,9 @@ grad_g(M, p) = [project(M, p, mI[:, i]) for i in 1:d] p0 = project(M2, [ones(2)..., zeros(d - 3)..., 0.1]) ``` -We directly start with recording the subsolvers Iteration. +We directly start with recording the sub solvers Iteration. We can specify what to record in the subsolver using the `sub_kwargs` -keyword argument with a `Symbol => value` pair. Here we specify to record the iteration and the cost in every subsolvers step. +keyword argument with a `Symbol => value` pair. Here we specify to record the iteration and the cost in every sub solvers step. Furthermore, we have to “collect” this recording after every sub solver run. This is done with the `:Subsolver` keyword in the main `record=` keyword. @@ -273,7 +273,7 @@ s3 = exact_penalty_method( ); ``` -Then the following displays also the reasons why each of the recorded subsolvers stopped and the corresponding cost +Then the following displays also the reasons why each of the recorded sub solvers stopped and the corresponding cost ```{julia} get_record(s3) diff --git a/tutorials/ImplementASolver.qmd b/tutorials/ImplementASolver.qmd index 51e0840cb5..cf65ca1040 100644 --- a/tutorials/ImplementASolver.qmd +++ b/tutorials/ImplementASolver.qmd @@ -4,7 +4,7 @@ author: "Ronny Bergmann" --- When you have used a few solvers from `Manopt.jl` for example like in the opening -tutorial [Get started: optimize!](https://manoptjl.org/stable/tutorials/Optimize!.html) +tutorial [Get started: optimize](https://manoptjl.org/stable/tutorials/Optimize!.html) you might come to the idea of implementing a solver yourself. After a short introduction of the algorithm we aim to implement, @@ -204,7 +204,7 @@ In practice, however, it is preferable to cache intermediate values like cost of Now we can just run the solver already. We take the same example as for the other tutorials -We first define our task, the Riemannian Center of Mass from the [Get started: optimize!](https://manoptjl.org/stable/tutorials/Optimize!.html) tutorial. +We first define our task, the Riemannian Center of Mass from the [Get started: optimize](https://manoptjl.org/stable/tutorials/Optimize!.html) tutorial. ```{julia} #| output: false diff --git a/tutorials/ImplementOwnManifold.qmd b/tutorials/ImplementOwnManifold.qmd index 93fc77ae8e..1373b6f49f 100644 --- a/tutorials/ImplementOwnManifold.qmd +++ b/tutorials/ImplementOwnManifold.qmd @@ -10,7 +10,7 @@ CurrentModule = Manopt ```` When you have used a few solvers from [`Manopt.jl`](https://manoptjl.org/) for example like in the opening -tutorial [🏔️ Get started: optimize!](https://manoptjl.org/stable/tutorials/Optimize!.html) +tutorial [🏔️ Get started with Manopt.jl](https://manoptjl.org/stable/tutorials/Optimize!.html) and also familiarized yourself with how to work with manifolds in general at [🚀 Get Started with `Manifolds.jl`](https://juliamanifolds.github.io/Manifolds.jl/stable/tutorials/getstarted.html), you might come across the point that you want to @@ -80,7 +80,7 @@ struct ScaledSphere <: AbstractManifold{ℝ} end ``` -We would like to compute a mean and/or median similar to [🏔️ Get started: optimize!](https://manoptjl.org/stable/tutorials/Optimize!.html). +We would like to compute a mean and/or median similar to [🏔️ Get started with Manopt.jl!](https://manoptjl.org/stable/tutorials/Optimize!.html). For given a set of points $q_1,\ldots,q_n$ we want to compute [Karcher:1977](@cite) ```math diff --git a/tutorials/InplaceGradient.qmd b/tutorials/InplaceGradient.qmd index cb63b1dc4a..a2b3174566 100644 --- a/tutorials/InplaceGradient.qmd +++ b/tutorials/InplaceGradient.qmd @@ -7,7 +7,7 @@ When it comes to time critical operations, a main ingredient in Julia is given b mutating functions, that is those that compute in place without additional memory allocations. In the following, we illustrate how to do this with `Manopt.jl`. -Let's start with the same function as in [Get started: optimize!](https://manoptjl.org/stable/tutorials/Optimize!.html) +Let's start with the same function as in [🏔️ Get started with Manopt.jl](https://manoptjl.org/stable/tutorials/Optimize!.html) and compute the mean of some points, only that here we use the sphere $\mathbb S^{30}$ and $n=800$ points. diff --git a/tutorials/Optimize.qmd b/tutorials/Optimize.qmd index 3a155ed82f..1e8b3bda74 100644 --- a/tutorials/Optimize.qmd +++ b/tutorials/Optimize.qmd @@ -1,5 +1,5 @@ --- -title: "🏔️ Get started: optimize." +title: "🏔️ Get started with Manopt.jl." author: Ronny Bergmann --- @@ -20,7 +20,7 @@ This can also be written as ``` where the aim is to compute the minimizer $p^*$ numerically. -As an example, consider the generalisation of the [(arithemtic) mean](https://en.wikipedia.org/wiki/Arithmetic_mean). +As an example, consider the generalisation of the [(arithmetic) mean](https://en.wikipedia.org/wiki/Arithmetic_mean). In the Euclidean case with $d∈\mathbb N$, that is for $n∈\mathbb N$ data points $y_1,\ldots,y_n ∈ ℝ^d$ the mean ```math From 5e7f23241802e3abdbabab5e2eb706beaa3bec03 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Tue, 4 Feb 2025 13:01:09 +0100 Subject: [PATCH 24/29] Update metadata. --- .vale.ini | 13 +++++++++---- .zenodo.json | 5 +++++ docs/src/solvers/index.md | 2 +- docs/styles/config/vocabularies/Manopt/accept.txt | 4 ++++ 4 files changed, 19 insertions(+), 5 deletions(-) diff --git a/.vale.ini b/.vale.ini index 3f4c7cfd64..054a57d418 100644 --- a/.vale.ini +++ b/.vale.ini @@ -13,7 +13,10 @@ jl = md BasedOnStyles = Vale, Google [docs/src/contributing.md] -BasedOnStyles = +BasedOnStyles = Vale, Google +Google.Will = false ; given format and really with intend a _will_ +Google.Headings = false ; some might jeally ahabe [] in their headers +Google.FirstPerson = false ; we pose a few contribution points as first-person questions [Changelog.md, CONTRIBUTING.md] BasedOnStyles = Vale, Google @@ -46,6 +49,8 @@ BasedOnStyles = Vale, Google TokenIgnores = (\$+[^\n$]+\$+) Google.We = false # For tutorials we want to address the user directly. -[docs/src/tutorials/*.md] - ; ignore since they are derived files -BasedOnStyles = +[docs/src/tutorials/*.md] ; actually .qmd for the first, second autogenerated +BasedOnStyles = Vale, Google +; ignore (1) math (2) ref and cite keys (3) code in docs (4) math in docs (5,6) indented blocks +TokenIgnores = (\$+[^\n$]+\$+) +Google.We = false # For tutorials we want to address the user directly. \ No newline at end of file diff --git a/.zenodo.json b/.zenodo.json index d4c2bec6ee..7e5fd36fe7 100644 --- a/.zenodo.json +++ b/.zenodo.json @@ -25,6 +25,11 @@ "name": "Riemer, Tom-Christian", "type": "ProjectMember" }, + { + "affiliation": "NTNU Trondheim", + "name": "Oddsen, Sander Engen", + "type": "ProjectMember" + }, { "name": "Schilly, Harald", "type": "Other" diff --git a/docs/src/solvers/index.md b/docs/src/solvers/index.md index f59b3faeb6..ba5746cbd6 100644 --- a/docs/src/solvers/index.md +++ b/docs/src/solvers/index.md @@ -99,7 +99,7 @@ For these you can use * [Steihaug-Toint Truncated Conjugate-Gradient Method](truncated_conjugate_gradient_descent.md) a solver for a constrained problem defined on a tangent space. -## Alphabetical list List of algorithms +## Alphabetical list of algorithms | Solver | Function | State | |:---------|:----------------|:---------| diff --git a/docs/styles/config/vocabularies/Manopt/accept.txt b/docs/styles/config/vocabularies/Manopt/accept.txt index afcf45d687..73e5f1fef9 100644 --- a/docs/styles/config/vocabularies/Manopt/accept.txt +++ b/docs/styles/config/vocabularies/Manopt/accept.txt @@ -77,9 +77,11 @@ Munkvold [Mm]ead [Nn]elder Nesterov +Nesterovs Newton nonmonotone nonpositive +[Nn]onsmooth [Pp]arametrising Parametrising [Pp]ock @@ -111,9 +113,11 @@ Stephansen Stokkenes [Ss]ubdifferential [Ss]ubgradient +[Ss]ubgradients subsampled [Ss]ubsolver summand +summands superlinear supertype th From 577dfb51573a9705122984ab9c8fc3890c14ad97 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Tue, 4 Feb 2025 13:17:07 +0100 Subject: [PATCH 25/29] Rearrange the order of names. --- docs/src/about.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/about.md b/docs/src/about.md index 79e2b25339..57f8968ee0 100644 --- a/docs/src/about.md +++ b/docs/src/about.md @@ -10,10 +10,10 @@ Thanks to the following contributors to `Manopt.jl`: * [Constantin Ahlmann-Eltze](https://const-ae.name) implemented the [gradient and differential `check` functions](helpers/checks.md) * [Renée Dornig](https://github.com/r-dornig) implemented the [particle swarm](solvers/particle_swarm.md), the [Riemannian Augmented Lagrangian Method](solvers/augmented_Lagrangian_method.md), the [Exact Penalty Method](solvers/exact_penalty_method.md), as well as the [`NonmonotoneLinesearch`](@ref). These solvers are also the first one with modular/exchangable sub solvers. * [Willem Diepeveen](https://www.maths.cam.ac.uk/person/wd292) implemented the [primal-dual Riemannian semismooth Newton](solvers/primal_dual_semismooth_Newton.md) solver. -* [Sander Engen Oddsen](https://github.com/oddsen) contributed to the implementation of the [LTMADS](solvers/mesh_adaptive_direct_search.md) solver. * [Hajg Jasa](https://www.ntnu.edu/employees/hajg.jasa) implemented the [convex bundle method](solvers/convex_bundle_method.md) and the [proximal bundle method](solvers/proximal_bundle_method.md) and a default subsolver each of them. * Even Stephansen Kjemsås contributed to the implementation of the [Frank Wolfe Method](solvers/FrankWolfe.md) solver. * Mathias Ravn Munkvold contributed most of the implementation of the [Adaptive Regularization with Cubics](solvers/adaptive-regularization-with-cubics.md) solver as well as its [Lanczos](@ref arc-Lanczos) subsolver +* [Sander Engen Oddsen](https://github.com/oddsen) contributed to the implementation of the [LTMADS](solvers/mesh_adaptive_direct_search.md) solver. * [Tom-Christian Riemer](https://www.tu-chemnitz.de/mathematik/wire/mitarbeiter.php) implemented the [trust regions](solvers/trust_regions.md) and [quasi Newton](solvers/quasi_Newton.md) solvers as well as the [truncated conjugate gradient descent](solvers/truncated_conjugate_gradient_descent.md) subsolver. * [Markus A. Stokkenes](https://www.linkedin.com/in/markus-a-stokkenes-b41bba17b/) contributed most of the implementation of the [Interior Point Newton Method](solvers/interior_point_Newton.md) as well as its default [Conjugate Residual](solvers/conjugate_residual.md) subsolver * [Manuel Weiss](https://scoop.iwr.uni-heidelberg.de/author/manuel-weiß/) implemented most of the [conjugate gradient update rules](@ref cg-coeffs) From 6ba7b9ab286c7543419d3476fc4953138f8332d2 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Tue, 4 Feb 2025 14:17:06 +0100 Subject: [PATCH 26/29] Update docs/src/references.bib Co-authored-by: Sander <148573547+oddsen@users.noreply.github.com> --- docs/src/references.bib | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/references.bib b/docs/src/references.bib index d1c9745755..0f1825ab52 100644 --- a/docs/src/references.bib +++ b/docs/src/references.bib @@ -326,7 +326,7 @@ @article{DiepeveenLellmann:2021 YEAR = {2021}, } @techreport{Dreisigmeyer:2007, - AUTHOR = {Dreisigmayer, David W.}, + AUTHOR = {Dreisigmeyer, David W.}, INSTITUTION = {Optimization Online}, TITLE = {Direct Search Alogirthms over Riemannian Manifolds}, URL = {https://optimization-online.org/?p=9134}, From 58f3b1a9aeb69146743964b9e8cb6b900bdcd220 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Tue, 4 Feb 2025 16:22:14 +0100 Subject: [PATCH 27/29] fix 2 more typos. --- .vale.ini | 8 +------- Changelog.md | 20 ++++++++++---------- 2 files changed, 11 insertions(+), 17 deletions(-) diff --git a/.vale.ini b/.vale.ini index 054a57d418..a10bd713a4 100644 --- a/.vale.ini +++ b/.vale.ini @@ -12,13 +12,7 @@ jl = md [docs/src/*.md] BasedOnStyles = Vale, Google -[docs/src/contributing.md] -BasedOnStyles = Vale, Google -Google.Will = false ; given format and really with intend a _will_ -Google.Headings = false ; some might jeally ahabe [] in their headers -Google.FirstPerson = false ; we pose a few contribution points as first-person questions - -[Changelog.md, CONTRIBUTING.md] +[{docs/src/contributing.md, Changelog.md, CONTRIBUTING.md}] BasedOnStyles = Vale, Google Google.Will = false ; given format and really with intend a _will_ Google.Headings = false ; some might jeally ahabe [] in their headers diff --git a/Changelog.md b/Changelog.md index 7ba5eba06c..68894f87d5 100644 --- a/Changelog.md +++ b/Changelog.md @@ -29,16 +29,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * The geodesic regression example, first because it is not correct, second because it should become part of ManoptExamples.jl once it is correct. -## [0.5.4] - December 11, 2024 +## [0.5.4] December 11, 2024 ### Added * An automated detection whether the tutorials are present if not an also no quarto run is done, an automated `--exclude-tutorials` option is added. * Support for ManifoldDiff 0.4 -* icons upfront external links when they link to another package or wikipedia. +* icons upfront external links when they link to another package or Wikipedia. -## [0.5.3] – October 18, 2024 +## [0.5.3] October 18, 2024 ### Added @@ -50,7 +50,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 It however would warn new users, that activate tutorial mode. * Start a `ManoptTestSuite` subpackage to store dummy types and common test helpers in. -## [0.5.2] – October 5, 2024 +## [0.5.2] October 5, 2024 ### Added @@ -61,7 +61,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * fix a few typos in the documentation * improved the documentation for the initial guess of [`ArmijoLinesearchStepsize`](https://manoptjl.org/stable/plans/stepsize/#Manopt.ArmijoLinesearch). -## [0.5.1] – September 4, 2024 +## [0.5.1] September 4, 2024 ### Changed @@ -71,7 +71,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * the `proximal_point` method. -## [0.5.0] – August 29, 2024 +## [0.5.0] August 29, 2024 This breaking update is mainly concerned with improving a unified experience through all solvers and some usability improvements, such that for example the different gradient update rules are easier to specify. @@ -151,7 +151,7 @@ In general we introduce a few factories, that avoid having to pass the manifold * bumped `Manifolds.jl`to version 0.10; this mainly means that any algorithm working on a product manifold and requiring `ArrayPartition` now has to explicitly do `using RecursiveArrayTools`. ### Fixed -* the `AverageGradientRule` filled its internal vector of gradients wrongly – or mixed it up in parallel transport. This is now fixed. +* the `AverageGradientRule` filled its internal vector of gradients wrongly or mixed it up in parallel transport. This is now fixed. ### Removed @@ -171,13 +171,13 @@ In general we introduce a few factories, that avoid having to pass the manifold * to update a stopping criterion in a solver state, replace the old `update_stopping_criterion!(state, :Val, v)` tat passed down to the stopping criterion by the explicit pass down with `set_parameter!(state, :StoppingCriterion, :Val, v)` -## [0.4.69] – August 3, 2024 +## [0.4.69] August 3, 2024 ### Changed * Improved performance of Interior Point Newton Method. -## [0.4.68] – August 2, 2024 +## [0.4.68] August 2, 2024 ### Added @@ -195,7 +195,7 @@ In general we introduce a few factories, that avoid having to pass the manifold * A `StopWhenRelativeResidualLess` for the `conjugate_residual` * A `StopWhenKKTResidualLess` for the `interior_point_newton` -## [0.4.67] – July 25, 2024 +## [0.4.67] July 25, 2024 ### Added From 20901fcca388e434797e19c4134e04329ded4f26 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Tue, 4 Feb 2025 17:35:35 +0100 Subject: [PATCH 28/29] Bring vale to zero errors. --- .vale.ini | 5 +- Changelog.md | 27 +- docs/src/solvers/index.md | 2 +- .../config/vocabularies/Manopt/accept.txt | 2 + tutorials/AutomaticDifferentiation.qmd | 4 +- tutorials/InplaceGradient.html | 703 ++++++++++++++++++ tutorials/Optimize.qmd | 2 +- 7 files changed, 726 insertions(+), 19 deletions(-) create mode 100644 tutorials/InplaceGradient.html diff --git a/.vale.ini b/.vale.ini index a10bd713a4..0aab5a6f57 100644 --- a/.vale.ini +++ b/.vale.ini @@ -43,8 +43,9 @@ BasedOnStyles = Vale, Google TokenIgnores = (\$+[^\n$]+\$+) Google.We = false # For tutorials we want to address the user directly. -[docs/src/tutorials/*.md] ; actually .qmd for the first, second autogenerated +[docs/src/tutorials/*.md] ; Can I somehow just deactivate these? BasedOnStyles = Vale, Google ; ignore (1) math (2) ref and cite keys (3) code in docs (4) math in docs (5,6) indented blocks TokenIgnores = (\$+[^\n$]+\$+) -Google.We = false # For tutorials we want to address the user directly. \ No newline at end of file +Google.We = false # For tutorials we want to address the user directly. +Google.Spacing = false # one reference uses this diff --git a/Changelog.md b/Changelog.md index 68894f87d5..c0e472feb1 100644 --- a/Changelog.md +++ b/Changelog.md @@ -1,11 +1,12 @@ # Changelog -All notable Changes to the Julia package `Manopt.jl` will be documented in this file. The file was started with Version `0.4`. +All notable Changes to the Julia package `Manopt.jl` are documented in this file. +The file was started with Version `0.4`. 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.5.6] unreleased +## [0.5.6] February 10, 20265 ### Added @@ -48,7 +49,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * stabilize `max_stepsize` to also work when `injectivity_radius` dos not exist. It however would warn new users, that activate tutorial mode. -* Start a `ManoptTestSuite` subpackage to store dummy types and common test helpers in. +* Start a `ManoptTestSuite` sub package to store dummy types and common test helpers in. ## [0.5.2] October 5, 2024 @@ -76,12 +77,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 This breaking update is mainly concerned with improving a unified experience through all solvers and some usability improvements, such that for example the different gradient update rules are easier to specify. -In general we introduce a few factories, that avoid having to pass the manifold to keyword arguments +In general this introduces a few factories, that avoid having to pass the manifold to keyword arguments ### Added * A `ManifoldDefaultsFactory` that postpones the creation/allocation of manifold-specific fields in for example direction updates, step sizes and stopping criteria. As a rule of thumb, internal structures, like a solver state should store the final type. Any high-level interface, like the functions to start solvers, should accept such a factory in the appropriate places and call the internal `_produce_type(factory, M)`, for example before passing something to the state. -* a `documentation_glossary.jl` file containing a glossary of often used variables in fields, arguments, and keywords, to print them in a unified manner. The same for usual sections, tex, and math notation that is often used within the doc-strings. +* a `documentation_glossary.jl` file containing a glossary of often used variables in fields, arguments, and keywords, to print them in a unified manner. The same for usual sections, text, and math notation that is often used within the doc-strings. ### Changed @@ -106,12 +107,12 @@ In general we introduce a few factories, that avoid having to pass the manifold * `HestenesStiefelCoefficient` is now called `HestenesStiefelCoefficientRule`. For the `HestenesStiefelCoefficient` the manifold as its first parameter is no longer necessary and the vector transport has been unified/moved to the `vector_transport_method=` keyword. * `LiuStoreyCoefficient` is now called `LiuStoreyCoefficientRule`. For the `LiuStoreyCoefficient` the manifold as its first parameter is no longer necessary and the vector transport has been unified/moved to the `vector_transport_method=` keyword. * `PolakRibiereCoefficient` is now called `PolakRibiereCoefficientRule`. For the `PolakRibiereCoefficient` the manifold as its first parameter is no longer necessary and the vector transport has been unified/moved to the `vector_transport_method=` keyword. - * the `SteepestDirectionUpdateRule` is now called `SteepestDescentCoefficientRule`. The `SteepestDescentCoefficient` is equivalent, but creates the new factory interims wise. + * the `SteepestDirectionUpdateRule` is now called `SteepestDescentCoefficientRule`. The `SteepestDescentCoefficient` is equivalent, but creates the new factory temporarily. * `AbstractGradientGroupProcessor` is now called `AbstractGradientGroupDirectionRule` - * the `StochasticGradient` is now called `StochasticGradientRule`. The `StochasticGradient` is equivalent, but creates the new factory interims wise, so that the manifold is not longer necessary. + * the `StochasticGradient` is now called `StochasticGradientRule`. The `StochasticGradient` is equivalent, but creates the new factory temporarily, so that the manifold is not longer necessary. * the `AlternatingGradient` is now called `AlternatingGradientRule`. - The `AlternatingGradient` is equivalent, but creates the new factory interims wise, so that the manifold is not longer necessary. -* `quasi_Newton` had a keyword `scale_initial_operator=` that was inconsistently declared (sometimes bool, sometimes real) and was unused. + The `AlternatingGradient` is equivalent, but creates the new factory temporarily, so that the manifold is not longer necessary. +* `quasi_Newton` had a keyword `scale_initial_operator=` that was inconsistently declared (sometimes boolean, sometimes real) and was unused. It is now called `initial_scale=1.0` and scales the initial (diagonal, unit) matrix within the approximation of the Hessian additionally to the $\frac{1}{\lVert g_k\rVert}$ scaling with the norm of the oldest gradient for the limited memory variant. For the full matrix variant the initial identity matrix is now scaled with this parameter. * Unify doc strings and presentation of keyword arguments * general indexing, for example in a vector, uses `i` @@ -144,7 +145,7 @@ In general we introduce a few factories, that avoid having to pass the manifold * `AdaptiveRegularizationState(M, sub_problem [, sub_state]; kwargs...)` replaces the (anyways unused) variant to only provide the objective; both `X` and `p` moved to keyword arguments. * `AugmentedLagrangianMethodState(M, objective, sub_problem; evaluation=...)` was added - * ``AugmentedLagrangianMethodState(M, objective, sub_problem, sub_state; evaluation=...)` now has `p=rand(M)` as keyword argument instead of being the second positional one + * `AugmentedLagrangianMethodState(M, objective, sub_problem, sub_state; evaluation=...)` now has `p=rand(M)` as keyword argument instead of being the second positional one * `ExactPenaltyMethodState(M, sub_problem; evaluation=...)` was added and `ExactPenaltyMethodState(M, sub_problem, sub_state; evaluation=...)` now has `p=rand(M)` as keyword argument instead of being the second positional one * `DifferenceOfConvexState(M, sub_problem; evaluation=...)` was added and `DifferenceOfConvexState(M, sub_problem, sub_state; evaluation=...)` now has `p=rand(M)` as keyword argument instead of being the second positional one * `DifferenceOfConvexProximalState(M, sub_problem; evaluation=...)` was added and `DifferenceOfConvexProximalState(M, sub_problem, sub_state; evaluation=...)` now has `p=rand(M)` as keyword argument instead of being the second positional one @@ -185,9 +186,9 @@ In general we introduce a few factories, that avoid having to pass the manifold * a `conjugate_residual` Algorithm to solve a linear system on a tangent space. * `ArmijoLinesearch` now allows for additional `additional_decrease_condition` and `additional_increase_condition` keywords to add further conditions to accept additional conditions when to accept an decreasing or increase of the stepsize. * add a `DebugFeasibility` to have a debug print about feasibility of points in constrained optimisation employing the new `is_feasible` function -* add a `InteriorPointCentralityCondition` check that can be added for step candidates within the line search of `interior_point_newton` +* add a `InteriorPointCentralityCondition` that can be added for step candidates within the line search of `interior_point_newton` * Add Several new functors - * the `LagrangianCost`, `LagrangianGradient`, `LagrangianHessian`, that based on a constrained objective allow to construct the hessian objective of its Lagrangian + * the `LagrangianCost`, `LagrangianGradient`, `LagrangianHessian`, that based on a constrained objective allow to construct the Hessian objective of its Lagrangian * the `CondensedKKTVectorField` and its `CondensedKKTVectorFieldJacobian`, that are being used to solve a linear system within `interior_point_newton` * the `KKTVectorField` as well as its `KKTVectorFieldJacobian` and ``KKTVectorFieldAdjointJacobian` * the `KKTVectorFieldNormSq` and its `KKTVectorFieldNormSqGradient` used within the Armijo line search of `interior_point_newton` @@ -241,7 +242,7 @@ In general we introduce a few factories, that avoid having to pass the manifold * Remodel `ConstrainedManifoldObjective` to store an `AbstractManifoldObjective` internally instead of directly `f` and `grad_f`, allowing also Hessian objectives therein and implementing access to this Hessian -* Fixed a bug that Lanczos produced NaNs when started exactly in a minimizer, since we divide by the gradient norm. +* Fixed a bug that Lanczos produced NaNs when started exactly in a minimizer, since the algorithm initially divides by the gradient norm. ### Deprecated diff --git a/docs/src/solvers/index.md b/docs/src/solvers/index.md index ba5746cbd6..6509e739e9 100644 --- a/docs/src/solvers/index.md +++ b/docs/src/solvers/index.md @@ -22,7 +22,7 @@ For derivative free only function evaluations of ``f`` are used. * [Nelder-Mead](NelderMead.md) a simplex based variant, that is using ``d+1`` points, where ``d`` is the dimension of the manifold. * [Particle Swarm](particle_swarm.md) 🫏 use the evolution of a set of points, called swarm, to explore the domain of the cost and find a minimizer. -* [Mesh adapive direct search](mesh_adaptive_direct_search.md) performs a mesh based exploration (poll) and search. +* [Mesh adaptive direct search](mesh_adaptive_direct_search.md) performs a mesh based exploration (poll) and search. * [CMA-ES](cma_es.md) uses a stochastic evolutionary strategy to perform minimization robust to local minima of the objective. ## First order diff --git a/docs/styles/config/vocabularies/Manopt/accept.txt b/docs/styles/config/vocabularies/Manopt/accept.txt index 73e5f1fef9..2ebb074b6b 100644 --- a/docs/styles/config/vocabularies/Manopt/accept.txt +++ b/docs/styles/config/vocabularies/Manopt/accept.txt @@ -34,6 +34,7 @@ eigen eigendecomposition elementwise Ehresmann +Engen Fenchel Ferreira Frank @@ -82,6 +83,7 @@ Newton nonmonotone nonpositive [Nn]onsmooth +Oddsen [Pp]arametrising Parametrising [Pp]ock diff --git a/tutorials/AutomaticDifferentiation.qmd b/tutorials/AutomaticDifferentiation.qmd index fe5da3abc9..3684eee52a 100644 --- a/tutorials/AutomaticDifferentiation.qmd +++ b/tutorials/AutomaticDifferentiation.qmd @@ -1,5 +1,5 @@ --- -title: "Using Automatic Differentiation in Manopt.jl" +title: "Using automatic differentiation in Manopt.jl" --- Since [Manifolds.jl](https://juliamanifolds.github.io/Manifolds.jl/latest/) 0.7, the support of automatic differentiation support has been extended. @@ -122,7 +122,7 @@ norm(M, p, X1 - X2) We obtain quite a good approximation of the gradient. -``## [2. Conversion of a Euclidean Gradient in the Embedding to a Riemannian Gradient of a (not Necessarily Isometrically) Embedded Manifold](@id EmbeddedGradient)``{=commonmark} +``## [2. Conversion of a Euclidean gradient in the embedding to a Riemannian Gradient of a (not Necessarily Isometrically) embedded manifold](@id EmbeddedGradient)``{=commonmark} Let $\tilde f: ℝ^m → ℝ$ be a function on the embedding of an $n$-dimensional manifold $\mathcal M \subset ℝ^m$and let $f: \mathcal M → ℝ$ denote the restriction of $\tilde f$ to the manifold $\mathcal M$. diff --git a/tutorials/InplaceGradient.html b/tutorials/InplaceGradient.html new file mode 100644 index 0000000000..f668af1857 --- /dev/null +++ b/tutorials/InplaceGradient.html @@ -0,0 +1,703 @@ + + + + + + + + + + +Speedup using in-place evaluation + + + + + + + + + + + + + + + + + + + + + + + +
+ +
+ +
+
+

Speedup using in-place evaluation

+
+ + + +
+ +
+
Author
+
+

Ronny Bergmann

+
+
+ + + +
+ + + +
+ + +

When it comes to time critical operations, a main ingredient in Julia is given by mutating functions, that is those that compute in place without additional memory allocations. In the following, we illustrate how to do this with Manopt.jl.

+

Let’s start with the same function as in 🏔️ Get started with Manopt.jl and compute the mean of some points, only that here we use the sphere \(\mathbb S^{30}\) and \(n=800\) points.

+

From the aforementioned example.

+

We first load all necessary packages.

+
+
using Manopt, Manifolds, Random, BenchmarkTools
+using ManifoldDiff: grad_distance, grad_distance!
+Random.seed!(42);
+
+

And setup our data

+
+
Random.seed!(42)
+m = 30
+M = Sphere(m)
+n = 800
+σ = π / 8
+p = zeros(Float64, m + 1)
+p[2] = 1.0
+data = [exp(M, p, σ * rand(M; vector_at=p)) for i in 1:n];
+
+
+

Classical definition

+

The variant from the previous tutorial defines a cost \(f(x)\) and its gradient \(\operatorname{grad}f(p)\) ““”

+
+
f(M, p) = sum(1 / (2 * n) * distance.(Ref(M), Ref(p), data) .^ 2)
+grad_f(M, p) = sum(1 / n * grad_distance.(Ref(M), data, Ref(p)))
+
+
grad_f (generic function with 1 method)
+
+
+

We further set the stopping criterion to be a little more strict. Then we obtain

+
+
sc = StopWhenGradientNormLess(3e-10)
+p0 = zeros(Float64, m + 1); p0[1] = 1/sqrt(2); p0[2] = 1/sqrt(2)
+m1 = gradient_descent(M, f, grad_f, p0; stopping_criterion=sc);
+
+

We can also benchmark this as

+
+
@benchmark gradient_descent($M, $f, $grad_f, $p0; stopping_criterion=$sc)
+
+
+
BenchmarkTools.Trial: 86 samples with 1 evaluation.
+ Range (min … max):  52.720 ms … 93.308 ms  ┊ GC (min … max):  8.27% … 11.47%
+ Time  (median):     55.064 ms              ┊ GC (median):    10.10%
+ Time  (mean ± σ):   58.153 ms ±  7.376 ms  ┊ GC (mean ± σ):  10.36% ±  1.44%
+   ▃█▅▅▄    ▁                                                  
+  ▅█████▅█▁▅█▇▁▁▅▁▁▁▁▁▅▅▇▁▁▅▇▁▁▇▅▁▅▁▁▁▁▅▅▁▁▁▁▁▁▁▁▁▁▁▅▁▁▁▁▁▁▁▅ ▁
+  52.7 ms      Histogram: log(frequency) by time      84.9 ms <
+ Memory estimate: 173.54 MiB, allocs estimate: 1167345.
+
+
+
+
+
+

In-place computation of the gradient

+

We can reduce the memory allocations by implementing the gradient to be evaluated in-place. We do this by using a functor. The motivation is twofold: on one hand, we want to avoid variables from the global scope, for example the manifold M or the data, being used within the function. Considering to do the same for more complicated cost functions might also be worth pursuing.

+

Here, we store the data (as reference) and one introduce temporary memory to avoid reallocation of memory per grad_distance computation. We get

+
+
struct GradF!{TD,TTMP}
+    data::TD
+    tmp::TTMP
+end
+function (grad_f!::GradF!)(M, X, p)
+    fill!(X, 0)
+    for di in grad_f!.data
+        grad_distance!(M, grad_f!.tmp, di, p)
+        X .+= grad_f!.tmp
+    end
+    X ./= length(grad_f!.data)
+    return X
+end
+
+

For the actual call to the solver, we first have to generate an instance of GradF! and tell the solver, that the gradient is provided in an InplaceEvaluation. We can further also use gradient_descent! to even work in-place of the initial point we pass.

+
+
grad_f2! = GradF!(data, similar(data[1]))
+m2 = deepcopy(p0)
+gradient_descent!(
+    M, f, grad_f2!, m2; evaluation=InplaceEvaluation(), stopping_criterion=sc
+);
+
+

We can again benchmark this

+
+
@benchmark gradient_descent!(
+    $M, $f, $grad_f2!, m2; evaluation=$(InplaceEvaluation()), stopping_criterion=$sc
+) setup = (m2 = deepcopy($p0))
+
+
+
BenchmarkTools.Trial: 135 samples with 1 evaluation.
+ Range (min … max):  35.592 ms … 59.467 ms  ┊ GC (min … max): 0.00% … 0.00%
+ Time  (median):     36.393 ms              ┊ GC (median):    0.00%
+ Time  (mean ± σ):   37.177 ms ±  3.086 ms  ┊ GC (mean ± σ):  0.64% ± 2.40%
+  ▄█▇▃▄                                                        
+  █████▇▆▇▅▅▆▅▁▅▁▁▁▁▁▁▁▁▁▁▆▁▁▁▁▁▅▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅ ▅
+  35.6 ms      Histogram: log(frequency) by time      57.8 ms <
+ Memory estimate: 3.59 MiB, allocs estimate: 6860.
+
+
+
+

which is faster by about a factor of 2 compared to the first solver-call. Note that the results m1 and m2 are of course the same.

+
+
distance(M, m1, m2)
+
+
4.8317610992693745e-11
+
+
+
+
+

Technical details

+

This tutorial is cached. It was last run on the following package versions.

+
+
+Code +
using Pkg
+Pkg.status()
+
+
+
Status `~/Repositories/Julia/Manopt.jl/tutorials/Project.toml`
+⌃ [6e4b80f9] BenchmarkTools v1.5.0
+⌃ [5ae59095] Colors v0.12.11
+⌃ [31c24e10] Distributions v0.25.113
+  [26cc04aa] FiniteDifferences v0.12.32
+  [7073ff75] IJulia v1.26.0
+  [8ac3fa9e] LRUCache v1.6.1
+⌅ [af67fdf4] ManifoldDiff v0.3.13
+⌃ [1cead3c2] Manifolds v0.10.7
+⌃ [3362f125] ManifoldsBase v0.15.22
+  [0fc0a36d] Manopt v0.5.5 `~/Repositories/Julia/Manopt.jl`
+  [91a5bcdd] Plots v1.40.9
+⌃ [731186ca] RecursiveArrayTools v3.27.4
+Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`
+
+
+
+
+Code +
using Dates
+now()
+
+
+
2025-02-04T17:30:59.608
+
+
+
+ +
+ + +
+ + + + + \ No newline at end of file diff --git a/tutorials/Optimize.qmd b/tutorials/Optimize.qmd index 1e8b3bda74..aac1741841 100644 --- a/tutorials/Optimize.qmd +++ b/tutorials/Optimize.qmd @@ -1,5 +1,5 @@ --- -title: "🏔️ Get started with Manopt.jl." +title: "🏔️ Get started with Manopt.jl" author: Ronny Bergmann --- From 8cd4880413c89d71de3c47bfb9e42ac94da15996 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Wed, 5 Feb 2025 10:03:40 +0100 Subject: [PATCH 29/29] Fix a few more typos. --- docs/src/tutorials/InplaceGradient.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/src/tutorials/InplaceGradient.md b/docs/src/tutorials/InplaceGradient.md index 47f64f447f..a6a503caab 100644 --- a/docs/src/tutorials/InplaceGradient.md +++ b/docs/src/tutorials/InplaceGradient.md @@ -5,7 +5,7 @@ When it comes to time critical operations, a main ingredient in Julia is given b mutating functions, that is those that compute in place without additional memory allocations. In the following, we illustrate how to do this with `Manopt.jl`. -Let’s start with the same function as in [Get started: optimize!](https://manoptjl.org/stable/tutorials/Optimize!.html) +Let’s start with the same function as in [🏔️ Get started with Manopt.jl](https://manoptjl.org/stable/tutorials/Optimize!.html) and compute the mean of some points, only that here we use the sphere $\mathbb S^{30}$ and $n=800$ points. @@ -32,7 +32,7 @@ p[2] = 1.0 data = [exp(M, p, σ * rand(M; vector_at=p)) for i in 1:n]; ``` -## Classical Definition +## Classical definition The variant from the previous tutorial defines a cost $f(x)$ and its gradient $\operatorname{grad}f(p)$ ““” @@ -63,13 +63,13 @@ We can also benchmark this as Time (median): 47.207 ms ┊ GC (median): 2.45% Time (mean ± σ): 47.364 ms ± 608.514 μs ┊ GC (mean ± σ): 2.53% ± 0.25% - ▄▇▅▇█▄▇ + ▄▇▅▇█▄▇ ▅▇▆████████▇▇▅▅▃▁▆▁▁▁▅▁▁▅▁▃▃▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅ ▃ 46.8 ms Histogram: frequency by time 50.2 ms < Memory estimate: 182.50 MiB, allocs estimate: 615822. -## In-place Computation of the Gradient +## In-place computation of the gradient We can reduce the memory allocations by implementing the gradient to be evaluated in-place. We do this by using a [functor](https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects). @@ -77,7 +77,7 @@ The motivation is twofold: on one hand, we want to avoid variables from the glob for example the manifold `M` or the `data`, being used within the function. Considering to do the same for more complicated cost functions might also be worth pursuing. -Here, we store the data (as reference) and one introduce temporary memory in order to avoid +Here, we store the data (as reference) and one introduce temporary memory to avoid reallocation of memory per `grad_distance` computation. We get ``` julia @@ -121,7 +121,7 @@ We can again benchmark this Time (median): 27.768 ms ┊ GC (median): 0.00% Time (mean ± σ): 28.504 ms ± 4.338 ms ┊ GC (mean ± σ): 0.60% ± 1.96% - ▂█▇▂ ▂ + ▂█▇▂ ▂ ▆▇████▆█▆▆▄▄▃▄▄▃▃▃▁▃▃▃▃▃▃▃▃▃▄▃▃▃▃▃▃▁▃▁▁▃▁▁▁▁▁▁▃▃▁▁▃▃▁▁▁▁▃▃▃ ▃ 27.4 ms Histogram: frequency by time 31.4 ms <