Fix check_point for SPD (#593)

* Fix check_point for SPD * increase accuracy, remove unnecessary checks
JuliaManifolds · Apr 28, 2023 · 73cfc7b · 73cfc7b · mateuszbaran · Apr 28, 2023
1 parent 8ccd5d8
commit 73cfc7b
Show file tree

Hide file tree

Showing 5 changed files with 14 additions and 6 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Manifolds"
 uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
 authors = ["Seth Axen <seth.axen@gmail.com>", "Mateusz Baran <mateuszbaran89@gmail.com>", "Ronny Bergmann <manopt@ronnybergmann.net>", "Antoine Levitt <antoine.levitt@gmail.com>"]
-version = "0.8.57"
+version = "0.8.58"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"

diff --git a/src/groups/rotation_action.jl b/src/groups/rotation_action.jl
@@ -97,10 +97,7 @@ function inverse_apply_diff(A::RotationActionOnVector{N,F,LeftAction}, a, p, X)
 end
 inverse_apply_diff(A::RotationActionOnVector{N,F,RightAction}, a, p, X) where {N,F} = a * X
 
-function optimal_alignment(A::RotationActionOnVector{N,T,LeftAction}, p, q) where {N,T}
-    is_point(A.manifold, p, true)
-    is_point(A.manifold, q, true)
-
+function optimal_alignment(::RotationActionOnVector{N,T,LeftAction}, p, q) where {N,T}
     Xmul = p * transpose(q)
     F = svd(Xmul)
     L = size(Xmul)[2]

diff --git a/src/manifolds/SymmetricPositiveDefinite.jl b/src/manifolds/SymmetricPositiveDefinite.jl
@@ -138,7 +138,7 @@ function check_point(M::SymmetricPositiveDefinite{N}, p; kwargs...) where {N}
             "The point $(p) does not lie on $(M) since its not a symmetric matrix:",
         )
     end
-    if !all(eigvals(p) .> 0)
+    if !isposdef(p)
         return DomainError(
             eigvals(p),
             "The point $p does not lie on $(M) since its not a positive definite matrix.",

diff --git a/src/manifolds/SymmetricPositiveDefiniteGeneralizedBuresWasserstein.jl b/src/manifolds/SymmetricPositiveDefiniteGeneralizedBuresWasserstein.jl
@@ -94,6 +94,8 @@ function exp!(
     m = M.metric.M
     Y = lyapc(p, m, -X) #lyap solves qpM + Mpq - X =0
     q .= p .+ X .+ m * Y * p * Y * m
+    # symmetrizing for better accuracy
+    copyto!(q, (q .+ q') ./ 2)
     return q
 end
 

diff --git a/test/manifolds/symmetric_positive_definite.jl b/test/manifolds/symmetric_positive_definite.jl
@@ -267,4 +267,13 @@ include("../utils.jl")
         @test ismissing(pS.sqrt)
         @test ismissing(pS.sqrt_inv)
     end
+
+    @testset "test BigFloat" begin
+        M = SymmetricPositiveDefinite(2)
+        p1 = BigFloat[
+            1.6590891025248637458133771360735408961772918701171875 -2.708777790960681386422947980463504791259765625e-07
+            -2.708777790960681386422947980463504791259765625e-07 1.6590893171834280028775765458703972399234771728515625
+        ]
+        @test is_point(M, p1)
+    end
 end