Stabilize SPD distance and a small fix to TangentBundle (#256)

* stabilize distance and introduce a proper inner for tangent bundle.

* Update src/manifolds/VectorBundle.jl

Co-authored-by: Mateusz Baran <[email protected]>

* remove an unnecessary case for the tangent bundle.

* bump version.

Co-authored-by: Mateusz Baran <[email protected]>

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.4.6"
+version = "0.4.7"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

src/manifolds/SymmetricPositiveDefiniteLinearAffine.jl

@@ -22,10 +22,10 @@ d_{\mathcal P(n)}(p,q)
 where $\operatorname{Log}$ denotes the matrix logarithm and
 $\lVert\cdot\rVert_{\mathrm{F}}$ denotes the matrix Frobenius norm.
 """
-function distance(M::SymmetricPositiveDefinite{N}, p, q) where {N}
+function distance(::SymmetricPositiveDefinite{N}, p, q) where {N}
     # avoid numerical instabilities in cholesky
     norm(p - q) < eps(eltype(p + q)) && return zero(eltype(p + q))
-    cq = cholesky(q)
+    cq = cholesky(Symmetric(q)) # to avoid numerical inaccuracies
     s = eigvals(Symmetric(cq.L \ p / cq.U))
     return any(s .<= eps()) ? zero(eltype(p)) : sqrt(sum(abs.(log.(s)) .^ 2))
 end
@@ -46,7 +46,7 @@ where $\operatorname{Exp}$ denotes to the matrix exponential.
 """
 exp(::SymmetricPositiveDefinite, ::Any...)
 
-function exp!(M::SymmetricPositiveDefinite{N}, q, p, X) where {N}
+function exp!(::SymmetricPositiveDefinite{N}, q, p, X) where {N}
     e = eigen(Symmetric(p))
     U = e.vectors
     S = e.values
@@ -73,7 +73,7 @@ Return a orthonormal basis `Ξ` as a vector of tangent vectors (of length
 with eigenvalues `κ` and where the direction `B.frame_direction` has curvature `0`.
 """
 function get_basis(
-    M::SymmetricPositiveDefinite{N},
+    ::SymmetricPositiveDefinite{N},
     p,
     B::DiagonalizingOrthonormalBasis,
 ) where {N}
@@ -95,7 +95,7 @@ function get_coordinates!(
     Y,
     p,
     X,
-    B::DefaultOrthonormalBasis,
+    ::DefaultOrthonormalBasis,
 ) where {N}
     dim = manifold_dimension(M)
     @assert size(Y) == (dim,)
@@ -115,7 +115,7 @@ function get_vector!(
     Y,
     p,
     X,
-    B::DefaultOrthonormalBasis,
+    ::DefaultOrthonormalBasis,
 ) where {N}
     dim = manifold_dimension(M)
     @assert size(X) == (div(N * (N + 1), 2),)
@@ -142,7 +142,7 @@ a [`MetricManifold`](@ref) with [`LinearAffineMetric`](@ref). The formula reads
 g_p(X,Y) = \operatorname{tr}(p^{-1} X p^{-1} Y),
 ````
 """
-function inner(M::SymmetricPositiveDefinite, p, X, Y)
+function inner(::SymmetricPositiveDefinite, p, X, Y)
     F = cholesky(Symmetric(p))
     return tr((F \ Symmetric(X)) * (F \ Symmetric(Y)))
 end

src/manifolds/VectorBundle.jl

@@ -76,7 +76,7 @@ CotangentBundleFibers(M::Manifold) = VectorBundleFibers(CotangentSpace, M)
         𝔽,
         TFiber<:VectorBundleFibers{<:VectorSpaceType,<:Manifold{𝔽}},
         TX,
-    } <: Manifold{𝔽}   
+    } <: Manifold{𝔽}
 
 A vector space at a point `p` on the manifold.
 This is modelled using [`VectorBundleFibers`](@ref) with only a vector-like part
@@ -633,8 +633,10 @@ function inner(B::VectorBundle, p, X, Y)
     px, Vx = submanifold_components(B.manifold, p)
     VXM, VXF = submanifold_components(B.manifold, X)
     VYM, VYF = submanifold_components(B.manifold, Y)
-    return inner(B.manifold, px, VXM, VYM) + inner(B.fiber, Vx, VXF, VYF)
+    return inner(B.manifold, px, VXM, VYM) +
+           inner(VectorSpaceAtPoint(B.fiber, px), Vx, VXF, VYF)
 end
+
 """
     inner(M::TangentSpaceAtPoint, p, X, Y)