ο‚›GitHub

Symmetric positive definite matrices

Manifolds.SymmetricPositiveDefinite β€” Type
SymmetricPositiveDefinite{T} = HermitianPositiveDefinite{ℝ,T}

The manifold of symmetric positive definite matrices, i.e.

\[\mathcal P(n) = \bigl\{ p ∈ ℝ^{nΓ—n}\ \big|\ p = p^{\mathrm{T}} \text{ and } a^\mathrm{T}pa > 0 \text{ for all } a ∈ ℝ^{n}\backslash\{0\} \bigr\}\]

The tangent space at $p ∈ \mathcal P(n)$ reads

\[ T_p\mathcal P(n) = \bigl\{ X \in \mathbb R^{nΓ—n} \big|\ X=X^\mathrm{T} \bigr\},\]

i.e. the set of symmetric matrices.

This manifold is the real case $𝔽 = ℝ$ of the HermitianPositiveDefinite matrices, $\mathcal P(n) = \mathcal H(n)$.

Constructor

SymmetricPositiveDefinite(n; parameter::Symbol=:type)

generates the manifold $\mathcal P(n) \subset ℝ^{nΓ—n}$

source

This manifold can – for example – be illustrated as ellipsoids: since the eigenvalues are all positive they can be taken as lengths of the axes of an ellipsoids while the directions are given by the eigenvectors.

An example set of data

The manifold can be equipped with different metrics

Common and metric independent functions

The symmetric positive definite matrices are the real case of the Hermitian positive definite matrices, and all functions that do not depend on the field are documented there. The following keep their real docstrings.

Base.rand β€” Method
rand(M::SymmetricPositiveDefinite; vector_at=nothing, tangent_distr=:Gaussian, Οƒ::Real=1/sqrt(n))

Generate a random symmetric positive definite matrix on the SymmetricPositiveDefinite manifold M, or, if vector_at is a point $p$, a random tangent vector at $p$. The parameter Οƒ has no effect on a point.

For tangent_distr = :Gaussian the coordinates of the tangent vector in an orthonormal basis of $T_p\mathcal P(n)$ are independent and normally distributed with standard deviation Οƒ. For tangent_distr = :Rician the tangent vector is $RR^{\mathrm{T}}$ with $R = L + \sqrt{Οƒ}U$, where $L$ is the lower triangular Cholesky factor of $p$ and $U$ is upper triangular with independent standard normal entries; here Οƒ defaults to $1/\lVert p\rVert_{\mathrm{F}}$ instead.

source

Default metric: the affine invariant metric

Manifolds.AffineInvariantMetric β€” Type
AffineInvariantMetric <: AbstractMetric

The linear affine metric is the metric for Hermitian, and in particular symmetric, positive definite matrices, that employs matrix logarithms and exponentials, which yields a linear and affine metric.

source

This metric is also the default metric, i.e. any call of the following functions with P=SymmetricPositiveDefinite(3) will result in MetricManifold(P,AffineInvariantMetric())and hence yield the formulae described in this section. The same formulas hold for the Hermitian positive definite matrices.

ManifoldDiff.riemannian_Hessian β€” Method
riemannian_Hessian(M::SymmetricPositiveDefinite, p, G, H, X)

The Riemannian Hessian can be computed as stated in Eq. (7.3) [Ngu23]. Let $\nabla f(p)$ denote the Euclidean gradient G, $\nabla^2 f(p)[X]$ the Euclidean Hessian H, and $\operatorname{sym}(X) = \frac{1}{2}\bigl(X^{\mathrm{T}}+X\bigr)$ the symmetrization operator. Then the formula reads

\[ \operatorname{Hess}f(p)[X] = p\operatorname{sym}(βˆ‡^2 f(p)[X])p + \operatorname{sym}\bigl( X\operatorname{sym}\bigl(βˆ‡ f(p)\bigr)p)\]

source
Manifolds.volume_density β€” Method
volume_density(::SymmetricPositiveDefinite, p, X)

Compute the volume density of the SymmetricPositiveDefinite manifold at p in direction X. See [CKA17], Section 6.2 for details. Note that metric in Manifolds.jl has a different scaling factor than the reference. It is the determinant of the differential of $\exp_p$ at X in orthonormal bases. With the eigenvalues $Ξ»_1,…,Ξ»_n$ of $p^{-1/2}Xp^{-1/2}$ it reads

\[ΞΈ_p(X) = \prod_{1 ≀ i < j ≀ n} \frac{\sinh\bigl(\lvert Ξ»_i-Ξ»_j\rvert/2\bigr)}{\lvert Ξ»_i-Ξ»_j\rvert/2},\]

where a factor with $Ξ»_i = Ξ»_j$ is $1$, see [dSLCY25], Proposition 4.3.

source
ManifoldsBase.change_metric β€” Method
change_metric(M::SymmetricPositiveDefinite, E::EuclideanMetric, p, X)

Given a tangent vector $X ∈ T_p\mathcal P(n)$ with respect to the EuclideanMetric g_E, this function changes into the AffineInvariantMetric (default) metric on the SymmetricPositiveDefinite M.

To be precise we are looking for $c\colon T_p\mathcal P(n) β†’ T_p\mathcal P(n)$ such that for all $Y,Z ∈ T_p\mathcal P(n)$ it holds

\[⟨Y,Z⟩ = \operatorname{tr}(YZ) = \operatorname{tr}(p^{-1}c(Y)p^{-1}c(Z)) = g_p(c(Y),c(Z)).\]

Since $g_p(Y,Z) = \operatorname{tr}(p^{-1}Yp^{-1}Z) = \operatorname{tr}\bigl((p^{-\frac{1}{2}}Yp^{-\frac{1}{2}})(p^{-\frac{1}{2}}Zp^{-\frac{1}{2}})\bigr)$ by the cyclic property of the trace is the Frobenius inner product of $p^{-\frac{1}{2}}Yp^{-\frac{1}{2}}$ and $p^{-\frac{1}{2}}Zp^{-\frac{1}{2}}$, see [PFA06], Section 3.3, the map $c(X) = p^{\frac{1}{2}}Xp^{\frac{1}{2}}$ is computed, which is symmetric whenever $X$ is.

source
ManifoldsBase.change_representer β€” Method
change_representer(M::SymmetricPositiveDefinite, E::EuclideanMetric, p, X)

Compute the representer of the linear function given by $X ∈ T_p\mathcal M$ with respect to the (default) AffineInvariantMetric on the SymmetricPositiveDefinite M. Here X represents the linear function on the tangent space at p with respect to the EuclideanMetric g_E.

To be precise we are looking for $Z∈T_p\mathcal P(n)$ such that for all $Y∈T_p\mathcal P(n)$ it holds

\[⟨X,Y⟩ = \operatorname{tr}(XY) = \operatorname{tr}(p^{-1}Zp^{-1}Y) = g_p(Z,Y)\]

and hence $Z = pXp$.

source
ManifoldsBase.distance β€” Method
distance(M::SymmetricPositiveDefinite, p, q)
distance(M::MetricManifold{SymmetricPositiveDefinite,AffineInvariantMetric}, p, q)

Compute the distance on the SymmetricPositiveDefinite manifold between p and q, as a MetricManifold with AffineInvariantMetric. The formula reads

\[d_{\mathcal P(n)}(p,q) = \lVert \operatorname{Log}(p^{-\frac{1}{2}}qp^{-\frac{1}{2}})\rVert_{\mathrm{F}}.,\]

where $\operatorname{Log}$ denotes the matrix logarithm and $\lVertβ‹…\rVert_{\mathrm{F}}$ denotes the matrix Frobenius norm.

source
ManifoldsBase.get_basis β€” Method
[Ξ,κ] = get_basis(M::SymmetricPositiveDefinite, p, B::DefaultOrthonormalBasis)
[Ξ,κ] = get_basis(M::MetricManifold{<:SymmetricPositiveDefinite,AffineInvariantMetric}, p, B::DefaultOrthonormalBasis)

Return a default ONB for the tangent space $T_p\mathcal P(n)$ of the SymmetricPositiveDefinite with respect to the AffineInvariantMetric.

\[ g_p(X,Y) = \operatorname{tr}(p^{-1} X p^{-1} Y),\]

The basis constructed here is based on the ONB for symmetric matrices constructed as follows. Let

\[\Delta_{i,j} = (a_{k,l})_{k,l=1}^n \quad \text{ with } a_{k,l} = \begin{cases} 1 & \mbox{ for } k=l \text{ if } i=j\\ \frac{1}{\sqrt{2}} & \mbox{ for } k=i, l=j \text{ or } k=j, l=i\\ 0 & \text{ else.} \end{cases}\]

which forms an ONB for the space of symmetric matrices.

We then form the ONB by

\[ \Xi_{i,j} = p^{\frac{1}{2}}\Delta_{i,j}p^{\frac{1}{2}},\qquad i=1,\ldots,n, j=i,\ldots,n.\]

source
ManifoldsBase.get_basis_diagonalizing β€” Method
[Ξ,κ] = get_basis_diagonalizing(M::SymmetricPositiveDefinite, p, B::DiagonalizingOrthonormalBasis)
[Ξ,κ] = get_basis_diagonalizing(M::MetricManifold{<:SymmetricPositiveDefinite,AffineInvariantMetric}, p, B::DiagonalizingOrthonormalBasis)

Return an orthonormal basis Ξ of the tangent space at p that diagonalizes the curvature tensor $R(u,v)w$, together with its eigenvalues κ. The basis is a vector of tangent vectors (of length manifold_dimension of M) on the MetricManifold of the SymmetricPositiveDefinite manifold M with AffineInvariantMetric, and the direction B.frame_direction $V$ has curvature 0.

The construction is based on an ONB for the symmetric matrices similar to get_basis(::SymmetricPositiveDefinite, p, ::DefaultOrthonormalBasis just that the ONB here is build from the eigen vectors of $p^{\frac{1}{2}}Vp^{\frac{1}{2}}$.

source
ManifoldsBase.get_vector β€” Method
get_vector(::SymmetricPositiveDefinite, p, c, ::DefaultOrthonormalBasis)

Using the basis from get_basis the vector reconstruction with respect to this ONB can be simplified to

\[ X = p^{\frac{1}{2}} \Biggl( \sum_{i=1,j=i}^n c_k \Delta_{i,j} \Biggr) p^{\frac{1}{2}}\]

where $k$ is the linearized index of the $i=1,\ldots,n, j=i,\ldots,n$.

source
ManifoldsBase.parallel_transport_to β€” Method
parallel_transport_to(M::SymmetricPositiveDefinite, p, X, q)
parallel_transport_to(M::MetricManifold{SymmetricPositiveDefinite,AffineInvariantMetric}, p, X, y)

Compute the parallel transport of X from the tangent space at p to the tangent space at q on the SymmetricPositiveDefinite as a MetricManifold with the AffineInvariantMetric. The formula reads

\[\mathcal P_{q←p}X = p^{\frac{1}{2}} \operatorname{Exp}\bigl( \frac{1}{2}p^{-\frac{1}{2}}\log_p(q)p^{-\frac{1}{2}} \bigr) p^{-\frac{1}{2}}X p^{-\frac{1}{2}} \operatorname{Exp}\bigl( \frac{1}{2}p^{-\frac{1}{2}}\log_p(q)p^{-\frac{1}{2}} \bigr) p^{\frac{1}{2}},\]

where $\operatorname{Exp}$ denotes the matrix exponential and log the logarithmic map on SymmetricPositiveDefinite (again with respect to the AffineInvariantMetric).

source
ManifoldsBase.sectional_curvature_min β€” Method
sectional_curvature_min(M::SymmetricPositiveDefinite)

Return the minimum sectional curvature of the SymmetricPositiveDefinite manifold with the AffineInvariantMetric.

This is 0 for $1Γ—1$ matrices, where the manifold has dimension one, and $-1/2$ otherwise: the sectional curvatures of this metric are at most 0, see sectional_curvature_max, and their absolute values are bounded by $1/2$ and by no smaller constant, see [CB23], Section 2, comment 4 after Theorem 2.7, p. 1445.

source

Bures-Wasserstein metric

ManifoldsBase.change_representer β€” Method
change_representer(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,BuresWassersteinMetric}, E::EuclideanMetric, p, X)

Compute the representer of the linear function given by $X ∈ T_p\mathcal M$ with respect to the BuresWassersteinMetric on the SymmetricPositiveDefinite M. Here X represents the linear function on the tangent space at p with respect to the EuclideanMetric g_E.

To be precise we are looking for $Z∈T_p\mathcal P(n)$ such that for all $Y∈T_p\mathcal P(n)$ it holds

\[⟨X,Y⟩ = \operatorname{tr}(XY) = ⟨Z,Y⟩_{\mathrm{BW}}\]

for all $Y$ and hence we get $Z$= 2(A+A^{\mathrm{T}})$with$A=Xp``.

source

Generalized Bures-Wasserstein metric

Manifolds.GeneralizedBuresWassersteinMetric β€” Type
GeneralizedBuresWassersteinMetric{T<:AbstractMatrix} <: RiemannianMetric

The generalized Bures Wasserstein metric for symmetric positive definite matrices, see [HMJG21].

This metric internally stores the symmetric positive definite matrix $B$ to generalise the metric, which is called $M$ in the cited paper.

source
ManifoldsBase.change_representer β€” Method
change_representer(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric}, E::EuclideanMetric, p, X)

Compute the representer of the linear function given by $X ∈ T_p\mathcal M$ with respect to the GeneralizedBuresWassersteinMetric with matrix $B$ on the SymmetricPositiveDefinite M. Here X represents the linear function on the tangent space at p with respect to the EuclideanMetric g_E.

To be precise we are looking for $Z∈T_p\mathcal P(n)$ such that for all $Y∈T_p\mathcal P(n)$ it holds

\[⟨X,Y⟩ = \operatorname{tr}(XY) = ⟨Z,Y⟩_{\mathrm{BW}}\]

for all $Y$ and hence we get $Z = 2pXB + 2BXp$.

source
ManifoldsBase.inner β€” Method
inner(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,<:GeneralizedBuresWassersteinMetric}, p, X, Y)

Compute the inner product on the SymmetricPositiveDefinite manifold M with respect to the GeneralizedBuresWassersteinMetric with matrix $B$, given in Eq. (3) of [HMJG23] by

\[ ⟨X,Y⟩ = \frac{1}{2}\operatorname{tr}(L_{p,B}(X)Y)\]

where $q=L_{p,B}(X)$ denotes the generalized Lyapunov operator, i.e. it solves $pqB + Bqp = X$, as defined below Eq. (3) there.

source

Log-Euclidean metric

Manifolds.LogEuclideanMetric β€” Type
LogEuclideanMetric <: RiemannianMetric

The LogEuclidean Metric consists of the Euclidean metric applied to all elements after mapping them into the Lie Algebra, i.e. performing a matrix logarithm beforehand.

source

Log-Cholesky metric

Base.exp β€” Method
exp(M::MetricManifold{SymmetricPositiveDefinite,LogCholeskyMetric}, p, X)

Compute the exponential map on the SymmetricPositiveDefinite M with LogCholeskyMetric from p into direction X. The formula reads

\[\exp_p X = (\exp_y W)(\exp_y W)^\mathrm{T}\]

where $\exp_xW$ is the exponential map on CholeskySpace, $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}$.

source
ManifoldsBase.distance β€” Method
distance(M::MetricManifold{SymmetricPositiveDefinite,LogCholeskyMetric}, p, q)

Compute the distance on the manifold of SymmetricPositiveDefinite matrices, i.e. between two symmetric positive definite matrices p and q with respect to the LogCholeskyMetric. The formula reads

\[d_{\mathcal P(n)}(p,q) = \sqrt{ \lVert ⌊ x βŒ‹ - ⌊ y βŒ‹ \rVert_{\mathrm{F}}^2 + \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, $βŒŠβ‹…βŒ‹$ denotes the strictly lower triangular matrix of its argument, and $\lVertβ‹…\rVert_{\mathrm{F}}$ the Frobenius norm.

source
ManifoldsBase.inner β€” Method
inner(M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric}, p, X, Y)

Compute the inner product of two matrices X, Y in the tangent space of p on the SymmetricPositiveDefinite manifold M, as a MetricManifold with LogCholeskyMetric. The formula reads

\[ g_p(X,Y) = ⟨a_z(X),a_z(Y)⟩_z,\]

where $βŸ¨β‹…,β‹…βŸ©_x$ denotes inner product on the CholeskySpace, $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}$

source
ManifoldsBase.parallel_transport_to β€” Method
vector_transport_to(
    M::MetricManifold{ℝ,<:SymmetricPositiveDefinite,LogCholeskyMetric},
    p,
    X,
    q,
    ::ParallelTransport,
)

Parallel transport the tangent vector X at p along the geodesic to q with respect to the SymmetricPositiveDefinite manifold M and LogCholeskyMetric. The parallel transport is based on the parallel transport on CholeskySpace: 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 of $W$ on CholeskySpace from $x$ to $y$, the formula here reads

\[\mathcal P_{q←p}X = yV^{\mathrm{T}} + Vy^{\mathrm{T}}.\]

source

Efficient representation

When a point p is used in several occasions, it might be beneficial to store the eigenvalues and vectors of p and optionally its square root and the inverse of the square root. The MatrixSqrtManifoldPoint can be used for exactly that; SPDPoint is deprecated in its favour.

Literature

[CKA17]
E.Β Chevallier, E.Β Kalunga and J.Β Angulo. Kernel Density Estimation on Spaces of Gaussian Distributions and Symmetric Positive Definite Matrices. SIAMΒ JournalΒ onΒ ImagingΒ Sciences 10, 191–215 (2017).
[CB23]
C.Β Criscitiello and N.Β Boumal. An Accelerated First-Order Method for Non-convex Optimization on Manifolds. FoundationsΒ ofΒ ComputationalΒ Mathematics 23, 1433–1509 (2023).
[HMJG21]
A.Β Han, B.Β Mishra, P.Β Jawanpuria and J.Β Gao. Learning with symmetric positive definite matrices via generalized Bures-Wasserstein geometry, arXivΒ Preprint (2021), arXiv:2110.10464.
[HMJG23]
A.Β Han, B.Β Mishra, P.Β Jawanpuria and J.Β Gao. Learning with Symmetric Positive Definite Matrices via Generalized Bures-Wasserstein Geometry. In: Geometric Science of Information, Vol.Β 14071 of Lecture Notes in Computer Science, edited by F.Β Nielsen and F.Β Barbaresco (Springer Nature Switzerland, 2023); pp.Β 405–415.
[HZ23]
Z.Β Huang and L.Β Zheng. Target detection based on generalized Bures–Wasserstein distance. EURASIPΒ JournalΒ onΒ AdvancesΒ inΒ SignalΒ Processing 2023 (2023).
[Lin19]
Z.Β Lin. Riemannian Geometry of Symmetric Positive Definite Matrices via Cholesky Decomposition. SIAMΒ JournalΒ onΒ MatrixΒ AnalysisΒ andΒ Applications 40, 1353–1370 (2019), arXiv:1908.09326.
[MMP18]
L.Β MalagΓ³, L.Β Montrucchio and G.Β Pistone. Wasserstein Riemannian geometry of Gaussian densities. InformationΒ Geometry 1, 137–179 (2018).
[Ngu23]
D.Β Nguyen. Operator-Valued Formulas for Riemannian Gradient and Hessian and Families of Tractable Metrics in Riemannian Optimization. JournalΒ ofΒ OptimizationΒ TheoryΒ andΒ Applications 198, 135–164 (2023), arXiv:2009.10159.
[PFA06]
X.Β Pennec, P.Β Fillard and N.Β Ayache. A Riemannian Framework for Tensor Computing. InternationalΒ JournalΒ ofΒ ComputerΒ Vision 66, 41–66 (2006).
[Ren11]
Q.Β Rentmeesters. A gradient method for geodesic data fitting on some symmetric Riemannian manifolds. In: IEEE Conference on Decision and Control and European Control Conference (2011); pp.Β 7141–7146.
[dSLCY25]
T.Β deΒ Surrel, F.Β Lotte, S.Β Chevallier and F.Β Yger. Wrapped Gaussian on the manifold of Symmetric Positive Definite Matrices. In: Proceedings of the 42nd International Conference on Machine Learning, Vol.Β 267 of Proceedings of Machine Learning Research (PMLR, 2025); pp.Β 12816–12839.
[TP23]
Y.Β Thanwerdas and X.Β Pennec. O(n)-invariant Riemannian metrics on SPD matrices. LinearΒ AlgebraΒ andΒ itsΒ Applications 661, 163–201 (2023).