Symplectic Grassmann

NEWS.md

@@ -5,6 +5,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [unreleased] 2023-12-x
+
+### Added
+
+* added the real symplectic Grassmann manifold `SymplecticGrassmann`
+
 ## [0.9.11] – 2023-12-27
 
 ### Fixed
@@ -37,27 +43,27 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 * Improved distribution of random vector generation for rotation matrices and complex circle.
 
-## [0.9.7] - 2023-11-14
+## [0.9.7] * 2023-11-14
 
 ### Fixed
 
 * Fixed `is_flat` for `CholeskySpace` and `SymmetricPositiveDefinite` with `LogCholeskyMetric` [https://github.com/JuliaManifolds/Manifolds.jl/issues/684](https://github.com/JuliaManifolds/Manifolds.jl/issues/684).
 
-## [0.9.6] - 2023-11-09
+## [0.9.6] * 2023-11-09
 
 ### Fixed
 
 * Fixed real coefficient basis for complex circle (an issue exposed by [https://github.com/JuliaManifolds/ManifoldsBase.jl/pull/173](https://github.com/JuliaManifolds/ManifoldsBase.jl/pull/173)).
 * Fixed `VeeOrthogonalBasis` test for non-real manifolds.
 
-## [0.9.5] - 2023-11-08
+## [0.9.5] * 2023-11-08
 
 ### Changed
 
 * `identity_element` now returns a complex matrix for unitary group.
 * `number_of_coordinates` is now exported.
 
-## [0.9.4] - 2023-11-06
+## [0.9.4] * 2023-11-06
 
 ### Added
 
@@ -72,27 +78,27 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 * Fixed issue with incorrect implementation of `apply_diff_group` in `GroupOperationAction` with left backward and right forward action [#669](https://github.com/JuliaManifolds/Manifolds.jl/issues/669).
 
-## [0.9.3] - 2023-10-28
+## [0.9.3] * 2023-10-28
 
 ### Added
 
 * Support for `BoundaryValueDiffEq` v5.
 
-## [0.9.2] - 2023-10-27
+## [0.9.2] * 2023-10-27
 
 ### Added
 
 * `rand(G; vector_at=Identity(G))` now works for translation, special orthogonal and special Euclidean groups `G` (issue [#665](https://github.com/JuliaManifolds/Manifolds.jl/issues/665)).
 * `get_embedding` now works for `GeneralUnitaryMultiplicationGroup`.
 * Github action that checks for NEWS.md changes.
 
-## [0.9.1] - 2023-10-25
+## [0.9.1] * 2023-10-25
 
 ### Added
 
 * a new retraction and its inverse for the fixed Rank Manifolds, the orthographic retraction.
 
-## [0.9.0] - 2023-10-24
+## [0.9.0] * 2023-10-24
 
 ### Added
 

docs/make.jl

@@ -113,6 +113,7 @@ makedocs(;
                 "Generalized Stiefel" => "manifolds/generalizedstiefel.md",
                 "Generalized Grassmann" => "manifolds/generalizedgrassmann.md",
                 "Grassmann" => "manifolds/grassmann.md",
+                "Hamiltonian" => "manifolds/hamiltonian.md",
                 "Hyperbolic space" => "manifolds/hyperbolic.md",
                 "Lorentzian manifold" => "manifolds/lorentz.md",
                 "Multinomial doubly stochastic matrices" => "manifolds/multinomialdoublystochastic.md",
@@ -134,6 +135,7 @@ makedocs(;
                 "SPD, fixed determinant" => "manifolds/spdfixeddeterminant.md",
                 "Symmetric positive semidefinite fixed rank" => "manifolds/symmetricpsdfixedrank.md",
                 "Symplectic" => "manifolds/symplectic.md",
+                "Symplectic Grassmann" => "manifolds/symplecticgrassmann.md",
                 "Symplectic Stiefel" => "manifolds/symplecticstiefel.md",
                 "Torus" => "manifolds/torus.md",
                 "Tucker" => "manifolds/tucker.md",

docs/src/manifolds/hamiltonian.md

@@ -0,0 +1,7 @@
+# Symmetric matrices
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["manifolds/Hamiltonian.jl"]
+Order = [:type, :function]
+```

docs/src/manifolds/symplecticgrassmann.md

@@ -0,0 +1,30 @@
+# (Real) Symplectic Grassmann
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["manifolds/SymplecticGrassmann.jl"]
+Order = [:type, :function]
+```
+
+## The (default) symplectic Stiefel representation
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["manifolds/SymplecticGrassmannStiefel.jl"]
+Order = [:type, :function]
+```
+
+## The symplectic projector representation
+
+```@autodocs
+Modules = [Manifolds]
+Pages = ["manifolds/SymplecticGrassmannProjector.jl"]
+Order = [:type, :function]
+```
+
+## Literature
+
+```@bibliography
+Pages = ["symplecticgrassmann.md"]
+Canonical=false
+```

docs/src/manifolds/symplecticstiefel.md

@@ -1,4 +1,4 @@
-# Symplectic Stiefel
+# (Real) Symplectic Stiefel
 
 The [`SymplecticStiefel`](@ref) manifold, denoted $\operatorname{SpSt}(2n, 2k)$,
 represents canonical symplectic bases of $2k$ dimensonal symplectic subspaces of $\mathbb{R}^{2n \times 2n}$.

src/Manifolds.jl

@@ -4,6 +4,7 @@
 module Manifolds
 
 import Base:
+    ^,
     angle,
     copyto!,
     convert,
@@ -444,8 +445,6 @@ include("manifolds/SymmetricPositiveDefiniteAffineInvariant.jl")
 include("manifolds/SymmetricPositiveDefiniteLogCholesky.jl")
 include("manifolds/SymmetricPositiveDefiniteLogEuclidean.jl")
 include("manifolds/SymmetricPositiveSemidefiniteFixedRank.jl")
-include("manifolds/Symplectic.jl")
-include("manifolds/SymplecticStiefel.jl")
 include("manifolds/Tucker.jl")
 #
 include("manifolds/ProbabilitySimplex.jl")
@@ -462,6 +461,13 @@ include("manifolds/KendallsShapeSpace.jl")
 # Introduce the quotient, Grassmann, only after Stiefel
 include("manifolds/Grassmann.jl")
 
+# Introduce Symplectic and so on manifolds only after Grassmann
+# Since that defines the StiefelPoint, StiefelTVector
+include("manifolds/Symplectic.jl")
+include("manifolds/Hamiltonian.jl") # Hamiltonian requires symplectic
+include("manifolds/SymplecticStiefel.jl")
+include("manifolds/SymplecticGrassmann.jl") # Requires SymplecticStiefel
+
 # Product or power based manifolds
 include("manifolds/Torus.jl")
 include("manifolds/Multinomial.jl")
@@ -659,6 +665,7 @@ export Euclidean,
     SPDFixedDeterminant,
     SymmetricPositiveSemidefiniteFixedRank,
     Symplectic,
+    SymplecticGrassmann,
     SymplecticStiefel,
     SymplecticMatrix,
     Torus,
@@ -682,7 +689,7 @@ export HyperboloidTVector,
     ProjectorTVector,
     StiefelTVector
 export AbstractNumbers, ℝ, ℂ, ℍ
-
+export Hamiltonian
 # decorator manifolds
 export AbstractDecoratorManifold
 export IsIsometricEmbeddedManifold, IsEmbeddedManifold, IsEmbeddedSubmanifold
@@ -775,6 +782,7 @@ export CachedBasis,
 export ComponentManifoldError, CompositeManifoldError
 # Functions on Manifolds
 export ×,
+    ^,
     action_side,
     allocate,
     allocate_result,
@@ -900,6 +908,8 @@ export ×,
     skewness,
     std,
     sym_rem,
+    symplectic_inverse,
+    symplectic_inverse!,
     symplectic_inverse_times,
     symplectic_inverse_times!,
     submanifold,

src/manifolds/Grassmann.jl

@@ -144,7 +144,7 @@ is_flat(M::Grassmann) = manifold_dimension(M) == 1
 @doc raw"""
     manifold_dimension(M::Grassmann)
 
-Return the dimension of the [`Grassmann(n,k,𝔽)`](@ref) manifold `M`, i.e.
+Return the dimension of the [`Grassmann`](@ref)`(n,k,𝔽)` manifold `M`, i.e.
 
 ````math
 \dim \operatorname{Gr}(n,k) = k(n-k) \dim_ℝ 𝔽,
@@ -195,7 +195,7 @@ function get_total_space(M::Grassmann{Tuple{Int,Int},𝔽}) where {𝔽}
 end
 
 #
-# Reprenter specific implementations in their corresponding subfiles
+# Representer specific implementations in their corresponding subfiles
 #
 include("GrassmannStiefel.jl")
 include("GrassmannProjector.jl")

src/manifolds/Hamiltonian.jl

@@ -0,0 +1,65 @@
+#
+# This file requires Symplectic to be defined, since it need the symplectic inverse A^+
+# This type is used in the Symplectic Grassmann
+@doc raw"""
+    Hamiltonian{T,S<:AbstractMatrix{<:T}} <: AbstractMatrix{T}
+
+A type to store an hamiltonien matrix, i.e. A square matrix matrix for which ``A^+ = -A`` where
+
+```math
+A^+ = J_{2n}A^{\mathrm{T}}J_{2n}, \qquad J_{2n} \begin{pmatrix} 0 & I_n\\-I_n & 0 \end{pmatrix},
+```
+
+and ``I_n`` denotes the ``n × n``
+"""
+struct Hamiltonian{T,S<:AbstractMatrix{<:T}} <: AbstractMatrix{T}
+    value::S
+    function Hamiltonian(A::S) where {T,S<:AbstractMatrix{<:T}}
+        n = div(size(A, 1), 2)
+        @assert size(A, 1) == 2 * n "The first dimension of A ($(size(A,1))) is not even"
+        @assert size(A, 2) == 2 * n "The matrix A is of size ($(size(A))), which is not square."
+        return new{T,S}(A)
+    end
+end
+# Avoid double wrapping / unwrap if that happened
+Hamiltonian(A::Hamiltonian) = Hamiltonian(A.value)
+# Conversion
+function Matrix(A::Hamiltonian)
+    return copy(A.value)
+end
+
+@doc raw"""
+    ^(A::Hamilonian, ::typeof(+))
+
+Compute the [`symplectic_inverse`](@ref) of a Hamiltonian (A)
+"""
+function ^(A::Hamiltonian, ::typeof(+))
+    return Hamiltonian(symplectic_inverse(A.value))
+end
+
+@doc raw"""
+    is_hamiltonian(A; kwargs...)
+
+Test whether a matrix `A` is hamiltonian.
+The test consists of verifying whether
+
+```math
+A^+ = -A
+```
+where ``A^+`` denotes the [`symplectic_inverse`](@ref) of `A`.
+
+The passed keyword arguments are passed on to the [`isapprox`](@ref)
+check within
+"""
+function is_hamiltonian(A::AbstractMatrix; kwargs...)
+    return isapprox(symplectic_inverse(A), -A; kwargs...)
+end
+function is_hamiltonian(A::Hamiltonian; kwargs...)
+    return isapprox((A^+).value, -A.value; kwargs...)
+end
+
+function show(io::IO, ::MIME"text/plain", A::Hamiltonian)
+    return print(io, "Hamiltonian($(A.value))")
+end
+
+size(A::Hamiltonian) = size(A.value)