# Optimization on the symplectic Stiefel manifold:
SR decomposition-based retraction and applications

Bin Gao  Nguyen Thanh Son<EMAIL_ADDRESS>Tatjana Stykel  State Key
Laboratory of Scientific and Engineering Computing, Academy of Mathematics and
Systems Science, Chinese Academy of Sciences, 100190 Beijing, China  Institut
für Mathematik $\,\&$ Centre for Advanced Analytics and Predictive Sciences,
Universität Augsburg, Universitätsstraße 12a, 86159 Augsburg, Germany
and Thai Nguyen University of Sciences, 24118 Thai Nguyen, Vietnam  Institut
für Mathematik $\,\&$ Centre for Advanced Analytics and Predictive Sciences,
Universität Augsburg, Universitätsstraße 12a, 86159 Augsburg, Germany

###### Abstract

Numerous problems in optics, quantum physics, stability analysis, and control
of dynamical systems can be brought to an optimization problem with matrix
variable subjected to the symplecticity constraint. As this constraint nicely
forms a so-called symplectic Stiefel manifold, Riemannian optimization is
preferred, because one can borrow ideas from unconstrained optimization
methods after preparing necessary geometric tools. Retraction is arguably the
most important one which decides the way iterates are updated given a search
direction. Two retractions have been constructed so far: one relies on the
Cayley transform and the other is designed using quasi-geodesic curves. In
this paper, we propose a new retraction which is based on an SR matrix
decomposition. We prove that its domain contains the open unit ball which is
essential in proving the global convergence of the associated gradient-based
optimization algorithm. Moreover, we consider three applications—symplectic
target matrix problem, symplectic eigenvalue computation, and symplectic model
reduction of Hamiltonian systems—with various examples. The extensive
numerical comparisons reveal the strengths of the proposed optimization
algorithm.

###### keywords:

Symplectic Stiefel manifold , Riemannian optimization, retraction, SR
decomposition, symplectic target problem, symplectic eigenvalue, Hamiltonian
systems, symplectic model reduction

###### MSC:

[2020] 15A23 , 32C25 , 65F15 , 65F99 , 65K05 , 65P10 , 90C30

##  1 Introduction

In scientific computing and physics, one often has to work with structured
matrices. One type of them is the set of symplectic matrices defined as

${\mathrm{Sp}(2k,2n)}:=\\{X\in\mathbb{R}^{2n\times 2k}\quad:\quad
X^{T}J_{2n}X=J_{2k}\\},$

where $J_{2n}=\left[\begin{smallmatrix}0&I_{n}\\\
-I_{n}&0\end{smallmatrix}\right]$ is the $2n\times 2n$ skew-symmetric Poisson
matrix and $I_{n}$ denotes the $n\times n$ identity matrix. This set has been
indicated to be a closed, unbounded, embedded submanifold of
${\mathbb{R}^{2n\times 2k}}$, termed as the _symplectic Stiefel manifold_ [1].
For $k=n$, this manifold forms the symplectic Lie group denoted by
${\mathrm{Sp}(2n)}$. Another view on the symplectic Stiefel manifold is to
consider it as a quotient manifold of the two Lie groups
${\mathrm{Sp}(2n)}/\mathrm{Sp}(2n-2k)$, see [2] for details. In this paper, we
consider a minimization problem with the symplecticity constraint given by

$\min_{X\in{\mathrm{Sp}(2k,2n)}}f(X)$ (1)

with a continuously differentiable cost function $f$. In order to solve such
an equality-constrained problem, one can apply, for example, a standard
penalty method or an augmented Lagrangian approach [3, Chapter 17]. However,
by exploiting the rich structure of ${\mathrm{Sp}(2k,2n)}$, Riemannian
optimization becomes preferred as in this setting, various well-known
unconstrained optimization methods in Euclidean spaces can be extended to the
case of nonlinear manifold.

To pursue this direction, a Riemannian structure of the symplectic Stiefel
manifold ${\mathrm{Sp}(2k,2n)}$ and necessary geometric tools are required.
The tangent spaces of ${\mathrm{Sp}(2k,2n)}$ have been completely
characterized in [1]. Moreover, a canonical-like metric and an Euclidean
metric have been introduced in [1] and [4], respectively, based on which the
normal spaces, the orthogonal projections onto the tangent and normal spaces
and the Riemannian gradients have been investigated. In order to perform a
search on ${\mathrm{Sp}(2k,2n)}$, we need a mapping that retracts the
Riemannian gradient from the tangent space at the current point to the
manifold. To this end, two retractions have been proposed: the first one is
based on the Cayley transform, henceforth referred to as the _Cayley
retraction_ , and the other is constructed to mimic the shape of a geodesic,
called the _quasi-geodesic retraction_. Based on these retractions, a non-
monotone line search algorithm has been developed in [1] and various tests
have been performed there to validate it. An equivalent expression for the
Cayley retraction on ${\mathrm{Sp}(2k,2n)}$ has also been derived in [2]
starting from a pseudo-Riemannian exponential.

In this paper, we propose a new retraction which is based on an SR matrix
decomposition [5, 6, 7]. Such a retraction is therefore referred to as the _SR
retraction_. We also establish in Proposition 3.5 and Theorem 3.7 that for
rectangular $2n\times 2k$ matrices, the SR decomposition almost always exists.
More importantly, we show in Theorem 3.8 that the update can be computed if
the length of the update with respect to the spectral matrix norm is smaller
than one. Like the existing retractions, the new one can be combined with any
metric which results in new Riemannian optimization schemes whose global
convergence can be shown in a way similar to [1, Theorem 5.7].

The minimization problem (1) appears in various areas of physics and
scientific computing. For example, the task of averaging optical transference
matrices, investigating beam dynamics, and optimal control of quantum gate can
be formulated as the minimization problem (1) with $k=n$, see, e.g., [8, 9,
10]. Further, stability analysis of weakly damped gyroscopic systems [11] can
be performed by solving (1) with the trace cost function which enables
determining so-called symplectic eigenvalues [12, 13, 14]. This fact also
motivates us to perform an extensive comparison of the gradient-based
optimization schemes on a wide variety of problems from different
applications. These consist of the symplectic target problem which arises in
optimal control of symplectic quantum gate, the computation of symplectic
eigenvalues of a symmetric positive-(semi)definite matrix which can be used
for stability analysis of gyroscopic systems, and model reduction of
Hamiltonian systems.

Especially, we would like to emphasize the application to structure-preserving
model reduction of Hamiltonian systems. Such systems have a specific structure
and they possess some underlying physical properties such as conservation of
energy, described by the Hamiltonian function, and conservation of mass which
should be preserved during the model reduction process. This can be achieved
by using the proper symplectic decomposition (PSD) model reduction approach
developed in [15]. It consists in determining a symplectic reduced basis
matrix from a set of snapshots, i.e., solutions at different time instances,
which minimizes the symplectic projection error in a least squares sense. This
leads to a nonlinear optimization problem of the form (1). Unlike the
orthogonal case, where a singular value decomposition (SVD) of the data matrix
can be used, an explicit solution to the problem with the simplecticity
constraint is still unknown. To simplify this problem, in [15], the optimal
candidates are limited to the ones that additionally have orthonormal columns.
The same condition is also used in a greedy approach for parametric
Hamiltonian systems in [16]. Alternatively, the additional orthonormality
condition is eliminated in [17], but the approach there is based on the _SVD-
like decomposition_ and therefore, as those in [15], its solution is in
general not guaranteed to be optimal although it can have this property for a
special class of Hamiltonian systems in a restrictive setting [18].

In contrast to these approaches, in this paper, thanks to the progress on
optimization techniques on the symplectic Stiefel manifold, we address the
symplectic model reduction problem in an optimal way. A similar approach was
recently considered in [19] in the framework of the symplectic Grassmann
manifold [2] but no model reduction errors were reported there. For nonlinear
Hamiltonian systems, the symplectic projection is combined with an
approximation of the nonlinear term computed by the discrete empirical
interpolation decomposition method (DEIM) proposed in [20] and its structure-
preserving variant [21]. In our experiments, we numerically evaluate the
optimization-based model reduction methods and compare them with other
structure-preserving reduction techniques developed for Hamiltonian systems.

After presenting the general notation, we organize the paper as follows. In
Section 2, we briefly review basic geometric concepts and facts for the
symplectic Riemannian manifold ${\mathrm{Sp}(2k,2n)}$. Notably, we present in
detail the canonical-like and Euclidean metrics, and the corresponding
formulations of the Riemannian gradient of the cost function. Section 3 is
devoted to retractions, which are indispensable in Riemannian optimization. We
start with reviewing the Cayley and quasi-geodesic retractions. Then, in the
main part of this section, we introduce a new retraction based on an SR
decomposition and discuss the existence conditions and computational issues.
Section 4 recalls the non-monotone line search algorithm based on the
Riemannian gradient for solving the minimization problem (1). In Section 5, we
present three applications with several test models and report on numerical
results illustrating the properties of different Riemannian optimization
schemes. Finally, the concluding remarks are given in Section 6.

Notation. We denote by ${\cal S}_{\mathrm{sym}}(m)$ and ${\cal
S}_{\mathrm{skew}}(m)$ the sets of all $m\times m$ real symmetric and skew-
symmetric matrices, respectively, and
$\operatorname*{skew}(A)=\tfrac{1}{2}(A-A^{T})$ stands for the skew-symmetric
part of a square matrix $A$. The determinant, the trace, and the image of a
matrix $A$ are denoted by $\det(A)$, $\operatorname*{tr}(A)$, and
$\operatorname*{im}(A)$, respectively. If $A$ is a complex matrix, $(A)$ and
$\operatorname*{Im}(A)$ denote the real and imaginary part of $A$,
respectively. For $A_{1},\ldots,A_{\ell}$ being square matrices, we denote by
$\operatorname*{diag}(A_{1},\ldots,A_{\ell})$ the block diagonal matrix. The
Frobenius and spectral matrix norms are denoted by $\|\cdot\|_{F}$ and
$\|\cdot\|_{2}$, respectively, and the Euclidean vector norm is denoted by
$\|\cdot\|$. For a matrix $M\in{\mathbb{R}^{2k\times 2k}}$ and $j=1,\ldots,k$,
we denote by $(M)_{1:2j,1:2j}$ the $2j\times 2j$ leading principal submatrix
of $M$. The dimension of a subspace $\,\mathcal{U}\subset\mathbb{R}^{n\times
m}$ is denoted by $\dim(\mathcal{U})$. Finally, for a function $h$ defined on
the Euclidean space ${\mathbb{R}^{2n\times 2k}}$, $\nabla h$ denotes the
standard Euclidean gradient of $h$.

##  2 Riemannian geometry of the symplectic Stiefel manifold

In this section, we briefly review the geometric structure of the symplectic
Stiefel manifold ${\mathrm{Sp}(2k,2n)}$ studied recently in [1, 4]. Let us
start with the result confirming that ${\mathrm{Sp}(2k,2n)}$ is a smooth
embedded submanifold of the Euclidean space $\mathbb{R}^{2n\times 2k}$ and it
has dimension $4nk-k(2k-1)$, see [1, Proposition 3.1]. Alternatively, one can
show that ${\mathrm{Sp}(2k,2n)}$ is diffeomorphic to a quotient space
$\mathrm{Sp}(2n)/\mathrm{Sp}(2n-2k)$ and, hence, it admits a structure of a
quotient manifold [2, Proposition 3.1]. Further, it has been shown in [1,
Proposition 3.3] that the _tangent space_ of ${\mathrm{Sp}(2k,2n)}$ at
$X\in{\mathrm{Sp}(2k,2n)}$ is given by

$\displaystyle{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}$
$\displaystyle=\\!\left\\{Z\in{\mathbb{R}^{2n\times
2k}}\;:\;Z^{T}J_{2n}X+X^{T}J_{2n}Z=0\right\\}$ (2)
$\displaystyle=\\!\left\\{XJ_{2k}W\\!+\\!J_{2n}X_{\perp}K\;:\;W\\!\in\\!\mathcal{S}_{\rm
sym}(2k),K\\!\in\\!\mathbb{R}^{(2n-2k)\times 2k}\right\\}\\!,$ (3)

where $X_{\perp}\in\mathbb{R}^{2n\times(2n-2k)}$ has full rank and satisfy
$X^{T}X_{\perp}=0$. The choice for $X_{\perp}$ is obviously not unique and has
a certain effect on numerical performance. It has been shown numerically that
$X_{\perp}$ with orthonormal columns, i.e., $X_{\perp}^{T}X_{\perp}=I$, is
preferred, see [1, Section 6] for details. In what follows, we will restrict
ourselves to this choice.

The Riemannian structure of the symplectic Stiefel manifold
${\mathrm{Sp}(2k,2n)}$ greatly depends on the metric. Let
$g_{X}:{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}\times{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}\to\mathbb{R}$
denote a Riemannian metric on ${\mathrm{Sp}(2k,2n)}$ at $X$. Note that
sometimes, the employed metrics does not depend on $X$. To simplify the
notation, we omit the subscript. The _normal space_ to ${\mathrm{Sp}(2k,2n)}$
at $X\in{\mathrm{Sp}(2k,2n)}$ with respect to $g$ is defined as

$\bigl{(}{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}\bigr{)}^{\perp}=\bigl{\\{}N\in\mathbb{R}^{2n\times
2k}\;:\;g(N,Z)=0\text{ for all
}Z\in{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}\bigr{\\}}.$

It is well known that any $Y\in{\mathbb{R}^{2n\times 2k}}$ can be decomposed
as

$Y=\mathcal{P}_{X}(Y)+\mathcal{P}_{X}^{\perp}(Y),$

where $\mathcal{P}_{X}$ and $\mathcal{P}_{X}^{\perp}$ denote the orthogonal
projections with respect to $g$ onto the tangent and normal spaces,
respectively.

The _Riemannian gradient_ of a differentiable function
$f:{\mathrm{Sp}(2k,2n)}\to\mathbb{R}$ at $X\in{\mathrm{Sp}(2k,2n)}$ with
respect to the metric $g$, denoted by $\operatorname*{grad}f(X)$, is defined
as the unique element of ${\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}$ that satisfies
the condition

$g\bigl{(}\operatorname*{grad}f(X),Z\bigr{)}=\mathrm{D}\bar{f}(X)[Z]\qquad\text{
for all }Z\in{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)},$

where $\bar{f}$ is a smooth extension of $f$ around $X$ in
$\mathbb{R}^{2n\times 2k}$, and $\mathrm{D}\bar{f}(X)$ denotes the Fréchet
derivative of $\bar{f}$ at $X$. Using [22, (3.37)], the Riemannian gradient
can be determined as

$\operatorname*{grad}f(X)=\mathcal{P}_{X}\bigl{(}\operatorname*{grad}\bar{f}(X)\bigr{)}.$
(4)

Next, we introduce the canonical-like and Euclidean metrics considered in [1,
4] and present particular geometric concepts related to these metrics.

###  2.1 Canonical-like metric

A class of metrics on ${\mathrm{Sp}(2k,2n)}$ has been introduced in [1]. For a
parameter $\rho>0$ and the tangent vectors
$Z_{i}=XJ_{2k}W_{i}+J_{2n}X_{\perp}K_{i}$ with $W_{i}\in{\cal
S}_{\mathrm{sym}}(2k)$ and $K_{i}\in\mathbb{R}^{(2n-2k)\times 2k}$ for
$i=1,2$, the _canonical-like metric_ is defined as

$\displaystyle g_{c,\rho}(Z_{1},Z_{2})$
$\displaystyle:=\frac{1}{\rho}\operatorname*{tr}(W_{1}^{T}W_{2})+\operatorname*{tr}(K_{1}^{T}K_{2}).$

The normal space with respect to this metric is then given by

$\bigl{(}{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}\bigr{)}_{c}^{\perp}=\bigl{\\{}XJ_{2k}\varOmega\enskip:\enskip\varOmega\in{\cal
S}_{\mathrm{skew}}(2k)\bigr{\\}}.$

Further, the corresponding orthogonal projections $\mathcal{P}_{X,c}$ and
$\mathcal{P}_{X,c}^{\perp}$ onto the tangent and normal spaces have the
following form

$\mathcal{P}_{X,c}(Y)=S_{X,Y}J_{2n}X,\qquad\mathcal{P}_{X,c}^{\perp}(Y)=XJ_{2k}\operatorname*{skew}(X^{T}J_{2n}^{T}Y),$

where $Y\in{\mathbb{R}^{2n\times 2k}}$, and

$S_{X,Y}=G_{X}Y(XJ_{2k})^{T}+XJ_{2k}(G_{X}Y)^{T},\quad
G_{X}=I_{2n}-\frac{1}{2}XJ_{2k}X^{T}J_{2n}^{T},$

see [1, Proposition 4.3]. Using (4), the Riemannian gradient of a function $f$
with respect to the metric $g_{c,\rho}$ can then be represented as

$\mathrm{grad}_{c}f(X)=\mathcal{P}_{X,c}(\mathrm{grad}_{c}\bar{f}(X))=S_{{X,\nabla\bar{f}}}J_{2n}X$

with

${S_{X,\nabla\bar{f}}}=H_{X}\nabla\bar{f}(X)(XJ_{2k})^{T}+XJ_{2k}(H_{X}\nabla\bar{f}(X))^{T},\quad
H_{X}=\frac{\rho}{2}XX^{T}+J_{2n}X_{\perp}X_{\perp}^{T}J_{2n}^{T},$

see [1, Proposition 4.5] for details.

###  2.2 Euclidean metric

Another metric on ${\mathrm{Sp}(2k,2n)}$ has been investigated in [4]. For
$Z_{i}=XJ_{2k}W_{i}+J_{2n}X_{\perp}K_{i}\in{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}$
with the matrices $W_{i}\in{\cal S}_{\mathrm{sym}}(2k)$ and
$K_{i}\in\mathbb{R}^{(2n-2k)\times 2k}$ for $i=1,2$, the _Euclidean metric_ is
defined as

$\displaystyle g_{e}(Z_{1},Z_{2})$
$\displaystyle:=\operatorname*{tr}(Z_{1}^{T}Z_{2})$
$\displaystyle~{}=\operatorname*{tr}(W_{1}^{T}J_{2k}^{T}X^{T}XJ_{2k}W_{2})+\operatorname*{tr}(K_{1}^{T}X^{T}_{\perp}X_{\perp}K_{2})$
$\displaystyle\quad\enskip+\operatorname*{tr}(W_{1}^{T}J_{2k}^{T}X^{T}J_{2n}{}X_{\perp}K_{2})+\operatorname*{tr}(K_{1}^{T}X^{T}_{\perp}J_{2n}^{T}XJ_{2k}W_{2}).$

Based on this metric, the normal space to ${\mathrm{Sp}(2k,2n)}$ can be
represented as

$\bigl{(}{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}\bigr{)}_{e}^{\perp}=\bigl{\\{}J_{2n}X\varOmega\enskip:\enskip\varOmega\in{\cal
S}_{\mathrm{skew}}(2k)\bigr{\\}},$ (5)

see [4, Proposition 1]. Further, using (3) and (5), we obtain the following
expressions for the orthogonal projections $\mathcal{P}_{X,e}$ and
$\mathcal{P}_{X,e}^{\perp}$ onto the tangent and normal spaces, respectively,
with respect to the Euclidean metric [4, Proposition 2]:

$\mathcal{P}_{X,e}(Y)=Y-J_{2n}X{\varOmega_{X,Y}},\qquad\mathcal{P}_{X,e}^{\perp}(Y)=J_{2n}X{\varOmega_{X,Y}},$

where $\varOmega_{X,Y}\in{\cal S}_{\mathrm{skew}}(2k)$ is the solution of the
Lyapunov equation

$X^{T}X\varOmega+\varOmega\,X^{T}X=2\operatorname*{skew}(X^{T}J_{2n}^{T}Y).$
(6)

The existence and uniqueness of the solution of the Lyapunov equation (6)
immediately follows from the fact that the coefficient matrix $X^{T}X$ is
symmetric and positive definite [23].

Since ${\mathrm{Sp}(2k,2n)}$ is endowed with the Euclidean metric $g_{e}$, the
corresponding Riemannian gradient $\mathrm{grad}_{e}f(X)$ can directly be
calculated by using (4) as follows

$\mathrm{grad}_{e}f(X)=\mathcal{P}_{X,e}(\mathrm{grad}_{e}\bar{f}(X))=\mathcal{P}_{X,e}(\nabla\bar{f}(X))=\nabla\bar{f}(X)-J_{2n}X\varOmega_{X,\nabla\bar{f}},$

where $\varOmega_{X,\nabla\bar{f}}\in{\cal S}_{\mathrm{skew}}(2k)$ solves the
Lyapunov equation

$X^{T}X\,\varOmega+\varOmega\,X^{T}X=2\operatorname*{skew}\left(X^{T}J_{2n}^{T}\nabla\bar{f}(X)\right).$

We summarize geometrical notions and their formulations for different metrics
on the symplectic Stiefel manifold ${\mathrm{Sp}(2k,2n)}$ in Table 2.2.

Table 1: Geometric objects on the symplectic Stiefel manifold
${\mathrm{Sp}(2k,2n)}$ for different metrics. Here,
$Z_{i}=XJ_{2k}W_{i}+J_{2n}X_{\perp}K_{i}$ with $W_{i}\in{\cal
S}_{\mathrm{sym}}(2k)$ and $K_{i}\in\mathbb{R}^{(2n-2k)\times 2k}$, and the
definitions of $S_{X,Y}$, $S_{X,\nabla\bar{f}}$, $\varOmega_{X,Y}$, and
$\varOmega_{X,\nabla\bar{f}}$ can be found in Section 2.

| Canonical-like | Euclidean
---|---|---
6pt.6pt. metric | $g(Z_{1},Z_{2})$ | $\frac{1}{\rho}\operatorname*{tr}(W_{1}^{T}W_{2})+\operatorname*{tr}(K_{1}^{T}K_{2})$ | $\operatorname*{tr}(Z_{1}^{T}Z_{2})$
normal space | $\bigl{(}{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}\bigr{)}^{\perp}$ | ${XJ_{2k}\varOmega:\varOmega\in{\cal S}_{\mathrm{skew}}(2k)}$ | ${J_{2n}X\varOmega:\varOmega\in{\cal S}_{\mathrm{skew}}(2k)}$
projection | $\mathcal{P}_{X}(Y)$ | $S_{X,Y}J_{2n}X$ | $Y-J_{2n}X{\varOmega_{X,Y}}$
1pt. | $\mathcal{P}_{X}^{\perp}(Y)$ | $XJ_{2k}\operatorname*{skew}(X^{T}J_{2n}^{T}Y)$ | $J_{2n}X{\varOmega_{X,Y}}$
gradient | $\operatorname*{grad}f(X)$ | $S_{X,\nabla\bar{f}}\,J_{2n}X$ | $\nabla\bar{f}(X)-J_{2n}X\varOmega_{X,\nabla\bar{f}}$

##  3 Retractions on ${\mathrm{Sp}(2k,2n)}$

Retractions play a key role in Riemannian optimization since they allow to
transfer data from a tangent space, which contains gradients of smooth
functions as its elements, to the manifold, which is the search space in the
optimization problem. The concept of retraction on the symplectic Stiefel
manifold ${\mathrm{Sp}(2k,2n)}$ can be defined as follows. Let

$\mathrm{T}{\mathrm{Sp}(2k,2n)}=\bigcup_{X\in{\mathrm{Sp}(2k,2n)}}{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}$

be the tangent bundle to ${\mathrm{Sp}(2k,2n)}$. A smooth mapping
$\mathcal{R}:\mathrm{T}{\mathrm{Sp}(2k,2n)}\to{\mathrm{Sp}(2k,2n)}$ is called
a  _retraction_ if for all $X\in{\mathrm{Sp}(2k,2n)}$, the restriction of
$\mathcal{R}$ to ${\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}$, denoted by
$\mathcal{R}_{X}$, satisfies the following properties:

  1. 1)

$\mathcal{R}_{X}(0_{X})=X$, where $0_{X}$ denotes the origin of
${\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}$;

  2. 2)

$\left.\tfrac{{\rm d}}{{\rm d}t}\mathcal{R}_{X}(tZ)\right|_{t=0}=Z$ for all
$Z\in{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}$.

The retraction provides a first-order approximation to the Riemannian
exponential map [24, Section 5.4] which is of crucial importance in
optimization algorithms on Riemannian manifolds. Since retractions are
computationally less expensive compared to the exponential map while retaining
the convergence properties of optimization schemes, they have attracted
extensive interest in research. There are different approaches for the
construction of retractions. In Subsections 3.1 and 3.2, we briefly review the
Cayley and quasi-geodesic retractions on ${\mathrm{Sp}(2k,2n)}$ first
introduced in [1] and then, in Subsection 3.3, we propose a new one based on
an SR decomposition.

###  3.1 Cayley retraction

For $Z\in{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}$, the _Cayley retraction_ on the
symplectic Stiefel manifold ${\mathrm{Sp}(2k,2n)}$ is defined by

$\mathcal{R}_{X}^{\rm
cay}(Z):=\Bigl{(}I_{2n}-\frac{1}{2}S_{X,Z}J_{2n}\Bigr{)}^{-1}\Bigl{(}I_{2n}+\frac{1}{2}S_{X,Z}J_{2n}\Bigr{)}X,$
(7)

where $S_{X,Z}=G_{X}Z(XJ_{2k})^{T}+XJ_{2k}(G_{X}Z)^{T}$ and
$G_{X}=I_{2n}-\frac{1}{2}XJ_{2k}X^{T}\\!J_{2n}^{T}$, see [1, Section 5.2]. It
exists if and only if $2$ is not an eigenvalue of the Hamiltonian matrix
$S_{X,Z}J_{2n}$. The computation of the Cayley retraction (7) involves a
matrix inverse of size $2n\times 2n$. In fact, it can be economically computed
by the Sherman–Morrison–Woodbury formula as discussed in [1, Proposition 5.5]
which requires the inversion of a $4k\times 4k$ matrix. This is certainly
advantageous since in most applications, $k$ is considerably smaller than $n$.
In [2, Proposition 5.2], an equivalent expression for the Cayley retraction

$\mathcal{R}_{X}^{\rm
cay}(Z)=-X+\bigl{(}H_{X,Z}+2X\bigr{)}\Bigl{(}\frac{1}{4}J_{2k}^{T}H_{X,Z}^{T}J_{2n}H_{X,Z}-\frac{1}{2}J_{2k}^{T}X^{T}J_{2n}Z+I_{2k}\Bigr{)}^{-1}$

with $H_{X,Z}=Z-XJ_{2k}^{T}X^{T}J_{2n}Z$ has been presented, which requires
solving a linear system with a $2k\times 2k$ matrix only. Note that this
expression coincides with that considered in [25, Lemma 3.1].

###  3.2 Quasi-geodesic retraction

Given $Z\in{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}$, the _quasi-geodesic
retraction_ is defined by

$\mathcal{R}^{\operatorname*{qgeo}}_{X}(Z):=\left[X,\;Z\right]\exp\left(\begin{bmatrix}-J_{2k}W&J_{2k}Z^{T}J_{2n}Z\\\
I_{2k}&-J_{2k}W\end{bmatrix}\right)\begin{bmatrix}I_{2k}\\\
0\end{bmatrix}\exp(J_{2k}W),$ (8)

where $W=X^{T}J_{2n}Z$ and $\exp(\cdot)$ denotes the matrix exponential [1,
Section 5.1]. This retraction is globally defined. Note that the calculation
of (8) requires computing two exponentials of matrices of size $4k\times 4k$
and $2k\times 2k$ which is dominating for relatively large $k$.

###  3.3 SR decomposition-based retraction

For some matrix manifolds, retractions can also be defined using related
matrix decompositions, e.g., [24]. In the case of the symplectic Stiefel
manifold ${\mathrm{Sp}(2k,2n)}$, the SR decomposition introduced first in [5]
appears to be very useful to define a new decomposition-based retraction.

To begin with, we consider the perfect shuffle permutation matrix

$P_{2k}=[e_{1},e_{3},\ldots,e_{2k-1},e_{2},\ldots,e_{2k}],$ (9)

where $e_{j}$, $j=1,\ldots,2k$, is the $j$-th canonical basis vector of
$\mathbb{R}^{2k}$. It can be shown by direct calculation that this matrix is
orthogonal and

$P_{2k}J_{2k}P_{2k}^{T}=\operatorname*{diag}(J_{2},\ldots,J_{2})=:\hat{J}_{2k},$
(10)

where $P_{2k}^{T}=[e_{1},e_{k+1},e_{2},e_{k+2},\ldots,e_{k},e_{2k}]$. Further,
we introduce a congruence matrix set

$T_{2k}(P_{2k})=\\{P_{2k}^{T}\hat{R}P_{2k}\enskip:\enskip\hat{R}\in\mathbb{R}^{2k\times
2k}\;\text{ is upper triangular}\\}.$

Then an  _SR decomposition_ of $A\in\mathbb{R}^{2n\times 2k}$ with $k\leq n$
is defined as

$A=SR,$ (11)

where $S\in{\mathrm{Sp}(2k,2n)}$ and $R\in T_{2k}(P_{2k})$. The existence of
such a decomposition for square matrices has been established in [6, Theorem
3.8]. This result can be adapted to rectangular matrices as follows.

######  Theorem 3.1.

Let $A\in\mathbb{R}^{2n\times 2k}$ have full column rank and let $P_{2k}$ be
as in (9). There exists an SR decomposition $A=SR$ with
$S\in{\mathrm{Sp}(2k,2n)}$ and $R\in T_{2k}(P_{2k})$ if and only if all
leading minors of even dimension of the matrix $P_{2k}A^{T}J_{2n}AP_{2k}^{T}$
are nonzero, i.e., $\mathrm{det}(P_{2k}A^{T}J_{2n}AP^{T}_{2k})_{1:2j,1:2j}\neq
0$ for $j=1,\dots,k$.

###### Proof.

Assume that $A$ has an SR decomposition $A=SR$ with $S\in{\mathrm{Sp}(2k,2n)}$
and $R\in T_{2k}(P_{2k})$. Then we obtain

$\displaystyle P_{2k}A^{T}J_{2n}AP_{2k}^{T}$
$\displaystyle=P_{2k}R^{T}S^{T}J_{2n}SRP_{2k}^{T}=P_{2k}R^{T}J_{2k}RP_{2k}^{T}=\hat{R}^{T}\hat{J}_{2k}\hat{R},$

where $\hat{R}=P_{2k}RP_{2k}^{T}$ is upper triangular and $\hat{J}_{2k}$ is as
in (10). Since $A$ has full column rank, $R$ and therefore $\hat{R}$ are
nonsingular. In view of these facts, for $j=1,\ldots,k$, the $2j\times 2j$
leading principal submatrices of $\hat{R}^{T}\hat{J}_{2k}\hat{R}$ have the
form

$(\hat{R}^{T}\hat{J}_{2k}\hat{R})_{1:2j,1:2j}=(\hat{R}^{T})_{1:2j,1:2j}(\hat{J}_{2k})_{1:2j,1:2j}(\hat{R})_{1:2j,1:2j}$

and, hence, they are, as a product of nonsingular matrices, nonsingular. Thus,
all leading minors of even dimension of $P_{2k}A^{T}J_{2n}AP_{2k}^{T}$ are
nonzero.

On the other hand, if all leading minors of even dimension of the skew-
symmetric matrix $P_{2k}A^{T}J_{2n}AP_{2k}^{T}$ are nonzero, then by [26,
Theorem 2.2] this matrix has the Cholesky-like decomposition

$P_{2k}A^{T}J_{2n}AP_{2k}^{T}=\hat{R}^{T}\hat{J}_{2k}\hat{R}$ (12)

with a nonsingular upper triangular matrix $\hat{R}$. In this case, we have

$A^{T}J_{2n}A=P_{2k}^{T}\hat{R}^{T}\hat{J}_{2k}\hat{R}P_{2k}=R^{T}J_{2k}R,$

where $R=P_{2k}^{T}\hat{R}P_{2k}\in T_{2k}(P_{2k})$ is nonsingular.
Furthermore, the matrix $S=AR^{-1}$ is symplectic due to
$S^{T}J_{2n}S=R^{-T}A^{T}J_{2n}AR^{-1}=J_{2k}$. Thus, $A=SR$ is an SR
decomposition. ∎

Similarly to the square case [6, Remark 3.9], the SR decomposition of
rectangular matrices is non-unique. Indeed, if $A$ is decomposed as in (11),
then for any nonsingular diagonal matrix $D\in\mathbb{R}^{k\times k}$,  A =
(Sdiag(D,D^-1))(diag(D^-1,D)R)  is also an SR decomposition of $A$. The
freedom of choice of the factors $S$ and $R$ in (11) is usually exploited to
improve numerical stability, see, e.g., [27, 28]. To serve the purpose of
constructing a retraction on ${\mathrm{Sp}(2k,2n)}$, we follow [29] and
restrict the factor $R$ in (11) to the matrix set

$\begin{array}[]{rl}T_{2k}^{0}(P_{2k})=\bigl{\\{}P_{2k}^{T}\hat{R}P_{2k}\enskip:&\hat{R}=[r_{ij}]\in\mathbb{R}^{2k\times
2k}\text{ is upper triangular with }r_{2j-1,2j}=0,\bigr{.}\\\
&\bigl{.}\,r_{2j-1,2j-1}>0,\text{ and }|r_{2j,2j}|=r_{2j-1,2j-1}\text{ for
}j=1,\ldots,k\bigr{\\}}.\end{array}$ (13)

This choice guarantees the uniqueness of the resulting SR decomposition [29].
Note that this fact can also be inferred from [26, Theorem 2.2], which
establishes the uniqueness of the Cholesky-like decomposition (12) with
$\hat{R}$ as in (13).

In order to compute the SR decomposition (11), we employ a symplectic
Gram–Schmidt algorithm developed in [7]. Based on Theorem 3.1, we investigate
the well-posedness of this algorithm which was not discussed in [7]. For ease
of explanation, our consideration is divided into three steps.

First, the SR decomposition of two-column matrices is needed. Given a matrix
$A=[a_{1},a_{2}]\in\mathbb{R}^{2n\times 2}$, we are looking for an elementary
SR (ESR) decomposition $A=SR$ with a $2\times 2$ upper triangular matrix

$R=\begin{bmatrix}r_{11}&r_{12}\\\ 0&r_{22}\end{bmatrix}$ (14)

and a two-column symplectic matrix $S=[s_{1},s_{2}]$. It is straightforward to
verify that $S=[s_{1},s_{2}]$ is symplectic if and only if
$s_{1}^{T}J_{2n}s_{2}=1$. Using this relation, we find that
$r_{11}r_{22}=a_{1}^{T}J_{2n}a_{2}$ and $r_{12}$ is arbitrary. By Theorem 3.1,
the existence of the ESR decomposition $A=SR$ is equivalent to the condition
$a_{1}^{T}J_{2n}a_{2}\not=0$ which is also known as the non-isotropy condition
for the subspace spanned by $a_{1}$ and $a_{2}$. Requiring $R\in
T_{2}^{0}(P_{2})$, we obtain

$r_{12}=0,\quad r_{11}=\sqrt{|a_{1}^{T}J_{2n}a_{2}|},\quad
r_{22}=\mbox{sign}(a_{1}^{T}J_{2n}a_{2})r_{11},$ (15)

where $\mbox{sign}(\cdot)$ denotes the sign of the corresponding value. The
resulting decomposition is referred to as the diagonal elementary SR (DESR)
decomposition which can easily be seen to be unique. Note that the choice (15)
corresponds to the version ESR4 in [27] which was proven to yield $R$ with a
minimal condition number. For convenience, we summarize the computation of the
DESR decomposition in Algorithm 1 and collect its properties in the following
lemma.

Algorithm 1 Diagonal elementary SR (DESR) decomposition

0: $A=[a_{1},a_{2}]\in\mathbb{R}^{2n\times 2}$.

0: $S=[s_{1},s_{2}]\in\mathrm{Sp}(2,2n)$ and
$R=\operatorname*{diag}(r_{11},r_{22})$ with $0<r_{11}=|r_{22}|$ such that
$A=SR$.

1: Compute $\omega=a_{1}^{T}J_{2n}a_{2}$.

2: if $\omega\neq 0$ then

3: Compute $\ r_{11}=\sqrt{|\omega|}$ and
$r_{22}=\mathrm{sign}(\omega)r_{11}$.

4: Compute $s_{1}=a_{1}/r_{11}$ and $s_{2}=a_{2}/r_{22}$.

5: else

6: Error: the DESR decomposition does not exist.

7: end if

######  Lemma 3.2.

For a matrix $A=[a_{1},a_{2}]\in\mathbb{R}^{2n\times 2}$ with
$a_{1}^{T}J_{2n}a_{2}\neq 0$, Algorithm 1 produces a unique DESR decomposition
$A=SR$. In particular, if $A$ is symplectic, then $S=A$ and $R=I_{2}$.

Second, for $A\in\mathbb{R}^{2n\times 2k}$, we consider a decomposition

$A=\hat{S}\hat{R},$ (16)

where $\hat{S}\in\mathbb{R}^{2n\times 2k}$ is a perfect shuffle permuted
symplectic (PSPS) matrix satisfying

$\hat{S}^{T}J_{2n}\hat{S}=\hat{J}_{2k},$ (17)

and

$\hat{R}=\begin{bmatrix}\hat{R}_{11}&\cdots&\hat{R}_{1k}\\\ &\ddots&\vdots\\\
&&\hat{R}_{kk}\end{bmatrix}$ (18)

is upper triangular with $\hat{R}_{ij}\in\mathbb{R}^{2\times 2}$ for $1\leq
i\leq j\leq k$. Additionally, the blocks $\hat{R}_{jj}$, $1\leq j\leq k$, are
assumed to be diagonal with diagonal elements ordered nonincreasingly and
having the same absolute value. Note that the matrix $R$ in (14) with entries
satisfying (15) has this structure. It follows from (10) and (17) that
$\hat{S}$ is PSPS if and only if $\hat{S}P_{2k}$ is symplectic, which
justifies the name PSPS. Moreover, if $\hat{S}$ is PSPS, the $2n\times 2$
blocks $\hat{S}_{j}$ of the matrix $\hat{S}=[\hat{S}_{1},\dots,\hat{S}_{k}]$
are symplectic. The following lemma establishes the existence and uniqueness
of the decomposition (16).

######  Lemma 3.3.

Let $A\in\mathbb{R}^{2n\times 2k}$ have full column rank. The decomposition
(16) with a PSPS matrix $\hat{S}$ and an upper triangular matrix $\hat{R}$ as
in (18) exists and is unique if and only if all leading minors of even
dimension of the matrix $A^{T}J_{2n}A$ are nonzero, i.e.,
$\mathrm{det}(A^{T}J_{2n}A)_{1:2j,1:2j}\neq 0$ for $j=1,\dots,k$.

###### Proof.

The necessity and sufficiency for the existence of the decomposition (16) can
be proved analogously to Theorem 3.1. The special choice of the diagonal
blocks of $\hat{R}$ guarantees the uniqueness. ∎

Partitioning $A=[A_{1},\ldots,A_{k}]$ with $A_{j}\in\mathbb{R}^{2n\times 2}$
makes it possible to use the DESR decompositions for computing the $2n\times
2$ symplectic blocks $\hat{S}_{j}$ of
$\hat{S}=[\hat{S}_{1},\ldots,\hat{S}_{k}]$ and the $2\times 2$ blocks
$\hat{R}_{jj}$ of $\hat{R}$ in (18). To this end, we consider

$\hat{J}_{2k}^{T}\hat{S}^{T}J_{2n}A=\hat{J}_{2k}^{T}\hat{S}^{T}J_{2n}\hat{S}\hat{R}=\hat{J}_{2k}^{T}\hat{J}_{2k}\hat{R}=\hat{R},$

where equalities follow from (16) and (17). According to the block partition
of $A$ and $\hat{S}$, the block diagonal structure of $\hat{J}_{2k}$, and the
block upper triangular structure of $\hat{R}$, we can work out a block
Gram–Schmidt-type algorithm for computing the decomposition (16), see
Algorithm 2.

Algorithm 2 Basic symplectic Gram–Schmidt algorithm

0: $A=[A_{1},\ldots,A_{k}]\in\mathbb{R}^{2n\times 2k}$ with
$A_{j}\in\mathbb{R}^{2n\times 2}$ for $j=1,\ldots,k$.

0: A PSPS matrix $\hat{S}=[\hat{S}_{1},\ldots,\hat{S}_{k}]$ and an upper
triangular matrix $\hat{R}$ as in (18) such that $A=\hat{S}\hat{R}$.

1: Compute the DESR decomposition $A_{1}=\hat{S}_{1}\hat{R}_{11}$ using
Algorithm 1.

2: for $j=2,\ldots,k$ do

3: for $i=1,\ldots,j-1$ do

4: Compute $\hat{R}_{ij}=J_{2}^{T}\hat{S}_{i}^{T}J_{2n}A_{j}$.

5: end for

6: Compute $W_{j}=A_{j}-\sum\limits_{i=1}^{j-1}\hat{S}_{i}\hat{R}_{ij}$.

7: Compute the DESR decomposition $W_{j}=\hat{S}_{j}\hat{R}_{jj}$ using
Algorithm 1.

8: end for

The following lemma provides the sufficient conditions for the existence of
the DESR decomposition of the matrices $A_{1}$ and $W_{j}$, $j=2,\ldots,k$, in
Algorithm 2.

######  Lemma 3.4.

Let $A=[A_{1},\ldots,A_{k}]$ with $A_{j}\in\mathbb{R}^{2n\times 2}$ for
$j=1,\ldots,k$ be such that all leading minors of even dimension of the matrix
$A^{T}J_{2n}A$ are nonzero. Then Algorithm 2 produces the decomposition (16)
with a PSPS matrix $\hat{S}$ and an upper triangular matrix $\hat{R}$ as in
(18) without breakdown. In particular, if $A$ is PSPS, then $\hat{S}=A$ and
$\hat{R}=I_{2k}$.

###### Proof.

We show by finite induction that the matrices $W_{1}=A_{1}$,
$W_{2},\ldots,W_{k}$ in Algorithm 2 have the DESR decomposition. By
assumption, the matrix $W_{1}^{T}J_{2n}W_{1}=A_{1}^{T}J_{2n}A_{1}$ is
nonsingular, and hence by Lemma 3.2, the matrix $W_{1}$ has the DESR
decomposition. Assume that $j-1$ steps in Algorithm 2 have been executed. For
simplicity reasons, we introduce the following matrices
$B_{2(j-1)}=[A_{1},\ldots,A_{j-1}]$,
$Q_{2(j-1)}=[\hat{S}_{1},\ldots,\hat{S}_{j-1}]$, and

$\displaystyle\hat{R}_{2(j-1)}$
$\displaystyle=\begin{bmatrix}\hat{R}_{11}&\cdots&\hat{R}_{1,j-1}\\\
&\ddots&\vdots\\\ &&\hat{R}_{j-1,j-1}\end{bmatrix}.$

Note that by construction, $\hat{R}_{2(j-1)}$ is nonsingular. Using the
expressions for $\hat{R}_{ij}$ and $W_{j}$ in steps 4 and 6 of Algorithm 2,
respectively, the matrix $W_{j}$ can be represented as
$W_{j}=A_{j}-Q_{2(j-1)}\hat{J}_{2(j-1)}^{T}Q_{2(j-1)}^{T}J_{2n}A_{j}$. Then we
have

$\displaystyle W_{j}^{T}J_{2n}W_{j}$
$\displaystyle=A_{j}^{T}J_{2n}A_{j}-2A_{j}^{T}J_{2n}^{T}Q_{2(j-1)}\hat{J}_{2(j-1)}Q_{2(j-1)}^{T}J_{2n}A_{j}$
$\displaystyle\quad+A_{j}^{T}J_{2n}^{T}Q_{2(j-1)}\hat{J}_{2(j-1)}Q_{2(j-1)}^{T}J_{2n}Q_{2(j-1)}\hat{J}_{2(j-1)}^{T}Q_{2(j-1)}^{T}J_{2n}A_{j}$
$\displaystyle=A_{j}^{T}J_{2n}A_{j}-A_{j}^{T}J_{2n}^{T}Q_{2(j-1)}\hat{J}_{2(j-1)}Q_{2(j-1)}^{T}J_{2n}A_{j}.$

It follows from $B_{2(j-1)}=Q_{2(j-1)}\hat{R}_{2(j-1)}$ that
Q_2(j-1)^TJ_2nA_j=^R_2(j-1)^-TB_2(j-1)^TJ_2nA_j.  Therefore,

$\displaystyle W_{j}^{T}J_{2n}W_{j}$
$\displaystyle=A_{j}^{T}J_{2n}A_{j}\\!-\\!A_{j}^{T}J_{2n}^{T}B_{2(j-1)}\hat{R}_{2(j-1)}^{-1}\hat{J}_{2(j-1)}\hat{R}_{2(j-1)}^{-T}B_{2(j-1)}^{T}J_{2n}A_{j}$
$\displaystyle=A_{j}^{T}J_{2n}A_{j}\\!+\\!\bigl{(}\\!A_{j}^{T}J_{2n}^{T}B_{2(j-1)}\bigr{)}\\!\bigl{(}\\!\hat{R}_{2(j-1)}^{T}\hat{J}_{2(j-1)}\hat{R}_{2(j-1)}\bigr{)}^{-1}\\!\bigl{(}\\!B_{2(j-1)}^{T}J_{2n}A_{j}\bigr{)}$
$\displaystyle=A_{j}^{T}J_{2n}A_{j}\\!+\\!\bigl{(}\\!A_{j}^{T}J_{2n}^{T}B_{2(j-1)}\bigr{)}\\!\bigl{(}B_{2(j-1)}^{T}J_{2n}B_{2(j-1)}\bigr{)}^{-1}\\!\bigl{(}\\!B_{2(j-1)}^{T}J_{2n}A_{j}\bigr{)}$

is the Schur complement of the block $B_{2(j-1)}^{T}J_{2n}B_{2(j-1)}$ of the
matrix

$\begin{bmatrix}B_{2(j-1)}^{T}J_{2n}B_{2(j-1)}&B_{2(j-1)}^{T}J_{2n}A_{j}\\\
-A_{j}^{T}J_{2n}^{T}B_{2(j-1)}^{T}&A_{j}^{T}J_{2n}A_{j}\end{bmatrix}=[A_{1},\ldots,A_{j}]^{T}J_{2n}[A_{1},\ldots,A_{j}].$

Since this matrix is nonsingular, the Schur complement $W_{j}^{T}J_{2n}W_{j}$
is also nonsingular. Thus, by Lemma 3.2, $W_{j}$ has the DESR decomposition.
The particular case for a PSPS matrix $A$ follows from the uniqueness of the
decomposition (16). ∎

Finally, once the basic symplectic Gram–Schmidt algorithm with the DESR
decompositions described above is well-defined, the computation of the SR
decomposition (11) requires just two extra permutation steps as presented in
Algorithm 3.

Algorithm 3 Symplectic Gram–Schmidt algorithm

0: $A=[A_{1},\ldots,A_{k}]\in\mathbb{R}^{2n\times 2k}$ with
$A_{j}\in\mathbb{R}^{2n\times 2}$ for $j=1,\ldots,k$.

0: $S\in{\mathrm{Sp}(2k,2n)}$ and $R\in T_{2k}^{0}(P_{2k})$ such that $A=SR$.

1: Compute the decomposition $AP^{T}_{2k}=\hat{S}\hat{R}$ by using Algorithm
2.

2: Compute $S=\hat{S}P_{2k}$.

3: Compute $R=P^{T}_{2k}\hat{R}P_{2k}$.

Using Lemma 3.4, we can establish the following properties of Algorithm 3.

######  Proposition 3.5.

Let $A\in\mathbb{R}^{2n\times 2k}$ and let $P_{2k}$ be as in (9). Assume that
all leading minors of even dimension of $P_{2k}A^{T}J_{2n}AP_{2k}^{T}$ are
nonzero. Then Algorithm 3 produces a unique SR decomposition (11) with
$S\in{\mathrm{Sp}(2k,2n)}$ and $R\in T_{2k}^{0}(P_{2k})$.

###### Proof.

The existence and uniqueness of the SR decomposition (11) provided by
Algorithm 3 follow from that of the decomposition $AP^{T}_{2k}=\hat{S}\hat{R}$
with a PSPS matrix $\hat{S}$ and an upper triangular matrix $\hat{R}$ as in
(18) which is guaranteed by Lemma 3.4. Therefore, the decomposition
$AP^{T}_{2k}=\hat{S}\hat{R}=\hat{S}P_{2k}P^{T}_{2k}\hat{R}$ yields that
$A=(\hat{S}P_{2k})(P^{T}_{2k}\hat{R}P_{2k})=SR$, where
$S=\hat{S}P_{2k}\in{\mathrm{Sp}(2k,2n)}$ and $R=P^{T}_{2k}\hat{R}P_{2k}\in
T_{2k}^{0}(P_{2k})$. ∎

######  Remark 3.6.

For clarity of theoretical discussion, we have presented here the basic
symplectic Gram–Schmidt algorithm only. In practice, however, similarly to the
Gram–Schmidt orthonormalization process [23], for a better numerical behavior,
a modified basic symplectic Gram–Schmidt procedure should be used. This
procedure leads to a modified symplectic Gram–Schmidt algorithm, see [7] for
more details, which is indeed employed in numerical experiments reported in
Section 5.

We are now ready to introduce a new retraction on ${\mathrm{Sp}(2k,2n)}$ which
is based on the SR decomposition.

######  Theorem 3.7.

Given $X\in{\mathrm{Sp}(2k,2n)}$ and a tangent vector
$Z\in{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}$, denote by $\mathrm{sf}(X+Z)$ the
factor $S\in{\mathrm{Sp}(2k,2n)}$ in the SR decomposition $X+Z=SR$ computed by
Algorithm 3. Then the mapping

$\mathcal{R}_{X}^{\rm SR}(Z)=\mathrm{sf}(X+Z)$ (19)

defines a retraction on ${\mathrm{Sp}(2k,2n)}$.

###### Proof.

We consider the mapping

$\begin{array}[]{rcl}\Theta:{\mathrm{Sp}(2k,2n)}\times
T_{2k}^{0}(P_{2k})&\rightarrow&\mathcal{U}\\\ (S,R)&\mapsto&SR,\end{array}$

where $\mathcal{U}\subset\mathbb{R}_{*}^{2n\times 2k}$ consists of matrices
having an SR decomposition and $\mathbb{R}_{*}^{2n\times 2k}$ denotes the set
of real $2n\times 2k$ matrices of full rank. It follows from Theorem 3.1 that
$\,\mathcal{U}$ is an open subset of $\mathbb{R}_{*}^{2n\times 2k}$. Moreover,
$\Theta$ admits a neutral element $I_{2k}\in T_{2k}^{0}(P_{2k})$ satisfying

$\Theta(S,I_{2k})=SI_{2k}=S\quad\mbox{ for all }S\in{\mathrm{Sp}(2k,2n)}.$

Since $\mathrm{dim}\bigl{(}{\mathrm{Sp}(2k,2n)}\bigr{)}=4nk-k(2k-1)$ and
$\mathrm{dim}\bigl{(}T_{2k}^{0}(P_{2k})\bigr{)}=k(2k-1)$, we obtain that
dim(Sp(2k,2n))+dim(T_2k^0(P_2k))=4nk=dim(R_*^2n×2k).  The mapping $\Theta$ is
smooth as it is just the matrix multiplication restricted to the submanifolds.
Furthermore, the inverse mapping $\Theta^{-1}$ is defined on the whole
$\,\mathcal{U}$. For any $A\in\mathcal{U}$, the first component $S$ of
$\Theta^{-1}(A)$ is obtained by Algorithm 3 applied to $A$. Since this
algorithm consists of basic mathematical operations only, it is smooth on
$\mathcal{U}$. The second component of $\Theta^{-1}(A)$ is determined as
$R=S^{+}A$, where $S^{+}=J_{2k}^{T}S^{T}J_{2n}$ is the symplectic inverse of
$S$. This means that $\Theta$ is indeed a diffeomorphism. Then by [24,
Proposition 4.1.2], the mapping $\mathcal{R}_{X}^{\rm SR}$ in (19) is a
retraction on ${\mathrm{Sp}(2k,2n)}$. ∎

For the global convergence of a Riemannian gradient-based optimization
algorithm that employs this retraction, its domain is crucial. Given
$X\in{\mathrm{Sp}(2k,2n)}$, using Theorems 3.1 and 3.7, we can deduce that the
set of the tangent vectors $Z$, for which the retraction $\mathcal{R}_{X}^{\rm
SR}(Z)$ does not exist, has measure zero. This fact is unfortunately not
enough to guarantee the convergence of the associated algorithm as it requires
that the retraction is locally well-defined around the origin $0_{X}$ in the
tangent space, i.e., the domain of the retraction should contain an open ball
centered at $0_{X}$ in the tangent space
${\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}$, see [1, Theorem 5.7]. For the SR
retraction developed here, the following theorem indicates that the statement
holds for the ball of radius one.

######  Theorem 3.8.

Let $X\in{\mathrm{Sp}(2k,2n)}$. If a tangent vector
$Z\in{\mathrm{T}_{X}}{\mathrm{Sp}(2k,2n)}$ satisfies $\|Z\|_{2}<1$, then $X+Z$
has an SR decomposition.

###### Proof.

We will show that all leading minors of even dimension of the matrix
$P_{2k}(X+Z)^{T}J_{2n}(X+Z)P_{2k}^{T}$ are nonzero. Taking (2), the
orthogonality of $P_{2k}$, and (10) into account, we obtain

$\displaystyle P_{2k}(X\\!+\\!Z)^{T}J_{2n}(X\\!+\\!Z)P_{2k}^{T}$
$\displaystyle=P_{2k}(X^{T}\\!J_{2n}X\\!+\\!Z^{T}\\!J_{2n}X\\!+\\!X^{T}\\!J_{2n}Z\\!+\\!Z^{T}\\!J_{2n}Z)P_{2k}^{T}$
$\displaystyle=P_{2k}(J_{2k}+Z^{T}\\!J_{2n}Z)P_{2k}^{T}$
$\displaystyle=P_{2k}J_{2k}P_{2k}^{T}\big{(}I_{2k}-(P_{2k}J_{2k}P_{2k}^{T})(P_{2k}Z^{T}\\!J_{2n}ZP_{2k}^{T})\big{)}$
$\displaystyle=\hat{J}_{2k}\big{(}I_{2k}-\hat{J}_{2k}(P_{2k}Z^{T}\\!J_{2n}ZP_{2k}^{T})\big{)}.$

Because of the special structure of $\hat{J}_{2k}$, the nonsingularity of

$\Big{(}P_{2k}(X+Z)^{T}J_{2n}(X+Z)P_{2k}^{T}\Big{)}_{1:2j,1:2j}\\!=\\!\Big{(}\hat{J}_{2k}\big{(}I_{2k}-\hat{J}_{2k}(P_{2k}Z^{T}\\!J_{2n}ZP_{2k}^{T})\big{)}\Big{)}_{1:2j,1:2j}$

is the same as that of
$\big{(}I_{2k}-\hat{J}_{2k}(P_{2k}Z^{T}\\!J_{2n}ZP_{2k}^{T})\big{)}_{1:2j,1:2j}$
for $j=1,\ldots,k$. The latter in turn can be shown to be true. Indeed, since
$\|J_{2n}\|_{2}=\|\hat{J}_{2k}\|_{2}=1$ and $P_{2k}$ is orthogonal, we have

$\displaystyle\Big{\|}\big{(}\hat{J}_{2k}(P_{2k}Z^{T}\\!J_{2n}ZP_{2k}^{T})\big{)}_{1:2j,1:2j}\Big{\|}_{2}$
$\displaystyle=\Big{\|}(\hat{J}_{2k})_{1:2j,1:2j}(P_{2k}Z^{T}\\!J_{2n}ZP_{2k}^{T})_{1:2j,1:2j}\Big{\|}_{2}$
$\displaystyle\leq\Big{\|}(P_{2k}Z^{T}\\!J_{2n}ZP_{2k}^{T})_{1:2j,1:2j}\Big{\|}_{2}$
$\displaystyle\leq\big{\|}P_{2k}Z^{T}\\!J_{2n}ZP_{2k}^{T}\big{\|}_{2}\leq\big{\|}Z\big{\|}_{2}^{2}<1$

for $j=1,\ldots,k$. Thus, in view of Theorem 3.1, the proof is complete. ∎

##  4 Riemannian gradient method with non-monotone line search

We now present the Riemannian gradient method with non-monotone line search
for solving the constrained optimization problem (1). Starting with an initial
guess $X_{0}\in{\mathrm{Sp}(2k,2n)}$, this method generates a sequence of
iterates $\\{X_{i}\\}$ using a search direction
$-\operatorname*{grad}f(X_{i})\in\mathrm{T}_{X_{i}}{\mathrm{Sp}(2k,2n)}$ as

$X_{i+1}=\mathcal{R}_{X_{i}}(-\tau_{i}\,\operatorname*{grad}f(X_{i})),$

where $\mathcal{R}_{X_{i}}$ is one of the retractions defined in Section 3,
and $\tau_{i}>0$ is an appropriate step size. We summarize the resulting
Riemannian gradient method combined with the alternating Barzilai–Borwein
strategy [30] for the step size in Algorithm 4.

Algorithm 4 Riemannian gradient method for the optimization problem (1)

0: The cost function $f$, metric $g$ and retraction $\mathcal{R}$ on
${\mathrm{Sp}(2k,2n)}$, initial guess $X_{0}\in{\mathrm{Sp}(2k,2n)}$,
$\gamma_{0}>0$, $0<\gamma_{\min}<\gamma_{\max}$, $\beta,\delta\in(0,1)$,
$\alpha\in[0,1]$, $q_{0}=1$, $c_{0}=f(X_{0})$.

0: Sequence of iterates $\\{X_{i}\\}$.

1: for $i=0,1,2,\dots$ do

2: Compute $Z_{i}=-\operatorname*{grad}f(X_{i})$.

3: if $i>0$ then

4:
$\gamma_{i}=\left\\{\begin{array}[]{ll}\dfrac{\|W_{i-1}\|_{F}^{2}}{\left|\operatorname*{tr}(W_{i-1}^{T}Y_{i-1})\right|}&\text{for
odd }i,\\\\[8.53581pt]
\dfrac{\left|\operatorname*{tr}(W_{i-1}^{T}Y_{i-1})\right|}{\|Y_{i-1}\|_{F}^{2}}&\text{for
even }i,\end{array}\right.$where $W_{i-1}=X_{i}-X_{i-1}$ and
$Y_{i-1}=Z_{i}-Z_{i-1}$.

5: end if

6: Calculate the trial step size
$\gamma_{i}=\max\bigl{(}\gamma_{\min},\min(\gamma_{i},\gamma_{\max})\bigr{)}$.

7: Find the smallest integer $\ell$ such that the non-monotone condition
$f\big{(}\mathcal{R}_{X_{i}}(\tau_{i}Z_{i})\big{)}\leq
c_{m}+\beta\,\tau_{i}\,g\big{(}\operatorname*{grad}f(X_{i}),Z_{i}\big{)}$
holds, where $\tau_{i}=\gamma_{i}\,\delta^{\ell}$.

8: Set $X_{i+1}=\mathcal{R}_{X_{i}}(\tau_{i}Z_{i})$.

9: Update $q_{i+1}=\alpha q_{i}+1$ and $\displaystyle{c_{i+1}=\frac{\alpha
q_{i}}{q_{i+1}}c_{i}+\frac{1}{q_{i+1}}f(X_{i+1})}$.

10: end for

It has been shown in [1, Theorem 5.7] that every accumulation point $X_{*}$ of
the sequence $\\{X_{i}\\}$ generated by Algorithm 4 is a critical point of the
cost function $f$ in (1), i.e., $\operatorname*{grad}f(X_{*})=0$, no matter
which metric $g$ and which retraction $\mathcal{R}$ are employed provided that
$0_{X_{*}}\in\mathrm{T}_{X_{*}}{\mathrm{Sp}(2k,2n)}$ is in the interior of the
domain of $\mathcal{R}$. Note that due to Theorem 3.8, the new retraction (19)
satisfies this condition.

##  5 Applications

Algorithm 4 is presented without any specification of metric nor retraction.
As discussed in Sections 2 and 3, we can take either the canonical-like (C) or
Euclidean (E) metric and one of the retractions based on Cayley transform
(Cayley), quasi-geodesics (QGeo) and SR decomposition (SR). It is important to
note that they are independently constructed and can freely be combined. As a
result, we obtain totally six optimization schemes which are referred to as
CayleyC, CayleyE, QGeoC, QGeoE, SRC, and SRE. In this section, we present
different applications for optimization on the symplectic Stiefel manifold
${\mathrm{Sp}(2k,2n)}$ and compare the numerical performance of different
optimization schemes applied to various problems and scenarios. Three
quantities of most interest are the values of the cost function, the norms of
the Riemannian gradient of this function, and the feasibility violation of the
approximate optimal solution. Depending on the problem, other performance and
error measures will also be investigated. The numerical experiments are
performed on a standard desktop with Intel(R) Core(TM) i9-11900K (at 3.50GHz,
16MB Cache, 32GB RAM) running MATLAB R2021b under Ubuntu 22.04. The code
including different optimization schemes is available from
https://github.com/opt-gaobin/spopt.

Various parameters in Algorithm 4 need to be set. Most of them are fixed for
all problems except for the maximal number of iterations niter and the
tolerance gtol in the stopping criterion
$\|\operatorname*{grad}f(X_{i})\|_{F}\leq\texttt{gtol}$. For the backtracking
search, we set $\beta=1\mathrm{e}{-4}$, $\delta=1\mathrm{e}{-1}$,
$\gamma_{0}=1\mathrm{e}{-3}$ (if it is not otherwise specified),
$\gamma_{\min}=1\mathrm{e}{-15}$, $\gamma_{\max}=1$ for the symplectic
eigenvalue computation problem and $\gamma_{\max}=1\mathrm{e}\\!+\\!5$ for the
rest. For the non-monotone condition, we choose $\alpha=0.85$. Further, we use
the canonical-like metric $g_{c,\rho}$ with $\rho=1/2$ as recommended in [1,
Subsection 6.2]. Note also that Theorem 3.8 provides a sufficient condition on
the size of $\|\tau_{i}Z_{i}\|_{2}$ in step 8 of Algorithm 4 in the case of
using the SR retraction. In our experiments, we never experienced a problem
with the SR decomposition even without any adjustment of $\tau_{i}$ for this
existence purpose.

###  5.1 Symplectic target matrix problem

In this problem, one has to solve the following minimization problem

$\min_{X\in{\mathrm{Sp}(2n)}}f(X):=\|X-W\|^{2}_{F},$ (20)

where $W\in{\mathrm{Sp}(2n)}$ is given. In the _optimal control of
(symplectic) quantum gate_ , it is shown that one can realize a target/ideal
quantum gate with symplectic transformations by minimizing the distance
between the real and ideal gates which in turn can be simplified to solving
(20), where $W$ is a matrix representation of the gate [10, 31]. A more
general problem than (20), where the cost function is given in the form of a
uniformly weighted sum of distances to several symplectic targets, has been
investigated in [32]. An extension of this problem to the symplectic Stiefel
manifold ${\mathrm{Sp}(2k,2n)}$ and a general target matrix
$W\in\mathbb{R}^{2n\times 2k}$, termed as the _nearest symplectic matrix
problem_ , has been considered as a test problem in [1]. Here, we examine two
test examples. In the first one, the SUM gate is used as the symplectic target
gate, which is given by

$W=\left[\begin{array}[]{c c c r}1&0&0&0\\\ 1&1&0&0\\\ 0&0&1&-1\\\
0&0&0&1\end{array}\right],$

see [31] and references therein. We run two experiments with different initial
guesses and the same maximal number of iterations $\texttt{niter}=500$ and
tolerance $\texttt{gtol}=1\mathrm{e}{-12}$. First, we choose $X_{0}=I_{4}$ as
the starting point. Figure 1 demonstrates that all six optimization schemes
converge to the global minimizer $W$. Notably, the ones with the Euclidean
metric get to the target after only a few steps while those with the
canonical-like metric need more iterations. Moreover, the symplecticity of
their iterates is truly attained which makes the corresponding plots disappear
in the right subfigure in Figure 1.

For another initial guess
$X_{0}=\operatorname*{diag}(1.728,-1.2,1/1.728,-1/1.2)$, the convergence
results and the feasibility violation are presented in Figure 2. One can see
that all optimization schemes converge not to the minimizer but to a critical
point which, in view of [31, Theorem 4.2], is a saddle point. The effect of
different metrics on the numerical behavior, though not considerably as the
first run, but still can be seen: the ones with the Euclidean metric
apparently converge faster.

In the second example, we employ artificial data in order to check different
scenarios and problem sizes. To this end, we choose the symplectic target
matrix $W=\left[\begin{smallmatrix}I&0\\\ V&I\end{smallmatrix}\right]$ and the
initial iterate $X_{0}=\left[\begin{smallmatrix}I&Y\\\
0&I\end{smallmatrix}\right]$ with randomly generated matrices $V,Y\in{\cal
S}_{\mathrm{sym}}(200)$, see [33, Lemma 2.1]. We run Algorithm 4 for
$\texttt{niter}=1000$ iterations with the tolerance
$\texttt{gtol}=1\mathrm{e}{-10}$. Figure 3 shows that the difference between
the optimization schemes based on the canonical-like and Euclidean metric is
even severer: the ones with the canonical-like metric do not provide
reasonable results after $1000$ iterations, while the others reach the
tolerance after less than $50$ iterations. Comparing the convergence results
for the latter in the lower row of Figure 3, we observe that the SRE scheme
has better performance than CayleyE and QGeoE. Replacing the initial guess
$X_{0}$ with another symplectic matrices listed in [33, Lemma 2.1], we obtain
the similar results.

These experiments apparently tell us that if the cost function expresses the
distance resulting from the Frobenius norm, which can also be thought of as
the Euclidean norm in the corresponding matrix space, the Riemannian
optimization schemes with the Euclidean metric are preferred. Moreover, even
in the convergent case, the computed critical point can be a saddle point
which is far from being satisfactory.

Figure 1: Optimal control of the SUM gate with convergence to the global
minimizer: a comparison of the cost function values (left), Riemannian
gradient norms (middle), and the feasibility violation (right) for different
optimization schemes.  Figure 2: Optimal control of the SUM gate with
convergence to a saddle point: a comparison of the cost function values
(left), Riemannian gradient norms (middle), and the feasibility violation
(right) for different optimization schemes.  Figure 3: A symplectic target
matrix problem: a comparison of the cost function values (left), Riemannian
gradient norms (middle), and the feasibility violation (right) for six
different optimization schemes (upper row) and three schemes with the
Euclidean metric (lower row).

###  5.2 Computing symplectic eigenvalues

Originated from Williamson’s work [34], for a $2n\times 2n$ symmetric
positive-definite (SPD) real matrix $A$, there exists a matrix
$S\in{\mathrm{Sp}(2n)}$ such that

$S^{T}AS=\begin{bmatrix}D&0\\\ 0&D\end{bmatrix},$ (21)

where $D=\operatorname*{diag}(d_{1},\ldots,d_{n})$ is a diagonal matrix with
positive diagonal entries. The right-hand side of (21) is termed as
_Williamson’s diagonal form_ of $A$, and the positive numbers
$d_{1},\ldots,d_{n}$ are referred to as _symplectic eigenvalues_ of $A$. For
ease of later argument, the symplectic eigenvalues are always numbered in the
nondecreasing order, i.e., $d_{1}\leq\cdots\leq d_{n}$. Note that the
symplectic eigenvalues of $A$ differ from the standard eigenvalues of $A$ but
strongly relate to the eigenvalues of the Hamiltonian matrix $J_{2n}A$ or the
symmetric/skew-symmetric matrix pencil $A-\lambda J_{2n}$, and the Hermitian
matrix $\mathrm{i}A^{1/2}J_{2n}A^{1/2}$ with $\mathrm{i}=\sqrt{-1}$ and
$A^{1/2}$ being the symmetric square root of $A$ or the Hermitian pencil
$A-\lambda\mathrm{i}J_{2n}$, which have been intensively investigated in [35,
36, 12, 37], to name a few. A  _pair of symplectic eigenvectors_
$u,v\in\mathbb{R}^{n}\backslash\\{0\\}$ associated with a symplectic
eigenvalue $d$ of $A$ are those that satisfy

$Au=d\,J_{2n}v,\qquad Av=-d\,J_{2n}u.$

Symplectic eigenvalues and eigenvectors can be numerically computed using a
symplectic Lanczos method via the connection with so-called positive-definite
Hamiltonian matrices [35] or by solving a trace minimization problem

$\min_{X\in{\mathrm{Sp}(2k,2n)}}f(X):=\operatorname*{tr}(X^{T}AX)$ (22)

using a Riemannian optimization method [13]. It has been shown in [38, 12]
that the minimal value of $f$ in (22) is twice the sum of $k$ smallest
symplectic eigenvalues of $A$. It is worth to note that this result can also
be derived based on the trace minimization theorems, called the Ky-Fan
theorem, for standard eigenvalues of real symmetric positive-semidefinite
pencils [39, 40] or its extended complex Hermitian version [41]. Recently, in
[14, 42], the notion of symplectic eigenvalues and its trace minimization
theorem have been extended to a special class of symmetric positive-
semidefinite (SPSD) matrices which have symplectic null space. Symplectic
eigenvalues find applications in quantum mechanics, optics, stability analysis
of gyroscopic systems, and in quantization process of superconducting networks
modeled by Hamiltonian systems, see [38, 43, 11, 44, 45, 46].

In the optimization approach, to compute the $k$ smallest symplectic
eigenvalues with $k$ pairs of symplectic eigenvectors for $1\leq k\leq n$, one
first finds a minimizer $X_{*}\in{\mathrm{Sp}(2k,2n)}$ of the cost function
$f$ in (22) and then diagonalizes $X_{*}^{T}AX_{*}$ by an orthosymplectic
matrix $K\in\mathbb{R}^{2k\times 2k}$ such that
$K^{T}X_{*}^{T}AX_{*}K=\operatorname*{diag}(\Lambda,\Lambda)$ with
$\Lambda=\operatorname*{diag}(\lambda_{1},\ldots,\lambda_{k})$. It has been
shown in [13] that $\lambda_{1},\ldots,\lambda_{k}$ and the pair of the $j$-th
and $(j+k)$-th columns of $X_{*}K$, $j=1,\ldots,k$, are the sought symplectic
eigenvalues and associated symplectic eigenvector pairs of $A$, respectively.

In the first example, in view of [14], we construct an SPSD matrix

$A=J_{2n}Q\begin{bmatrix}D&0\\\ 0&D\end{bmatrix}(J_{2n}Q)^{T}$

with a symplectic null space by setting  D = diag(0,…,0,m+1,m+2,…,n-1,n)  with
$0<m<k<n$ and $Q=KL(n/5,1.2,-\sqrt{n/5})$, where $L(l,c,d)$ is the (non-
orthogonal) symplectic Gauss transformation

$L(l,c,d)=\begin{bmatrix}I_{l-2}&&&&&&&\\\ &c&&&&&d&\\\ &&c&&&d&&\\\
&&&I_{n-l}&&&&\\\ &&&&I_{l-2}&&&\\\ &&&&&c^{-1}&&\\\ &&&&&&c^{-1}&\\\
&&&&&&&I_{n-l}\\\ \end{bmatrix}$

and $K\\!=\\!\left[\\!\begin{smallmatrix}\enskip\;(U)&\operatorname*{Im}(U)\\\
-\operatorname*{Im}(U)&(U)\end{smallmatrix}\\!\right]{\color[rgb]{0.62,0.0,0.77}\definecolor[named]{pgfstrokecolor}{rgb}{0.62,0.0,0.77}\\!}{\in{\mathrm{Sp}(2n)}}$
with a randomly generated unitary matrix $U\\!\in\\!\mathbb{C}^{n\times n}$.
Obviously, $A$ is of rank $2(n-m)$ and its null space has dimension $2m$. In
our experiments, we take $n=1000$, $k=5$, and $m=2$. In this case, the $5$
smallest symplectic eigenvalues of $A$ are $0,0,3,4,5$. The accuracy of the
computed symplectic eigenvalues $\tilde{d}_{j}$, $j=1,\ldots,k$, is verified
by the $l_{1}$-norm error defined as $\sum_{j=1}^{k}|d_{j}-\tilde{d}_{j}|$.
For all six optimization schemes, we run Algorithm 4 with
$\texttt{niter}=5000$ iterations and the tolerance
$\texttt{gtol}=1\mathrm{e}-12$. In Table 2, we present the computed symplectic
eigenvalues, the $l_{1}$-norm errors, and the errors in the minimal value of
the cost function, i.e., ${\rm err}_{f}=|f(X_{5000})-24|$. One can see that
the optimization schemes reach the same accuracy as the MATLAB function eigs
applied to the Hermitian pencil $A-\lambda\mathrm{i}J_{2n}$ with the default
tolerance $1\mathrm{e}{-14}$. Certainly, as less direct, these schemes are
more computationally expensive. Nevertheless, an obviously seen advantage over
eigs is that they provide purely real symplectic eigenvalues. In Figure 4, we
also present the history of the iterations consisting of the errors in the
minimal value of the cost function, the norms of its Riemannian gradient, and
the feasibility violation.

Figure 4: Symplectic eigenvalue computation for an SPSD matrix: a comparison
of the cost function values (left), Riemannian gradient norms (middle), and
the feasibility violation (right) for different optimization schemes.

Regarding efficiency comparison, in Table 3, we report the average time
consumed in one iteration step by different optimization schemes with the
aforementioned setting and with $k=n=1000$. One can see that the schemes using
the Cayley retraction are fastest in most cases. Moreover, compared to the
schemes with the canonical-like metric, those based on the Euclidean metric
tend to be slightly faster when $k\ll n$ but are obviously slower when $k$
approaches $n$. The main reason is the difference in the orthogonal
projections for computing the Riemannian gradients, c.f., Table 2.2: while for
the canonical-like metric, one has to multiply matrices of size scalable with
$2n$, a Lyapunov matrix equation with coefficients of size $2k\times 2k$ is
required to solve for the Euclidean metric. Therefore, for $k\ll n$, solving
such an equation can be much faster while when $k=n$, the cost is definitely
more expensive than matrix-matrix multiplication.

Table 2: Symplectic eigenvalue computation for an SPSD matrix: 5 computed smallest symplectic eigenvalues, $l_{1}$-errors, and the errors in the minimal value of the cost function.  scheme | Canonical-like (C) | Euclidean (E)
---|---|---
symplectic eigenvalues |  $\;l_{1}$-error | $\quad{\rm err}_{f}\quad$ | symplectic eigenvalues |  $\;l_{1}$-error | $\quad{\rm err}_{f}\quad$
| 0.000000000000000 |  |  | 0.000000000000000 |  |
| 0.000000000000082 |  |  | 0.000000000000049 |  |
Cayley | 2.999999999999981 | $1.42\mathrm{e}{-13}$ | $8.53\mathrm{e}{-14}$ | 3.000000000000000 | $1.14\mathrm{e}{-13}$ | $2.06\mathrm{e}{-13}$
| 3.999999999999988 |  |  | 3.999999999999991 |  |
| 5.000000000000029 |  |  | 5.000000000000055 |  |
| 0.000000000000000 |  |  | 0.000000000000000 |  |
| 0.000000000000057 |  |  | 0.000000000000063 |  |
QGeo | 2.999999999999885 | $3.71\mathrm{e}{-13}$ | $4.23\mathrm{e}{-13}$ | 2.999999999999929 | $1.86\mathrm{e}{-13}$ | $8.52\mathrm{e}{-14}$
| 3.999999999999878 |  |  | 3.999999999999998 |  |
| 4.999999999999924 |  |  | 5.000000000000051 |  |
| 0.000000000000000 |  |  | 0.000000000000000 |  |
| 0.000000000000066 |  |  | 0.000000000000056 |  |
SR | 2.999999999999979 | $1.06\mathrm{e}{-13}$ | $1.67\mathrm{e}{-13}$ | 2.999999999999991 | $1.09\mathrm{e}{-13}$ | $7.83\mathrm{e}{-14}$
| 3.999999999999993 |  |  | 3.999999999999986 |  |
| 5.000000000000012 |  |  | 5.000000000000031 |  |
| $\qquad(18+60\,\mathrm{i})\mathrm{e}{-15}\;\,$ |  |  |  |  |
| $\qquad(16+62\,\mathrm{i})\mathrm{e}{-15}\;\,$ |  |  |  |  |
eigs | $3+((15-\enskip 3\,\mathrm{i})\mathrm{e}{-15})$ | $2.19\mathrm{e}{-13}$ |  |  |  |
| $4-((14-11\,\mathrm{i})\mathrm{e}{-15})$ |  |  |  |  |
| $5+((59-\enskip 8\,\mathrm{i})\mathrm{e}{-15})$ |  |  |  |  |
Table 3: Symplectic eigenvalue computation for an SPSD matrix: time consumed in seconds per step by different optimization schemes. |  |  | Cayley | QGeo | SR
---|---|---|---|---|---
$n=1000$ |  $k=5$ | Canonical-like (C) |  8.50$\mathrm{e}{-3}$ |  9.38$\mathrm{e}{-3}$ |  9.21$\mathrm{e}{-3}$
Euclidean (E) |  8.75$\mathrm{e}{-3}$ |  7.90$\mathrm{e}{-3}$ |  7.72$\mathrm{e}{-3}$
$k=1000$ | Canonical-like (C) |  2.31$\mathrm{e}{+0}$ |  8.01$\mathrm{e}{+0}$ |  9.88$\mathrm{e}{+0}$
Euclidean (E) |  2.63$\mathrm{e}{+1}$ |  3.64$\mathrm{e}{+1}$ |  3.18$\mathrm{e}{+1}$

Regarding efficiency comparison, in Table 3, we report the average time
consumed in one iteration step by different optimization schemes with the
aforementioned setting and with $k=n=1000$. One can see that the schemes using
the Cayley retraction are fastest in most cases. Moreover, compared to the
schemes with the canonical-like metric, those based on the Euclidean metric
tend to be slightly faster when $k\ll n$ but are obviously slower when $k$
approaches $n$. The main reason is the difference in the orthogonal
projections for computing the Riemannian gradients, c.f., Table 2.2: while for
the canonical-like metric, one has to multiply matrices of size scalable with
$2n$, a Lyapunov matrix equation with coefficients of size $2k\times 2k$ is
required to solve for the Euclidean metric. Therefore, for $k\ll n$, solving
such an equation can be much faster while when $k=n$, the cost is definitely
more expensive than matrix-matrix multiplication.

In the second example, we consider a mechanical system which is used in
analysis of vibration and frequency response of wire saws [47]. For such a
system, the SPD matrix is given by $A=J_{2n}H$ with the Hamiltonian matrix

$H=\begin{bmatrix}-\frac{1}{2}GM^{-1}&\frac{1}{4}GM^{-1}G-K\\\
M^{-1}&-\frac{1}{2}M^{-1}G\end{bmatrix},$

where $M,G,K$ are the mass, damping, and stiffness matrices of the underlying
mechanical system. Following the setting in [13], we are lead to the
minimization problem (22) with $n=2000$ and $k=5$. As we do not know the exact
symplectic eigenvalues, in Figure 5, we only present the history comparison of
the six optimization schemes when Algorithm 4 is run with
$\texttt{niter}=100000$ iterations and the tolerance
$\texttt{gtol}=1\mathrm{e}{-6}$. We see that the cost function seems to
stagnate even though the Riemannian gradients are far from being zero.

Figure 5: Symplectic eigenvalue computation for a wire saw model: a comparison
of the cost function values (left), Riemannian gradient norms (middle), and
the feasibility violation (right) for different optimization schemes.

Though with examples in this subsection, there is no considerable difference
in the final values of the cost function, we can easily observe that the
feasibility violation varies quite a bit for different schemes. As expected,
those using the SR retraction apparently maintain the symplecticity constraint
the best, while the schemes based on the quasi-geodesic retraction perform the
worst. This is due to the fact that the symplecticity is retrieved by the SR
decomposition at every iteration. This feature prevents the accumulative
errors in the constraint, which can be seen in the quasi-geodesic case. This
argument is actually independent of the cost function and therefore holds true
for any problem. Moreover, these tests also reconfirm an observation on the
constraint violation for the optimization schemes based on the Cayley and
quasi-geodesic retractions made in [4].

###  5.3 Symplectic model reduction of Hamiltonian systems

The third minimization problem arises in structure-preserving model reduction
of nonlinear Hamiltonian systems given by

$\dot{x}=J_{2n}\nabla_{x}H(x),\qquad x(0)=x_{0},$ (23)

where $x\in\mathbb{R}^{2n}$ is the state vector, $x_{0}\in\mathbb{R}^{2n}$ is
the initial vector, $H:\mathbb{R}^{2n}\to\mathbb{R}$ is the continuously
differentiable Hamiltonian function describing the internal energy of the
system, $\nabla_{x}H(x)$ is the Euclidean gradient of $H$ with respect to $x$,
and $J_{2n}$ is the structure matrix which describes the interconnection of
energy storage elements. For such systems, the Hamiltonian $H$ is a first
integral since it remains invariant along the solution of (23). Hamiltonian
systems arise in a wide range of applications including mechanical systems,
molecular dynamics, network design and electromagnetic field simulation, e.g.,
[48, 49, 50].

The goal of model reduction is to approximate the full-order model (FOM) (23)
by a reduced-order model (ROM) which preserves the Hamiltonian structure

$\dot{\tilde{x}}=J_{2k}\nabla_{\tilde{x}}\tilde{H}(\tilde{x}),\qquad\tilde{x}(0)=\tilde{x}_{0},$
(24)

where $\tilde{x}\in\mathbb{R}^{2k}$ is the reduced state, and
$\tilde{H}:\mathbb{R}^{2k}\to\mathbb{R}$ is the reduced Hamiltonian with $k\ll
n$. Note that the preservation of the Hamiltonian structure during the
reduction process ensures the conservation of energy of the reduced-order
model. In the last decade, structure- and energy-preserving model reduction of
Hamiltonian systems has attracted a lot of attention and several model
reduction methods have been developed for such systems [15, 16, 51, 17, 52,
53, 54].

One of the most popular model reduction approaches for nonlinear Hamiltonian
systems is the PSD method presented first in [15]. It is based on collecting
the snapshots $A=[x(t_{1}),\ldots,x(t_{s})]\in\mathbb{R}^{2n\times s}$ of the
Hamiltonian system (23) and determining a projection matrix
$U\in{\mathrm{Sp}(2k,2n)}$ which solves the following constrained minimization
problem

$\min_{X\in{\mathrm{Sp}(2k,2n)}}f(X):=\|A-XX^{+}A\|_{F}^{2},$ (25)

where $X^{+}=J_{2k}^{T}\,X^{T}J_{2n}$ is a symplectic inverse of $X$. Although
the existence of a minimizer of (25) is still unknown, the cost function is
smooth, bounded from below, and with a monotonically decreasing line search
strategy and an experimental convergence recognized by, e.g., the distance
between two consecutive iterates, one might hope that the value of the cost
function at the computed point is close to the infimum. The solution of (23)
is then approximated by $x\approx U\tilde{x}$, where
$\tilde{x}\in\mathbb{R}^{2k}$ is a coordinate vector of the approximation with
respect to the basis $U$. Replacing $x$ with $U\tilde{x}$ and multiplying the
resulting equation from the left with $U^{+}$, we obtain the reduced-order
system (24) with the reduced initial vector $\tilde{x}_{0}=U^{+}x_{0}$ and the
reduced Hamiltonian $\tilde{H}(\tilde{x})=H(U\tilde{x})$. It has been shown in
[15] that the error in the Hamiltonian given by $\Delta
H(t)=H(x(t))-\tilde{H}(\tilde{x}(t))$ is constant for all $t\in\mathbb{R}$.
Moreover, if $x_{0}\in\operatorname*{im}(U)$, then $\Delta H(t)\equiv 0$. This
implies that the reduced-order Hamiltonian system (24) preserves the energy.
Another important property of the symplectic projection is the preservation of
the stability of equilibrium points, see [15, 16] for details.

Due to the non-convexity and unboundedness of the feasible set
${\mathrm{Sp}(2k,2n)}$, the minimization problem (25) was considered to be
difficult to solve. By imposing additional orthogonality constraint
$X^{T}X=I$, different algorithms have been developed in [15] for computing
suboptimal solutions. A similar constraint is also used in [16] for a greedy
algorithm applied to parametric Hamiltonian systems. Another approach for
generating a non-orthonormal symplectic basis matrix has been presented in
[17] which is based on the SVD-like decomposition [55] of the snapshot matrix.
Unlike the standard SVD, where the optimality of the approximation is well
known, the theoretical results derived in [17] does not infer a similar claim.

Here, we employ the Riemannian optimization algorithm for structure-preserving
model reduction of Hamiltonian systems. The key potential advantage of this
approach is the possibility to reach an optimal solution instead of a
suboptimal one. Note that since the cost function $f$ in (25) satisfies the
homogeneity property $f(XS)\\!=\\!f(X)$ for any $S\in\mathrm{Sp}(2k)$, one
could also compute a symplectic projection matrix by solving the minimization
problem on the symplectic Grassmann manifold [2, 19]. A comparison of these
two Riemannian optimization approaches in the context of PSD is out of the
scope of this paper.

To compute symplectic reduced bases using Algorithm 4, initialization is
required. To reduce the risk of convergence to a saddle point as illustrated
in Subsection 5.1, we take the symplectic matrix $X_{0}$ produced by the
cotangent lift method [15]. It is a block diagonal orthosymplectic matrix
$X_{0}=\operatorname*{diag}(\hat{X},\hat{X})$, where the columns of
$\hat{X}\in\mathbb{R}^{n\times k}$ are the left singular vectors of the
combined snapshot matrix $\left[[I_{n},\,0]A,\;[0,\,I_{n}]A\right]$. As we
start with a suboptimal solution, it is expected that the trial step size
$\gamma_{i}$ should generally be small, and therefore, we set
$\gamma_{0}=1\mathrm{e}{-8}$. As will be shown below, the optimization method
always improves this result and yields smaller approximation errors in model
reduction.

We also compare our optimization-based model reduction methods with other
existing reduction techniques for Hamiltonian systems. Among the PSD
approaches proposed in [15], we choose the cotangent lift method due to its
good performance, where the balance between the simplicity and the accuracy is
taken into account. Moreover, whenever possible, we include the PSD SVD-like
decomposition method from [17] in our comparison experiments.

As we are specially interested in the energy conservation in the course of
model reduction of Hamiltonian systems, we employ the Crank–Nicolson
integration method with a constant time step size $h_{t}$. It has been shown
in [56] that this method is of second order and that it delivers a discrete
solution whose mass and energy are conserved. The Crank–Nicolson scheme, being
implicit, requires the numerical solution of nonlinear systems in each time
step. For this purpose, we use in our experiments the MATLAB function fsolve.

For nonlinear dynamical systems, one has to additionally approximate the
nonlinear term in order to maintain the benefit of smaller order of the ROM
(24). This task for nonlinear Hamiltonian systems is even more challenging
since the approximation must be done in such a way that it does not destroy
the Hamiltonian structure which guarantees the energy conservation. For this
purpose, different approaches can be used.

  * 1)

The first one is the PSD-DEIM (also termed as SDEIM in [15]), which is a
combination of the PSD proposed in [15] and the DEIM developed in [20].
Assuming that the Hamiltonian function in (23) has the form

$H(x)=\frac{1}{2}x^{T}Mx+h(x)$

with an SPD matrix $M\in\mathbb{R}^{2n\times 2n}$ and a nonlinear function
$h:\mathbb{R}^{2n}\to\mathbb{R}$, the gradient
$\nabla_{\tilde{x}}\tilde{H}(\tilde{x})$ in the ROM (24) is then approximated
by

$\nabla_{\tilde{x}}\tilde{H}(\tilde{x})\approx
U^{T}MU\tilde{x}+U^{T}V(P^{T}V)^{-1}P^{T}\nabla_{x}h(U\tilde{x}),$ (26)

where $V\in\mathbb{R}^{2n\times m}$ is a DEIM basis matrix and
$P=[e_{i_{1}},\ldots,e_{i_{m}}]$ is a selector matrix associated with an index
set $\\{i_{1},\ldots,i_{m}\\}$ determined by a greedy procedure applied to
$V$.

  * 2)

The second approach is the structure-preserving method developed for port-
Hamiltonian systems in [21]. In this case, one uses the approximation

$\nabla_{\tilde{x}}\tilde{H}(\tilde{x})\approx
U^{T}MU\tilde{x}+U^{T}V(P^{T}V)^{-1}P^{T}\nabla_{x}h(P(V^{T}P)^{-1}V^{T}U\tilde{x}).$
(27)

The resulting model reduction method is referred to as structure-preserving
PSD-DEIM.

The key point of the approximations (26) and (27) is that only a small number
of selected components of the nonlinear term are needed to be evaluated. In
our implementation, this number $m$ is approximately set to $2.5k$. One can
observe that the PSD-DEIM approach (26) focuses more on the task of
approximation and actually the resulting ROM deviates from being a Hamiltonian
system. Nevertheless, the rate of this deviation is shown to be bounded and
this upper bound depends on the approximation quality of the nonlinear term
[15, Theorem 5.1]. On the contrary, the approximation (27) seems to devote
more to the structure preservation task. Indeed, it provides a reduced-order
Hamiltonian model with an approximate Hamiltonian function

$\tilde{H}(\tilde{x})=\frac{1}{2}\tilde{x}^{T}U^{T}MU\tilde{x}+h(P(V^{T}P)^{-1}V^{T}U\tilde{x}).$

In our opinion, however, the approximation $h(P(V^{T}P)^{-1}V^{T}U\tilde{x})$
can unexpectedly cause large errors in the Hamiltonian function and the state.
In Subsection 5.3.5, we present a comparison of these two approaches. In our
experiments, for model reduction of nonlinear-systems, we use the PSD-DEIM
method unless stated otherwise.

To verify the approximation properties of ROMs, we consider the relative
errors in the state vector and the energy given by

${\rm
RE}_{x}=\dfrac{\|x-\tilde{x}\|_{L^{2}(0,T)}}{\|x\|_{L^{2}(0,T)}},\qquad\quad{\rm
RE}_{H}=\dfrac{\|H(x(\cdot))-\tilde{H}(\tilde{x}(\cdot))\|_{L^{2}(0,T)}}{\|H(x(\cdot))\|_{L^{2}(0,T)}},$
(28)

respectively, where $\|\cdot\|_{L^{2}(0,T)}$ denotes the $L^{2}$-norm of the
corresponding function. Moreover, for each model, we choose one ROM of certain
dimension and present the relative state vector error
$\|x(t)-\tilde{x}(t)\|/\mu_{[0,T]}(\|x\|)$, where $\mu_{[0,T]}(\|x\|)$ denotes
the mean of the function $\|x(t)\|$ over the time interval $[0,T]$, and the
relative energy error $|H(x(t))-\tilde{H}(\tilde{x}(t))|/|H(x(0))|$ versus the
time. We also expose the so-called _average accelerating factor_ (a.a.f.) of
MOR which is the ratio of the simulation time for the FOM (23) over the
average of those for the ROMs (24) computed by different model reduction
methods.

In the following, we present the results of various numerical experiments. In
each example, we first introduce a model in the form of PDEs and its
associated energy. This model is then spatially discretized using an
appropriate method to give a finite-dimensional Hamiltonian system of the form
(23) with the corresponding semi-discrete energy. Furthermore, physical and
numerical parameters are set up.

In the preliminary tests, we have compared different optimization schemes
applied to all models similarly to Subsection 5.2. Since the symplecticity
constraint is important in preserving the Hamiltonian structure, we choose the
ones that maintain this constraint better. Moreover, among those, we also want
to see the possible difference resulting from the use of different metrics.
Therefore, for simplicity, we restrict our comparisons in Subsections
5.3.1–5.3.4 to the CayleyC and SRE schemes. In addition, general discussion is
given in Subsection 5.3.5.

####  5.3.1 Linear wave equation

As a first example, we consider a 1D linear wave equation subject to the
periodic boundary conditions

$\begin{array}[]{rclcl}z_{tt}&=&c^{2}\,z_{\xi\xi}&\quad\text{ in
}&(0,T)\times(a,b),\\\ z(t,a)&=&z(t,b)&\quad\text{ in }&(0,T),\\\
z(0,\xi)&=&z_{0}(\xi),\enskip z_{t}(0,\xi)=z_{1}(\xi)&\quad\text{ in
}&(a,b),\end{array}$ (29)

where constant $c>0$ is the wave speed and $z$ is the unknown function of time
and space, see [15]. The Hamiltonian function associated with (29) is given by

$\mathcal{H}(z)=\int_{a}^{b}\frac{1}{2}z_{t}^{2}+\frac{c^{2}}{2}z_{\xi}^{2}\>{\rm
d}\xi.$ (30)

Using a finite difference discretization method on a uniform spatial grid
$\xi_{j}=a+jh_{\xi}$, $j=0,...,n$, with the spatial step size
$h_{\xi}=(b-a)/n$ and setting $q_{j}(t)=z(t,\xi_{j})$,
$p_{j}(t)=z_{t}(t,\xi_{j})$, $q=[q_{1},\ldots,q_{n}]^{T}$, and
$p=[p_{1},\ldots,p_{n}]^{T}$, we obtain the linear Hamiltonian system

$\displaystyle\begin{bmatrix}\dot{q}\\\ \dot{p}\end{bmatrix}$
$\displaystyle=J_{2n}\begin{bmatrix}-c^{2}D_{\xi\xi}&0\\\
0&I\end{bmatrix}\begin{bmatrix}q\\\ p\end{bmatrix}$ (31)

with the quadratic Hamiltonian function

$H(q,p)=\sum_{j=1}^{n}\frac{p_{j}^{2}}{2}+\frac{c^{2}(q_{j+1}-q_{j})^{2}}{2h_{\xi}^{2}}+\frac{c^{2}(q_{j}-q_{j-1})^{2}}{2h_{\xi}^{2}}=\frac{1}{2}p^{T}p-\frac{c^{2}}{2}q^{T}D_{\xi\xi}\,q,$

where $D_{\xi\xi}$ is the three-point finite difference matrix approximating
the second-order spatial differential operator. Note that $h_{\xi}H(q,p)$
gives the spatial discretization of the Hamiltonian function (30). In our
tests, we chose $c=0.1$, $a=0$, $b=1$, $h_{\xi}=0.002$, $T=50$, $h_{t}=0.01$,
and the initial functions $z_{1}(\xi)\equiv 0$ and
$z_{0}(\xi)=\phi(10|\xi-\tfrac{1}{2}|)$ with the cubic spline function

$\phi(\eta)=\left\\{\begin{array}[]{lcl}1-\frac{3}{2}\eta^{2}+\frac{3}{4}\eta^{3}&\text{
if }&0\leq\eta\leq 1,\\\\[2.84526pt] \frac{1}{4}(2-\eta)^{3}&\text{ if
}&1<\eta\leq 2,\\\\[2.84526pt] 0&\text{ if }&2<\eta.\end{array}\right.$

The resulting system (31) has dimension $2n=1000$. For the PSD data, we
extract $s=500$ snapshots from the Crank–Nicolson solution of the FOM and
compute the ROMs of dimension $2k$ with $k=10,20,40,80$. The relative errors
defined in (28) and a.a.f. for the ROMs computed by the cotangent lift
(CotLift), the Riemannian optimization-based methods CayleyC and SRE with
$\texttt{iter}=1000$ iterations, and the SVD-like reduction method (SVD-like)
are presented in Table 4. Furthermore, Figure 6 shows the changes of the
relative state and energy errors with time for the ROM of dimension $2k=80$.

Table 4: Linear wave equation, $n=500$: model reduction results (the italic font means a warning of ill-conditioning received during the simulation of the reduced-order model). $k$ | CotLift | CayleyC | SRE | SVD-like | a.a.f.
---|---|---|---|---|---
$\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$ | $\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$ | $\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$ | $\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$
10 |  4.91$\mathrm{e}{-2}$ | $1.45\mathrm{e}{-12}$ |  4.90$\mathrm{e}{-2}$ | $1.46\mathrm{e}{-12}$ |  4.90$\mathrm{e}{-2}$ | $1.12\mathrm{e}{-12}$ |  1.78$\mathrm{e}{-1}$ | $1.28\mathrm{e}{-2}$ | 73
20 |  1.63$\mathrm{e}{-2}$ | $2.26\mathrm{e}{-12}$ |  1.61$\mathrm{e}{-2}$ | $2.74\mathrm{e}{-12}$ |  1.60$\mathrm{e}{-2}$ | $9.88\mathrm{e}{-13}$ |  7.07$\mathrm{e}{-2}$ | $1.54\mathrm{e}{-4}$ | 48
40 |  9.79$\mathrm{e}{-3}$ | $4.22\mathrm{e}{-12}$ |  9.15$\mathrm{e}{-3}$ | $7.32\mathrm{e}{-12}$ |  9.05$\mathrm{e}{-3}$ | $1.10\mathrm{e}{-12}$ | _5.71e $-$2_ | _1.23e $-$4_ | 36
80 |  6.28$\mathrm{e}{-3}$ | $1.56\mathrm{e}{-12}$ |  6.23$\mathrm{e}{-3}$ | $1.23\mathrm{e}{-11}$ |  6.23$\mathrm{e}{-3}$ | $2.86\mathrm{e}{-12}$ |  7.92$\mathrm{e}{-2}$ | $1.96\mathrm{e}{-5}$ |  3
Figure 6: Linear wave equation, $n=500$, $k=40$: relative errors in the state
vector (left) and the Hamiltonian function (right).

####  5.3.2 Sine-Gordon equation

The second example is the sine-Gordon equation with Dirichlet boundary
conditions

$\begin{array}[]{rclcl}z_{tt}&=&z_{\xi\xi}-\sin(z)&\quad\text{ in
}&(0,T)\times(a,b),\\\ z(t,a)&=&\phi_{a}(t),\quad
z(t,b)=\phi_{b}(t),&\quad\text{ in }&(0,T),\\\ z(0,\xi)&=&z_{0}(\xi),\enskip
z_{t}(0,\xi)=z_{1}(\xi)&\quad\text{ in }&(a,b).\end{array}$ (32)

Such an equation appears, for example, in differential geometry and in a wide
variety of physical applications. For the initial conditions

$z_{0}(\xi)=4\arctan\biggl{(}\exp\Bigl{(}\frac{\xi-\xi_{0}}{\sqrt{1-v^{2}}}\Bigr{)}\biggr{)},\qquad
z_{1}(\xi)=\frac{-4\,v\,\exp\bigl{(}\frac{\xi-\xi_{0}}{\sqrt{1-v^{2}}}\bigr{)}}{\sqrt{1-v^{2}}\Bigl{(}1+\exp\bigl{(}\frac{2(\xi-\xi_{0})}{\sqrt{1-v^{2}}}\bigr{)}\Bigr{)}},$

the sine-Gordon equation (32) admits the exact solitary solution

$z(t,\xi)=4\arctan\biggl{(}\exp\Bigl{(}\frac{\xi-\xi_{0}-vt}{\sqrt{1-v^{2}}}\Bigr{)}\biggr{)},$

where $0<v<1$ is the velocity of the solitary wave and $\xi_{0}>0$.

The finite difference method on a uniform spatial grid $\xi_{j}=jh_{\xi}$ for
$j=0,...,n+1$, with the spatial step size $h_{\xi}=(b-a)/(n+1)$ leads to the
nonlinear Hamiltonian system

$\begin{bmatrix}\dot{q}\\\
\dot{p}\end{bmatrix}=J_{2n}\left(\begin{bmatrix}-D_{\xi\xi}&0\\\
0&I\end{bmatrix}\begin{bmatrix}q\\\ p\end{bmatrix}+\begin{bmatrix}\psi(q)\\\
0\end{bmatrix}\right),$ (33)

where $D_{\xi\xi}$ is the three-point finite difference matrix,
$q=[q_{1},\ldots,q_{n}]^{T}$ with $q_{j}(t)=z(t,\xi_{j})$,
$p=[p_{1},\ldots,p_{n}]^{T}$ with $p_{j}(t)=z_{t}(t,\xi_{j})$, and

$\psi(q)=\bigl{[}\sin(q_{1})-\frac{\phi_{a}}{h_{\xi}^{2}},\sin(q_{2}),\ldots,\sin(q_{n-1}),\sin(q_{n})-\frac{\phi_{b}}{h_{\xi}^{2}}\bigr{]}^{T}.$

Note that $D_{\xi\xi}$ in (33) slightly differs from that in the semi-
discretized wave equation (31) due to the different boundary conditions. The
Hamiltonian function of (33) is given by

$H(q,p)=\frac{1}{2}p^{T}p-\frac{1}{2}q^{T}D_{\xi\xi}q+\sum_{i=1}^{n}\bigl{(}1\\!-\\!\cos(q_{i})\bigr{)}+\frac{\phi_{a}^{2}}{2h_{\xi}^{2}}-\frac{\phi_{a}q_{1}}{h_{\xi}^{2}}+\frac{\phi_{b}^{2}}{2h_{\xi}^{2}}-\frac{\phi_{b}q_{n}}{h_{\xi}^{2}}.$

In our experiments, we take $v=0.2$, $\xi_{0}=10$, $a=0$, $b=50$,
$h_{\xi}=0.025$, $n=1001$, $T=90$, and $h_{t}=0.05$. After the FOM simulation,
$s=450$ snapshots are uniformly extracted from the solution for constructing
the reduced bases of size $2k$ with $k=11,13,15,17$ by using the proposed
optimization schemes with $\texttt{iter}=6000$ iterations. Model reduction
results with detailed setting values are given in Table 5 and Figure 7.

Table 5: Sine-Gordon equation, $n=1001$: model reduction results. $k$ | CotLift | CayleyC | SRE | a.a.f.
---|---|---|---|---
$\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$ | $\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$ | $\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$
11 | $1.43\mathrm{e}{-2}$ | $4.11\mathrm{e}{-2}$ | $1.70\mathrm{e}{-2}$ | $9.70\mathrm{e}{+0}$ | $9.38\mathrm{e}{-3}$ | $1.37\mathrm{e}{-2}$ | 118
13 | $8.32\mathrm{e}{-3}$ | $1.69\mathrm{e}{-2}$ | $4.21\mathrm{e}{-3}$ | $3.97\mathrm{e}{-3}$ | $6.78\mathrm{e}{-3}$ | $3.72\mathrm{e}{-3}$ | 100
15 | $5.65\mathrm{e}{-3}$ | $9.12\mathrm{e}{-3}$ | $3.60\mathrm{e}{-3}$ | $7.41\mathrm{e}{-4}$ | $4.97\mathrm{e}{-3}$ | $2.56\mathrm{e}{-3}$ |  92
17 | $3.34\mathrm{e}{-3}$ | $4.57\mathrm{e}{-3}$ | $3.05\mathrm{e}{-3}$ | $1.00\mathrm{e}{-3}$ | $3.43\mathrm{e}{-3}$ | $1.62\mathrm{e}{-3}$ |  86
Figure 7: Sine-Gordon equation, $n=1001,k=13$: relative errors in the state
vector (left) and the Hamiltonian function (right).

####  5.3.3 Schrödinger equation

The third test is a 1D nonlinear Schrödinger equation with periodic boundary
conditions

$\begin{array}[]{rclcl}\mathrm{i}\,z_{t}&=&-z_{\xi\xi}-\varepsilon|z|^{2}z&\quad\text{
in }&(0,T)\times(-\frac{L}{2},\frac{L}{2}),\\\
z(t,-\frac{L}{2})&=&z(t,\frac{L}{2})&\quad\text{ in }&(0,T),\\\
z(0,\xi)&=&z_{0}(\xi)&\quad\text{ in }&(-\frac{L}{2},\frac{L}{2}),\end{array}$
(34)

where $z$ is the unknown complex valued wave function, $\varepsilon>0$, and
the initial condition has the form

$z_{0}(\xi)=\frac{\sqrt{2}}{\cosh(\xi-\xi_{0})}\exp\Bigl{(}\mathrm{i}\frac{c\,(\xi-\xi_{0})}{2}\Bigr{)}$

with the wave speed $c>0$ and $\xi_{0}\geq 0$. By introducing the real and
imaginary part of the wave function $z=z_{1}+\mathrm{i}\,z_{2}$, we can turn
the Schrödinger equation (34) into an infinite-dimensional Hamiltonian system

$\displaystyle\begin{split}(z_{1})_{t}&=-(z_{2})_{\xi\xi}-\varepsilon(z_{1}^{2}+z_{2}^{2})z_{2},\\\
(z_{2})_{t}&=\enskip\,\,(z_{1})_{\xi\xi}+\varepsilon(z_{1}^{2}+z_{2}^{2})z_{1},\end{split}$
(35)

with the associated Hamiltonian function

$\mathcal{H}(z_{1},z_{2})=\int_{0}^{L}\frac{1}{2}(z_{1})_{\xi}^{2}+\frac{1}{2}(z_{2})_{\xi}^{2}-\frac{\varepsilon}{4}(z_{1}^{2}+z_{2}^{2})^{2}\,{\rm
d}\xi.$

A spatial discretization of (35) by using the finite difference method on a
uniform spatial grid $\xi_{j}=jh_{\xi}$ with $h_{\xi}=L/n$ leads to the
nonlinear Hamiltonian system

$\displaystyle\begin{split}\begin{bmatrix}\dot{q}\\\
\dot{p}\end{bmatrix}&=J_{2n}\begin{bmatrix}-D_{\xi\xi}-G(q,p)&0\\\
0&-D_{\xi\xi}-G(q,p)\end{bmatrix}\begin{bmatrix}q\\\
p\end{bmatrix},\end{split}$ (36)

where $q=[q_{1},\ldots,q_{n}]^{T}$ with $q_{j}(t)=z_{1}(t,\xi_{j})$,
$p=[p_{1},\ldots,p_{n}]^{T}$ with $p_{j}(t)=z_{2}(t,\xi_{j})$, $D_{\xi\xi}$ is
the three-point finite difference matrix, and  G(q,p)=εdiag(q_1^2 + p_1^2,…,
q_n^2 + p_n^2).  The Hamiltonian function of (36) has the form

$H(q,p)=\sum_{i=1}^{n}\biggl{(}\frac{q_{i}^{2}-q_{i}q_{i-1}}{h_{\xi}^{2}}+\frac{p_{i}^{2}-p_{i}p_{i-1}}{h_{\xi}^{2}}-\frac{\varepsilon}{4}(p_{i}^{2}+q_{i}^{2})^{2}\biggr{)}.$

Numerical setting for simulation is similar to that in [16]:
$\varepsilon=1.0932$, $L=2\pi/0.11$, $c=1$, $\xi_{0}=0$, $h_{\xi}=0.2231$,
$n=1024$, $T=30$, and $h_{t}=0.01$. From the FOM solution at 3000 time
instances, we employ $s=750$ of them for computing the ROMs of dimension $2k$
with $k=95,100,105,110$. We run the optimization schemes with ${\tt iter}=500$
iterations. Similarly, we report the model reduction results in Table 6 and
Figure 8.

Table 6: Schrödinger equation, $n=1024$: model reduction results. $k$ | CotLift | CayleyC | SRE | a.a.f.
---|---|---|---|---
$\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$ | $\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$ | $\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$
95 | $7.62\mathrm{e}{-2}$ | $1.22\mathrm{e}{-1}$ | $8.06\mathrm{e}{-2}$ | $1.35\mathrm{e}{-1}$ | $5.88\mathrm{e}{-2}$ | $7.62\mathrm{e}{-2}$ | 1.74
100 | $6.77\mathrm{e}{-2}$ | $9.67\mathrm{e}{-2}$ | $1.98\mathrm{e}{-2}$ | $8.46\mathrm{e}{-3}$ | $2.94\mathrm{e}{-2}$ | $1.85\mathrm{e}{-2}$ | 1.52
105 | $4.48\mathrm{e}{-2}$ | $4.30\mathrm{e}{-2}$ | $4.71\mathrm{e}{-2}$ | $4.86\mathrm{e}{-2}$ | $2.14\mathrm{e}{-2}$ | $9.77\mathrm{e}{-3}$ | 1.36
110 | $3.60\mathrm{e}{-2}$ | $2.78\mathrm{e}{-2}$ | $3.95\mathrm{e}{-2}$ | $3.41\mathrm{e}{-2}$ | $3.42\mathrm{e}{-2}$ | $2.55\mathrm{e}{-2}$ | 1.35
Figure 8: Schrödinger equation, $n=1024$, $k=100$: relative errors in the
state vector (left) and the Hamiltonian function (right).

####  5.3.4 Vlasov equation

Finally, we consider a 1D Vlasov equation [57]

$z_{t}+v^{T}z_{\xi}-E(\xi)^{T}z_{v}=0,$ (37)

where $z=z(t,\xi,v)$ is the particle density function, $E(\xi)=-\tfrac{\rm
d}{{\rm d}\xi}\phi(\xi)$ is the electric field with a potential function
$\phi(\xi)$. In the particle-in-cell method [58], one assumes that there are
totally $n$ particles whose position and velocity at time $t$ are denoted,
respectively, by $q_{j}(t)$ and $p_{j}(t)$, $j=1,\ldots,n$. Approximating $z$
by

$\tilde{z}(t,\xi,v)=\sum_{j=1}^{n}w_{j}\delta(\xi-
q_{j}(t))\delta(v-p_{j}(t)),$

where $\delta$ is the Dirac function and $w_{j}$ is the $j$-th weight, and
requiring that the zeroth and first-order moments in $\xi$ and $v$ of (37) are
satisfied for this approximation, we obtain the equations

$\dot{q}_{j}=p_{j},\quad\dot{p}_{j}=-E(q_{j}),\quad j=1,\dots,n.$ (38)

Introducing $q=[q_{1}^{T},\ldots,q_{n}^{T}]^{T}$ and
$p=[p_{1}^{T},\ldots,p_{n}^{T}]^{T}$, equations (38) can shortly be written as
the Hamiltonian system (23) with $x=[q^{T},\,p^{T}]^{T}$ and the Hamiltonian
function

$H(q,p)=\sum_{j=1}^{n}\biggl{(}\frac{1}{2}p_{i}^{T}p_{i}-\phi(q_{i})\biggr{)}.$

We use the same numerical setting as in [57]: $t\in(0,0.2]$, $\xi\in[0,1]$,
the periodic boundary condition and the initial condition

$z(0,\xi,v)=\frac{1+\varepsilon\cos(2\pi\xi)}{\sqrt{2\pi}(a+1)}\left(\exp\Big{(}\\!-\frac{v^{2}}{2}\Big{)}+\frac{a}{\sigma}\,\exp\Big{(}\\!-\frac{(v-v_{0})^{2}}{2\sigma^{2}}\Big{)}\right)$
(39)

with $\varepsilon=0.3$, $a=0.3$, $v_{0}=4$, and $\sigma=1$. We assume that
there are $n\\!=\\!1000$ particles in the system. The time interval is equally
divided into subintervals of length $h_{t}=0.0001$. The initial values $q(0)$
and $p(0)$ are generated as random variables from the distribution111Tristan
Ursell (2020). Generate Random Numbers from a 2D Discrete Distribution, MATLAB
Central File Exchange. Retrieved April 29, 2020.
(https://www.mathworks.com/matlabcentral/fileexchange/35797-generate-random-
numbers-from-a-2d-discrete-distribution) (39). The electric field is chosen as
$E(\xi)=3\cos(4\pi\xi)$. Among $2000$ time instances of the FOM solution,
$s=400$ snapshots are used for constructing the ROMs of dimension $2k$ with
$k=6,8,10,12$. In the optimization schemes, we choose ${\tt iter}=1000$. Table
7 and Figure 9 show the obtained model reduction results.

Table 7: Vlasov equation, $n=1001$: model reduction results. $k$ | CotLift | CayleyC | SRE | SVD-like | a.a.f.
---|---|---|---|---|---
$\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$ | $\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$ | $\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$ | $\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$
$6$ | $1.76\mathrm{e}{-3}$ | $1.32\mathrm{e}{-\,\ 6}$ | $1.65\mathrm{e}{-3}$ | $2.40\mathrm{e}{-\ \,7}$ | $1.10\mathrm{e}{-3}$ | $6.67\mathrm{e}{-\,\ 7}$ | $2.94\mathrm{e}{-4}$ | $4.37\mathrm{e}{-\ \,9}$ | $13$
8 | $5.04\mathrm{e}{-4}$ | $6.53\mathrm{e}{-\,\ 9}$ | $4.75\mathrm{e}{-4}$ | $1.94\mathrm{e}{-\,\ 8}$ | $5.11\mathrm{e}{-4}$ | $5.11\mathrm{e}{-\,\ 9}$ | $7.61\mathrm{e}{-5}$ | $7.56\mathrm{e}{-11}$ | $12$
10 | $2.23\mathrm{e}{-4}$ | $4.68\mathrm{e}{-10}$ | $2.22\mathrm{e}{-4}$ | $6.73\mathrm{e}{-10}$ | $1.19\mathrm{e}{-4}$ | $1.08\mathrm{e}{-\ \,9}$ | $5.43\mathrm{e}{-4}$ | $6.00\mathrm{e}{-\ \,9}$ | $11$
12 | $6.32\mathrm{e}{-5}$ | $5.33\mathrm{e}{-11}$ | $6.30\mathrm{e}{-5}$ | $5.26\mathrm{e}{-11}$ | $6.31\mathrm{e}{-5}$ | $5.20\mathrm{e}{-11}$ | $2.05\mathrm{e}{-2}$ | $1.46\mathrm{e}{-\,\ 3}$ | $10$
Figure 9: Vlasov equation, $n=1001$, $k=10$: relative errors in the state
vector (left) and the Hamiltonian function (right).

####  5.3.5 Discussion

First, for different test models, the relative errors in the state
($\mathrm{RE}_{x}$) generally decrease when the ROMs computed by using the
optimization schemes are getting larger, see Figure 10 for a comparison.
However, the error in the energy ($\mathrm{RE}_{H}$) sometimes does not behave
the same. A difference in the linear and nonlinear cases can be observed.
Recalling that theoretically, if $x_{0}\in\mbox{im}(U)$, then the error in the
Hamiltonian function should be zero. However, in practice, due to different
errors, most considerably caused by the time integration solver, the computed
error in the Hamiltonian function is nonzero and apparently independent of the
dimension of the ROMs, see, e.g., Table 4 for the linear wave equation.
Meanwhile, for nonlinear models, this error additionally suffers from the
approximation of the nonlinear term. In view of [15, Theorem 5.1], we can
generally state that the more DEIM modes in approximating the nonlinear term
we use, which is the number of columns of $V$ in (26) and (27), the smaller
error in the energy we obtain. In our implementation, this number is nailed to
about $2.5k$ and therefore explains the monotone decrease of the error versus
reduced order in the case of nonlinear Hamiltonian systems.

Figure 10: Relative errors $\mathrm{RE}_{x}$ in the state vector for the ROMs
of different dimensions.

Second, as the cotangent lift method provides only a suboptimal solution to
the minimization problem (25), the optimization-based methods almost always
improve the reduction results by delivering slightly smaller errors in both
state vector and energy function. In the optimization aspect, this fact can be
viewed as a result of the decrease in the value of the cost function $f$.
Figure 11 reports the history of the values of the cost function for the wave
and sine-Gordon equations. One can observe that $f(X_{i})$ is decreasing in
Algorithm 4 starting from a suboptimal solution generated by the cotangent
lift method, leading to a reduction in errors. However, the correlation
between the amounts of decrease in the values of the cost function and the
resulting model reduction error $\mathrm{RE}_{x}$ is unclear. Furthermore,
note that the amount of improvement is in general rather modest. On the one
hand, the starting point is suboptimal which means it is already quite a good
approximation. On the other hand, the optimization method here is of first
order, i.e., only first-order information of the cost function is employed.
This class of methods is known to be slow when approaching the limit.
Moreover, the final model reduction error depends also on other factors, e.g.,
the approximation of the nonlinear term, which also attenuates the
improvement.

Figure 11: Values of the cost function for the wave and sine-Gordon equations
provided by the CayleyC and SRE optimization schemes.

Third, unlike linear Hamiltonian systems, for which the simulation time of
ROMs linearly depends on the reduced order, that of nonlinear systems is
difficult to predict as it also strongly depends on the structure of the
nonlinear term. Moreover, if implicit integration schemes, like the
Crank–Nicolson method in our case, are used for simulation, then the
convergence properties of the employed nonlinear solver may significantly
impact the computational time. In view of this, we can explain the low
accelerating factor for the Schrödinger equation compared to other models.
Indeed, the nonlinear term $\big{[}(G(q,p)q)^{T}\;(G(q,p)p)^{T}\big{]}^{T}$ in
the FOM (36) allows a vectorization in computing which results in quite a fast
simulation of the FOM even in large dimension while that of the ROM does not
have such a structure any more. As a consequence, even though the PSD-DEIM
method is exploited, i.e., only a moderate number ($m$ instead of $2n$) of
components of the nonlinear term are evaluated, the simulation of the ROM is
still quite time consuming, especially in designing a selection strategy via a
selector matrix $P$. Note also that for the Schrödinger model, the ratio of
the order reduction, $k/n$, is largest, which also contributes to the observed
fact.

Fourth, the SVD-like approach does not work reliably for the sine-Gordon and
Schrödinger equations with the reported setting which makes us hesitate to
present the result. Here, we used the available code from [17] and it is
claimed there that the basis can deviate from being symplectic in some cases.
This situation is unfortunately experienced in our experiments for these
models, and also for the wave equation when $k=40$. It is most probably the
reason why the simulation of the computed ROMs suffers from ill-conditioning
and results in large errors. Nevertheless, the SVD-like reduction method
delivers ROMs with the smallest errors in some situations for the Vlasov
equation. Further investigations are required to get a full picture of this
approach.

Finally, presenting Table 8 with the relative errors in the state vector and
energy for the reduced-order Vlasov models with different approximated
nonlinear terms given in (26) and (27), we would like to convey the message
that the structure preservation must be encompassed by an adequate
approximation to have expectedly satisfactory result. The ROM with (27) is
obviously Hamiltonian but a considerable error might have been caused during
the approximation of the nonlinear term. As a consequence, the overall quality
of the ROM obtained by the structure-preserving PSD-DEIM, including the energy
preservation, is apparently not as good as that computed by the PSD-DEIM,
which is not exactly structure-preserving.

Table 8: Vlasov equation, $n=1001$, $k=6$: comparison of the PSD-DEIM and structure-preserving PSD-DEIM methods. PSD-DEIM | CotLift | CayleyC | SRE | SVD-like
---|---|---|---|---
$\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$ | $\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$ | $\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$ | $\mathrm{RE}_{x}$ | $\mathrm{RE}_{H}$
(26) | $1.76\mathrm{e}{-3}$ | $1.32\mathrm{e}{-6}$ | $1.65\mathrm{e}{-3}$ | $2.40\mathrm{e}{-7}$ | $1.10\mathrm{e}{-3}$ | $6.67\mathrm{e}{-7}$ | $2.94\mathrm{e}{-4}$ | $4.37\mathrm{e}{-9}$
(27) | $1.08\mathrm{e}{-2}$ | $6.35\mathrm{e}{-5}$ | $2.62\mathrm{e}{-3}$ | $6.34\mathrm{e}{-5}$ | $3.08\mathrm{e}{-3}$ | $6.35\mathrm{e}{-5}$ | $1.22\mathrm{e}{-2}$ | $6.34\mathrm{e}{-5}$

##  6 Conclusion

We have proposed a new retraction on the symplectic Stiefel manifold which is
based on an SR decomposition. Its domain contains the unit ball and thus
results in the global convergence of the corresponding Riemannian gradient-
based optimization method.

Various applications and examples have also been presented for validating and
comparing the optimization schemes derived by combining different metrics and
retractions. Numerical results showed that, depending on the problem and the
setting, choice for metric and retraction must be taken with care to get
better result. Especially, schemes that use the SR retraction maintains the
symplecticity constraint the best. Running the same number of iteration, when
$k$ is considerably smaller than $n$, there is not much difference in the time
consumed by the schemes based on either the canonical-like metric or the
Euclidean metric. However, if $k$ is approaching $n$, the schemes with
canonical-like metric are faster. The numerical results also suggest that in
this case, the Cayley retraction is favorably combined with the canonical-like
metric as this combination tends to be faster than the others.

In addition, as an accumulation point of the optimization iterates can be a
saddle point, deep investigation on the structure of the set of critical
points of the cost function is always helpful in assuring to compute a
minimizer. Finally, the fact that in many cases, especially in the symplectic
model reduction problem, the presented optimization methods can only slightly
improve the established model reduction methods, urges an investigation for
faster optimization methods such as conjugate gradient method and Newton
method.

## Acknowledgment

Part of this work was done when Bin Gao and Nguyen Thanh Son were with ICTEAM
Institute, UCLouvain and supported by the Fonds de la Recherche Scientifique –
FNRS and the Fonds Wetenschappelijk Onderzoek – Vlaanderen under EOS Project
no 30468160. The authors also would like to thank P.-A. Absil for helpful
discussions.

## References

  * [1] B. Gao, N. T. Son, P.-A. Absil, T. Stykel, Riemannian optimization on the symplectic Stiefel manifold, SIAM J. Optim. 31 (2) (2021) 1546–1575.  doi:10.1137/20M1348522.
  * [2] T. Bendokat, R. Zimmermann, The real symplectic Stiefel and Grassmann manifolds: metrics, geodesics and applications, Preprint arXiv:2108.12447 [math.DG] (2021).  arXiv:2108.12447.
  * [3] J. Nocedal, S. J. Wright, Numerical Optimization, 2nd Edition, Springer-Verlag, New York, 2006.  doi:10.1007/978-0-387-40065-5.
  * [4] B. Gao, N. T. Son, P.-A. Absil, T. Stykel, Geometry of the symplectic Stiefel manifold endowed with the Euclidean metric, in: F. Nielsen, F. Barbaresco (Eds.), Geometric Science of Information: GSI 2021, Vol. 12829 of Lecture Notes in Computer Science, Springer Nature, Cham, Switzerland, 2021, pp. 789–796.  doi:10.1007/978-3-030-80209-7_85.
  * [5] J. Della-Dora, Numerical linear algorithms and group theory, Linear Algebra Appl. 10 (3) (1975) 267–283.  doi:10.1016/0024-3795(75)90074-9.
  * [6] A. Bunse-Gerstner, Matrix factorizations for symplectic QR-like methods, Linear Algebra Appl. 83 (1986) 49–77.  doi:10.1016/0024-3795(86)90265-X.
  * [7] A. Salam, On theoretical and numerical aspects of symplectic Gram–Schmidt-like algorithms, Numer. Algor. 39 (4) (2005) 437–462.  doi:10.1007/s11075-005-0963-2.
  * [8] A. J. Draft, F. Neri, G. Rangarajan, D. R. Douglas, L. M. Healy, R. D. Ryne, Lie algebraic treatment of linear and nonlinear beam dynamics, Ann. Rev. Nuclear Particle Sci. 38 (1988) 455–496.  doi:10.1146/annurev.ns.38.120188.002323.
  * [9] W. Harris, The average eye, Opthalmic Physiol. Optics 24 (6) (2004) 580–585.  doi:10.1111/j.1475-1313.2004.00239.x.
  * [10] R.-B. Wu, R. Chakrabarti, H. Rabitz, Optimal control theory for continuous-variable quantum gates, Phys. Rev. A 77 (2008) 052303.  doi:10.1103/PhysRevA.77.052303.
  * [11] P. Lancaster, Stability of linear gyroscopic systems: A review, Linear Algebra Appl. 439 (2013) 686–706.  doi:10.1016/j.laa.2012.12.026.
  * [12] R. Bhatia, T. Jain, On the symplectic eigenvalues of positive definite matrices, J. Math. Phys. 56 (11) (2015) 112201.  doi:10.1063/1.4935852.
  * [13] N. T. Son, P.-A. Absil, B. Gao, T. Stykel, Computing symplectic eigenpairs of symmetric positive-definite matrices via trace minimization and Riemannian optimization., SIAM J. Matrix Anal. Appl. 42 (4) (2021) 1732–1757.  doi:10.1137/21M1390621.
  * [14] N. T. Son, T. Stykel, Symplectic eigenvalues of positive-semidefinite matrices and the trace minimization theorem, Electron. J. Linear Algebra 38 (2022) 607–616.  doi:10.13001/ela.2022.7351.
  * [15] L. Peng, K. Mohseni, Symplectic model reduction of Hamiltonian systems, SIAM J. Sci. Comput. 38 (1) (2016) A1–A27.  doi:10.1137/140978922.
  * [16] B. Afkham, J. Hesthaven, Structure preserving model of parametric Hamiltonian systems, SIAM J. Sci. Comput. 39 (6) (2017) A2616–A2644.
  * [17] P. Buchfink, A. Bhatt, B. Haasdonk, Symplectic model order reduction with non-orthonormal bases, Math. Comput. Appl. 24 (2) (2019) Article ID 43.  doi:10.3390/mca24020043.
  * [18] P. Buchfink, S. Glas, B. Haasdonk, Optimal bases for symplectic model order reduction of canonizable linear Hamiltonian systems, IFAC-PapersOnLine 55 (20) (2022) 463–468.  doi:10.1016/j.ifacol.2022.09.138.
  * [19] T. Bendokat, R. Zimmermann, Geometric optimization for structure-preserving model reduction of Hamiltonian systems, IFAC-PapersOnLine 55 (20) (2022) 457–462.  doi:10.1016/j.ifacol.2022.09.137.
  * [20] S. Chaturantabut, D. Sorensen, Nonlinear model reduction via discrete empirical interpolation, SIAM J. Sci. Comput. 32 (5) (2010) 2737–2764.  doi:10.1137/090766498.
  * [21] S. Chaturantabut, C. Beattie, S. Gugercin, Structure-preserving model reduction for nonlinear port-Hamiltonian systems, SIAM J. Sci. Comput. 38 (5) (2016) B837–B865.  doi:10.1137/15M1055085.
  * [22] P.-A. Absil, R. Mahony, R. Sepulchre, Riemannian geometry of Grassmann manifolds with a view on algorithmic computations, Acta Appl. Math 8 (2004) 199–220.  doi:10.1023/B:ACAP.0000013855.14971.91.
  * [23] G. Golub, C. V. Loan, Matrix Computations. 4th ed, The Johns Hopkins University Press, Baltimore, London, 2013.
  * [24] P.-A. Absil, R. Mahony, R. Sepulchre, Optimization Algorithms on Matrix Manifolds, Princeton University Press, Princeton, NJ, 2008.
  * [25] H. Oviedo, R. Herrera, An efficient retraction mapping for the symplectic Stiefel manifold, Preprint 8478 (2021).
URL https://optimization-online.org/2021/07/8478

  * [26] P. Benner, R. Byers, H. Faßbender, V. Mehrmann, D. Watkins, Cholesky-like factorizations of skew-symmetric matrices, Electron. Trans. Numer. Anal. 11 (2000) 85–93.
  * [27] H. Faßbender, R. Rozložník, On the conditioning of factors in the SR decomposition, Linear Algebra Appl. 505 (2016) 224–244.  doi:0.1016/j.laa.2016.04.038.
  * [28] A. Salam, E. Al-Aidarous, A. E. Farouka, Optimal symplectic Householder transformations for SR decomposition, Linear Algebra Appl. 429 (5-6) (2008) 1334–1353.  doi:10.1016/j.laa.2008.02.029.
  * [29] V. Mehrmann, Der SR-Algorithmus zur Berechnung der Eigenwerte einer Matrix, Diplomarbeit, Universität Bielefeld (1979).
  * [30] J. Barzilai, J. M. Borwein, Two-point step size gradient methods, IMA J. Numer. Anal. 8 (1) (1988) 141–148.  doi:10.1093/imanum/8.1.141.
  * [31] R.-B. Wu, R. Chakrabarti, H. Rabitz, Critical landscape topology for optimization on the symplectic group, J. Optim. Theory Appl. 145 (2010) 387–406.  doi:10.1007/s10957-009-9641-1.
  * [32] S. Fiori, A Riemannian steepest descent approach over the inhomogeneous symplectic group: application to the averaging of linear optical systems, Appl. Math. Comput. 283 (2016) 251–264.  doi:10.1016/j.amc.2016.02.018.
  * [33] F. Dopico, C. Johnson, Parametrization of the matrix symplectic group and applications, SIAM J. Matrix Anal. Appl. 31 (2) (2009) 650–673.  doi:10.1137/060678221.
  * [34] J. Williamson, On the algebraic problem concerning the normal forms of linear dynamical systems, Am. J. Math. 58 (1) (1936) 141–163.
  * [35] P. Amodio, On the computation of few eigenvalues of positive definite Hamiltonian matrices, Future Generation Computer Systems 22 (4) (2006) 403–411.  doi:10.1016/j.future.2004.11.027.
  * [36] P. Lancaster, R. Rodman, Canonical forms for Hermitian matrix pairs under strict equivalence and congruence, SIAM Rev. 47 (2006) 407–443.  doi:10.1137/S003614450444556X.
  * [37] T. Jain, H. K. Mishra, Derivatives of symplectic eigenvalues and a Lidskii type theorem, Canad. J. Math. 74 (2022) 457–485.  doi:10.4153/S0008414X2000084X.
  * [38] T. Hiroshima, Additivity and multiplicativity properties of some Gaussian channels for Gaussian inputs, Phys. Rev. A 73 (2006) 012330.  doi:10.1103/PhysRevA.73.012330.
  * [39] J. Kovač-Striko, K. Veselić, Trace minimization and definiteness of symmetric pencils, Linear Algebra Appl. 216 (1995) 139–158.  doi:10.1016/0024-3795(93)00126-K.
  * [40] I. Nakić, K. Veselić, Wielandt and Ky-Fan theorem for matrix pairs, Linear Algebra Appl. 369 (2003) 77–93.  doi:10.1016/S0024-3795(02)00733-4.
  * [41] X. Liang, R. C. Li, Z. Bai, Trace minimization principles for positive semi-definite pencils, Linear Algebra Appl. 438 (2013) 3085–3106.  doi:10.1016/j.laa.2012.12.003.
  * [42] X. Liang, L. Wang, L.-H. Zhang, R.-C. Li, On generalizing trace minimization principles, Linear Algebra Appl. 656 (2023) 483–509.  doi:10.1016/j.laa.2022.10.012.
  * [43] P. Benner, H. Fassbender, M. Stoll., Solving large-scale quadratic eigenvalue problems with Hamiltonian eigenstructure using a structure-preserving Krylov subspace method, Electron. Trans. Numer. Anal. 29 (2008) 212–229.
  * [44] F. J. Buijnsters, A. Fasolino, M. I. Katsnelson, Zero modes in magnetic systems: General theory and an efficient computational scheme, Phys. Rev. B 89 (2014) 174433.  doi:10.1103/PhysRevB.89.174433.
  * [45] M. Krbek, T. Tyc, J. Vlach, Inequalities for quantum marginal problems with continuous variables, J. Math. Phys. 55 (2014) 062201.  doi:10.1063/1.4880198.
  * [46] I. L. Egusquiza, A. Parra-Rodziguez, Algebraic canonical quantization of lumped superconducting networks, Phys. Rev. B 106 (2022) 024510.  doi:10.1103/PhysRevB.106.024510.
  * [47] S. Wei, I. Kao, Vibration analysis of wire and frequency response in the modern wiresaw manufacturing process, J. Sound Vib. 231 (5) (2000) 2383–1395.  doi:10.1006/jsvi.1999.247.
  * [48] J. Marsden, T. Ratiu, Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems, Vol. 17 of Texts in Applied Mathematics, Springer-Verlag, New York, NY, 1999.  doi:10.1007/978-0-387-21792-5.
  * [49] V. Duindam, A. Macchelli, S. Stramigioli, H. Bruyninckx, Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach, Springer, Berlin, Heidelberg, 2009.  doi:10.1007/978-3-642-03196-0.
  * [50] A. van der Schaft, D. Jeltsema, Port-Hamiltonian systems theory: An introductory overview, Foundations and Trends in Systems and Control 1 (2-3) (2014) 173–378.  doi:10.1561/2600000002.
  * [51] B. Afkham, A. Bhatt, B. Haasdonk, J. Hesthaven, Symplectic model-reduction with a weighted inner product, Preprint arXiv:1803.07799 [math.NA] (2018).  arXiv:1803.07799.
  * [52] C. Pagliantini, Dynamical reduced basis methods for Hamiltonian systems, Numer. Math. 148 (2021) 409–448.  doi:10.1007/s00211-021-01211-w.
  * [53] H. Sharma, Z. Wang, B. Kramer, Hamiltonian operator inference: Physics-preserving learning of reduced-order models for canonical Hamiltonian systems, Physica D: Nonlinear Phenomena (2021) 133122doi:10.1016/j.physd.2021.133122.
  * [54] J. S. Hesthaven, C. Pagliantini, N. Ripamonti, Rank-adaptive structure-preserving model order reduction of Hamiltonian systems, ESAIM: M2AN 56 (2022) 617–650.  doi:10.1051/m2an/2022013.
  * [55] H. Xu, An SVD-like matrix decomposition and its applications, Linear Algebra Appl. 368 (2003) 1–24.  doi:10.1016/S0024-3795(03)00370-7.
  * [56] H. Li, Y. Wang, Q. Sheng, An energy-preserving Crank–Nicolson Galerkin method for Hamiltonian partial differential equations, Numer. Methods Partial Differential Equations 32 (5) (2016) 1485–1504.  doi:10.1002/num.22062.
  * [57] T. M. Tyranowski, M. Kraus, Symplectic model reduction methods for the Vlasov equation, Preprint arXiv:1910.06026v1 [physics.comp-ph] (2019).  arXiv:1910.06026v1.
  * [58] G. Lapenta, Particle-based simulation of plasmas, in: G. Colonna, A. D’Angola (Eds.), Plasma Modeling: Methods and Applications, IOP Publishing, 2016, pp. 4–1 to 4–37.  doi:10.1088/978-0-7503-1200-4ch4.