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
366 stars 52 forks source link

Why does exponential map on the Sphere return points outside of the Sphere? #715

Closed Nimrais closed 5 months ago

Nimrais commented 5 months ago

I am using Manifolds.jl v0.9.15.

The following is confusing for me

julia> using Manifolds

julia> sp = Sphere(1)
Sphere(1, ℝ)

julia> p = rand(sp)
2-element Vector{Float64}:
 0.8953358144881157
 0.4453917144434801

julia> exp(sp, p, [1, 0])
2-element Vector{Float64}:
 1.3252229899021544
 0.24064617032837635

julia> norm(exp(sp, p, [1, 0]))
1.3468951526599684

What am I doing wrong? I thought that the exponential map should return points that are staying on the manifold.

kellertuer commented 5 months ago

Your point p is not on the sphere, it is not of norm 1. The behaviour of the exponential map for non-points and/or non-tangents is undefined.

I first thought your point is not unit norm, but it actually is. Still, your tangent vector has to be orthogonal to p, and that is not the case. But I can actually quite well explain the behaviour. The resulting point stays on the same sphere with the radius as p has.

See

julia> is_vector(Sphere(1), p, [1,0]; error=:error)
ERROR: DomainError with 0.8953358144881157:
The vector [1, 0] is not a tangent vector to [0.8953358144881157, 0.4453917144434801] on Sphere(1, ℝ), since it is not orthogonal in the embedding.

For non-points I know the exponential map (with tangent vectors) stays on the same radius as norm(p), for non-tangent vectors I do not know what exp does, but since it is not defined, it is quite often the case that the point is not on the sphere.

kellertuer commented 5 months ago

To add: If you take a tangent vector

julia> X = [-p[2], p[1]]
2-element Vector{Float64}:
 -0.4453917144434801
  0.8953358144881157

julia> is_vector(Sphere(1), p, X, true, true)
true

you also get

julia> q = exp(Sphere(1), p, X)
2-element Vector{Float64}:
 0.10896780051622529
 0.9940452798794712

julia> is_point(Sphere(1), q)
true
Nimrais commented 5 months ago

Ahh, yes thanks for the fast answer, I see the issue now!

kellertuer commented 5 months ago

Great :) We specifically provided the two functions is_point and is_vector to help a bit with checks.

On the other hand we do not always check all input. You can activate that (to narrow down an error for example) by wrapping your manifold in a ValidationManifold https://juliamanifolds.github.io/ManifoldsBase.jl/stable/manifolds/#A-manifold-for-validation, that actually calls the above two functions on inputs and results within most functions of Manifolds.jl