diff --git a/NEWS.md b/NEWS.md index 383156383..944bc2b8f 100644 --- 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 index 9e3f537a5..3422c26e1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Manifolds" uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann ", "Antoine Levitt "] -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 index c600d00bd..f0cad921f 100644 --- 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 index b96cb58a8..d96fbefc5 100644 --- 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 index d90397b9d..14ded4494 100644 --- a/src/manifolds/GrassmannStiefel.jl +++ b/src/manifolds/GrassmannStiefel.jl @@ -32,9 +32,9 @@ 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) @@ -42,6 +42,7 @@ default_vector_transport_method(::Grassmann, ::Type{<:StiefelPoint}) = Projectio 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 index 9f6827b27..c300e8388 100644 --- 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