Adjust injectivity radii for (inverse) retractions on the Sphere (#900)

NEWS.md

@@ -5,6 +5,15 @@ All notable changes to ´Manifolds.jl´ will be documented in this file.
 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.11.30] Unreleased
+
+### Added
+
+* `injectivity_radius(M, p, m)` and `injectivity_radius(M, m)` where `m` is an `AbstractInverseRetractionMethod` and `M` is a `Sphere` are introduced.
+* `injectivity_radius(M, m)` where `M` is a `Sphere` and `m` is a `LogarithmicInverseRetraction` or a `ProjectionInverseRetraction`.
+* `bijectivity_radius(M, p, m)` and `bijectivity_radius(M, m)` are defined for `Sphere`s.
+* `Manifolds.Test.test_manifold` now tests injectivity radii for inverse retractions.
+
 ## [0.11.29] Unreleased
 
 ### Fixed

ext/ManifoldsTestExt/ManifoldsTestExt.jl

@@ -63,7 +63,7 @@ Possible entries of the `expectations` dictionary are
 * for `injectivity_radius`
     - you can provide the global test for just the function,
     - you can provide the radius for a specific point with `(injectivity_radius, p)`
-    - you can provide a global or local one for retractions as well usinr `(injectivity_radius, rm)` and `injectivity_radius, p, rm`, respectively
+    - you can provide a global or local one for retractions as well using `(injectivity_radius, rm)` and `injectivity_radius, p, rm`, respectively
 * for `get_basis`, the key is a tuple of the function and the basis, e.g. `(get_basis, B) => ...` to the expected basis
 * for `get_coordinates` the key is a tuple of the function and the basis, e.g. `(get_coordinates, B) => c`
 * for `get_vector` the key is a tuple of the function, the coordinate vector, and the basis, e.g. `(get_vector, c, B) => X`
@@ -315,7 +315,18 @@ function Manifolds.Test.test_manifold(M::AbstractManifold, properties::Dict, exp
                     name = "injectivity_radius(M, p, $rm)", # shorten name within large suite
                 )
             end
-
+            for irm in inverse_retraction_methods
+                ismissing(irm) && continue
+                expected_irm = get_expectation(expectations, (injectivity_radius, points[1], irm))
+                expected_irm_global = get_expectation(expectations, (injectivity_radius, irm))
+                Manifolds.Test.test_injectivity_radius(
+                    M, points[1];
+                    expected_value = expected_irm,
+                    expected_global_value = expected_irm_global,
+                    retraction_method = irm,
+                    name = "injectivity_radius(M, p, $irm)",
+                )
+            end
         end
         if (inner in functions) && !ismissing(vector)
             expected_inner = get_expectation(expectations, inner)

src/Manifolds.jl

@@ -308,6 +308,7 @@ using ManifoldsBase:
     SoftmaxRetraction,
     SoftmaxInverseRetraction,
     StabilizedRetraction,
+    StabilizedInverseRetraction,
     StopForwardingType,
     TangentSpace,
     TangentSpaceType,
@@ -756,7 +757,10 @@ export AbstractInverseRetractionMethod,
     PolarLightInverseRetraction,
     ProjectionInverseRetraction,
     ShootingInverseRetraction,
-    SoftmaxInverseRetraction
+    SoftmaxInverseRetraction,
+    StabilizedInverseRetraction
+# Bijectivity radii are exported here to be tested, but this should be removed if bijectivity radii get moved to ManifoldsBase.
+export bijectivity_radius
 # Estimation methods for median and mean
 export AbstractApproximationMethod,
     GradientDescentEstimation,

src/manifolds/Sphere.jl

@@ -286,7 +286,7 @@ function get_vector_orthonormal!(M::AbstractSphere{ℝ}, Y, p, c, ::RealNumbers)
 end
 
 _doc_injectivity_radius_sphere = raw"""
-    injectivity_radius(M::AbstractSphere[, p, ::ExponentialRetraction])
+    injectivity_radius(M::AbstractSphere[, p, ::StabilizedRetraction])
 
 Return the injectivity radius for the [`AbstractSphere`](@ref) `M`, which is globally ``π``.
 """
@@ -301,12 +301,11 @@ _doc_injectivity_radius_sphere_projection = raw"""
     injectivity_radius(M::Sphere, p, ::ProjectionRetraction)
 
 Return the injectivity radius for the [`ProjectionRetraction`](@extref `ManifoldsBase.ProjectionRetraction`) on the
-[`AbstractSphere`](@ref), which is globally ``\frac{π}{2}``.
+[`AbstractSphere`](@ref), which is globally ``\infty``.
 """
 
 @doc "$(_doc_injectivity_radius_sphere_projection)"
 injectivity_radius(::AbstractSphere, ::ProjectionRetraction)
-
 @doc "$(_doc_injectivity_radius_sphere_projection)"
 injectivity_radius(::AbstractSphere, p, ::ProjectionRetraction)
 
@@ -317,8 +316,8 @@ end
 function injectivity_radius(M::AbstractSphere, p, m::AbstractRetractionMethod)
     return _injectivity_radius(M, p, m)
 end
-_injectivity_radius(::AbstractSphere, ::ExponentialRetraction) = π
-_injectivity_radius(::AbstractSphere, ::ProjectionRetraction) = π / 2
+_injectivity_radius(::AbstractSphere, ::StabilizedRetraction) = π
+_injectivity_radius(::AbstractSphere, ::ProjectionRetraction) = Inf
 
 @doc raw"""
     inverse_retract(M::AbstractSphere, p, q, ::ProjectionInverseRetraction)
@@ -337,6 +336,62 @@ function inverse_retract_project!(::AbstractSphere, X, p, q)
     return (X .= q ./ real(dot(p, q)) .- p)
 end
 
+_doc_injectivity_radius_sphere_inverse_projection = raw"""
+    injectivity_radius(M::AbstractSphere, p, ::ProjectionInverseRetraction)
+
+Return the injectivity radius for the [`ProjectionInverseRetraction`](@extref `ManifoldsBase.ProjectionInverseRetraction`) on the [`AbstractSphere`](@ref), which is the largest geodesic distance ``ξ`` such that the inverse projection retraction from `p` is invertible whenever ``d_{𝕊^{d-1}}(p,q) ≤ ξ``. This is globally ``\frac{\pi}{2}``.
+"""
+
+@doc "$(_doc_injectivity_radius_sphere_inverse_projection)"
+injectivity_radius(::AbstractSphere, ::ProjectionInverseRetraction)
+
+@doc "$(_doc_injectivity_radius_sphere_inverse_projection)"
+injectivity_radius(::AbstractSphere, p, ::ProjectionInverseRetraction)
+
+function injectivity_radius(M::AbstractSphere, m::AbstractInverseRetractionMethod)
+    return _injectivity_radius(M, m)
+end
+function injectivity_radius(M::AbstractSphere, p, m::AbstractInverseRetractionMethod)
+    return _injectivity_radius(M, p, m)
+end
+
+function _injectivity_radius(M::AbstractSphere, m::LogarithmicInverseRetraction)
+    return injectivity_radius(M, StabilizedRetraction())
+end
+function _injectivity_radius(M::AbstractSphere, p, m::LogarithmicInverseRetraction)
+    return injectivity_radius(M, p, StabilizedRetraction())
+end
+_injectivity_radius(::AbstractSphere, ::ProjectionInverseRetraction) = π / 2
+function _injectivity_radius(M::AbstractSphere, p, m::ProjectionInverseRetraction)
+    return _injectivity_radius(M, m)
+end
+
+_doc_bijectivity_radius_sphere = raw"""
+    bijectivity_radius(M::AbstractSphere[, p, m::Union{AbstractRetractionMethod, AbstractInverseRetractionMethod}])
+
+Return the bijectivity radius for the abstract sphere `M` endowed with (inverse) retraction method `m`, which is the smallest radius applicable to both the tangent spaces and the Riemannian metric for which the (inverse) retraction `m` has a well-defined inverse. This radius is computed as the minimum between the injectivity radius of `M` endowed with (inverse) retraction method `m` (at `p`), and the injectivity radius of the inverse of `m` (at `p`).
+"""
+
+@doc "$(_doc_bijectivity_radius_sphere)"
+bijectivity_radius(::AbstractSphere, ::Union{AbstractRetractionMethod, AbstractInverseRetractionMethod})
+@doc"$(_doc_bijectivity_radius_sphere)"
+bijectivity_radius(::AbstractSphere, p, ::Union{AbstractRetractionMethod, AbstractInverseRetractionMethod})
+
+function bijectivity_radius(M::AbstractSphere, m::Union{AbstractRetractionMethod, AbstractInverseRetractionMethod})
+    if m ∈ (StabilizedRetraction(), StabilizedInverseRetraction())
+        return injectivity_radius(M)
+    elseif m ∈ (ProjectionRetraction(), ProjectionInverseRetraction())
+        return injectivity_radius(M, ProjectionInverseRetraction())
+    end
+end
+function bijectivity_radius(M::AbstractSphere, p, m::Union{AbstractRetractionMethod, AbstractInverseRetractionMethod})
+    if m ∈ (StabilizedRetraction(), StabilizedInverseRetraction())
+        return injectivity_radius(M, p)
+    elseif m ∈ (ProjectionRetraction(), ProjectionInverseRetraction())
+        return injectivity_radius(M, p, ProjectionInverseRetraction())
+    end
+end
+
 """
     is_flat(M::AbstractSphere)
 

test/manifolds/test_sphere.jl

@@ -33,7 +33,7 @@ using ManifoldDiff
             :InvalidVectors => [p],
             :NormalVectors => [V],
             :Covectors => [ξ],
-            :RetractionMethods => [ExponentialRetraction(), ProjectionRetraction()],
+            :RetractionMethods => [StabilizedRetraction(), ProjectionRetraction()],
             :VectorTransportMethods => [ParallelTransport(), SchildsLadderTransport(), PoleLadderTransport()]
         ),
         # Expectations
@@ -48,7 +48,9 @@ using ManifoldDiff
             (get_vectors, DefaultOrthonormalBasis()) => :Orthonormal,
             (get_vectors, DefaultOrthogonalBasis()) => :Orthogonal,
             injectivity_radius => π,
-            (injectivity_radius, ProjectionRetraction()) => π / 2,
+            (injectivity_radius, ProjectionRetraction()) => Inf,
+            (injectivity_radius, LogarithmicInverseRetraction()) => π,
+            (injectivity_radius, ProjectionInverseRetraction()) => π / 2,
             is_default_metric => EuclideanMetric(),
             log => X, norm => π / 4,
             parallel_transport_to => parallel_transport_to(M, p, X, q),