Skip to content

Commit

Permalink
Fix embedding functions.
Browse files Browse the repository at this point in the history
  • Loading branch information
kellertuer committed Jan 21, 2024
1 parent b94a8e4 commit e10d353
Show file tree
Hide file tree
Showing 3 changed files with 45 additions and 27 deletions.
27 changes: 0 additions & 27 deletions src/manifolds/SymplecticGrassmann.jl
Original file line number Diff line number Diff line change
Expand Up @@ -73,33 +73,6 @@ end
# Define Stiefel as the array fallback
ManifoldsBase.@default_manifold_fallbacks SymplecticGrassmann StiefelPoint StiefelTVector value value

@doc raw"""
inner(::SymplecticGrassmann, p, X, Y)
Compute the Riemannian inner product ``g^{\mathrm{SpGr}}_p(X,Y)``, where ``p``
is a point on the [`SymplecticStiefel`](@ref) manifold and ``X,Y \in \mathrm{Hor}_p^π\operatorname{SpSt}(2n,2k)``
are horizontal tangent vectors. The formula reads according to Proposition Lemma 4.8 [BendokatZimmermann:2021](@cite).
```math
g^{\mathrm{SpGr}}_p(X,Y) = \operatorname{tr}\bigl(
(p^{\mathrm{T}}p)^{-1}X^{\mathrm{T}}(I_{2n} - pp^+)Y
\bigr),
```
where ``I_{2n}`` denotes the identity matrix and ``(⋅)^+`` the [`symplectic_inverse`](@ref).
"""
function inner(M::SymplecticGrassmann, p, X, Y)
n, k = get_parameter(M.size)
J = SymplecticElement(p, X, Y) # in BZ21 also J
# Procompute lu(p'p) since we solve a^{-1}* 3 times
a = lu(p' * p) # note that p'p is symmetric, thus so is its inverse c=a^{-1}
# we split the original trace into two one with I -> (X'Yc)
# 1) we permute X' and Y c to c^{\mathrm{T}}Y^{\mathrm{T}}X = a\(Y'X) (avoids a large interims matrix)
# 2) the second we compute as c (X'p)(p^+Y) since both brackets are the smaller matrices
return tr(a \ (Y' * X)) - tr(
a \ ((X' * p) * symplectic_inverse_times(SymplecticStiefel(2 * n, 2 * k), p, Y)),
)
end

@doc raw"""
manifold_dimension(::SymplecticGrassmann)
Expand Down
18 changes: 18 additions & 0 deletions src/manifolds/SymplecticGrassmannProjector.jl
Original file line number Diff line number Diff line change
Expand Up @@ -70,3 +70,21 @@ function check_vector(
)
end
end

embed!(::SymplecticGrassmann, q, p::ProjectorPoint) = copyto!(q, p.value)
function embed!(::SymplecticGrassmann, Y, p::ProjectorPoint, X::ProjectorTVector)
return copyto!(Y, X.value)

Check warning on line 76 in src/manifolds/SymplecticGrassmannProjector.jl

View check run for this annotation

Codecov / codecov/patch

src/manifolds/SymplecticGrassmannProjector.jl#L74-L76

Added lines #L74 - L76 were not covered by tests
end
embed(::SymplecticGrassmann, p::ProjectorPoint) = p.value
embed(::SymplecticGrassmann, p::ProjectorPoint, X::ProjectorTVector) = X.value

Check warning on line 79 in src/manifolds/SymplecticGrassmannProjector.jl

View check run for this annotation

Codecov / codecov/patch

src/manifolds/SymplecticGrassmannProjector.jl#L78-L79

Added lines #L78 - L79 were not covered by tests

function get_embedding(

Check warning on line 81 in src/manifolds/SymplecticGrassmannProjector.jl

View check run for this annotation

Codecov / codecov/patch

src/manifolds/SymplecticGrassmannProjector.jl#L81

Added line #L81 was not covered by tests
::SymplecticGrassmann{TypeParameter{Tuple{n,k}},𝔽},
p::ProjectorPoint,
) where {n,k,𝔽}
return Euclidean(2n, 2n; field=𝔽)

Check warning on line 85 in src/manifolds/SymplecticGrassmannProjector.jl

View check run for this annotation

Codecov / codecov/patch

src/manifolds/SymplecticGrassmannProjector.jl#L85

Added line #L85 was not covered by tests
end
function get_embedding(M::SymplecticGrassmann{Tuple{Int,Int},𝔽}, ::ProjectorPoint) where {𝔽}
n, _ = get_parameter(M.size)
return Euclidean(2n, 2n; field=𝔽, parameter=:field)

Check warning on line 89 in src/manifolds/SymplecticGrassmannProjector.jl

View check run for this annotation

Codecov / codecov/patch

src/manifolds/SymplecticGrassmannProjector.jl#L87-L89

Added lines #L87 - L89 were not covered by tests
end
27 changes: 27 additions & 0 deletions src/manifolds/SymplecticGrassmannStiefel.jl
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,33 @@ function get_embedding(M::SymplecticGrassmann{Tuple{Int,Int},𝔽}) where {𝔽}
return SymplecticStiefel(2n, 2k, 𝔽; parameter=:field)
end

@doc raw"""
inner(::SymplecticGrassmann, p, X, Y)
Compute the Riemannian inner product ``g^{\mathrm{SpGr}}_p(X,Y)``, where ``p``
is a point on the [`SymplecticStiefel`](@ref) manifold and ``X,Y \in \mathrm{Hor}_p^π\operatorname{SpSt}(2n,2k)``
are horizontal tangent vectors. The formula reads according to Proposition Lemma 4.8 [BendokatZimmermann:2021](@cite).
```math
g^{\mathrm{SpGr}}_p(X,Y) = \operatorname{tr}\bigl(
(p^{\mathrm{T}}p)^{-1}X^{\mathrm{T}}(I_{2n} - pp^+)Y
\bigr),
```
where ``I_{2n}`` denotes the identity matrix and ``(⋅)^+`` the [`symplectic_inverse`](@ref).
"""
function inner(M::SymplecticGrassmann, p, X, Y)
n, k = get_parameter(M.size)
J = SymplecticElement(p, X, Y) # in BZ21 also J

Check warning on line 81 in src/manifolds/SymplecticGrassmannStiefel.jl

View check run for this annotation

Codecov / codecov/patch

src/manifolds/SymplecticGrassmannStiefel.jl#L79-L81

Added lines #L79 - L81 were not covered by tests
# Procompute lu(p'p) since we solve a^{-1}* 3 times
a = lu(p' * p) # note that p'p is symmetric, thus so is its inverse c=a^{-1}

Check warning on line 83 in src/manifolds/SymplecticGrassmannStiefel.jl

View check run for this annotation

Codecov / codecov/patch

src/manifolds/SymplecticGrassmannStiefel.jl#L83

Added line #L83 was not covered by tests
# we split the original trace into two one with I -> (X'Yc)
# 1) we permute X' and Y c to c^{\mathrm{T}}Y^{\mathrm{T}}X = a\(Y'X) (avoids a large interims matrix)
# 2) the second we compute as c (X'p)(p^+Y) since both brackets are the smaller matrices
return tr(a \ (Y' * X)) - tr(

Check warning on line 87 in src/manifolds/SymplecticGrassmannStiefel.jl

View check run for this annotation

Codecov / codecov/patch

src/manifolds/SymplecticGrassmannStiefel.jl#L87

Added line #L87 was not covered by tests
a \ ((X' * p) * symplectic_inverse_times(SymplecticStiefel(2 * n, 2 * k), p, Y)),
)
end

@doc raw"""
inverse_retract(::SymplecticGrassmann, p, q, ::CayleyInverseRetraction)
inverse_retract!(::SymplecticGrassmann, q, p, X, ::CayleyInverseRetraction)
Expand Down

0 comments on commit e10d353

Please sign in to comment.