From b1c53a117dbb50659a6932d49fdeaa69e2c0bd7f Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Tue, 31 Jan 2023 22:01:36 +0100 Subject: [PATCH 01/37] Even faster storage handling --- src/plans/conjugate_gradient_plan.jl | 298 ++++++++++++---------- src/plans/solver_state.jl | 146 +++++++++-- src/plans/stepsize.jl | 24 +- src/solvers/conjugate_gradient_descent.jl | 7 +- test/solvers/test_conjugate_gradient.jl | 17 +- 5 files changed, 309 insertions(+), 183 deletions(-) diff --git a/src/plans/conjugate_gradient_plan.jl b/src/plans/conjugate_gradient_plan.jl index cc13f33428..9cdf017a5a 100644 --- a/src/plans/conjugate_gradient_plan.jl +++ b/src/plans/conjugate_gradient_plan.jl @@ -1,3 +1,22 @@ + +struct DirectionUpdateRuleStorage{TC<:DirectionUpdateRule,TStorage<:StoreStateAction} + coefficient::TC + storage::TStorage +end + +function DirectionUpdateRuleStorage( + M::AbstractManifold, dur::DirectionUpdateRule, p0=rand(M), X0=zero_vector(M, p0) +) + ursp = update_rule_storage_points(dur) + ursv = update_rule_storage_vectors(dur) + sa = StoreStateAction( + Symbol[], + NamedTuple{ursp}(map(x -> copy(M, p0), ursp)), + NamedTuple{ursv}(map(x -> copy(M, p0, X0), ursv)), + ) + return DirectionUpdateRuleStorage{typeof(dur),typeof(sa)}(dur, sa) +end + @doc raw""" ConjugateGradientState <: AbstractGradientSolverState @@ -23,7 +42,7 @@ mutable struct ConjugateGradientDescentState{ P, T, F, - TCoeff<:DirectionUpdateRule, + TCoeff<:DirectionUpdateRuleStorage, TStepsize<:Stepsize, TStop<:StoppingCriterion, TRetr<:AbstractRetractionMethod, @@ -43,7 +62,7 @@ mutable struct ConjugateGradientDescentState{ p::P, sC::StoppingCriterion, s::Stepsize, - dC::DirectionUpdateRule, + dC::DirectionUpdateRuleStorage, retr::AbstractRetractionMethod=ExponentialRetraction(), vtr::AbstractVectorTransportMethod=ParallelTransport(), initial_gradient::T=zero_vector(M, p), @@ -67,7 +86,7 @@ function ConjugateGradientDescentState( p::P, sC::StoppingCriterion, s::Stepsize, - dU::DirectionUpdateRule, + dU::DirectionUpdateRuleStorage, retr::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), vtr::AbstractVectorTransportMethod=default_vector_transport_method(M, typeof(p)), initial_gradient::T=zero_vector(M, p), @@ -100,26 +119,24 @@ Construct the conjugate descent coefficient update rule, a new storage is create > R. Fletcher, __Practical Methods of Optimization vol. 1: Unconstrained Optimization__ > John Wiley & Sons, New York, 1987. doi [10.1137/1024028](https://doi.org/10.1137/1024028) """ -mutable struct ConjugateDescentCoefficient <: DirectionUpdateRule - storage::StoreStateAction - function ConjugateDescentCoefficient( - a::StoreStateAction=StoreStateAction((:Iterate, :Gradient)) - ) - return new(a) - end -end -function (u::ConjugateDescentCoefficient)( +struct ConjugateDescentCoefficient <: DirectionUpdateRule end + +update_rule_storage_points(::ConjugateDescentCoefficient) = (:Iterate,) +update_rule_storage_vectors(::ConjugateDescentCoefficient) = (:Gradient,) + +function (u::DirectionUpdateRuleStorage{ConjugateDescentCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_storage(u.storage, :Iterate) || !has_storage(u.storage, :Gradient) - update_storage!(u.storage, cgs) # if not given store current as old + if !has_point_storage(u.storage, :Iterate) || !has_tangent_storage(u.storage, :Gradient) + update_storage!(u.storage, amp, cgs) # if not given store current as old return 0.0 end - p_old = get_storage(u.storage, :Iterate) - X_old = get_storage(u.storage, :Gradient) - update_storage!(u.storage, cgs) - return inner(M, cgs.p, cgs.X, cgs.X) / inner(M, p_old, -cgs.δ, X_old) + p_old = get_point_storage(u.storage, :Iterate) + X_old = get_tangent_storage(u.storage, :Gradient) + coef = inner(M, cgs.p, cgs.X, cgs.X) / inner(M, p_old, -cgs.δ, X_old) + update_storage!(u.storage, amp, cgs) + return coef end @doc raw""" @@ -157,36 +174,36 @@ default vector transport and a new storage is created by default. > SIAM J. Optim., 10 (1999), pp. 177–182. > doi: [10.1137/S1052623497318992](https://doi.org/10.1137/S1052623497318992) """ -mutable struct DaiYuanCoefficient{TVTM<:AbstractVectorTransportMethod} <: - DirectionUpdateRule +struct DaiYuanCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM - storage::StoreStateAction - function DaiYuanCoefficient( - t::AbstractVectorTransportMethod=ParallelTransport(), - a::StoreStateAction=StoreStateAction((:Iterate, :Gradient, :δ)), - ) - return new{typeof(t)}(t, a) + function DaiYuanCoefficient(t::AbstractVectorTransportMethod=ParallelTransport()) + return new{typeof(t)}(t) end end -function (u::DaiYuanCoefficient)( + +update_rule_storage_points(::DaiYuanCoefficient) = (:Iterate,) +update_rule_storage_vectors(::DaiYuanCoefficient) = (:Gradient, :δ) + +function (u::DirectionUpdateRuleStorage{<:DaiYuanCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_storage(u.storage, :Iterate) || - !has_storage(u.storage, :Gradient) || - !has_storage(u.storage, :δ) - update_storage!(u.storage, cgs) # if not given store current as old + if !has_point_storage(u.storage, :Iterate) || + !has_tangent_storage(u.storage, :Gradient) || + !has_tangent_storage(u.storage, :δ) + update_storage!(u.storage, amp, cgs) # if not given store current as old return 0.0 end - p_old = get_storage(u.storage, :Iterate) - X_old = get_storage(u.storage, :Gradient) - δ_old = get_storage(u.storage, :δ) - update_storage!(u.storage, cgs) + p_old = get_point_storage(u.storage, :Iterate) + X_old = get_tangent_storage(u.storage, :Gradient) + δ_old = get_tangent_storage(u.storage, :δ) - gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.transport_method) + gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) ν = cgs.X - gradienttr #notation y from [HZ06] - δtr = vector_transport_to(M, p_old, δ_old, cgs.p, u.transport_method) - return inner(M, cgs.p, cgs.X, cgs.X) / inner(M, p_old, δtr, ν) + δtr = vector_transport_to(M, p_old, δ_old, cgs.p, u.coefficient.transport_method) + coef = inner(M, cgs.p, cgs.X, cgs.X) / inner(M, p_old, δtr, ν) + update_storage!(u.storage, amp, cgs) + return coef end @doc raw""" @@ -214,25 +231,23 @@ Construct the Fletcher Reeves coefficient update rule, a new storage is created > Comput. J., 7 (1964), pp. 149–154. > doi: [10.1093/comjnl/7.2.149](http://dx.doi.org/10.1093/comjnl/7.2.149) """ -mutable struct FletcherReevesCoefficient <: DirectionUpdateRule - storage::StoreStateAction - function FletcherReevesCoefficient( - a::StoreStateAction=StoreStateAction((:Iterate, :Gradient)) - ) - return new(a) - end -end -function (u::FletcherReevesCoefficient)( +struct FletcherReevesCoefficient <: DirectionUpdateRule end + +update_rule_storage_points(::FletcherReevesCoefficient) = (:Iterate,) +update_rule_storage_vectors(::FletcherReevesCoefficient) = (:Gradient,) + +function (u::DirectionUpdateRuleStorage{FletcherReevesCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_storage(u.storage, :Iterate) || !has_storage(u.storage, :Gradient) - update_storage!(u.storage, cgs) # if not given store current as old + if !has_point_storage(u.storage, :Iterate) || !has_tangent_storage(u.storage, :Gradient) + update_storage!(u.storage, amp, cgs) # if not given store current as old end - p_old = get_storage(u.storage, :Iterate) - X_old = get_storage(u.storage, :Gradient) - update_storage!(u.storage, cgs) - return inner(M, cgs.p, cgs.X, cgs.X) / inner(M, p_old, X_old, X_old) + p_old = get_point_storage(u.storage, :Iterate) + X_old = get_tangent_storage(u.storage, :Gradient) + coef = inner(M, cgs.p, cgs.X, cgs.X) / inner(M, p_old, X_old, X_old) + update_storage!(u.storage, amp, cgs) + return coef end @doc raw""" @@ -274,32 +289,31 @@ default vector transport and a new storage is created by default. mutable struct HagerZhangCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM - storage::StoreStateAction - function HagerZhangCoefficient( - t::AbstractVectorTransportMethod=ParallelTransport(), - a::StoreStateAction=StoreStateAction((:Iterate, :Gradient, :δ)), - ) - return new{typeof(t)}(t, a) + function HagerZhangCoefficient(t::AbstractVectorTransportMethod=ParallelTransport()) + return new{typeof(t)}(t) end end -function (u::HagerZhangCoefficient)( + +update_rule_storage_points(::HagerZhangCoefficient) = (:Iterate,) +update_rule_storage_vectors(::HagerZhangCoefficient) = (:Gradient, :δ) + +function (u::DirectionUpdateRuleStorage{<:HagerZhangCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_storage(u.storage, :Iterate) || - !has_storage(u.storage, :Gradient) || - !has_storage(u.storage, :δ) - update_storage!(u.storage, cgs) # if not given store current as old + if !has_point_storage(u.storage, :Iterate) || + !has_tangent_storage(u.storage, :Gradient) || + !has_tangent_storage(u.storage, :δ) + update_storage!(u.storage, amp, cgs) # if not given store current as old return 0.0 end - p_old = get_storage(u.storage, :Iterate) - X_old = get_storage(u.storage, :Gradient) - δ_old = get_storage(u.storage, :δ) - update_storage!(u.storage, cgs) + p_old = get_point_storage(u.storage, :Iterate) + X_old = get_tangent_storage(u.storage, :Gradient) + δ_old = get_tangent_storage(u.storage, :δ) - gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.transport_method) + gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) ν = cgs.X - gradienttr #notation y from [HZ06] - δtr = vector_transport_to(M, p_old, δ_old, cgs.p, u.transport_method) + δtr = vector_transport_to(M, p_old, δ_old, cgs.p, u.coefficient.transport_method) denom = inner(M, cgs.p, δtr, ν) νknormsq = inner(M, cgs.p, ν, ν) β = @@ -308,7 +322,9 @@ function (u::HagerZhangCoefficient)( # Numerical stability from Manopt / Hager-Zhang paper ξn = norm(M, cgs.p, cgs.X) η = -1 / (ξn * min(0.01, norm(M, p_old, X_old))) - return max(β, η) + coef = max(β, η) + update_storage!(u.storage, amp, cgs) + return coef end @doc raw""" @@ -346,35 +362,38 @@ See also [`conjugate_gradient_descent`](@ref) > J. Research Nat. Bur. Standards, 49 (1952), pp. 409–436. > doi: [10.6028/jres.049.044](http://dx.doi.org/10.6028/jres.049.044) """ -mutable struct HestenesStiefelCoefficient{TVTM<:AbstractVectorTransportMethod} <: - DirectionUpdateRule +struct HestenesStiefelCoefficient{TVTM<:AbstractVectorTransportMethod} <: + DirectionUpdateRule transport_method::TVTM - storage::StoreStateAction function HestenesStiefelCoefficient( - transport_method::AbstractVectorTransportMethod=ParallelTransport(), - storage_action::StoreStateAction=StoreStateAction((:Iterate, :Gradient, :δ)), + transport_method::AbstractVectorTransportMethod=ParallelTransport() ) - return new{typeof(transport_method)}(transport_method, storage_action) + return new{typeof(transport_method)}(transport_method) end end -function (u::HestenesStiefelCoefficient)( + +update_rule_storage_points(::HestenesStiefelCoefficient) = (:Iterate,) +update_rule_storage_vectors(::HestenesStiefelCoefficient) = (:Gradient, :δ) + +function (u::DirectionUpdateRuleStorage{<:HestenesStiefelCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_storage(u.storage, :Iterate) || - !has_storage(u.storage, :Gradient) || - !has_storage(u.storage, :δ) - update_storage!(u.storage, cgs) # if not given store current as old + if !has_point_storage(u.storage, :Iterate) || + !has_tangent_storage(u.storage, :Gradient) || + !has_tangent_storage(u.storage, :δ) + update_storage!(u.storage, amp, cgs) # if not given store current as old return 0.0 end - p_old = get_storage(u.storage, :Iterate) - X_old = get_storage(u.storage, :Gradient) - δ_old = get_storage(u.storage, :δ) - update_storage!(u.storage, cgs) - gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.transport_method) - δtr = vector_transport_to(M, p_old, δ_old, cgs.p, u.transport_method) + p_old = get_point_storage(u.storage, :Iterate) + X_old = get_tangent_storage(u.storage, :Gradient) + δ_old = get_tangent_storage(u.storage, :δ) + + gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) + δtr = vector_transport_to(M, p_old, δ_old, cgs.p, u.coefficient.transport_method) ν = cgs.X - gradienttr #notation from [HZ06] β = inner(M, cgs.p, cgs.X, ν) / inner(M, cgs.p, δtr, ν) + update_storage!(u.storage, amp, cgs) return max(0, β) end @@ -414,33 +433,33 @@ default vector transport and a new storage is created by default. > J. Optim. Theory Appl., 69 (1991), pp. 129–137. > doi: [10.1007/BF00940464](https://doi.org/10.1007/BF00940464) """ -mutable struct LiuStoreyCoefficient{TVTM<:AbstractVectorTransportMethod} <: - DirectionUpdateRule +struct LiuStoreyCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM - storage::StoreStateAction - function LiuStoreyCoefficient( - t::AbstractVectorTransportMethod=ParallelTransport(), - a::StoreStateAction=StoreStateAction((:Iterate, :Gradient, :δ)), - ) - return new{typeof(t)}(t, a) + function LiuStoreyCoefficient(t::AbstractVectorTransportMethod=ParallelTransport()) + return new{typeof(t)}(t) end end -function (u::LiuStoreyCoefficient)( + +update_rule_storage_points(::LiuStoreyCoefficient) = (:Iterate,) +update_rule_storage_vectors(::LiuStoreyCoefficient) = (:Gradient, :δ) + +function (u::DirectionUpdateRuleStorage{<:LiuStoreyCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_storage(u.storage, :Iterate) || - !has_storage(u.storage, :Gradient) || - !has_storage(u.storage, :δ) - update_storage!(u.storage, cgs) # if not given store current as old + if !has_point_storage(u.storage, :Iterate) || + !has_tangent_storage(u.storage, :Gradient) || + !has_tangent_storage(u.storage, :δ) + update_storage!(u.storage, amp, cgs) # if not given store current as old end - p_old = get_storage(u.storage, :Iterate) - X_old = get_storage(u.storage, :Gradient) - δ_old = get_storage(u.storage, :δ) - update_storage!(u.storage, cgs) - gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.transport_method) + p_old = get_point_storage(u.storage, :Iterate) + X_old = get_tangent_storage(u.storage, :Gradient) + δ_old = get_tangent_storage(u.storage, :δ) + gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) ν = cgs.X - gradienttr # notation y from [HZ06] - return inner(M, cgs.p, cgs.X, ν) / inner(M, p_old, -δ_old, X_old) + coef = inner(M, cgs.p, cgs.X, ν) / inner(M, p_old, -δ_old, X_old) + update_storage!(u.storage, amp, cgs) + return coef end @doc raw""" @@ -485,31 +504,30 @@ See also [`conjugate_gradient_descent`](@ref) > USSR Comp. Math. Math. Phys., 9 (1969), pp. 94–112. > doi: [10.1016/0041-5553(69)90035-4](https://doi.org/10.1016/0041-5553(69)90035-4) """ -mutable struct PolakRibiereCoefficient{TVTM<:AbstractVectorTransportMethod} <: - DirectionUpdateRule +struct PolakRibiereCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM - storage::StoreStateAction - function PolakRibiereCoefficient( - t::AbstractVectorTransportMethod=ParallelTransport(), - a::StoreStateAction=StoreStateAction((:Iterate, :Gradient)), - ) - return new{typeof(t)}(t, a) + function PolakRibiereCoefficient(t::AbstractVectorTransportMethod=ParallelTransport()) + return new{typeof(t)}(t) end end -function (u::PolakRibiereCoefficient)( + +update_rule_storage_points(::PolakRibiereCoefficient) = (:Iterate,) +update_rule_storage_vectors(::PolakRibiereCoefficient) = (:Gradient,) + +function (u::DirectionUpdateRuleStorage{<:PolakRibiereCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_storage(u.storage, :Iterate) || !has_storage(u.storage, :Gradient) - update_storage!(u.storage, cgs) # if not given store current as old + if !has_point_storage(u.storage, :Iterate) || !has_tangent_storage(u.storage, :Gradient) + update_storage!(u.storage, amp, cgs) # if not given store current as old end - p_old = get_storage(u.storage, :Iterate) - X_old = get_storage(u.storage, :Gradient) - update_storage!(u.storage, cgs) + p_old = get_point_storage(u.storage, :Iterate) + X_old = get_tangent_storage(u.storage, :Gradient) - gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.transport_method) + gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) ν = cgs.X - gradienttr β = inner(M, cgs.p, cgs.X, ν) / inner(M, p_old, X_old, X_old) + update_storage!(u.storage, amp, cgs) return max(0, β) end @@ -522,7 +540,11 @@ hence return an update ``β = 0`` for all [`ConjugateGradientDescentState`](@ref See also [`conjugate_gradient_descent`](@ref) """ struct SteepestDirectionUpdateRule <: DirectionUpdateRule end -function (u::SteepestDirectionUpdateRule)( + +update_rule_storage_points(::SteepestDirectionUpdateRule) = () +update_rule_storage_vectors(::SteepestDirectionUpdateRule) = () + +function (u::DirectionUpdateRuleStorage{SteepestDirectionUpdateRule})( ::DefaultManoptProblem, ::ConjugateGradientDescentState, i ) return 0.0 @@ -577,7 +599,6 @@ mutable struct ConjugateGradientBealeRestart{ DUR<:DirectionUpdateRule,VT<:AbstractVectorTransportMethod,F } <: DirectionUpdateRule direction_update::DUR - storage::StoreStateAction threshold::F vector_transport_method::VT function ConjugateGradientBealeRestart( @@ -585,30 +606,39 @@ mutable struct ConjugateGradientBealeRestart{ threshold=0.2; manifold::AbstractManifold=DefaultManifold(), vector_transport_method::V=default_vector_transport_method(manifold), - a::StoreStateAction=StoreStateAction((:Iterate, :Gradient, :δ)), ) where {D<:DirectionUpdateRule,V<:AbstractVectorTransportMethod} return new{D,V,typeof(threshold)}( - direction_update, a, threshold, vector_transport_method + direction_update, threshold, vector_transport_method ) end end -function (u::ConjugateGradientBealeRestart)( + +function update_rule_storage_points(dur::ConjugateGradientBealeRestart) + return update_rule_storage_points(dur.direction_update) +end +function update_rule_storage_vectors(dur::ConjugateGradientBealeRestart) + return update_rule_storage_vectors(dur.direction_update) +end + +function (u::DirectionUpdateRuleStorage{<:ConjugateGradientBealeRestart})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_storage(u.storage, :Iterate) || !has_storage(u.storage, :Gradient) - update_storage!(u.storage, cgs) # if not given store current as old + if !has_point_storage(u.storage, :Iterate) || !has_tangent_storage(u.storage, :Gradient) + update_storage!(u.storage, amp, cgs) # if not given store current as old end - p_old = get_storage(u.storage, :Iterate) - X_old = get_storage(u.storage, :Gradient) + p_old = get_point_storage(u.storage, :Iterate) + X_old = get_tangent_storage(u.storage, :Gradient) # call actual rule β = u.direction_update(amp, cgs, i) - # update storage only after that in case they share - update_storage!(u.storage, cgs) denom = norm(M, cgs.p, cgs.X) - Xoldpk = vector_transport_to(M, p_old, X_old, cgs.p, u.vector_transport_method) + Xoldpk = vector_transport_to( + M, p_old, X_old, cgs.p, u.coefficient.vector_transport_method + ) num = inner(M, cgs.p, cgs.X, Xoldpk) + # update storage only after that in case they share + update_storage!(u.storage, amp, cgs) return (num / denom) > u.threshold ? zero(β) : β end diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index 8bc454cef1..b4a120ac9a 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -211,15 +211,40 @@ Initialize the Functor to a set of keys, where the dictionary is initialized to be empty. Further, `once` determines whether more that one update per iteration are effective, otherwise only the first update is stored, all others are ignored. """ -mutable struct StoreStateAction <: AbstractStateAction - values::Dict{Symbol,<:Any} - keys::NTuple{N,Symbol} where {N} +mutable struct StoreStateAction{TPS,TXS,TPI,TTI} <: AbstractStateAction + values::Dict{Symbol,Any} + keys::Vector{Symbol} # for values + point_values::TPS + tangent_values::TXS + point_init::TPI + tangent_init::TTI once::Bool last_stored::Int function StoreStateAction( - keys::NTuple{N,Symbol} where {N}=NTuple{0,Symbol}(), once=true + general_keys::Vector{Symbol}=Symbol[], + point_values=NamedTuple(), + tangent_values=NamedTuple(), + once=true, ) - return new(Dict{Symbol,Any}(), keys, once, -1) + point_init = NamedTuple{keys(point_values)}(map(u -> false, keys(point_values))) + tangent_init = NamedTuple{keys(tangent_values)}( + map(u -> false, keys(tangent_values)) + ) + return new{ + typeof(point_values), + typeof(tangent_values), + typeof(point_init), + typeof(tangent_init), + }( + Dict{Symbol,Any}(), + general_keys, + point_values, + tangent_values, + point_init, + tangent_init, + once, + -1, + ) end end function (a::StoreStateAction)( @@ -227,45 +252,67 @@ function (a::StoreStateAction)( ) #update values (maybe only once) if !a.once || a.last_stored != i - for key in a.keys - if key === :Iterate - M = get_manifold(amp) - a.values[key] = copy(M, get_iterate(s)) - elseif key === :Gradient - M = get_manifold(amp) - p = get_iterate(s) - a.values[key] = copy(M, p, get_gradient(s)) - elseif hasproperty(s, key) - a.values[key] = deepcopy(getproperty(s, key)) - end - end + update_storage!(a, amp, s) end return a.last_stored = i end """ - get_storage(a,key) + get_storage(a::AbstractStateAction, key::Symbol) + +Return the internal value of the [`AbstractStateAction`](@ref) `a` at the +`Symbol` `key`. +""" +get_storage(a::AbstractStateAction, key::Symbol) = a.values[key] + +""" + get_point_storage(a::AbstractStateAction, key::Symbol) + +Return the internal value of the [`AbstractStateAction`](@ref) `a` at the +`Symbol` `key` that represents a point. +""" +get_point_storage(a::AbstractStateAction, key::Symbol) = a.point_values[key] + +""" + get_tangent_storage(a::AbstractStateAction, key::Symbol) + +Return the internal value of the [`AbstractStateAction`](@ref) `a` at the +`Symbol` `key` that represents a tangent vector. +""" +get_tangent_storage(a::AbstractStateAction, key::Symbol) = a.tangent_values[key] -return the internal value of the [`AbstractStateAction`](@ref) `a` at the +""" + has_storage(a::AbstractStateAction, key::Symbol) + +Return whether the [`AbstractStateAction`](@ref) `a` has a value stored at the +`Symbol` `key`. +""" +has_storage(a::AbstractStateAction, key::Symbol) = haskey(a.values, key) + +""" + has_point_storage(a::AbstractStateAction, key::Symbol) + +Return whether the [`AbstractStateAction`](@ref) `a` has a point value stored at the `Symbol` `key`. """ -get_storage(a::AbstractStateAction, key) = a.values[key] +has_point_storage(a::AbstractStateAction, key::Symbol) = a.point_init[key] """ - get_storage(a,key) + has_tangent_storage(a::AbstractStateAction, key::Symbol) -return whether the [`AbstractStateAction`](@ref) `a` has a value stored at the +Return whether the [`AbstractStateAction`](@ref) `a` has a point value stored at the `Symbol` `key`. """ -has_storage(a::AbstractStateAction, key) = haskey(a.values, key) +has_tangent_storage(a::AbstractStateAction, key::Symbol) = a.tangent_init[key] """ - update_storage!(a, s) + update_storage!(a::AbstractStateAction, s::AbstractManoptSolverState) -update the [`AbstractStateAction`](@ref) `a` internal values to the ones given on +Update the [`AbstractStateAction`](@ref) `a` internal values to the ones given on the [`AbstractManoptSolverState`](@ref) `s`. """ function update_storage!(a::AbstractStateAction, s::AbstractManoptSolverState) + error("FFDF") for key in a.keys if key === :Iterate a.values[key] = deepcopy(get_iterate(s)) @@ -278,8 +325,55 @@ function update_storage!(a::AbstractStateAction, s::AbstractManoptSolverState) return a.keys end +function _storage_key_true(nt::NamedTuple) + return map(key -> NamedTuple{(key,),Tuple{Bool}}(true), keys(nt)) +end + +""" + update_storage!(a::AbstractStateAction, amp::AbstractManoptProblem, s::AbstractManoptSolverState) + +Update the [`AbstractStateAction`](@ref) `a` internal values to the ones given on +the [`AbstractManoptSolverState`](@ref) `s`. +Optimized using the information from `amp` +""" +function update_storage!( + a::AbstractStateAction, amp::AbstractManoptProblem, s::AbstractManoptSolverState +) + for key in a.keys + if key === :Iterate + a.values[key] = deepcopy(get_iterate(s)) + elseif key === :Gradient + a.values[key] = deepcopy(get_gradient(s)) + else + a.values[key] = deepcopy(getproperty(s, key)) + end + end + + M = get_manifold(amp) + + pt_kts = _storage_key_true(a.point_values) + map(keys(a.point_values), pt_kts) do key, kt + if key === :Iterate + copyto!(M, a.point_values[key], get_iterate(s)) + else + copyto!(M, a.point_values[key], getproperty(s, key)) + end + a.point_init = merge(a.point_init, kt) + end + tv_kts = _storage_key_true(a.tangent_values) + map(keys(a.tangent_values), tv_kts) do key, kt + if key === :Gradient + copyto!(M, a.tangent_values[key], get_gradient(s)) + else + copyto!(M, a.tangent_values[key], getproperty(s, key)) + end + a.tangent_init = merge(a.tangent_init, kt) + end + return a.keys +end + """ - update_storage!(a, d) + update_storage!(a::AbstractStateAction, d::Dict{Symbol,<:Any}) Update the [`AbstractStateAction`](@ref) `a` internal values to the ones given in the dictionary `d`. The values are merged, where the values from `d` are preferred. diff --git a/src/plans/stepsize.jl b/src/plans/stepsize.jl index d010663437..f7a5ebe4cf 100644 --- a/src/plans/stepsize.jl +++ b/src/plans/stepsize.jl @@ -266,24 +266,24 @@ function linesearch_backtrack( f0=F(x); stop_step=0.0, ) where {TF,T} - xNew = retract(M, x, s * η, retr) + xNew = retract(M, x, η, s, retr) fNew = F(xNew) search_dir_inner = real(inner(M, x, η, gradFx)) extended = false while fNew < f0 + decrease * s * search_dir_inner # increase extended = true s = s / contract - retract!(M, xNew, x, s * η, retr) + retract!(M, xNew, x, η, s, retr) fNew = F(xNew) end if extended s *= contract # undo last increase - retract!(M, xNew, x, s * η, retr) + retract!(M, xNew, x, η, s, retr) fNew = F(xNew) end while fNew > f0 + decrease * s * search_dir_inner # decrease s = contract * s - retract!(M, xNew, x, s * η, retr) + retract!(M, xNew, x, η, s, retr) fNew = F(xNew) (s < stop_step) && break end @@ -631,7 +631,7 @@ function (a::WolfePowellLinesearch)( s_minus = step f0 = get_cost(mp, cur_p) - p_new = retract(M, cur_p, step * η, a.retraction_method) + p_new = retract(M, cur_p, η, step, a.retraction_method) fNew = get_cost(mp, p_new) η_xNew = vector_transport_to(M, cur_p, η, p_new, a.vector_transport_method) if fNew > f0 + a.c1 * step * real(inner(M, get_iterate(ams), η, get_gradient(ams))) @@ -640,7 +640,7 @@ function (a::WolfePowellLinesearch)( ) && (s_minus > 10^(-9)) # decrease s_minus = s_minus * 0.5 step = s_minus - retract!(M, p_new, get_iterate(ams), step * η, a.retraction_method) + retract!(M, p_new, get_iterate(ams), η, step, a.retraction_method) fNew = get_cost(mp, p_new) end s_plus = 2.0 * s_minus @@ -656,18 +656,18 @@ function (a::WolfePowellLinesearch)( (s_plus < max_step_increase)# increase s_plus = s_plus * 2.0 step = s_plus - retract!(M, p_new, get_iterate(ams), step * η, a.retraction_method) + retract!(M, p_new, get_iterate(ams), η, step, a.retraction_method) fNew = get_cost(mp, p_new) end s_minus = s_plus / 2.0 end end - retract!(M, p_new, get_iterate(ams), s_minus * η, a.retraction_method) + retract!(M, p_new, get_iterate(ams), η, s_minus, a.retraction_method) vector_transport_to!(M, η_xNew, get_iterate(ams), η, p_new, a.vector_transport_method) while real(inner(M, p_new, get_gradient(mp, p_new), η_xNew)) < a.c2 * real(inner(M, get_iterate(ams), η, get_gradient(ams))) step = (s_minus + s_plus) / 2 - retract!(M, p_new, get_iterate(ams), step * η, a.retraction_method) + retract!(M, p_new, get_iterate(ams), η, step, a.retraction_method) fNew = get_cost(mp, p_new) if fNew <= f0 + a.c1 * step * real(inner(M, get_iterate(ams), η, get_gradient(ams))) s_minus = step @@ -675,7 +675,7 @@ function (a::WolfePowellLinesearch)( s_plus = step end abs(s_plus - s_minus) <= a.linesearch_stopsize && break - retract!(M, p_new, get_iterate(ams), s_minus * η, a.retraction_method) + retract!(M, p_new, get_iterate(ams), η, s_minus, a.retraction_method) vector_transport_to!( M, η_xNew, get_iterate(ams), η, p_new, a.vector_transport_method ) @@ -775,7 +775,7 @@ function (a::WolfePowellBinaryLinesearch)( β = Inf t = 1.0 f0 = get_cost(amp, get_iterate(ams)) - xNew = retract(M, get_iterate(ams), t * η, a.retraction_method) + xNew = retract(M, get_iterate(ams), η, t, a.retraction_method) fNew = get_cost(amp, xNew) η_xNew = vector_transport_to(M, get_iterate(ams), η, xNew, a.vector_transport_method) gradient_new = get_gradient(amp, xNew) @@ -790,7 +790,7 @@ function (a::WolfePowellBinaryLinesearch)( (!nAt && nWt) && (α = t) # A(t) holds but W(t) fails t = isinf(β) ? 2 * α : (α + β) / 2 # Update trial point - retract!(M, xNew, get_iterate(ams), t * η, a.retraction_method) + retract!(M, xNew, get_iterate(ams), η, t, a.retraction_method) fNew = get_cost(amp, xNew) gradient_new = get_gradient(amp, xNew) vector_transport_to!( diff --git a/src/solvers/conjugate_gradient_descent.jl b/src/solvers/conjugate_gradient_descent.jl index 0ac0ccb952..0278766bf5 100644 --- a/src/solvers/conjugate_gradient_descent.jl +++ b/src/solvers/conjugate_gradient_descent.jl @@ -107,7 +107,7 @@ function conjugate_gradient_descent!( p, stopping_criterion, stepsize, - coefficient, + DirectionUpdateRuleStorage(M, coefficient, p), retraction_method, vector_transport_method, X, @@ -127,10 +127,11 @@ function step_solver!(amp::AbstractManoptProblem, cgs::ConjugateGradientDescentS M = get_manifold(amp) p_old = copy(M, cgs.p) current_stepsize = get_stepsize(amp, cgs, i, cgs.δ) - retract!(M, cgs.p, cgs.p, current_stepsize * cgs.δ, cgs.retraction_method) + retract!(M, cgs.p, cgs.p, cgs.δ, current_stepsize, cgs.retraction_method) get_gradient!(amp, cgs.X, cgs.p) cgs.β = cgs.coefficient(amp, cgs, i) vector_transport_to!(M, cgs.δ, p_old, cgs.δ, cgs.p, cgs.vector_transport_method) - cgs.δ .= -cgs.X .+ cgs.β * cgs.δ + cgs.δ .*= cgs.β + cgs.δ .-= cgs.X return cgs end diff --git a/test/solvers/test_conjugate_gradient.jl b/test/solvers/test_conjugate_gradient.jl index 39a7eaec9f..cea104d468 100644 --- a/test/solvers/test_conjugate_gradient.jl +++ b/test/solvers/test_conjugate_gradient.jl @@ -1,3 +1,4 @@ +using Revise using Manopt, Manifolds, ManifoldsBase, Test using LinearAlgebra: Diagonal, dot, eigvals, eigvecs @@ -18,12 +19,12 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs δ2 = [0.5, 2.0] diff = grad_2 - grad_1 - dU = SteepestDirectionUpdateRule() + dU = Manopt.DirectionUpdateRuleStorage(M, SteepestDirectionUpdateRule()) cgs = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm, zero_vector(M, x0)) @test cgs.coefficient(dmp, cgs, 1) == 0 @test default_stepsize(M, typeof(cgs)) isa ArmijoLinesearch - dU = ConjugateDescentCoefficient() + dU = Manopt.DirectionUpdateRuleStorage(M, ConjugateDescentCoefficient()) cgs = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm, zero_vector(M, x0)) cgs.X = grad_1 cgs.δ = δ1 @@ -33,7 +34,7 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs cgs.δ = δ2 @test cgs.coefficient(dmp, cgs, 2) == dot(grad_2, grad_2) / dot(-δ2, grad_1) - dU = DaiYuanCoefficient() + dU = Manopt.DirectionUpdateRuleStorage(M, DaiYuanCoefficient()) cgs = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm) cgs.X = grad_1 cgs.δ = δ1 @@ -43,7 +44,7 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs cgs.δ = δ2 @test cgs.coefficient(dmp, cgs, 2) == dot(grad_2, grad_2) / dot(δ2, grad_2 - grad_1) - dU = FletcherReevesCoefficient() + dU = Manopt.DirectionUpdateRuleStorage(M, FletcherReevesCoefficient()) cgs = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm) cgs.X = grad_1 cgs.δ = δ1 @@ -53,7 +54,7 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs cgs.δ = δ2 @test cgs.coefficient(dmp, cgs, 2) == dot(grad_2, grad_2) / dot(grad_1, grad_1) - dU = HagerZhangCoefficient() + dU = Manopt.DirectionUpdateRuleStorage(M, HagerZhangCoefficient()) cgs = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm) cgs.X = grad_1 cgs.δ = δ1 @@ -66,7 +67,7 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs @test cgs.coefficient(dmp, cgs, 2) == dot(diff, grad_2) / denom - 2 * ndiffsq * dot(δ1, grad_2) / denom^2 - dU = HestenesStiefelCoefficient() + dU = Manopt.DirectionUpdateRuleStorage(M, HestenesStiefelCoefficient()) O = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm) O.X = grad_1 O.δ = δ1 @@ -75,7 +76,7 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs O.δ = δ2 @test O.coefficient(dmp, O, 2) == dot(diff, grad_2) / dot(δ1, diff) - dU = LiuStoreyCoefficient() + dU = Manopt.DirectionUpdateRuleStorage(M, LiuStoreyCoefficient()) O = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm) O.X = grad_1 O.δ = δ1 @@ -84,7 +85,7 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs O.δ = δ2 @test O.coefficient(dmp, O, 2) == -dot(grad_2, diff) / dot(δ1, grad_1) - dU = PolakRibiereCoefficient() + dU = Manopt.DirectionUpdateRuleStorage(M, PolakRibiereCoefficient()) O = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm) O.X = grad_1 O.δ = δ1 From cf192bd3e33faca95eb8ebe015761c1509c28065 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Tue, 31 Jan 2023 22:06:17 +0100 Subject: [PATCH 02/37] minor improvement --- src/solvers/conjugate_gradient_descent.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/conjugate_gradient_descent.jl b/src/solvers/conjugate_gradient_descent.jl index 0278766bf5..1a5e650f78 100644 --- a/src/solvers/conjugate_gradient_descent.jl +++ b/src/solvers/conjugate_gradient_descent.jl @@ -107,7 +107,7 @@ function conjugate_gradient_descent!( p, stopping_criterion, stepsize, - DirectionUpdateRuleStorage(M, coefficient, p), + DirectionUpdateRuleStorage(M, coefficient, p, X), retraction_method, vector_transport_method, X, From 2a319d8a2b5970bcdf9e5083343cbd0028fe806b Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Tue, 31 Jan 2023 22:15:59 +0100 Subject: [PATCH 03/37] cleanup and asserts --- src/plans/solver_state.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index b4a120ac9a..f188a319ff 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -312,7 +312,6 @@ Update the [`AbstractStateAction`](@ref) `a` internal values to the ones given o the [`AbstractManoptSolverState`](@ref) `s`. """ function update_storage!(a::AbstractStateAction, s::AbstractManoptSolverState) - error("FFDF") for key in a.keys if key === :Iterate a.values[key] = deepcopy(get_iterate(s)) @@ -345,7 +344,7 @@ function update_storage!( elseif key === :Gradient a.values[key] = deepcopy(get_gradient(s)) else - a.values[key] = deepcopy(getproperty(s, key)) + qa.values[key] = deepcopy(getproperty(s, key)) end end @@ -356,7 +355,7 @@ function update_storage!( if key === :Iterate copyto!(M, a.point_values[key], get_iterate(s)) else - copyto!(M, a.point_values[key], getproperty(s, key)) + copyto!(M, a.point_values[key], getproperty(s, key)::typeof(a.point_values[key])) end a.point_init = merge(a.point_init, kt) end @@ -365,7 +364,7 @@ function update_storage!( if key === :Gradient copyto!(M, a.tangent_values[key], get_gradient(s)) else - copyto!(M, a.tangent_values[key], getproperty(s, key)) + copyto!(M, a.tangent_values[key], getproperty(s, key)::typeof(a.tangent_values[key])) end a.tangent_init = merge(a.tangent_init, kt) end From dba15cbcadfb187df3a6f614f02ef43267a0f4fb Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 1 Feb 2023 16:14:55 +0100 Subject: [PATCH 04/37] fixing things and addressing code review --- src/plans/conjugate_gradient_plan.jl | 19 +++++++++++-------- src/plans/debug.jl | 12 ++++++------ src/plans/record.jl | 8 ++++---- src/plans/solver_state.jl | 18 +++++++++++------- src/plans/stepsize.jl | 14 ++++++++++---- src/plans/stopping_criterion.jl | 12 ++++++++---- src/solvers/conjugate_gradient_descent.jl | 2 +- test/plans/test_conjugate_gradient_plan.jl | 2 ++ test/plans/test_debug.jl | 8 ++++---- test/plans/test_primal_dual_plan.jl | 4 ++-- test/plans/test_record.jl | 2 +- test/solvers/test_conjugate_gradient.jl | 17 ++++++++--------- 12 files changed, 68 insertions(+), 50 deletions(-) diff --git a/src/plans/conjugate_gradient_plan.jl b/src/plans/conjugate_gradient_plan.jl index 9cdf017a5a..22dfee2596 100644 --- a/src/plans/conjugate_gradient_plan.jl +++ b/src/plans/conjugate_gradient_plan.jl @@ -62,20 +62,21 @@ mutable struct ConjugateGradientDescentState{ p::P, sC::StoppingCriterion, s::Stepsize, - dC::DirectionUpdateRuleStorage, + dC::DirectionUpdateRule, retr::AbstractRetractionMethod=ExponentialRetraction(), vtr::AbstractVectorTransportMethod=ParallelTransport(), initial_gradient::T=zero_vector(M, p), ) where {P,T} + coef = DirectionUpdateRuleStorage(M, dC, p, initial_gradient) βT = allocate_result_type(M, ConjugateGradientDescentState, (p, initial_gradient)) - cgs = new{P,T,βT,typeof(dC),typeof(s),typeof(sC),typeof(retr),typeof(vtr)}() + cgs = new{P,T,βT,typeof(coef),typeof(s),typeof(sC),typeof(retr),typeof(vtr)}() cgs.p = p cgs.X = initial_gradient cgs.δ = copy(M, p, initial_gradient) cgs.stop = sC cgs.retraction_method = retr cgs.stepsize = s - cgs.coefficient = dC + cgs.coefficient = coef cgs.vector_transport_method = vtr cgs.β = zero(βT) return cgs @@ -86,7 +87,7 @@ function ConjugateGradientDescentState( p::P, sC::StoppingCriterion, s::Stepsize, - dU::DirectionUpdateRuleStorage, + dU::DirectionUpdateRule, retr::AbstractRetractionMethod=default_retraction_method(M, typeof(p)), vtr::AbstractVectorTransportMethod=default_vector_transport_method(M, typeof(p)), initial_gradient::T=zero_vector(M, p), @@ -614,10 +615,12 @@ mutable struct ConjugateGradientBealeRestart{ end function update_rule_storage_points(dur::ConjugateGradientBealeRestart) - return update_rule_storage_points(dur.direction_update) + dur_p = update_rule_storage_points(dur.direction_update) + return :Iterate in dur_p ? dur_p : (:Iterate, dur_p...) end function update_rule_storage_vectors(dur::ConjugateGradientBealeRestart) - return update_rule_storage_vectors(dur.direction_update) + dur_X = update_rule_storage_vectors(dur.direction_update) + return :Gradient in dur_X ? dur_X : (:Gradient, dur_X...) end function (u::DirectionUpdateRuleStorage{<:ConjugateGradientBealeRestart})( @@ -631,7 +634,7 @@ function (u::DirectionUpdateRuleStorage{<:ConjugateGradientBealeRestart})( X_old = get_tangent_storage(u.storage, :Gradient) # call actual rule - β = u.direction_update(amp, cgs, i) + β = u.coefficient.direction_update(amp, cgs, i) denom = norm(M, cgs.p, cgs.X) Xoldpk = vector_transport_to( @@ -640,5 +643,5 @@ function (u::DirectionUpdateRuleStorage{<:ConjugateGradientBealeRestart})( num = inner(M, cgs.p, cgs.X, Xoldpk) # update storage only after that in case they share update_storage!(u.storage, amp, cgs) - return (num / denom) > u.threshold ? zero(β) : β + return (num / denom) > u.coefficient.threshold ? zero(β) : β end diff --git a/src/plans/debug.jl b/src/plans/debug.jl index cad9dc940e..de1ace8f49 100644 --- a/src/plans/debug.jl +++ b/src/plans/debug.jl @@ -139,7 +139,7 @@ debug for the amount of change of the iterate (stored in `get_iterate(o)` of the during the last iteration. See [`DebugEntryChange`](@ref) for the general case # Keyword Parameters -* `storage` – (`StoreStateAction( (:Iterate,) )`) – (eventually shared) the storage of the previous action +* `storage` – (`StoreStateAction( [:Gradient] )`) – (eventually shared) the storage of the previous action * `prefix` – (`"Last Change:"`) prefix of the debug output (ignored if you set `format`) * `io` – (`stdout`) default steream to print the debug to. * `format` - ( `"$prefix %f"`) format to print the output using an sprintf format. @@ -154,7 +154,7 @@ mutable struct DebugChange{TInvRetr<:AbstractInverseRetractionMethod} <: DebugAc storage::StoreStateAction invretr::TInvRetr function DebugChange(; - storage::StoreStateAction=StoreStateAction((:Iterate,)), + storage::StoreStateAction=StoreStateAction([:Iterate]), io::IO=stdout, prefix::String="Last Change: ", format::String="$(prefix)%f", @@ -183,7 +183,7 @@ debug for the amount of change of the gradient (stored in `get_gradient(o)` of t during the last iteration. See [`DebugEntryChange`](@ref) for the general case # Keyword Parameters -* `storage` – (`StoreStateAction( (:Gradient,) )`) – (eventually shared) the storage of the previous action +* `storage` – (`StoreStateAction( [:Gradient] )`) – (eventually shared) the storage of the previous action * `prefix` – (`"Last Change:"`) prefix of the debug output (ignored if you set `format`) * `io` – (`stdout`) default steream to print the debug to. * `format` - ( `"$prefix %f"`) format to print the output using an sprintf format. @@ -193,7 +193,7 @@ mutable struct DebugGradientChange <: DebugAction format::String storage::StoreStateAction function DebugGradientChange(; - storage::StoreStateAction=StoreStateAction((:Gradient,)), + storage::StoreStateAction=StoreStateAction([:Iterate, :Gradient]), io::IO=stdout, prefix::String="Last Change: ", format::String="$(prefix)%f", @@ -379,11 +379,11 @@ mutable struct DebugEntryChange <: DebugAction function DebugEntryChange( f::Symbol, d; - storage::StoreStateAction=StoreStateAction((f,)), + storage::StoreStateAction=StoreStateAction([f]), prefix::String="Change of $f:", format::String="$prefix%s", io::IO=stdout, - initial_value::T where {T}=NaN, + initial_value::Any=NaN, ) if !isa(initial_value, Number) || !isnan(initial_value) #set initial value update_storage!(storage, Dict(f => initial_value)) diff --git a/src/plans/record.jl b/src/plans/record.jl index a8dc996684..308e55db4f 100644 --- a/src/plans/record.jl +++ b/src/plans/record.jl @@ -329,7 +329,7 @@ mutable struct RecordChange{TInvRetr<:AbstractInverseRetractionMethod} <: Record storage::StoreStateAction invretr::TInvRetr function RecordChange( - a::StoreStateAction=StoreStateAction((:Iterate,)); + a::StoreStateAction=StoreStateAction([:Iterate]); manifold::AbstractManifold=DefaultManifold(1), invretr::AbstractInverseRetractionMethod=default_inverse_retraction_method( manifold @@ -339,7 +339,7 @@ mutable struct RecordChange{TInvRetr<:AbstractInverseRetractionMethod} <: Record end function RecordChange( p, - a::StoreStateAction=StoreStateAction((:Iterate,)); + a::StoreStateAction=StoreStateAction([:Iterate]); manifold::AbstractManifold=DefaultManifold(1), invretr::AbstractInverseRetractionMethod=default_inverse_retraction_method( manifold, typeof(p) @@ -403,11 +403,11 @@ mutable struct RecordEntryChange <: RecordAction field::Symbol distance::Any storage::StoreStateAction - function RecordEntryChange(f::Symbol, d, a::StoreStateAction=StoreStateAction((f,))) + function RecordEntryChange(f::Symbol, d, a::StoreStateAction=StoreStateAction([f])) return new(Float64[], f, d, a) end function RecordEntryChange( - v::T where {T}, f::Symbol, d, a::StoreStateAction=StoreStateAction((f,)) + v::T where {T}, f::Symbol, d, a::StoreStateAction=StoreStateAction([f]) ) update_storage!(a, Dict(f => v)) return new(Float64[], f, d, a) diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index f188a319ff..50577388f8 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -222,9 +222,9 @@ mutable struct StoreStateAction{TPS,TXS,TPI,TTI} <: AbstractStateAction last_stored::Int function StoreStateAction( general_keys::Vector{Symbol}=Symbol[], - point_values=NamedTuple(), - tangent_values=NamedTuple(), - once=true, + point_values::NamedTuple=NamedTuple(), + tangent_values::NamedTuple=NamedTuple(), + once::Bool=true, ) point_init = NamedTuple{keys(point_values)}(map(u -> false, keys(point_values))) tangent_init = NamedTuple{keys(tangent_values)}( @@ -344,7 +344,7 @@ function update_storage!( elseif key === :Gradient a.values[key] = deepcopy(get_gradient(s)) else - qa.values[key] = deepcopy(getproperty(s, key)) + a.values[key] = deepcopy(getproperty(s, key)) end end @@ -355,7 +355,9 @@ function update_storage!( if key === :Iterate copyto!(M, a.point_values[key], get_iterate(s)) else - copyto!(M, a.point_values[key], getproperty(s, key)::typeof(a.point_values[key])) + copyto!( + M, a.point_values[key], getproperty(s, key)::typeof(a.point_values[key]) + ) end a.point_init = merge(a.point_init, kt) end @@ -364,7 +366,9 @@ function update_storage!( if key === :Gradient copyto!(M, a.tangent_values[key], get_gradient(s)) else - copyto!(M, a.tangent_values[key], getproperty(s, key)::typeof(a.tangent_values[key])) + copyto!( + M, a.tangent_values[key], getproperty(s, key)::typeof(a.tangent_values[key]) + ) end a.tangent_init = merge(a.tangent_init, kt) end @@ -380,5 +384,5 @@ the dictionary `d`. The values are merged, where the values from `d` are preferr function update_storage!(a::AbstractStateAction, d::Dict{Symbol,<:Any}) merge!(a.values, d) # update keys - return a.keys = Tuple(keys(a.values)) + return a.keys = collect(keys(a.values)) end diff --git a/src/plans/stepsize.jl b/src/plans/stepsize.jl index f7a5ebe4cf..334ea394a9 100644 --- a/src/plans/stepsize.jl +++ b/src/plans/stepsize.jl @@ -379,7 +379,10 @@ and vector transport are set to the default ones on `M`, repsectively. The constructors return the functor to perform nonmonotone line search. """ mutable struct NonmonotoneLinesearch{ - TRM<:AbstractRetractionMethod,VTM<:AbstractVectorTransportMethod,T<:AbstractVector + TRM<:AbstractRetractionMethod, + VTM<:AbstractVectorTransportMethod, + T<:AbstractVector, + TSSA<:StoreStateAction, } <: Linesearch retraction_method::TRM vector_transport_method::VTM @@ -390,7 +393,7 @@ mutable struct NonmonotoneLinesearch{ initial_stepsize::Float64 old_costs::T strategy::Symbol - storage::StoreStateAction + storage::TSSA linesearch_stopsize::Float64 function NonmonotoneLinesearch( M::AbstractManifold=DefaultManifold(2); @@ -405,7 +408,7 @@ mutable struct NonmonotoneLinesearch{ min_stepsize::Float64=1e-3, max_stepsize::Float64=1e3, strategy::Symbol=:direct, - storage::StoreStateAction=StoreStateAction((:Iterate, :Gradient)), + storage::StoreStateAction=StoreStateAction([:Iterate, :Gradient]), linesearch_stopsize::Float64=0.0, ) if strategy ∉ [:direct, :inverse, :alternating] @@ -436,7 +439,10 @@ mutable struct NonmonotoneLinesearch{ throw(DomainError(memory_size, "The memory_size has to be greater than zero.")) end return new{ - typeof(retraction_method),typeof(vector_transport_method),Vector{Float64} + typeof(retraction_method), + typeof(vector_transport_method), + Vector{Float64}, + typeof(storage), }( retraction_method, vector_transport_method, diff --git a/src/plans/stopping_criterion.jl b/src/plans/stopping_criterion.jl index 57d86a27d4..ddfcaf3f0d 100644 --- a/src/plans/stopping_criterion.jl +++ b/src/plans/stopping_criterion.jl @@ -131,19 +131,23 @@ initialize the stopping criterion to a threshold `ε` using the default. You can also provide an inverse_retraction_method for the `distance` or a manifold to use its default inverse retraction. """ -mutable struct StopWhenChangeLess{IRT} <: StoppingCriterion +mutable struct StopWhenChangeLess{ + IRT<:AbstractInverseRetractionMethod,TSSA<:StoreStateAction +} <: StoppingCriterion threshold::Float64 reason::String - storage::StoreStateAction + storage::TSSA inverse_retraction::IRT end function StopWhenChangeLess( ε::Float64; - storage::StoreStateAction=StoreStateAction((:Iterate,)), + storage::StoreStateAction=StoreStateAction([:Iterate]), manifold::AbstractManifold=DefaultManifold(3), inverse_retraction_method::IRT=default_inverse_retraction_method(manifold), ) where {IRT<:AbstractInverseRetractionMethod} - return StopWhenChangeLess{IRT}(ε, "", storage, inverse_retraction_method) + return StopWhenChangeLess{IRT,typeof(storage)}( + ε, "", storage, inverse_retraction_method + ) end function (c::StopWhenChangeLess)(mp::AbstractManoptProblem, s::AbstractManoptSolverState, i) diff --git a/src/solvers/conjugate_gradient_descent.jl b/src/solvers/conjugate_gradient_descent.jl index 1a5e650f78..70147e5ea2 100644 --- a/src/solvers/conjugate_gradient_descent.jl +++ b/src/solvers/conjugate_gradient_descent.jl @@ -107,7 +107,7 @@ function conjugate_gradient_descent!( p, stopping_criterion, stepsize, - DirectionUpdateRuleStorage(M, coefficient, p, X), + coefficient, retraction_method, vector_transport_method, X, diff --git a/test/plans/test_conjugate_gradient_plan.jl b/test/plans/test_conjugate_gradient_plan.jl index b4674bd18b..058c3507ee 100644 --- a/test/plans/test_conjugate_gradient_plan.jl +++ b/test/plans/test_conjugate_gradient_plan.jl @@ -2,6 +2,8 @@ using Manopt, Manifolds, Test struct DummyCGCoeff <: DirectionUpdateRule end (u::DummyCGCoeff)(p, s, i) = 0.2 +Manopt.update_rule_storage_points(::DummyCGCoeff) = () +Manopt.update_rule_storage_vectors(::DummyCGCoeff) = () @testset "Test Restart CG" begin M = Euclidean(2) diff --git a/test/plans/test_debug.jl b/test/plans/test_debug.jl index 2ebbffddfb..9d2343293a 100644 --- a/test/plans/test_debug.jl +++ b/test/plans/test_debug.jl @@ -49,18 +49,18 @@ end DebugEntry(:p; prefix="x:", io=io)(mp, st, -1) @test String(take!(io)) == "" # Change of Iterate - a2 = DebugChange(; storage=StoreStateAction((:Iterate,)), prefix="Last: ", io=io) + a2 = DebugChange(; storage=StoreStateAction([:Iterate]), prefix="Last: ", io=io) a2(mp, st, 0) # init st.p = [3.0, 2.0] a2(mp, st, 1) a2inv = DebugChange(; - storage=StoreStateAction((:Iterate,)), + storage=StoreStateAction([:Iterate]), prefix="Last: ", io=io, invretr=PolarInverseRetraction(), ) a2mani = DebugChange(; - storage=StoreStateAction((:Iterate,)), + storage=StoreStateAction([:Iterate]), prefix="Last: ", io=io, manifold=TestPolarManifold(), @@ -71,7 +71,7 @@ end @test String(take!(io)) == "Last: 1.000000" # Change of Gradient a3 = DebugGradientChange(; - storage=StoreStateAction((:Gradient,)), prefix="Last: ", io=io + storage=StoreStateAction([:Gradient]), prefix="Last: ", io=io ) a3(mp, st, 0) # init st.X = [1.0, 0.0] diff --git a/test/plans/test_primal_dual_plan.jl b/test/plans/test_primal_dual_plan.jl index 8e4c5edd81..c8dfae4f48 100644 --- a/test/plans/test_primal_dual_plan.jl +++ b/test/plans/test_primal_dual_plan.jl @@ -175,7 +175,7 @@ using Manopt, Manifolds, ManifoldsBase, Test @test_throws DomainError dual_residual(p_exact, o_err, p_old, ξ_old, n_old) end @testset "Debug prints" begin - a = StoreStateAction((:Iterate, :X, :n, :m)) + a = StoreStateAction([:Iterate, :X, :n, :m]) update_storage!(a, Dict(:Iterate => p_old, :X => ξ_old, :n => n_old, :m => copy(m))) io = IOBuffer() @@ -255,7 +255,7 @@ using Manopt, Manifolds, ManifoldsBase, Test @test startswith(s, "Primal Residual:") end @testset "Records" begin - a = StoreStateAction((:Iterate, :X, :n, :m)) + a = StoreStateAction([:Iterate, :X, :n, :m]) update_storage!(a, Dict(:Iterate => p_old, :X => ξ_old, :n => n_old, :m => copy(m))) io = IOBuffer() diff --git a/test/plans/test_record.jl b/test/plans/test_record.jl index c350cee495..77b255b2ef 100644 --- a/test/plans/test_record.jl +++ b/test/plans/test_record.jl @@ -127,7 +127,7 @@ using Manifolds, Manopt, Test, ManifoldsBase, Dates # RecordEntryChange set_iterate!(gds, M, p) e = RecordEntryChange(:p, (p, o, x, y) -> distance(get_manifold(p), x, y)) - @test update_storage!(e.storage, gds) == (:p,) + @test update_storage!(e.storage, gds) == [:p] e2 = RecordEntryChange(dmp, :p, (p, o, x, y) -> distance(get_manifold(p), x, y)) @test e.field == e2.field e(dmp, gds, 1) diff --git a/test/solvers/test_conjugate_gradient.jl b/test/solvers/test_conjugate_gradient.jl index cea104d468..39a7eaec9f 100644 --- a/test/solvers/test_conjugate_gradient.jl +++ b/test/solvers/test_conjugate_gradient.jl @@ -1,4 +1,3 @@ -using Revise using Manopt, Manifolds, ManifoldsBase, Test using LinearAlgebra: Diagonal, dot, eigvals, eigvecs @@ -19,12 +18,12 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs δ2 = [0.5, 2.0] diff = grad_2 - grad_1 - dU = Manopt.DirectionUpdateRuleStorage(M, SteepestDirectionUpdateRule()) + dU = SteepestDirectionUpdateRule() cgs = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm, zero_vector(M, x0)) @test cgs.coefficient(dmp, cgs, 1) == 0 @test default_stepsize(M, typeof(cgs)) isa ArmijoLinesearch - dU = Manopt.DirectionUpdateRuleStorage(M, ConjugateDescentCoefficient()) + dU = ConjugateDescentCoefficient() cgs = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm, zero_vector(M, x0)) cgs.X = grad_1 cgs.δ = δ1 @@ -34,7 +33,7 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs cgs.δ = δ2 @test cgs.coefficient(dmp, cgs, 2) == dot(grad_2, grad_2) / dot(-δ2, grad_1) - dU = Manopt.DirectionUpdateRuleStorage(M, DaiYuanCoefficient()) + dU = DaiYuanCoefficient() cgs = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm) cgs.X = grad_1 cgs.δ = δ1 @@ -44,7 +43,7 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs cgs.δ = δ2 @test cgs.coefficient(dmp, cgs, 2) == dot(grad_2, grad_2) / dot(δ2, grad_2 - grad_1) - dU = Manopt.DirectionUpdateRuleStorage(M, FletcherReevesCoefficient()) + dU = FletcherReevesCoefficient() cgs = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm) cgs.X = grad_1 cgs.δ = δ1 @@ -54,7 +53,7 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs cgs.δ = δ2 @test cgs.coefficient(dmp, cgs, 2) == dot(grad_2, grad_2) / dot(grad_1, grad_1) - dU = Manopt.DirectionUpdateRuleStorage(M, HagerZhangCoefficient()) + dU = HagerZhangCoefficient() cgs = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm) cgs.X = grad_1 cgs.δ = δ1 @@ -67,7 +66,7 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs @test cgs.coefficient(dmp, cgs, 2) == dot(diff, grad_2) / denom - 2 * ndiffsq * dot(δ1, grad_2) / denom^2 - dU = Manopt.DirectionUpdateRuleStorage(M, HestenesStiefelCoefficient()) + dU = HestenesStiefelCoefficient() O = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm) O.X = grad_1 O.δ = δ1 @@ -76,7 +75,7 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs O.δ = δ2 @test O.coefficient(dmp, O, 2) == dot(diff, grad_2) / dot(δ1, diff) - dU = Manopt.DirectionUpdateRuleStorage(M, LiuStoreyCoefficient()) + dU = LiuStoreyCoefficient() O = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm) O.X = grad_1 O.δ = δ1 @@ -85,7 +84,7 @@ using LinearAlgebra: Diagonal, dot, eigvals, eigvecs O.δ = δ2 @test O.coefficient(dmp, O, 2) == -dot(grad_2, diff) / dot(δ1, grad_1) - dU = Manopt.DirectionUpdateRuleStorage(M, PolakRibiereCoefficient()) + dU = PolakRibiereCoefficient() O = ConjugateGradientDescentState(M, x0, sC, s, dU, retr, vtm) O.X = grad_1 O.δ = δ1 From 284782a56f72aa6848bb2e46fe57f56f266e557f Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 1 Feb 2023 21:37:46 +0100 Subject: [PATCH 05/37] some general improvements and docs --- src/plans/conjugate_gradient_plan.jl | 9 +++--- src/plans/solver_state.jl | 38 ++++++++++++++++++++--- src/plans/stopping_criterion.jl | 2 +- src/solvers/conjugate_gradient_descent.jl | 4 +-- src/solvers/quasi_Newton.jl | 9 +++--- test/plans/test_debug.jl | 6 ++-- 6 files changed, 50 insertions(+), 18 deletions(-) diff --git a/src/plans/conjugate_gradient_plan.jl b/src/plans/conjugate_gradient_plan.jl index 22dfee2596..8d5bce2c6c 100644 --- a/src/plans/conjugate_gradient_plan.jl +++ b/src/plans/conjugate_gradient_plan.jl @@ -5,14 +5,13 @@ struct DirectionUpdateRuleStorage{TC<:DirectionUpdateRule,TStorage<:StoreStateAc end function DirectionUpdateRuleStorage( - M::AbstractManifold, dur::DirectionUpdateRule, p0=rand(M), X0=zero_vector(M, p0) + M::AbstractManifold, dur::DirectionUpdateRule, p=rand(M), X=zero_vector(M, p) ) ursp = update_rule_storage_points(dur) ursv = update_rule_storage_vectors(dur) + # StoreStateAction makes a copy sa = StoreStateAction( - Symbol[], - NamedTuple{ursp}(map(x -> copy(M, p0), ursp)), - NamedTuple{ursv}(map(x -> copy(M, p0, X0), ursv)), + Symbol[], NamedTuple{ursp}(map(x -> p, ursp)), NamedTuple{ursv}(map(x -> X, ursv)) ) return DirectionUpdateRuleStorage{typeof(dur),typeof(sa)}(dur, sa) end @@ -49,6 +48,7 @@ mutable struct ConjugateGradientDescentState{ TVTM<:AbstractVectorTransportMethod, } <: AbstractGradientSolverState p::P + p_old::P X::T δ::T β::F @@ -71,6 +71,7 @@ mutable struct ConjugateGradientDescentState{ βT = allocate_result_type(M, ConjugateGradientDescentState, (p, initial_gradient)) cgs = new{P,T,βT,typeof(coef),typeof(s),typeof(sC),typeof(retr),typeof(vtr)}() cgs.p = p + cgs.p_old = copy(M, p) cgs.X = initial_gradient cgs.δ = copy(M, p, initial_gradient) cgs.stop = sC diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index 50577388f8..08af978954 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -194,9 +194,24 @@ iteration, i.e. acts on `(p,o,i)`, where `p` is a [`AbstractManoptProblem`](@ref # Fields * `values` – a dictionary to store interims values based on certain `Symbols` -* `keys` – an `NTuple` of `Symbols` to refer to fields of `AbstractManoptSolverState` +* `keys` – a `Vector` of `Symbols` to refer to fields of `AbstractManoptSolverState` +* `point_values` – a `NamedTuple` of mutable values of points on a manifold to be stored in + `StoreStateAction`. Manifold is later determined by `AbstractManoptProblem` passed + to `update_storage!`. +* `point_init` – a `NamedTuple` of boolean values indicating whether a point in + `point_values` with matching key has been already initialized to a value. When it is + false, it corresponds to a general value not being stored for the key present in the + vector `keys`. +* `tangent_values` – a `NamedTuple` of mutable values of tangent vectors on a manifold to be + stored in `StoreStateAction`. Manifold is later determined by `AbstractManoptProblem` + passed to `update_storage!`. It is not specified at which point the vectors are tangent + but for storage it should not matter. +* `vector_init` – a `NamedTuple` of boolean values indicating whether a tangent vector in + `tangent_values` with matching key has been already initialized to a value. When it is + false, it corresponds to a general value not being stored for the key present in the + vector `keys`. * `once` – whether to update the internal values only once per iteration -* `lastStored` – last iterate, where this `AbstractStateAction` was called (to determine `once` +* `lastStored` – last iterate, where this `AbstractStateAction` was called (to determine `once`) # Constructiors @@ -205,11 +220,18 @@ iteration, i.e. acts on `(p,o,i)`, where `p` is a [`AbstractManoptProblem`](@ref Initialize the Functor to an (empty) set of keys, where `once` determines whether more that one update per iteration are effective - AbstractStateAction(keys, once=true]) + function StoreStateAction( + general_keys::Vector{Symbol}=Symbol[], + point_values::NamedTuple=NamedTuple(), + tangent_values::NamedTuple=NamedTuple(), + once::Bool=true, + )) Initialize the Functor to a set of keys, where the dictionary is initialized to be empty. Further, `once` determines whether more that one update per iteration are effective, otherwise only the first update is stored, all others are ignored. +Make a copy of points and tangent vectors passed to `point_values` and `tangent_values` +for later storage respective fields. """ mutable struct StoreStateAction{TPS,TXS,TPI,TTI} <: AbstractStateAction values::Dict{Symbol,Any} @@ -230,6 +252,12 @@ mutable struct StoreStateAction{TPS,TXS,TPI,TTI} <: AbstractStateAction tangent_init = NamedTuple{keys(tangent_values)}( map(u -> false, keys(tangent_values)) ) + point_values_copy = NamedTuple{keys(point_values)}( + map(u -> copy(point_values[u]), keys(point_values)) + ) + tangent_values_copy = NamedTuple{keys(tangent_values)}( + map(u -> copy(tangent_values[u]), keys(tangent_values)) + ) return new{ typeof(point_values), typeof(tangent_values), @@ -238,8 +266,8 @@ mutable struct StoreStateAction{TPS,TXS,TPI,TTI} <: AbstractStateAction }( Dict{Symbol,Any}(), general_keys, - point_values, - tangent_values, + point_values_copy, + tangent_values_copy, point_init, tangent_init, once, diff --git a/src/plans/stopping_criterion.jl b/src/plans/stopping_criterion.jl index ddfcaf3f0d..baf4ead078 100644 --- a/src/plans/stopping_criterion.jl +++ b/src/plans/stopping_criterion.jl @@ -121,7 +121,7 @@ For the storage a [`StoreStateAction`](@ref) is used StopWhenChangeLess( ε::Float64; - storage::StoreStateAction=StoreStateAction((:Iterate,)), + storage::StoreStateAction=StoreStateAction([:Iterate]), manifold::AbstractManifold=DefaultManifold(3), inverse_retraction_method::IRT=default_inverse_retraction_method(manifold) ) diff --git a/src/solvers/conjugate_gradient_descent.jl b/src/solvers/conjugate_gradient_descent.jl index 70147e5ea2..8b27ccca9e 100644 --- a/src/solvers/conjugate_gradient_descent.jl +++ b/src/solvers/conjugate_gradient_descent.jl @@ -125,12 +125,12 @@ function initialize_solver!(amp::AbstractManoptProblem, cgs::ConjugateGradientDe end function step_solver!(amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i) M = get_manifold(amp) - p_old = copy(M, cgs.p) + copyto!(M, cgs.p_old, cgs.p) current_stepsize = get_stepsize(amp, cgs, i, cgs.δ) retract!(M, cgs.p, cgs.p, cgs.δ, current_stepsize, cgs.retraction_method) get_gradient!(amp, cgs.X, cgs.p) cgs.β = cgs.coefficient(amp, cgs, i) - vector_transport_to!(M, cgs.δ, p_old, cgs.δ, cgs.p, cgs.vector_transport_method) + vector_transport_to!(M, cgs.δ, cgs.p_old, cgs.δ, cgs.p, cgs.vector_transport_method) cgs.δ .*= cgs.β cgs.δ .-= cgs.X return cgs diff --git a/src/solvers/quasi_Newton.jl b/src/solvers/quasi_Newton.jl index 7b202c42d9..e28dde7ceb 100644 --- a/src/solvers/quasi_Newton.jl +++ b/src/solvers/quasi_Newton.jl @@ -269,18 +269,19 @@ function step_solver!(mp::AbstractManoptProblem, qns::QuasiNewtonState, iter) η = qns.direction_update(mp, qns) α = qns.stepsize(mp, qns, iter, η) p_old = copy(M, get_iterate(qns)) - retract!(M, qns.p, qns.p, α * η, qns.retraction_method) + αη = α * η + retract!(M, qns.p, qns.p, η, α, qns.retraction_method) β = locking_condition_scale( - M, qns.direction_update, p_old, α * η, qns.p, qns.vector_transport_method + M, qns.direction_update, p_old, αη, qns.p, qns.vector_transport_method ) vector_transport_to!( - M, qns.sk, p_old, α * η, qns.p, get_update_vector_transport(qns.direction_update) + M, qns.sk, p_old, αη, qns.p, get_update_vector_transport(qns.direction_update) ) vector_transport_to!( M, qns.X, p_old, qns.X, qns.p, get_update_vector_transport(qns.direction_update) ) new_grad = get_gradient(mp, qns.p) - qns.yk = new_grad / β - qns.X + qns.yk = new_grad ./ β .- qns.X update_hessian!(qns.direction_update, mp, qns, p_old, iter) return qns end diff --git a/test/plans/test_debug.jl b/test/plans/test_debug.jl index 9d2343293a..9a0405a8ad 100644 --- a/test/plans/test_debug.jl +++ b/test/plans/test_debug.jl @@ -48,8 +48,10 @@ end @test String(take!(io)) == "x: $p" DebugEntry(:p; prefix="x:", io=io)(mp, st, -1) @test String(take!(io)) == "" - # Change of Iterate - a2 = DebugChange(; storage=StoreStateAction([:Iterate]), prefix="Last: ", io=io) + # Change of Iterate and recording a custom field + a2 = DebugChange(; + storage=StoreStateAction([:Iterate], (; p=p)), prefix="Last: ", io=io + ) a2(mp, st, 0) # init st.p = [3.0, 2.0] a2(mp, st, 1) From e9b3d17fd42c82611beba621c842db58db3ea64c Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 1 Feb 2023 21:50:20 +0100 Subject: [PATCH 06/37] a bit more docs and type bounds --- src/plans/solver_state.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index 08af978954..5b15ca1141 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -213,6 +213,11 @@ iteration, i.e. acts on `(p,o,i)`, where `p` is a [`AbstractManoptProblem`](@ref * `once` – whether to update the internal values only once per iteration * `lastStored` – last iterate, where this `AbstractStateAction` was called (to determine `once`) +To handle the general storage, use `get_storage` and `has_storage`. For the point storage +use `get_point_storage` and `has_point_storage`. For tangent vector storage use +`get_tangent_storage` and `has_tangent_storage`. Point and tangent storage have been +optimized to be more efficient + # Constructiors AbstractStateAction([keys=(), once=true]) @@ -233,7 +238,9 @@ are effective, otherwise only the first update is stored, all others are ignored Make a copy of points and tangent vectors passed to `point_values` and `tangent_values` for later storage respective fields. """ -mutable struct StoreStateAction{TPS,TXS,TPI,TTI} <: AbstractStateAction +mutable struct StoreStateAction{ + TPS<:NamedTuple,TXS<:NamedTuple,TPI<:NamedTuple,TTI<:NamedTuple +} <: AbstractStateAction values::Dict{Symbol,Any} keys::Vector{Symbol} # for values point_values::TPS @@ -338,6 +345,8 @@ has_tangent_storage(a::AbstractStateAction, key::Symbol) = a.tangent_init[key] Update the [`AbstractStateAction`](@ref) `a` internal values to the ones given on the [`AbstractManoptSolverState`](@ref) `s`. + +Warning: it does not update point and tangent vector storage. """ function update_storage!(a::AbstractStateAction, s::AbstractManoptSolverState) for key in a.keys From 85dc5f8d2ddd913c9d1a7de7f5966070653ab768 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 1 Feb 2023 21:59:52 +0100 Subject: [PATCH 07/37] unify constructors for CG coefficients --- src/plans/conjugate_gradient_plan.jl | 57 +++++++++++++++++----------- 1 file changed, 35 insertions(+), 22 deletions(-) diff --git a/src/plans/conjugate_gradient_plan.jl b/src/plans/conjugate_gradient_plan.jl index 8d5bce2c6c..cf8f9cc93f 100644 --- a/src/plans/conjugate_gradient_plan.jl +++ b/src/plans/conjugate_gradient_plan.jl @@ -63,8 +63,8 @@ mutable struct ConjugateGradientDescentState{ sC::StoppingCriterion, s::Stepsize, dC::DirectionUpdateRule, - retr::AbstractRetractionMethod=ExponentialRetraction(), - vtr::AbstractVectorTransportMethod=ParallelTransport(), + retr::AbstractRetractionMethod=default_retraction_method(M), + vtr::AbstractVectorTransportMethod=default_vector_transport_method(M), initial_gradient::T=zero_vector(M, p), ) where {P,T} coef = DirectionUpdateRuleStorage(M, dC, p, initial_gradient) @@ -163,9 +163,9 @@ Then the coefficient reads See also [`conjugate_gradient_descent`](@ref) # Constructor - DaiYuanCoefficient( - t::AbstractVectorTransportMethod=ParallelTransport(), - a::StoreStateAction=(), + function DaiYuanCoefficient( + M::AbstractManifold=DefaultManifold(2); + t::AbstractVectorTransportMethod=default_vector_transport_method(M) ) Construct the Dai Yuan coefficient update rule, where the parallel transport is the @@ -178,7 +178,10 @@ default vector transport and a new storage is created by default. """ struct DaiYuanCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM - function DaiYuanCoefficient(t::AbstractVectorTransportMethod=ParallelTransport()) + function DaiYuanCoefficient( + M::AbstractManifold=DefaultManifold(2); + t::AbstractVectorTransportMethod=default_vector_transport_method(M), + ) return new{typeof(t)}(t) end end @@ -275,9 +278,9 @@ This method includes a numerical stability proposed by those authors. See also [`conjugate_gradient_descent`](@ref) # Constructor - HagerZhangCoefficient( - t::AbstractVectorTransportMethod=ParallelTransport(), - a::StoreStateAction=(), + function HagerZhangCoefficient( + M::AbstractManifold = DefaultManifold(2); + t::AbstractVectorTransportMethod=default_vector_transport_method(M) ) Construct the Hager Zhang coefficient update rule, where the parallel transport is the @@ -291,7 +294,10 @@ default vector transport and a new storage is created by default. mutable struct HagerZhangCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM - function HagerZhangCoefficient(t::AbstractVectorTransportMethod=ParallelTransport()) + function HagerZhangCoefficient( + M::AbstractManifold=DefaultManifold(2); + t::AbstractVectorTransportMethod=default_vector_transport_method(M), + ) return new{typeof(t)}(t) end end @@ -349,9 +355,9 @@ Then the update reads where ``P_{a\gets b}(⋅)`` denotes a vector transport from the tangent space at ``a`` to ``b``. # Constructor - HestenesStiefelCoefficient( - t::AbstractVectorTransportMethod=ParallelTransport(), - a::StoreStateAction=() + function HestenesStiefelCoefficient( + M::AbstractManifold = DefaultManifold(2); + transport_method::AbstractVectorTransportMethod=default_vector_transport_method(M) ) Construct the Heestens Stiefel coefficient update rule, where the parallel transport is the @@ -368,7 +374,8 @@ struct HestenesStiefelCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM function HestenesStiefelCoefficient( - transport_method::AbstractVectorTransportMethod=ParallelTransport() + M::AbstractManifold=DefaultManifold(2); + transport_method::AbstractVectorTransportMethod=default_vector_transport_method(M), ) return new{typeof(transport_method)}(transport_method) end @@ -422,9 +429,9 @@ Then the coefficient reads See also [`conjugate_gradient_descent`](@ref) # Constructor - LiuStoreyCoefficient( - t::AbstractVectorTransportMethod=ParallelTransport(), - a::StoreStateAction=() + function LiuStoreyCoefficient( + M::AbstractManifold = DefaultManifold(2); + t::AbstractVectorTransportMethod=default_vector_transport_method(M) ) Construct the Lui Storey coefficient update rule, where the parallel transport is the @@ -437,7 +444,10 @@ default vector transport and a new storage is created by default. """ struct LiuStoreyCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM - function LiuStoreyCoefficient(t::AbstractVectorTransportMethod=ParallelTransport()) + function LiuStoreyCoefficient( + M::AbstractManifold=DefaultManifold(2); + t::AbstractVectorTransportMethod=default_vector_transport_method(M), + ) return new{typeof(t)}(t) end end @@ -486,9 +496,9 @@ Then the update reads # Constructor - PolakRibiereCoefficient( - t::AbstractVectorTransportMethod=ParallelTransport(), - a::StoreStateAction=() + function PolakRibiereCoefficient( + M::AbstractManifold=DefaultManifold(2); + t::AbstractVectorTransportMethod=default_vector_transport_method(M) ) Construct the PolakRibiere coefficient update rule, where the parallel transport is the @@ -508,7 +518,10 @@ See also [`conjugate_gradient_descent`](@ref) """ struct PolakRibiereCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM - function PolakRibiereCoefficient(t::AbstractVectorTransportMethod=ParallelTransport()) + function PolakRibiereCoefficient( + M::AbstractManifold=DefaultManifold(2); + t::AbstractVectorTransportMethod=default_vector_transport_method(M), + ) return new{typeof(t)}(t) end end From 5f4059a59cab4b1c8c41713628a786cb35bb04ae Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Fri, 3 Feb 2023 18:52:23 +0100 Subject: [PATCH 08/37] Update one docstring. --- src/solvers/truncated_conjugate_gradient_descent.jl | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/solvers/truncated_conjugate_gradient_descent.jl b/src/solvers/truncated_conjugate_gradient_descent.jl index f3ef396171..bad0078145 100644 --- a/src/solvers/truncated_conjugate_gradient_descent.jl +++ b/src/solvers/truncated_conjugate_gradient_descent.jl @@ -164,17 +164,13 @@ mehtod is larger than the trust-region radius, i.e. $\Vert η_{k}^{*} \Vert_x # Fields * `reason` – stores a reason of stopping if the stopping criterion has one be reached, see [`get_reason`](@ref). -* `storage` – stores the necessary parameters `η, δ, residual` to check the - criterion. # Constructor - StopWhenTrustRegionIsExceeded([a]) + StopWhenTrustRegionIsExceeded() initialize the StopWhenTrustRegionIsExceeded functor to indicate to stop after -the norm of the next iterate is greater than the trust-region radius using the -[`StoreStateAction`](@ref) `a`, which is initialized to store -`:η, :δ, :residual` by default. +the norm of the next iterate is greater than the trust-region radius. # See also From 49f003138d9ef724c373d6eb7857edb134ecf0f9 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sun, 5 Feb 2023 20:58:49 +0100 Subject: [PATCH 09/37] simplify initialization marking --- src/plans/solver_state.jl | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index 5b15ca1141..65610e18c0 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -361,10 +361,6 @@ function update_storage!(a::AbstractStateAction, s::AbstractManoptSolverState) return a.keys end -function _storage_key_true(nt::NamedTuple) - return map(key -> NamedTuple{(key,),Tuple{Bool}}(true), keys(nt)) -end - """ update_storage!(a::AbstractStateAction, amp::AbstractManoptProblem, s::AbstractManoptSolverState) @@ -387,8 +383,7 @@ function update_storage!( M = get_manifold(amp) - pt_kts = _storage_key_true(a.point_values) - map(keys(a.point_values), pt_kts) do key, kt + map(keys(a.point_values)) do key if key === :Iterate copyto!(M, a.point_values[key], get_iterate(s)) else @@ -396,10 +391,10 @@ function update_storage!( M, a.point_values[key], getproperty(s, key)::typeof(a.point_values[key]) ) end - a.point_init = merge(a.point_init, kt) end - tv_kts = _storage_key_true(a.tangent_values) - map(keys(a.tangent_values), tv_kts) do key, kt + a.point_init = NamedTuple{keys(a.point_values)}(map(u -> true, keys(a.point_values))) + + map(keys(a.tangent_values)) do key if key === :Gradient copyto!(M, a.tangent_values[key], get_gradient(s)) else @@ -407,8 +402,11 @@ function update_storage!( M, a.tangent_values[key], getproperty(s, key)::typeof(a.tangent_values[key]) ) end - a.tangent_init = merge(a.tangent_init, kt) end + a.tangent_init = NamedTuple{keys(a.tangent_values)}( + map(u -> true, keys(a.tangent_values)) + ) + return a.keys end From 5983e5aee70a279af3732316863a2af12113973e Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Mon, 6 Feb 2023 18:35:44 +0100 Subject: [PATCH 10/37] unified storage accessors with fallbacks --- src/plans/conjugate_gradient_plan.jl | 76 +++++++++++++++------------- src/plans/debug.jl | 9 +++- src/plans/solver_state.jl | 65 +++++++++++++++++++----- 3 files changed, 100 insertions(+), 50 deletions(-) diff --git a/src/plans/conjugate_gradient_plan.jl b/src/plans/conjugate_gradient_plan.jl index cf8f9cc93f..d245c4cbdf 100644 --- a/src/plans/conjugate_gradient_plan.jl +++ b/src/plans/conjugate_gradient_plan.jl @@ -130,12 +130,13 @@ function (u::DirectionUpdateRuleStorage{ConjugateDescentCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_point_storage(u.storage, :Iterate) || !has_tangent_storage(u.storage, :Gradient) + if !has_storage(u.storage, PointStorageKey(:Iterate)) || + !has_storage(u.storage, TangentStorageKey(:Gradient)) update_storage!(u.storage, amp, cgs) # if not given store current as old return 0.0 end - p_old = get_point_storage(u.storage, :Iterate) - X_old = get_tangent_storage(u.storage, :Gradient) + p_old = get_storage(u.storage, PointStorageKey(:Iterate)) + X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) coef = inner(M, cgs.p, cgs.X, cgs.X) / inner(M, p_old, -cgs.δ, X_old) update_storage!(u.storage, amp, cgs) return coef @@ -193,15 +194,15 @@ function (u::DirectionUpdateRuleStorage{<:DaiYuanCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_point_storage(u.storage, :Iterate) || - !has_tangent_storage(u.storage, :Gradient) || - !has_tangent_storage(u.storage, :δ) + if !has_storage(u.storage, PointStorageKey(:Iterate)) || + !has_storage(u.storage, TangentStorageKey(:Gradient)) || + !has_storage(u.storage, TangentStorageKey(:δ)) update_storage!(u.storage, amp, cgs) # if not given store current as old return 0.0 end - p_old = get_point_storage(u.storage, :Iterate) - X_old = get_tangent_storage(u.storage, :Gradient) - δ_old = get_tangent_storage(u.storage, :δ) + p_old = get_storage(u.storage, PointStorageKey(:Iterate)) + X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) + δ_old = get_storage(u.storage, TangentStorageKey(:δ)) gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) ν = cgs.X - gradienttr #notation y from [HZ06] @@ -245,11 +246,12 @@ function (u::DirectionUpdateRuleStorage{FletcherReevesCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_point_storage(u.storage, :Iterate) || !has_tangent_storage(u.storage, :Gradient) + if !has_storage(u.storage, PointStorageKey(:Iterate)) || + !has_storage(u.storage, TangentStorageKey(:Gradient)) update_storage!(u.storage, amp, cgs) # if not given store current as old end - p_old = get_point_storage(u.storage, :Iterate) - X_old = get_tangent_storage(u.storage, :Gradient) + p_old = get_storage(u.storage, PointStorageKey(:Iterate)) + X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) coef = inner(M, cgs.p, cgs.X, cgs.X) / inner(M, p_old, X_old, X_old) update_storage!(u.storage, amp, cgs) return coef @@ -309,15 +311,15 @@ function (u::DirectionUpdateRuleStorage{<:HagerZhangCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_point_storage(u.storage, :Iterate) || - !has_tangent_storage(u.storage, :Gradient) || - !has_tangent_storage(u.storage, :δ) + if !has_storage(u.storage, PointStorageKey(:Iterate)) || + !has_storage(u.storage, TangentStorageKey(:Gradient)) || + !has_storage(u.storage, TangentStorageKey(:δ)) update_storage!(u.storage, amp, cgs) # if not given store current as old return 0.0 end - p_old = get_point_storage(u.storage, :Iterate) - X_old = get_tangent_storage(u.storage, :Gradient) - δ_old = get_tangent_storage(u.storage, :δ) + p_old = get_storage(u.storage, PointStorageKey(:Iterate)) + X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) + δ_old = get_storage(u.storage, TangentStorageKey(:δ)) gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) ν = cgs.X - gradienttr #notation y from [HZ06] @@ -388,15 +390,15 @@ function (u::DirectionUpdateRuleStorage{<:HestenesStiefelCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_point_storage(u.storage, :Iterate) || - !has_tangent_storage(u.storage, :Gradient) || - !has_tangent_storage(u.storage, :δ) + if !has_storage(u.storage, PointStorageKey(:Iterate)) || + !has_storage(u.storage, TangentStorageKey(:Gradient)) || + !has_storage(u.storage, TangentStorageKey(:δ)) update_storage!(u.storage, amp, cgs) # if not given store current as old return 0.0 end - p_old = get_point_storage(u.storage, :Iterate) - X_old = get_tangent_storage(u.storage, :Gradient) - δ_old = get_tangent_storage(u.storage, :δ) + p_old = get_storage(u.storage, PointStorageKey(:Iterate)) + X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) + δ_old = get_storage(u.storage, TangentStorageKey(:δ)) gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) δtr = vector_transport_to(M, p_old, δ_old, cgs.p, u.coefficient.transport_method) @@ -459,14 +461,14 @@ function (u::DirectionUpdateRuleStorage{<:LiuStoreyCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_point_storage(u.storage, :Iterate) || - !has_tangent_storage(u.storage, :Gradient) || - !has_tangent_storage(u.storage, :δ) + if !has_storage(u.storage, PointStorageKey(:Iterate)) || + !has_storage(u.storage, TangentStorageKey(:Gradient)) || + !has_storage(u.storage, TangentStorageKey(:δ)) update_storage!(u.storage, amp, cgs) # if not given store current as old end - p_old = get_point_storage(u.storage, :Iterate) - X_old = get_tangent_storage(u.storage, :Gradient) - δ_old = get_tangent_storage(u.storage, :δ) + p_old = get_storage(u.storage, PointStorageKey(:Iterate)) + X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) + δ_old = get_storage(u.storage, TangentStorageKey(:δ)) gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) ν = cgs.X - gradienttr # notation y from [HZ06] coef = inner(M, cgs.p, cgs.X, ν) / inner(M, p_old, -δ_old, X_old) @@ -533,11 +535,12 @@ function (u::DirectionUpdateRuleStorage{<:PolakRibiereCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_point_storage(u.storage, :Iterate) || !has_tangent_storage(u.storage, :Gradient) + if !has_storage(u.storage, PointStorageKey(:Iterate)) || + !has_storage(u.storage, TangentStorageKey(:Gradient)) update_storage!(u.storage, amp, cgs) # if not given store current as old end - p_old = get_point_storage(u.storage, :Iterate) - X_old = get_tangent_storage(u.storage, :Gradient) + p_old = get_storage(u.storage, PointStorageKey(:Iterate)) + X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) ν = cgs.X - gradienttr @@ -641,11 +644,12 @@ function (u::DirectionUpdateRuleStorage{<:ConjugateGradientBealeRestart})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i ) M = get_manifold(amp) - if !has_point_storage(u.storage, :Iterate) || !has_tangent_storage(u.storage, :Gradient) + if !has_storage(u.storage, PointStorageKey(:Iterate)) || + !has_storage(u.storage, TangentStorageKey(:Gradient)) update_storage!(u.storage, amp, cgs) # if not given store current as old end - p_old = get_point_storage(u.storage, :Iterate) - X_old = get_tangent_storage(u.storage, :Gradient) + p_old = get_storage(u.storage, PointStorageKey(:Iterate)) + X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) # call actual rule β = u.coefficient.direction_update(amp, cgs, i) diff --git a/src/plans/debug.jl b/src/plans/debug.jl index de1ace8f49..a6ac3b997d 100644 --- a/src/plans/debug.jl +++ b/src/plans/debug.jl @@ -171,7 +171,12 @@ function (d::DebugChange)(mp::AbstractManoptProblem, st::AbstractManoptSolverSta (i > 0) && Printf.format( d.io, Printf.Format(d.format), - distance(M, get_iterate(st), get_storage(d.storage, :Iterate), d.invretr), + distance( + M, + get_iterate(st), + get_storage(d.storage, PointStorageKey(:Iterate)), + d.invretr, + ), ) d.storage(mp, st, i) return nothing @@ -208,7 +213,7 @@ function (d::DebugGradientChange)( (i > 0) && Printf.format( d.io, Printf.Format(d.format), - distance(M, get_gradient(st), get_storage(d.storage, :Gradient)), + distance(M, get_gradient(st), get_storage(d.storage, TangentStorageKey(:Gradient))), ) d.storage(pm, st, i) return nothing diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index 65610e18c0..d5bd77319f 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -171,6 +171,24 @@ end get_solver_result(s::AbstractManoptSolverState, ::Val{false}) = get_iterate(s) get_solver_result(s::AbstractManoptSolverState, ::Val{true}) = get_solver_result(s.state) +""" + struct PointStorageKey{key} end + +Refer to point storage of [`StoreStateAction`](@ref) in `get_storage` and `has_storage` +functions +""" +struct PointStorageKey{key} end +PointStorageKey(key::Symbol) = PointStorageKey{key}() + +""" + struct TangentStorageKey{key} end + +Refer to tangent storage of [`StoreStateAction`](@ref) in `get_storage` and `has_storage` +functions +""" +struct TangentStorageKey{key} end +TangentStorageKey(key::Symbol) = TangentStorageKey{key}() + # # Common Actions for decorated AbstractManoptSolverState # @@ -213,10 +231,9 @@ iteration, i.e. acts on `(p,o,i)`, where `p` is a [`AbstractManoptProblem`](@ref * `once` – whether to update the internal values only once per iteration * `lastStored` – last iterate, where this `AbstractStateAction` was called (to determine `once`) -To handle the general storage, use `get_storage` and `has_storage`. For the point storage -use `get_point_storage` and `has_point_storage`. For tangent vector storage use -`get_tangent_storage` and `has_tangent_storage`. Point and tangent storage have been -optimized to be more efficient +To handle the general storage, use `get_storage` and `has_storage` with keys as `Symbol`s. +For the point storage use `PointStorageKey`. For tangent vector storage use +`TangentStorageKey`. Point and tangent storage have been optimized to be more efficient. # Constructiors @@ -301,20 +318,32 @@ Return the internal value of the [`AbstractStateAction`](@ref) `a` at the get_storage(a::AbstractStateAction, key::Symbol) = a.values[key] """ - get_point_storage(a::AbstractStateAction, key::Symbol) + get_storage(a::AbstractStateAction, ::PointStorageKey{key}) where {key} Return the internal value of the [`AbstractStateAction`](@ref) `a` at the `Symbol` `key` that represents a point. """ -get_point_storage(a::AbstractStateAction, key::Symbol) = a.point_values[key] +function get_storage(a::AbstractStateAction, ::PointStorageKey{key}) where {key} + if haskey(a.point_values, key) + return a.point_values[key] + else + return get_storage(a, key) + end +end """ - get_tangent_storage(a::AbstractStateAction, key::Symbol) + get_storage(a::AbstractStateAction, ::TangentStorageKey{key}) where {key} Return the internal value of the [`AbstractStateAction`](@ref) `a` at the `Symbol` `key` that represents a tangent vector. """ -get_tangent_storage(a::AbstractStateAction, key::Symbol) = a.tangent_values[key] +function get_storage(a::AbstractStateAction, ::TangentStorageKey{key}) where {key} + if haskey(a.tangent_values, key) + return a.tangent_values[key] + else + return get_storage(a, key) + end +end """ has_storage(a::AbstractStateAction, key::Symbol) @@ -325,20 +354,32 @@ Return whether the [`AbstractStateAction`](@ref) `a` has a value stored at the has_storage(a::AbstractStateAction, key::Symbol) = haskey(a.values, key) """ - has_point_storage(a::AbstractStateAction, key::Symbol) + has_storage(a::AbstractStateAction, ::PointStorageKey{key}) where {key} Return whether the [`AbstractStateAction`](@ref) `a` has a point value stored at the `Symbol` `key`. """ -has_point_storage(a::AbstractStateAction, key::Symbol) = a.point_init[key] +function has_storage(a::AbstractStateAction, ::PointStorageKey{key}) where {key} + if haskey(a.point_init, key) + return a.point_init[key] + else + return has_storage(a, key) + end +end """ - has_tangent_storage(a::AbstractStateAction, key::Symbol) + has_storage(a::AbstractStateAction, ::TangentStorageKey{key}) where {key} Return whether the [`AbstractStateAction`](@ref) `a` has a point value stored at the `Symbol` `key`. """ -has_tangent_storage(a::AbstractStateAction, key::Symbol) = a.tangent_init[key] +function has_storage(a::AbstractStateAction, ::TangentStorageKey{key}) where {key} + if haskey(a.tangent_init, key) + return a.tangent_init[key] + else + return has_storage(a, key) + end +end """ update_storage!(a::AbstractStateAction, s::AbstractManoptSolverState) From d95a1c9b8bbed95c41c264587fc4c259165065c4 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 8 Feb 2023 12:29:32 +0100 Subject: [PATCH 11/37] new StoreStateAction constructor --- src/plans/conjugate_gradient_plan.jl | 44 +++++++++--------- src/plans/solver_state.jl | 52 +++++++++++++--------- src/plans/stepsize.jl | 2 +- test/plans/test_conjugate_gradient_plan.jl | 4 +- test/plans/test_debug.jl | 4 +- test/plans/test_record.jl | 2 +- 6 files changed, 59 insertions(+), 49 deletions(-) diff --git a/src/plans/conjugate_gradient_plan.jl b/src/plans/conjugate_gradient_plan.jl index d245c4cbdf..4438b17f15 100644 --- a/src/plans/conjugate_gradient_plan.jl +++ b/src/plans/conjugate_gradient_plan.jl @@ -10,9 +10,7 @@ function DirectionUpdateRuleStorage( ursp = update_rule_storage_points(dur) ursv = update_rule_storage_vectors(dur) # StoreStateAction makes a copy - sa = StoreStateAction( - Symbol[], NamedTuple{ursp}(map(x -> p, ursp)), NamedTuple{ursv}(map(x -> X, ursv)) - ) + sa = StoreStateAction(M, Symbol[], ursp, ursv; p_init=p, X_init=X) return DirectionUpdateRuleStorage{typeof(dur),typeof(sa)}(dur, sa) end @@ -123,8 +121,8 @@ Construct the conjugate descent coefficient update rule, a new storage is create """ struct ConjugateDescentCoefficient <: DirectionUpdateRule end -update_rule_storage_points(::ConjugateDescentCoefficient) = (:Iterate,) -update_rule_storage_vectors(::ConjugateDescentCoefficient) = (:Gradient,) +update_rule_storage_points(::ConjugateDescentCoefficient) = Tuple{:Iterate} +update_rule_storage_vectors(::ConjugateDescentCoefficient) = Tuple{:Gradient} function (u::DirectionUpdateRuleStorage{ConjugateDescentCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i @@ -187,8 +185,8 @@ struct DaiYuanCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdat end end -update_rule_storage_points(::DaiYuanCoefficient) = (:Iterate,) -update_rule_storage_vectors(::DaiYuanCoefficient) = (:Gradient, :δ) +update_rule_storage_points(::DaiYuanCoefficient) = Tuple{:Iterate} +update_rule_storage_vectors(::DaiYuanCoefficient) = Tuple{:Gradient,:δ} function (u::DirectionUpdateRuleStorage{<:DaiYuanCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i @@ -239,8 +237,8 @@ Construct the Fletcher Reeves coefficient update rule, a new storage is created """ struct FletcherReevesCoefficient <: DirectionUpdateRule end -update_rule_storage_points(::FletcherReevesCoefficient) = (:Iterate,) -update_rule_storage_vectors(::FletcherReevesCoefficient) = (:Gradient,) +update_rule_storage_points(::FletcherReevesCoefficient) = Tuple{:Iterate} +update_rule_storage_vectors(::FletcherReevesCoefficient) = Tuple{:Gradient} function (u::DirectionUpdateRuleStorage{FletcherReevesCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i @@ -304,8 +302,8 @@ mutable struct HagerZhangCoefficient{TVTM<:AbstractVectorTransportMethod} <: end end -update_rule_storage_points(::HagerZhangCoefficient) = (:Iterate,) -update_rule_storage_vectors(::HagerZhangCoefficient) = (:Gradient, :δ) +update_rule_storage_points(::HagerZhangCoefficient) = Tuple{:Iterate} +update_rule_storage_vectors(::HagerZhangCoefficient) = Tuple{:Gradient,:δ} function (u::DirectionUpdateRuleStorage{<:HagerZhangCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i @@ -383,8 +381,8 @@ struct HestenesStiefelCoefficient{TVTM<:AbstractVectorTransportMethod} <: end end -update_rule_storage_points(::HestenesStiefelCoefficient) = (:Iterate,) -update_rule_storage_vectors(::HestenesStiefelCoefficient) = (:Gradient, :δ) +update_rule_storage_points(::HestenesStiefelCoefficient) = Tuple{:Iterate} +update_rule_storage_vectors(::HestenesStiefelCoefficient) = Tuple{:Gradient,:δ} function (u::DirectionUpdateRuleStorage{<:HestenesStiefelCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i @@ -454,8 +452,8 @@ struct LiuStoreyCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpd end end -update_rule_storage_points(::LiuStoreyCoefficient) = (:Iterate,) -update_rule_storage_vectors(::LiuStoreyCoefficient) = (:Gradient, :δ) +update_rule_storage_points(::LiuStoreyCoefficient) = Tuple{:Iterate} +update_rule_storage_vectors(::LiuStoreyCoefficient) = Tuple{:Gradient,:δ} function (u::DirectionUpdateRuleStorage{<:LiuStoreyCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i @@ -528,8 +526,8 @@ struct PolakRibiereCoefficient{TVTM<:AbstractVectorTransportMethod} <: Direction end end -update_rule_storage_points(::PolakRibiereCoefficient) = (:Iterate,) -update_rule_storage_vectors(::PolakRibiereCoefficient) = (:Gradient,) +update_rule_storage_points(::PolakRibiereCoefficient) = Tuple{:Iterate} +update_rule_storage_vectors(::PolakRibiereCoefficient) = Tuple{:Gradient} function (u::DirectionUpdateRuleStorage{<:PolakRibiereCoefficient})( amp::AbstractManoptProblem, cgs::ConjugateGradientDescentState, i @@ -559,8 +557,8 @@ See also [`conjugate_gradient_descent`](@ref) """ struct SteepestDirectionUpdateRule <: DirectionUpdateRule end -update_rule_storage_points(::SteepestDirectionUpdateRule) = () -update_rule_storage_vectors(::SteepestDirectionUpdateRule) = () +update_rule_storage_points(::SteepestDirectionUpdateRule) = Tuple{} +update_rule_storage_vectors(::SteepestDirectionUpdateRule) = Tuple{} function (u::DirectionUpdateRuleStorage{SteepestDirectionUpdateRule})( ::DefaultManoptProblem, ::ConjugateGradientDescentState, i @@ -631,13 +629,13 @@ mutable struct ConjugateGradientBealeRestart{ end end -function update_rule_storage_points(dur::ConjugateGradientBealeRestart) +@inline function update_rule_storage_points(dur::ConjugateGradientBealeRestart) dur_p = update_rule_storage_points(dur.direction_update) - return :Iterate in dur_p ? dur_p : (:Iterate, dur_p...) + return :Iterate in dur_p.parameters ? dur_p : Tuple{:Iterate,dur_p.parameters...} end -function update_rule_storage_vectors(dur::ConjugateGradientBealeRestart) +@inline function update_rule_storage_vectors(dur::ConjugateGradientBealeRestart) dur_X = update_rule_storage_vectors(dur.direction_update) - return :Gradient in dur_X ? dur_X : (:Gradient, dur_X...) + return :Gradient in dur_X.parameters ? dur_X : Tuple{:Gradient,dur_X.parameters...} end function (u::DirectionUpdateRuleStorage{<:ConjugateGradientBealeRestart})( diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index d5bd77319f..a67d45a64a 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -242,6 +242,21 @@ For the point storage use `PointStorageKey`. For tangent vector storage use Initialize the Functor to an (empty) set of keys, where `once` determines whether more that one update per iteration are effective + function StoreStateAction( + M::AbstractManifold, + dictionary_symbols::Vector{Symbol}, + ::Type{TPS}, + ::Type{TTS}; + p_init=rand(M), + X_init=zero_vector(M, p_init), + once=true, + ) where {TPS<:Tuple,TTS<:Tuple} + +Initialize the general storage keys to `dictionary_symbols`, point storage keys to `TPS` and +tangent vector storage tu `TTS`. For example you may call +`StorageStateAction(M, Symbol[], Tuple{:Iterate}, Tuple{:Gradient})` to create efficient +storage for point representing iterate and tangent vector representing gradient. + function StoreStateAction( general_keys::Vector{Symbol}=Symbol[], point_values::NamedTuple=NamedTuple(), @@ -299,6 +314,22 @@ mutable struct StoreStateAction{ ) end end +@inline function StoreStateAction( + M::AbstractManifold, + dictionary_symbols::Vector{Symbol}, + ::Type{TPS}, + ::Type{TTS}; + p_init=rand(M), + X_init=zero_vector(M, p_init), + once=true, +) where {TPS<:Tuple,TTS<:Tuple} + TPS_tuple = Tuple(TPS.parameters) + TTS_tuple = Tuple(TTS.parameters) + point_values = NamedTuple{TPS_tuple}(map(_ -> p_init, TPS_tuple)) + tangent_values = NamedTuple{TTS_tuple}(map(_ -> X_init, TTS_tuple)) + return StoreStateAction(dictionary_symbols, point_values, tangent_values, once) +end + function (a::StoreStateAction)( amp::AbstractManoptProblem, s::AbstractManoptSolverState, i::Int ) @@ -381,27 +412,6 @@ function has_storage(a::AbstractStateAction, ::TangentStorageKey{key}) where {ke end end -""" - update_storage!(a::AbstractStateAction, s::AbstractManoptSolverState) - -Update the [`AbstractStateAction`](@ref) `a` internal values to the ones given on -the [`AbstractManoptSolverState`](@ref) `s`. - -Warning: it does not update point and tangent vector storage. -""" -function update_storage!(a::AbstractStateAction, s::AbstractManoptSolverState) - for key in a.keys - if key === :Iterate - a.values[key] = deepcopy(get_iterate(s)) - elseif key === :Gradient - a.values[key] = deepcopy(get_gradient(s)) - else - a.values[key] = deepcopy(getproperty(s, key)) - end - end - return a.keys -end - """ update_storage!(a::AbstractStateAction, amp::AbstractManoptProblem, s::AbstractManoptSolverState) diff --git a/src/plans/stepsize.jl b/src/plans/stepsize.jl index 334ea394a9..07f5e07dea 100644 --- a/src/plans/stepsize.jl +++ b/src/plans/stepsize.jl @@ -473,7 +473,7 @@ function (a::NonmonotoneLinesearch)( p_old = get_storage(a.storage, :Iterate) X_old = get_storage(a.storage, :Gradient) end - update_storage!(a.storage, s) + update_storage!(a.storage, mp, s) return a( get_manifold(mp), get_iterate(s), diff --git a/test/plans/test_conjugate_gradient_plan.jl b/test/plans/test_conjugate_gradient_plan.jl index 058c3507ee..1129d6add2 100644 --- a/test/plans/test_conjugate_gradient_plan.jl +++ b/test/plans/test_conjugate_gradient_plan.jl @@ -2,8 +2,8 @@ using Manopt, Manifolds, Test struct DummyCGCoeff <: DirectionUpdateRule end (u::DummyCGCoeff)(p, s, i) = 0.2 -Manopt.update_rule_storage_points(::DummyCGCoeff) = () -Manopt.update_rule_storage_vectors(::DummyCGCoeff) = () +Manopt.update_rule_storage_points(::DummyCGCoeff) = Tuple{} +Manopt.update_rule_storage_vectors(::DummyCGCoeff) = Tuple{} @testset "Test Restart CG" begin M = Euclidean(2) diff --git a/test/plans/test_debug.jl b/test/plans/test_debug.jl index 9a0405a8ad..5198b3ac9d 100644 --- a/test/plans/test_debug.jl +++ b/test/plans/test_debug.jl @@ -50,7 +50,9 @@ end @test String(take!(io)) == "" # Change of Iterate and recording a custom field a2 = DebugChange(; - storage=StoreStateAction([:Iterate], (; p=p)), prefix="Last: ", io=io + storage=StoreStateAction(M, Symbol[], Tuple{:Iterate}, Tuple{}; p_init=p), + prefix="Last: ", + io=io, ) a2(mp, st, 0) # init st.p = [3.0, 2.0] diff --git a/test/plans/test_record.jl b/test/plans/test_record.jl index 77b255b2ef..925e2ca844 100644 --- a/test/plans/test_record.jl +++ b/test/plans/test_record.jl @@ -127,7 +127,7 @@ using Manifolds, Manopt, Test, ManifoldsBase, Dates # RecordEntryChange set_iterate!(gds, M, p) e = RecordEntryChange(:p, (p, o, x, y) -> distance(get_manifold(p), x, y)) - @test update_storage!(e.storage, gds) == [:p] + @test update_storage!(e.storage, dmp, gds) == [:p] e2 = RecordEntryChange(dmp, :p, (p, o, x, y) -> distance(get_manifold(p), x, y)) @test e.field == e2.field e(dmp, gds, 1) From 2a6de87af698d8670e41e9f7c3a7d6a30d80a20a Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Thu, 9 Feb 2023 18:43:13 +0100 Subject: [PATCH 12/37] storage tests --- test/plans/test_storage.jl | 22 ++++++++++++++++++++++ test/runtests.jl | 1 + 2 files changed, 23 insertions(+) create mode 100644 test/plans/test_storage.jl diff --git a/test/plans/test_storage.jl b/test/plans/test_storage.jl new file mode 100644 index 0000000000..c198d0200d --- /dev/null +++ b/test/plans/test_storage.jl @@ -0,0 +1,22 @@ +using Test, Manopt, ManifoldsBase + +@testset "StoreStateAction" begin + M = ManifoldsBase.DefaultManifold(2) + p = [4.0, 2.0] + st = GradientDescentState( + M, p; stopping_criterion=StopAfterIteration(20), stepsize=ConstantStepsize(M) + ) + f(M, q) = distance(M, q, p) .^ 2 + grad_f(M, q) = -2 * log(M, q, p) + mp = DefaultManoptProblem(M, ManifoldGradientObjective(f, grad_f)) + + a = StoreStateAction(M, [:p, :X], Tuple{}, Tuple{}) + + @test !has_storage(a, Manopt.PointStorageKey(:p)) + @test !has_storage(a, Manopt.TangentStorageKey(:X)) + update_storage!(a, mp, st) + @test has_storage(a, Manopt.PointStorageKey(:p)) + @test has_storage(a, Manopt.TangentStorageKey(:X)) + @test get_storage(a, Manopt.PointStorageKey(:p)) == p + @test get_storage(a, Manopt.PointStorageKey(:X)) == [0.0, 0.0] +end diff --git a/test/runtests.jl b/test/runtests.jl index 76ed279ead..e72a1168e8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,7 @@ using Manifolds, Manopt, ManifoldsBase, Test include("plans/test_state.jl") include("plans/test_conjugate_gradient_plan.jl") include("plans/test_debug.jl") + include("plans/test_storage.jl") include("plans/test_cache_plan.jl") include("plans/test_nelder_mead_plan.jl") include("plans/test_nonmutating.jl") From a13f9202047512772ea3b7204bfacb43eaa2b1b3 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Thu, 9 Feb 2023 21:03:33 +0100 Subject: [PATCH 13/37] even more tests --- test/plans/test_storage.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/test/plans/test_storage.jl b/test/plans/test_storage.jl index c198d0200d..288d96cb0a 100644 --- a/test/plans/test_storage.jl +++ b/test/plans/test_storage.jl @@ -18,5 +18,14 @@ using Test, Manopt, ManifoldsBase @test has_storage(a, Manopt.PointStorageKey(:p)) @test has_storage(a, Manopt.TangentStorageKey(:X)) @test get_storage(a, Manopt.PointStorageKey(:p)) == p - @test get_storage(a, Manopt.PointStorageKey(:X)) == [0.0, 0.0] + @test get_storage(a, Manopt.TangentStorageKey(:X)) == [0.0, 0.0] + + a2 = StoreStateAction(M, Symbol[], Tuple{:p}, Tuple{:X}) + @test !has_storage(a2, Manopt.PointStorageKey(:p)) + @test !has_storage(a2, Manopt.TangentStorageKey(:X)) + update_storage!(a2, mp, st) + @test has_storage(a2, Manopt.PointStorageKey(:p)) + @test has_storage(a2, Manopt.TangentStorageKey(:X)) + @test get_storage(a2, Manopt.PointStorageKey(:p)) == p + @test get_storage(a2, Manopt.TangentStorageKey(:X)) == [0.0, 0.0] end From 18136a003c54e122529854d298165572ff9ad251 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Wed, 15 Feb 2023 17:21:09 +0100 Subject: [PATCH 14/37] adapt to scaled gradients and some minor changes --- src/plans/stepsize.jl | 4 ++-- test/functions/test_bezier.jl | 8 ++++---- test/solvers/test_alternating_gradient.jl | 6 ++++-- 3 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/plans/stepsize.jl b/src/plans/stepsize.jl index 07f5e07dea..3206092673 100644 --- a/src/plans/stepsize.jl +++ b/src/plans/stepsize.jl @@ -248,7 +248,7 @@ perform a linesearch for * an initial stepsize `s` usually called ``γ`` * a sufficient `decrease` * a `contract`ion factor ``σ`` -* a `retr`action, which defaults to the `ExponentialRetraction()` +* a `retr`action, which defaults to the `default_retraction_method(M)` * a search direction ``η = -\operatorname{grad}F(x)`` * an offset, ``f_0 = F(x)`` * a keyword `stop_step` as a minimal step size when to stop @@ -261,7 +261,7 @@ function linesearch_backtrack( s, decrease, contract, - retr::AbstractRetractionMethod=ExponentialRetraction(), + retr::AbstractRetractionMethod=default_retraction_method(M), η::T=-gradFx, f0=F(x); stop_step=0.0, diff --git a/test/functions/test_bezier.jl b/test/functions/test_bezier.jl index 0840a71fff..5719515889 100644 --- a/test/functions/test_bezier.jl +++ b/test/functions/test_bezier.jl @@ -48,7 +48,7 @@ using Manopt, Manifolds, Test z = zero_vector(Mp, Bvec) distance(Mp, grad_acceleration_bezier(M, Bvec, degrees, T), z) @test norm(Mp, Bvec, grad_acceleration_bezier(M, Bvec, degrees, T) - z) ≈ 0 atol = - 10^(-12) + 2e-12 d = [pT, exp(M, pC, [0.3, 0.0, 0.0]), pB] λ = 3.0 @@ -61,15 +61,15 @@ using Manopt, Manifolds, Test # when the data are the junctions @test norm( Mp, Bvec, grad_L2_acceleration_bezier(M, Bvec, degrees, T, λ, [pT, pC, pB]) - z - ) ≈ 0 atol = 10^(-12) + ) ≈ 0 atol = 2e-12 z[4][1] = -0.9 @test norm(Mp, Bvec, grad_L2_acceleration_bezier(M, Bvec, degrees, T, λ, d) - z) ≈ 0 atol = - 10^(-12) + 2e-12 # when the data is weighted with zero @test cost_L2_acceleration_bezier(M, Bvec, degrees, T, 0.0, d) ≈ 0 atol = 10^(-10) z[4][1] = 0.0 @test norm(Mp, Bvec, grad_L2_acceleration_bezier(M, Bvec, degrees, T, 0.0, d) - z) ≈ - 0 atol = 10^(-12) + 0 atol = 2e-12 end @testset "de Casteljau variants" begin M = Sphere(2) diff --git a/test/solvers/test_alternating_gradient.jl b/test/solvers/test_alternating_gradient.jl index 430b28347c..55db588dee 100644 --- a/test/solvers/test_alternating_gradient.jl +++ b/test/solvers/test_alternating_gradient.jl @@ -5,8 +5,10 @@ using Manopt, Manifolds, Test N = M × M data = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]] function f(N, p) - return 1 / 2 * - (distance(N[1], p[N, 1], data[1])^2 + distance(N[2], p[N, 2], data[2])^2) + return 1 / 2 * ( + distance(N[1], p[N, Val(1)], data[1])^2 + + distance(N[2], p[N, Val(2)], data[2])^2 + ) end grad_f1(N, p) = -log(N[1], p[N, 1], data[1]) grad_f1!(N, X, p) = (X .= -log(N[1], p[N, 1], data[1])) From 5d5bc8199d495effb5f54b8fe51e27533631c06f Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Thu, 16 Feb 2023 10:56:03 +0100 Subject: [PATCH 15/37] Update DebugChange to a unified field name scheme and tries to use the new storage. --- docs/Project.toml | 2 +- src/plans/debug.jl | 36 ++++++++++++++++++++++++++++-------- test/plans/test_debug.jl | 12 ++++++------ 3 files changed, 35 insertions(+), 15 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 684e95481a..91bc5162d6 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -21,4 +21,4 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" BenchmarkTools = "1.3" Documenter = "0.27" Manifolds = "0.8.27" -ManifoldsBase = "0.13" +ManifoldsBase = "0.13, 0.14" diff --git a/src/plans/debug.jl b/src/plans/debug.jl index a6ac3b997d..52ce9b8613 100644 --- a/src/plans/debug.jl +++ b/src/plans/debug.jl @@ -152,18 +152,38 @@ mutable struct DebugChange{TInvRetr<:AbstractInverseRetractionMethod} <: DebugAc io::IO format::String storage::StoreStateAction - invretr::TInvRetr - function DebugChange(; - storage::StoreStateAction=StoreStateAction([:Iterate]), + inverse_retraction_method::TInvRetr + function DebugChange( + M::AbstractManifold=DefaultManifold(); + storage::Union{Nothing,StoreStateAction}=nothing, io::IO=stdout, prefix::String="Last Change: ", format::String="$(prefix)%f", - manifold::AbstractManifold=DefaultManifold(1), - invretr::AbstractInverseRetractionMethod=default_inverse_retraction_method( - manifold + manifold::Union{Nothing,AbstractManifold}=nothing, + invretr::Union{Nothing,AbstractInverseRetractionMethod}=nothing, + inverse_retraction_method::AbstractInverseRetractionMethod=default_inverse_retraction_method( + M ), ) - return new{typeof(invretr)}(io, format, storage, invretr) + irm = inverse_retraction_method + # Deprecated, remove in Manopt 0.5 + if !isnothing(manifold) + @warn "The `manifold` keyword is deprecated, use the first positional argument `M`. This keyword for now sets `inverse_retracion_method`." + irm = default_inverse_retraction_method(manifold) + end + if !isnothing(invretr) + @warn "invretr keyword is deprecated, use `inverse_retraction_method`, which this one overrides for now." + irm = invretr + end + # Is this how this is intended? + if isnothing(storage) + if M isa DefaultManifold + storage = StoreStateAction([:Iterate]) + else + storage = StoreStateAction(M, Symbol[], Tuple{:Iterate}) + end + end + return new{typeof(irm)}(io, format, storage, irm) end end function (d::DebugChange)(mp::AbstractManoptProblem, st::AbstractManoptSolverState, i) @@ -175,7 +195,7 @@ function (d::DebugChange)(mp::AbstractManoptProblem, st::AbstractManoptSolverSta M, get_iterate(st), get_storage(d.storage, PointStorageKey(:Iterate)), - d.invretr, + d.inverse_retraction_method, ), ) d.storage(mp, st, i) diff --git a/test/plans/test_debug.jl b/test/plans/test_debug.jl index 5198b3ac9d..d0c416f4f8 100644 --- a/test/plans/test_debug.jl +++ b/test/plans/test_debug.jl @@ -61,17 +61,17 @@ end storage=StoreStateAction([:Iterate]), prefix="Last: ", io=io, - invretr=PolarInverseRetraction(), + inverse_retraction_method=PolarInverseRetraction(), ) - a2mani = DebugChange(; + a2mani = DebugChange( + TestPolarManifold(); storage=StoreStateAction([:Iterate]), prefix="Last: ", io=io, - manifold=TestPolarManifold(), ) - @test a2inv.invretr === PolarInverseRetraction() - @test a2mani.invretr === PolarInverseRetraction() - @test a2.invretr === LogarithmicInverseRetraction() + @test a2inv.inverse_retraction_method === PolarInverseRetraction() + @test a2mani.inverse_retraction_method === PolarInverseRetraction() + @test a2.inverse_retraction_method === LogarithmicInverseRetraction() @test String(take!(io)) == "Last: 1.000000" # Change of Gradient a3 = DebugGradientChange(; From 71ce59627bf82592666a013d21b63904105bdff3 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Thu, 16 Feb 2023 12:02:05 +0100 Subject: [PATCH 16/37] addressing review --- src/plans/conjugate_gradient_plan.jl | 68 ++++++++++++++-------------- 1 file changed, 33 insertions(+), 35 deletions(-) diff --git a/src/plans/conjugate_gradient_plan.jl b/src/plans/conjugate_gradient_plan.jl index 4438b17f15..58177f1116 100644 --- a/src/plans/conjugate_gradient_plan.jl +++ b/src/plans/conjugate_gradient_plan.jl @@ -5,12 +5,15 @@ struct DirectionUpdateRuleStorage{TC<:DirectionUpdateRule,TStorage<:StoreStateAc end function DirectionUpdateRuleStorage( - M::AbstractManifold, dur::DirectionUpdateRule, p=rand(M), X=zero_vector(M, p) + M::AbstractManifold, + dur::DirectionUpdateRule, + p_init=rand(M), + X_init=zero_vector(M, p_init), ) ursp = update_rule_storage_points(dur) ursv = update_rule_storage_vectors(dur) # StoreStateAction makes a copy - sa = StoreStateAction(M, Symbol[], ursp, ursv; p_init=p, X_init=X) + sa = StoreStateAction(M, Symbol[], ursp, ursv; p_init=p_init, X_init=X_init) return DirectionUpdateRuleStorage{typeof(dur),typeof(sa)}(dur, sa) end @@ -177,13 +180,13 @@ default vector transport and a new storage is created by default. """ struct DaiYuanCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM - function DaiYuanCoefficient( - M::AbstractManifold=DefaultManifold(2); - t::AbstractVectorTransportMethod=default_vector_transport_method(M), - ) + function DaiYuanCoefficient(t::AbstractVectorTransportMethod) return new{typeof(t)}(t) end end +function DaiYuanCoefficient(M::AbstractManifold=DefaultManifold(2)) + return DaiYuanCoefficient(default_vector_transport_method(M)) +end update_rule_storage_points(::DaiYuanCoefficient) = Tuple{:Iterate} update_rule_storage_vectors(::DaiYuanCoefficient) = Tuple{:Gradient,:δ} @@ -278,10 +281,8 @@ This method includes a numerical stability proposed by those authors. See also [`conjugate_gradient_descent`](@ref) # Constructor - function HagerZhangCoefficient( - M::AbstractManifold = DefaultManifold(2); - t::AbstractVectorTransportMethod=default_vector_transport_method(M) - ) + function HagerZhangCoefficient(t::AbstractVectorTransportMethod) + function HagerZhangCoefficient(M::AbstractManifold = DefaultManifold(2)) Construct the Hager Zhang coefficient update rule, where the parallel transport is the default vector transport and a new storage is created by default. @@ -294,13 +295,14 @@ default vector transport and a new storage is created by default. mutable struct HagerZhangCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM - function HagerZhangCoefficient( - M::AbstractManifold=DefaultManifold(2); - t::AbstractVectorTransportMethod=default_vector_transport_method(M), - ) + + function HagerZhangCoefficient(t::AbstractVectorTransportMethod) return new{typeof(t)}(t) end end +function HagerZhangCoefficient(M::AbstractManifold=DefaultManifold(2)) + return HagerZhangCoefficient(default_vector_transport_method(M)) +end update_rule_storage_points(::HagerZhangCoefficient) = Tuple{:Iterate} update_rule_storage_vectors(::HagerZhangCoefficient) = Tuple{:Gradient,:δ} @@ -355,10 +357,8 @@ Then the update reads where ``P_{a\gets b}(⋅)`` denotes a vector transport from the tangent space at ``a`` to ``b``. # Constructor - function HestenesStiefelCoefficient( - M::AbstractManifold = DefaultManifold(2); - transport_method::AbstractVectorTransportMethod=default_vector_transport_method(M) - ) + function HestenesStiefelCoefficient(transport_method::AbstractVectorTransportMethod) + function HestenesStiefelCoefficient(M::AbstractManifold = DefaultManifold(2)) Construct the Heestens Stiefel coefficient update rule, where the parallel transport is the default vector transport and a new storage is created by default. @@ -373,13 +373,13 @@ See also [`conjugate_gradient_descent`](@ref) struct HestenesStiefelCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM - function HestenesStiefelCoefficient( - M::AbstractManifold=DefaultManifold(2); - transport_method::AbstractVectorTransportMethod=default_vector_transport_method(M), - ) - return new{typeof(transport_method)}(transport_method) + function HestenesStiefelCoefficient(t::AbstractVectorTransportMethod) + return new{typeof(t)}(t) end end +function HestenesStiefelCoefficient(M::AbstractManifold=DefaultManifold(2)) + return HestenesStiefelCoefficient(default_vector_transport_method(M)) +end update_rule_storage_points(::HestenesStiefelCoefficient) = Tuple{:Iterate} update_rule_storage_vectors(::HestenesStiefelCoefficient) = Tuple{:Gradient,:δ} @@ -429,10 +429,8 @@ Then the coefficient reads See also [`conjugate_gradient_descent`](@ref) # Constructor - function LiuStoreyCoefficient( - M::AbstractManifold = DefaultManifold(2); - t::AbstractVectorTransportMethod=default_vector_transport_method(M) - ) + function LiuStoreyCoefficient(t::AbstractVectorTransportMethod) + function LiuStoreyCoefficient(M::AbstractManifold = DefaultManifold(2)) Construct the Lui Storey coefficient update rule, where the parallel transport is the default vector transport and a new storage is created by default. @@ -444,13 +442,13 @@ default vector transport and a new storage is created by default. """ struct LiuStoreyCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM - function LiuStoreyCoefficient( - M::AbstractManifold=DefaultManifold(2); - t::AbstractVectorTransportMethod=default_vector_transport_method(M), - ) + function LiuStoreyCoefficient(t::AbstractVectorTransportMethod) return new{typeof(t)}(t) end end +function LiuStoreyCoefficient(M::AbstractManifold=DefaultManifold(2)) + return LiuStoreyCoefficient(default_vector_transport_method(M)) +end update_rule_storage_points(::LiuStoreyCoefficient) = Tuple{:Iterate} update_rule_storage_vectors(::LiuStoreyCoefficient) = Tuple{:Gradient,:δ} @@ -518,13 +516,13 @@ See also [`conjugate_gradient_descent`](@ref) """ struct PolakRibiereCoefficient{TVTM<:AbstractVectorTransportMethod} <: DirectionUpdateRule transport_method::TVTM - function PolakRibiereCoefficient( - M::AbstractManifold=DefaultManifold(2); - t::AbstractVectorTransportMethod=default_vector_transport_method(M), - ) + function PolakRibiereCoefficient(t::AbstractVectorTransportMethod) return new{typeof(t)}(t) end end +function PolakRibiereCoefficient(M::AbstractManifold=DefaultManifold(2)) + return PolakRibiereCoefficient(default_vector_transport_method(M)) +end update_rule_storage_points(::PolakRibiereCoefficient) = Tuple{:Iterate} update_rule_storage_vectors(::PolakRibiereCoefficient) = Tuple{:Gradient} From e7270fb6d8504389d8608ea1723841551a12c671 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Thu, 16 Feb 2023 18:27:47 +0100 Subject: [PATCH 17/37] Set p_init and X_init as keyword arguments. --- Project.toml | 6 +++--- src/plans/conjugate_gradient_plan.jl | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 497483fd1f..b07f860371 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Manopt" uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5" authors = ["Ronny Bergmann "] -version = "0.4.5" +version = "0.4.8" [deps] ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" @@ -26,9 +26,9 @@ ColorSchemes = "3.5.0" ColorTypes = "0.9.1, 0.10, 0.11" Colors = "0.11.2, 0.12" DataStructures = "0.17, 0.18" -ManifoldDiff = "0.2.1" +ManifoldDiff = "0.2, 0.3" Manifolds = "0.8.43" -ManifoldsBase = "0.13.28" +ManifoldsBase = "0.13.28, 0.14" Requires = "0.5, 1" StaticArrays = "0.12, 1.0" julia = "1.8" diff --git a/src/plans/conjugate_gradient_plan.jl b/src/plans/conjugate_gradient_plan.jl index 58177f1116..78ac5e6025 100644 --- a/src/plans/conjugate_gradient_plan.jl +++ b/src/plans/conjugate_gradient_plan.jl @@ -6,7 +6,7 @@ end function DirectionUpdateRuleStorage( M::AbstractManifold, - dur::DirectionUpdateRule, + dur::DirectionUpdateRule; p_init=rand(M), X_init=zero_vector(M, p_init), ) @@ -68,7 +68,7 @@ mutable struct ConjugateGradientDescentState{ vtr::AbstractVectorTransportMethod=default_vector_transport_method(M), initial_gradient::T=zero_vector(M, p), ) where {P,T} - coef = DirectionUpdateRuleStorage(M, dC, p, initial_gradient) + coef = DirectionUpdateRuleStorage(M, dC; p_init=p, X_init=initial_gradient) βT = allocate_result_type(M, ConjugateGradientDescentState, (p, initial_gradient)) cgs = new{P,T,βT,typeof(coef),typeof(s),typeof(sC),typeof(retr),typeof(vtr)}() cgs.p = p From 471382a6b5731771f7c9730f13a13a07da3fe903 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sat, 18 Feb 2023 20:57:08 +0100 Subject: [PATCH 18/37] new StoreStateAction constructor --- src/plans/conjugate_gradient_plan.jl | 5 +++-- src/plans/debug.jl | 4 ++-- src/plans/solver_state.jl | 15 +++++++++++++++ 3 files changed, 20 insertions(+), 4 deletions(-) diff --git a/src/plans/conjugate_gradient_plan.jl b/src/plans/conjugate_gradient_plan.jl index 5c7dac0a17..7c32b11198 100644 --- a/src/plans/conjugate_gradient_plan.jl +++ b/src/plans/conjugate_gradient_plan.jl @@ -13,7 +13,9 @@ function DirectionUpdateRuleStorage( ursp = update_rule_storage_points(dur) ursv = update_rule_storage_vectors(dur) # StoreStateAction makes a copy - sa = StoreStateAction(M, Symbol[], ursp, ursv; p_init=p_init, X_init=X_init) + sa = StoreStateAction( + M; store_points=ursp, store_vectors=ursv, p_init=p_init, X_init=X_init + ) return DirectionUpdateRuleStorage{typeof(dur),typeof(sa)}(dur, sa) end @@ -610,7 +612,6 @@ The default threshold is chosen as `0.2` as recommended in [^Powell1977]. threshold=0.2; manifold::AbstractManifold = DefaultManifold(), vector_transport_method::V=default_vector_transport_method(manifold), - a::StoreStateAction=StoreStateAction((:Iterate, :Gradient, :δ)), ) [^Beale1972]: diff --git a/src/plans/debug.jl b/src/plans/debug.jl index 41404b3c0c..8edab7ba45 100644 --- a/src/plans/debug.jl +++ b/src/plans/debug.jl @@ -216,9 +216,9 @@ mutable struct DebugChange{IR<:AbstractInverseRetractionMethod} <: DebugAction # Is this how this is intended? if isnothing(storage) if M isa DefaultManifold - storage = StoreStateAction([:Iterate]) + storage = StoreStateAction(M; store_fields=[:Iterate]) else - storage = StoreStateAction(M, Symbol[], Tuple{:Iterate}) + storage = StoreStateAction(M; store_points=Tuple{:Iterate}) end end return new{typeof(irm)}(io, format, storage, irm) diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index d58861aa3e..30be629c79 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -331,6 +331,21 @@ end tangent_values = NamedTuple{TTS_tuple}(map(_ -> X_init, TTS_tuple)) return StoreStateAction(dictionary_symbols, point_values, tangent_values, once) end +@inline function StoreStateAction( + M::AbstractManifold; + store_fields::Vector{Symbol}=Symbol[], + store_points::Type{TPS}=Tuple{}, + store_vectors::Type{TTS}=Tuple{}, + p_init=rand(M), + X_init=zero_vector(M, p_init), + once=true, +) where {TPS<:Tuple,TTS<:Tuple} + TPS_tuple = Tuple(TPS.parameters) + TTS_tuple = Tuple(TTS.parameters) + point_values = NamedTuple{TPS_tuple}(map(_ -> p_init, TPS_tuple)) + tangent_values = NamedTuple{TTS_tuple}(map(_ -> X_init, TTS_tuple)) + return StoreStateAction(store_fields, point_values, tangent_values, once) +end function (a::StoreStateAction)( amp::AbstractManoptProblem, s::AbstractManoptSolverState, i::Int From 96f91994eb98a632d28dfde68bec6490131c504e Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 19 Feb 2023 14:35:17 +0100 Subject: [PATCH 19/37] remove the old complicated constructor. Trying one new constructor as well, not succeeding. --- src/plans/debug.jl | 17 +++++++++++++---- src/plans/solver_state.jl | 15 --------------- test/plans/test_debug.jl | 13 ++++++++++--- test/plans/test_storage.jl | 4 ++-- test/test_deprecated.jl | 6 ++++++ 5 files changed, 31 insertions(+), 24 deletions(-) create mode 100644 test/test_deprecated.jl diff --git a/src/plans/debug.jl b/src/plans/debug.jl index 8edab7ba45..25ef10028b 100644 --- a/src/plans/debug.jl +++ b/src/plans/debug.jl @@ -197,8 +197,8 @@ mutable struct DebugChange{IR<:AbstractInverseRetractionMethod} <: DebugAction io::IO=stdout, prefix::String="Last Change: ", format::String="$(prefix)%f", - manifold::Union{Nothing,AbstractManifold}=nothing, - invretr::Union{Nothing,AbstractInverseRetractionMethod}=nothing, + manifold::Union{Nothing,<:AbstractManifold}=nothing, + invretr::Union{Nothing,<:AbstractInverseRetractionMethod}=nothing, inverse_retraction_method::AbstractInverseRetractionMethod=default_inverse_retraction_method( M ), @@ -415,13 +415,22 @@ mutable struct DebugGradientChange{VTR<:AbstractVectorTransportMethod} <: DebugA storage::StoreStateAction vector_transport_method::VTR function DebugGradientChange( - M::AbstractManifold=DefaultManifold(2); - storage::StoreStateAction=StoreStateAction([:Gradient, :Iterate]), + M::AbstractManifold=DefaultManifold(); + storage::Union{Nothing,StoreStateAction}=nothing, io::IO=stdout, prefix::String="Last Change: ", format::String="$(prefix)%f", vector_transport_method::VTR=default_vector_transport_method(M), ) where {VTR<:AbstractVectorTransportMethod} + if isnothing(storage) + if M isa DefaultManifold + storage = StoreStateAction(M; store_fields=[:Iterate, :Gradient]) + else + storage = StoreStateAction( + M; store_points=Tuple{:Iterate}, store_tangents=Tuple{:Gradient} + ) + end + end return new{VTR}(io, format, storage, vector_transport_method) end end diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index 30be629c79..ad3897cc3b 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -316,21 +316,6 @@ mutable struct StoreStateAction{ ) end end -@inline function StoreStateAction( - M::AbstractManifold, - dictionary_symbols::Vector{Symbol}, - ::Type{TPS}, - ::Type{TTS}; - p_init=rand(M), - X_init=zero_vector(M, p_init), - once=true, -) where {TPS<:Tuple,TTS<:Tuple} - TPS_tuple = Tuple(TPS.parameters) - TTS_tuple = Tuple(TTS.parameters) - point_values = NamedTuple{TPS_tuple}(map(_ -> p_init, TPS_tuple)) - tangent_values = NamedTuple{TTS_tuple}(map(_ -> X_init, TTS_tuple)) - return StoreStateAction(dictionary_symbols, point_values, tangent_values, once) -end @inline function StoreStateAction( M::AbstractManifold; store_fields::Vector{Symbol}=Symbol[], diff --git a/test/plans/test_debug.jl b/test/plans/test_debug.jl index 2b50ad79a5..748faf7f8c 100644 --- a/test/plans/test_debug.jl +++ b/test/plans/test_debug.jl @@ -1,4 +1,4 @@ -using Manopt, Test, ManifoldsBase, Dates +using Manopt, Test, ManifoldsBase, Dates, Manifolds struct TestPolarManifold <: AbstractManifold{ℝ} end @@ -56,7 +56,7 @@ struct TestDebugAction <: DebugAction end @test String(take!(io)) == "" # Change of Iterate and recording a custom field a2 = DebugChange(; - storage=StoreStateAction(M, Symbol[], Tuple{:Iterate}, Tuple{}; p_init=p), + storage=StoreStateAction(M; store_points=Tuple{:Iterate}, p_init=p), prefix="Last: ", io=io, ) @@ -64,7 +64,7 @@ struct TestDebugAction <: DebugAction end st.p = [3.0, 2.0] a2(mp, st, 1) a2inv = DebugChange(; - storage=StoreStateAction([:Iterate]), + storage=StoreStateAction(M; store_fields=[:Iterate]), prefix="Last: ", io=io, inverse_retraction_method=PolarInverseRetraction(), @@ -305,16 +305,23 @@ struct TestDebugAction <: DebugAction end @test repr(d3) == "DebugGroup([$(d1), $(d2)])" ts = "[ $(Manopt.status_summary(d1)), $(Manopt.status_summary(d2)) ]" @test Manopt.status_summary(d3) == ts + d4 = DebugEvery(d1, 4) @test repr(d4) == "DebugEvery($(d1), 4, true)" @test Manopt.status_summary(d4) === "[$(d1), 4]" + ts2 = "DebugChange(; format=\"Last Change: %f\", inverse_retraction=LogarithmicInverseRetraction())" @test repr(DebugChange()) == ts2 @test Manopt.status_summary(DebugChange()) == "(:Change, \"Last Change: %f\")" + # check that a nondefault manifold works as well - not sure how to test this then + d = DebugChange(Euclidean(2)) + @test repr(DebugCost()) == "DebugCost(; format=\"F(x): %f\")" @test Manopt.status_summary(DebugCost()) == "(:Cost, \"F(x): %f\")" + @test repr(DebugDivider("|")) == "DebugDivider(; divider=\"|\")" @test Manopt.status_summary(DebugDivider("a")) == "\"a\"" + @test repr(DebugEntry(:a)) == "DebugEntry(:a; format=\"a: %s\")" end end diff --git a/test/plans/test_storage.jl b/test/plans/test_storage.jl index 288d96cb0a..ee1a3af6d1 100644 --- a/test/plans/test_storage.jl +++ b/test/plans/test_storage.jl @@ -10,7 +10,7 @@ using Test, Manopt, ManifoldsBase grad_f(M, q) = -2 * log(M, q, p) mp = DefaultManoptProblem(M, ManifoldGradientObjective(f, grad_f)) - a = StoreStateAction(M, [:p, :X], Tuple{}, Tuple{}) + a = StoreStateAction(M; store_fields=[:p, :X]) @test !has_storage(a, Manopt.PointStorageKey(:p)) @test !has_storage(a, Manopt.TangentStorageKey(:X)) @@ -20,7 +20,7 @@ using Test, Manopt, ManifoldsBase @test get_storage(a, Manopt.PointStorageKey(:p)) == p @test get_storage(a, Manopt.TangentStorageKey(:X)) == [0.0, 0.0] - a2 = StoreStateAction(M, Symbol[], Tuple{:p}, Tuple{:X}) + a2 = StoreStateAction(M; store_points=Tuple{:p}, store_vectors=Tuple{:X}) @test !has_storage(a2, Manopt.PointStorageKey(:p)) @test !has_storage(a2, Manopt.TangentStorageKey(:X)) update_storage!(a2, mp, st) diff --git a/test/test_deprecated.jl b/test/test_deprecated.jl new file mode 100644 index 0000000000..761ae0c406 --- /dev/null +++ b/test/test_deprecated.jl @@ -0,0 +1,6 @@ +using Manopt, ManifoldsBase, Test + +@testset "test deprecated definitions still work" begin + @test_logs (:warn,) DebugChange(; invretr=LogarithmicInverseRetraction()) + @test_logs (:warn,) DebugChange(; manifold=ManifoldsBase.DefaultManifold()) +end From 1b4c1084de4314ac432b49f377e1aa94d9c645ed Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 19 Feb 2023 15:02:55 +0100 Subject: [PATCH 20/37] Adapt RecordChange. --- src/plans/debug.jl | 11 ++++------- src/plans/record.jl | 30 +++++++++++++++++++++++------- test/plans/test_record.jl | 2 +- test/test_deprecated.jl | 1 + 4 files changed, 29 insertions(+), 15 deletions(-) diff --git a/src/plans/debug.jl b/src/plans/debug.jl index 25ef10028b..3f36f0a88c 100644 --- a/src/plans/debug.jl +++ b/src/plans/debug.jl @@ -171,7 +171,7 @@ end # Special single ones # @doc raw""" - DebugChange() + DebugChange(M=DefaultManifold()) debug for the amount of change of the iterate (stored in `get_iterate(o)` of the [`AbstractManoptSolverState`](@ref)) during the last iteration. See [`DebugEntryChange`](@ref) for the general case @@ -181,9 +181,7 @@ during the last iteration. See [`DebugEntryChange`](@ref) for the general case * `prefix` – (`"Last Change:"`) prefix of the debug output (ignored if you set `format`) * `io` – (`stdout`) default steream to print the debug to. * `format` - ( `"$prefix %f"`) format to print the output using an sprintf format. -* `manifold` (`DefaultManifold(1)`) manifold whose default inverse retraction should be used - for approximating the distance. -* `invretr` - (`default_inverse_retraction_method(manifold)`) the inverse retraction to be +* `inverse_retraction_method` - (`default_inverse_retraction_method(M)`) the inverse retraction to be used for approximating distance. """ mutable struct DebugChange{IR<:AbstractInverseRetractionMethod} <: DebugAction @@ -197,8 +195,8 @@ mutable struct DebugChange{IR<:AbstractInverseRetractionMethod} <: DebugAction io::IO=stdout, prefix::String="Last Change: ", format::String="$(prefix)%f", - manifold::Union{Nothing,<:AbstractManifold}=nothing, - invretr::Union{Nothing,<:AbstractInverseRetractionMethod}=nothing, + manifold::Union{Nothing,AbstractManifold}=nothing, + invretr::Union{Nothing,AbstractInverseRetractionMethod}=nothing, inverse_retraction_method::AbstractInverseRetractionMethod=default_inverse_retraction_method( M ), @@ -213,7 +211,6 @@ mutable struct DebugChange{IR<:AbstractInverseRetractionMethod} <: DebugAction @warn "invretr keyword is deprecated, use `inverse_retraction_method`, which this one overrides for now." irm = invretr end - # Is this how this is intended? if isnothing(storage) if M isa DefaultManifold storage = StoreStateAction(M; store_fields=[:Iterate]) diff --git a/src/plans/record.jl b/src/plans/record.jl index 56a41baaf8..03500f16d4 100644 --- a/src/plans/record.jl +++ b/src/plans/record.jl @@ -368,20 +368,36 @@ during the last iteration. * `inverse_retraction_method` - (`default_inverse_retraction_method(manifold, p)`) the inverse retraction to be used for approximating distance. -# Additional constructor keyword parameters -* `manifold` (`DefaultManifold(1)`) manifold whose default inverse retraction should be used -for approximating the distance. +# Constructor + + RecordChange(M=DefaultManifold();) + +with the above fields as keywords. For the `DefaultManifold` only the field storage is used. +Providing the actual manifold moves the default storage to the efficient point storage. """ mutable struct RecordChange{TInvRetr<:AbstractInverseRetractionMethod} <: RecordAction recorded_values::Vector{Float64} storage::StoreStateAction inverse_retraction_method::TInvRetr function RecordChange( - a::StoreStateAction=StoreStateAction([:Iterate]); - manifold::AbstractManifold=DefaultManifold(1), - inverse_retraction_method::IRT=default_inverse_retraction_method(manifold), + M::AbstractManifold=DefaultManifold(); + storage::Union{Nothing,StoreStateAction}=nothing, + manifold::Union{Nothing,AbstractManifold}=nothing, + inverse_retraction_method::IRT=default_inverse_retraction_method(M), ) where {IRT<:AbstractInverseRetractionMethod} - return new{IRT}(Vector{Float64}(), a, inverse_retraction_method) + irm = inverse_retraction_method + if !isnothing(manifold) + @warn "The `manifold` keyword is deprecated, use the first positional argument `M`. This keyword for now sets `inverse_retracion_method`." + irm = default_inverse_retraction_method(manifold) + end + if isnothing(storage) + if M isa DefaultManifold + storage = StoreStateAction(M; store_fields=[:Iterate]) + else + storage = StoreStateAction(M; store_points=Tuple{:Iterate}) + end + end + return new{typeof(irm)}(Vector{Float64}(), storage, irm) end function RecordChange( p, diff --git a/test/plans/test_record.jl b/test/plans/test_record.jl index 6bdbca8f43..3fc69421c1 100644 --- a/test/plans/test_record.jl +++ b/test/plans/test_record.jl @@ -123,7 +123,7 @@ using Manifolds, Manopt, Test, ManifoldsBase, Dates @test e.recorded_values == [1.0] # no p0 -> assume p is the first iterate dinvretr = RecordChange(; inverse_retraction_method=PolarInverseRetraction()) - dmani = RecordChange(; manifold=Symplectic(2)) + dmani = RecordChange(Symplectic(2)) @test dinvretr.inverse_retraction_method === PolarInverseRetraction() @test dmani.inverse_retraction_method === CayleyInverseRetraction() @test d.inverse_retraction_method === LogarithmicInverseRetraction() diff --git a/test/test_deprecated.jl b/test/test_deprecated.jl index 761ae0c406..59cf30ac64 100644 --- a/test/test_deprecated.jl +++ b/test/test_deprecated.jl @@ -3,4 +3,5 @@ using Manopt, ManifoldsBase, Test @testset "test deprecated definitions still work" begin @test_logs (:warn,) DebugChange(; invretr=LogarithmicInverseRetraction()) @test_logs (:warn,) DebugChange(; manifold=ManifoldsBase.DefaultManifold()) + @test_logs (:warn,) RecordChange(; manifold=ManifoldsBase.DefaultManifold()) end From 9b93cc5535374c2df7a6a007d6a496bb0d3f74f9 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sun, 19 Feb 2023 15:10:59 +0100 Subject: [PATCH 21/37] StoreStateAction accepts vectors of symbols in kwargs --- src/plans/solver_state.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index ad3897cc3b..8808463dcb 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -319,14 +319,22 @@ end @inline function StoreStateAction( M::AbstractManifold; store_fields::Vector{Symbol}=Symbol[], - store_points::Type{TPS}=Tuple{}, - store_vectors::Type{TTS}=Tuple{}, + store_points::Union{Type{TPS},Vector{Symbol}}=Tuple{}, + store_vectors::Union{Type{TTS},Vector{Symbol}}=Tuple{}, p_init=rand(M), X_init=zero_vector(M, p_init), once=true, ) where {TPS<:Tuple,TTS<:Tuple} - TPS_tuple = Tuple(TPS.parameters) - TTS_tuple = Tuple(TTS.parameters) + if store_points isa Vector{Symbol} + TPS_tuple = tuple(store_points...) + else + TPS_tuple = Tuple(TPS.parameters) + end + if store_vectors isa Vector{Symbol} + TTS_tuple = tuple(store_vectors...) + else + TTS_tuple = Tuple(TTS.parameters) + end point_values = NamedTuple{TPS_tuple}(map(_ -> p_init, TPS_tuple)) tangent_values = NamedTuple{TTS_tuple}(map(_ -> X_init, TTS_tuple)) return StoreStateAction(store_fields, point_values, tangent_values, once) From 2acf3f051ff69490b0797007d6367c725ddebb9d Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 19 Feb 2023 15:35:01 +0100 Subject: [PATCH 22/37] adapt the remaining ones. --- src/plans/stepsize.jl | 15 ++++++++++++--- src/plans/stopping_criterion.jl | 15 ++++++++++++++- 2 files changed, 26 insertions(+), 4 deletions(-) diff --git a/src/plans/stepsize.jl b/src/plans/stepsize.jl index e979277ee6..7443934433 100644 --- a/src/plans/stepsize.jl +++ b/src/plans/stepsize.jl @@ -377,7 +377,7 @@ and ``γ`` is the sufficient decrease parameter ``∈(0,1)``. We can then find t * `max_stepsize` – (`1e3`) upper bound for the Barzilai-Borwein step size greater than min_stepsize * `retraction_method` – (`ExponentialRetraction()`) the rectraction to use * `strategy` – (`direct`) defines if the new step size is computed using the direct, indirect or alternating strategy -* `storage` – (`x`, `gradient`) a [`StoreStateAction`](@ref) to store `old_x` and `old_gradient`, the x-value and corresponding gradient of the previous iteration +* `storage` – (for `:Iterate` and `:Gradient`) a [`StoreStateAction`](@ref) * `stepsize_reduction` – (`0.5`) step size reduction factor contained in the interval (0,1) * `sufficient_decrease` – (`1e-4`) sufficient decrease parameter contained in the interval (0,1) * `vector_transport_method` – (`ParallelTransport()`) the vector transport method to use @@ -413,7 +413,7 @@ mutable struct NonmonotoneLinesearch{ storage::TSSA linesearch_stopsize::Float64 function NonmonotoneLinesearch( - M::AbstractManifold=DefaultManifold(2); + M::AbstractManifold=DefaultManifold(); initial_stepsize::Float64=1.0, retraction_method::AbstractRetractionMethod=default_retraction_method(M), vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method( @@ -425,7 +425,7 @@ mutable struct NonmonotoneLinesearch{ min_stepsize::Float64=1e-3, max_stepsize::Float64=1e3, strategy::Symbol=:direct, - storage::StoreStateAction=StoreStateAction([:Iterate, :Gradient]), + storage::Union{Nothing,StoreStateAction}=nothing, linesearch_stopsize::Float64=0.0, ) if strategy ∉ [:direct, :inverse, :alternating] @@ -455,6 +455,15 @@ mutable struct NonmonotoneLinesearch{ if memory_size <= 0 throw(DomainError(memory_size, "The memory_size has to be greater than zero.")) end + if isnothing(storage) + if (M isa DefaultManifold) || (rand(M) isa Number) + storage = StoreStateAction(M; store_fields=[:Iterate, :Gradient]) + else + storage = StoreStateAction( + M; store_points=Tuple{:Iterate}, store_vectors=Tuple{:Gradient} + ) + end + end return new{ typeof(retraction_method), typeof(vector_transport_method), diff --git a/src/plans/stopping_criterion.jl b/src/plans/stopping_criterion.jl index a8497505c2..727f2794f1 100644 --- a/src/plans/stopping_criterion.jl +++ b/src/plans/stopping_criterion.jl @@ -175,12 +175,25 @@ mutable struct StopWhenChangeLess{ inverse_retraction::IRT at_iteration::Int end +function StopWhenChangeLess( + M::AbstractManifold, + ε::Float64; + storage::StoreStateAction=StoreStateAction(M; store_points=Tuple{:Iterate}), + inverse_retraction_method::IRT=default_inverse_retraction_method(M), +) where {IRT<:AbstractInverseRetractionMethod} + return StopWhenChangeLess{IRT,typeof(storage)}( + ε, "", storage, inverse_retraction_method, 0 + ) +end function StopWhenChangeLess( ε::Float64; storage::StoreStateAction=StoreStateAction([:Iterate]), - manifold::AbstractManifold=DefaultManifold(3), + manifold::AbstractManifold=DefaultManifold(), inverse_retraction_method::IRT=default_inverse_retraction_method(manifold), ) where {IRT<:AbstractInverseRetractionMethod} + if !(manifold isa DefaultManifold) + @warn "The `manifold` keyword is deprecated, use the first positional argument `M`. This keyword for now sets `inverse_retracion_method`." + end return StopWhenChangeLess{IRT,typeof(storage)}( ε, "", storage, inverse_retraction_method, 0 ) From 0d00e83eb44ac3657b5a532a02627a52ebda0a0e Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 19 Feb 2023 15:43:54 +0100 Subject: [PATCH 23/37] Update Documentation of StoreStateAction. --- src/plans/solver_state.jl | 42 +++++++++++++-------------------------- 1 file changed, 14 insertions(+), 28 deletions(-) diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index 8808463dcb..228e5c44d6 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -239,38 +239,24 @@ For the point storage use `PointStorageKey`. For tangent vector storage use # Constructiors - AbstractStateAction([keys=(), once=true]) + StoreStateAction(M) -Initialize the Functor to an (empty) set of keys, where `once` determines -whether more that one update per iteration are effective +## Keyword arguments - function StoreStateAction( - M::AbstractManifold, - dictionary_symbols::Vector{Symbol}, - ::Type{TPS}, - ::Type{TTS}; - p_init=rand(M), - X_init=zero_vector(M, p_init), - once=true, - ) where {TPS<:Tuple,TTS<:Tuple} - -Initialize the general storage keys to `dictionary_symbols`, point storage keys to `TPS` and -tangent vector storage tu `TTS`. For example you may call -`StorageStateAction(M, Symbol[], Tuple{:Iterate}, Tuple{:Gradient})` to create efficient -storage for point representing iterate and tangent vector representing gradient. +* `store_fields` (`Symbol[]`) +* `store_points` (`Symbol[]`) +* `store_vectors` (`Symbol[]`) - function StoreStateAction( - general_keys::Vector{Symbol}=Symbol[], - point_values::NamedTuple=NamedTuple(), - tangent_values::NamedTuple=NamedTuple(), - once::Bool=true, - )) +as vectors of symbols each referring to fields of the state (lower case symbols) +or semantic ones (upper case). + +* `p_init` (`rand(M)`) +* `X_init` (`zero_vector(M, p_init)`) + +are used to initialize the point and vector storages, change these if you use other +types (than the default) for your points/vectors on `M`. -Initialize the Functor to a set of keys, where the dictionary is initialized to -be empty. Further, `once` determines whether more that one update per iteration -are effective, otherwise only the first update is stored, all others are ignored. -Make a copy of points and tangent vectors passed to `point_values` and `tangent_values` -for later storage respective fields. +* `once` (`true`) whether to update internal storage only once per iteration or on every update call """ mutable struct StoreStateAction{ TPS<:NamedTuple,TXS<:NamedTuple,TPI<:NamedTuple,TTI<:NamedTuple From 400bc3a6971e4701d7dba24030627fdc6525cb8c Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 19 Feb 2023 15:56:33 +0100 Subject: [PATCH 24/37] includes deprecated tests. --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 5665d324a8..28da4963d6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -54,4 +54,5 @@ using Manifolds, Manopt, ManifoldsBase, Test include("solvers/test_trust_regions.jl") include("solvers/test_trust_regions_hessian_update.jl") end + include("test_deprecated.jl") end From ee6dce90136d4d03fbd9cc48aa1f0d08b5b0e854 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sun, 19 Feb 2023 18:13:54 +0100 Subject: [PATCH 25/37] fix storage for circle --- src/plans/solver_state.jl | 58 +++++++++++++++++++++++++++---- test/plans/test_storage.jl | 71 ++++++++++++++++++++++++-------------- 2 files changed, 97 insertions(+), 32 deletions(-) diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index 228e5c44d6..1d11e8f85d 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -202,6 +202,22 @@ for example within the [`DebugSolverState`](@ref) or within the [`RecordSolverSt """ abstract type AbstractStateAction end +""" + _storage_copy_point(M::AbstractManifold, p) + +Make a copy of point `p` from manifold `M` for storage in [`StoreStateAction`](@ref). +""" +_storage_copy_point(M::AbstractManifold, p) = copy(M, p) +_storage_copy_point(::AbstractManifold, p::Number) = fill(p) + +""" + _storage_copy_vector(M::AbstractManifold, X) + +Make a copy of tangent vector `X` from manifold `M` for storage in [`StoreStateAction`](@ref). +""" +_storage_copy_vector(M::AbstractManifold, X) = copy(M, SA_F64[], X) +_storage_copy_vector(::AbstractManifold, X::Number) = fill(X) + @doc raw""" StoreStateAction <: AbstractStateAction @@ -259,7 +275,7 @@ types (than the default) for your points/vectors on `M`. * `once` (`true`) whether to update internal storage only once per iteration or on every update call """ mutable struct StoreStateAction{ - TPS<:NamedTuple,TXS<:NamedTuple,TPI<:NamedTuple,TTI<:NamedTuple + TPS_asserts,TXS_assert,TPS<:NamedTuple,TXS<:NamedTuple,TPI<:NamedTuple,TTI<:NamedTuple } <: AbstractStateAction values::Dict{Symbol,Any} keys::Vector{Symbol} # for values @@ -273,21 +289,24 @@ mutable struct StoreStateAction{ general_keys::Vector{Symbol}=Symbol[], point_values::NamedTuple=NamedTuple(), tangent_values::NamedTuple=NamedTuple(), - once::Bool=true, + once::Bool=true; + M::AbstractManifold=DefaultManifold(), ) point_init = NamedTuple{keys(point_values)}(map(u -> false, keys(point_values))) tangent_init = NamedTuple{keys(tangent_values)}( map(u -> false, keys(tangent_values)) ) point_values_copy = NamedTuple{keys(point_values)}( - map(u -> copy(point_values[u]), keys(point_values)) + map(u -> _storage_copy_point(M, point_values[u]), keys(point_values)) ) tangent_values_copy = NamedTuple{keys(tangent_values)}( - map(u -> copy(tangent_values[u]), keys(tangent_values)) + map(u -> _storage_copy_vector(M, tangent_values[u]), keys(tangent_values)) ) return new{ typeof(point_values), typeof(tangent_values), + typeof(point_values_copy), + typeof(tangent_values_copy), typeof(point_init), typeof(tangent_init), }( @@ -323,7 +342,28 @@ end end point_values = NamedTuple{TPS_tuple}(map(_ -> p_init, TPS_tuple)) tangent_values = NamedTuple{TTS_tuple}(map(_ -> X_init, TTS_tuple)) - return StoreStateAction(store_fields, point_values, tangent_values, once) + return StoreStateAction(store_fields, point_values, tangent_values, once; M=M) +end + +@generated function extract_type_from_namedtuple(::Type{nt}, ::Val{key}) where {nt,key} + for i in 1:length(nt.parameters[1]) + if nt.parameters[1][i] === key + return nt.parameters[2].parameters[i] + end + end + return Any +end + +function _store_point_assert_type( + ::StoreStateAction{TPS_asserts,TXS_assert}, key::Val +) where {TPS_asserts,TXS_assert} + return extract_type_from_namedtuple(TPS_asserts, key) +end + +function _store_tangent_assert_type( + ::StoreStateAction{TPS_asserts,TXS_assert}, key::Val +) where {TPS_asserts,TXS_assert} + return extract_type_from_namedtuple(TXS_assert, key) end function (a::StoreStateAction)( @@ -435,7 +475,9 @@ function update_storage!( copyto!(M, a.point_values[key], get_iterate(s)) else copyto!( - M, a.point_values[key], getproperty(s, key)::typeof(a.point_values[key]) + M, + a.point_values[key], + getproperty(s, key)::_store_point_assert_type(a, Val(key)), ) end end @@ -446,7 +488,9 @@ function update_storage!( copyto!(M, a.tangent_values[key], get_gradient(s)) else copyto!( - M, a.tangent_values[key], getproperty(s, key)::typeof(a.tangent_values[key]) + M, + a.tangent_values[key], + getproperty(s, key)::_store_tangent_assert_type(a, Val(key)), ) end end diff --git a/test/plans/test_storage.jl b/test/plans/test_storage.jl index ee1a3af6d1..62cc255fb2 100644 --- a/test/plans/test_storage.jl +++ b/test/plans/test_storage.jl @@ -1,31 +1,52 @@ -using Test, Manopt, ManifoldsBase +using Test, Manopt, ManifoldsBase, Manifolds @testset "StoreStateAction" begin - M = ManifoldsBase.DefaultManifold(2) - p = [4.0, 2.0] - st = GradientDescentState( - M, p; stopping_criterion=StopAfterIteration(20), stepsize=ConstantStepsize(M) - ) - f(M, q) = distance(M, q, p) .^ 2 - grad_f(M, q) = -2 * log(M, q, p) - mp = DefaultManoptProblem(M, ManifoldGradientObjective(f, grad_f)) + @testset "manifold $M" for M in [ManifoldsBase.DefaultManifold(2), Circle()] + if M isa Circle + p = 0.4 + p_fast = fill(p) + X_zero = 0.0 + X_zero_fast = fill(0.0) + else + p = [4.0, 2.0] + p_fast = p + X_zero = [0.0, 0.0] + X_zero_fast = X_zero + end - a = StoreStateAction(M; store_fields=[:p, :X]) + st = GradientDescentState( + M, p; stopping_criterion=StopAfterIteration(20), stepsize=ConstantStepsize(M) + ) + f(M, q) = distance(M, q, p) .^ 2 + grad_f(M, q) = -2 * log(M, q, p) + mp = DefaultManoptProblem(M, ManifoldGradientObjective(f, grad_f)) - @test !has_storage(a, Manopt.PointStorageKey(:p)) - @test !has_storage(a, Manopt.TangentStorageKey(:X)) - update_storage!(a, mp, st) - @test has_storage(a, Manopt.PointStorageKey(:p)) - @test has_storage(a, Manopt.TangentStorageKey(:X)) - @test get_storage(a, Manopt.PointStorageKey(:p)) == p - @test get_storage(a, Manopt.TangentStorageKey(:X)) == [0.0, 0.0] + a = StoreStateAction(M; store_fields=[:p, :X]) - a2 = StoreStateAction(M; store_points=Tuple{:p}, store_vectors=Tuple{:X}) - @test !has_storage(a2, Manopt.PointStorageKey(:p)) - @test !has_storage(a2, Manopt.TangentStorageKey(:X)) - update_storage!(a2, mp, st) - @test has_storage(a2, Manopt.PointStorageKey(:p)) - @test has_storage(a2, Manopt.TangentStorageKey(:X)) - @test get_storage(a2, Manopt.PointStorageKey(:p)) == p - @test get_storage(a2, Manopt.TangentStorageKey(:X)) == [0.0, 0.0] + @test !has_storage(a, Manopt.PointStorageKey(:p)) + @test !has_storage(a, Manopt.TangentStorageKey(:X)) + update_storage!(a, mp, st) + @test has_storage(a, Manopt.PointStorageKey(:p)) + @test has_storage(a, Manopt.TangentStorageKey(:X)) + @test get_storage(a, Manopt.PointStorageKey(:p)) == p + @test get_storage(a, Manopt.TangentStorageKey(:X)) == X_zero + + a2 = StoreStateAction(M; store_points=Tuple{:p}, store_vectors=Tuple{:X}) + @test !has_storage(a2, Manopt.PointStorageKey(:p)) + @test !has_storage(a2, Manopt.TangentStorageKey(:X)) + update_storage!(a2, mp, st) + @test has_storage(a2, Manopt.PointStorageKey(:p)) + @test has_storage(a2, Manopt.TangentStorageKey(:X)) + @test get_storage(a2, Manopt.PointStorageKey(:p)) == p_fast + @test get_storage(a2, Manopt.TangentStorageKey(:X)) == X_zero_fast + + a3 = StoreStateAction(M; store_points=[:p], store_vectors=[:X]) + @test !has_storage(a3, Manopt.PointStorageKey(:p)) + @test !has_storage(a3, Manopt.TangentStorageKey(:X)) + update_storage!(a3, mp, st) + @test has_storage(a3, Manopt.PointStorageKey(:p)) + @test has_storage(a3, Manopt.TangentStorageKey(:X)) + @test get_storage(a3, Manopt.PointStorageKey(:p)) == p_fast + @test get_storage(a3, Manopt.TangentStorageKey(:X)) == X_zero_fast + end end From f2b8a913a55a700a9982665063efcfa1ce1fa904 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 19 Feb 2023 18:18:34 +0100 Subject: [PATCH 26/37] improve test coverage, fixes a bug. --- src/plans/debug.jl | 2 +- src/plans/solver_state.jl | 5 +++++ test/plans/test_debug.jl | 3 +++ test/plans/test_stepsize.jl | 3 ++- test/plans/test_stopping_criteria.jl | 2 ++ test/plans/test_storage.jl | 6 +++++- test/solvers/test_exact_penalty.jl | 6 +++--- test/test_deprecated.jl | 3 ++- 8 files changed, 23 insertions(+), 7 deletions(-) diff --git a/src/plans/debug.jl b/src/plans/debug.jl index 3f36f0a88c..70179e2b3b 100644 --- a/src/plans/debug.jl +++ b/src/plans/debug.jl @@ -424,7 +424,7 @@ mutable struct DebugGradientChange{VTR<:AbstractVectorTransportMethod} <: DebugA storage = StoreStateAction(M; store_fields=[:Iterate, :Gradient]) else storage = StoreStateAction( - M; store_points=Tuple{:Iterate}, store_tangents=Tuple{:Gradient} + M; store_points=[:Iterate], store_vectors=[:Gradient] ) end end diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index 228e5c44d6..e44b3c5d19 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -239,6 +239,11 @@ For the point storage use `PointStorageKey`. For tangent vector storage use # Constructiors + StoreStateAction(s::Vector{Symbol}) + +This is equivalent as providing `s` to the keyword `store_fields`, just that here, no manifold +is necessay for the construciton. + StoreStateAction(M) ## Keyword arguments diff --git a/test/plans/test_debug.jl b/test/plans/test_debug.jl index 748faf7f8c..e930806587 100644 --- a/test/plans/test_debug.jl +++ b/test/plans/test_debug.jl @@ -205,6 +205,9 @@ struct TestDebugAction <: DebugAction end dgc_s = "DebugGradientChange(; format=\"Last Change: %f\", vector_transport_method=ParallelTransport())" @test repr(dgc) == dgc_s @test Manopt.status_summary(dgc) == "(:GradientChange, \"Last Change: %f\")" + # Faster storage + dgc2 = DebugGradientChange(Euclidean(2)) + @test repr(dgc2) == dgc_s end @testset "Debug Warnings" begin diff --git a/test/plans/test_stepsize.jl b/test/plans/test_stepsize.jl index 2a4ee11b0c..2d2c843c9b 100644 --- a/test/plans/test_stepsize.jl +++ b/test/plans/test_stepsize.jl @@ -9,7 +9,8 @@ using Manopt, Manifolds, Test s2 = NonmonotoneLinesearch() @test startswith(repr(s2), "NonmonotoneLinesearch() with keyword arguments\n") - + s2b = NonmonotoneLinesearch(Euclidean(2)) # with manifold -> faster storage + @test startswith(repr(s2b), "NonmonotoneLinesearch() with keyword arguments\n") s3 = WolfePowellBinaryLinesearch() @test startswith(repr(s3), "WolfePowellBinaryLinesearch(DefaultManifold(), ") # no stepsize yet so repr and summary are the same diff --git a/test/plans/test_stopping_criteria.jl b/test/plans/test_stopping_criteria.jl index 7b5ed6120c..1752907ea4 100644 --- a/test/plans/test_stopping_criteria.jl +++ b/test/plans/test_stopping_criteria.jl @@ -69,6 +69,8 @@ end b = StopWhenChangeLess(1e-6) sb = "StopWhenChangeLess(1.0e-6)\n $(Manopt.status_summary(b))" @test repr(b) == sb + b2 = StopWhenChangeLess(Euclidean(), 1e-6) # second constructor + @test repr(b2) == sb c = StopWhenGradientNormLess(1e-6) sc = "StopWhenGradientNormLess(1.0e-6)\n $(Manopt.status_summary(c))" @test repr(c) == sc diff --git a/test/plans/test_storage.jl b/test/plans/test_storage.jl index ee1a3af6d1..ca02baa6eb 100644 --- a/test/plans/test_storage.jl +++ b/test/plans/test_storage.jl @@ -20,7 +20,7 @@ using Test, Manopt, ManifoldsBase @test get_storage(a, Manopt.PointStorageKey(:p)) == p @test get_storage(a, Manopt.TangentStorageKey(:X)) == [0.0, 0.0] - a2 = StoreStateAction(M; store_points=Tuple{:p}, store_vectors=Tuple{:X}) + a2 = StoreStateAction(M; store_points=[:p], store_vectors=[:X]) @test !has_storage(a2, Manopt.PointStorageKey(:p)) @test !has_storage(a2, Manopt.TangentStorageKey(:X)) update_storage!(a2, mp, st) @@ -28,4 +28,8 @@ using Test, Manopt, ManifoldsBase @test has_storage(a2, Manopt.TangentStorageKey(:X)) @test get_storage(a2, Manopt.PointStorageKey(:p)) == p @test get_storage(a2, Manopt.TangentStorageKey(:X)) == [0.0, 0.0] + a2b = StoreStateAction(M; store_points=Tuple{:p}, store_vectors=Tuple{:X}) + @test keys(a2.point_values) == keys(a2b.point_values) + @test keys(a2.tangent_values) == keys(a2b.tangent_values) + @test keys(a2.keys) == keys(a2b.keys) end diff --git a/test/solvers/test_exact_penalty.jl b/test/solvers/test_exact_penalty.jl index c9295f0deb..097be5c356 100644 --- a/test/solvers/test_exact_penalty.jl +++ b/test/solvers/test_exact_penalty.jl @@ -2,9 +2,9 @@ using Manopt, ManifoldsBase, Manifolds, Test using LinearAlgebra: I, tr @testset "Test REPM with a nonneg. PCA" begin - d = 20 + d = 4 M = Sphere(d - 1) - S = [ones(4)..., zeros(d - 4)...] + S = [ones(2)..., zeros(d - 2)...] v0 = project(M, S) Z = v0 * v0' f(M, p) = -tr(transpose(p) * Z * p) / 2 @@ -12,7 +12,7 @@ using LinearAlgebra: I, tr g(M, p) = -p # i.e. p ≥ 0 mI = -Matrix{Float64}(I, d, d) grad_g(M, p) = [project(M, p, mI[:, i]) for i in 1:d] - x0 = project(M, [ones(3)..., zeros(d - 3)...]) + x0 = project(M, [ones(2)..., zeros(d - 1)..., 0.1]) sol_lse = exact_penalty_method(M, f, grad_f, x0; g=g, grad_g=grad_g) sol_lqh = exact_penalty_method( M, f, grad_f, x0; g=g, grad_g=grad_g, smoothing=LinearQuadraticHuber() diff --git a/test/test_deprecated.jl b/test/test_deprecated.jl index 59cf30ac64..2937c767bc 100644 --- a/test/test_deprecated.jl +++ b/test/test_deprecated.jl @@ -1,7 +1,8 @@ -using Manopt, ManifoldsBase, Test +using Manopt, Manifolds, ManifoldsBase, Test @testset "test deprecated definitions still work" begin @test_logs (:warn,) DebugChange(; invretr=LogarithmicInverseRetraction()) @test_logs (:warn,) DebugChange(; manifold=ManifoldsBase.DefaultManifold()) @test_logs (:warn,) RecordChange(; manifold=ManifoldsBase.DefaultManifold()) + @test_logs (:warn,) StopWhenChangeLess(1e-9; manifold=Euclidean()) end From 013e3c6ebbdadf9571d1237af5b08d54113924f3 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 19 Feb 2023 18:22:13 +0100 Subject: [PATCH 27/37] stabilize a test. --- test/solvers/test_exact_penalty.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/solvers/test_exact_penalty.jl b/test/solvers/test_exact_penalty.jl index 097be5c356..c6cfca8198 100644 --- a/test/solvers/test_exact_penalty.jl +++ b/test/solvers/test_exact_penalty.jl @@ -12,7 +12,7 @@ using LinearAlgebra: I, tr g(M, p) = -p # i.e. p ≥ 0 mI = -Matrix{Float64}(I, d, d) grad_g(M, p) = [project(M, p, mI[:, i]) for i in 1:d] - x0 = project(M, [ones(2)..., zeros(d - 1)..., 0.1]) + x0 = project(M, [ones(2)..., zeros(d - 3)..., 0.1]) sol_lse = exact_penalty_method(M, f, grad_f, x0; g=g, grad_g=grad_g) sol_lqh = exact_penalty_method( M, f, grad_f, x0; g=g, grad_g=grad_g, smoothing=LinearQuadraticHuber() From 136b0a09bc6261eba5d3748bc25336fb9251d4ce Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 19 Feb 2023 19:06:58 +0100 Subject: [PATCH 28/37] Adds a changelog. --- Changelog.md | 78 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100644 Changelog.md diff --git a/Changelog.md b/Changelog.md new file mode 100644 index 0000000000..b64c53fe55 --- /dev/null +++ b/Changelog.md @@ -0,0 +1,78 @@ +# Changelog + +All notable Changes to the Julia package `Manopt.jl` will be 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.4.8] + +### Added + +* started this Changelog (checking the last few patches backwards) + +## [0.4.7] - 14/02/2023 + +### Changed + +* Bump [compat] entry of ManifoldDiff to also include 0.3 + +## [0.4.6] - 03/02/2023 + +### Fixed + +* Fixed a few stopping criteria even indicated to stop before the algorithm started. + +## [0.4.5] - 24/01/2023 + +### Changed + +* the new default functions that include `p` are used where possible +* a first step towards faster storage handling + +## [0.4.4] - 20/01/2023 + +### Added + +* Introduce `ConjugateGradientBealeRestart` to allow CG restarts using Beale‘s rule + +### Fixed + +* fix a type in `HestenesStiefelCoefficient` + + +## [0.4.3] - 17/01/2023 + +### Fixed + +* the CG coefficient `β` can now be complex +* fix a bug in `grad_distance` + +## [0.4.2] - 16/01/2023 + +### Changed + +* the usage of `inner` in linesearch methods, such that they work well with complex manifolds as well + + +## [0.4.1] - 15/01/2023 + +### Fixed + +* a `max_stepsize` per manifold to avoid leaving the injectivity radius, which it also defaults to + +## {0.4.0] - 10/01/2023 + +### Added + +* Dependency on `ManifoldDiff.jl` and a start of moving actual derivatives, differentials, and gradients there. +* `AbstractManifoldObjective` to store the objective within the `AbstractManoptProblem` +* Introduce a `CostGrad` structure to store a function that computes the cost and gradient within one function. + +### Changed + +* `AbstractManoptProblem` replaces `Problem` +* the problem now contains a +* `AbstractManoptSolverState` replaces `Options` +* `random_point(M)` is replaced by `rand(M)` from `ManifoldsBase.jl +* `random_tangent(M, p)` is replaced by `rand(M; vector_at=p)` \ No newline at end of file From 2fb8a32995b8ea19e3bdaa7d63bc8ede4c49cb70 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sun, 19 Feb 2023 19:22:30 +0100 Subject: [PATCH 29/37] faster storage update; test one edge case --- src/plans/solver_state.jl | 8 +++++--- test/plans/test_storage.jl | 12 ++++++++++++ 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index 3db9cf0dfd..6b12e80374 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -474,8 +474,7 @@ function update_storage!( end M = get_manifold(amp) - - map(keys(a.point_values)) do key + @inline function update_points(key) if key === :Iterate copyto!(M, a.point_values[key], get_iterate(s)) else @@ -486,9 +485,10 @@ function update_storage!( ) end end + map(update_points, keys(a.point_values)) a.point_init = NamedTuple{keys(a.point_values)}(map(u -> true, keys(a.point_values))) - map(keys(a.tangent_values)) do key + @inline function update_tangent(key) if key === :Gradient copyto!(M, a.tangent_values[key], get_gradient(s)) else @@ -499,6 +499,8 @@ function update_storage!( ) end end + map(update_tangent, keys(a.tangent_values)) + a.tangent_init = NamedTuple{keys(a.tangent_values)}( map(u -> true, keys(a.tangent_values)) ) diff --git a/test/plans/test_storage.jl b/test/plans/test_storage.jl index 9ed1a67595..1211a6a178 100644 --- a/test/plans/test_storage.jl +++ b/test/plans/test_storage.jl @@ -44,6 +44,11 @@ using Test, Manopt, ManifoldsBase, Manifolds @test keys(a2.tangent_values) == keys(a2b.tangent_values) @test keys(a2.keys) == keys(a2b.keys) + # make sure fast storage is actually fast + @test (@allocated update_storage!(a2, mp, st)) == 0 + @test (@allocated has_storage(a2, Manopt.PointStorageKey(:p))) == 0 + @test (@allocated get_storage(a2, Manopt.PointStorageKey(:p))) == 0 + a3 = StoreStateAction(M; store_points=[:p], store_vectors=[:X]) @test !has_storage(a3, Manopt.PointStorageKey(:p)) @test !has_storage(a3, Manopt.TangentStorageKey(:X)) @@ -52,5 +57,12 @@ using Test, Manopt, ManifoldsBase, Manifolds @test has_storage(a3, Manopt.TangentStorageKey(:X)) @test get_storage(a3, Manopt.PointStorageKey(:p)) == p_fast @test get_storage(a3, Manopt.TangentStorageKey(:X)) == X_zero_fast + + # make sure fast storage is actually fast + @test (@allocated update_storage!(a3, mp, st)) == 0 + @test (@allocated has_storage(a3, Manopt.PointStorageKey(:p))) == 0 + @test (@allocated get_storage(a3, Manopt.PointStorageKey(:p))) == 0 end + + @test Manopt.extract_type_from_namedtuple(typeof((; a=10, b='a')), Val(:c)) === Any end From 934181645c456c6327cf6fbe80754df1e6cf3431 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Sun, 19 Feb 2023 19:26:37 +0100 Subject: [PATCH 30/37] Update Changelog. --- Changelog.md | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/Changelog.md b/Changelog.md index b64c53fe55..e7eb034817 100644 --- a/Changelog.md +++ b/Changelog.md @@ -9,7 +9,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added -* started this Changelog (checking the last few patches backwards) +* a `status_summary` that displays the main parameters within several structures of Manopt, + most prominently a solver state + +### Changed + +* Improved storage performance by introducing separate named tuples for points and vectors +* changed the `show` methods of `AbstractManoptSolverState`s to display their `state_summary +* Move tutorials to be rendered with Quarto into the documentation. ## [0.4.7] - 14/02/2023 From 568c51fa0f7e4695e3844b8d669c1dffd2fd22eb Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Sun, 19 Feb 2023 20:55:22 +0100 Subject: [PATCH 31/37] small cleanup --- src/plans/conjugate_gradient_plan.jl | 48 +++++++++--------- src/plans/primal_dual_plan.jl | 24 ++++----- src/plans/record.jl | 22 ++++---- src/plans/solver_state.jl | 76 +++++++++++++--------------- src/plans/stepsize.jl | 2 +- test/plans/test_storage.jl | 20 ++++---- 6 files changed, 94 insertions(+), 98 deletions(-) diff --git a/src/plans/conjugate_gradient_plan.jl b/src/plans/conjugate_gradient_plan.jl index 7c32b11198..6ff40f91b9 100644 --- a/src/plans/conjugate_gradient_plan.jl +++ b/src/plans/conjugate_gradient_plan.jl @@ -134,12 +134,12 @@ function (u::DirectionUpdateRuleStorage{ConjugateDescentCoefficient})( ) M = get_manifold(amp) if !has_storage(u.storage, PointStorageKey(:Iterate)) || - !has_storage(u.storage, TangentStorageKey(:Gradient)) + !has_storage(u.storage, VectorStorageKey(:Gradient)) update_storage!(u.storage, amp, cgs) # if not given store current as old return 0.0 end p_old = get_storage(u.storage, PointStorageKey(:Iterate)) - X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) + X_old = get_storage(u.storage, VectorStorageKey(:Gradient)) coef = inner(M, cgs.p, cgs.X, cgs.X) / inner(M, p_old, -cgs.δ, X_old) update_storage!(u.storage, amp, cgs) return coef @@ -199,14 +199,14 @@ function (u::DirectionUpdateRuleStorage{<:DaiYuanCoefficient})( ) M = get_manifold(amp) if !has_storage(u.storage, PointStorageKey(:Iterate)) || - !has_storage(u.storage, TangentStorageKey(:Gradient)) || - !has_storage(u.storage, TangentStorageKey(:δ)) + !has_storage(u.storage, VectorStorageKey(:Gradient)) || + !has_storage(u.storage, VectorStorageKey(:δ)) update_storage!(u.storage, amp, cgs) # if not given store current as old return 0.0 end p_old = get_storage(u.storage, PointStorageKey(:Iterate)) - X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) - δ_old = get_storage(u.storage, TangentStorageKey(:δ)) + X_old = get_storage(u.storage, VectorStorageKey(:Gradient)) + δ_old = get_storage(u.storage, VectorStorageKey(:δ)) gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) ν = cgs.X - gradienttr #notation y from [HZ06] @@ -254,11 +254,11 @@ function (u::DirectionUpdateRuleStorage{FletcherReevesCoefficient})( ) M = get_manifold(amp) if !has_storage(u.storage, PointStorageKey(:Iterate)) || - !has_storage(u.storage, TangentStorageKey(:Gradient)) + !has_storage(u.storage, VectorStorageKey(:Gradient)) update_storage!(u.storage, amp, cgs) # if not given store current as old end p_old = get_storage(u.storage, PointStorageKey(:Iterate)) - X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) + X_old = get_storage(u.storage, VectorStorageKey(:Gradient)) coef = inner(M, cgs.p, cgs.X, cgs.X) / inner(M, p_old, X_old, X_old) update_storage!(u.storage, amp, cgs) return coef @@ -321,14 +321,14 @@ function (u::DirectionUpdateRuleStorage{<:HagerZhangCoefficient})( ) M = get_manifold(amp) if !has_storage(u.storage, PointStorageKey(:Iterate)) || - !has_storage(u.storage, TangentStorageKey(:Gradient)) || - !has_storage(u.storage, TangentStorageKey(:δ)) + !has_storage(u.storage, VectorStorageKey(:Gradient)) || + !has_storage(u.storage, VectorStorageKey(:δ)) update_storage!(u.storage, amp, cgs) # if not given store current as old return 0.0 end p_old = get_storage(u.storage, PointStorageKey(:Iterate)) - X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) - δ_old = get_storage(u.storage, TangentStorageKey(:δ)) + X_old = get_storage(u.storage, VectorStorageKey(:Gradient)) + δ_old = get_storage(u.storage, VectorStorageKey(:δ)) gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) ν = cgs.X - gradienttr #notation y from [HZ06] @@ -401,14 +401,14 @@ function (u::DirectionUpdateRuleStorage{<:HestenesStiefelCoefficient})( ) M = get_manifold(amp) if !has_storage(u.storage, PointStorageKey(:Iterate)) || - !has_storage(u.storage, TangentStorageKey(:Gradient)) || - !has_storage(u.storage, TangentStorageKey(:δ)) + !has_storage(u.storage, VectorStorageKey(:Gradient)) || + !has_storage(u.storage, VectorStorageKey(:δ)) update_storage!(u.storage, amp, cgs) # if not given store current as old return 0.0 end p_old = get_storage(u.storage, PointStorageKey(:Iterate)) - X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) - δ_old = get_storage(u.storage, TangentStorageKey(:δ)) + X_old = get_storage(u.storage, VectorStorageKey(:Gradient)) + δ_old = get_storage(u.storage, VectorStorageKey(:δ)) gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) δtr = vector_transport_to(M, p_old, δ_old, cgs.p, u.coefficient.transport_method) @@ -473,13 +473,13 @@ function (u::DirectionUpdateRuleStorage{<:LiuStoreyCoefficient})( ) M = get_manifold(amp) if !has_storage(u.storage, PointStorageKey(:Iterate)) || - !has_storage(u.storage, TangentStorageKey(:Gradient)) || - !has_storage(u.storage, TangentStorageKey(:δ)) + !has_storage(u.storage, VectorStorageKey(:Gradient)) || + !has_storage(u.storage, VectorStorageKey(:δ)) update_storage!(u.storage, amp, cgs) # if not given store current as old end p_old = get_storage(u.storage, PointStorageKey(:Iterate)) - X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) - δ_old = get_storage(u.storage, TangentStorageKey(:δ)) + X_old = get_storage(u.storage, VectorStorageKey(:Gradient)) + δ_old = get_storage(u.storage, VectorStorageKey(:δ)) gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) ν = cgs.X - gradienttr # notation y from [HZ06] coef = inner(M, cgs.p, cgs.X, ν) / inner(M, p_old, -δ_old, X_old) @@ -550,11 +550,11 @@ function (u::DirectionUpdateRuleStorage{<:PolakRibiereCoefficient})( ) M = get_manifold(amp) if !has_storage(u.storage, PointStorageKey(:Iterate)) || - !has_storage(u.storage, TangentStorageKey(:Gradient)) + !has_storage(u.storage, VectorStorageKey(:Gradient)) update_storage!(u.storage, amp, cgs) # if not given store current as old end p_old = get_storage(u.storage, PointStorageKey(:Iterate)) - X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) + X_old = get_storage(u.storage, VectorStorageKey(:Gradient)) gradienttr = vector_transport_to(M, p_old, X_old, cgs.p, u.coefficient.transport_method) ν = cgs.X - gradienttr @@ -661,11 +661,11 @@ function (u::DirectionUpdateRuleStorage{<:ConjugateGradientBealeRestart})( ) M = get_manifold(amp) if !has_storage(u.storage, PointStorageKey(:Iterate)) || - !has_storage(u.storage, TangentStorageKey(:Gradient)) + !has_storage(u.storage, VectorStorageKey(:Gradient)) update_storage!(u.storage, amp, cgs) # if not given store current as old end p_old = get_storage(u.storage, PointStorageKey(:Iterate)) - X_old = get_storage(u.storage, TangentStorageKey(:Gradient)) + X_old = get_storage(u.storage, VectorStorageKey(:Gradient)) # call actual rule β = u.coefficient.direction_update(amp, cgs, i) diff --git a/src/plans/primal_dual_plan.jl b/src/plans/primal_dual_plan.jl index 557c4991b1..691b852d8e 100644 --- a/src/plans/primal_dual_plan.jl +++ b/src/plans/primal_dual_plan.jl @@ -644,7 +644,7 @@ mutable struct DebugDualResidual <: DebugAction format::String storage::StoreStateAction function DebugDualResidual(; - storage::StoreStateAction=StoreStateAction((:Iterate, :X, :n)), + storage::StoreStateAction=StoreStateAction([:Iterate, :X, :n]), io::IO=stdout, prefix="Dual Residual: ", format="$prefix%s", @@ -653,7 +653,7 @@ mutable struct DebugDualResidual <: DebugAction end function DebugDualResidual( initial_values::Tuple{P,T,Q}; - storage::StoreStateAction=StoreStateAction((:Iterate, :X, :n)), + storage::StoreStateAction=StoreStateAction([:Iterate, :X, :n]), io::IO=stdout, prefix="Dual Residual: ", format="$prefix%s", @@ -705,7 +705,7 @@ mutable struct DebugPrimalResidual <: DebugAction format::String storage::StoreStateAction function DebugPrimalResidual(; - storage::StoreStateAction=StoreStateAction((:Iterate, :X, :n)), + storage::StoreStateAction=StoreStateAction([:Iterate, :X, :n]), io::IO=stdout, prefix="Primal Residual: ", format="$prefix%s", @@ -714,7 +714,7 @@ mutable struct DebugPrimalResidual <: DebugAction end function DebugPrimalResidual( values::Tuple{P,T,Q}; - storage::StoreStateAction=StoreStateAction((:Iterate, :X, :n)), + storage::StoreStateAction=StoreStateAction([:Iterate, :X, :n]), io::IO=stdout, prefix="Primal Residual: ", format="$prefix%s", @@ -764,7 +764,7 @@ mutable struct DebugPrimalDualResidual <: DebugAction format::String storage::StoreStateAction function DebugPrimalDualResidual(; - storage::StoreStateAction=StoreStateAction((:Iterate, :X, :n)), + storage::StoreStateAction=StoreStateAction([:Iterate, :X, :n]), io::IO=stdout, prefix="PD Residual: ", format="$prefix%s", @@ -773,7 +773,7 @@ mutable struct DebugPrimalDualResidual <: DebugAction end function DebugPrimalDualResidual( values::Tuple{P,T,Q}; - storage::StoreStateAction=StoreStateAction((:Iterate, :X, :n)), + storage::StoreStateAction=StoreStateAction([:Iterate, :X, :n]), io::IO=stdout, prefix="PD Residual: ", format="$prefix%s", @@ -811,7 +811,7 @@ Print the change of the primal variable by using [`DebugChange`](@ref), see their constructors for detail. """ function DebugPrimalChange(; - storage::StoreStateAction=StoreStateAction((:Iterate,)), + storage::StoreStateAction=StoreStateAction([:Iterate]), prefix="Primal Change: ", kwargs..., ) @@ -847,7 +847,7 @@ mutable struct DebugDualChange <: DebugAction format::String storage::StoreStateAction function DebugDualChange(; - storage::StoreStateAction=StoreStateAction((:X, :n)), + storage::StoreStateAction=StoreStateAction([:X, :n]), io::IO=stdout, prefix="Dual Change: ", format="$prefix%s", @@ -856,7 +856,7 @@ mutable struct DebugDualChange <: DebugAction end function DebugDualChange( values::Tuple{T,P}; - storage::StoreStateAction=StoreStateAction((:X, :n)), + storage::StoreStateAction=StoreStateAction([:X, :n]), io::IO=stdout, prefix="Dual Change: ", format="$prefix%s", @@ -897,13 +897,13 @@ This method is further set display `o.n`. DebugDualBaseIterate(; kwargs...) = DebugEntry(:n; kwargs...) """ - DebugDualChange(; storage=StoreStateAction((:X)), io::IO=stdout) + DebugDualChange(; storage=StoreStateAction([:n]), io::IO=stdout) Print the change of the dual base variable by using [`DebugEntryChange`](@ref), see their constructors for detail, on `o.n`. """ function DebugDualBaseChange(; - storage::StoreStateAction=StoreStateAction((:n)), prefix="Dual Base Change:", kwargs... + storage::StoreStateAction=StoreStateAction([:n]), prefix="Dual Base Change:", kwargs... ) return DebugEntryChange( :n, @@ -925,7 +925,7 @@ This method is further set display `o.m`. DebugPrimalBaseIterate(opts...; kwargs...) = DebugEntry(:m, opts...; kwargs...) """ - DebugPrimalBaseChange(a::StoreStateAction=StoreStateAction((:m)),io::IO=stdout) + DebugPrimalBaseChange(a::StoreStateAction=StoreStateAction([:m]),io::IO=stdout) Print the change of the primal base variable by using [`DebugEntryChange`](@ref), see their constructors for detail, on `o.n`. diff --git a/src/plans/record.jl b/src/plans/record.jl index 03500f16d4..e96ab3ff8c 100644 --- a/src/plans/record.jl +++ b/src/plans/record.jl @@ -375,9 +375,11 @@ during the last iteration. with the above fields as keywords. For the `DefaultManifold` only the field storage is used. Providing the actual manifold moves the default storage to the efficient point storage. """ -mutable struct RecordChange{TInvRetr<:AbstractInverseRetractionMethod} <: RecordAction +mutable struct RecordChange{ + TInvRetr<:AbstractInverseRetractionMethod,TStorage<:StoreStateAction +} <: RecordAction recorded_values::Vector{Float64} - storage::StoreStateAction + storage::TStorage inverse_retraction_method::TInvRetr function RecordChange( M::AbstractManifold=DefaultManifold(); @@ -397,7 +399,7 @@ mutable struct RecordChange{TInvRetr<:AbstractInverseRetractionMethod} <: Record storage = StoreStateAction(M; store_points=Tuple{:Iterate}) end end - return new{typeof(irm)}(Vector{Float64}(), storage, irm) + return new{typeof(irm),typeof(storage)}(Vector{Float64}(), storage, irm) end function RecordChange( p, @@ -408,7 +410,7 @@ mutable struct RecordChange{TInvRetr<:AbstractInverseRetractionMethod} <: Record ), ) where {IRT<:AbstractInverseRetractionMethod} update_storage!(a, Dict(:Iterate => p)) - return new{IRT}(Vector{Float64}(), a, inverse_retraction_method) + return new{IRT,typeof(a)}(Vector{Float64}(), a, inverse_retraction_method) end end function (r::RecordChange)(amp::AbstractManoptProblem, s::AbstractManoptSolverState, i::Int) @@ -474,19 +476,17 @@ record a certain entries change during iterates * `distance` – function (p,o,x1,x2) to compute the change/distance between two values of the entry * `storage` – a [`StoreStateAction`](@ref) to store (at least) `getproperty(o, d.field)` """ -mutable struct RecordEntryChange <: RecordAction +mutable struct RecordEntryChange{TStorage<:StoreStateAction} <: RecordAction recorded_values::Vector{Float64} field::Symbol distance::Any - storage::StoreStateAction + storage::TStorage function RecordEntryChange(f::Symbol, d, a::StoreStateAction=StoreStateAction([f])) - return new(Float64[], f, d, a) + return new{typeof(a)}(Float64[], f, d, a) end - function RecordEntryChange( - v::T where {T}, f::Symbol, d, a::StoreStateAction=StoreStateAction([f]) - ) + function RecordEntryChange(v, f::Symbol, d, a::StoreStateAction=StoreStateAction([f])) update_storage!(a, Dict(f => v)) - return new(Float64[], f, d, a) + return new{typeof(a)}(Float64[], f, d, a) end end function (r::RecordEntryChange)( diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index 6b12e80374..0b5b86394d 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -183,13 +183,13 @@ struct PointStorageKey{key} end PointStorageKey(key::Symbol) = PointStorageKey{key}() """ - struct TangentStorageKey{key} end + struct VectorStorageKey{key} end Refer to tangent storage of [`StoreStateAction`](@ref) in `get_storage` and `has_storage` functions """ -struct TangentStorageKey{key} end -TangentStorageKey(key::Symbol) = TangentStorageKey{key}() +struct VectorStorageKey{key} end +VectorStorageKey(key::Symbol) = VectorStorageKey{key}() # # Common Actions for decorated AbstractManoptSolverState @@ -238,12 +238,12 @@ iteration, i.e. acts on `(p,o,i)`, where `p` is a [`AbstractManoptProblem`](@ref `point_values` with matching key has been already initialized to a value. When it is false, it corresponds to a general value not being stored for the key present in the vector `keys`. -* `tangent_values` – a `NamedTuple` of mutable values of tangent vectors on a manifold to be +* `vector_values` – a `NamedTuple` of mutable values of tangent vectors on a manifold to be stored in `StoreStateAction`. Manifold is later determined by `AbstractManoptProblem` passed to `update_storage!`. It is not specified at which point the vectors are tangent but for storage it should not matter. * `vector_init` – a `NamedTuple` of boolean values indicating whether a tangent vector in - `tangent_values` with matching key has been already initialized to a value. When it is + `vector_values` with matching key has been already initialized to a value. When it is false, it corresponds to a general value not being stored for the key present in the vector `keys`. * `once` – whether to update the internal values only once per iteration @@ -251,7 +251,7 @@ iteration, i.e. acts on `(p,o,i)`, where `p` is a [`AbstractManoptProblem`](@ref To handle the general storage, use `get_storage` and `has_storage` with keys as `Symbol`s. For the point storage use `PointStorageKey`. For tangent vector storage use -`TangentStorageKey`. Point and tangent storage have been optimized to be more efficient. +`VectorStorageKey`. Point and tangent storage have been optimized to be more efficient. # Constructiors @@ -285,48 +285,46 @@ mutable struct StoreStateAction{ values::Dict{Symbol,Any} keys::Vector{Symbol} # for values point_values::TPS - tangent_values::TXS + vector_values::TXS point_init::TPI - tangent_init::TTI + vector_init::TTI once::Bool last_stored::Int function StoreStateAction( general_keys::Vector{Symbol}=Symbol[], point_values::NamedTuple=NamedTuple(), - tangent_values::NamedTuple=NamedTuple(), + vector_values::NamedTuple=NamedTuple(), once::Bool=true; M::AbstractManifold=DefaultManifold(), ) point_init = NamedTuple{keys(point_values)}(map(u -> false, keys(point_values))) - tangent_init = NamedTuple{keys(tangent_values)}( - map(u -> false, keys(tangent_values)) - ) + vector_init = NamedTuple{keys(vector_values)}(map(u -> false, keys(vector_values))) point_values_copy = NamedTuple{keys(point_values)}( map(u -> _storage_copy_point(M, point_values[u]), keys(point_values)) ) - tangent_values_copy = NamedTuple{keys(tangent_values)}( - map(u -> _storage_copy_vector(M, tangent_values[u]), keys(tangent_values)) + vector_values_copy = NamedTuple{keys(vector_values)}( + map(u -> _storage_copy_vector(M, vector_values[u]), keys(vector_values)) ) return new{ typeof(point_values), - typeof(tangent_values), + typeof(vector_values), typeof(point_values_copy), - typeof(tangent_values_copy), + typeof(vector_values_copy), typeof(point_init), - typeof(tangent_init), + typeof(vector_init), }( Dict{Symbol,Any}(), general_keys, point_values_copy, - tangent_values_copy, + vector_values_copy, point_init, - tangent_init, + vector_init, once, -1, ) end end -@inline function StoreStateAction( +@noinline function StoreStateAction( M::AbstractManifold; store_fields::Vector{Symbol}=Symbol[], store_points::Union{Type{TPS},Vector{Symbol}}=Tuple{}, @@ -346,8 +344,8 @@ end TTS_tuple = Tuple(TTS.parameters) end point_values = NamedTuple{TPS_tuple}(map(_ -> p_init, TPS_tuple)) - tangent_values = NamedTuple{TTS_tuple}(map(_ -> X_init, TTS_tuple)) - return StoreStateAction(store_fields, point_values, tangent_values, once; M=M) + vector_values = NamedTuple{TTS_tuple}(map(_ -> X_init, TTS_tuple)) + return StoreStateAction(store_fields, point_values, vector_values, once; M=M) end @generated function extract_type_from_namedtuple(::Type{nt}, ::Val{key}) where {nt,key} @@ -365,7 +363,7 @@ function _store_point_assert_type( return extract_type_from_namedtuple(TPS_asserts, key) end -function _store_tangent_assert_type( +function _store_vector_assert_type( ::StoreStateAction{TPS_asserts,TXS_assert}, key::Val ) where {TPS_asserts,TXS_assert} return extract_type_from_namedtuple(TXS_assert, key) @@ -404,14 +402,14 @@ function get_storage(a::AbstractStateAction, ::PointStorageKey{key}) where {key} end """ - get_storage(a::AbstractStateAction, ::TangentStorageKey{key}) where {key} + get_storage(a::AbstractStateAction, ::VectorStorageKey{key}) where {key} Return the internal value of the [`AbstractStateAction`](@ref) `a` at the -`Symbol` `key` that represents a tangent vector. +`Symbol` `key` that represents a vector vector. """ -function get_storage(a::AbstractStateAction, ::TangentStorageKey{key}) where {key} - if haskey(a.tangent_values, key) - return a.tangent_values[key] +function get_storage(a::AbstractStateAction, ::VectorStorageKey{key}) where {key} + if haskey(a.vector_values, key) + return a.vector_values[key] else return get_storage(a, key) end @@ -440,14 +438,14 @@ function has_storage(a::AbstractStateAction, ::PointStorageKey{key}) where {key} end """ - has_storage(a::AbstractStateAction, ::TangentStorageKey{key}) where {key} + has_storage(a::AbstractStateAction, ::VectorStorageKey{key}) where {key} Return whether the [`AbstractStateAction`](@ref) `a` has a point value stored at the `Symbol` `key`. """ -function has_storage(a::AbstractStateAction, ::TangentStorageKey{key}) where {key} - if haskey(a.tangent_init, key) - return a.tangent_init[key] +function has_storage(a::AbstractStateAction, ::VectorStorageKey{key}) where {key} + if haskey(a.vector_init, key) + return a.vector_init[key] else return has_storage(a, key) end @@ -488,22 +486,20 @@ function update_storage!( map(update_points, keys(a.point_values)) a.point_init = NamedTuple{keys(a.point_values)}(map(u -> true, keys(a.point_values))) - @inline function update_tangent(key) + @inline function update_vector(key) if key === :Gradient - copyto!(M, a.tangent_values[key], get_gradient(s)) + copyto!(M, a.vector_values[key], get_gradient(s)) else copyto!( M, - a.tangent_values[key], - getproperty(s, key)::_store_tangent_assert_type(a, Val(key)), + a.vector_values[key], + getproperty(s, key)::_store_vector_assert_type(a, Val(key)), ) end end - map(update_tangent, keys(a.tangent_values)) + map(update_vector, keys(a.vector_values)) - a.tangent_init = NamedTuple{keys(a.tangent_values)}( - map(u -> true, keys(a.tangent_values)) - ) + a.vector_init = NamedTuple{keys(a.vector_values)}(map(u -> true, keys(a.vector_values))) return a.keys end diff --git a/src/plans/stepsize.jl b/src/plans/stepsize.jl index 7443934433..ec6f0f8532 100644 --- a/src/plans/stepsize.jl +++ b/src/plans/stepsize.jl @@ -456,7 +456,7 @@ mutable struct NonmonotoneLinesearch{ throw(DomainError(memory_size, "The memory_size has to be greater than zero.")) end if isnothing(storage) - if (M isa DefaultManifold) || (rand(M) isa Number) + if M isa DefaultManifold storage = StoreStateAction(M; store_fields=[:Iterate, :Gradient]) else storage = StoreStateAction( diff --git a/test/plans/test_storage.jl b/test/plans/test_storage.jl index 1211a6a178..6f6125bf51 100644 --- a/test/plans/test_storage.jl +++ b/test/plans/test_storage.jl @@ -24,24 +24,24 @@ using Test, Manopt, ManifoldsBase, Manifolds a = StoreStateAction(M; store_fields=[:p, :X]) @test !has_storage(a, Manopt.PointStorageKey(:p)) - @test !has_storage(a, Manopt.TangentStorageKey(:X)) + @test !has_storage(a, Manopt.VectorStorageKey(:X)) update_storage!(a, mp, st) @test has_storage(a, Manopt.PointStorageKey(:p)) - @test has_storage(a, Manopt.TangentStorageKey(:X)) + @test has_storage(a, Manopt.VectorStorageKey(:X)) @test get_storage(a, Manopt.PointStorageKey(:p)) == p - @test get_storage(a, Manopt.TangentStorageKey(:X)) == X_zero + @test get_storage(a, Manopt.VectorStorageKey(:X)) == X_zero a2 = StoreStateAction(M; store_points=Tuple{:p}, store_vectors=Tuple{:X}) @test !has_storage(a2, Manopt.PointStorageKey(:p)) - @test !has_storage(a2, Manopt.TangentStorageKey(:X)) + @test !has_storage(a2, Manopt.VectorStorageKey(:X)) update_storage!(a2, mp, st) @test has_storage(a2, Manopt.PointStorageKey(:p)) - @test has_storage(a2, Manopt.TangentStorageKey(:X)) + @test has_storage(a2, Manopt.VectorStorageKey(:X)) @test get_storage(a2, Manopt.PointStorageKey(:p)) == p_fast - @test get_storage(a2, Manopt.TangentStorageKey(:X)) == X_zero_fast + @test get_storage(a2, Manopt.VectorStorageKey(:X)) == X_zero_fast a2b = StoreStateAction(M; store_points=Tuple{:p}, store_vectors=Tuple{:X}) @test keys(a2.point_values) == keys(a2b.point_values) - @test keys(a2.tangent_values) == keys(a2b.tangent_values) + @test keys(a2.vector_values) == keys(a2b.vector_values) @test keys(a2.keys) == keys(a2b.keys) # make sure fast storage is actually fast @@ -51,12 +51,12 @@ using Test, Manopt, ManifoldsBase, Manifolds a3 = StoreStateAction(M; store_points=[:p], store_vectors=[:X]) @test !has_storage(a3, Manopt.PointStorageKey(:p)) - @test !has_storage(a3, Manopt.TangentStorageKey(:X)) + @test !has_storage(a3, Manopt.VectorStorageKey(:X)) update_storage!(a3, mp, st) @test has_storage(a3, Manopt.PointStorageKey(:p)) - @test has_storage(a3, Manopt.TangentStorageKey(:X)) + @test has_storage(a3, Manopt.VectorStorageKey(:X)) @test get_storage(a3, Manopt.PointStorageKey(:p)) == p_fast - @test get_storage(a3, Manopt.TangentStorageKey(:X)) == X_zero_fast + @test get_storage(a3, Manopt.VectorStorageKey(:X)) == X_zero_fast # make sure fast storage is actually fast @test (@allocated update_storage!(a3, mp, st)) == 0 From ede0d6cd8240f503eb6ea9ce7e3e6195631f5a20 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Mon, 20 Feb 2023 11:37:37 +0100 Subject: [PATCH 32/37] fix nonmutating storage --- src/plans/solver_state.jl | 31 +++++++++++++++++++++++++------ src/plans/stepsize.jl | 7 ++++--- test/plans/test_storage.jl | 20 ++++++++++---------- 3 files changed, 39 insertions(+), 19 deletions(-) diff --git a/src/plans/solver_state.jl b/src/plans/solver_state.jl index 0b5b86394d..9a25a0b4e5 100644 --- a/src/plans/solver_state.jl +++ b/src/plans/solver_state.jl @@ -202,13 +202,22 @@ for example within the [`DebugSolverState`](@ref) or within the [`RecordSolverSt """ abstract type AbstractStateAction end +mutable struct StorageRef{T} + x::T +end + +function Base.copyto!(sr::StorageRef, new_x) + sr.x = copy(new_x) + return sr +end + """ _storage_copy_point(M::AbstractManifold, p) Make a copy of point `p` from manifold `M` for storage in [`StoreStateAction`](@ref). """ _storage_copy_point(M::AbstractManifold, p) = copy(M, p) -_storage_copy_point(::AbstractManifold, p::Number) = fill(p) +_storage_copy_point(::AbstractManifold, p::Number) = StorageRef(p) """ _storage_copy_vector(M::AbstractManifold, X) @@ -216,7 +225,7 @@ _storage_copy_point(::AbstractManifold, p::Number) = fill(p) Make a copy of tangent vector `X` from manifold `M` for storage in [`StoreStateAction`](@ref). """ _storage_copy_vector(M::AbstractManifold, X) = copy(M, SA_F64[], X) -_storage_copy_vector(::AbstractManifold, X::Number) = fill(X) +_storage_copy_vector(::AbstractManifold, X::Number) = StorageRef(X) @doc raw""" StoreStateAction <: AbstractStateAction @@ -393,9 +402,14 @@ get_storage(a::AbstractStateAction, key::Symbol) = a.values[key] Return the internal value of the [`AbstractStateAction`](@ref) `a` at the `Symbol` `key` that represents a point. """ -function get_storage(a::AbstractStateAction, ::PointStorageKey{key}) where {key} +@inline function get_storage(a::AbstractStateAction, ::PointStorageKey{key}) where {key} if haskey(a.point_values, key) - return a.point_values[key] + val = a.point_values[key] + if val isa StorageRef + return val.x + else + return val + end else return get_storage(a, key) end @@ -407,9 +421,14 @@ end Return the internal value of the [`AbstractStateAction`](@ref) `a` at the `Symbol` `key` that represents a vector vector. """ -function get_storage(a::AbstractStateAction, ::VectorStorageKey{key}) where {key} +@inline function get_storage(a::AbstractStateAction, ::VectorStorageKey{key}) where {key} if haskey(a.vector_values, key) - return a.vector_values[key] + val = a.vector_values[key] + if val isa StorageRef + return val.x + else + return val + end else return get_storage(a, key) end diff --git a/src/plans/stepsize.jl b/src/plans/stepsize.jl index ec6f0f8532..1db2718a32 100644 --- a/src/plans/stepsize.jl +++ b/src/plans/stepsize.jl @@ -491,13 +491,14 @@ function (a::NonmonotoneLinesearch)( η=-get_gradient(mp, get_iterate(s)); kwargs..., ) - if !has_storage(a.storage, :Iterate) || !has_storage(a.storage, :Gradient) + if !has_storage(a.storage, PointStorageKey(:Iterate)) || + !has_storage(a.storage, VectorStorageKey(:Gradient)) p_old = get_iterate(s) X_old = get_gradient(mp, p_old) else #fetch - p_old = get_storage(a.storage, :Iterate) - X_old = get_storage(a.storage, :Gradient) + p_old = get_storage(a.storage, PointStorageKey(:Iterate)) + X_old = get_storage(a.storage, VectorStorageKey(:Gradient)) end update_storage!(a.storage, mp, s) return a( diff --git a/test/plans/test_storage.jl b/test/plans/test_storage.jl index 6f6125bf51..ebc495757b 100644 --- a/test/plans/test_storage.jl +++ b/test/plans/test_storage.jl @@ -4,14 +4,10 @@ using Test, Manopt, ManifoldsBase, Manifolds @testset "manifold $M" for M in [ManifoldsBase.DefaultManifold(2), Circle()] if M isa Circle p = 0.4 - p_fast = fill(p) X_zero = 0.0 - X_zero_fast = fill(0.0) else p = [4.0, 2.0] - p_fast = p X_zero = [0.0, 0.0] - X_zero_fast = X_zero end st = GradientDescentState( @@ -37,8 +33,8 @@ using Test, Manopt, ManifoldsBase, Manifolds update_storage!(a2, mp, st) @test has_storage(a2, Manopt.PointStorageKey(:p)) @test has_storage(a2, Manopt.VectorStorageKey(:X)) - @test get_storage(a2, Manopt.PointStorageKey(:p)) == p_fast - @test get_storage(a2, Manopt.VectorStorageKey(:X)) == X_zero_fast + @test get_storage(a2, Manopt.PointStorageKey(:p)) == p + @test get_storage(a2, Manopt.VectorStorageKey(:X)) == X_zero a2b = StoreStateAction(M; store_points=Tuple{:p}, store_vectors=Tuple{:X}) @test keys(a2.point_values) == keys(a2b.point_values) @test keys(a2.vector_values) == keys(a2b.vector_values) @@ -47,7 +43,9 @@ using Test, Manopt, ManifoldsBase, Manifolds # make sure fast storage is actually fast @test (@allocated update_storage!(a2, mp, st)) == 0 @test (@allocated has_storage(a2, Manopt.PointStorageKey(:p))) == 0 - @test (@allocated get_storage(a2, Manopt.PointStorageKey(:p))) == 0 + if M isa ManifoldsBase.DefaultManifold + @test (@allocated get_storage(a2, Manopt.PointStorageKey(:p))) == 0 + end a3 = StoreStateAction(M; store_points=[:p], store_vectors=[:X]) @test !has_storage(a3, Manopt.PointStorageKey(:p)) @@ -55,13 +53,15 @@ using Test, Manopt, ManifoldsBase, Manifolds update_storage!(a3, mp, st) @test has_storage(a3, Manopt.PointStorageKey(:p)) @test has_storage(a3, Manopt.VectorStorageKey(:X)) - @test get_storage(a3, Manopt.PointStorageKey(:p)) == p_fast - @test get_storage(a3, Manopt.VectorStorageKey(:X)) == X_zero_fast + @test get_storage(a3, Manopt.PointStorageKey(:p)) == p + @test get_storage(a3, Manopt.VectorStorageKey(:X)) == X_zero # make sure fast storage is actually fast @test (@allocated update_storage!(a3, mp, st)) == 0 @test (@allocated has_storage(a3, Manopt.PointStorageKey(:p))) == 0 - @test (@allocated get_storage(a3, Manopt.PointStorageKey(:p))) == 0 + if M isa ManifoldsBase.DefaultManifold + @test (@allocated get_storage(a3, Manopt.PointStorageKey(:p))) == 0 + end end @test Manopt.extract_type_from_namedtuple(typeof((; a=10, b='a')), Val(:c)) === Any From 60623873d5fbf45473f313162d00c1f7c47d1f50 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Mon, 20 Feb 2023 11:51:39 +0100 Subject: [PATCH 33/37] use more storage keys --- src/plans/debug.jl | 4 ++-- src/plans/record.jl | 4 ++-- src/plans/stopping_criterion.jl | 4 ++-- src/solvers/particle_swarm.jl | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/plans/debug.jl b/src/plans/debug.jl index 70179e2b3b..dcf0bf56cd 100644 --- a/src/plans/debug.jl +++ b/src/plans/debug.jl @@ -436,8 +436,8 @@ function (d::DebugGradientChange)( ) if i > 0 M = get_manifold(pm) - p_old = get_storage(d.storage, :Iterate) - X_old = get_storage(d.storage, :Gradient) + p_old = get_storage(d.storage, PointStorageKey(:Iterate)) + X_old = get_storage(d.storage, VectorStorageKey(:Gradient)) p = get_iterate(st) X = get_gradient(st) l = norm( diff --git a/src/plans/record.jl b/src/plans/record.jl index e96ab3ff8c..98fa111d59 100644 --- a/src/plans/record.jl +++ b/src/plans/record.jl @@ -417,11 +417,11 @@ function (r::RecordChange)(amp::AbstractManoptProblem, s::AbstractManoptSolverSt M = get_manifold(amp) record_or_reset!( r, - if has_storage(r.storage, :Iterate) + if has_storage(r.storage, PointStorageKey(:Iterate)) distance( M, get_iterate(s), - get_storage(r.storage, :Iterate), + get_storage(r.storage, PointStorageKey(:Iterate)), r.inverse_retraction_method, ) else diff --git a/src/plans/stopping_criterion.jl b/src/plans/stopping_criterion.jl index 727f2794f1..83571500dd 100644 --- a/src/plans/stopping_criterion.jl +++ b/src/plans/stopping_criterion.jl @@ -203,9 +203,9 @@ function (c::StopWhenChangeLess)(mp::AbstractManoptProblem, s::AbstractManoptSol c.reason = "" c.at_iteration = 0 end - if has_storage(c.storage, :Iterate) + if has_storage(c.storage, PointStorageKey(:Iterate)) M = get_manifold(mp) - x_old = get_storage(c.storage, :Iterate) + x_old = get_storage(c.storage, PointStorageKey(:Iterate)) d = distance(M, get_iterate(s), x_old, c.inverse_retraction) if d < c.threshold && i > 0 c.reason = "The algorithm performed a step with a change ($d) less than $(c.threshold).\n" diff --git a/src/solvers/particle_swarm.jl b/src/solvers/particle_swarm.jl index df861c36df..895eca272a 100644 --- a/src/solvers/particle_swarm.jl +++ b/src/solvers/particle_swarm.jl @@ -298,8 +298,8 @@ get_solver_result(s::ParticleSwarmState) = s.g # but also lives in the power manifold on M, so we have to adapt StopWhenChangeless # function (c::StopWhenChangeLess)(mp::AbstractManoptProblem, s::ParticleSwarmState, i) - if has_storage(c.storage, :Iterate) - x_old = get_storage(c.storage, :Iterate) + if has_storage(c.storage, PointStorageKey(:Iterate)) + x_old = get_storage(c.storage, PointStorageKey(:Iterate)) n = length(s.x) d = distance( PowerManifold(get_manifold(mp), NestedPowerRepresentation(), n), s.x, x_old From 5f3e98db42294f666007a8f5f597c1674208e528 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Mon, 20 Feb 2023 19:25:15 +0100 Subject: [PATCH 34/37] add a few stati I missed in the last PR (Sorrey). --- src/plans/gradient_plan.jl | 15 ++++++++++++++- src/plans/record.jl | 2 +- test/plans/test_debug.jl | 11 +++++++++++ test/plans/test_record.jl | 5 ++++- 4 files changed, 30 insertions(+), 3 deletions(-) diff --git a/src/plans/gradient_plan.jl b/src/plans/gradient_plan.jl index 591ae135ee..25ffbb3bfd 100644 --- a/src/plans/gradient_plan.jl +++ b/src/plans/gradient_plan.jl @@ -470,6 +470,10 @@ function (d::DebugGradient)(::AbstractManoptProblem, s::AbstractManoptSolverStat Printf.format(d.io, Printf.Format(d.format), get_gradient(s)) return nothing end +function show(io::IO, dg::DebugGradient) + return print(io, "DebugGradient(; format=\"$(dg.format)\")") +end +status_summary(dg::DebugGradient) = "(:Gradient, \"$(dg.format)\")" @doc raw""" DebugGradientNorm <: DebugAction @@ -508,6 +512,10 @@ function (d::DebugGradientNorm)( ) return nothing end +function show(io::IO, dgn::DebugGradientNorm) + return print(io, "DebugGradientNorm(; format=\"$(dgn.format)\")") +end +status_summary(dgn::DebugGradientNorm) = "(:GradientNorm, \"$(dgn.format)\")" @doc raw""" DebugStepsize <: DebugAction @@ -538,7 +546,10 @@ function (d::DebugStepsize)( Printf.format(d.io, Printf.Format(d.format), get_last_stepsize(p, s, i)) return nothing end - +function show(io::IO, ds::DebugStepsize) + return print(io, "DebugStepsize(; format=\"$(ds.format)\")") +end +status_summary(ds::DebugStepsize) = "(:Stepsize, \"$(ds.format)\")" # # Records # @@ -562,6 +573,7 @@ function (r::RecordGradient{T})( ) where {T} return record_or_reset!(r, get_gradient(s), i) end +show(io::IO, ::RecordGradient{T}) where {T} = print(io, "RecordGradient{$T}()") @doc raw""" RecordGradientNorm <: RecordAction @@ -578,6 +590,7 @@ function (r::RecordGradientNorm)( M = get_manifold(mp) return record_or_reset!(r, norm(M, get_iterate(ast), get_gradient(ast)), i) end +show(io::IO, ::RecordGradientNorm) = print(io, "RecordGradientNorm()") @doc raw""" RecordStepsize <: RecordAction diff --git a/src/plans/record.jl b/src/plans/record.jl index 98fa111d59..213b644dbd 100644 --- a/src/plans/record.jl +++ b/src/plans/record.jl @@ -434,7 +434,7 @@ function (r::RecordChange)(amp::AbstractManoptProblem, s::AbstractManoptSolverSt end function show(io::IO, rc::RecordChange) return print( - io, "DebugChange(; inverse_retraction_method=$(rc.inverse_retraction_method))" + io, "RecordChange(; inverse_retraction_method=$(rc.inverse_retraction_method))" ) end status_summary(rc::RecordChange) = ":Change" diff --git a/test/plans/test_debug.jl b/test/plans/test_debug.jl index e930806587..e7d0cfa014 100644 --- a/test/plans/test_debug.jl +++ b/test/plans/test_debug.jl @@ -326,5 +326,16 @@ struct TestDebugAction <: DebugAction end @test Manopt.status_summary(DebugDivider("a")) == "\"a\"" @test repr(DebugEntry(:a)) == "DebugEntry(:a; format=\"a: %s\")" + + @test repr(DebugStepsize()) == "DebugStepsize(; format=\"s:%s\")" + @test Manopt.status_summary(DebugStepsize()) == "(:Stepsize, \"s:%s\")" + + @test repr(DebugGradientNorm()) == "DebugGradientNorm(; format=\"|grad f(p)|:%s\")" + dgn_s = "(:GradientNorm, \"|grad f(p)|:%s\")" + @test Manopt.status_summary(DebugGradientNorm()) == dgn_s + + @test repr(DebugGradient()) == "DebugGradient(; format=\"grad f(p):%s\")" + dg_s = "(:Gradient, \"grad f(p):%s\")" + @test Manopt.status_summary(DebugGradient()) == dg_s end end diff --git a/test/plans/test_record.jl b/test/plans/test_record.jl index 3fc69421c1..8440699a46 100644 --- a/test/plans/test_record.jl +++ b/test/plans/test_record.jl @@ -110,7 +110,7 @@ using Manifolds, Manopt, Test, ManifoldsBase, Dates @test c2[:It1] == [10, 20] # RecordChange d = RecordChange() - sd = "DebugChange(; inverse_retraction_method=LogarithmicInverseRetraction())" + sd = "RecordChange(; inverse_retraction_method=LogarithmicInverseRetraction())" @test repr(d) == sd @test Manopt.status_summary(d) == ":Change" d(dmp, gds, 1) @@ -211,4 +211,7 @@ using Manifolds, Manopt, Test, ManifoldsBase, Dates # stop after 20 so 21 hits h3(dmp, gds, 20) @test length(h3.recorded_values) == 1 + @test repr(RecordGradientNorm()) == "RecordGradientNorm()" + # since only the type is stored we get + @test repr(RecordGradient(zeros(3))) == "RecordGradient{Vector{Float64}}()" end From eba6978330a5aafd2c584c02c172e19aadddade6 Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Tue, 21 Feb 2023 17:02:07 +0100 Subject: [PATCH 35/37] a few small quasi Newton optimizations --- src/Manopt.jl | 2 ++ src/plans/quasi_newton_plan.jl | 9 +++++---- src/solvers/quasi_Newton.jl | 14 ++++++++------ 3 files changed, 15 insertions(+), 10 deletions(-) diff --git a/src/Manopt.jl b/src/Manopt.jl index 349a0c1260..af3f8d64ad 100644 --- a/src/Manopt.jl +++ b/src/Manopt.jl @@ -77,6 +77,8 @@ using ManifoldsBase: default_retraction_method, default_vector_transport_method, distance, + embed_project, + embed_project!, exp, exp!, geodesic, diff --git a/src/plans/quasi_newton_plan.jl b/src/plans/quasi_newton_plan.jl index e6622c3c50..dde7fd8341 100644 --- a/src/plans/quasi_newton_plan.jl +++ b/src/plans/quasi_newton_plan.jl @@ -478,14 +478,15 @@ function (d::QuasiNewtonLimitedMemoryDirectionUpdate{InverseBFGS})(mp, st) for i in m:-1:1 d.ρ[i] = 1 / inner(M, p, d.memory_s[i], d.memory_y[i]) # 1 sk 2 yk d.ξ[i] = inner(M, p, d.memory_s[i], r) * d.ρ[i] - r .= r .- d.ξ[i] .* d.memory_y[i] + r .-= d.ξ[i] .* d.memory_y[i] end r .= 1 / (d.ρ[m] * norm(M, p, last(d.memory_y))^2) .* r for i in 1:m - r .= r .+ (d.ξ[i] - d.ρ[i] * inner(M, p, d.memory_y[i], r)) .* d.memory_s[i] + r .+= (d.ξ[i] - d.ρ[i] * inner(M, p, d.memory_y[i], r)) .* d.memory_s[i] end - d.project && project!(M, r, p, r) - return -r + d.project && embed_project!(M, r, p, r) + r .*= -1 + return r end @doc raw""" diff --git a/src/solvers/quasi_Newton.jl b/src/solvers/quasi_Newton.jl index 2dd73d5430..fed3cd1a43 100644 --- a/src/solvers/quasi_Newton.jl +++ b/src/solvers/quasi_Newton.jl @@ -49,6 +49,7 @@ mutable struct QuasiNewtonState{ retraction_method::RTR stepsize::S stop::SC + p_old::P vector_transport_method::VT end function QuasiNewtonState( @@ -86,6 +87,7 @@ function QuasiNewtonState( retraction_method, stepsize, stopping_criterion, + copy(M, p), vector_transport_method, ) end @@ -287,21 +289,21 @@ function step_solver!(mp::AbstractManoptProblem, qns::QuasiNewtonState, iter) qns.X = get_gradient(mp, qns.p) η = qns.direction_update(mp, qns) α = qns.stepsize(mp, qns, iter, η) - p_old = copy(M, get_iterate(qns)) + copyto!(M, qns.p_old, get_iterate(qns)) αη = α * η retract!(M, qns.p, qns.p, η, α, qns.retraction_method) β = locking_condition_scale( - M, qns.direction_update, p_old, αη, qns.p, qns.vector_transport_method + M, qns.direction_update, qns.p_old, αη, qns.p, qns.vector_transport_method ) vector_transport_to!( - M, qns.sk, p_old, αη, qns.p, get_update_vector_transport(qns.direction_update) + M, qns.sk, qns.p_old, αη, qns.p, get_update_vector_transport(qns.direction_update) ) vector_transport_to!( - M, qns.X, p_old, qns.X, qns.p, get_update_vector_transport(qns.direction_update) + M, qns.X, qns.p_old, qns.X, qns.p, get_update_vector_transport(qns.direction_update) ) new_grad = get_gradient(mp, qns.p) - qns.yk = new_grad ./ β .- qns.X - update_hessian!(qns.direction_update, mp, qns, p_old, iter) + qns.yk .= new_grad ./ β .- qns.X + update_hessian!(qns.direction_update, mp, qns, qns.p_old, iter) return qns end From a0c5bbbe8ceb8383a4fc41ec4ad21c4fcf6b86de Mon Sep 17 00:00:00 2001 From: Mateusz Baran Date: Tue, 21 Feb 2023 17:43:35 +0100 Subject: [PATCH 36/37] revert p_old storage in quasi Newton --- src/solvers/quasi_Newton.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/solvers/quasi_Newton.jl b/src/solvers/quasi_Newton.jl index fed3cd1a43..f986cd2e0a 100644 --- a/src/solvers/quasi_Newton.jl +++ b/src/solvers/quasi_Newton.jl @@ -49,7 +49,6 @@ mutable struct QuasiNewtonState{ retraction_method::RTR stepsize::S stop::SC - p_old::P vector_transport_method::VT end function QuasiNewtonState( @@ -87,7 +86,6 @@ function QuasiNewtonState( retraction_method, stepsize, stopping_criterion, - copy(M, p), vector_transport_method, ) end @@ -289,21 +287,21 @@ function step_solver!(mp::AbstractManoptProblem, qns::QuasiNewtonState, iter) qns.X = get_gradient(mp, qns.p) η = qns.direction_update(mp, qns) α = qns.stepsize(mp, qns, iter, η) - copyto!(M, qns.p_old, get_iterate(qns)) + p_old = copy(M, get_iterate(qns)) αη = α * η retract!(M, qns.p, qns.p, η, α, qns.retraction_method) β = locking_condition_scale( - M, qns.direction_update, qns.p_old, αη, qns.p, qns.vector_transport_method + M, qns.direction_update, p_old, αη, qns.p, qns.vector_transport_method ) vector_transport_to!( - M, qns.sk, qns.p_old, αη, qns.p, get_update_vector_transport(qns.direction_update) + M, qns.sk, p_old, αη, qns.p, get_update_vector_transport(qns.direction_update) ) vector_transport_to!( - M, qns.X, qns.p_old, qns.X, qns.p, get_update_vector_transport(qns.direction_update) + M, qns.X, p_old, qns.X, qns.p, get_update_vector_transport(qns.direction_update) ) new_grad = get_gradient(mp, qns.p) qns.yk .= new_grad ./ β .- qns.X - update_hessian!(qns.direction_update, mp, qns, qns.p_old, iter) + update_hessian!(qns.direction_update, mp, qns, p_old, iter) return qns end From f4172a09b162730198a056b7e1abea472e66091d Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Tue, 21 Feb 2023 18:33:03 +0100 Subject: [PATCH 37/37] Finalize 0.4.8 --- Changelog.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Changelog.md b/Changelog.md index e7eb034817..177bfe6c30 100644 --- a/Changelog.md +++ b/Changelog.md @@ -5,7 +5,7 @@ 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.8] +## [0.4.8] - 21/02/2023 ### Added