Symmetric positive definite matrices
Manifolds.SymmetricPositiveDefinite
β TypeSymmetricPositiveDefinite{N} <: AbstractDecoratorManifold{π½}
The manifold of symmetric positive definite matrices, i.e.
\[\mathcal P(n) = \bigl\{ p β β^{n Γ n}\ \big|\ a^\mathrm{T}pa > 0 \text{ for all } a β β^{n}\backslash\{0\} \bigr\}\]
The tangent space at $T_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,
Constructor
SymmetricPositiveDefinite(n)
generates the manifold $\mathcal P(n) \subset β^{n Γ n}$
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.
The manifold can be equipped with different metrics
Common and metric independent functions
Base.rand
β Methodrand(M::SymmetricPositiveDefinite; Ο::Real=1)
Generate a random symmetric positive definite matrix on the SymmetricPositiveDefinite
manifold M
.
ManifoldsBase.check_point
β Methodcheck_point(M::SymmetricPositiveDefinite, p; kwargs...)
checks, whether p
is a valid point on the SymmetricPositiveDefinite
M
, i.e. is a matrix of size (N,N)
, symmetric and positive definite. The tolerance for the second to last test can be set using the kwargs...
.
ManifoldsBase.check_vector
β Methodcheck_vector(M::SymmetricPositiveDefinite, p, X; kwargs... )
Check whether X
is a tangent vector to p
on the SymmetricPositiveDefinite
M
, i.e. atfer check_point
(M,p)
, X
has to be of same dimension as p
and a symmetric matrix, i.e. this stores tangent vetors as elements of the corresponding Lie group. The tolerance for the last test can be set using the kwargs...
.
ManifoldsBase.injectivity_radius
β Methodinjectivity_radius(M::SymmetricPositiveDefinite[, p])
injectivity_radius(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}[, p])
injectivity_radius(M::MetricManifold{SymmetricPositiveDefinite,LogCholeskyMetric}[, p])
Return the injectivity radius of the SymmetricPositiveDefinite
. Since M
is a Hadamard manifold with respect to the LinearAffineMetric
and the LogCholeskyMetric
, the injectivity radius is globally $β$.
ManifoldsBase.manifold_dimension
β Methodmanifold_dimension(M::SymmetricPositiveDefinite)
returns the dimension of SymmetricPositiveDefinite
M
$=\mathcal P(n), n β β$, i.e.
\[\dim \mathcal P(n) = \frac{n(n+1)}{2}.\]
ManifoldsBase.project
β Methodproject(M::SymmetricPositiveDefinite, p, X)
project a matrix from the embedding onto the tangent space $T_p\mathcal P(n)$ of the SymmetricPositiveDefinite
matrices, i.e. the set of symmetric matrices.
ManifoldsBase.representation_size
β Methodrepresentation_size(M::SymmetricPositiveDefinite)
Return the size of an array representing an element on the SymmetricPositiveDefinite
manifold M
, i.e. $n Γ n$, the size of such a symmetric positive definite matrix on $\mathcal M = \mathcal P(n)$.
ManifoldsBase.zero_vector
β Methodzero_vector(M::SymmetricPositiveDefinite,x)
returns the zero tangent vector in the tangent space of the symmetric positive definite matrix x
on the SymmetricPositiveDefinite
manifold M
.
Default metric: the linear affine metric
Manifolds.LinearAffineMetric
β TypeLinearAffineMetric <: AbstractMetric
The linear affine metric is the metric for symmetric positive definite matrices, that employs matrix logarithms and exponentials, which yields a linear and affine metric.
This metric is also the default metric, i.e. any call of the following functions with P=SymmetricPositiveDefinite(3)
will result in MetricManifold(P,LinearAffineMetric())
and hence yield the formulae described in this seciton.
Base.exp
β Methodexp(M::SymmetricPositiveDefinite, p, X)
exp(M::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric}, p, X)
Compute the exponential map from p
with tangent vector X
on the SymmetricPositiveDefinite
M
with its default MetricManifold
having the LinearAffineMetric
. The formula reads
\[\exp_p X = p^{\frac{1}{2}}\operatorname{Exp}(p^{-\frac{1}{2}} X p^{-\frac{1}{2}})p^{\frac{1}{2}},\]
where $\operatorname{Exp}$ denotes to the matrix exponential.
Base.log
β Methodlog(M::SymmetricPositiveDefinite, p, q)
log(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}, p, q)
Compute the logarithmic map from p
to q
on the SymmetricPositiveDefinite
as a MetricManifold
with LinearAffineMetric
. The formula reads
\[\log_p q = p^{\frac{1}{2}}\operatorname{Log}(p^{-\frac{1}{2}}qp^{-\frac{1}{2}})p^{\frac{1}{2}},\]
where $\operatorname{Log}$ denotes to the matrix logarithm.
Manifolds.change_metric
β Methodchange_metric(M::SymmetricPositiveDefinite{n}, 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 LinearAffineMetric
(default) metric on the SymmetricPositiveDefinite
M
.
To be precise we are looking for $c\colon T_p\mathcal P(n) \to 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(Z),c(Y))\]
and hence $c(X) = pX$ is computed.
Manifolds.change_representer
β Methodchange_representer(M::SymmetricPositiveDefinite, E::EuclideanMetric, p, X)
Given a tangent vector $X β T_p\mathcal M$ representing a linear function on the tangent space at p
with respect to the EuclideanMetric
g_E
, this is turned into the representer with respect to the (default) metric, the LinearAffineMetric
on the SymmetricPositiveDefinite
M
.
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$.
ManifoldsBase.distance
β Methoddistance(M::SymmetricPositiveDefinite, p, q)
distance(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}, p, q)
Compute the distance on the SymmetricPositiveDefinite
manifold between p
and q
, as a MetricManifold
with LinearAffineMetric
. 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\cdot\rVert_{\mathrm{F}}$ denotes the matrix Frobenius norm.
ManifoldsBase.get_basis
β Method[Ξ,ΞΊ] = get_basis(M::SymmetricPositiveDefinite, p, B::DefaultOrthonormalBasis)
[Ξ,ΞΊ] = get_basis(M::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric}, p, B::DefaultOrthonormalBasis)
Return a default ONB for the tangent space $T_p\mathcal P(n)$ of the SymmetricPositiveDefinite
with respect to the LinearAffineMetric
.
\[ 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.\]
ManifoldsBase.get_basis_diagonalizing
β Method[Ξ,ΞΊ] = get_basis_diagonalizing(M::SymmetricPositiveDefinite, p, B::DiagonalizingOrthonormalBasis)
[Ξ,ΞΊ] = get_basis_diagonalizing(M::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric}, p, B::DiagonalizingOrthonormalBasis)
Return a orthonormal basis Ξ
as a vector of tangent vectors (of length manifold_dimension
of M
) in the tangent space of p
on the MetricManifold
of SymmetricPositiveDefinite
manifold M
with LinearAffineMetric
that diagonalizes the curvature tensor $R(u,v)w$ with eigenvalues ΞΊ
and where 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}}$.
ManifoldsBase.get_coordinates
β Methodget_coordinates(::SymmetricPositiveDefinite, p, X, ::DefaultOrthonormalBasis)
Using the basis from get_basis
the coordinates with respect to this ONB can be simplified to
\[ c_k = \mathrm{tr}(p^{-\frac{1}{2}}\Delta_{i,j} X)\]
where $k$ is trhe linearized index of the $i=1,\ldots,n, j=i,\ldots,n$.
ManifoldsBase.get_vector
β Methodget_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$.
ManifoldsBase.inner
β Methodinner(M::SymmetricPositiveDefinite, p, X, Y)
inner(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}, p, X, Y)
Compute the inner product of X
, Y
in the tangent space of p
on the SymmetricPositiveDefinite
manifold M
, as a MetricManifold
with LinearAffineMetric
. The formula reads
\[g_p(X,Y) = \operatorname{tr}(p^{-1} X p^{-1} Y),\]
ManifoldsBase.parallel_transport_to
β Methodparallel_transport_to(M::SymmetricPositiveDefinite, p, X, q)
parallel_transport_to(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}, 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 LinearAffineMetric
. 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 LinearAffineMetric
).
Bures-Wasserstein metric
Manifolds.BuresWassersteinMetric
β TypeBurresWassertseinMetric <: AbstractMetric
The Bures Wasserstein metric for symmetric positive definite matrices[MalagoMontruccioPistone2018].
Base.exp
β Methodexp(::MatricManifold{β,SymmetricPositiveDefinite,BuresWassersteinMetric}, p, X)
Compute the exponential map on SymmetricPositiveDefinite
with respect to the BuresWassersteinMetric
given by
\[ \exp_p(X) = p+X+L_p(X)pL_p(X)\]
where $q=L_p(X)$ denotes the Lyapunov operator, i.e. it solves $pq + qp = X$.
Base.log
β Methodlog(::MatricManifold{SymmetricPositiveDefinite,BuresWassersteinMetric}, p, q)
Compute the logarithmic map on SymmetricPositiveDefinite
with respect to the BuresWassersteinMetric
given by
\[ \log_p(q) = (pq)^{\frac{1}{2}} + (qp)^{\frac{1}{2}} - 2 p\]
where $q=L_p(X)$ denotes the Lyapunov operator, i.e. it solves $pq + qp = X$.
Manifolds.change_representer
β Methodchange_representer(M::MetricManifold{β,SymmetricPositiveDefinite,BuresWassersteinMetric}, E::EuclideanMetric, p, X)
Given a tangent vector $X β T_p\mathcal M$ representing a linear function on the tangent space at p
with respect to the EuclideanMetric
g_E
, this is turned into the representer with respect to the (default) metric, the BuresWassersteinMetric
on the SymmetricPositiveDefinite
M
.
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``.
ManifoldsBase.distance
β Methoddistance(::MatricManifold{SymmetricPositiveDefinite,BuresWassersteinMetric}, p, q)
Compute the distance with respect to the BuresWassersteinMetric
on SymmetricPositiveDefinite
matrices, i.e.
\[d(p,q) = \operatorname{tr}(p) + \operatorname{tr}(q) - 2\operatorname{tr}\Bigl( (p^{\frac{1}{2}}qp^{\frac{1}{2}} \bigr)^\frac{1}{2} \Bigr),\]
where the last trace can be simplified (by rotating the matrix products in the trace) to $\operatorname{tr}(pq)$.
ManifoldsBase.inner
β Methodinner(::MetricManifold{β,SymmetricPositiveDefinite,BuresWassersteinMetric}, p, X, Y)
Compute the inner product SymmetricPositiveDefinite
with respect to the BuresWassersteinMetric
given by
\[ β¨X,Yβ© = \frac{1}{2}\operatorname{tr}(L_p(X)Y)\]
where $q=L_p(X)$ denotes the Lyapunov operator, i.e. it solves $pq + qp = X$.
Generalized Bures-Wasserstein metric
Manifolds.GeneralizedBuresWassersteinMetric
β TypeGeneralizedBurresWassertseinMetric{T<:AbstractMatrix} <: AbstractMetric
The generalized Bures Wasserstein metric for symmetric positive definite matrices, see[HanMishraJawanpuriaGao2021].
This metric internally stores the symmetric positive definite matrix $M$ to generalise the metric, where the name also follows the mentioned preprint.
Base.exp
β Methodexp(::MatricManifold{β,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, X)
Compute the exponential map on SymmetricPositiveDefinite
with respect to the GeneralizedBuresWassersteinMetric
given by
\[ \exp_p(X) = p+X+\mathcal ML_{p,M}(X)pML_{p,M}(X)\]
where $q=L_{M,p}(X)$ denotes the generalized Lyapunov operator, i.e. it solves $pqM + Mqp = X$.
Base.log
β Methodlog(::MatricManifold{SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, q)
Compute the logarithmic map on SymmetricPositiveDefinite
with respect to the BuresWassersteinMetric
given by
\[ \log_p(q) = M(M^{-1}pM^{-1}q)^{\frac{1}{2}} + (qM^{-1}pM^{-1})^{\frac{1}{2}}M - 2 p.\]
Manifolds.change_representer
β Methodchange_representer(M::MetricManifold{β,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, E::EuclideanMetric, p, X)
Given a tangent vector $X β T_p\mathcal M$ representing a linear function on the tangent space at p
with respect to the EuclideanMetric
g_E
, this is turned into the representer with respect to the (default) metric, the GeneralizedBuresWassersteinMetric
on the SymmetricPositiveDefinite
M
.
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 = 2pXM + 2MXp$.
ManifoldsBase.distance
β Methoddistance(::MatricManifold{SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, q)
Compute the distance with respect to the BuresWassersteinMetric
on SymmetricPositiveDefinite
matrices, i.e.
\[d(p,q) = \operatorname{tr}(M^{-1}p) + \operatorname{tr}(M^{-1}q) - 2\operatorname{tr}\bigl( (p^{\frac{1}{2}}M^{-1}qM^{-1}p^{\frac{1}{2}} \bigr)^{\frac{1}{2}},\]
ManifoldsBase.inner
β Methodinner(::MetricManifold{β,SymmetricPositiveDefinite,GeneralizedBuresWassersteinMetric}, p, X, Y)
Compute the inner product SymmetricPositiveDefinite
with respect to the GeneralizedBuresWassersteinMetric
given by
\[ β¨X,Yβ© = \frac{1}{2}\operatorname{tr}(L_{p,M}(X)Y)\]
where $q=L_{M,p}(X)$ denotes the generalized Lyapunov operator, i.e. it solves $pqM + Mqp = X$.
Log-Euclidean metric
Manifolds.LogEuclideanMetric
β TypeLogEuclideanMetric <: 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.
ManifoldsBase.distance
β Methoddistance(M::MetricManifold{SymmetricPositiveDefinite{N},LogEuclideanMetric}, p, q)
Compute the distance on the SymmetricPositiveDefinite
manifold between p
and q
as a MetricManifold
with LogEuclideanMetric
. The formula reads
\[ d_{\mathcal P(n)}(p,q) = \lVert \operatorname{Log} p - \operatorname{Log} q \rVert_{\mathrm{F}}\]
where $\operatorname{Log}$ denotes the matrix logarithm and $\lVert\cdot\rVert_{\mathrm{F}}$ denotes the matrix Frobenius norm.
Log-Cholesky metric
Manifolds.LogCholeskyMetric
β TypeLogCholeskyMetric <: RiemannianMetric
The Log-Cholesky metric imposes a metric based on the Cholesky decomposition as introduced by [Lin2019].
Base.exp
β Methodexp(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 $(\cdot)_\frac{1}{2}$ denotes the lower triangular matrix with the diagonal multiplied by $\frac{1}{2}$.
Base.log
β Methodlog(M::MetricManifold{SymmetricPositiveDefinite,LogCholeskyMetric}, p, q)
Compute the logarithmic map on SymmetricPositiveDefinite
M
with respect to the LogCholeskyMetric
emanating from p
to q
. The formula can be adapted from the CholeskySpace
as
\[\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 of $q$ and the just mentioned logarithmic map is the one on CholeskySpace
.
ManifoldsBase.distance
β Methoddistance(M::MetricManifold{SymmetricPositiveDefinite,LogCholeskyMetric}, p, q)
Compute the distance on the manifold of SymmetricPositiveDefinite
nmatrices, 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, $β\cdotβ$ denbotes the strictly lower triangular matrix of its argument, and $\lVert\cdot\rVert_{\mathrm{F}}$ the Frobenius norm.
ManifoldsBase.inner
β Methodinner(M::MetricManifold{LogCholeskyMetric,β,SymmetricPositiveDefinite}, 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 $β¨\cdot,\cdotβ©_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 $(\cdot)_\frac{1}{2}$ denotes the lower triangular matrix with the diagonal multiplied by $\frac{1}{2}$
ManifoldsBase.parallel_transport_to
β Methodvector_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 $(\cdot)_\frac{1}{2}$ denotes the lower triangular matrix with the diagonal multiplied by $\frac{1}{2}$. With $V$ the parallel transport on CholeskySpace
from $x$ to $y$. The formula hear reads
\[\mathcal P_{qβp}X = yV^{\mathrm{T}} + Vy^{\mathrm{T}}.\]
Statistics
Statistics.mean
β Methodmean(
M::SymmetricPositiveDefinite,
x::AbstractVector,
[w::AbstractWeights,]
method = GeodesicInterpolation();
kwargs...,
)
Compute the Riemannian mean
of x
using GeodesicInterpolation
.
Literature
- MalagoMontruccioPistone2018
MalagΓ², L., Montrucchio, L., Pistone, G.: Wasserstein Riemannian geometry of Gaussian densities. Information Geometry, 1, pp. 137β179, 2018. doi: 10.1007/s41884-018-0014-4
- HanMishraJawanpuriaGao2021
Han, A., Mishra, B., Jawanpuria, P., Gao, J.: Generalized Bures-Wasserstein geometry for positive definite matrices. arXiv: 2110.10464.
- Lin2019
Lin, Zenhua: "Riemannian Geometry of Symmetric Positive Definite Matrices via Cholesky Decomposition", arXiv: 1908.09326.