Rework high level interfaces

Changelog.md

@@ -5,12 +5,9 @@ All notable Changes to the Julia package `Manopt.jl` will be documented in this
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.4.64] unreleased
 
-## [0.4.63] June 4, 2024
-
-### Changed
-
-* Fixed a bug that Lanczos produced NaNs when started exactly in a minimizer, since we divide by the gradient norm.
+* Fix ambiguities that occurred due to point being nonmutating.
 
 ## [0.4.63] May 11, 2024
 

src/plans/gradient_plan.jl

@@ -221,13 +221,6 @@ function get_subgradient!(
     return get_gradient!(M, X, agmo, p)
 end
 
-function _to_mutating_gradient(grad_f, evaluation::AllocatingEvaluation)
-    return grad_f_(M, p) = [grad_f(M, p[])]
-end
-function _to_mutating_gradient(grad_f, evaluation::InplaceEvaluation)
-    return grad_f_(M, X, p) = (X .= [grad_f(M, p[])])
-end
-
 @doc raw"""
     get_gradient(agst::AbstractGradientSolverState)
 
@@ -458,7 +451,6 @@ mutable struct Nesterov{P,R<:Real} <: DirectionUpdateRule
     shrinkage::Function
     inverse_retraction_method::AbstractInverseRetractionMethod
 end
-Nesterov(M::AbstractManifold, p::Number; kwargs...) = Nesterov(M, [p]; kwargs...)
 function Nesterov(
     M::AbstractManifold,
     p::P;
@@ -469,7 +461,8 @@ function Nesterov(
         M, P
     ),
 ) where {P,T}
-    return Nesterov{P,T}(γ, μ, copy(M, p), shrinkage, inverse_retraction_method)
+    p_ = _ensure_mutating_variable(p)
+    return Nesterov{typeof(p_),T}(γ, μ, copy(M, p_), shrinkage, inverse_retraction_method)
 end
 function (n::Nesterov)(mp::AbstractManoptProblem, s::AbstractGradientSolverState, i)
     M = get_manifold(mp)

src/plans/objective.jl

@@ -75,6 +75,54 @@ function ReturnManifoldObjective(
     return ReturnManifoldObjective{E,O2,O1}(o)
 end
 
+function _ensure_mutating_cost(cost, q::Number)
+    return isnothing(cost) ? cost : (M, p) -> cost(M, p[])
+end
+function _ensure_mutating_cost(cost, p)
+    return cost
+end
+
+function _ensure_mutating_gradient(grad_f, p, evaluation::AbstractEvaluationType)
+    return grad_f
+end
+function _ensure_mutating_gradient(grad_f, q::Number, evaluation::AllocatingEvaluation)
+    return isnothing(grad_f) ? grad_f : (M, p) -> [grad_f(M, p[])]
+end
+function _ensure_mutating_gradient(grad_f, q::Number, evaluation::InplaceEvaluation)
+    return isnothing(grad_f) ? grad_f : (M, X, p) -> (X .= [grad_f(M, p[])])
+end
+
+function _ensure_mutating_hessian(hess_f, p, evaluation::AbstractEvaluationType)
+    return hess_f
+end
+function _ensure_mutating_hessian(hess_f, q::Number, evaluation::AllocatingEvaluation)
+    return isnothing(hess_f) ? hess_f : (M, p, X) -> [hess_f(M, p[], X[])]
+end
+function _ensure_mutating_hessian(hess_f, q::Number, evaluation::InplaceEvaluation)
+    return isnothing(hess_f) ? hess_f : (M, Y, p, X) -> (Y .= [hess_f(M, p[], X[])])
+end
+
+function _ensure_mutating_prox(prox_f, p, evaluation::AbstractEvaluationType)
+    return prox_f
+end
+function _ensure_mutating_prox(prox_f, q::Number, evaluation::AllocatingEvaluation)
+    return isnothing(prox_f) ? prox_f : (M, λ, p) -> [prox_f(M, λ, p[])]
+end
+function _ensure_mutating_prox(prox_f, q::Number, evaluation::InplaceEvaluation)
+    return isnothing(prox_f) ? prox_f : (M, q, λ, p) -> (q .= [prox_f(M, λ, p[])])
+end
+
+_ensure_mutating_variable(p) = p
+_ensure_mutating_variable(q::Number) = [q]
+"""
+    _ensure_matching_output(p, q)
+    _ensure_matching_output(e, q, s)
+
+If p is a number and q is a vector, return q[] (a number) otherwise q
+"""
+_ensure_matching_output(::T, q::Vector{T}) where {T} = length(q) == 1 ? q[] : q
+_ensure_matching_output(p, q) = q
+
 """
     dispatch_objective_decorator(o::AbstractManoptSolverState)
 

src/solvers/DouglasRachford.jl

@@ -189,30 +189,15 @@ function DouglasRachford(
     parallel=0,
     kwargs...,
 ) where {TF}
-    N, f_, (prox1, prox2), parallel_, p0 = parallel_to_alternating_DR(
-        M, f, proxes_f, p, parallel, evaluation
+    p_ = _ensure_mutating_variable(p)
+    f_ = _ensure_mutating_cost(f, p)
+    proxes_f_ = [_ensure_mutating_prox(prox_f, p, evaluation) for prox_f in proxes_f]
+    N, f__, (prox1, prox2), parallel_, q = parallel_to_alternating_DR(
+        M, f_, proxes_f_, p_, parallel, evaluation
     )
     mpo = ManifoldProximalMapObjective(f_, (prox1, prox2); evaluation=evaluation)
-    return DouglasRachford(N, mpo, p0; evaluation=evaluation, parallel=parallel_, kwargs...)
-end
-function DouglasRachford(
-    M::AbstractManifold,
-    f::TF,
-    proxes_f::Vector{<:Any},
-    p::Number;
-    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
-    kwargs...,
-) where {TF}
-    q = [p]
-    f_(M, p) = f(M, p[])
-    if evaluation isa AllocatingEvaluation
-        proxes_f_ = [(M, λ, p) -> [pf(M, λ, p[])] for pf in proxes_f]
-    else
-        proxes_f_ = [(M, q, λ, p) -> (q .= [pf(M, λ, p[])]) for pf in proxes_f]
-    end
-    rs = DouglasRachford(M, f_, proxes_f_, q; evaluation=evaluation, kwargs...)
-    #return just a number if  the return type is the same as the type of q
-    return (typeof(q) == typeof(rs)) ? rs[] : rs
+    rs = DouglasRachford(N, mpo, q; evaluation=evaluation, parallel=parallel_, kwargs...)
+    return _ensure_matching_output(p, rs)
 end
 function DouglasRachford(
     M::AbstractManifold, mpo::O, p; kwargs...

src/solvers/FrankWolfe.jl

@@ -198,24 +198,12 @@ function Frank_Wolfe_method(
     evaluation::AbstractEvaluationType=AllocatingEvaluation(),
     kwargs...,
 )
-    mgo = ManifoldGradientObjective(f, grad_f; evaluation=evaluation)
-    return Frank_Wolfe_method(M, mgo, p; evaluation=evaluation, kwargs...)
-end
-function Frank_Wolfe_method(
-    M::AbstractManifold,
-    f,
-    grad_f,
-    p::Number;
-    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
-    kwargs...,
-)
-    # redefine initial point
-    q = [p]
-    f_(M, p) = f(M, p[])
-    grad_f_ = _to_mutating_gradient(grad_f, evaluation)
-    rs = Frank_Wolfe_method(M, f_, grad_f_, q; evaluation=evaluation, kwargs...)
-    #return just a number if  the return type is the same as the type of q
-    return (typeof(q) == typeof(rs)) ? rs[] : rs
+    p_ = _ensure_mutating_variable(p)
+    f_ = _ensure_mutating_cost(f, p)
+    grad_f_ = _ensure_mutating_gradient(grad_f, p, evaluation)
+    mgo = ManifoldGradientObjective(f_, grad_f_; evaluation=evaluation)
+    rs = Frank_Wolfe_method(M, mgo, p_; evaluation=evaluation, kwargs...)
+    return _ensure_matching_output(p, rs)
 end
 function Frank_Wolfe_method(
     M::AbstractManifold, mgo::O, p; kwargs...

src/solvers/NelderMead.jl

@@ -31,21 +31,19 @@ end
 function NelderMeadSimplex(M::AbstractManifold)
     return NelderMeadSimplex([rand(M) for i in 1:(manifold_dimension(M) + 1)])
 end
-function NelderMeadSimplex(M::AbstractManifold, p::Number, B::AbstractBasis; kwargs...)
-    return NelderMeadSimplex(M, [p], B; kwargs...)
-end
 function NelderMeadSimplex(
     M::AbstractManifold,
     p,
     B::AbstractBasis=DefaultOrthonormalBasis();
     a::Real=0.025,
     retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)),
 )
+    p_ = _ensure_mutating_variable(p)
     M_dim = manifold_dimension(M)
     vecs = [
-        get_vector(M, p, [ifelse(i == j, a, zero(a)) for i in 1:M_dim], B) for j in 0:M_dim
+        get_vector(M, p_, [ifelse(i == j, a, zero(a)) for i in 1:M_dim], B) for j in 0:M_dim
     ]
-    pts = map(X -> retract(M, p, X, retraction_method), vecs)
+    pts = map(X -> retract(M, p_, X, retraction_method), vecs)
     return NelderMeadSimplex(pts)
 end
 

src/solvers/adaptive_regularization_with_cubics.jl

@@ -257,32 +257,13 @@ function adaptive_regularization_with_cubics(
     evaluation::AbstractEvaluationType=AllocatingEvaluation(),
     kwargs...,
 ) where {TF,TDF,THF}
-    mho = ManifoldHessianObjective(f, grad_f, Hess_f; evaluation=evaluation)
-    return adaptive_regularization_with_cubics(M, mho, p; evaluation=evaluation, kwargs...)
-end
-function adaptive_regularization_with_cubics(
-    M::AbstractManifold,
-    f::TF,
-    grad_f::TDF,
-    Hess_f::THF,
-    p::Number;
-    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
-    kwargs...,
-) where {TF,TDF,THF}
-    q = [p]
-    f_(M, p) = f(M, p[])
-    Hess_f_ = Hess_f
-    if evaluation isa AllocatingEvaluation
-        grad_f_ = (M, p) -> [grad_f(M, p[])]
-        Hess_f_ = (M, p, X) -> [Hess_f(M, p[], X[])]
-    else
-        grad_f_ = (M, X, p) -> (X .= [grad_f(M, p[])])
-        Hess_f_ = (M, Y, p, X) -> (Y .= [Hess_f(M, p[], X[])])
-    end
-    rs = adaptive_regularization_with_cubics(
-        M, f_, grad_f_, Hess_f_, q; evaluation=evaluation, kwargs...
-    )
-    return (typeof(q) == typeof(rs)) ? rs[] : rs
+    p_ = _ensure_mutating_variable(p)
+    f_ = _ensure_mutating_cost(f, p)
+    grad_f_ = _ensure_mutating_gradient(grad_f, p, evaluation)
+    Hess_f_ = _ensure_mutating_hessian(Hess_f, p, evaluation)
+    mho = ManifoldHessianObjective(f_, grad_f_, Hess_f_; evaluation=evaluation)
+    rs = adaptive_regularization_with_cubics(M, mho, p_; evaluation=evaluation, kwargs...)
+    return _ensure_matching_output(p, rs)
 end
 function adaptive_regularization_with_cubics(M::AbstractManifold, f, grad_f; kwargs...)
     return adaptive_regularization_with_cubics(M, f, grad_f, rand(M); kwargs...)

src/solvers/augmented_Lagrangian_method.jl

@@ -267,43 +267,27 @@ function augmented_Lagrangian_method(
     grad_h=nothing,
     kwargs...,
 ) where {TF,TGF}
-    q = copy(M, p)
+    p_ = _ensure_mutating_variable(p)
+    f_ = _ensure_mutating_cost(f, p)
+    grad_f_ = _ensure_mutating_gradient(grad_f, p, evaluation)
+    g_ = _ensure_mutating_cost(g, p)
+    grad_g_ = _ensure_mutating_gradient(grad_g, p, evaluation)
+    h_ = _ensure_mutating_cost(h, p)
+    grad_h_ = _ensure_mutating_gradient(grad_h, p, evaluation)
+
     cmo = ConstrainedManifoldObjective(
-        f, grad_f, g, grad_g, h, grad_h; evaluation=evaluation
+        f_, grad_f_, g_, grad_g_, h_, grad_h_; evaluation=evaluation
     )
-    return augmented_Lagrangian_method!(M, cmo, q; evaluation=evaluation, kwargs...)
+    rs = augmented_Lagrangian_method(M, cmo, p_; evaluation=evaluation, kwargs...)
+    return _ensure_matching_output(p, rs)
 end
 function augmented_Lagrangian_method(
-    M::AbstractManifold, cmo::O, p=rand(M); kwargs...
+    M::AbstractManifold, cmo::O, p=rand(M);
+    kwargs...
 ) where {O<:Union{ConstrainedManifoldObjective,AbstractDecoratedManifoldObjective}}
     q = copy(M, p)
     return augmented_Lagrangian_method!(M, cmo, q; kwargs...)
 end
-function augmented_Lagrangian_method(
-    M::AbstractManifold,
-    f::TF,
-    grad_f::TGF,
-    p::Number;
-    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
-    g=nothing,
-    grad_g=nothing,
-    grad_h=nothing,
-    h=nothing,
-    kwargs...,
-) where {TF,TGF}
-    q = [p]
-    f_(M, p) = f(M, p[])
-    grad_f_ = _to_mutating_gradient(grad_f, evaluation)
-    g_ = isnothing(g) ? nothing : (M, p) -> g(M, p[])
-    grad_g_ = isnothing(grad_g) ? nothing : _to_mutating_gradient(grad_g, evaluation)
-    h_ = isnothing(h) ? nothing : (M, p) -> h(M, p[])
-    grad_h_ = isnothing(grad_h) ? nothing : _to_mutating_gradient(grad_h, evaluation)
-    cmo = ConstrainedManifoldObjective(
-        f_, grad_f_, g_, grad_g_, h_, grad_h_; evaluation=evaluation
-    )
-    rs = augmented_Lagrangian_method(M, cmo, q; evaluation=evaluation, kwargs...)
-    return (typeof(q) == typeof(rs)) ? rs[] : rs
-end
 
 @doc raw"""
     augmented_Lagrangian_method!(M, f, grad_f, p=rand(M); kwargs...)

src/solvers/conjugate_gradient_descent.jl

@@ -94,24 +94,12 @@ end
 function conjugate_gradient_descent(
     M::AbstractManifold, f::TF, grad_f::TDF, p; evaluation=AllocatingEvaluation(), kwargs...
 ) where {TF,TDF}
-    mgo = ManifoldGradientObjective(f, grad_f; evaluation=evaluation)
-    return conjugate_gradient_descent(M, mgo, p; evaluation=evaluation, kwargs...)
-end
-function conjugate_gradient_descent(
-    M::AbstractManifold,
-    f::TF,
-    grad_f::TDF,
-    p::Number;
-    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
-    kwargs...,
-) where {TF,TDF}
-    # redefine initial point
-    q = [p]
-    f_(M, p) = f(M, p[])
-    grad_f_ = _to_mutating_gradient(grad_f, evaluation)
-    rs = conjugate_gradient_descent(M, f_, grad_f_, q; evaluation=evaluation, kwargs...)
-    #return just a number if  the return type is the same as the type of q
-    return (typeof(q) == typeof(rs)) ? rs[] : rs
+    p_ = _ensure_mutating_variable(p)
+    f_ = _ensure_mutating_cost(f, p)
+    grad_f_ = _ensure_mutating_gradient(grad_f, p, evaluation)
+    mgo = ManifoldGradientObjective(f_, grad_f_; evaluation=evaluation)
+    rs = conjugate_gradient_descent(M, mgo, p_; evaluation=evaluation, kwargs...)
+    return _ensure_matching_output(p, rs)
 end
 function conjugate_gradient_descent(
     M::AbstractManifold, mgo::O, p=rand(M); kwargs...

src/solvers/cyclic_proximal_point.jl

@@ -57,27 +57,12 @@ function cyclic_proximal_point(
     evaluation::AbstractEvaluationType=AllocatingEvaluation(),
     kwargs...,
 )
-    mpo = ManifoldProximalMapObjective(f, proxes_f; evaluation=evaluation)
-    return cyclic_proximal_point(M, mpo, p; evaluation=evaluation, kwargs...)
-end
-function cyclic_proximal_point(
-    M::AbstractManifold,
-    f,
-    proxes_f::Union{Tuple,AbstractVector},
-    p::Number;
-    evaluation::AbstractEvaluationType=AllocatingEvaluation(),
-    kwargs...,
-)
-    q = [p]
-    f_(M, p) = f(M, p[])
-    if evaluation isa AllocatingEvaluation
-        proxes_f_ = [(M, λ, p) -> [pf(M, λ, p[])] for pf in proxes_f]
-    else
-        proxes_f_ = [(M, q, λ, p) -> (q .= [pf(M, λ, p[])]) for pf in proxes_f]
-    end
-    rs = cyclic_proximal_point(M, f_, proxes_f_, q; evaluation=evaluation, kwargs...)
-    #return just a number if  the return type is the same as the type of q
-    return (typeof(q) == typeof(rs)) ? rs[] : rs
+    p_ = _ensure_mutating_variable(p)
+    f_ = _ensure_mutating_cost(f, p)
+    proxes_f_ = [_ensure_mutating_prox(prox_f, p, evaluation) for prox_f in proxes_f]
+    mpo = ManifoldProximalMapObjective(f_, proxes_f_; evaluation=evaluation)
+    rs = cyclic_proximal_point(M, mpo, p_; evaluation=evaluation, kwargs...)
+    return _ensure_matching_output(p, rs)
 end
 function cyclic_proximal_point(
     M::AbstractManifold, mpo::O, p; kwargs...