JuliaManifolds / Manopt.jl

🏔️Manopt. jl – Optimization on Manifolds in Julia
http://manoptjl.org
Other
314 stars 40 forks source link

Failed calling of the retraction of submanifolds in PowerManifold #330

Closed yuwenchen95 closed 9 months ago

yuwenchen95 commented 9 months ago

I tried to optimize a quadratic function over a PowerManifold whose submanifolds are a cluster of Stiefel manifolds.

###############################################################
# Riemmanian manifold algorithms
###############################################################
using Manifolds, ManifoldsBase, Manopt
using SparseArrays, LinearAlgebra

n = 100
d = 3
dim = 300
C = sprand(300,300,0.1)
dim = d*n

rnk = Int(ceil(sqrt(2*dim)))
M = PowerManifold(Stiefel(rnk,d),n) 
function F(M::AbstractPowerManifold,x)
    xred = reshape(x,rnk,d*n)
    return tr(xred*C*xred')
end

###############################################################
# in-place gradient
###############################################################
struct grad!{TD,TTMP}
    C::TD
    tmp::TTMP
end

function (gradF!::grad!)(M,y,x)
    xred = reshape(x,rnk,d*n)
    yred = reshape(y,rnk,d*n)
    mul!(yred,xred,gradF!.C,2.,0.)
    Threads.@threads for i = 1:n
        id = Threads.threadid()

        cur_x = @view x[:,:,i]
        cur_y = @view y[:,:,i]

        dotxy = @view gradF!.tmp[:,:,id]
        mul!(dotxy,cur_x',cur_y)
        @. dotxy = (dotxy + dotxy')/2.
        cur_y .-= cur_x*dotxy
    end

    return y
end

gradF2! = grad!(C,zeros(d,d,Threads.nthreads()))

###############################################################################
# 1st-order method
###############################################################################
# start to solve
x0 = rand(rnk,d,n)
for i = 1:n
    cur_x = @view x0[:,:,i]
    (Ui,Si,VTi) = LAPACK.gesdd!('S',cur_x)
    mul!(cur_x,Ui,VTi)
end

ManifoldsBase.default_retraction_method(M::Stiefel) = PolarRetraction()
ManifoldsBase.default_retraction_method(M::PowerManifold) = ManifoldsBase.default_retraction_method(M.manifold)
ManifoldsBase.default_inverse_retraction_method(M::Stiefel) = PolarInverseRetraction()
ManifoldsBase.default_inverse_retraction_method(M::PowerManifold) = ManifoldsBase.default_inverse_retraction_method(M.manifold)

dg=[:Iteration, " | ", :IterativeTime, " | ", :Cost, "\n", :Stop]
rd=[:p, :IterativeTime, :Cost]
ro=true
st = ArmijoLinesearch(M)
sc = StopAfterIteration(50) | StopWhenGradientNormLess(1e-10)

m = deepcopy(x0)
rgd = gradient_descent!(
    M, F, gradF2!, m; evaluation=InplaceEvaluation(), stepsize=st, stopping_criterion=sc,
    record = rd, return_state=true#, debug= dg
)

However, it returns the following error:

julia> include("test_powermanifolds.jl")
ERROR: LoadError: MethodError: no method matching retract_polar!(::PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ::Array{Float64, 3}, ::Array{Float64, 3}, ::Array{Float64, 3})

Closest candidates are:
  retract_polar!(::AbstractManifold, ::Any, ::Any, ::Any, ::Any)
   @ ManifoldsBase C:\Users\ddt00\.julia\packages\ManifoldsBase\RQVnN\src\retractions.jl:904
  retract_polar!(::Flag, ::Any, ::Any, ::Any, ::Number)
   @ Manifolds C:\Users\ddt00\.julia\packages\Manifolds\FajGI\src\manifolds\FlagStiefel.jl:185
  retract_polar!(::Stiefel, ::Any, ::Any, ::Any, ::Number)
   @ Manifolds C:\Users\ddt00\.julia\packages\Manifolds\FajGI\src\manifolds\Stiefel.jl:496
  ...

Stacktrace:
  [1] retract_polar!(M::PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, q::Array{Float64, 3}, p::Array{Float64, 3}, X::Array{Float64, 3}, t::Float64)
    @ ManifoldsBase C:\Users\ddt00\.julia\packages\ManifoldsBase\RQVnN\src\retractions.jl:905
  [2] _retract!(M::PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, q::Array{Float64, 3}, p::Array{Float64, 3}, X::Array{Float64, 3}, t::Float64, ::PolarRetraction; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ManifoldsBase C:\Users\ddt00\.julia\packages\ManifoldsBase\RQVnN\src\retractions.jl:796
  [3] _retract!(M::PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, q::Array{Float64, 3}, p::Array{Float64, 3}, X::Array{Float64, 3}, t::Float64, ::PolarRetraction)
    @ ManifoldsBase C:\Users\ddt00\.julia\packages\ManifoldsBase\RQVnN\src\retractions.jl:795
  [4] retract!(M::PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, q::Array{Float64, 3}, p::Array{Float64, 3}, X::Array{Float64, 3}, t::Float64, m::PolarRetraction)
    @ ManifoldsBase C:\Users\ddt00\.julia\packages\ManifoldsBase\RQVnN\src\retractions.jl:771
  [5] _retract(M::PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, p::Array{Float64, 3}, X::Array{Float64, 3}, t::Float64, m::PolarRetraction)
    @ ManifoldsBase C:\Users\ddt00\.julia\packages\ManifoldsBase\RQVnN\src\retractions.jl:744
  [6] retract(M::PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, p::Array{Float64, 3}, X::Array{Float64, 3}, t::Float64, m::PolarRetraction)
    @ ManifoldsBase C:\Users\ddt00\.julia\packages\ManifoldsBase\RQVnN\src\retractions.jl:736
  [7] linesearch_backtrack(M::PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, f::Manopt.var"#104#105"{DefaultManoptProblem{PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}}}, p::Array{Float64, 3}, grad_f_at_p::Array{Float64, 3}, s::Float64, decrease::Float64, contract::Float64, retr::PolarRetraction, η::Array{Float64, 3}, f0::Float64; stop_when_stepsize_less::Float64, stop_when_stepsize_exceeds::Float64, stop_increasing_at_step::Int64, stop_decreasing_at_step::Int64)
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\plans\stepsize.jl:368
  [8] (::ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)})(mp::DefaultManoptProblem{PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}}, s::GradientDescentState{Array{Float64, 3}, Array{Float64, 3}, StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, IdentityUpdateRule, PolarRetraction}, i::Int64, η::Array{Float64, 3}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\plans\stepsize.jl:288
  [9] (::ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)})(mp::DefaultManoptProblem{PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}}, s::GradientDescentState{Array{Float64, 3}, Array{Float64, 3}, StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, IdentityUpdateRule, PolarRetraction}, i::Int64, η::Array{Float64, 3})
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\plans\stepsize.jl:279
 [10] (::ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)})(mp::DefaultManoptProblem{PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}}, s::GradientDescentState{Array{Float64, 3}, Array{Float64, 3}, StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, IdentityUpdateRule, PolarRetraction}, i::Int64)
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\plans\stepsize.jl:279
 [11] _get_stepsize(amp::DefaultManoptProblem{PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}}, ams::GradientDescentState{Array{Float64, 3}, Array{Float64, 3}, StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, IdentityUpdateRule, PolarRetraction}, ::Val{false}, vars::Int64; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\plans\stepsize.jl:1170
 [12] _get_stepsize(amp::DefaultManoptProblem{PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}}, ams::GradientDescentState{Array{Float64, 3}, Array{Float64, 3}, StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, IdentityUpdateRule, PolarRetraction}, ::Val{false}, vars::Int64)
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\plans\stepsize.jl:1163
 [13] get_stepsize(amp::DefaultManoptProblem{PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}}, ams::GradientDescentState{Array{Float64, 3}, Array{Float64, 3}, StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, IdentityUpdateRule, PolarRetraction}, vars::Int64; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\plans\stepsize.jl:1152
 [14] get_stepsize(amp::DefaultManoptProblem{PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}}, ams::GradientDescentState{Array{Float64, 3}, Array{Float64, 3}, StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, IdentityUpdateRule, PolarRetraction}, vars::Int64)
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\plans\stepsize.jl:1149
 [15] (::IdentityUpdateRule)(mp::DefaultManoptProblem{PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}}, s::GradientDescentState{Array{Float64, 3}, Array{Float64, 3}, StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, IdentityUpdateRule, PolarRetraction}, i::Int64)
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\solvers\gradient_descent.jl:81
 [16] step_solver!(p::DefaultManoptProblem{PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}}, s::GradientDescentState{Array{Float64, 3}, Array{Float64, 3}, StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, IdentityUpdateRule, PolarRetraction}, i::Int64)
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\solvers\gradient_descent.jl:267
 [17] step_solver!(amp::DefaultManoptProblem{PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}}, rss::RecordSolverState{GradientDescentState{Array{Float64, 3}, Array{Float64, 3}, StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, IdentityUpdateRule, PolarRetraction}, NamedTuple{(:Iteration,), Tuple{RecordGroup}}}, i::Int64)
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\solvers\record_solver.jl:20
 [18] step_solver!(p::DefaultManoptProblem{PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}}, s::Manopt.ReturnSolverState{RecordSolverState{GradientDescentState{Array{Float64, 3}, Array{Float64, 3}, StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, IdentityUpdateRule, PolarRetraction}, NamedTuple{(:Iteration,), Tuple{RecordGroup}}}}, i::Int64)
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\solvers\solver.jl:149
 [19] solve!(p::DefaultManoptProblem{PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}}, s::Manopt.ReturnSolverState{RecordSolverState{GradientDescentState{Array{Float64, 3}, Array{Float64, 3}, StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, IdentityUpdateRule, PolarRetraction}, NamedTuple{(:Iteration,), Tuple{RecordGroup}}}})
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\solvers\solver.jl:136
 [20] gradient_descent!(M::PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, mgo::ManifoldGradientObjective{InplaceEvaluation, typeof(F), grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}}, p::Array{Float64, 3}; retraction_method::PolarRetraction, stepsize::ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, stopping_criterion::StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, debug::Vector{Any}, direction::IdentityUpdateRule, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:record, :return_state), Tuple{Vector{Symbol}, Bool}}})
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\solvers\gradient_descent.jl:256
 [21] gradient_descent!(M::PowerManifold{ℝ, Stiefel{ManifoldsBase.TypeParameter{Tuple{25, 3}}, ℝ}, Tuple{Int64}, ArrayPowerRepresentation}, f::Function, grad_f::grad!{SparseMatrixCSC{Float64, Int64}, Array{Float64, 3}}, p::Array{Float64, 3}; evaluation::InplaceEvaluation, kwargs::Base.Pairs{Symbol, Any, NTuple{4, Symbol}, NamedTuple{(:stepsize, :stopping_criterion, :record, :return_state), Tuple{ArmijoLinesearch{PolarRetraction, typeof(Manopt.armijo_initial_guess)}, StopWhenAny{Tuple{StopAfterIteration, StopWhenGradientNormLess}}, Vector{Symbol}, Bool}}})
    @ Manopt C:\Users\ddt00\.julia\packages\Manopt\LJaqp\src\solvers\gradient_descent.jl:229
 [22] top-level scope
    @ D:\code\paper_code\paper_admm_bm-main\test_powermanifolds.jl:73
 [23] include(fname::String)
    @ Base.MainInclude .\client.jl:478
 [24] top-level scope
    @ REPL[1]:1
in expression starting at D:\code\paper_code\paper_admm_bm-main\test_powermanifolds.jl:73

which seems that I need to define the retract_polar!() myself according to this page? I would expect the function can detect the submanifold (Stiefel) and call the retract function for Stiefel manifolds.

Btw, It's a bit confusing to me since I am running the same code I wrote about 1 year ago and there is no error like this one before. The packages related to manifold computation seem to change a lot since then.

kellertuer commented 9 months ago

Hi, thanks for using Manopt.jl – your error is indeed more related to Manifolds.jl, where both Stiefel and the PowerManifold are defined (Manopt.jl depends only on ManifoldsBase)

May I first ask which version are you using? We made retractions on PowerManifold recently a bit more easy to use, maybe we missed something.

I have two small remarks on the code: 1) If you like, we have a shorthand for power manifold, you can also write M = Stiefel(rnk,d)^100

2) Why are you defining

ManifoldsBase.default_retraction_method(M::Stiefel) = PolarRetraction()
ManifoldsBase.default_retraction_method(M::PowerManifold) = ManifoldsBase.default_retraction_method(M.manifold)
ManifoldsBase.default_inverse_retraction_method(M::Stiefel) = PolarInverseRetraction()
ManifoldsBase.default_inverse_retraction_method(M::PowerManifold) = ManifoldsBase.default_inverse_retraction_method(M.manifold)

usually it is better to use the retraction_method= keyword of the solver to set the retraction. Your power manifold definitions are actually the default, see for example https://github.com/JuliaManifolds/ManifoldsBase.jl/blob/81d2d93978527120dc95bc628a08f3197a3c6c9c/src/PowerManifold.jl#L458-L463.

But I can recreate your error and take a look (since I am surprised by that myself).

Concerning the packages changed a bit, we had one breaking change in Manifolds (Oct ‘23 version 0.9) in the last 2 years (v0.8 was in May 2022, which the Power Manifold is from. So yes there was a breaking change, but that should not be related to the power manifold. We sure introduced quite a few features over the last year, but we only in very rare cases do braking changes – the version 0l9 introduced more stable and faster storage of parameters in manifolds.

For Manopt.jl the last breaking change was in January, and that was mainly since we migrate function definitions (costs, gradients,...) either to ManifoldDiff.jl or to ManoptExamples.jl, since Manopt.jl should – at some point – only focus on the algorithms and algorithms elements (step sizes, stopping criteria,...); but also here we try to minimise the amount of breaking changes (since we are at 0.4, we had 3 since the beginning in 2016). Of course that is sometimes a tradeoff, but we try to as few as possible breaking changes, the one in January was to introduce the subsolvers properly I think.

edit: if you really want to be sure your 1-year-old code still runs, run it on the same versions as back then. Yes we introduced a small bug in the last change (see below), but that one as also a breaking change – so your old code might also just really break on a breaking change. That is what breaking changes do.

kellertuer commented 9 months ago

Ha. Found it!

We changed exactly two things in ManifoldsBase 0.15, and you it a combination of these that we missed ;)

and in this rework we also removed an one PowerRetraction we had for a while. And there we missed one, let‘s say “dispatch-path”. The PR linked directly above fixes this. If you are in a hurry, you can import and define

function ManifoldsBase.retract!(
    M::AbstractPowerManifold,
    q,
    p,
    X,
    t::Number,
    m::AbstractRetractionMethod = default_retraction_method(M, typeof(p)),
)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        retract!(
            M.manifold,
            _write(M, rep_size, q, i),
            _read(M, rep_size, p, i),
            _read(M, rep_size, X, i),
            t,
            m,
        )
    end
    return q
end

or you wait for ManifoldsBase.jl 0.15.4, which is just missing a test coverage check for now. Or use ManifoldsBase.jl 0.14.x together with Manifolds.jl 0.8.x until then :)

Thanks for pointing this out.

yuwenchen95 commented 9 months ago

Ha. Found it!

We changed exactly two things in ManifoldsBase 0.15, and you it a combination of these that we missed ;)

  • retractions follow since then the “allocate-then-dispatch” mode we do in the rest of the packages for quite some time (it reduced the code by half, and it is now easier to define new retractions)
  • Power manifolds moved to ManifoldsBase (so one can use them without loading all of Manifolds.jl).

and in this rework we also removed an one PowerRetraction we had for a while. And there we missed one, let‘s say “dispatch-path”. The PR linked directly above fixes this. If you are in a hurry, you can import and define

function ManifoldsBase.retract!(
    M::AbstractPowerManifold,
    q,
    p,
    X,
    t::Number,
    m::AbstractRetractionMethod = default_retraction_method(M, typeof(p)),
)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        retract!(
            M.manifold,
            _write(M, rep_size, q, i),
            _read(M, rep_size, p, i),
            _read(M, rep_size, X, i),
            t,
            m,
        )
    end
    return q
end

or you wait for ManifoldsBase.jl 0.15.4, which is just missing a test coverage check for now.

Thanks for pointing this out.

Hi,

I tried the fixing branch and it works now. Thank you so much for the quick update!

kellertuer commented 9 months ago

Great that that fixed everything already :)

And sure, the retractions are hopefully relatively easy to use – the code behind is a bit tricky, since we try to be flexible and efficient. So that one does not see bugs in that directly, might (but hopefully does not too often) happen. The same for the hopefully relatively easy to use power manifold.