From af8136937161b66cc96b4a5e975c9c838017b1ee Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Tue, 17 Nov 2020 16:14:44 +0100 Subject: [PATCH] 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 * remove an unnecessary case for the tangent bundle. * bump version. Co-authored-by: Mateusz Baran --- Project.toml | 2 +- .../SymmetricPositiveDefiniteLinearAffine.jl | 14 +++++++------- src/manifolds/VectorBundle.jl | 6 ++++-- 3 files changed, 12 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index a224b5425c..a2ac3c5b48 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Manifolds" uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann ", "Antoine Levitt "] -version = "0.4.6" +version = "0.4.7" [deps] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" diff --git a/src/manifolds/SymmetricPositiveDefiniteLinearAffine.jl b/src/manifolds/SymmetricPositiveDefiniteLinearAffine.jl index 55384c0c7f..2e8d10bb6d 100644 --- a/src/manifolds/SymmetricPositiveDefiniteLinearAffine.jl +++ b/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 diff --git a/src/manifolds/VectorBundle.jl b/src/manifolds/VectorBundle.jl index 2d0cf3564b..45fb056fa2 100644 --- a/src/manifolds/VectorBundle.jl +++ b/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)