Implement and document parallel transport on Grassmann. (#731)

* Implement and document PT in Grassmann. * Dispatch the allocating one correctly (in theory). * Finish PT on Grassmann * Simplify code. --------- Co-authored-by: Mateusz Baran <[email protected]>
JuliaManifolds · Jun 18, 2024 · 676b0f5 · 676b0f5 · kellertuer · Jun 18, 2024
1 parent 7785932
commit 676b0f5
Show file tree

Hide file tree

Showing 6 changed files with 97 additions and 28 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -5,6 +5,17 @@ All notable changes to this project 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.9.20] – 2024-06-17
+
+### Added
+
+* implemented parallel transport on the Grassmann manifold with respect to Stiefel representation
+
+### Changed
+
+* since now all exp/log/parallel transport are available for all representations of `Grassmann`,
+  these are now also set as defaults, since they are more exact.
+
 ## [0.9.19] – 2024-06-12
 
 ### Changed

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>", "Antoine Levitt <[email protected]>"]
-version = "0.9.19"
+version = "0.9.20"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

diff --git a/ext/ManifoldsTestExt/tests_general.jl b/ext/ManifoldsTestExt/tests_general.jl
@@ -477,8 +477,8 @@ function test_manifold(
         test_default_vector_transport && Test.@testset "default vector transport" begin
             v1t1 = vector_transport_to(M, pts[1], X1, pts32)
             v1t2 = vector_transport_direction(M, pts[1], X1, X2)
-            Test.@test is_vector(M, pts32, v1t1; atol=tvatol)
-            Test.@test is_vector(M, pts32, v1t2; atol=tvatol)
+            Test.@test is_vector(M, pts32, v1t1; atol=tvatol, error=:warn)
+            Test.@test is_vector(M, pts32, v1t2; atol=tvatol, error=:warn)
             Test.@test isapprox(M, pts32, v1t1, v1t2)
             Test.@test isapprox(M, pts[1], vector_transport_to(M, pts[1], X1, pts[1]), X1)
 
@@ -506,8 +506,10 @@ function test_manifold(
                 pts32 = retract(M, pts[1], X2, rtr_m)
                 test_to && (v1t1 = vector_transport_to(M, pts[1], X1, pts32, vtm))
                 test_dir && (v1t2 = vector_transport_direction(M, pts[1], X1, X2, vtm))
-                test_to && Test.@test is_vector(M, pts32, v1t1, true; atol=tvatol)
-                test_dir && Test.@test is_vector(M, pts32, v1t2, true; atol=tvatol)
+                test_to &&
+                    Test.@test is_vector(M, pts32, v1t1; atol=tvatol, error=:warn)
+                test_dir &&
+                    Test.@test is_vector(M, pts32, v1t2; atol=tvatol, error=:warn)
                 (test_to && test_dir) &&
                     Test.@test isapprox(M, pts32, v1t1, v1t2, atol=tvatol)
                 test_to && Test.@test isapprox(

diff --git a/src/manifolds/Grassmann.jl b/src/manifolds/Grassmann.jl
@@ -214,7 +214,7 @@ to a projector representation of said subspace, i.e. compute the [`canonical_pro
 for
 
 ```math
-  π^{\mathrm{SG}}(p) = pp^{\mathrm{T)}.
+  π^{\mathrm{SG}}(p) = pp^{\mathrm{T}}.
 ```
 """
 convert(::Type{ProjectorPoint}, p::AbstractMatrix) = ProjectorPoint(p * p')
@@ -232,24 +232,19 @@ for
 convert(::Type{ProjectorPoint}, p::StiefelPoint) = ProjectorPoint(p.value * p.value')
 
 """
+    default_retraction_method(M::Grassmann)
+    default_retraction_method(M::Grassmann, ::Type{StiefelPoint})
     default_retraction_method(M::Grassmann, ::Type{ProjectorPoint})
 
 Return `ExponentialRetraction` as the default on the [`Grassmann`](@ref) manifold
-with projection matrices
-"""
-default_retraction_method(::Grassmann, ::Type{ProjectorPoint}) = ExponentialRetraction()
+for both representations.
 """
-    default_retraction_method(M::Grassmann)
-    default_retraction_method(M::Grassmann, ::Type{StiefelPoint})
+default_retraction_method(::Grassmann) = ExponentialRetraction()
 
-Return `PolarRetracion` as the default on the [`Grassmann`](@ref) manifold
-with projection matrices
-"""
-default_retraction_method(::Grassmann) = PolarRetraction()
 """
     default_vector_transport_method(M::Grassmann)
 
-Return the `ProjectionTransport` as the default vector transport method
-for the [`Grassmann`](@ref) manifold.
+Return the default vector transport method for the [`Grassmann`](@ref) manifold,
+which is `ParallelTransport``()`.
 """
-default_vector_transport_method(::Grassmann) = ProjectionTransport()
+default_vector_transport_method(::Grassmann) = ParallelTransport()
diff --git a/src/manifolds/GrassmannStiefel.jl b/src/manifolds/GrassmannStiefel.jl
@@ -32,16 +32,17 @@ ManifoldsBase.@default_manifold_fallbacks (Stiefel{<:Any,ℝ}) StiefelPoint Stie
 ManifoldsBase.@default_manifold_fallbacks Grassmann StiefelPoint StiefelTVector value value
 
 function default_vector_transport_method(::Grassmann, ::Type{<:AbstractArray})
-    return ProjectionTransport()
+    return ParallelTransport()
 end
-default_vector_transport_method(::Grassmann, ::Type{<:StiefelPoint}) = ProjectionTransport()
+default_vector_transport_method(::Grassmann, ::Type{<:StiefelPoint}) = ParallelTransport()
 
 @doc raw"""
     distance(M::Grassmann, p, q)
 
 Compute the Riemannian distance on [`Grassmann`](@ref) manifold `M```= \mathrm{Gr}(n,k)``.
 
 The distance is given by
+
 ````math
 d_{\mathrm{Gr}(n,k)}(p,q) = \operatorname{norm}(\log_p(q)).
 ````
@@ -182,6 +183,54 @@ function log!(M::Grassmann, X, p, q)
     return X
 end
 
+@doc raw"""
+    parallel_transport_direction(M::Grassmann, p, X, Y)
+
+Compute the parallel transport of ``X \in   T_p\mathcal M`` along the
+geodesic starting in direction ``\dot γ (0) = Y``.
+
+ Let ``Y = USV`` denote the SVD decomposition of ``Y``.
+Then the parallel transport is given by the formula according to Equation (8.5) (p. 171) [AbsilMahonySepulchre:2008](@cite) as
+
+```math
+\mathcal P_{p,Y} X = -pV \sin(S)U^{\mathrm{T}}X + U\cos(S)U^{\mathrm{T}}X + (I-UU^{\mathrm{T}})X
+```
+
+where the sine and cosine applied to the diagonal matrix ``S`` are meant to be elementwise
+"""
+parallel_transport_direction(M::Grassmann, p, X, Y)
+
+# Hook into default since here we have direction first
+function parallel_transport_direction(M::Grassmann, p, X, Y)
+    Z = zero_vector(M, exp(M, p, X))
+    return parallel_transport_direction!(M, Z, p, X, Y)
+end
+
+function parallel_transport_direction!(M::Grassmann, Z, p, X, Y)
+    d = svd(Y)
+    return copyto!(
+        M,
+        Z,
+        p,
+        (-p * d.V .* sin.(d.S') + d.U .* cos.(d.S')) * (d.U' * X) + (I - d.U * d.U') * X,
+    )
+end
+
+@doc raw"""
+    parallel_transport_to(M::Grassmann, p, X, q)
+
+Compute the parallel transport of ``X ∈  T_p\mathcal M`` along the
+geodesic connecting ``p`` to ``q``.
+
+This method uses the [logarithmic map](@ref log(::Grassmann, ::Any...)) and the [parallel transport in that direction](@ref parallel_transport_direction(M::Grassmann, p, X, Y)).
+"""
+parallel_transport_to(M::Grassmann, p, X, q)
+
+function parallel_transport_to!(M::Grassmann, Z, p, X, q)
+    Y = log(M, p, q)
+    return parallel_transport_direction!(M, Z, p, X, Y)
+end
+
 @doc raw"""
     project(M::Grassmann, p)
 
@@ -329,7 +378,10 @@ end
 
 Compute the value of Riemann tensor on the real [`Grassmann`](@ref) manifold.
 The formula reads [Rentmeesters:2011](@cite)
-``R(X,Y)Z = (XY^\mathrm{T} - YX^\mathrm{T})Z + Z(Y^\mathrm{T}X - X^\mathrm{T}Y)``.
+
+```math
+R(X,Y)Z = (XY^\mathrm{T} - YX^\mathrm{T})Z + Z(Y^\mathrm{T}X - X^\mathrm{T}Y).
+```
 """
 riemann_tensor(::Grassmann{<:Any,ℝ}, p, X, Y, Z)
 
@@ -373,6 +425,14 @@ function uniform_distribution(M::Grassmann{<:Any,ℝ}, p)
     return ProjectedPointDistribution(M, d, (M, q, p) -> (q .= svd(p).U), p)
 end
 
+# switch order and not dispatch on the _to variant
+function vector_transport_direction(M::Grassmann, p, X, Y, ::ParallelTransport)
+    return parallel_transport_direction(M, p, X, Y)
+end
+function vector_transport_direction!(M::Grassmann, Z, p, X, Y, ::ParallelTransport)
+    return parallel_transport_direction!(M, Z, p, X, Y)
+end
+
 @doc raw"""
     vector_transport_to(M::Grassmann, p, X, q, ::ProjectionTransport)
 

diff --git a/test/manifolds/grassmann.jl b/test/manifolds/grassmann.jl
@@ -9,11 +9,11 @@ include("../header.jl")
             @test manifold_dimension(M) == 2
             @test !is_flat(M)
             @test is_flat(Grassmann(2, 1))
-            @test default_retraction_method(M) == PolarRetraction()
-            @test default_retraction_method(M, typeof(zeros(3, 2))) == PolarRetraction()
+            @test default_retraction_method(M) == ExponentialRetraction()
+            @test default_retraction_method(M, typeof(zeros(3, 2))) ==
+                  ExponentialRetraction()
             @test default_retraction_method(M, ProjectorPoint) == ExponentialRetraction()
-            @test default_retraction_method(M) == PolarRetraction()
-            @test default_vector_transport_method(M) == ProjectionTransport()
+            @test default_vector_transport_method(M) == ParallelTransport()
             @test get_total_space(M) == Stiefel(3, 2, ℝ)
             @test get_orbit_action(M) ==
                   Manifolds.RowwiseMultiplicationAction(M, Orthogonal(2))
@@ -82,7 +82,8 @@ include("../header.jl")
                 test_injectivity_radius=false,
                 test_project_tangent=true,
                 test_project_point=true,
-                test_default_vector_transport=false,
+                test_default_vector_transport=true,
+                vector_transport_methods=[ParallelTransport(), ProjectionTransport()],
                 point_distributions=[Manifolds.uniform_distribution(M, pts[1])],
                 test_vee_hat=false,
                 test_rand_point=true,
@@ -148,8 +149,8 @@ include("../header.jl")
         @testset "default_* functions" begin
             p = [1.0 0.0; 0.0 1.0; 0.0 0.0]
             pS = StiefelPoint(p)
-            @test default_vector_transport_method(M, typeof(p)) == ProjectionTransport()
-            @test default_vector_transport_method(M, typeof(pS)) == ProjectionTransport()
+            @test default_vector_transport_method(M, typeof(p)) == ParallelTransport()
+            @test default_vector_transport_method(M, typeof(pS)) == ParallelTransport()
         end
     end