{'fraction_non_alphanumeric': 0.059442054687696526, 'fraction_numerical': 0.024098810152692877, 'mean_word_length': 4.244591029023747, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The real symplectic Stiefel and Grassmann manifolds: metrics, geodesics and applications

The real symplectic Stiefel manifold is the manifold of symplectic bases of symplectic subspaces of a fixed dimension. It features in a large variety of applications in physics and engineering. In this work, we study this manifold with the goal of providing theory and matrix-based numerical tools fit for basic data processing. Geodesics are fundamental for data processing. However, these are so far unknown. Pursuing a Lie group approach, we close this gap and derive efficiently computable formulas for the geodesics both with respect to a natural pseudo-Riemannian metric and a novel Riemannian metric. In addition, we provide efficiently computable and invertible retractions. Moreover, we introduce the real symplectic Grassmann manifold, i.e., the manifold of symplectic subspaces. Again, we derive efficient formulas for pseudo-Riemannian and Riemannian geodesics and invertible retractions. The findings are illustrated by numerical experiments, where we consider optimization via gradient descent on both manifolds and compare the proposed methods with the state of the art. In particular, we treat the 'nearest symplectic matrix' problem and the problem of optimal data representation via a low-rank symplectic subspace. The latter task is associated with the problem of finding a 'proper symplectic decomposition', which is important in structure-preserving model order reduction of Hamiltonian systems.

Symplectic structures feature in a large variety of applications in physics and engineering, most prominently in Hamiltonian mechanics [3]. Hamiltonian systems are used in applications ranging from molecular dynamics to celestial mechanics, see [16] and references therein. The symplectic Stiefel manifold is of special importance to optimization problems of the form see [14] and references therein, since it allows to tackle such constrained optimization problems on R 2n×2k as unconstrained optimization problems on SpSt(2n, 2k). Fields of applications include the symplectic eigenvalue problem [7,24] and projection-based structure-preserving model order reduction for Hamiltonian systems. Here, an optimization problem of the form (1.1) appears as the central task of computing a so-called proper symplectic decomposition [2,22,9].
Riemannian optimization requires that we have explicit formulas for essential geometric quantities at hand, as well as efficient algorithms for practical computations.
It is understood that the inner geometry of the symplectic Stiefel manifold depends on the chosen metric.
Related work and state of the art. Optimization on the real symplectic group is considered in [12] with respect to a pseudo-Riemannian metric and in [26,8] with respect to a left-invariant Riemannian metric. Quotients of the real symplectic group relating to the real symplectic Stiefel and Grassmann manifolds are briefly introduced in [23, Subsection 2.1]. To the best of the authors' knowledge, the first treatment of the real symplectic Stiefel manifold with a view on numerical applications is [14]. The optimization algorithm developed there forms the state of the art. The same team of authors compared this method with optimization with respect to a Riemannian metric that stems from a Euclidean metric in the later work [13].
Main original contributions. Starting from the classical real symplectic group, equipped with a bi-invariant pseudo-Riemannian metric, we use a Lie group approach to investigate the real symplectic Stiefel manifold. This original approach allows us to exploit quotient manifold results from semi-Riemannian geometry [21] and enables us to derive the first closed-form expressions for the corresponding pseudo-Riemannian geodesics on SpSt(2n, 2k). Complementary to the pseudo-Riemannian approach, we also introduce a Riemannian metric and derive closed-form expressions for the corresponding Riemannian geodesics. In view of optimization tasks, we provide a formula for the gradient associated with the Riemannian metric and efficiently computable and invertible retractions that approximate the pseudo-Riemannian geodesics.
Moreover, we initiate a study of the manifold of symplectic subspaces that are spanned by symplectic bases, which we term the real symplectic Grassmann manifold in analogy to the classical Stiefel and Grassmann manifolds. Continuing the quotient manifold approach, we derive corresponding formulas for pseudo-Riemannian and Riemannian geodesics and retractions. We promote symplectic subspaces as the main objects of interest in structure-preserving model order reduction of parameterized Hamiltonian systems.
We illustrate the theoretical findings by means of numerical examples. More precisely, we investigate the numerical feasibility of the proposed methods, we tackle the nearest symplectic matrix problem on the real symplectic Stiefel manifold and compute the optimal symplectic subspace containing a given data matrix on the real symplectic Grassmann manifold. The latter problem is directly associated with finding a proper symplectic decomposition in the context of structure preserving model reduction. We juxtapose the methods' performance with the state of the art.
Organisation. Section 2 reviews basic facts on the real symplectic group and states its geodesics associated with a natural bi-invariant pseudo-Riemannian metric and a right-invariant Riemannian metric. In Section 3 we investigate the real symplectic Stiefel manifold, where we cover basic geometry, Riemannian and pseudo-Riemannian metrics and their geodesics as well as the Riemannian gradient. Section 4 introduces the real symplectic Grassmann manifold as a quotient space and provides formulas for the inherited metrics and geodesics and the Riemannian gradient. Suitable retractions fit for replacing the actual geodesics in efficient implementations are given in Section 5. Numerical experiments are contained in Section 6 and Section 7 concludes on the paper.
2. The real symplectic group. Symplectic vectors spaces are the objects of interest for the (local) study of Hamiltonian systems. An introduction can be found in [3]. By definition, a real symplectic vector space is a real vector space V together with a nondegenerate, skew-symmetric bilinear form ω : V × V → R. This means ω is bilinear and fulfills 1. ω(x, y) = 0 for all y ∈ V implies x = 0 (nondegenerate), 2. ω(x, y) = −ω(y, x) (skew-symmetric). As a standard result, such a V is even-dimensional. For any subspace U ⊂ V, the symplectic form ω allows to the define the symplectic complement Since in general U ⊥ ∩ U = {0}, four special cases of subspaces are classified. A subspace U of (V, ω) is called The last case, a symplectic subspace, means that ω restricts to a symplectic form on U, i.e., (U, ω| U ) is a symplectic space in itself.
The linear Darboux theorem [3] states that for any two symplectic vector spaces of the same dimension, there is a linear isomorphism between them preserving the symplectic form. We can therefore restrict our considerations to the standard symplectic vector space (R 2n , ω 0 ), where the standard symplectic form is defined as where n ∈ N and I n is the n × n identity matrix. Note that J T 2n = −J 2n = J −1 2n . The real symplectic group is the matrix Lie group of transformations which leave the standard symplectic form invariant. It has been studied for example in [12,26,8] with a view on applications, and in [3] from a more abstract point of view.
Define for any matrix A ∈ R 2n×2k the symplectic inverse [22] A + := J T 2k A T J 2n .
The real symplectic group is then defined as For any M ∈ Sp(2n, R), and any x, y ∈ R 2n , it holds that The corresponding Lie algebra is given by the Hamiltonian matrices Accordingly, the tangent space of the real symplectic group at M is given by translation by M , i.e.

Pseudo-Riemannian metric.
Similarly to [12], we define a bi-invariant pseudo-Riemannian metric on Sp(2n, R) by h M : The factor 1 2 is introduced for convenience. If X 1 = M Ω 1 and X 2 = M Ω 2 , then i.e. ·, · M is exactly − 1 2 times the pseudo-Riemannian metric defined in [12] and therefore − 1 2 times the Khvedelidze-Mladenov metric [17] on the general linear group. By properties of the trace and the symplectic inverse, it can be immediately verified that the pseudo-Riemannian metric defined in this way is bi-invariant. Therefore, making use of [21,Proposition 11.9], the (pseudo-Riemannian) geodesics are given by the one-parameter subgroups , we introduce a right-invariant Riemannian metric in anticipation of the quotient structure that is considered in the upcoming Section 3. We also derive the corresponding gradient and geodesics.
The mapping g M : defines point-wise a right-invariant Riemannian metric on the real symplectic group Sp(2n, R). The right-invariance follows from the fact that for every N ∈ Sp(2n, R), . The Riemannian gradient for this metric is given as follows: Let f : Sp(2n, R) → R be differentiable and let ∇f M be the Euclidean gradient of a continuous extension of f to an open subset of R 2n×2n around M ∈ Sp(2n, R), evaluated at M . Then the Riemannian gradient of f at M (with respect to the metric g M ) is This follows from the fact that by definition grad g f (M ) is the unique tangent vector at M such that g M (grad g f (M ), X) = df M (X) = tr((∇f M ) T X) holds for all X ∈ T M Sp(2n, R). By making use of the fact that M + X = −X + M , grad g f (M ) solves this equation, and grad g f (M ) ∈ T M (Sp(2n, R)) follows from grad g f (M )M + = −M (grad g f (M )) + . We can derive the Riemannian geodesics corresponding to the Riemannian metric (2.2) analogously to [25,Proposition 4.2], where the Riemannian geodesics corresponding to a right-invariant metric on the general linear group GL(n) were derived. Proposition 2.1. Let M ∈ Sp(2n, R) and X ∈ T M Sp(2n, R). The Riemannian geodesic γ with γ(0) = M andγ(0) = X for the Riemannian metric (2.2) is given by [25,Proposition 4.2] can be transferred straightforwardly to this setting.
3. The real symplectic Stiefel manifold. The real symplectic Stiefel manifold is defined as It contains the matrices U ∈ R 2n×2k , whose column vectors form symplectic bases for the 2k-dimensional symplectic subspaces of (R 2n , ω 0 ) and was treated in [14,13,24]. Note the formal similarity with the (compact) Stiefel manifold St(n, k) = U ∈ R n×k U T U = I k . As a novelty, and in contrast to the aforementioned references, we will pursue a Lie group-based approach to study the real symplectic Stiefel manifold. We will furthermore introduce a new pseudo-Riemannian and a new Riemannian metric and derive the geodesics for both.
Denote the projection onto the first k columns of a matrix, when multiplied from the right, by and the projection onto the first k and the (n + 1)th to the (n + k)th column by Our first goal is to recognize the real symplectic Stiefel manifold as a quotient of the real symplectic group. To this end, we introduce the following canonical projection: Proposition 3.1. The real symplectic Stiefel manifold is diffeomorphic to the quotient It has dimension dim(SpSt(2n, 2k)) = (4n − 2k + 1)k.
Proof. The set SpSt(2n, 2k) is the orbit of E under the group action of Sp(2n, R) that is induced by left-multiplication. The stabilizer of this group action is isomorphic to Sp(2(n − k), R). From [19,Theorem 21.20], it follows that SpSt(2n, 2k) has a unique smooth manifold structure for which the group action is smooth. It furthermore follows that the dimension of the real symplectic Stiefel manifold is dim(SpSt(2n, 2k)) = dim(Sp(2n, R)) − dim(Sp(2(n − k), R)) = (4n − 2k + 1)k, Fig. 1. Visualization of the quotient structure of the real symplectic Grassmann and Stiefel manifold with respect to the symplectic group. Any point P ∈ SpGr(2n, 2k) has an equivalence class in SpSt(2n, 2k) as its pre-image under ρ, visualized by the blue line through U ∈ SpSt(2n, 2k). This equivalence class in turn has again an equivalence class of symplectic matrices in Sp(2n, R) as its pre-image under π, visualized by the blue area around M . The equivalence class of a single point U ∈ SpSt(2n, 2k) is of lower dimension, visualized by the red line throughM ∈ Sp(2n, R).
The quotient manifold structure of the real symplectic Stiefel manifold (and of the real symplectic Grassmann manifold, which is to be discussed later on) with the symplectic group as the associated total space is visualized in Figure 1.
Note that the tangent space parametrization [14, (3.8b)] is similar to (3.3), but the chosen complement there is not necessarily a symplectic complement.
3.1. Pseudo-Riemannian metric on SpSt(2n, 2k). According to our quotient Lie group approach, SpSt(2n, 2k) inherits a pseudo-Riemannian metric from the pseudo-Riemannian metric (2.1) on the total space Sp(2n, R) in a natural way, by making use of horizontal lifts. A big advantage of this construction is that the corresponding geodesics can then be obtained via the projection of horizontal geodesics, i.e., geodesics with horizontal tangent vectors on Sp(2n, R) [21,Corollary 7.46].
Eventually, we will exploit abstract results of semi-Riemannian geometry [21, §11] for determining the geodesics. To enable this, we show next that SpSt(2n, 2k) ∼ = Sp(2n, R)/ stab E is reductive, and naturally reductive with respect to the pseudo-Riemannian metric h of (2.1), see [21,Definition 11.21 & 11.23] for an explanation of these terms.
Lemma 3.4. The real symplectic Stiefel manifold SpSt(2n, 2k) is reductive.With respect to the pseudo-Riemannian metric h, it is naturally reductive.
Proof. In view of (3.4) and the fact that Ver π sp(2n, R) is isomorphic to the Lie algebra of Sp(2(n − k), R) ∼ = stab E , we need to show that the complementary subspace Hor π,h sp(2n, R) is Ad(stab E )-invariant in order to establish reductiveness, where Ad M (Ω) = M ΩM + , see [21, §11, p. 303]. In fact, for every M ∈ stab E (and therefore M + ∈ stab E ) and every Ω ∈ Hor π,h sp(2n, R), it holds that By Proposition (3.3), it follows that Ad M (Ω) ∈ Hor π,h sp(2n, R), which means that The fact that SpSt(2n, 2k) is naturally reductive with respect to h follows from a direct calculation, by making use of the fact that the projection from sp(2n, R) onto Recall that the Lie algebra sp(2n, R) = T I Sp(2n, R) is the tangent space at the identity. By left translation, every tangent space T M Sp(2n, R) can be split into a vertical and horizontal part, and U = M E. It follows therefore that the application of dπ M to Hor π,h M Sp(2n, R) gives (3.3). This implies that T U SpSt(2n, 2k) can also be parameterized as , Ω ∈ Hor π,h sp(2n, R) .

Proposition 3.5 (Alternative tangent vector parameterization). Every tangent vector
with Ω ∈ Hor π,h sp(2n, R), is unique. It can be calculated via Proof. This can be seen by making use of U = M E and Proposition 3.3.
Via horizontal lifts, a pseudo-Riemannian metric on the real symplectic Stiefel manifold can be defined as follows: For two tangent vectors ∆ 1 , ∆ 2 ∈ T U SpSt(2n, 2k) and U = M E, calculateΩ(U, ∆ i ) according to (3.7), i = 1, 2. The horizontal lift to Hor π,h M Sp(2n, R) is then given by where and the fact that Ω i fulfills Proposition 3.3. In the following, we exploit that pseudo-Riemannian submersions [21,Definition 7.44] are particularly useful to find the geodesics in a quotient space, given that the geodesics in the associated total space are known, see [21,Corollary 7.46].
With respect to the metric h SpSt U , π is a pseudo-Riemannian submersion.
Proof. A direct calculation shows the identities of (3.8). Since h is a pseudo-Riemannian metric, h SpSt U is a pseudo-Riemannian metric as well. The projection π is then a pseudo-Riemannian submersion by Lemma 3.4 and [21, Lemma 11.24].
The geodesics, calculated via the exponential mapping with respect to h SpSt U , can be found via the projection of the exponential mapping in Sp(2n, R).
In the form of (3.9), the exponential mapping depends on the matrix exponential of a 2n × 2n matrix. For tangent vectors ∆ ∈ T U SpSt(2n, 2k) with ∆ + (I 2n − U U + )∆ invertible, we can reduce the computational complexity to 4k × 4k. This is rendered possible by the fact that for X, Y ∈ R n×k with Y T X ∈ R k×k non-singular, we have from [10, Prop. 3] that Proof. We start from (3.9) withΩ(U, ∆) according to (3.7). Introducing Note that for the calculation of (3.11) we don't need the invertibility of H + H, and one can check that the right hand side is always an element in the symplectic Stiefel manifold SpSt(2n, 2k). We can therefore always apply (3.11) to calculate a curve.
Remark 3.9. The simplified formula for the symplectic Stiefel exponential (3.11) is similar to the so-called quasigeodesic retraction defined in [14, Lemma 5.1], which in our notation is given as The two curves are however not identical. Note also the structural similarity with the formula for the Euclidean Stiefel geodesics of [11, Section 2.2.2].

Right-invariant
Riemannian metric on SpSt(2n, 2k). The real symplectic Stiefel manifold may be equipped with different Riemannian metrics. The so-called canonical-like metric has been studied in [14], while [13] considers a restriction of the Euclidean metric. To the best of the authors' knowledge, the geodesics for these metrics are unknown. Complementary to the aforementioned Riemannian metrics, we use the Riemannian metric g M of (2.2) on Sp(2n, R) to introduce a third Riemannian metric on SpSt(2n, 2k) via a horizontal lift, which allows us to find the corresponding geodesics. This metric is invariant under the group action of Sp(2k, R) from the right and induces therefore a Riemannian metric on the symplectic subspaces that will be considered in Subsection 4.2.
We begin by splitting the tangent space T M Sp(2n, R) at M ∈ Sp(2n, R) into a vertical and a horizontal part with respect to g M and the projection π from (3.2), As the vertical part is defined as the kernel of dπ M , it is the same as in (3.6). The horizontal part, however, is different, since it is now given as the orthogonal complement of Ver π M Sp(2n, R) with respect to the metric g M of (2.2). This yields Here, P = J T 2n U U + J 2n and U = π(M ) = M E. Equation (3.13) can be established as follows: Any horizontal tangent vector X =ΩM , withΩ ∈ sp(2n, R), fulfills for all DefineΩ := M + J 2nΩ J 2n M . Then 0 = 1 2 tr(ΩΩ) for all Ω ∈ Ver π sp(2n, R) implies thatΩ fulfills (3.5). By making use of E = M + M E = M + U , it follows that Conversely, ifΩ fulfills the above equation, then it follows thatΩ := M + J 2nΩ J 2n M fulfills (3.5) and thereforeΩM ∈ Hor π,g M Sp(2n, R). As usual, for U = π(M ) we can identify the tangent space T U SpSt(2n, 2k) with the horizontal space Hor π,g M Sp(2n, R). Any ∆ ∈ T U SpSt(2n, 2k) is of the form ∆ = π(ΩM ) =ΩU for someΩM ∈ Hor π,g M Sp(2n, R), and we can find the horizontal lift This follows from the facts that , where the last equation follows from a straightforward calculation.
The right action of Sp(2(n−k), It is furthermore isometric, by right-invariance of (2.2). From [20,Theorem 2.28], it follows that there is a unique Riemannian metric on SpSt(2n, 2k) such that π is a Riemannian submersion. This Riemannian metric is given via the horizontal lift.
Proposition 3.10. The Riemannian metric on SpSt(2n, 2k), for which π is a Riemannian submersion, is right-invariant and given point-wise by Proof. The Riemannian submersion property and right-invariance hold by the definition of g SpSt U via the horizontal lift. The second equality follows from the combination of (2.2), (3.14) and (3.15).
The Riemannian gradient of a function f : SpSt(2n, 2k) → R with respect to g SpSt is given by where ∇f (U ) denotes the Euclidean gradient of a smooth extension of f around U ∈ SpSt(2n, 2k) in R 2n×2k at U . This holds because U + grad g f (U ) = −(grad g f (U )) + U , which implies grad g f (U ) ∈ T U SpSt(2n, 2k), and because grad g f (U ) solves for all ∆ ∈ T U SpSt(2n, 2k). Riemannian geodesics on Sp(2n, R) with a horizontal tangent vector at every point project to Riemannian geodesics on SpSt(2n, 2k) by [21,Corollary 7.46]. (Mind that the referenced result is stated in the pseudo-Riemannian setting, but also holds true in the Riemannian case.) We show that Riemannian geodesics on Sp(2n, R) with initial horizontal tangent vector have a horizontal tangent vector throughout.
We are now ready to state the Riemannian geodesics on SpSt(2n, 2k) with respect to the Riemannian metric g SpSt from (3.16).
Proof. This follows directly from the preceding discussion and the definition of the horizontal lift.
Equation (3.18) is formulated with 2n × 2n-matrices, but may in practical calculations be reduced to work with tall, skinny 2n × 8k matrices and matrix exponentials of an 8k × 8k and a 4k × 4k matrix, respectively. To this end, defineĀ ∈ so(2k, R) bȳ With∆ := UĀ +H ∈ T U SpSt(2n, 2k) it holds thatΩ(∆) = Y X T , where This follows from (3.15) and solvingΩ(∆) T = (I 2n − 1 2 U U + )∆U + −U∆ + (I 2n − 1 2 U U + ) for∆, i.e.,Ā = U +Ω (∆) T U andH = (I 2n − U U + )Ω(∆) T U . Furthermore, definê Proposition 3.13. With notation as above, it holds that Proof. First, note that sinceΩ(∆) = Y X T , it holds that exp m (Ω(∆) T )U = exp m (XY T )U . We make use of (3.10), which implies 4. The real symplectic Grassmann manifold. Similar to the usual Grassmann manifold [4,11] of linear subspaces of a fixed dimension, we define the real symplectic Grassmann manifold as the manifold of symplectic subspaces of dimension 2k of (R 2n , ω 0 ). This must not be confused with the Lagrangian Grassmannian, the manifold of Lagrangian subspaces, which is also referred to as the symplectic Grassmann manifold by some authors. The quotient manifold approach which we use is similar to the course of action in [5]. As in the case of linear subspaces, we identify a symplectic subspace with the associated symplectic projection onto it.   Because every U ∈ SpSt(2n, 2k) has a representation U = M E, M ∈ Sp(2n, R), it is sufficient to show that every P ∈ SpGr(2n, 2k) is equal to P = U U + for some U ∈ SpSt(2n, 2k). This fact is established as follows: Since P J 2n P T is skew-symmetric, it features a 'Schur-like decomposition' [27, eq. (5)] of the form where Q ∈ O(2n) is a real orthogonal matrix. Moreover Σ = diag(σ 1 , . . . , σ k ), where σ i > 0 for all i = 1, . . . , k, because rank(P ) = 2k [27,Proposition 3]. From P + = P = , with I 2n,2k as in (3.1), it furthermore holds that P = U U + . The fact that U + U = I 2k , i.e., U ∈ SpSt(2n, 2k), follows from P 2 = P . The other inclusion, i.e., φ(M, E 0 ) ∈ SpGr(2n, 2k) for all M ∈ Sp(2n, R) is immediate. The stabilizer of the group action φ(·, E 0 ) is given by The manifold structure now follows from [19,Theorem 21.20]. The real symplectic Grassmann manifold is by [19,Theorem 21.18] diffeomorphic to the homogeneous space The dimension of SpGr(2n, 2k) is obtained via the standard formula Note that the real symplectic Grassmann manifold SpGr(2n, 2k) has the same dimension as the Grassmann manifold Gr(2n, 2k). The manifolds are not the same however, since not every 2k dimensional subspace of R 2n is also a symplectic subspace of the standard symplectic space (R 2n , ω 0 ).
Let M ∈ Sp(2n, R) such that P = M E 0 M + and define U = M E ∈ SpSt(2n, 2k) and U s = M E s . Then P = U U + , and it follows from the preceding proposition that T P SpGr(2n, 2k) can be identified with Hor ρ,h U SpSt(2n, 2k) via the horizontal lift. For Γ ∈ T P SpGr(2n, 2k), this horizontal lift is explicitly given by Γ hor U = ΓU ∈ Hor ρ,h U SpSt(2n, 2k), as can be seen by the fact that there is Ω ∈ sp(2n, R) such that Γ = [Ω, U U + ] and By making use of the horizontal lift to SpSt(2n, 2k), we can define a pseudo-Riemannian metric on the real symplectic Grassmann manifold SpGr(2n, 2k). Proof. Similar to Proposition 3.6.
The connection between SpSt(2n, 2k) and SpGr(2n, 2k) allows us to state the following decomposition of real symplectic Stiefel matrices.
As the metric g SpGr of (4.6) is defined via a horizontal lift, we obtain the associated geodesics by projection.

Riemannian metric on
Again, we split the tangent space T U SpSt(2n, 2k) at U ∈ SpSt(2n, 2k) into a vertical part with respect to ρ and a horizontal part with respect to ρ and g SpSt U from (3.16). The former yields (4.3), as the vertical space is independent of the metric. The latter gives This follows from g SpSt U U A, (U H + − HU + ) T U = 0 for all U A ∈ Ver ρ U SpSt(2n, 2k) and all H ∈ R 2n×2k with U + H = 0, and by counting degrees of freedom. For any ∆ = (U H + − HU + ) T U ∈ Hor ρ,g U SpSt(2n, 2k), the corresponding H (which is not to be confused with (I 2n − U U + )∆ here) can by calculated via We can identify the horizontal space Hor ρ,g U SpSt(2n, 2k) with the tangent space T ρ(U ) SpGr(2n, 2k) and define a Riemannian metric on SpGr(2n, 2k) via the restriction of the Riemannian metric g SpSt of (3.16) to the horizontal spaces.
For Γ i ∈ T U U + SpGr(2n, 2k), let (Γ i ) hor,g U ∈ Hor ρ,g U SpSt(2n, 2k) be the horizontal lift to Hor ρ,g U SpSt(2n, 2k), i.e. dρ U ((Γ i ) hor,g U ) = Γ i . The mapping defines pointwise a Riemannian metric. We are not aware of an explicit mapping to calculate the horizontal lift with respect to g SpSt for a given Γ ∈ T U U + SpGr(2n, 2k). Nevertheless, one can directly work with symplectic Stiefel representatives and horizontal tangent vectors, i.e., with 2n × 2k-matrices.   Sp(2k, R). We assume that f can (locally) be extended to a smooth function on 2n × 2k-matrices, for convenience again denoted by f . The Riemannian gradient of f with respect to g SpGr is given by where ∇f U denotes the Euclidean gradient of a smooth extension of f around U in R 2n×2k . This follows from [1, Equation (3.39)] and g SpSt U (grad g f (U ), ∆) = df U (∆) = tr((∇f U ) T ∆) for all ∆ ∈ Hor ρ,g U SpSt(2n, 2k), as well as the fact that grad g f (U ) ∈ Hor ρ,g U SpSt(2n, 2k). Proposition 4.9. Let U ∈ SpSt(2n, 2k) and ∆ ∈ Hor ρ,g U SpSt(2n, 2k). The lifted symplectic Grassmann geodesic from U in direction ∆ is given by whereΩ is given by (3.15).
Proof. We need to show that the tangent vector d dt Exp SpGr U (t∆) is horizontal for every t. Then, the claim follows from [21,Cor. 7.46]. Since ∆ ∈ Hor ρ,g U SpSt(2n, 2k) is equivalent to J 2n U J T 2k U TΩ (∆)J 2n U J T 2k U T = 0, the proof follows in the same fashion as the one of Lemma 3.11.

Retractions and computational issues.
Calculating the matrix exponential of an n × n matrix is computationally expensive if n is large. Furthermore, numerical experiments show that even though the matrix exponential of a Hamiltonian matrix is theoretically guaranteed to yield a symplectic matrix as an output, this is not necessarily the case in practice, where one needs to rely on numerical tools to compute the standard matrix exponential. While there are specialized algorithms for the matrix exponential of a Hamiltonian matrix [18], there is another alternative: The Cayley map. In this section, we propose the use of the Cayley map for approximating the pseudo-Riemannian geodesics in order to define retractions on the symplectic Stiefel and Grassmann manifold. Furthermore, these retractions turn out to be invertible in closed form on both manifolds, which can for example be used for interpolation and optimization purposes and for defining local coordinates. In the experiments of Section 6, the Cayley-based retraction turns out to be computationally cheaper and to retain the manifold structure to a much higher numerical accuracy.
A retraction [1] on a smooth manifold M with tangent bundle T M is a smooth mapping R : T M → M such that for any x ∈ M , The Cayley transformation cay(X) := (I n + X)(I n − X) −1 , X ∈ R n×n is widely used as a standard approximation of the matrix exponential exp m (2X). Of special interest in the present context is the property that cay maps from sp(2n, R) to Sp(2n, R) [3]. This was also exploited in [12]. The inverse of the Cayley transform is given by [3] cay −1 (M ) = (M − I n )(I n + M ) −1 .

Cayley retraction on the real symplectic Stiefel manifold.
Replacing the matrix exponential in the pseudo-Riemannian exponential (3.9) on the symplectic Stiefel manifold with the Cayley transform leads to the Cayley retraction defined in [14,Definition 5.2]. Yet note that the Cayley retraction in the aforementioned reference was found unaware of the pseudo-Riemannian geodesics by transferring the Cayley retraction on the classical Stiefel manifold St(n, k) to the symplectic case.
In [14,Proposition 5.5], it was proposed to use the Sherman-Morrison-Woodbury formula where A ∈ R n×n , X, Y ∈ R n×k , to reduce the matrix inverse in (5.1) from 2n × 2n to 4k × 4k. We show that we can even reduce it to a matrix inversion of dimensions 2k × 2k.
Proposition 5.2. Let U ∈ SpSt(2n, 2k) and ∆ ∈ T U SpSt(2n, 2k). Define A := U + ∆ and H := ∆ − U A. Then Proof. For t = 0, the equality holds. In the following, assume t = 0. Similarly to the proof of Proposition 3.8, define X = 1 Then againΩ(U, ∆) = XY T as in (3.7), and It holds that Block-matrix inversion via the Schur complement yields Putting everything together, we obtain which shows the claim.
Unlike the pseudo-Riemannian exponential (3.11) or the Riemannian exponential (3.19), we can invert the Cayley retraction (5.2) in closed form. Apart from interpolation, this facilitates the calculation of local coordinates on SpSt(2n, 2k). SpSt(2n, 2k). If (I 2k +U + V ) −1 and (I 2k +V + U ) −1 exist, it holds for and it holds that which shows the claim.

Cayley retraction on the real symplectic Grassmann manifold.
With the quotient manifold approach to the symplectic Stiefel manifold and the definition of the symplectic Grassmann manifold, we can show an additional property of R SpSt : It maps horizontal tangent vectors (with respect to the pseudo-Riemannian metric h SpSt from (3.8)) to curves with horizontal tangent vectors everywhere. We can therefore use it to calculate approximations of the pseudo-Riemannian symplectic Grassmann geodesics lifted to the symplectic Stiefel manifold.
which implies the claimed result.
We can also directly invert R SpGr P (Γ) on symplectic Stiefel representatives.
Proposition 5.7. Let U, V ∈ SpSt(2n, 2k). If and are well-defined, it holds that The difference between the connecting curves from Proposition 5.3 and Proposition 5.7 is visualized in Figure 2.
6. Numerical Experiments. In this section, we study the feasibility of different retractions on SpSt(2n, 2k) and investigate optimization problems via gradient descent on SpSt(2n, 2k) and SpGr(2n, 2k), respectively. All experiments are conducted with MATLAB version R2021a on a laptop with Ubuntu 18.04, an Intel® Core™ i7-8850H CPU and 16GB RAM. We generate random Hamiltonian matrices via Ω = A B C −A T , where A, B and C are generated by randn(n,n), and then B and C are symmetrized. For reproducability, all random matrices are constructed with the random stream s = RandStream('mt19937ar'). 6.1. Feasibility of different retractions on the symplectic Stiefel manifold. We compare the numerical feasibility of the Riemannian geodesic (3.19) with respect to g SpSt , the Cayley-retraction (5.2), the pseudo-Riemannian geodesic (3.11) and the quasi-geodesic retraction (3.12) from [14]. To this end, we generate a (pseudo) random point on SpSt(2n, 2k) via U = cay(Ω)E, where Ω ∈ sp(2n, R) is scaled to ||Ω|| F = 1. We furthermore generate a (pseudo) random tangent vector ∆ ∈ T U SpSt(2n, 2k), also scaled to ||∆|| F = 1. For the chosen retractions R, we calculate U (t) = R U (t∆) with t ∈ [0, 10 3 ] and plot the feasibility ||U (t) + U (t) − I 2n || F . The average over 10 runs is shown in Figure 3 for n = 1000, k = 20 (left) and n = 1000, k = 200 (right). It can be seen that the Riemannian geodesic, the pseudo-Riemannian geodesic and the quasi-geodesic retraction, which all rely on the matrix exponential, fail numerically to stay on SpSt(2n, 2k) for tangent vectors of Frobenius-norm larger than O(10 2 ). The Cayley-Retraction, while less feasible at some points, fulfills the manifold condition up to an error of about 10 −8 for tangent vectors of any tested size on SpSt(2000, 40) and up to an error of about 10 −4 on SpSt(2000, 400).
6.2. Gradient descent on the real symplectic Stiefel manifold. We tackle an academic instance of the 'nearest symplectic matrix' problem min U ∈SpSt(2n,2k) ||U − A|| 2 F via a Riemannian gradient descent. For this, we set A = randn(2*n,2*k) and then normalize A = A/norm(A,2), as in [14]. The initial point for starting the optimization procedure is set to be U0 = cay(X/2)E, where X ∈ sp(2n, R) is a random 2n × 2n Hamiltonian matrix, scaled by X = X/norm(X,'fro'). As gradient descent algorithm we use [14,Alg. 1] with monotone line search and stopping criterion [14, Eqs. (6.1) and (6.2)]. For the reader's convenience, we restate the procedure here as Algorithm 6.1 in the precise form in which we use it. As the trial step size γ k , we use the alternating BB method γ ABB We compare gradient descent for the Riemannian metric g SpSt with geodesic step- ping and Cayley stepping, respectively, with gradient descent from [14] with Cayley stepping. For the gradient descent according to [14], we choose the optimal settings stated in this reference, i.e., the canonical-like metric g ρ with ρ = 1 2 and gradient (I), according to [14,Subsection 6.2.2]. In the actual implementation of all methods included in this comparison, care has been taken that the action of large matrices like J 2n and I 2n is applied directly, so that these matrices are never formed explicitly. Figure 4 displays the objective function value versus the iteration count (left) and the convergence history according to the gradient norm (right), respectively. For comparison purposes, all methods are run for a fixed number of 60 iterations. It can be seen that the algorithms deliver similar results in regard of the convergence by iterations, depending on the chosen tolerance. The run time however differs: In Table 1, we compare the three methods and state the average iterations and run time until numerical convergence over 10 runs. We furthermore denote the relative deviation from the respective minimum over all three methods after convergence. It can be seen that for SpSt (2000,40) and SpSt(2000, 400), gradient descent with Cayley stepping is the fastest method regarding run time, while Geodesic descent is the slowest. For SpSt(2000, 400), the run time for geodesic stepping increases drastically, since (3.19) requires the matrix exponential of both a 8k × 8k and a 4k × 4k matrix.
In Figure 5, we compare the convergence over time for one optimizer run on SpSt (2000, 40). For each step, the run time is measured over one full iteration of the outer for-loop in lines 1 to 24 in Algorithm 6.1. It can be seen that gradient descent with respect to the Riemannian metric g SpSt with Cayley stepping converges the fastest in terms of the run time. The iteration count for Cayley and geodesic stepping with respect to g SpSt are comparable.

12:
for h min ≤ h ≤ h max do 13: Break as converged.

22:
end if 23: end if 24: end for Output: Iterates U k . Table 1 Numerical performance for the cases considered in Section 6.2, taking averages over 10 runs. The (pseudo-)random data are generated for n = 1000, and either k = 20 or k = 200, respectively. The minimum is the respective minimum over all three methods after convergence.   for representing a given data matrix S ∈ R 2n×2n , i.e., This problem is associated with computing a proper symplectic decomposition, a task that is central in Hamiltonian model order reduction [22]. Here, we work in an academic setting, where the target matrix S is generated as a random symplectic 2) Analogous to Figure 5 but with the target matrix A scaled to ||A|| 2 = 2. Here, gρ denotes the canonical-like metric from [14] with Cayley stepping.
subspace representative plus an error term, i.e.
where A ∈ SpSt(2n, 2k) is a random symplectic Stiefel matrix found in the same manner as the initial point U 0 , and E is a random n × n-matrix, divided by its 2-norm. The parameters for the gradient descent algorithm are the same as in Subsection 6.2. The resulting average of the function value and the convergence history over 10 runs with a fixed number of 40 iterations is shown in Figure 7. It can be seen that gradient descent for all methods produces similar results in regards of the iteration count. For the gradient descent from [14] and for the gradient descent according to g SpSt with Cayley stepping, we ignore the quotient structure and treat (6.1) as a minimization problem on SpSt(2n, 2k). The run time and iteration count of the methods is compared in Table 2, similarly to Subsection 6.2. We also compare the convergence over run time for a single optimizer run in Figure 8. It can be seen that gradient descent with Cayley stepping according to g SpSt or g SpGr converges fastest and both methods perform comparable to the method of [14]. As is to be expected, for k = 200, geodesic stepping is again considerably slower. Note however that for all methods, processing the 2n × 2n input matrix S requires a high base level run time.
Remark 6.1. It is also possible to tackle optimization problems with pseudo-Riemannian methods [12,15]. For our experiments however, we achieved better result with the Riemannian methods. Nevertheless, pseudo-Riemannian optimization might prove beneficial in some settings.

Conclusion.
We introduced a novel pseudo-Riemannian framework for the real symplectic Stiefel manifold SpSt(2n, 2k). In analogy to the classical Stiefel and Grassmann manifolds, we introduced the real symplectic Grassmann manifold SpGr(2n, 2k). For a natural pseudo-Riemannian metric, we derived the corresponding geodesics. With the formulas at hand, we explained the Cayley retraction as an Here, gρ denotes the canonical-like metric from [14] with Cayley stepping.  approximation of the pseudo-Riemannian geodesics and found an efficiently computable expression for the retraction, which turned out to be invertible in closed form.
Secondly, we introduced a new Riemannian framework for both SpSt(2n, 2k) and SpGr(2n, 2k), coming from a right-invariant Riemannian metric on Sp(2n, R), and derived the corresponding Riemannian geodesics. Since to the best of the authors' knowledge, the Riemannian geodesics for no other Riemannian metric on SpSt(2n, 2k) are known, this opens up new possibilities for theoretical studies and applications.
In the experiments, we showed that gradient descent with the Riemannian geode- Table 2 Numerical performance for the cases considered in Section 6.3, taking averages over 10 runs. The (pseudo-)random data are generated for n = 1000, and either k = 20 or k = 200, respectively. The minimum is the respective minimum over all five methods after convergence. sics or optimized Cayley retraction outperforms the state-of-the-art method from [14] in some cases and delivers comparable results in others. Cayley stepping with respect to the Riemannian metric g SpSt on SpSt(2n, 2k) converges in general the fastest among all methods, regarding the run time.
The invertible retractions provide local coordinates on the manifolds SpSt(2n, 2k) and SpGr(2n, 2k), respectively. This renders it possible to apply tangent space methods, e.g. for interpolation purposes. A potential area of application of such tangent space interpolation is parametric model order reduction of Hamiltonian systems. The proposed coordinate transformations allow to approach this problem analogously to parametric model order reduction of general dynamical systems [6,28].