Add horizontal and vertical component calculation to Kendall's shape …

…space (#582)

* Add horizontal and vertical component calculation to Kendall's shape space.

* one more exception

* fix in-place vertical_component!

* fix exp on KSS

* docs for horizontal_component on KSS

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.8.50"
+version = "0.8.51"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"

src/Manifolds.jl

@@ -307,7 +307,7 @@ using ManifoldDiff:
 import ManifoldDiff: riemannian_gradient, riemannian_gradient!
 
 using Markdown: @doc_str
-using MatrixEquations: lyapc
+using MatrixEquations: lyapc, sylvc
 using Quaternions: Quaternions
 using Random
 using RecipesBase
@@ -727,6 +727,8 @@ export ×,
     grad_euclidean_to_manifold!,
     hat,
     hat!,
+    horizontal_component,
+    horizontal_component!,
     horizontal_lift,
     horizontal_lift!,
     identity_element,
@@ -812,6 +814,8 @@ export ×,
     vector_transport_to!,
     vee,
     vee!,
+    vertical_component,
+    vertical_component!,
     zero_vector,
     zero_vector!
 # Lie group types & functions

src/manifolds/KendallsShapeSpace.jl

@@ -66,7 +66,8 @@ See [^Guigui2021] for discussion about its computation.
 exp(M::KendallsShapeSpace, p, X)
 
 function exp!(M::KendallsShapeSpace, q, p, X)
-    return exp!(get_embedding(M), q, p, X)
+    Xh = horizontal_component(M, p, X)
+    return exp!(get_embedding(M), q, p, Xh)
 end
 
 embed(::KendallsShapeSpace, p) = p
@@ -82,6 +83,28 @@ function get_embedding(::KendallsShapeSpace{N,K}) where {N,K}
     return KendallsPreShapeSpace(N, K)
 end
 
+"""
+    horizontal_component(::KendallsShapeSpace, p, X)
+
+Compute the horizontal component of tangent vector `X` at `p` on [`KendallsShapeSpace`](@ref)
+`M`. See [^Guigui2021], Section 2.3 for details.
+"""
+horizontal_component(::KendallsShapeSpace, p, X)
+
+function horizontal_component!(::KendallsShapeSpace, Y, p, X)
+    B = p * transpose(p)
+    C = X * transpose(p) - p * transpose(X)
+    A = sylvc(B, B, C)
+    Y .= X .- A * p
+    return Y
+end
+
+function inner(M::KendallsShapeSpace, p, X, Y)
+    Xh = horizontal_component(M, p, X)
+    Yh = horizontal_component(M, p, Y)
+    return inner(get_embedding(M), p, Xh, Yh)
+end
+
 function _isapprox(M::KendallsShapeSpace, p, X, Y; atol=sqrt(max_eps(X, Y)), kwargs...)
     return isapprox(norm(M, p, X - Y), 0; atol=atol, kwargs...)
 end
@@ -114,11 +137,26 @@ end
 @doc raw"""
     manifold_dimension(M::KendallsShapeSpace)
 
-Return the dimension of the [`KendallsShapeSpace`](@ref) manifold `M`. The dimension is given by
-``n(k - 1) - 1 - n(n - 1)/2``.
+Return the dimension of the [`KendallsShapeSpace`](@ref) manifold `M`. The dimension is
+given by ``n(k - 1) - 1 - n(n - 1)/2`` in the typical case where ``k \geq n+1``, and
+``(k + 1)(k - 2) / 2`` otherwise, unless ``k`` is equal to 1, in which case the dimension
+is 0. See [^Kendall1984] for a discussion of the over-dimensioned case.
 """
 function manifold_dimension(::KendallsShapeSpace{n,k}) where {n,k}
-    return n * (k - 1) - 1 - div(n * (n - 1), 2)
+    if k < n + 1 # over-dimensioned case
+        if k == 1
+            return 0
+        else
+            return div((k + 1) * (k - 2), 2)
+        end
+    else
+        return n * (k - 1) - 1 - div(n * (n - 1), 2)
+    end
+end
+
+function norm(M::KendallsShapeSpace, p, X)
+    Xh = horizontal_component(M, p, X)
+    return norm(get_embedding(M), p, Xh)
 end
 
 function project!(M::KendallsShapeSpace, q, p)

src/manifolds/QuotientManifold.jl

@@ -145,9 +145,9 @@ get_orbit_action(::AbstractManifold)
     horizontal_lift(N::AbstractManifold, q, X)
     horizontal_lift(::QuotientManifold{𝔽,MT<:AbstractManifold{𝔽},NT<:AbstractManifold}, p, X) where {𝔽}
 
-Given a point `q` such that ``p=π(q)`` is a point on a quotient manifold `M`
-(implicitly given for the first case) and a tangent vector `X` this method
-computes a tangent vector `Y` on the horizontal space of ``T_q\mathcal N``,
+Given a point `q` in total space of quotient manifold `N` such that ``p=π(q)`` is a point on
+a quotient manifold `M` (implicitly given for the first case) and a tangent vector `X` this
+method computes a tangent vector `Y` on the horizontal space of ``T_q\mathcal N``,
 i.e. the subspace that is orthogonal to the kernel of ``Dπ(q)``.
 """
 function horizontal_lift(N::AbstractManifold, q, X)
@@ -156,14 +156,42 @@ function horizontal_lift(N::AbstractManifold, q, X)
 end
 
 @doc raw"""
-    horizontal_lift(N, q, X)
-    horizontal_lift(QuotientManifold{M,N}, p, X)
+    horizontal_lift!(N, Y, q, X)
+    horizontal_lift!(QuotientManifold{M,N}, Y, p, X)
 
 Compute the [`horizontal_lift`](@ref) of `X` from ``T_p\mathcal M``, ``p=π(q)``.
 to ``T_q\mathcal N` in place of `Y`.
 """
 horizontal_lift!(N::AbstractManifold, Y, q, X)
 
+"""
+    horizontal_component(N::AbstractManifold, p, X)    
+
+Compute the horizontal component of tangent vector `X` at point `p`
+in the total space of quotient manifold `N`.
+"""
+function horizontal_component(N::AbstractManifold, p, X)
+    Y = allocate_result(N, horizontal_component, X, p)
+    return horizontal_component!(N, Y, p, X)
+end
+
 function Base.show(io::IO, M::QuotientManifold)
     return print(io, "QuotientManifold($(M.manifold), $(M.total_space))")
 end
+
+"""
+    vertical_component(N::AbstractManifold, p, X)    
+
+Compute the vertical component of tangent vector `X` at point `p`
+in the total space of quotient manifold `N`.
+"""
+function vertical_component(N::AbstractManifold, p, X)
+    return X - horizontal_component(N, p, X)
+end
+
+function vertical_component!(N::AbstractManifold, Y, p, X)
+    horizontal_component!(N, Y, p, X)
+    Y .*= -1
+    Y .+= X
+    return Y
+end

test/manifolds/shape_space.jl

@@ -56,8 +56,36 @@ end
         0.3248027612629014 0.440253011955812 -0.7650557732187135
         0.26502337825226757 -0.06175142812400016 -0.20327195012826738
     ]
+    X1 = [
+        0.6090792159558263 -0.02523987621672985 -0.5838393397390964
+        0.4317628895706799 0.12108361184633629 -0.5528465014170161
+    ]
+    X1h = [
+        0.5218590427922166 0.05525104866717821 -0.5771100914593948
+        0.3319078589730016 0.2777009756923593 -0.6096088346653609
+    ]
+    X1v = [
+        0.08722017316360964 -0.08049092488390806 -0.006729248279701561
+        0.09985503059767825 -0.156617363846023 0.05676233324834479
+    ]
+    @testset "tangent vector components" begin
+        @test isapprox(M, p1, horizontal_component(M, p1, X1), X1h)
+        @test isapprox(M, p1, vertical_component(M, p1, X1), X1v)
+        Y = similar(X1)
+        vertical_component!(M, Y, p1, X1)
+        @test isapprox(M, p1, Y, X1v)
+        @test norm(M, p1, X1v) < 1e-16
+        @test abs(norm(M, p1, X1) - norm(M, p1, X1h)) < 1e-16
+    end
+
     @test_throws ManifoldDomainError is_point(M, [1 0 1; 1 -1 0], true)
     @test_throws ManifoldDomainError is_vector(M, p1, [1 0 1; 1 -1 0], true)
+
+    @testset "exp/distance/norm" begin
+        q1 = exp(M, p1, X1)
+        @test distance(M, p1, q1) ≈ norm(M, p1, X1)
+    end
+
     test_manifold(
         M,
         [p1, p2, p3];
@@ -71,4 +99,10 @@ end
         test_rand_tvector=true,
         rand_tvector_atol_multiplier=5,
     )
+    @testset "degenerate cases" begin
+        Md3_2 = KendallsShapeSpace(3, 2)
+        Md2_1 = KendallsShapeSpace(2, 1)
+        @test manifold_dimension(Md3_2) == 0
+        @test manifold_dimension(Md2_1) == 0
+    end
 end