Bases for CholeskySpace and LogCholesky metric (#780)

* Bases for CholeskySpace and LogCholesky metric

* add tests and address review

* rename function

NEWS.md

@@ -5,7 +5,12 @@ 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.10.12] - unreleased
+## [0.10.12] - 2025-01-10
+
+### Added
+
+* Orthonormal bases for `CholeskySpace` and `LogCholesky` metric for `SymmetricPositiveDefinite`.
+* `rand` for `CholeskySpace`.
 
 ### Changed
 

docs/src/misc/internals.md

@@ -6,6 +6,7 @@ This page documents the internal types and methods of `Manifolds.jl`'s that migh
 
 ```@docs
 Manifolds.eigen_safe
+Manifolds.get_parameter_type
 Manifolds.isnormal
 Manifolds.log_safe
 Manifolds.log_safe!

src/Manifolds.jl

@@ -370,6 +370,18 @@ include("utils.jl")
 
 include("manifold_fallbacks.jl")
 
+"""
+    get_parameter_type(M::AbstractManifold)
+
+Get `parameter` argument of the constructor of manifold `M`.
+Returns either `:field` or `:type`.
+
+# See also
+
+`get_parameter`, `TypeParameter`
+"""
+get_parameter_type(::AbstractManifold)
+
 """
     projected_distribution(M::AbstractManifold, d, [p=rand(d)])
 

src/manifolds/CholeskySpace.jl

@@ -2,7 +2,7 @@
     CholeskySpace{T} <: AbstractManifold{ℝ}
 
 The manifold of lower triangular matrices with positive diagonal and
-a metric based on the cholesky decomposition. The formulae for this manifold
+a metric based on the Cholesky decomposition. The formulae for this manifold
 are for example summarized in Table 1 of [Lin:2019](@cite).
 
 # Constructor
@@ -121,11 +121,38 @@ function exp!(::CholeskySpace, q, p, X)
     return q
 end
 
+function get_coordinates_orthonormal!(M::CholeskySpace, Xⁱ, p, X, ::RealNumbers)
+    n = get_parameter(M.size)[1]
+    view(Xⁱ, 1:n) .= diag(X)
+    xi_ind = n + 1
+    for i in 1:n
+        for j in (i + 1):n
+            Xⁱ[xi_ind] = X[j, i]
+            xi_ind += 1
+        end
+    end
+    return Xⁱ
+end
+
+function get_vector_orthonormal!(M::CholeskySpace, X, p, Xⁱ, ::RealNumbers)
+    n = get_parameter(M.size)[1]
+    fill!(X, 0)
+    view(X, diagind(X)) .= view(Xⁱ, 1:n) .* diag(p)
+    xi_ind = n + 1
+    for i in 1:n
+        for j in (i + 1):n
+            X[j, i] = Xⁱ[xi_ind]
+            xi_ind += 1
+        end
+    end
+    return X
+end
+
 @doc raw"""
     inner(M::CholeskySpace, p, X, Y)
 
 Compute the inner product on the [`CholeskySpace`](@ref) `M` at the
-lower triangular matric with positive diagonal `p` and the two tangent vectors
+lower triangular matrix with positive diagonal `p` and the two tangent vectors
 `X`,`Y`, i.e they are both lower triangular matrices with arbitrary diagonal.
 The formula reads
 
@@ -228,6 +255,38 @@ function parallel_transport_to!(::CholeskySpace, Y, p, X, q)
     return copyto!(Y, strictlyLowerTriangular(p) + Diagonal(diag(q) .* diag(X) ./ diag(p)))
 end
 
+function Random.rand!(
+    rng::AbstractRNG,
+    M::CholeskySpace,
+    pX;
+    vector_at=nothing,
+    σ::Real=one(eltype(pX)) /
+            (vector_at === nothing ? 1 : norm(convert(AbstractMatrix, vector_at))),
+    tangent_distr=:Gaussian,
+)
+    N = get_parameter(M.size)[1]
+    if vector_at === nothing
+        p_spd = rand(
+            rng,
+            SymmetricPositiveDefinite(N; parameter=:field);
+            σ=σ,
+            tangent_distr=tangent_distr,
+        )
+        pX .= cholesky(p_spd).L
+    else
+        p_spd = vector_at * vector_at'
+        X_spd = rand(
+            rng,
+            SymmetricPositiveDefinite(N; parameter=:field);
+            vector_at=p_spd,
+            σ=σ,
+            tangent_distr=tangent_distr,
+        )
+        pX .= spd_to_cholesky(p_spd, X_spd)[2]
+    end
+    return pX
+end
+
 @doc raw"""
     zero_vector(M::CholeskySpace, p)
 

src/manifolds/SymmetricPositiveDefinite.jl

@@ -236,6 +236,9 @@ function get_embedding(M::SymmetricPositiveDefinite{Tuple{Int}})
     return Euclidean(n, n; field=ℝ, parameter=:field)
 end
 
+get_parameter_type(::SymmetricPositiveDefinite{<:TypeParameter}) = :type
+get_parameter_type(::SymmetricPositiveDefinite{Tuple{Int}}) = :field
+
 @doc raw"""
     injectivity_radius(M::SymmetricPositiveDefinite[, p])
     injectivity_radius(M::MetricManifold{SymmetricPositiveDefinite,AffineInvariantMetric}[, p])

src/manifolds/SymmetricPositiveDefiniteLogCholesky.jl

@@ -31,13 +31,17 @@ d_{\mathcal P(n)}(p,q) = \sqrt{
  + \lVert \log(\operatorname{diag}(x)) - \log(\operatorname{diag}(y))\rVert_{\mathrm{F}}^2 }\ \ ,
 ````
 
-where ``x`` and ``y`` are the cholesky factors of ``p`` and ``q``, respectively,
+where ``x`` and ``y`` are the Cholesky factors of ``p`` and ``q``, respectively,
 ``⌊⋅⌋`` denbotes the strictly lower triangular matrix of its argument,
 and ``\lVert⋅\rVert_{\mathrm{F}}`` the Frobenius norm.
 """
 function distance(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, p, q)
     N = get_parameter(M.manifold.size)[1]
-    return distance(CholeskySpace(N), cholesky(p).L, cholesky(q).L)
+    return distance(
+        CholeskySpace(N; parameter=get_parameter_type(M.manifold)),
+        cholesky(p).L,
+        cholesky(q).L,
+    )
 end
 
 @doc raw"""
@@ -50,7 +54,7 @@ Compute the exponential map on the [`SymmetricPositiveDefinite`](@ref) `M` with
 \exp_p X = (\exp_y W)(\exp_y W)^\mathrm{T}
 ````
 
-where ``\exp_xW`` is the exponential map on [`CholeskySpace`](@ref), ``y`` is the cholesky
+where ``\exp_xW`` is the exponential map on [`CholeskySpace`](@ref), ``y`` is the Cholesky
 decomposition of ``p``, ``W = y(y^{-1}Xy^{-\mathrm{T}})_\frac{1}{2}``,
 and ``(⋅)_\frac{1}{2}``
 denotes the lower triangular matrix with the diagonal multiplied by ``\frac{1}{2}``.
@@ -60,7 +64,7 @@ exp(::MetricManifold{ℝ,SymmetricPositiveDefinite,LogCholeskyMetric}, ::Any...)
 function exp!(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, q, p, X)
     N = get_parameter(M.manifold.size)[1]
     (y, W) = spd_to_cholesky(p, X)
-    z = exp(CholeskySpace(N), y, W)
+    z = exp(CholeskySpace(N; parameter=get_parameter_type(M.manifold)), y, W)
     return copyto!(q, z * z')
 end
 function exp!(
@@ -73,6 +77,35 @@ function exp!(
     return exp!(M, q, p, t * X)
 end
 
+function get_coordinates_orthonormal!(
+    M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric},
+    Xⁱ,
+    p,
+    X,
+    rn::RealNumbers,
+)
+    N = get_parameter(M.manifold.size)[1]
+    MC = CholeskySpace(N; parameter=get_parameter_type(M.manifold))
+    (y, W) = spd_to_cholesky(p, X)
+    get_coordinates_orthonormal!(MC, Xⁱ, y, W, rn)
+    return Xⁱ
+end
+
+function get_vector_orthonormal!(
+    M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric},
+    X,
+    p,
+    Xⁱ,
+    rn::RealNumbers,
+)
+    N = get_parameter(M.manifold.size)[1]
+    MC = CholeskySpace(N; parameter=get_parameter_type(M.manifold))
+    y = cholesky(p).L
+    get_vector_orthonormal!(MC, X, y, Xⁱ, rn)
+    tangent_cholesky_to_tangent_spd!(p, X)
+    return X
+end
+
 @doc raw"""
     inner(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, p, X, Y)
 
@@ -85,15 +118,15 @@ a [`MetricManifold`](@ref) with [`LogCholeskyMetric`](@ref). The formula reads
 ````
 
 where ``⟨⋅,⋅⟩_x`` denotes inner product on the [`CholeskySpace`](@ref),
-``z`` is the cholesky factor of ``p``,
+``z`` is the Cholesky factor of ``p``,
 ``a_z(W) = z (z^{-1}Wz^{-\mathrm{T}})_{\frac{1}{2}}``, and ``(⋅)_\frac{1}{2}``
 denotes the lower triangular matrix with the diagonal multiplied by ``\frac{1}{2}``
 """
 function inner(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, p, X, Y)
     N = get_parameter(M.manifold.size)[1]
     (z, Xz) = spd_to_cholesky(p, X)
     (z, Yz) = spd_to_cholesky(p, z, Y)
-    return inner(CholeskySpace(N), z, Xz, Yz)
+    return inner(CholeskySpace(N; parameter=get_parameter_type(M.manifold)), z, Xz, Yz)
 end
 
 """
@@ -113,7 +146,7 @@ The formula can be adapted from the [`CholeskySpace`](@ref) as
 ````math
 \log_p q = xW^{\mathrm{T}} + Wx^{\mathrm{T}},
 ````
-where ``x`` is the cholesky factor of ``p`` and ``W=\log_x y`` for ``y`` the cholesky factor
+where ``x`` is the Cholesky factor of ``p`` and ``W=\log_x y`` for ``y`` the Cholesky factor
 of ``q`` and the just mentioned logarithmic map is the one on [`CholeskySpace`](@ref).
 """
 log(::MetricManifold{ℝ,SymmetricPositiveDefinite,LogCholeskyMetric}, ::Any...)
@@ -122,7 +155,7 @@ function log!(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetri
     N = get_parameter(M.manifold.size)[1]
     x = cholesky(p).L
     y = cholesky(q).L
-    log!(CholeskySpace(N), X, x, y)
+    log!(CholeskySpace(N; parameter=get_parameter_type(M.manifold)), X, x, y)
     return tangent_cholesky_to_tangent_spd!(x, X)
 end
 
@@ -138,7 +171,7 @@ end
 Parallel transport the tangent vector `X` at `p` along the geodesic to `q` with respect to
 the [`SymmetricPositiveDefinite`](@ref) manifold `M` and [`LogCholeskyMetric`](@ref).
 The parallel transport is based on the parallel transport on [`CholeskySpace`](@ref):
-Let ``x`` and ``y`` denote the cholesky factors of `p` and `q`, respectively and
+Let ``x`` and ``y`` denote the Cholesky factors of `p` and `q`, respectively and
 ``W = x(x^{-1}Xx^{-\mathrm{T}})_\frac{1}{2}``, where ``(⋅)_\frac{1}{2}`` denotes the lower
 triangular matrix with the diagonal multiplied by ``\frac{1}{2}``. With ``V`` the parallel
 transport on [`CholeskySpace`](@ref) from ``x`` to ``y``. The formula hear reads
@@ -164,6 +197,12 @@ function parallel_transport_to!(
     N = get_parameter(M.manifold.size)[1]
     y = cholesky(q).L
     (x, W) = spd_to_cholesky(p, X)
-    parallel_transport_to!(CholeskySpace(N), Y, x, W, y)
+    parallel_transport_to!(
+        CholeskySpace(N; parameter=get_parameter_type(M.manifold)),
+        Y,
+        x,
+        W,
+        y,
+    )
     return tangent_cholesky_to_tangent_spd!(y, Y)
 end

test/manifolds/cholesky_space.jl

@@ -33,6 +33,8 @@ include("../header.jl")
             test_vee_hat=false,
             exp_log_atol_multiplier=8.0,
             test_inplace=true,
+            test_rand_point=true,
+            test_rand_tvector=true,
         )
     end
     @testset "Test Error cases in is_point and is_vector" begin

test/manifolds/symmetric_positive_definite.jl

@@ -27,7 +27,7 @@ include("../header.jl")
     TEST_STATIC_SIZED && push!(types, MMatrix{3,3,Float64,9})
 
     for M in metrics
-        basis_types = if (M == M1 || M == M2)
+        basis_types = if (M == M1 || M == M2 || M == M3)
             (DefaultOrthonormalBasis(),)
         else
             ()
@@ -312,6 +312,7 @@ include("../header.jl")
         M = SymmetricPositiveDefinite(3; parameter=:field)
         @test typeof(get_embedding(M)) === Euclidean{Tuple{Int,Int},ℝ}
         @test repr(M) == "SymmetricPositiveDefinite(3; parameter=:field)"
+        @test Manifolds.get_parameter_type(M) === :field
     end
 
     @testset "Curvature" begin