Introduce the manifold of symmetric positive definite matrices with two metrics

Project.toml

@@ -14,6 +14,7 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SimpleTraits = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+UnsafeArrays = "c4a57d5a-5b31-53a6-b365-19f8c011fbd6"
 
 [extras]
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
@@ -22,4 +23,4 @@ ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["Test", "DoubleFloats", "ForwardDiff", "ReverseDiff"]

src/Manifolds.jl

@@ -256,7 +256,7 @@ end
 
 Inner product of tangent vectors `v` and `w` at point `x` from manifold `M`.
 """
-inner(M::Manifold, x, v, w) = error("inner: Inner product not implemented on a $(typeof(M)) for input point $(typeof(x)) and tangent vectors $(typeof(v)) and $(typeof(w)).")
+inner(M::Manifold, x, v, w) = error("Inner product not implemented on a $(typeof(M)) for input point $(typeof(x)) and tangent vectors $(typeof(v)) and $(typeof(w)).")
 
 """
     norm(M::Manifold, x, v)
@@ -372,11 +372,20 @@ function shortest_geodesic(M::Manifold, x, y, T::AbstractVector)
     return geodesic(M, x, log(M, x, y), T)
 end
 
-vector_transport!(M::Manifold, vto, x, v, y) = project_tangent!(M, vto, x, v)
+abstract type VectorTransportMethod end
 
-function vector_transport(M::Manifold, x, v, y)
+struct ParallelTransport <: VectorTransportMethod end
+struct VectorProjection <: VectorTransportMethod end
+
+vector_transport!(M::Manifold, vto, x, v, y) = vector_transport!(M,vto,x,v,y,::ParallelTransport)
+
+vector_transport!(M::Manifold, vto, x, v, y, method::VectorTransportMethod) = error("No vector transport method $(typeof(method)) on $(typeof(M)) for two points $(typeof(x)) and $(typeof(y)) and a tangent vector $(typeof(v)).")
+
+vector_transport!(M::Manifold, vto, x, v, y, ::VectorProjection) = project_tangent!(M, vto, x, v)
+
+function vector_transport(M::Manifold, x, v, y, method::VectorTransportMethod=ParallelTransport)
     vto = similar_result(M, vector_transport, v, x, y)
-    vector_transport!(M, vto, x, y, v)
+    vector_transport!(M, vto, x, y, v, method)
     return vto
 end
 
@@ -465,6 +474,7 @@ include("Metric.jl")
 include("Euclidean.jl")
 include("Rotations.jl")
 include("Sphere.jl")
+include("SymmetricPositiveDefinite.jl")
 include("ProjectedDistribution.jl")
 
 export Manifold,

src/SymmetricPositiveDefinite.jl

@@ -0,0 +1,176 @@
+using LinearAlgebra: svd, eig, eigen
+
+@doc doc"""
+    SymmetricPositiveDefinite{N} <: Manifold
+
+The manifold of symmetric positive definite matrices, i.e.
+
+```math
+\mathcal P(n) \coloneqq
+\bigl\{
+    x \in \mathbb R^{n\\times n} :
+    \xi^\mathrm{T}x\xi > 0 \text{ for all } \xi \in \mathbb R^{n}\backslash\{0\}
+\bigr\}
+```
+
+# Constructor
+
+    SymmetricPositiveDefinite(n)
+
+generates the $\mathcal P(n) \subset \mathbb R^{n\times n}$
+"""
+struct SymmetricPositiveDefinite{N} <: Manifold end
+SymmetricPositiveDefinite(n::Int) = SymmetricPositiveDefinite{n}()
+
+@doc doc"""
+    manifold_dimension(::SymmetricPositiveDefinite{N})
+
+returns the dimension of the manifold [`SymmetricPositiveDefinite`](@ref) $\mathcal P(n)$, N\in \mathbb N$, i.e.
+```math
+    \frac{n(n+1)}{2}    
+```
+"""
+@generated manifold_dimension(::SymmetricPositiveDefinite{N}) where {N} = N*(N+1)/2
+
+@doc doc"""
+
+"""
+struct LinearAffineMetric <: Metric end
+@traitimpl HasMetric{SymmetricPositiveDefinite, LogEuclidean}
+
+@doc doc"""
+
+"""
+struct LogEuclideanMetric <: Metric end
+@traitimpl HasMetric{SymmetricPositiveDefinite, LogEuclidean}
+
+distance(P::SymmetricPositiveDefinite{N},x,y) = distance(MetricManifold(S,LinearAffineMetric)},x,y)
+function distance(P::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric},x,y)
+    s = real.( eigen( x,y ).values )
+    return any(s .<= eps() ) ? 0 : sqrt(  sum( abs.(log.(s)).^2 )  )
+end
+
+inner(P::SymmetricPositiveDefinite{N}, x, w, v) = inner(MetricManifold(S,LinearAffineMetric),x,w,v)
+function inner(P::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric},x,w,v)
+	svd1 = svd(x)
+	U = svd1.U
+	S = svd1.S
+	SInv = Diagonal( 1 ./ S )
+	return tr( w * U * SInv * transpose(U) * v * U * SInv * transpose(U) )
+end
+
+norm(P::SymmetricPositiveDefinite{N},x,v) = sqrt( inner(P,x,v,v) )
+norm(P::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric},x,v) = sqrt( inner(P,x,v,v) )
+
+function exp!(P::SymmetricPositiveDefinite{N},y,x,v) = exp!(MetricManifold(SymmetricPositiveDefinite,LinearAffineMetric),y,x,v)
+function exp!(P::Metricmanifold{SymmetricPositiveDefinite{N},LinearAffineMetric}, y, x, v)
+	svd1 = svd(x)
+	U = svd1.U
+	S = svd1.S
+    Ssqrt = Diagonal( sqrt.(S) )
+    SsqrtInv = Diagonal( 1 ./ sqrt.(S) )
+    xSqrt = U*Ssqrt*transpose(U);
+    xSqrtInv = U*SsqrtInv*transpose(U)
+    T = xSqrtInv * (t.*ξ.value) * xSqrtInv
+    eig1 = eigen(0.5*( T + transpose(T) ) ) # numerical stabilization
+   	Se = Diagonal( exp.(eig1.values) )
+    Ue = eig1.vectors
+    y = xSqrt*Ue*Se*transpose(Ue)*xSqrt
+    y = 0.5*( y + transpose(y) ) ) # numerical stabilization
+    return y
+end
+
+function log!(P::SymmetricPositiveDefinite{N}, v, x, y) = exp!(MetricManifold(SymmetricPositiveDefinite,LinearAffineMetric),y,x,v)
+function log!(P::Metricmanifold{SymmetricPositiveDefinite{N},LinearAffineMetric}, v, x, y)
+    svd1 = svd( x )
+	U = svd1.U
+	S = svd1.S
+    Ssqrt = Diagonal( sqrt.(S) )
+    SsqrtInv = Diagonal( 1 ./ sqrt.(S) )
+    xSqrt = U*Ssqrt*transpose(U)
+    xSqrtInv = U*SsqrtInv*transpose(U)
+    T = xSqrtInv * getValue(y) * xSqrtInv
+	svd2 = svd(0.5*(T+transpose(T)))
+	Se = Diagonal( log.(max.(svd2.S,eps()) ) )
+	Ue = svd2.U
+	v = xSqrt * Ue*Se*transpose(Ue) * xSqrt
+	v = 0.5*( v + transpose(v) ) )
+    return v
+end
+
+function vector_transport!(M::SymmetricPositiveDefinite,vto, x, v, y, ::ParallelTransport)
+  if norm(x-y)<1e-13
+    vto = v
+    return vto
+  end
+  svd1 = svd(x)
+  U = svd1.U
+  S = svd1.S
+  Ssqrt = sqrt.(S)
+  SsqrtInv = Diagonal( 1 ./ Ssqrt )
+  Ssqrt = Diagonal( Ssqrt )
+  xSqrt = U*Ssqrt*transpose(U)
+  xSqrtInv = U*SsqrtInv*transpose(U)
+  tv = xSqrtInv * v * xSqrtInv
+  ty = xSqrtInv * y * xSqrtInv
+  svd2 = svd( 0.5*( ty + transpose(ty) ) )
+  Se = Diagonal( log.(svd2.S) )
+  Ue = svd2.U
+  ty2 = Ue*Se*transpose(Ue)
+  eig1 = eigen(  0.5 * (ty2 + transpose(ty2) )  )
+  Sf = Diagonal( exp.(eig1.values) )
+  Uf = eig1.vectors
+  vto = xSqrt*Uf*Sf*transpose(Uf)*(0.5*(tv+transpose(tv)))*Uf*Sf*transpose(Uf)*xSqrt
+  return vto
+end
+
+injectivity_radius(P::SymmetricPositiveDefinite, args...) = Infπ
+
+zero_tangent_vector(P::SymmetricPositiveDefinite, x) = zero(x)
+zero_tangent_vector(P::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric},x) = zero(x)
+zero_tangent_vector(P::MetricManifold{SymmetricPositiveDefinite{N},LogEuclidean},x) = zero(x)
+
+zero_tangent_vector!(P::SymmetricPositiveDefinite, v, x) = (v .= zero(x))
+zero_tangent_vector!(P::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric},v, x) = (v .= zero(x))
+zero_tangent_vector!(P::MetricManifold{SymmetricPositiveDefinite{N},LogEuclidean},v, x) = (v .= zero(x))
+
+"""
+    is_manifold_point(S,x; kwargs...)
+
+checks, whether `x` is a valid point on the [`SymmetricPositiveDefinite{N}`](@ref) `P`, 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...`.
+"""
+function is_manifold_point(P::SymmetricPositiveDefinite{N},x; kwargs...) where {N}
+    if size(x) != (N,N)
+        throw(DomainError(size(x),"The point $(x) does not lie on $P, since its size is not $( (N,N) )."))
+    end
+    if !isapprox(norm(x-transpose(x)), 0.; kwargs...)
+        throw(DomainError(norm(x), "The point $x does not lie on the sphere $P since its not a symmetric matrix:"))
+    end
+    if ! all( eigvals(x) > 0 )
+        throw(DomainError(norm(x), "The point $x does not lie on the sphere $P since its not a positive definite matrix."))
+    end
+    return true
+end
+
+"""
+    is_tangent_vector(S,x,v; kwargs... )
+
+checks whether `v` is a tangent vector to `x` on the [`Sphere`](@ref) `S`, i.e.
+atfer [`is_manifold_point`](@ref)`(S,x)`, `v` has to be of same dimension as `x`
+and orthogonal to `x`.
+The tolerance for the last test can be set using the ´kwargs...`.
+"""
+function is_tangent_vector(P::SymmetricPositiveDefinite{N},x,v; kwargs...) where N
+    is_manifold_point(P,x)
+    if size(v) != (N,N)
+        throw(DomainError(size(v),
+            "The vector $(v) is not a tangent to a point on $S since its size does not match $( (N,N) )."))
+    end
+    if !isapprox(norm(v-transpose(v)), 0.; kwargs...)
+        throw(DomainError(size(v),
+            "The vector $(v) is not a tangent to a point on $S (represented as an element of the Lie algebrasince its not symmetric."))
+    end
+    return true
+end