Symplectic Stiefel

The SymplecticStiefel manifold, denoted $\operatorname{SpSt}(2n, 2k)$, represents canonical symplectic bases of $2k$ dimensonal symplectic subspaces of $\mathbb{R}^{2n \times 2n}$. This means that the columns of each element $p \in \operatorname{SpSt}(2n, 2k) \subset \mathbb{R}^{2n \times 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 \mathbb{F}^{2k} \times \mathbb{F}^{2k} \rightarrow \mathbb{F}$, given by $\omega_{2k}(x, y) = x^T Q_{2k} y$ for elements $x, y \in \mathbb{F}^{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 \mathbb{F}^{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 Symplectic manifold, which is also known as the symplectic group.

Manifolds.SymplecticStiefel โ€” Type
SymplecticStiefel{n, k, ๐”ฝ} <: AbstractEmbeddedManifold{๐”ฝ, DefaultIsometricEmbeddingType}

The symplectic Stiefel manifold consists of all $2n ร— 2k, \; n \geq k$ matrices satisfying the requirement

\[\operatorname{SpSt}(2n, 2k, โ„) = \bigl\{ p โˆˆ โ„^{2n ร— 2n} \, \big| \, p^TQ_{2n}p = Q_{2k} \bigr\},\]


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

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

\[ \begin{align*} T_p\operatorname{SpSt}(2n, 2k) = \{&X \in \mathbb{R}^{2n \times 2k} \;|\; p^{T}Q_{2n}X + X^{T}Q_{2n}p = 0 \}, \\ = \{&X = pฮฉ + p^sB \;|\; ฮฉ โˆˆ โ„^{2k ร— 2k}, ฮฉ^+ = -ฮฉ, \\ &\; p^s โˆˆ \operatorname{SpSt}(2n, 2(n- k)), B โˆˆ โ„^{2(n-k) ร— 2k}, \}, \end{align*}\]

where $ฮฉ \in \mathfrak{sp}(2n,F)$ is Hamiltonian and $p^s$ means the symplectic complement of $p$ s.t. $p^{+}p^{s} = 0$.


SymplecticStiefel(2n::Int, 2k::Int, field::AbstractNumbers=โ„)
    -> SymplecticStiefel{div(2n, 2), div(2k, 2), field}()

Generate the (real-valued) symplectic Stiefel manifold of $2n \times 2k$ matrices which span a $2k$ dimensional symplectic subspace of $โ„^{2n \times 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 \in โ„^{2n ร— 2k}$.

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

Compute the exponential mapping

\[ \operatorname{exp}\colon T\operatorname{SpSt}(2n, 2k) \rightarrow \operatorname{SpSt}(2n, 2k)\]

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

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

\[ \bar{\Omega} = X (p^Tp)^{-1}p^T + Q_{2n}p(p^Tp)^{-1}X^T(I_{2n} - Q_{2n}^Tp(p^Tp)^{-1}p^TQ_{2n})Q_{2n} \in โ„^{2n \times 2n},\]

where $Q_{2n}$ is the SymplecticMatrix. Using this expression for $X$, the exponential mapping can be computed as

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

where $\operatorname{Exp}(\cdot)$ denotes the matrix exponential.

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

To this end, first define

\[\bar{A} = Q_{2k}p^TX(p^Tp)^{-1}Q_{2k} + (p^Tp)^{-1}X^T(p - Q_{2n}^Tp(p^Tp)^{-1}Q_{2k}) \in โ„^{2k \times 2k},\]


\[\bar{H} = (I_{2n} - pp^+)Q_{2n}X(p^Tp)^{-1}Q_{2k} \in โ„^{2n \times 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] \in โ„^{2n \times 4k},\]


\[ ฮป = \left[Q_{2n}^TpQ_{2k} \quad \left(\bar{\Delta}^+\left(I_{2n} - \frac{1}{2}pp^+\right)\right)^T\right] \in โ„^{2n \times 4k}.\]

With the above defined matrices it holds that $\bar{\Omega} = ฮปฮณ^T$. As a last preliminary step, concatenate $ฮณ$ and $ฮป$ to define the matrices $ฮ“ = [ฮป \quad -ฮณ] \in โ„^{2n \times 8k}$ and $ฮ› = [ฮณ \quad ฮป] \in โ„^{2n \times 8k}$.

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

\[ \operatorname{exp}_p(X) = ฮ“ \operatorname{Exp}(ฮ›ฮ“^T) \begin{bmatrix} 0_{4k} \\ I_{4k} \end{bmatrix} \operatorname{Exp}(ฮปฮณ^T) \begin{bmatrix} 0_{2k} \\ I_{2k} \end{bmatrix}.\]

which only requires computing the matrix exponentials of $ฮ›ฮ“^T \in โ„^{8k \times 8k}$ and $ฮปฮณ^T \in โ„^{4k \times 4k}$.

Base.inv โ€” Method
inv(::SymplecticStiefel{n, k}, A)
inv!(::SymplecticStiefel{n, k}, 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},\; A_{i, j} \in โ„^{2n ร— 2k}\]

the symplectic inverse is defined as:

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


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

For any $p \in \operatorname{SpSt}(2n, 2k)$ we have that $p^{+}p = I_{2k}$.

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

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

Base.rand โ€” Method
rand(M::SymplecticStiefel; vector_at=nothing,
    hamiltonian_norm=(vector_at === nothing ? 1/2 : 1.0))

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

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

To generate a random tangent vector in $T_p\operatorname{SpSt}(2n, 2k)$ this code exploits the second tangent vector space parametrization of SymplecticStiefel, showing that any $X \in T_p\operatorname{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 \in \mathfrak{sp}(2n,F)$ with Frobenius norm of hamiltonian_norm before returning $X = pฮฉ_X$.

Manifolds.canonical_projection โ€” Method
canonical_projection(::SymplecticStiefel, p_Sp)
canonical_projection!(::SymplecticStiefel{n,k}, p, p_Sp)

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

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

Manifolds.gradient โ€” Method
gradient(::SymplecticStiefel, f, p, backend::RiemannianProjectionBackend)
gradient!(::SymplecticStiefel, f, X, p, backend::RiemannianProjectionBackend)

Compute the gradient of $f\colon \operatorname{SpSt}(2n, 2k) \rightarrow โ„$ at $p \in \operatorname{SpSt}(2n, 2k)$.

This function first computes the embedding gradient $โˆ‡f(p) \in โ„^{2n \times 2k}$ using the AbstractRiemannianDiffBackend in the RiemannianProjectionBackend. Then it transforms the embedding gradient to the unique tangent vector space element $\text{grad}f(p) \in T_p\operatorname{SpSt}(2n, 2k)$ which fulfills the variational equation

\[ g_p(\text{grad}f(p), X) = \text{D}f(p) = \langle โˆ‡f(p), X \rangle \quad\forall\; X \in T_p\operatorname{SpSt}(2n, 2k).\]

The manifold gradient $\text{grad}f(p)$ is computed from $โˆ‡f(p)$ as

\[ \text{grad}f(p) = โˆ‡f(p)p^Tp + Q_{2n}pโˆ‡f(p)^TQ_{2n}p,\]

where $Q_{2n}$ is the SymplecticMatrix.

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 \in \operatorname{SpSt}(2n, 2k)$, multiplied with $q \in \operatorname{SpSt}(2n, 2k)$. That is, this function efficiently computes $p^+q = (Q_{2k}p^TQ_{2n})q \in โ„^{2k \times 2k}$, where $Q_{2n}, Q_{2k}$ are the SymplecticMatrix of sizes $2n \times 2n$ and $2k \times 2k$ respectively.

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

ManifoldsBase.check_point โ€” Method
check_point(M::SymplecticStiefel, p; kwargs...)

Check whether p is a valid point on the SymplecticStiefel, $\operatorname{SpSt}(2n, 2k)$ manifold. That is, the point has the right AbstractNumbers type and $p^{+}p$ is (approximately) the identity, where for $A \in \mathbb{R}^{2n \times 2k}$, $A^{+} = Q_{2k}^TA^TQ_{2n}$ is the symplectic inverse, with

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

The tolerance can be set with kwargs... (e.g. atol = 1.0e-14).

ManifoldsBase.check_vector โ€” Method
check_vector(M::Symplectic, p, X; kwargs...)

Checks whether X is a valid tangent vector at p on the SymplecticStiefel, $\operatorname{SpSt}(2n, 2k)$ manifold. First recall the definition of the symplectic inverse for $A \in \mathbb{R}^{2n \times 2k}$, $A^{+} = Q_{2k}^TA^TQ_{2n}$ is the symplectic inverse, with

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

The we check that $H = p^{+}X \in ๐”ค_{2k}$, where $๐”ค$ is the Lie Algebra of the symplectic group $\operatorname{Sp}(2k)$, characterized as [Bendokat2021],

\[ ๐”ค_{2k} = \{H \in โ„^{2k \times 2k} \;|\; H^+ = -H \}.\]

The tolerance can be set with kwargs... (e.g. atol = 1.0e-14).

ManifoldsBase.inner โ€” Method
inner(M::SymplecticStiefel{n, k}, p, X. Y)

Compute the Riemannian inner product $g^{\operatorname{SpSt}}$ at $p \in \operatorname{SpSt}$ between tangent vectors $X, X \in T_p\operatorname{SpSt}$. Given by Proposition 3.10 in [Bendokat2021].

\[g^{\operatorname{SpSt}}_p(X, Y) = \operatorname{tr}\left(X^T\left(I_{2n} - \frac{1}{2}Q_{2n}^Tp(p^Tp)^{-1}p^TQ_{2n}\right)Y(p^Tp)^{-1}\right).\]

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^{\operatorname{SpSt}}(q)$ such that the Cayley Retraction from $p$ along $X$ lands at $q$, i.e. $\mathcal{R}_p(X) = q$ [Bendokat2021].

First, recall the definition the standard symplectic matrix

\[Q = \begin{bmatrix} 0 & I \\ -I & 0 \end{bmatrix}\]

as well as the symplectic inverse of a matrix $A$, $A^{+} = Q^T A^T Q$.

For $p, q โˆˆ \operatorname{SpSt}(2n, 2k, โ„)$ then, we can define the inverse cayley retraction as long as the following matrices exist.

\[ U = (I + p^+ q)^{-1} \in โ„^{2k \times 2k}, \quad V = (I + q^+ p)^{-1} \in โ„^{2k \times 2k}.\]

If that is the case, the inverse cayley retration at $p$ applied to $q$ is

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

ManifoldsBase.manifold_dimension โ€” Method
manifold_dimension(::SymplecticStiefel{n, k})

Returns the dimension of the symplectic Stiefel manifold embedded in $โ„^{2n \times 2k}$, i.e. [Bendokat2021]

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

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

Given a point $p \in \operatorname{SpSt}(2n, 2k)$, project an element $A \in \mathbb{R}^{2n \times 2k}$ onto the tangent space $T_p\operatorname{SpSt}(2n, 2k)$ relative to the euclidean metric of the embedding $\mathbb{R}^{2n \times 2k}$.

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

\[ \operatorname{min}_{X \in \mathbb{R}^{2n \times 2k}} \frac{1}{2}||X - A||^2, \quad \text{s.t.}\; h(X)\colon= X^T Q p + p^T Q X = 0,\]

where $h : \mathbb{R}^{2n \times 2k} \rightarrow \operatorname{skew}(2k)$ defines the restriction of $X$ onto the tangent space $T_p\operatorname{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, computed inplace of q from p along X.

Given a point $p \in \operatorname{SpSt}(2n, 2k)$, every tangent vector $X \in T_p\operatorname{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) \in โ„^{2n \times 2n},\]

as shown in Proposition 3.5 of [Bendokat2021]. Using this representation of $X$, the Cayley retraction on $\operatorname{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 \times 2n$ matrix inverse in the expression above can be reduced down to inverting a $2k \times 2k$ matrix due to Proposition 5.2 of [Bendokat2021].

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}.\]

It is this expression we compute inplace of q.



  • Bendokat2021

    Bendokat, Thomas and Zimmermann, Ralf: The real symplectic Stiefel and Grassmann manifolds: metrics, geodesics and applications arXiv preprint arXiv:2108.12447, 2021 (