Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix check_point for SPD #593

Merged
merged 2 commits into from
Apr 28, 2023
Merged

Fix check_point for SPD #593

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4144489
Fix check_point for SPD
mateuszbaran Apr 27, 2023
d6ebea0
increase accuracy, remove unnecessary checks
mateuszbaran Apr 28, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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"
Expand Down
5 changes: 1 addition & 4 deletions src/groups/rotation_action.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
2 changes: 1 addition & 1 deletion src/manifolds/SymmetricPositiveDefinite.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.",
Expand Down
2 changes: 2 additions & 0 deletions src/manifolds/SymmetricPositiveDefiniteGeneralizedBuresWasserstein.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Comment on lines +97 to +98
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How is this increasing accuracy? How do you know the error in q is "symmetric", and reduced by symmetrization? Also, I wonder if this is optimal with respect to memory consumption? Perhaps something like

temp = m * Y * p * Y * m
temp .= temp .+ p .+ X
copyto!(q, temp)
# or
# copyto!(q, Symmetric(temp))
# or 
# copyto!(q, (temp .+ temp') ./ 2) # which first allocates the result and then copies it to `q`
# or, starting from Julia v1.10
# hermitianpart!(q, temp)

?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, it's more accurate than the previous implementation in the sense that the result is numerically closer to the manifold. I should've made a better comment there.

Maybe taking either the upper or lower part instead of average would be a better idea, I didn't check that. I didn't care much about memory consumption here because previous lines are already quite expensive.

return q
end

Expand Down
9 changes: 9 additions & 0 deletions test/manifolds/symmetric_positive_definite.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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