Bases for CholeskySpace and LogCholesky metric

NEWS.md

@@ -5,6 +5,13 @@ 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
+
+### Added
+
+* Orthonormal bases for `CholeskySpace` and `LogCholesky` metric for `SymmetricPositiveDefinite`.
+* `rand` for `CholeskySpace`.
+
 ## [0.10.11] - 2025-01-02
 
 ### Added

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.10.11"
+version = "0.10.12"
 
 [deps]
 Einsum = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0"

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 @@
     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 @@
     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 @@
     return Euclidean(n, n; field=ℝ, parameter=:field)
 end
 
+get_parameter_arg(::SymmetricPositiveDefinite{<:TypeParameter}) = :type
+get_parameter_arg(::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_arg(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_arg(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_arg(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_arg(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_arg(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_arg(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_arg(M.manifold)),
+        Y,
+        x,
+        W,
+        y,
+    )
     return tangent_cholesky_to_tangent_spd!(y, Y)
 end

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
             ()