JuliaManifolds / Manifolds.jl

Manifolds.jl provides a library of manifolds aiming for an easy-to-use and fast implementation.
https://juliamanifolds.github.io/Manifolds.jl
MIT License
368 stars 53 forks source link

retract_project not implemented for ProductManifold #659

Closed blegat closed 11 months ago

blegat commented 11 months ago
  MethodError: no method matching retract_project!(::ProductManifold{ℝ, Tuple{Sphere{1, ℝ}, Euclidean{Tuple{0}, ℝ}}}, ::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, ::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, ::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, ::Float64)

  Closest candidates are:
    retract_project!(::Spectrahedron, ::Any, ::Any, ::Any, ::Number)
     @ Manifolds ~/.julia/packages/Manifolds/KpLKJ/src/manifolds/Spectrahedron.jl:149
    retract_project!(::Stiefel, ::Any, ::Any, ::Any, ::Number)
     @ Manifolds ~/.julia/packages/Manifolds/KpLKJ/src/manifolds/StiefelEuclideanMetric.jl:182
    retract_project!(::GeneralizedGrassmann, ::Any, ::Any, ::Any, ::Number)
     @ Manifolds ~/.julia/packages/Manifolds/KpLKJ/src/manifolds/GeneralizedGrassmann.jl:372
    ...

  Stacktrace:
    [1] retract_project(M::ProductManifold{ℝ, Tuple{Sphere{1, ℝ}, Euclidean{Tuple{0}, ℝ}}}, p::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, X::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, t::Float64; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
      @ ManifoldsBase ~/.julia/packages/ManifoldsBase/vKGuo/src/retractions.jl:1179
    [2] retract_project(M::ProductManifold{ℝ, Tuple{Sphere{1, ℝ}, Euclidean{Tuple{0}, ℝ}}}, p::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, X::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, t::Float64)
      @ ManifoldsBase ~/.julia/packages/ManifoldsBase/vKGuo/src/retractions.jl:1177
    [3] _retract(M::ProductManifold{ℝ, Tuple{Sphere{1, ℝ}, Euclidean{Tuple{0}, ℝ}}}, p::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, X::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, t::Float64, ::ProjectionRetraction; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
      @ ManifoldsBase ~/.julia/packages/ManifoldsBase/vKGuo/src/retractions.jl:1110
    [4] _retract(M::ProductManifold{ℝ, Tuple{Sphere{1, ℝ}, Euclidean{Tuple{0}, ℝ}}}, p::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, X::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, t::Float64, ::ProjectionRetraction)
      @ ManifoldsBase ~/.julia/packages/ManifoldsBase/vKGuo/src/retractions.jl:1109
    [5] retract(M::ProductManifold{ℝ, Tuple{Sphere{1, ℝ}, Euclidean{Tuple{0}, ℝ}}}, p::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, X::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, t::Float64, m::ProjectionRetraction)
      @ ManifoldsBase ~/.julia/packages/ManifoldsBase/vKGuo/src/retractions.jl:1088
    [6] retract(::ManifoldsBase.EmptyTrait, M::ProductManifold{ℝ, Tuple{Sphere{1, ℝ}, Euclidean{Tuple{0}, ℝ}}}, p::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, X::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, t::Float64, m::ProjectionRetraction)
      @ ManifoldsBase ~/.julia/packages/ManifoldsBase/vKGuo/src/nested_trait.jl:346
    [7] retract (repeats 2 times)
      @ ~/.julia/packages/ManifoldsBase/vKGuo/src/nested_trait.jl:313 [inlined]
    [8] retract
      @ ~/.julia/packages/ManifoldsBase/vKGuo/src/nested_trait.jl:306 [inlined]
    [9] linesearch_backtrack(M::ProductManifold{ℝ, Tuple{Sphere{1, ℝ}, Euclidean{Tuple{0}, ℝ}}}, f::Manopt.var"#104#105"{Manopt.DefaultManoptProblem{ProductManifold{ℝ, Tuple{Sphere{1, ℝ}, Euclidean{Tuple{0}, ℝ}}}, Manopt.ManifoldGradientObjective{Manopt.AllocatingEvaluation, Macaulay.var"#eval_f_cb#96"{Matrix{Float64}}, Macaulay.var"#eval_grad_f_cb#97"{Macaulay.var"#partition#101"{Int64}, Matrix{Float64}}}}}, p::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, grad_f_at_p::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, s::Float64, decrease::Float64, contract::Float64, retr::ProjectionRetraction, η::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, f0::Float64; stop_when_stepsize_less::Float64, stop_when_stepsize_exceeds::Float64, stop_increasing_at_step::Int64, stop_decreasing_at_step::Int64)
      @ Manopt ~/.julia/packages/Manopt/A4a09/src/plans/stepsize.jl:368
   [10] (::Manopt.ArmijoLinesearch{ProjectionRetraction, typeof(Manopt.armijo_initial_guess)})(mp::Manopt.DefaultManoptProblem{ProductManifold{ℝ, Tuple{Sphere{1, ℝ}, Euclidean{Tuple{0}, ℝ}}}, Manopt.ManifoldGradientObjective{Manopt.AllocatingEvaluation, Macaulay.var"#eval_f_cb#96"{Matrix{Float64}}, Macaulay.var"#eval_grad_f_cb#97"{Macaulay.var"#partition#101"{Int64}, Matrix{Float64}}}}, s::Manopt.GradientDescentState{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Manopt.StopWhenAny{Tuple{Manopt.StopAfterIteration, Manopt.StopWhenGradientNormLess}}, Manopt.ArmijoLinesearch{ProjectionRetraction, typeof(Manopt.armijo_initial_guess)}, Manopt.IdentityUpdateRule, ProjectionRetraction}, i::Int64, η::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
      @ Manopt ~/.julia/packages/Manopt/A4a09/src/plans/stepsize.jl:288
   [11] (::Manopt.ArmijoLinesearch{ProjectionRetraction, typeof(Manopt.armijo_initial_guess)})(mp::Manopt.DefaultManoptProblem{ProductManifold{ℝ, Tuple{Sphere{1, ℝ}, Euclidean{Tuple{0}, ℝ}}}, Manopt.ManifoldGradientObjective{Manopt.AllocatingEvaluation, Macaulay.var"#eval_f_cb#96"{Matrix{Float64}}, Macaulay.var"#eval_grad_f_cb#97"{Macaulay.var"#partition#101"{Int64}, Matrix{Float64}}}}, s::Manopt.GradientDescentState{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Manopt.StopWhenAny{Tuple{Manopt.StopAfterIteration, Manopt.StopWhenGradientNormLess}}, Manopt.ArmijoLinesearch{ProjectionRetraction, typeof(Manopt.armijo_initial_guess)}, Manopt.IdentityUpdateRule, ProjectionRetraction}, i::Int64, η::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}})
      @ Manopt ~/.julia/packages/Manopt/A4a09/src/plans/stepsize.jl:279
   [12] (::Manopt.ArmijoLinesearch{ProjectionRetraction, typeof(Manopt.armijo_initial_guess)})(mp::Manopt.DefaultManoptProblem{ProductManifold{ℝ, Tuple{Sphere{1, ℝ}, Euclidean{Tuple{0}, ℝ}}}, Manopt.ManifoldGradientObjective{Manopt.AllocatingEvaluation, Macaulay.var"#eval_f_cb#96"{Matrix{Float64}}, Macaulay.var"#eval_grad_f_cb#97"{Macaulay.var"#partition#101"{Int64}, Matrix{Float64}}}}, s::Manopt.GradientDescentState{ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Manopt.StopWhenAny{Tuple{Manopt.StopAfterIteration, Manopt.StopWhenGradientNormLess}}, Manopt.ArmijoLinesearch{ProjectionRetraction, typeof(Manopt.armijo_initial_guess)}, Manopt.IdentityUpdateRule, ProjectionRetraction}, i::Int64)
      @ Manopt ~/.julia/packages/Manopt/A4a09/src/plans/stepsize.jl:27
kellertuer commented 11 months ago

This should be resolved once https://github.com/JuliaManifolds/ManifoldsBase.jl/pull/169 is merged and the product manifold. Then the product manifold is already defined in ManifoldsBase.jl.

I am indeed a bit surprised that this error appears. Without context / versions it's also not so easy to comment, but I had hoped we passed all retractions on (to each factor of the product already), but maybe we missed a retraction type until now (and I have not yet run into this).

blegat commented 11 months ago

It worked by with the default retraction and I got the error when getting ManifoldsBase.ProjectionRetraction. Feeld free to close then if it's going to be resolved by https://github.com/JuliaManifolds/ManifoldsBase.jl/pull/169

kellertuer commented 11 months ago

Ah, no it is not resolved by the linked PR, since the problem lies somewhere else. Let me illustrate.

If you do

using Manifolds
M = Sphere(2) × ℝ^3

(approximating your idea) you get a

ProductManifold with 2 submanifolds:
 Sphere(2, ℝ)
 Euclidean(3; field = ℝ)

Now note that the default retraction looks like (we could surely improve on the printing)

julia> default_retraction_method(M)
ProductRetraction{Tuple{ExponentialRetraction, ExponentialRetraction}}((ExponentialRetraction(), ExponentialRetraction()))

Since on Product manifolds you always have to provide a ProductRetraction as well, so what you need is not a ProjectionRetraction but for example

rtm = ProductRetraction(ProjectionRetraction(), ExponentialRetraction())

which I already allowed myself to adapt in the second component, so this does the projection based retraction on the sphere and still good-old + on the Euclidean part.

We currently do not have the “automatism” to map a retraction type (like ProjectionRetraction) to an n-copies-vector of itself passed to ProductRetraction.

You are probably correct that besides above show method also the docs of ProductRetraction could be improved.

edit: My confusion came from the fact that I used much more often the PowerManifold, but since there we have basically copies of a single manifold, there is no need for a PowerRetraction there, since all are the same anyways. For the product manifold this is not the case. There might even be products of manifolds, where no retraction exists on both manifolds. We can of course still discuss whether a single retraction should be passed down (to each factor of the product), but the error message of that – if it fails – might be misleading.

mateuszbaran commented 11 months ago

I think we could just implement the projection retraction on product manifold, this sounds fine.

kellertuer commented 11 months ago

But then I would vote for passing down any other retraction as well (i.e. like for power)?

mateuszbaran commented 11 months ago

IIRC there were some ambiguity issues that make passing arbitrary retractions on ProductManifold hard to do right. Passing a few particular ones should be fine though.

kellertuer commented 11 months ago

Ah, right, though I do not remember whether that stayed after the 3-layer structure. We can just try when first merging the product manifold PR in ManifoldsBase and then checking on the retraction-restructure PR :)

kellertuer commented 11 months ago

This now already works on ManifoldsBase main branch (i.e. that retract(M, p, X, ProjectionRetraciont()) works on product manifolds element wise, same for all other retractions). This will be part of the forthcoming release ManifoldsBase.jl 0.15), so I'll already close this issue.