(Real) Symplectic Stiefel

The SymplecticStiefel manifold, denoted $\operatorname{SpSt}(2n, 2k)$, represents canonical symplectic bases of $2k$ dimensonal symplectic subspaces of $ℝ^{2nΓ—2n}$. This means that the columns of each element $p \in \operatorname{SpSt}(2n, 2k) \subset ℝ^{2nΓ—2k}$ constitute a canonical symplectic basis of $\operatorname{span}(p)$. The canonical symplectic form is a non-degenerate, bilinear, and skew symmetric map $\omega_{2k}\colon 𝔽^{2k}×𝔽^{2k} β†’ 𝔽$, given by $\omega_{2k}(x, y) = x^T Q_{2k} y$ for elements $x, y \in 𝔽^{2k}$, with

\[ Q_{2k} = \begin{bmatrix} 0_k & I_k \\ -I_k & 0_k \end{bmatrix}.\]

Specifically given an element $p \in \operatorname{SpSt}(2n, 2k)$ we require that

\[ \omega_{2n} (p x, p y) = x^T(p^TQ_{2n}p)y = x^TQ_{2k}y = \omega_{2k}(x, y) \;\forall\; x, y \in 𝔽^{2k},\]

leading to the requirement on $p$ that $p^TQ_{2n}p = Q_{2k}$. In the case that $k = n$, this manifold reduces to the SymplecticMatrices manifold, which is also known as the symplectic group.

Manifolds.SymplecticStiefel β€” Type
SymplecticStiefel{T,𝔽} <: AbstractEmbeddedManifold{𝔽, DefaultIsometricEmbeddingType}

The symplectic Stiefel manifold consists of all $2nΓ—2k, n β‰₯ k$ matrices satisfying the requirement

\[\mathrm{SpSt}(2n, 2k, ℝ) := \bigl\{ p ∈ ℝ^{2nΓ—2n} \ \big| \ p^{\mathrm{T}}J_{2n}p = J_{2k} \bigr\},\]

where $J_{2n}$ denotes the SymplecticElement

\[J_{2n} = \begin{bmatrix} 0_n & I_n \\ -I_n & 0_n \end{bmatrix}.\]

The symplectic Stiefel tangent space at $p$ can be parametrized as [BZ21]

\[\begin{align*} T_p\mathrm{SpSt}(2n, 2k) &= \{X ∈ ℝ^{2nΓ—2k} ∣ p^{T}J_{2n}X + X^{T}J_{2n}p = 0 \}, \\ &= \{X = pΞ© + p^sB \mid Ξ© ∈ ℝ^{2kΓ—2k}, Ξ©^+ = -Ξ©, \\ &\quad\qquad p^s ∈ \mathrm{SpSt}(2n, 2(n- k)), B ∈ ℝ^{2(n-k)Γ—2k}, \}, \end{align*}\]

where $Ω ∈ \mathfrak{sp}(2n,F)$ is Hamiltonian and $p^s$ means the symplectic complement of $p$ s.t. $p^{+}p^{s} = 0$. Here $p^+$ denotes the symplectic_inverse.

You can also use StiefelPoint and StiefelTVector with this manifold, they are equivalent to using arrays.


SymplecticStiefel(2n::Int, 2k::Int, field::AbstractNumbers=ℝ; parameter::Symbol=:type)

Generate the (real-valued) symplectic Stiefel manifold of $2nΓ—2k$ matrices which span a $2k$ dimensional symplectic subspace of $ℝ^{2nΓ—2n}$. The constructor for the SymplecticStiefel manifold accepts the even column dimension $2n$ and an even number of columns $2k$ for the real symplectic Stiefel manifold with elements $p ∈ ℝ^{2nΓ—2k}$.

Base.exp β€” Method
exp(::SymplecticStiefel, p, X)
exp!(M::SymplecticStiefel, q, p, X)

Compute the exponential mapping

\[ \exp\colon T\mathrm{SpSt}(2n, 2k) β†’ \mathrm{SpSt}(2n, 2k)\]

at a point $p ∈ \mathrm{SpSt}(2n, 2k)$ in the direction of $X ∈ T_p\mathrm{SpSt}(2n, 2k)$.

The tangent vector $X$ can be written in the form $X = \bar{\Omega}p$ [BZ21], with

\[ \bar{\Omega} = X (p^{\mathrm{T}}p)^{-1}p^{\mathrm{T}} + J_{2n}p(p^{\mathrm{T}}p)^{-1}X^{\mathrm{T}}(I_{2n} - J_{2n}^{\mathrm{T}}p(p^{\mathrm{T}}p)^{-1}p^{\mathrm{T}}J_{2n})J_{2n} ∈ ℝ^{2nΓ—2n},\]

where $J_{2n} = \begin{bmatrix} 0_n & I_n \\ -I_n & 0_n \end{bmatrix}$ denotes the SymplecticElement.

Using this expression for $X$, the exponential mapping can be computed as

\[ \exp_p(X) = \operatorname{Exp}([\bar{\Omega} - \bar{\Omega}^{\mathrm{T}}]) \operatorname{Exp}(\bar{\Omega}^{\mathrm{T}})p,\]

where $\operatorname{Exp}(β‹…)$ denotes the matrix exponential.

Computing the above mapping directly however, requires taking matrix exponentials of two $2nΓ—2n$ matrices, which is computationally expensive when $n$ increases. Therefore we instead follow [BZ21] who express the above exponential mapping in a way which only requires taking matrix exponentials of an $8kΓ—8k$ matrix and a $4kΓ—4k$ matrix.

To this end, first define

\[\bar{A} = J_{2k}p^{\mathrm{T}}X(p^{\mathrm{T}}p)^{-1}J_{2k} + (p^{\mathrm{T}}p)^{-1}X^{\mathrm{T}}(p - J_{2n}^{\mathrm{T}}p(p^{\mathrm{T}}p)^{-1}J_{2k}) ∈ ℝ^{2kΓ—2k},\]


\[\bar{H} = (I_{2n} - pp^+)J_{2n}X(p^{\mathrm{T}}p)^{-1}J_{2k} ∈ ℝ^{2nΓ—2k}.\]

We then let $\bar{\Delta} = p\bar{A} + \bar{H}$, and define the matrices

\[ Ξ³ = \left[\left(I_{2n} - \frac{1}{2}pp^+\right)\bar{\Delta} \quad -p \right] ∈ ℝ^{2nΓ—4k},\]


\[ Ξ» = \left[J_{2n}^{\mathrm{T}}pJ_{2k} \quad \left(\bar{\Delta}^+\left(I_{2n} - \frac{1}{2}pp^+\right)\right)^{\mathrm{T}}\right] ∈ ℝ^{2nΓ—4k}.\]

With the above defined matrices it holds that $\bar{\Omega} = λγ^{\mathrm{T}}$. As a last preliminary step, concatenate $Ξ³$ and $Ξ»$ to define the matrices $Ξ“ = [Ξ» \quad -Ξ³] ∈ ℝ^{2nΓ—8k}$ and $Ξ› = [Ξ³ \quad Ξ»] ∈ ℝ^{2nΓ—8k}$.

With these matrix constructions done, we can compute the exponential mapping as

\[ \exp_p(X) = Ξ“ \operatorname{Exp}(ΛΓ^{\mathrm{T}}) \begin{bmatrix} 0_{4k} \\ I_{4k} \end{bmatrix} \operatorname{Exp}(λγ^{\mathrm{T}}) \begin{bmatrix} 0_{2k} \\ I_{2k} \end{bmatrix}.\]

which only requires computing the matrix exponentials of $ΛΓ^{\mathrm{T}} ∈ ℝ^{8kΓ—8k}$ and $λγ^{\mathrm{T}} ∈ ℝ^{4kΓ—4k}$.

Base.inv β€” Method
inv(::SymplecticStiefel, A)
inv!(::SymplecticStiefel, q, p)

Compute the symplectic inverse $A^+$ of matrix $A ∈ ℝ^{2nΓ—2k}$. Given a matrix

\[A ∈ ℝ^{2nΓ—2k},\quad A = \begin{bmatrix} A_{1, 1} & A_{1, 2} \\ A_{2, 1} & A_{2, 2} \end{bmatrix}, \quad A_{i, j} ∈ ℝ^{2nΓ—2k}\]

the symplectic inverse is defined as:

\[A^{+} := J_{2k}^{\mathrm{T}} A^{\mathrm{T}} J_{2n},\]

where $J_{2n} = \begin{bmatrix} 0_n & I_n \\ -I_n & 0_n \end{bmatrix}$ denotes the SymplecticElement.

The symplectic inverse of a matrix A can be expressed explicitly as:

\[A^{+} = \begin{bmatrix} A_{2, 2}^{\mathrm{T}} & -A_{1, 2}^{\mathrm{T}} \\[1.2mm] -A_{2, 1}^{\mathrm{T}} & A_{1, 1}^{\mathrm{T}} \end{bmatrix}.\]

Base.rand β€” Method
rand(M::SymplecticStiefel; vector_at=nothing, Οƒ = 1.0)

Generate a random point $p ∈ \mathrm{SpSt}(2n, 2k)$ or a random tangent vector $X ∈ T_p\mathrm{SpSt}(2n, 2k)$ if vector_at is set to a point $p ∈ \mathrm{Sp}(2n)$.

A random point on $\mathrm{SpSt}(2n, 2k)$ is found by first generating a random point on the symplectic manifold $\mathrm{Sp}(2n)$, and then projecting onto the Symplectic Stiefel manifold using the canonical_project $Ο€_{\mathrm{SpSt}(2n, 2k)}$. That is, $p = Ο€_{\mathrm{SpSt}(2n, 2k)}(p_{\mathrm{Sp}})$.

To generate a random tangent vector in $T_p\mathrm{SpSt}(2n, 2k)$ this code exploits the second tangent vector space parametrization of SymplecticStiefel, that any $X ∈ T_p\mathrm{SpSt}(2n, 2k)$ can be written as $X = pΞ©_X + p^sB_X$. To generate random tangent vectors at $p$ then, this function sets $B_X = 0$ and generates a random Hamiltonian matrix $Ξ©_X ∈ \mathfrak{sp}(2n,F)$ with Frobenius norm of Οƒ before returning $X = pΞ©_X$.

ManifoldDiff.riemannian_gradient β€” Method
X = riemannian_gradient(::SymplecticStiefel, f, p, Y; embedding_metric::EuclideanMetric=EuclideanMetric())
riemannian_gradient!(::SymplecticStiefel, f, X, p, Y; embedding_metric::EuclideanMetric=EuclideanMetric())

Compute the riemannian gradient X of f on SymplecticStiefel at a point p, provided that the gradient of the function $\tilde f$, which is f continued into the embedding is given by Y. The metric in the embedding is the Euclidean metric.

The manifold gradient X is computed from Y as

\[ X = Yp^{\mathrm{T}}p + J_{2n}pY^{\mathrm{T}}J_{2n}p,\]

where $J_{2n} = \begin{bmatrix} 0_n & I_n \\ -I_n & 0_n \end{bmatrix}$ denotes the SymplecticElement.

Manifolds.canonical_project β€” Method
canonical_project(::SymplecticStiefel, p_Sp)
canonical_project!(::SymplecticStiefel, p, p_Sp)

Define the canonical projection from $\mathrm{Sp}(2n, 2n)$ onto $\mathrm{SpSt}(2n, 2k)$, by projecting onto the first $k$ columns and the $n + 1$'th onto the $n + k$'th columns [BZ21].

It is assumed that the point $p$ is on $\mathrm{Sp}(2n, 2n)$.

Manifolds.symplectic_inverse_times β€” Method
symplectic_inverse_times(::SymplecticStiefel, p, q)
symplectic_inverse_times!(::SymplecticStiefel, A, p, q)

Directly compute the symplectic inverse of $p ∈ \mathrm{SpSt}(2n, 2k)$, multiplied with $q ∈ \mathrm{SpSt}(2n, 2k)$. That is, this function efficiently computes $p^+q = (J_{2k}p^{\mathrm{T}}J_{2n})q ∈ ℝ^{2kΓ—2k}$, where $J_{2n}, J_{2k}$ are the SymplecticElement of sizes $2nΓ—2n$ and $2kΓ—2k$ respectively.

This function performs this common operation without allocating more than a $2kΓ—2k$ matrix to store the result in, or in the case of the in-place function, without allocating memory at all.

ManifoldsBase.check_vector β€” Method
check_vector(M::SymplecticMatrices, p, X; kwargs...)

Checks whether X is a valid tangent vector at p on the SymplecticStiefel, $\mathrm{SpSt}(2n, 2k)$ manifold.

The check consists of verifying that $H = p^{+}X ∈ 𝔀_{2k}$, where $𝔀$ is the Lie Algebra of the symplectic group $\mathrm{Sp}(2k)$, that is the set of [HamiltonianMatrices])(@ref), where $(β‹…)^+$ denotes the symplectic_inverse.

ManifoldsBase.inner β€” Method
inner(M::SymplecticStiefel, p, X. Y)

Compute the Riemannian inner product $g^{\mathrm{SpSt}}$ at $p ∈ \mathrm{SpSt}$ of tangent vectors $Y, X ∈ T_p\mathrm{SpSt}$. Given by Proposition 3.10 in [BZ21].

\[g^{\mathrm{SpSt}}_p(X, Y) = \operatorname{tr}\Bigl( X^{\mathrm{T}}\bigl( I_{2n} - \frac{1}{2}J_{2n}^{\mathrm{T}} p(p^{\mathrm{T}}p)^{-1}p^{\mathrm{T}}J_{2n} \bigr) Y (p^{\mathrm{T}}p)^{-1}\Bigr).\]

ManifoldsBase.inverse_retract β€” Method
inverse_retract(::SymplecticStiefel, p, q, ::CayleyInverseRetraction)
inverse_retract!(::SymplecticStiefel, q, p, X, ::CayleyInverseRetraction)

Compute the Cayley Inverse Retraction $X = \mathcal{L}_p^{\mathrm{SpSt}}(q)$ such that the Cayley Retraction from $p$ along $X$ lands at $q$, i.e. $\mathcal{R}_p(X) = q$ [BZ21].

For $p, q ∈ \mathrm{SpSt}(2n, 2k, ℝ)$ we can define the inverse cayley retraction as long as the following matrices exist.

\[ U = (I + p^+ q)^{-1} ∈ ℝ^{2kΓ—2k}, \quad V = (I + q^+ p)^{-1} ∈ ℝ^{2kΓ—2k},\]

where $(β‹…)^+$ denotes the symplectic_inverse.

THen the inverse retraction reads

\[\mathcal{L}_p^{\mathrm{Sp}}(q) = 2p\bigl(V - U\bigr) + 2\bigl((p + q)U - p\bigr) ∈ T_p\mathrm{Sp}(2n).\]

ManifoldsBase.manifold_dimension β€” Method

Returns the dimension of the symplectic Stiefel manifold embedded in $ℝ^{2nΓ—2k}$, i.e. [BZ21]

\[ \operatorname{dim}(\mathrm{SpSt}(2n, 2k)) = (4n - 2k + 1)k.\]

ManifoldsBase.project β€” Method
project(::SymplecticStiefel, p, A)
project!(::SymplecticStiefel, Y, p, A)

Given a point $p ∈ \mathrm{SpSt}(2n, 2k)$, project an element $A ∈ ℝ^{2nΓ—2k}$ onto the tangent space $T_p\mathrm{SpSt}(2n, 2k)$ relative to the euclidean metric of the embedding $ℝ^{2nΓ—2k}$.

That is, we find the element $X ∈ T_p\mathrm{SpSt}(2n, 2k)$ which solves the constrained optimization problem

\[ \displaystyle\operatorname{min}_{X ∈ ℝ^{2nΓ—2k}} \frac{1}{2}||X - A||^2, \quad \text{s.t.}\; h(X) := X^{\mathrm{T}} J p + p^{\mathrm{T}} J X = 0,\]

where $h : ℝ^{2nΓ—2k} β†’ \operatorname{skew}(2k)$ defines the restriction of $X$ onto the tangent space $T_p\mathrm{SpSt}(2n, 2k)$.

ManifoldsBase.retract β€” Method
retract(::SymplecticStiefel, p, X, ::CayleyRetraction)
retract!(::SymplecticStiefel, q, p, X, ::CayleyRetraction)

Compute the Cayley retraction on the Symplectic Stiefel manifold, from p along X (computed inplace of q).

Given a point $p ∈ \mathrm{SpSt}(2n, 2k)$, every tangent vector $X ∈ T_p\mathrm{SpSt}(2n, 2k)$ is of the form $X = \tilde{\Omega}p$, with

\[ \tilde{\Omega} = \left(I_{2n} - \frac{1}{2}pp^+\right)Xp^+ - pX^+\left(I_{2n} - \frac{1}{2}pp^+\right) ∈ ℝ^{2nΓ—2n},\]

as shown in Proposition 3.5 of [BZ21]. Using this representation of $X$, the Cayley retraction on $\mathrm{SpSt}(2n, 2k)$ is defined pointwise as

\[ \mathcal{R}_p(X) = \operatorname{cay}\left(\frac{1}{2}\tilde{\Omega}\right)p.\]

The operator $\operatorname{cay}(A) = (I - A)^{-1}(I + A)$ is the Cayley transform.

However, the computation of an $2nΓ—2n$ matrix inverse in the expression above can be reduced down to inverting a $2kΓ—2k$ matrix due to Proposition 5.2 of [BZ21].

Let $A = p^+X$ and $H = X - pA$. Then an equivalent expression for the Cayley retraction defined pointwise above is

\[ \mathcal{R}_p(X) = -p + (H + 2p)(H^+H/4 - A/2 + I_{2k})^{-1}.\]

This expression is computed inplace of q.



T.Β Bendokat and R.Β Zimmermann. The real symplectic Stiefel and Grassmann manifolds: metrics, geodesics and applications, arXivΒ Preprint,Β 2108.12447 (2021), arXiv:2108.12447.