Introduce the manifold of symmetric positive definite matrices with two metrics

docs/make.jl

@@ -11,7 +11,8 @@ makedocs(
             "Basic manifolds" => [
                 "Euclidean" => "manifolds/euclidean.md",
                 "Rotations" => "manifolds/rotations.md",
-                "Sphere" => "manifolds/sphere.md"
+                "Sphere" => "manifolds/sphere.md",
+                "Symmetric Positive Definite" => "manifolds/symmetricpositivedefinite.md"
             ],
             "Combined manifolds" => [
                 "Power manifold" => "manifolds/power.md",

docs/src/manifolds/euclidean.md

@@ -1,7 +1,16 @@
-# Sphere
+# Euclidean Space
+
+The Euclidean space $\mathbb R^n$ is a simple model space, since it has
+curvature constantly zero everywhere and hence nearly all operations simplify.
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["Euclidean.jl"]
+Order = [:type]
+```
 
 ```@autodocs
 Modules = [Manifolds]
 Pages = ["Euclidean.jl"]
-Order = [:type, :function]
+Order = [:function]
 ```

docs/src/manifolds/symmetricpositivedefinite.md

@@ -0,0 +1,76 @@
+# Symmetric Positive Definite Matrices
+
+The symmetric positive definite matrices
+
+```math
+\mathcal P(n) = \bigl\{ A \in \mathbb R^{n\times n}\ \big|\ 
+A = A^{\mathrm{T}} \text{ and }
+x^{\mathrm{T}}Ax > 0 \text{ for } 0\neq x \in\mathbb R^n \bigr\}
+```
+
+```@docs
+SymmetricPositiveDefinite
+```
+
+can -- for example -- be illustrated as ellipsoids:  since the eigen values 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](../assets/images/SPDSignal.png)
+
+The manifold can be equipped with different metrics
+
+## Common and Metric Independent functions
+```@docs
+injectivity_radius(::SymmetricPositiveDefinite{N},a::Vararg{Any,N} where N) where N
+is_manifold_point(::SymmetricPositiveDefinite{N},x; kwargs...) where N
+is_tangent_vector(::SymmetricPositiveDefinite{N},x,v; kwargs...) where N
+manifold_dimension(::SymmetricPositiveDefinite{N}) where N 
+representation_size(::SymmetricPositiveDefinite) 
+zero_tangent_vector(::SymmetricPositiveDefinite{N},x) where N
+zero_tangent_vector!(::SymmetricPositiveDefinite{N}, v, x) where N
+```
+
+## Default Metric
+```@docs
+distance(P::SymmetricPositiveDefinite{N},x,y) where N
+exp!(P::SymmetricPositiveDefinite{N}, y, x, v) where N
+inner(P::SymmetricPositiveDefinite{N}, x, w, v) where N
+log!(P::SymmetricPositiveDefinite{N}, v, x, y) where N
+tangent_orthonormal_basis(P::SymmetricPositiveDefinite{N},x,v) where N
+vector_transport_to!(P::SymmetricPositiveDefinite{N},vto, x, v, y, m::AbstractVectorTransportMethod) where N
+```
+
+## Linear Affine Metric
+
+```@docs
+LinearAffineMetric
+```
+
+This metric is also the default metric, i.e.
+any call of the following functions with
+`SymmetricPositiveDefinite(3)` will result in
+`MetricManifold(P,LinearAffineMetric())`and hence yield the formulae described in this seciton.
+
+```@docs
+distance(P::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric},x,y) where N
+exp!(P::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric}, y, x, v) where N
+inner(P::MetricManifold{SymmetricPositiveDefinite{N}, LinearAffineMetric}, x, w, v) where N
+log!(P::MetricManifold{SymmetricPositiveDefinite{N}, LinearAffineMetric}, v, x, y) where N
+tangent_orthonormal_basis(M::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric},x,v) where N
+vector_transport_to!(M::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric}, vto, x, v, y, ::ParallelTransport) where N
+```
+
+## Log Euclidean Metric
+
+```@docs
+LogEuclideanMetric
+```
+
+And we obtain the following functions
+
+```@docs
+distance(P::MetricManifold{SymmetricPositiveDefinite{N},LogEuclideanMetric},x,y) where N
+```
+
+### Literature

src/ArrayManifold.jl

@@ -130,19 +130,21 @@ function zero_tangent_vector(M::ArrayManifold, x; kwargs...)
     return w
 end
 
-function vector_transport_to(M::ArrayManifold, x, v, y)
+function vector_transport_to(M::ArrayManifold, x, v, y, m::AbstractVectorTransportMethod)
     return vector_transport_to(M.manifold,
                                array_value(x),
                                array_value(v),
-                               array_value(y))
+                               array_value(y),
+                               m)
 end
 
-function vector_transport_to!(M::ArrayManifold, vto, x, v, y)
+function vector_transport_to!(M::ArrayManifold, vto, x, v, y, m::AbstractVectorTransportMethod)
     return vector_transport_to!(M.manifold,
                                 array_value(vto),
                                 array_value(x),
                                 array_value(v),
-                                array_value(y))
+                                array_value(y),
+                                m)
 end
 
 function is_manifold_point(M::ArrayManifold, x::MPoint; kwargs...)

src/CholeskySpace.jl

@@ -0,0 +1,97 @@
+using LinearAlgebra: diagm, diag, eigen, eigvals, eigvecs, Symmetric, Diagonal, factorize, tr, norm, cholesky, LowerTriangular
+
+@doc doc"""
+    CholeskySpace{N} <: Manifold
+
+the manifold of lower triangular matrices with positive diagonal and
+a metric based on the cholesky decomposition. The formulae for this manifold
+are for example summarized in Table 1 of
+> Lin, Zenhua: Riemannian Geometry of Symmetric Positive Definite Matrices via
+> Cholesky Decomposition, arXiv: 1908.09326
+"""
+struct CholeskySpace{N} <: Manifold end
+CholeskySpace(n::Int) = CholeskySpace{N}()
+
+@generated representation_size(::CholeskySpace{N}) where N = (N,N)
+
+@generated manifold_dimension(::CholeskySpace{N}) where N = N*(N+1)/2
+
+@doc doc"""
+    inner(M,x,v,w)
+
+computes the inner product on the [`CholeskySpace`](@ref) `M` at the
+lower triangular matric with positive diagonal `x` and the two tangent vectors
+`v`,`w`, i.e they are both lower triangular matrices with arbitrary diagonal.
+The formula reads
+
+````math
+    g_{x}(v,w) = \sum_{i>j} v_{ij}w_{ij} + \sum_{j=1}^m v_{ii}w_{ii}x_{ii}^{-2}
+````
+"""
+inner(::CholeskySpace{N},x,v,w) where N = sum(LowerTriangular(v).*LowerTriangular(w)) + sum(diag(v).*diag(w)./( diag(x).^2 ))
+
+@doc doc"""
+    distance(M,x,y)
+
+computes the Riemannian distance on the [`CholeskySpace`](@ref) `M` between two
+matrices `x`, `y` that are lower triangular with positive diagonal. The formula
+reads
+
+````math
+d_{\mathcal M}(x,y) = \sqrt{
+\sum_{i>j} (x_{ij}-y_{ij})^2 + 
+\sum_{j=1}^m (\log x_{jj} - \log_{y_jj})^2
+}
+````
+"""
+distance(::CholeskySpace{N},x,y) where N = sqrt(
+  sum( (LowerTriangular(x) - LowerTriangular(y)).^2 ) + sum( (log.(diag(x)) - log.(diag(y))).^2 )
+)
+
+@doc doc"""
+    exp!(M,y,x,v)
+
+compute the exponential map on the [`CholeskySpace`](@ref) `M` eminating from
+the lower triangular matrix with positive diagonal `x` towards the lower triangular
+matrx `v` and return the result in `y`. The formula reads
+
+````math
+\exp_x v = \lfloor x \rfloor + \lfloor v \rfloor
++\operatorname{diag}(x)\operatorname{diag}(x)\exp{ \operatorname{diag}(v)\operatorname{diag}(x)^{-1}}
+````
+where $\lfloor x\rfloor$ denotes the lower triangular matrix of $x$ and
+$\opertorname{diag}(x)$ the diagonal matrix of $x$
+"""
+function exp!(::CholeskySpace{N},y,x,v) where N
+    y .= LowerTriangular(x) + LowerTriangular(v) + Diagonal(x)*Diagonal(exp.(diag(v)./diag(x)))
+    return y
+end
+@doc doc"""
+    log!(M,v,x,y)
+
+compute the exponential map on the [`CholeskySpace`](@ref) `M` eminating from
+the lower triangular matrix with positive diagonal `x` towards the lower triangular
+matrx `v` and return the result in `y`. The formula reads
+
+````math
+\exp_x v = \lfloor x \rfloor - \lfloor y \rfloor
++\operatorname{diag}(x)\log{ \operatorname{diag}(y)\operatorname{diag}(x)^{-1}}
+````
+where $\lfloor x\rfloor$ denotes the lower triangular matrix of $x$ and
+$\opertorname{diag}(x)$ the diagonal matrix of $x$
+"""
+function log!(::CholeskySpace{N},v,x,y) where N
+    v .= LowerTriangular(y) + LowerTriangular(x) + Diagonal(x)*Diagonal(log.(diag(y)./diag(x)))
+    return v
+end
+
+@doc doc"""
+    zero_tangent_vector!(M,v,x)
+
+returns the zero tangent vector on the [`CholeskySpace`](@ref) `M` at `x` in
+the variable `v`.
+"""
+function zero_tangent_vector!(M::CholeskySpace{N},v,x) where N
+    fill!(v,0)
+    return v
+end

src/Euclidean.jl

@@ -9,20 +9,38 @@ Euclidean vector space $\mathbb R^n$.
 
 generates the $n$-dimensional vector space $\mathbb R^n$.
 
-   Euclidean(m, n)
+    Euclidean(n₁,n₂,...,nᵢ)
 
-generates the $mn$-dimensional vector space $\mathbb R^{m \times n}$, whose
-elements are interpreted as $m \times n$ matrices.
+generates the $n_1n_2\cdot\ldots n_i$-dimensional vector space $\mathbb R^{n_1, n_2, \ldots, n_i}$, whose
+elements are interpreted as $n_1 \times,n_2\times\cdots\times n_i$ arrays, e.g. for
+two parameters as matrices.
 """
 struct Euclidean{T<:Tuple} <: Manifold where {T} end
+Euclidean(n::Vararg{Int,N}) where N = Euclidean{Tuple{n...}}()
 
-Euclidean(n::Int) = Euclidean{Tuple{n}}()
-Euclidean(m::Int, n::Int) = Euclidean{Tuple{m,n}}()
+"""
+    representation_size(M::Euclidean{T})
 
+returns the array dimensions required to represent an element on the
+[`Euclidean`](@ref) manifold `M`, i.e. the vector of all array dimensions.
+"""
 @generated representation_size(::Euclidean{T}) where {T} = Tuple(T.parameters...)
 
+"""
+    manifold_dimension(M::Euclidean{T})
+
+returns the manifold dimension of the [`Euclidean`](@ref) manifold `M`, i.e.
+the product of all array dimensions.
+"""
 @generated manifold_dimension(::Euclidean{T}) where {T} = *(T.parameters...)
 
+"""
+    EuclideanMetric <: RiemannianMetric
+
+a general type for any manifold that employs the Euclidean Metric, for example
+the [`Euclidean`](@ref) manifold itself, or the [`Sphere`](@ref), where every
+tangent space (as a plane in the embedding) uses this metric (in the embedding).
+"""
 struct EuclideanMetric <: RiemannianMetric end
 
 @traitimpl HasMetric{Euclidean,EuclideanMetric}
@@ -38,31 +56,104 @@ log_local_metric_density(M::MetricManifold{<:Manifold,EuclideanMetric}, x) = zer
 @inline inner(::Euclidean, x, v, w) = dot(v, w)
 @inline inner(::MetricManifold{<:Manifold,EuclideanMetric}, x, v, w) = dot(v, w)
 
+"""
+    distance(M::Euclidean,x,y)
+
+compute the Euclidean distance between two points on the [`Euclidean`](@ref)
+manifold `M`, i.e. for vectors it's just the norm of the difference, for matrices
+and higher order arrays, the matrix and ternsor Frobenius norm, respectively.
+"""
 distance(::Euclidean, x, y) = norm(x-y)
+"""
+    norm(M::Euclidean,x,v)
+
+compute the norm of a tangent vector `v` at `x` on the [`Euclidean`](@ref)
+manifold `M`, i.e. since every tangent space can be identified with `M` itself
+in this case, just the (Frobenius) norm of `v`.
+"""
 norm(::Euclidean, x, v) = norm(v)
 norm(::MetricManifold{<:Manifold,EuclideanMetric}, x, v) = norm(v)
 
+"""
+    exp!(M::Euclidean,y, x, v)
+
+compute the exponential map on the [`Euclidean`](@ref) manifold `M` from `x`
+in direction `v` in place, i.e. in `y`, which in this case is just
+````math
+    y = x + v
+````
+"""
 exp!(M::Euclidean, y, x, v) = (y .= x .+ v)
+"""
+    log!(M::Euclidean, v, x, y)
 
+compute the logarithmic map on the [`Euclidean`](@ref) manifold `M` from `x`
+tpo `y`, stored in the direction `v` in place, i.e. in `v`, which in this case is just
+````math
+    v = y - x
+````
+"""
 log!(M::Euclidean, v, x, y) = (v .= y .- x)
+"""
+    zero_tangent_vector!(M::Euclidean, v, x)
 
+compute a zero vector in the tangent space of `x` on the [`Euclidean`](@ref)
+manifold `M`, whcih here is just a zero filled array the same size as the 
+in place parameter `v` which is assumed to have the same `size` as `x`.    
+"""
 function zero_tangent_vector!(M::Euclidean, v, x)
     fill!(v, 0)
     return v
 end
 
+"""
+    project_point!(M::Euclidean, x)
+
+project an arbitrary point `x` onto the [`Euclidean`](@ref) manifold `M`, which
+is of course just the identity map.
+"""
 project_point!(M::Euclidean, x) = x
 
+"""
+    project_tangent!(M::Euclidean, w, x, v)
+
+project an arbitrary vector `v` into the tangent space of a point `x` on the
+[`Euclidean`](@ref) manifold `M`, which is just the identity, since any tangent
+space of `M` can be identified with all of `M`.
+"""
 function project_tangent!(M::Euclidean, w, x, v)
     w .= v
     return w
 end
+"""
+    vector_transport_to!(M::Euclidean, vto, x, v, y, ::ParallelTransport)
+
+parallel transport the vector `v` from the tangent space at `x` to the
+tangent space at `y` on the [`Euclidean`](@ref) manifold `M`,
+i.e. the in place `w` is just set to `v`.
+"""
+function vector_transport_to!(M::Euclidean, vto, x, v, y, ::ParallelTransport)
+    vto .= v
+    return vto
+end
+"""
+    flat!(M::Euclidean, v, x, w)
 
+since cotangent and tangent vectors can directly be identified in the [`Euclidean`](@ref)
+case, this yields just the identity for a tangent vector `w` in the tangent space
+of `x` on `M`. The result is returned also in place in `v`.
+"""
 function flat!(M::Euclidean, v::FVector{CotangentSpaceType}, x, w::FVector{TangentSpaceType})
     copyto!(v.data, w.data)
     return v
 end
+"""
+    sharp!(M::Euclidean, v, x, w)
 
+since cotangent and tangent vectors can directly be identified in the [`Euclidean`](@ref)
+case, this yields just the identity for a cotangent vector `w` in the tangent space
+of `x` on `M`. The result is returned also in place in `v`.
+"""
 function sharp!(M::Euclidean, v::FVector{TangentSpaceType}, x, w::FVector{CotangentSpaceType})
     copyto!(v.data, w.data)
     return v