additions

# Symplectic Autoencoders for Model Reduction of Hamiltonian Systems

Benedikt Brantner, Michael Kraus
<EMAIL_ADDRESS><EMAIL_ADDRESS>
Max-Planck-Institut für Plasmaphysik
Boltzmannstraße 2, 85748 Garching, Deutschland
and
Technische Universität München, Zentrum Mathematik
Boltzmannstraße 3, 85748 Garching, Deutschland

###### Abstract

Many applications, such as optimization, uncertainty quantification and
inverse problems, require repeatedly performing simulations of large-
dimensional physical systems for different choices of parameters. This can be
prohibitively expensive.

In order to save computational cost, one can construct surrogate models by
expressing the system in a low-dimensional basis, obtained from training data.
This is referred to as model reduction. Past investigations have shown that,
when performing model reduction of Hamiltonian systems, it is crucial to
preserve the symplectic structure associated with the system in order to
ensure long-term numerical stability.

Up to this point structure-preserving reductions have largely been limited to
linear transformations. We propose a new neural network architecture in the
spirit of autoencoders, which are established tools for dimension reduction
and feature extraction in data science, to obtain more general mappings.

In order to train the network, a non-standard gradient descent approach is
applied that leverages the differential-geometric structure emerging from the
network design. The new architecture is shown to significantly outperform
existing designs in accuracy.

###### Contents

  1. 1 Introduction
  2. 2 Hamiltonian Systems
    1. 2.1 Infinite-dimensional Hamiltonian systems
  3. 3 Reduced Basis Methods
    1. 3.1 POD/PSD
      1. 3.1.1 Induced differential equations
      2. 3.1.2 Example of failing POD
    2. 3.2 Autoencoders
      1. 3.2.1 Feedforward Neural Networks
  4. 4 SympNets
    1. 4.1 Linear/Activation-type SympNets $\Psi_{LA}$
    2. 4.2 Gradient-type SympNets $\Psi_{G}$
  5. 5 Proposed Architecture
  6. 6 Manifold Optimization for Neural Networks
    1. 6.1 Optimization on Riemannian manifolds and retractions
    2. 6.2 Generalizing Adam
  7. 7 Experiments
    1. 7.1 Initial conditions
  8. 8 Summary
  9. A Riemannian gradients for the Stiefel and the Symplectic Stiefel manifold

##  1 Introduction

Real-world physical systems are often described by parametrized partial
differential equations (PDEs). In order to solve these, they are ususally
approximated by very high-dimensional ordinary differential equations (ODEs).
These discretized systems are referred to as full order models (FOM).
Applications such as optimization, uncertainty quantification and inverse
problems require performing repeated simulations of the FOM which can quickly
become prohibitively expensive. It has become common practice to construct
reduced order models (ROM) in a data-driven way to alleviate this
computational bottleneck. Ideally, the ROM is designed to capture the FOM
dynamics as accurately as possible (see [2]).

ROMs always consist of two mappings: a reduction $\mathcal{P}$ from the high-
to the low-dimensional space, and a reconstruction $\mathcal{R}$ from the low-
to the high-dimensional space. Once we have $\mathcal{P}$ and $\mathcal{R}$,
we can naturally obtain an ODE on the low-dimensional space by means of
Galerkin projection. This is ideally much cheaper to solve than the ODE
describing the FOM.

The most-common of these data-driven model reduction techniques is “proper
orthogonal decomposition” (POD111This is also known as principal component
analysis (PCA)., see [25]), in which both $\mathcal{P}$ and $\mathcal{R}$ are
orthonormal matrices. POD uses singular value decomposition (SVD) to identify
dominant global modes in the data matrix on which it is performed.

Despite its ubiquity POD has two limitations: First, it is strictly linear.
This means that problems with slowly-decaying singular values (i.e. advection-
dominated problems) require a large amount of modes (see [26, 7]).

For this reason it is becoming increasingly more common to construct
$\mathcal{P}$ and $\mathcal{R}$ using non-linear mappings. If neural networks
are used, this amounts to using autoencoders (see [26, 11]). Autoencoders are
neural networks that represent high-fidelity data using a low-dimensional
feature space. For reduced order modelling this feature space is the low-
dimensional space on which the dynamics evolve.

The second limitation of POD (this is also true for practically all
applications of autoencoders) is that it is oblivious to the physical
structure of the problem it is applied to.

This neglect can have grave consequences for the accuracy of the dynamics of
the reduced system. An extreme example of this is shown in section 3.1.2. Of
particular importance among these physical properties is symplecticity (see
e.g. [3]), which is in most cases equivalent to the system being Hamiltonian.
If the high-fidelity dynamics is Hamiltonian, an alternative to POD exists
that takes the Hamiltonian structure of the system into account: Proper
Symplectic Decomposition (PSD, [30, 33]). PSD, however, is affected by the
other problem that POD suffers from: It is strictly linear.

In this work we attempt to combine the advantages of autoencoders and PSD by
constructing non-linear mappings that also take symplecticity into account
(i.e. the Hamiltonian structure of the system). This is achieved by composing
an existing neural network architecture, called symplectic neural networks
(SympNets, [20]), with PSD-like matrices. These PSD-like matrices form a
manifold. Training this new neural network will then require techniques from
differential geometry and manifold optimization.

The remainder of this paper is structured as follows: sections 2, 3 and 4
discuss basic notions of Hamiltonian dynamics, reduced basis methods and
SympNets respectively. Section 5 gives a description of the proposed
architecture. Section 6 discusses details of the implementation and section 7
shows results of experiments performed with the new neural networks.

##  2 Hamiltonian Systems

The following offers a brief review of Hamiltonian dynamics. For more details
consult e.g. [6, 3].

A canonical Hamiltonian system on $\mathbb{R}^{2N}$ is defined through a
Hamiltonian function $H:\mathbb{R}^{2N}\to\mathbb{R}$. The coordinates on
$\mathbb{R}^{2N}$ are denoted by $(q^{1},\ldots,q^{N},p_{1},\ldots,p_{N})$.

From $H$ we obtain a Hamiltonian vector field by first taking the gradient of
the function and then applying the matrix222In theory this can be an arbitrary
skew-symmetric non-degenerate matrix. But since one can always find a linear
transformation from such a matrix onto $\mathbb{J}_{2N}$ one typically just
chooses the latter. Further we will refer to $\mathbb{J}$ as the symplectic
potential. It is sometimes also called Poisson tensor. $\mathbb{J}_{2N}$:

$\mathbb{J}_{2N}=\begin{pmatrix}\mathbb{O}_{N}&\mathbb{I}_{N}\\\
-\mathbb{I}_{N}&\mathbb{O}_{N}\end{pmatrix}.$ (1)

The Hamiltonian vector field hence takes the following form:

$V_{H}(q,p)=\mathbb{J}_{2N}\nabla{}H(q,p),$ (2)

or equivalently:

$\dot{q}=\frac{\partial{}H}{\partial{}p},\quad\dot{p}=-\frac{\partial{}H}{\partial{}q}.$
(3)

The flow $\varphi^{H}_{t}(\cdot)$ of this differential equation has two very
important properties: it preserves the Hamiltonian $H$ and the symplectic
potential $\mathbb{J}_{2N}$, i.e.

  1. 1.

$H\circ\varphi^{H}_{t}(z)=H(z)$ and

  2. 2.

$(\nabla\varphi^{H}_{t}(z))^{T}\mathbb{J}_{2N}\nabla\varphi^{H}_{t}(z)=\mathbb{J}_{2N}$.

There are both theoretical and empirical arguments showing that considering
this structure of the differential equations is extremely important when
treating Hamiltonian systems, because it greatly restricts the space of
possible solutions (see [16, 28] for details on how symplecticity restricts
the space of possible mappings). It has proven very important to consider this
structure when designing numerical schemes (see e.g. [17, 32, 27]).

###  2.1 Infinite-dimensional Hamiltonian systems

Analogous to Hamiltonian ODEs, we can also formulate infinite-dimensional
canonical Hamiltonian systems that are PDEs. Given a domain $\Omega$ we define
two (time-dependent) functions on it: $q:\mathcal{I}\times\Omega\to\mathbb{R}$
and $p:\mathcal{I}\times\Omega\to\mathbb{R}$. Here a Hamiltonian is a
functional that takes $q$ and $p$ as input:

$\mathcal{H}:C^{\infty}(\mathcal{I}\times\Omega)\times{}C^{\infty}(\mathcal{I}\times\Omega)\to\mathbb{R}.$
(4)

The resulting PDE takes a similar form to the case when we deal with
Hamiltonian ODEs:

$\partial_{t}{}q(t,\xi;\mu)=\frac{\delta{}\mathcal{H}}{\delta{}p},\quad\partial_{t}{}p(t,\xi;\mu)=-\frac{\delta{}\mathcal{H}}{\delta{}q}.$
(5)

An example of this is the “linear wave equation” (see [30, 9]) with
Hamiltonian

$\mathcal{H}(q,p;\mu):=\frac{1}{2}\int_{\Omega}\mu^{2}(\partial_{\xi}{}q(t,\xi;\mu))^{2}+p(t,\xi;\mu)^{2}d\xi.$
(6)

If we discretize the Hamiltonian directly with e.g. a structure-preserving
finite difference discretization we again end up with a finite-dimensional
Hamiltonian system (but typically in very high dimensions since it is the
discretization of a PDE). See section 7 for an example of this.

##  3 Reduced Basis Methods

The following gives a quick review of data-driven reduced basis methods
(RBMs). RBMs can alleviate the cost of solving parametric PDEs (PPDEs) by
leveraging machine learning techniques. The task here is the following:

Assume we are given a PPDE

$F(u(\mu);\mu)=0,\quad u(\mu)\in{}V,$ (7)

where $V$ is some infinite-dimensional Banach space. For solving this PDE we
have to first discretize $V$ in some form (finite element, finite volume,
finite difference, particle method, …); we denote this discretized space by
$V_{h}$ and refer to the induced differential equation on it as the full order
model (FOM). In order to achieve sensible simulation results we usually have
to choose the dimension of $V_{h}$ very large, which makes solving the
discretized PPDE very expensive. In practice, we typically also have to solve
the PPDE for many different parameter instances $\mu$, which can lead to
prohibitively high computational costs. Even though the solutions for
different parameter instances $\mu\in\mathbb{P}$ will be qualitatively
similar, we usually have to run the simulation from scratch every time we
change the parameters. Dealing with data and qualitatively similar solutions
(also backed up by the concept of the solution manifold333The solution
manifold is the space on which the solutions of the PDE evolve:
$\mathcal{M}:=\\{u(\mu):\mu\in\mathbb{P}\\}$.) this is an ideal setting to
apply machine learning techniques.

RBMs leverage the parametric structure of the system by constructing a low-
dimensional manifold and a corresponding embedding of the manifold in $V_{h}$
that approximate the dynamics in the high-dimensional space as accurately as
possible [11, 25]. This manifold and the embedding are determined in a data-
driven way based on a set of FOM solutions (this process is explained below).
In the following, the dimension of the big system is denoted by $2N$ and that
of the small system by $2n$. The small system is referred to as the reduced
order model (ROM).

The ODE in the high-dimensional space $V_{h}=\mathbb{R}^{2N}$, i.e. the FOM,
is denoted by:

$\begin{cases}\dot{\mathbf{u}}(t;\mu)=\mathbf{f}(t,\mathbf{u}(t;\mu);\mu),\qquad
t\in(0,T),\\\ \mathbf{u}(0;\mu)=\mathbf{u}_{0}(\mu),\end{cases}$ (8)

where $\mu$ is a set of parameters (potentially also including initial
conditions) and $\mathbf{u}(\cdot;\mu)$ is the solution of the ODE for
parameters $\mu$ that evolves on $V_{h}=\mathbb{R}^{2N}$.

The goal of reduced basis methods is the construction of a reduction
$\mathcal{P}:\mathbb{R}^{2N}\to\mathbb{R}^{2n}$ and a reconstruction (or lift)
$\mathcal{R}:\mathbb{R}^{2n}\to\mathbb{R}^{2N}$. The low-dimensional manifold
embedded in $V_{h}$, which was mentioned before, is
$\mathcal{R}(\mathbb{R}^{2n})$ in this case.

In practice, we want to construct $\mathcal{P}$ and $\mathcal{R}$ based on a
set of previous FOM simulations to make successive simulations cheaper. The
data coming from these previous simulations is represented in form of a
snapshot matrix

$M=[\mathbf{u}(t_{0};\mu_{i}),\ldots,\mathbf{u}(t_{T};\mu_{i})]_{i=1,\ldots,n_{\mathbb{P}}},$
(9)

where $t_{0}$ and $t_{f}$ are the initial and final time of the simulation.
With this the problem can be cast into the following form:

$\min_{\mathcal{P},\mathcal{R}}d(M,\mathcal{R}\circ\mathcal{P}(M)),$ (10)

where $d(\cdot,\cdot)$ is an appropriate metric (this can in principle be any
metric, but will be $L_{2}$ from now on). We refer to equation (10) as the
reconstruction error.

Having found the mappings $\mathcal{P}$ and $\mathcal{R}$ one has two choices
for a vector field on the reduced space:

$\displaystyle\mathbf{f}^{\mathrm{red}}(t,z;\mu)$
$\displaystyle=(\nabla_{\mathcal{R}(z)}\mathcal{P})(\mathbf{f}(t,\mathcal{R}(z);\mu))\text{
and }$ (11) $\displaystyle\mathbf{f}^{\mathrm{red}}(t,z;\mu)$
$\displaystyle=(\nabla_{z}\mathcal{R})^{+}(\mathbf{f}(t,\mathcal{R}(z);\mu)),$
(12)

where $(\cdot)^{+}$ indicates the pseudo-inverse. The second one of these
equations is a “manifold Galerkin projection” (see [26]) and gives the most
consistent vector field. We will therefore use the second equation to get
reduced dynamics from now on.

For Hamiltonian model reduction, both of these mappings have to be symplectic
in the sense of Section 2, i.e.,
$(\nabla_{(q,p)}\mathcal{P})\mathbb{J}_{2N}(\nabla_{(q,p)}\mathcal{P})^{T}=\mathbb{J}_{2n}$
and
$(\nabla_{(\tilde{q},\tilde{p})}\mathcal{R})^{T}\mathbb{J}_{2N}(\nabla_{(\tilde{q},\tilde{p})}\mathcal{R})=\mathbb{J}_{2n}$
where $(q,p)$ are the coordinates in the full space $\mathbb{R}^{2N}$ and
$(\tilde{q},\tilde{p})$ are the coordinates in the reduced space
$\mathbb{R}^{2n}$.

###  3.1 POD/PSD

The most widely-used model order reduction technique is Proper Orthogonal
Decomposition (POD, see [25]). This provides the ideal solution for the case
when $\mathcal{P}$ and $\mathcal{R}$ are assumed to be orthonormal
matrices444More precisely: $\mathcal{R}$ is assumed orthonormal and the
transpose of $\mathcal{P}$ is assumed to be orthonormal, i.e. $\mathcal{R}=A$
and $\mathcal{P}=B^{T}$ with $A^{T}A=B^{T}B=\mathbb{I}_{2n}\Longleftrightarrow
A,B\in{}St(2n,2N)$. and an $L_{2}$ norm is used. That a matrix is orthonormal
here means that $A^{T}A=\mathbb{I}_{2n}$. These orthonormal matrices form a
manifold, called the “Stiefel manifold” (see [17, 8]). The solution then just
requires the computation of an SVD for $M$: $M=U\Sigma{}V$. For the POD we set
$\mathcal{P}:\mathbb{R}^{N}\to\mathbb{R}^{n}$ and
$\mathcal{R}:\mathbb{R}^{n}\to\mathbb{R}^{N}$.

######  Proposition 1.

If $\mathcal{R}\equiv A$ is assumed to be a matrix and
$\mathcal{P}=\mathcal{R}^{T}$, then $A$ is necessarily orthonormal. It is
given by the first $n$ columns of the matrix $U$ in the SVD of $M$.

###### Proof.

Consider the SVD $M=U\Sigma{}V^{T}$. Then the reconstruction error reads
$||M-\mathcal{R}\circ\mathcal{P}(M)||_{2}=||U\Sigma{}V^{T}-AA^{T}U\Sigma{}V^{T}||_{2}=||(I-AA^{T})U\Sigma{}V^{T}||_{2}$.
Here $(I-AA^{T})$ is a matrix with rank $N-n$ (assuming that $A$ has full
rank). The minimum is achieved when $(I-AA^{T})$ removes the part of $U$ that
contributes the most to the error, i.e. the one for which $\Sigma$ has the
highest eigenvalues; this implies $A=[u_{1},\ldots,u_{n}]$, i.e. the first $n$
columns of $U$. ∎

POD is hence just SVD applied to a snapshot matrix. If the system, represented
by $M$, has additional structure (e.g. symplecticity), POD will be entirely
oblivious to this. Proper Symplectic Decomposition (PSD, [30]), on the other
hand, encodes the symplectic structure of the system into the reduced basis.
In this case the pseudo inverse of $A=\mathcal{R}$ is $A^{+}$, the symplectic
inverse: $A^{+}:=\mathbb{J}_{2n}A^{T}\mathbb{J}_{2N}^{T}$. This yields the
optimization problem:

$\min_{A\in{}Sp(2n,2N)}||M-AA^{+}M||_{2}$ (13)

where we look for the minimum on the symplectic Stiefel manifold
$Sp(2n,2N):=\\{A\in\mathbb{R}^{2N\times
2n}:A^{T}\mathbb{J}_{2N}A=\mathbb{J}_{2n}\\}$.

Here we cannot perform an SVD (or a symplectic analog) in the way it was used
in the proof of Theorem 1, which means we have to rely on techniques like
(Riemannian) optimization to find a global minimum on
$Sp(2n,2N)$555Technically speaking there is an equivalent of SVD where an
arbitrary matrix admits a decomposition $M=QDS^{-1}$ where $Q$ is orthogonal
and $S$ is symplectic (see [34]). This decomposition can, however, not be used
in the way it was done in the proof of Theorem 1 for SVD.. Some details of
optimizing on manifolds will be discussed in Section 6. Associated with
optimization on $Sp(2n,2N)$ is another difficulty: even though in principle
one could find a minimum of this optimization problem, this is numerically
highly unstable as the symplectic group and the symplectic Stiefel manifold
are, other than the orthonormal group and the Stiefel manifold, non-compact.

A potential workaround is enforcing an additional constraint such that $A$ is
a cotangent lift. The matrices we optimize over then take the following form:

$A=\begin{pmatrix}\Phi&\mathbb{O}_{N\times{}n}\\\
\mathbb{O}_{N\times{}n}&\Phi\end{pmatrix},\quad\Phi\in{St}(n,N).$ (14)

These matrices are easily verified to lie in $Sp(2n,2N)\cap{}St(2n,2N)$.
Computing matrices of the form as in Equation (14) again just requires an
application of SVD and this is referred to as proper symplectic decomposition
(PSD666Also not that PSD-like matrices shown in equation (14) are naturally
isomorphic to the Stiefel manifold $St(n,N)$., see [30]). So the difference
between POD and PSD is that for POD we find the minimum of the reconstruction
error on all orthonormal matrices, whereas for PSD we find it on matrices as
in Equation (14); i.e. symplecticity is encoded into the structure. Different
approaches to construct elements from $Sp(2n,2N)$, not just the cotangent
lift, to minimize the reconstruction error (see Equation (10)) are discussed
in [13]. In this reference, it is also observed that the cotangent lift works
remarkably well, also when compared to different, more involved techniques.
Henceforth we will use the terms cotangent lift and PSD interchangeably.

####  3.1.1 Induced differential equations

Both POD and PSD induce differential equations on the reduced system through
equation (12). In the case of POD, this is $\dot{z}=A^{T}V(Az)$ and for PSD it
is $\dot{z}=A^{+}V(Az)$, where $z$ are the reduced coordinates and $V$ is the
vector field. Here, $A^{+}:=\mathbb{J}_{2n}A^{T}\mathbb{J}_{2N}^{T}$ is the
symplectic inverse that was mentioned above.

If $V\equiv{}V_{H}$ is a Hamiltonian vector field, then the reduced system in
the PSD case will also be Hamiltonian:

$\dot{z}=A^{+}V_{H}(Az)=\mathbb{J}_{2n}A^{T}\mathbb{J}_{2N}^{T}\mathbb{J}_{2N}\nabla{}H(Az)=\mathbb{J}_{2n}A^{T}\nabla{}H(Az),$
(15)

where the Hamiltonian in the reduced system is now simply $z\mapsto{}H(Az)$.

####  3.1.2 Example of failing POD

An example where POD can go horribly wrong is the simple Hamiltonian system:

$H(\mathbf{q},\mathbf{p})=\sum_{i=1}^{2}\frac{p_{i}^{2}}{2}+\phi_{i}\frac{q_{i}^{2}}{2},$
(16)

with $\phi_{1}=0.05$, $\phi_{2}=\pi$ and
$((\mathbf{q}^{(0)})^{T},(\mathbf{p}^{(0)})^{T})=(q^{(0)}_{1},q^{(0)}_{2},p^{(0)}_{1},p_{2}^{(0)})=(0,0,1,3)$.
This is a system of two independent harmonic oscillators whose frequencies are
in an irrational relation to each other; the resulting orbit will thus wind
densely around an invariant 2-torus that is embedded in $\mathbb{R}^{4}$.

If one now takes the snapshot matrix

$M:=\begin{bmatrix}\mathbf{q}^{(0)}&\mathbf{q}^{(1)}&\cdots&\mathbf{q}^{(T)}\\\
\mathbf{p}^{(0)}&\mathbf{p}^{(1)}&\cdots&\mathbf{p}^{(T)}\end{bmatrix},$ (17)

and reduces the system to $\mathbb{R}^{2}$ (i.e. $N=2$ and $n=1$), then the
lifts/reductions have the following form777These are the matrices obtained in
the limit for $hT\to\infty$ and $h\to 0$, where $h$ is the stepsize.:

$A_{\mathrm{POD}}=\begin{pmatrix}1&0\\\ 0&0\\\ 0&0\\\
0&1\end{pmatrix}\quad\text{ and }\quad A_{\mathrm{PSD}}=\begin{pmatrix}1&0\\\
0&0\\\ 0&1\\\ 0&0\end{pmatrix}.$ (18)

This is because the $L_{2}$ distance between $(q_{1},0,0,p_{2})$ and
$(q_{1},q_{2},p_{1},p_{2})$ is the smallest possible one if only two modes are
considered. For the PSD, the ideal “reduced system” is $(q_{1},0,p_{1},0)$,
i.e. the first one of the two circles on the torus
$\mathbb{T}^{2}\simeq\mathbb{S}^{1}\times\mathbb{S}^{1}$.

To make this more clear, consider the full-order dynamics:

$\begin{pmatrix}\mathbf{q}^{(t)}\\\
\mathbf{p}^{(t)}\end{pmatrix}=\begin{pmatrix}q_{1}^{(t)}\\\ q_{2}^{(t)}\\\
p_{1}^{(t)}\\\
p_{2}^{(t)}\end{pmatrix}=\begin{pmatrix}p_{1}^{(0)}\sin(\sqrt{\phi_{1}}t)/\sqrt{\phi_{1}}\\\
p_{2}^{(0)}\sin(\sqrt{\phi_{2}}t)/\sqrt{\phi_{2}}\\\
p_{1}^{(0)}\cos(\sqrt{\phi_{1}}t)\\\
p_{2}^{(0)}\cos(\sqrt{\phi_{2}}t)\end{pmatrix}.$ (19)

The “dominant” of these dimensions is the first one, since its prefactor is
$p_{1}^{(0)}/\sqrt{\phi_{1}}=1/\sqrt{0.05}\approx 4.47$; the second biggest is
the forth dimension with prefactor $p_{2}^{(0)}=3$.

Since a 2-torus can be embedded in $\mathbb{R}^{3}$, a visualization of the
two reductions is possible. Figures 1(a) and 1(b) show
$A_{\mathrm{POD}}A_{\mathrm{POD}}^{T}M$ and
$A_{\mathrm{PSD}}A_{\mathrm{PSD}}^{+}M$ respectively (where $M$ is the data
matrix). The embedding for the two plots
$\mathbb{T}^{2}\hookrightarrow\mathbb{R}^{3}$ is
$(q_{1},q_{2},p_{1},p_{2})\mapsto(q_{1}+q_{2},p_{1},p_{2})$.

(a) POD: full order model (blue); reduced dynamics (orange).

(b) PSD: full order model (blue); reduced dynamics (orange).

In Figure 1(a), the POD picks the “wrong” two modes, i.e., they are not
corresponding $q$ and $p$ coordinates but form a Lissajous curve (orange) in
the $xz$-plane instead of a circle.

The difference between POD and PSD becomes even more apparent when looking at
the induced dynamics. The differential equation for the POD is:

$A_{\mathrm{POD}}^{T}V_{H}(A_{\mathrm{POD}}z)=\begin{pmatrix}1&0&0&0\\\
0&0&0&1\end{pmatrix}\begin{pmatrix}0\\\ p_{2}\\\ -\phi_{1}q_{1}\\\
0\end{pmatrix}=\begin{pmatrix}0\\\ 0\end{pmatrix},$ (20)

i.e. the induced trajectory is the constant one. On the other hand, for the
PSD, this is:

$A_{\mathrm{PSD}}^{+}V_{H}(A_{\mathrm{PSD}}z)=\begin{pmatrix}1&0&0&0\\\
0&0&1&0\end{pmatrix}\begin{pmatrix}p_{1}\\\ 0\\\ -\phi_{1}q_{1}\\\
0\end{pmatrix}=\begin{pmatrix}p_{1}\\\ -\phi_{1}q_{1}\end{pmatrix},$ (21)

i.e. the dynamics are correctly reduced to the first one of the two harmonic
oscillators. Another example of a failing POD for the Vlasov equation is shown
in [33].

###  3.2 Autoencoders

Whereas the previous section addressed that POD does not encode the
Hamiltonian dynamics in its reduced model, this section will deal with how to
remove the constraint that it is strictly linear.

The term “autoencoder” refers to a pair of neural networks, the encoder
$\mathcal{P}=\Psi^{\mathrm{enc}}$ and the decoder
$\mathcal{R}=\Psi^{\mathrm{dec}}$, that are “trained” to minimize a
reconstruction error $d(M,\Psi^{\mathrm{dec}}\circ\Psi^{\mathrm{enc}}(M))$ for
a given dataset $M$ (see [14, Chapter 14]). If $\Psi^{\mathrm{enc}}$ and
$\Psi^{\mathrm{dec}}$ are assumed to be linear (i.e. matrices), one again
arrives at POD, which may therefore be viewed as a simple, linear form of an
autoencoder. Originally, autoencoders were devised to perform feature
extraction, but have also successfully been applied to image generation (see
[4]). The following gives a short review of feedforward neural networks, the
simplest instance of a general neural network that can be used to construct
autoencoders. This will be used to introduce our neural network architecture
in Section 5.

####  3.2.1 Feedforward Neural Networks

Feedforward neural networks offer a method to build a basis for general
functions. In e.g. [18] it is shown that feedforward neural networks can
approximate any continuous function in a certain topology:

######  Theorem 1 (Universal approximation theorem for feedforward neural
networks [18]).

The set of two-layer neural networks is uniformly dense on compacta in
$C^{0}(\mathbb{R}^{N})$, the set of continuous functions
$\mathbb{R}^{N}\to\mathbb{R}$ if the activation function is non-polynomial.

The theorem, of course, also holds for neural networks that have more than two
layers. A layer of a feedforward neural network consists of an application of
an affine transformation, i.e., $x\mapsto Ax+b$ (where
$A\in\mathbb{R}^{m\times{}n}$ and $b\in\mathbb{R}^{m}$) and a
nonlinearity888Popular choices for nonlinearities include the nonlinear
rectified unit (RELU) $x\mapsto\mathrm{max}(0,x)$ and
$\mathrm{tanh}(x)=(1-e^{-x})/(1+e^{-x})$. $\sigma$. The last layer of a
feedforward neural network is typically just a linear transformation without
the nonlinearity. Theorem 1 states that the mapping $x\mapsto
v^{T}\sigma(Ax+b)$ (with $v\in\mathbb{R}^{m}$) can approximate any function
(on a compact set) if $m$ is chosen big enough.

The problem of finding the parameters $A$, $b$ and $v$ that parametrize the
function is usually solved in an iterative manner that is referred to as the
training of the neural network. The training aims to minimize a loss function
that maps the weights $\theta:=(A,b,v)$ to the non-negative real numbers
$\mathbb{R}^{+}:=[0,\infty)$. For the autoencoder example with feedforward
neural networks as encoder and decoder, the loss function is the
reconstruction error:

$L:\theta\equiv(\theta_{1},\theta_{2})\mapsto
d(M,\Psi^{\mathrm{dec}}_{\theta_{1}}\circ\Psi^{\mathrm{enc}}_{\theta_{2}}(M)).$
(22)

The simplest form of training is gradient descent with a constant learning
rate $\eta$, where the optimization is performed iteratively:
$\theta^{i+1}\leftarrow\theta^{i}-\eta\nabla_{\theta^{i}}L$. More involved
optimization techniques are discussed in Section 6.

Section 4 gives a brief discussion of SympNets. These are a special type of
feedforward neural network that can approximate any symplectic map. To this
end the space of maps that a neural network can approximate is restricted to
symplectic maps.

##  4 SympNets

SympNets (see [20]) offer a way of constructing general symplectic maps with
simple building blocks. The collections of networks which can be constructed
with these building blocks are called $\Psi_{LA}$ (sympnets made up of linear
and activation modules) and $\Psi_{G}$ (sympnets made up of gradient modules).
$\Psi_{LA}$ and $\Psi_{G}$ are dense in $\mathcal{SP}^{r}(\mathbb{R}^{2N})$,
the set of $r$-differentiable symplectic maps (for arbitrary $r$), in some
topology (similarly to theorem 1). More concretely:

######  Theorem 2 (Universal approximation theorem for SympNets [20]).

For any $r\in\mathbb{N^{*}}$, $\Psi_{LA}$ and $\Psi_{G}$ are both
$r$-uniformly dense on compacta in $\mathcal{SP}^{r}(\mathbb{R}^{2N})$ if the
activation function $\sigma$ is $r$-finite.

$\sigma$ being $r$-finite means that $0<\int|D^{r}(\sigma(x))|dx<\infty$.
$r$-uniformly dense on compacta means that for each
$g\in\mathcal{SP}^{r}(\mathbb{R}^{2N})$, $K$ compact, $\epsilon>0$ there
exists $\psi_{g,K,\epsilon}\in\Psi_{LA}/\Psi_{G}$ such that
$||g-\psi_{g,K,\epsilon}||_{r}<\epsilon$, where $||\cdot||_{r}$ is the Sobolev
norm.

The components making up $\Psi_{LA}$ and $\Psi_{G}$ will now briefly be
reviewed.

###  4.1 Linear/Activation-type SympNets $\Psi_{LA}$

Linear SympNet layers consist of upper/lower triangular matrices of the
following form:

$\begin{pmatrix}\mathbb{I}_{N}&S\\\
\mathbb{O}_{N}&\mathbb{I}_{N}\end{pmatrix}\text{ and
}\begin{pmatrix}\mathbb{I}_{N}&\mathbb{O}_{N}\\\
S^{\prime}&\mathbb{I}_{N}\end{pmatrix},$ (23)

where $S$ and $S^{\prime}$ are symmetric matrices. These matrices are by
construction symplectic.

In [21] it was shown that every symplectic matrix (i.e. an element of
$Sp(2N,2N)\equiv{}Sp(2N)$) can be represented by 5 of these linear layers.

Activation layers form the remaining component(s) of $\Psi_{LA}$. They take
the following form:

$\begin{pmatrix}q\\\
p\end{pmatrix}\mapsto\begin{pmatrix}q+\mathrm{diag}(a)\sigma(p)\\\
p\end{pmatrix}\text{ and }\begin{pmatrix}q\\\
p\end{pmatrix}\mapsto\begin{pmatrix}q\\\
p+\mathrm{diag}(a^{\prime})\sigma(q)\end{pmatrix},$ (24)

where $a\in\mathbb{R}^{n}$ and $[\mathrm{diag}(a)]_{ij}=a_{i}\delta_{ij}$.
$\Psi_{LA}$ is then the set of all possible combinations of the layers
presented in equations (23) and (24).

###  4.2 Gradient-type SympNets $\Psi_{G}$

The third type of SympNet layers are gradient layers. They have the following
form:

$\begin{pmatrix}q\\\
p\end{pmatrix}\mapsto\begin{pmatrix}q+K^{T}\mathrm{diag}(a)\sigma(Kp+b)\\\
p\end{pmatrix}\text{ and }\begin{pmatrix}q\\\
p\end{pmatrix}\mapsto\begin{pmatrix}q\\\ p+K^{\prime
T}\mathrm{diag}(a^{\prime})\sigma(K^{\prime}q+b)\end{pmatrix},$ (25)

where $K\in\mathbb{R}^{L\times{}N}$ and $a,b\in\mathbb{R}^{L}$ ith $L$
arbitrary. $\Psi_{G}$ is any combination of layers of the type presented in
equation (25).

The symplectic autoencoders presented in this work use SympNets and PSD-like
matrices to construct mappings that preserve symplecticity and also change
dimension. The first mapping in equation (25) will be refered to as GradientQ
and the second as GradientP as they change the $q$ and the $p$ component
respectively.

##  5 Proposed Architecture

The following illustrates the proposed architecture:

[width=.8]autoencoders3

It is a composition of SympNets and PSD-like matrices:

$\displaystyle\Psi^{\mathrm{enc}}$
$\displaystyle=\psi^{L}_{\mathrm{symp}}\circ(A_{PSD}^{L})^{+}\circ\cdots\circ\psi^{2}_{\mathrm{symp}}\circ(A_{PSD}^{1})^{+}\circ\psi^{1}_{\mathrm{symp}}\text{
and}$ (26) $\displaystyle\Psi^{\mathrm{dec}}$
$\displaystyle=\tilde{\psi}^{\tilde{L}}_{\mathrm{symp}}\circ\tilde{A}_{PSD}^{\tilde{L}}\circ\cdots\circ\tilde{\psi}^{2}_{\mathrm{symp}}\circ\tilde{A}_{PSD}^{1}\circ\tilde{\psi}^{1}_{\mathrm{symp}}.$
(27)

The SympNet layers here perform a symplectic preprocessing step before the
linear, symplectic reduction is employed with the PSD-like layers. The PSD-
like layers are matrices of the same form as in equation (14).

This then yields a neural network that is capable of approximating any
canonical symplectic mapping, also mappings that change the dimension of the
system - this is not possible with traditional SympNets.

One has to be careful with training this network since some of the weights
have to satisfy the symplecticity constraint. In order to train the SympNet
part of the network, standard general-purpose algorithms can be employed as
was also done in the original SympNet paper [20]. Optimizing the elements of
the Stiefel manifold(s) is more complicated, as the condition
$\Phi^{T}\Phi=\mathbb{I}_{n}$ has to be enforced in each training step. This
naturally leads to manifold optimization (see [5, 12, 1]).

The Stiefel manifold $St(n,N)\simeq\\{\text{$PSD$-like
matrices}\\}\subset\mathbb{R}^{2N\times 2n}$ is a homogeneous space for which
a generalization of common optimizers exist (see [8]). In addition the set
$St(n,N)$ is compact (whereas $Sp(2n,2N)$, for example, is not); this makes
optimization on $St(n,N)$ relatively easy.

Note that according to equation (12), if the jacobian of the reconstruction
$\Psi^{\mathrm{dec}}$ is symplectic, then one obtains a Hamiltonian system
with Hamiltonian $H\circ\Psi^{\mathrm{dec}}$ on the reduced space. This can be
achieved with the new architecture.

The following section summarizes the most important aspects of optimization on
manifolds relevant for neural networks.

##  6 Manifold Optimization for Neural Networks

This section discusses mathematical and computational aspects relevant for the
implementation of the new neural network architecture. This requires
generalizing existing optimizers to the manifolds that appear in the network,
i.e. $St(n,N)$. This manifold, as well as $Sp(2n,2N)$, are homogeneous spaces,
a special subclass of manifolds that make generalizing common optimizers
(stochastic gradient descent with momentum, Adam, etc.) possible. Details
regarding this are presented in [8]. For the purposes of this paper we limit
the discussion to standard first-order optimization on general Riemannian
manifolds, i.e. generalized gradient descent, and mention Adam only briefly.

###  6.1 Optimization on Riemannian manifolds and retractions

Here a generalization of the update rule $Y^{t+1}\leftarrow
Y^{t}-\eta\nabla_{Y^{t}}L$ from Euclidean spaces to manifolds is discussed. In
practice this involves two steps: the computation of a Riemannian gradient and
a retraction that (approximately) solves the geodesic equation.

Assume we are given a Riemannian manifold $(\mathcal{M},g)$ and an objective
function $L$ that we want to minimize on $\mathcal{M}$. Standard gradient
descent does this by repeatedly solving the geodesic equation. This is a
canonical second-order ODE associated with all Riemannian manifolds. This ODE
is obtained by applying the variational principle (see [24]) to the following
functional:

$\mathcal{L}(\gamma)=\int_{0}^{1}\sqrt{g_{\gamma(t)}(\dot{\gamma}(t),\dot{\gamma}(t))}.$
(28)

The geodesic equation, the result of this variational principle, is then
usually written in the following form:

$\frac{d^{2}\gamma^{\lambda}}{dt^{2}}+\Gamma^{\lambda}_{\mu\nu}\frac{d\gamma^{\mu}}{dt}\frac{d\gamma^{\nu}}{dt}=0,$
(29)

with $\Gamma^{\lambda}_{\mu\nu}$ begin the Christoffel symbols of the second
kind (see [6, Chapter 5.12]); the initial conditions are $\gamma(0)=Y$ and
$\dot{\gamma}(0)=W$. It is worth mentioning that we do not directly work with
equation (29) and Christoffel symbols, but exploit the special structure of
the manifolds we deal with (as homogeneous spaces) to make solving the
geodesic equation computationally more tractable. For optimizing neural
networks with weights on manifolds we now need a value for $W$ (i.e. an
initial condition) and a method to solve or approximate equation (29).

First we address how to obtain the value for $W$, i.e. the initial condition
for the second order ODE. If we are optimizing neural networks with weights on
vector spaces, $W$ is the result of automatic differentiation (AD, [19, 29])
multiplied with the negative of a learning rate, i.e. $W=-\eta\nabla{}L$. If
we are instead dealing with more general manifolds, we have to perform another
operation to obtain an element of the tangent space
$W\in{}T_{Y^{t}}\mathcal{M}$. This additional operation is called computing
the Riemannian gradient:

######  Definition 1.

The Riemannian gradient of a function $L$ on a Riemannian manifold
$(\mathcal{M},g)$ is an element of the tangent space
$\mathrm{grad}_{Y}L\in{}T_{Y}\mathcal{M}$ such that
$\forall{}V\in{}T_{Y}\mathcal{M}:\langle{}dL,V\rangle=g_{Y}(\mathrm{grad}_{Y}L,V)$.

For $St(n,N)$ and $Sp(2n,2N)$ the Riemannian gradient can be obtained through
a simple operation on the output of an AD routine999The output of the AD
routine is called $\nabla_{Y}L$ here, i.e. the Euclidean gradient of $L$ at
$Y$.:

$\displaystyle\mathrm{grad}_{Y}L=$
$\displaystyle\nabla_{Y}L-Y(\nabla_{Y}L)^{T}Y$ for $Y\in$ $\displaystyle
St(n,N),$ and (30) $\displaystyle\mathrm{grad}_{Y}L=$
$\displaystyle\nabla_{Y}L(Y^{T}Y)+\mathbb{J}_{2N}Y((\nabla_{Y}L)^{T}\mathbb{J}_{2N}Y)$
for $Y\in$ $\displaystyle Sp(2n,2N).$ (31)

Also see appendix A for a discussion of these expressions and [1, 5, 12, 23]
for similar computations. After having obtained the Riemannian gradient of $L$
at $Y^{t}$ one can solve or approximate the geodesic equation. For $St(n,N)$
and $Sp(2n,2N)$ there exist closed-form solutions for the geodesic equation,
see [10, equation (2.25)] and [5, equation (3.18)]. Alternatively one could
approximate the solution of the geodesic with so-called retractions, one of
the most popular of these is the Cayley transform (see e.g. [17, chaber IV]
and [1, chapter 4.1]). The solution of equation (29) with $\gamma(0)=Y$ and
$\dot{\gamma}(0)=V\in{}T_{Y}\mathcal{M}$ will be called
$\gamma_{\mathcal{M}}(t;V)$.

To sum up, the update rule $Y^{t+1}\leftarrow{}Y^{t}-\eta\nabla_{Y^{t}}L$ on
regular vector spaces is replaced with:

$Y^{t+1}\leftarrow\gamma_{\mathcal{M}}(1;-\eta\mathrm{grad}_{Y}L).$ (32)

So we first compute the Euclidean gradient $\nabla_{Y}L$ with an AD routine,
then obtain the Riemannian gradient with equation (30) and then use this as an
initial condition to solve (or approximate) the geodesic equation in (29).

###  6.2 Generalizing Adam

Modern training of neural networks does not simply perform gradient descent,
but uses optimization algorithms that include first and second moments (called
$s$ and $r$ in the algorithm below), which we also refer to as the cache. One
of the most successful optimization algorithms is Adam [22]:

1:step size $\eta$; exponential decay rates $\beta_{1},\beta_{2}\in[0,1)$;
constant for numerical stabilization $\delta$; first and second moments
$s^{t-1}$, $r^{t-1}$;

2:compute $G_{\theta}\leftarrow\nabla_{\theta}{}L(\theta^{t-1})$,

3:update
$s^{t}\leftarrow\frac{\beta_{1}-\beta_{1}^{t}}{1-\beta_{1}^{t}}s^{t-1}+\frac{1-\beta_{1}}{1-\beta_{1}^{t}}G_{\theta}$,

4:update
$r^{t}\leftarrow\frac{\beta_{2}-\beta_{2}^{t}}{1-\beta_{2}^{t}}r^{t-1}+\frac{1-\beta_{2}}{1-\beta_{2}^{t}}G_{\theta}\odot{}G_{\theta}$,

5:update
$\theta^{t}\leftarrow\theta^{t-1}-\eta\frac{s^{t}}{\sqrt{r^{t}+\delta}}$
$\triangleright$ operations are applied element-wise.

Algorithm 1 Pseudocode for a single gradient step with Adam.

The Adam algorithm (and other neural network optimizers) take as input the
current neural network parameters $\theta^{t-1}$ and a cache (in the Adam
algorithm above this is $s^{t-1}$ and $r^{t-1}$) and updates both of them.
Some aspects of this algorithm are easy to generalize to manifolds, other are
not. The following offers a short description of this, breaking the alorithm
into three main parts:

  1. 1.

Computing a first order derivative of the loss function with respect to the
network parameters. For usual neural network weights this is the output of an
AD routine (the Euclidean gradient). This can be easily generalized to
manifolds by computing the Riemannian gradient (see equation (30)).

  2. 2.

Updating the cache and computing a final velocity. In the vector space case
this final velocity is $v^{t}=-\eta{}s^{t}/(r^{t}+\delta)$. This is the step
that is very difficult to generalize to manifolds since many operations (such
as the Hadamard product) in the Adam algorithm lack a clear physical
interpretation.

  3. 3.

Performing the update. For the vector space case this is simply adding the
velocity onto the neural network weights, i.e.
$\theta^{t}\leftarrow\theta^{t-1}+v^{t}$. For manifolds this step consists of
solving (or approximating) the geodesic equation (29).

In geometric terms, step 2 above (i.e. he Hadamard product and associated
operations) depends on the chosen coordinate system and is inherently non-
geometric. A change of coordinates produces a different result. A workaround
for this exists by leveraging the goemetry of homogeneous spaces, a category
to which the Stiefel manifold and the symplectic Stiefel manifold belong.
Working on these we can identify a global tangent space representation
$\mathfrak{g}^{\mathrm{hor}}$. This is the space on which the cache is defined
and all the Adam operations that lack a clear geometric interpretation are
performed, i.e. the operations $\odot$, $\sqrt{(\cdot)}$, $/$ and
$(\cdot)+\delta$ in algorithm 1. The full generalization of Adam to
homogeneous spaces is discussed in [8]. This is also the algorithm that we use
for training the networks in the following section.

##  7 Experiments

We demonstrate the new method on the example of the linear wave equation (see
equation (6)) with $\xi\in\Omega:=(-1/2,1/2)$ and
$\mu\in\mathbb{P}:=[5/12,2/3]$ as a choice for domain and parameters
(similarly to [9]). We additionaly specify the following initial and boundary
conditions:

$\displaystyle q_{0}(\omega;\mu):=$ $\displaystyle q(0,\omega;\mu)$ (33)
$\displaystyle p(0,\omega;\mu)=\partial_{t}q(0,\xi;\mu)=$
$\displaystyle-\mu\partial_{\omega}{}q_{0}(\xi;\mu)$ (34) $\displaystyle
q(t,\omega;\mu)=$ $\displaystyle 0,\text{ for }\omega\in\partial\Omega.$ (35)

The precise shape of $q_{0}(\cdot;\cdot)$ is described below.

As with any other PDE, the wave equation can also be discretized to obtain a
ODE which can be solved numerically.

If we discretize $\mathcal{H}$ directly, to obtain a Hamiltonian on a finite-
dimensional vector space, we get a Hamiltonian ODE101010This conserves the
Hamiltonian structure of the system.:

$\mathcal{H}_{h}(z)=\sum_{i=1}^{\tilde{N}}\frac{\Delta{}x}{2}\bigg{[}p_{i}^{2}+\mu^{2}\frac{(q_{i}-q_{i-1})^{2}+(q_{i+1}-q_{i})^{2}}{2\Delta{}x}\bigg{]}=\frac{\Delta{}x}{2}p^{T}p+z^{T}Kz,$
(36)

where the components of the matrix $K$ are of the form:

$k_{ij}=\begin{cases}\frac{\mu^{2}}{4\Delta{}x}&\text{if
$(i,j)\in\\{(0,0),(\tilde{N}+1,\tilde{N}+1)\\}$ },\\\
-\frac{\mu^{2}}{2\Delta{}x}&\text{if $(i,j)=(1,0)$ or
$(i,j)=(\tilde{N},\tilde{N}+1)$}\\\ \frac{3\mu^{2}}{4\Delta{}x}&\text{if
$(i,j)\in\\{(1,1),(\tilde{N},\tilde{N})\\}$}\\\
\frac{\mu^{2}}{\Delta{}x}&\text{if $i=j$ and
$i\in\\{2,\ldots,(\tilde{N}-2)\\}$}\\\ -\frac{\mu^{2}}{2\Delta{}x}&\text{if
$|i-j|=1$ and $i,j\notin\\{0,\tilde{N}+1\\}$}\\\ 0&\text{else}.\end{cases}$
(37)

and
$q=(q_{0},q_{1},\cdots,q_{\tilde{N}})^{T},\,p=(p_{0},p_{1},\cdots,p_{\tilde{N}})^{T}$.

$\tilde{N}$ is the number of grid points, which leads to a symplectic vector
space of dimension $2N:=2(\tilde{N}+2)$.

The vector field of the FOM is described by (see [30]):

$\frac{dz}{dt}=\mathbb{J}_{d}\nabla_{z}\mathcal{H}_{h}=\mathbb{J}_{d}\begin{bmatrix}\Delta{}x\mathbb{I}&\mathbb{O}\\\
\mathbb{O}&K+K^{T}\end{bmatrix}z,\quad\mathbb{J}_{d}=\frac{\mathbb{J}_{2N}}{\Delta{}x}.$
(38)

The parameter-dependence of this system is encoder into the matrix $K$. Note
that, as opposed to equation (8), this system is not time-dependent. In the
following section we describe how we obtain the initial conditions.

###  7.1 Initial conditions

The initial conditions are based on the following third-degree spline (also
used in [9]):

$h(s)=\begin{cases}1-\frac{3}{2}s^{2}+\frac{3}{4}s^{3}&\text{if }0\leq s\leq
1\\\ \frac{1}{4}(2-s)^{3}&\text{if }1<s\leq 2\\\ 0&\text{else.}\end{cases}$
(39)

Plotted on the relevant domain it takes the following shape:

[width=.4]third_degree_spline

Figure 2: The third-degree spline from equation (39). This is used to model
all initial conditions.

Taking the above function $h(s)$ as a starting point, the initial conditions
for the linear wave equations are modelled with

$q_{0}(\omega;\mu)=h(s(\omega,\mu)).$ (40)

For the construction of $s(\cdot;\cdot)$ we consider the following factors:

  1. 1.

the analytic solution for equation (6) is
$q(t,\omega;\mu)=q_{0}(\omega-\mu{}t;\mu)$; i.e. the solutions of the linear
wave equation will travel with speed $\mu$. We should make sure that the wave
does not touch the right boundary of the domain for any of the values
$t\in\mathcal{I}$. So the peak should be sharper for higher values of $\mu$ as
the wave will travel faster.

  2. 2.

the wave should start at the left boundary of the domain in order to cover it
as much as possible.

Based on this we end up with the following choice of parametrized initial
conditions:

$q_{0}(\mu)(\omega)=h(s(\omega,\mu)),\quad
s(\omega,\mu)=20\mu\left|\omega+\frac{\mu}{2}\right|.$ (41)

Three initial conditions and their time evolutions are shown in figure 3. As
was required, we can see that the peak gets sharper and moves to the left as
we increase the parameter $\mu$; the curves also get a good coverage of the
domain $\Omega$.

Figure 3: Numerical solutions of equation (6) for three different values of
$\mu$. Note that the shape of the curves as well as the speed with which they
travel changes based on the parameter.

The wave equation has a slowly-decaying Kolmogorov $n$-width (see e.g. [15]),
which means linear methods like PSD will perform poorly.

For discretizing the system we choose $\tilde{N}$ in equation (36) to be
$\tilde{N}=128$, resulting in a FOM dimension of $2N=2(\tilde{N}+2)=260$; here
we get two additional degrees of freedom for the $q$ and $p$ coordinates
because we also have to consider the end-points (boundary conditions). We
integrate the system in equation (38) with implicit midpoint, which is a
general-purpose symplectic integrator (see [17]).

We train on data coming from 20 different choices of the parameter $\mu$ that
are equally spaced in the interval $\mathbb{P}=[5/12,2/3]$; the training is
done for 100 epochs. The precise architecture is shown in figure 4.

[width=.8]architecture

Figure 4: The architecture we use for all experiments. GradientQ and GradientP
are the gradient layers (see section 4) that change the $q$ and $p$ component
respectively. The syntax is the same used as in GeometricMachineLearning.jl.

This system was already treated with so-called “weakly-symplectic
autoencoders” ([9]). These do not enforce symplecticity through the neural
network architecture, as we do it here, but attempt to take symplecticity into
account through an extra term in the loss function (PiNN approach, see [31]).
This has obvious pitfalls as there are no guarantees on symplecticity which
lead to unpredictable behaviour in the online stage that is also discussed in
[9].

As was done in [9], we compare two errors: the projection error and the
reduction error; these are:

$e_{\mathrm{proj}}(\mu):=\sqrt{\frac{\sum_{t=0}^{K}||\mathbf{x}^{(t)}(\mu)-\mathcal{R}\circ\mathcal{P}(\mathbf{x}^{(t)}(\mu))||^{2}}{\sum_{t=0}^{K}||\mathbf{x}^{(t)}(\mu)||^{2}}},$
(42)

and

$e_{\mathrm{red}}(\mu):=\sqrt{\frac{\sum_{t=0}^{K}||\mathbf{x}^{(t)}(\mu)-\mathcal{R}(\mathbf{x}^{(t)}_{r}(\mu))||^{2}}{\sum_{t=0}^{K}||\mathbf{x}^{(t)}(\mu)||^{2}}},$
(43)

where $||\cdot||$ is the $L_{2}$ norm (one could also optimize for different
norms), $\mathbf{x}^{(t)}$ is the solution of the FOM at point $t$ and
$\mathbf{x}^{(t)}_{r}$ is the solution of the ROM (in the reduced basis) at
point $t$. This is obtained by integraing the reduced system (second vector
field in equation (12)). The projection error computes how well a reduced
basis, represented by the reduction $\mathcal{P}$ and the reconstruction
$\mathcal{R}$, can represent the data with which it is built. The reduction
error on the other hand measures how far the reduced system diverges from the
full-order system during integration (online stage).

(a) $\mu=0.47$

(b) $\mu=0.51$

(c) $\mu=0.55$

(d) $\mu=0.625$

Figure 5: Projection and reduction errors for the different choices of the
parameter $\mu$.

Training of the neural networks is done with
GeometricMachineLearning.jl111111https://github.com/JuliaGNI/GeometricMachineLearning.jl.
This contains almost all relevent components for optimizing the neural network
shown in figure 4, including the generalization of common optimizers to
manifolds. (It does however not include an automatic differentiation routine.
This is taken from the Julia package
Zygote.jl121212https://github.com/FluxML/Zygote.jl.) As optimizer we pick the
Adam optimzier (see [22]) with the values $\eta=0.001$, $\rho_{1}=0.9$,
$\rho_{2}=0.99$ and $\delta=10^{-8}$. The optimization for the PSD-like layers
is done with the generalized optimizer introduced in [8]. The data and network
parameters are all of type Float64.

The encoder and decoder are evaluated for four different choices of the
parameter $\mu$: 0.470, 0.510, 0.550 and 0.625; none of these appear in the
training set. The results are shown in figure 5 including a comparison to PSD.

The symplectic autoencoders almost consistently outperform the PSD for all of
these parameters, exept for $\mu=0.625$. We also observe a monotonous decrease
in the reduction error as we increase the dimension of the reduced space. This
was not observed for the “weakly-symplectic autoencoders” in [9], where
increasing the dimension of the reduced space did not necessarily imply an
improvement in the reduction error.

##  8 Summary

We introduced a novel neural network architecture that, for the first time,
can map between spaces of different dimension while simultaneously preserving
symplecticity. In order to achieve this we rely on manifold optimization and a
generalized Adam optimizer.

The advantages of the new neural network design was demonstrated by comparing
it to proper symplectic decomposition (PSD). PSD was up to this point the only
data-driven method that can map between spaces of different dimension while
also preserving symplecticity. PSD is however, unlike the neural network
introduced in this paper, strictly limited to linear mappings.

## Acknowledgements

We would like to thank Tobias Blickhan for valuable discussions and Tatjana
Stykel for reading an early version of the manuscript and providing valuable
remarks that greatly improved the exposition of this work.

##  Appendix A Riemannian gradients for the Stiefel and the Symplectic Stiefel
manifold

Here we check that the definitions for the Riemannian gradients in equation
(30) are the correct ones. The Riemannian metrics for $St(n,N)$ and
$St(2n,2N)$ are:

$\displaystyle g_{Y}(V_{1},V_{2})$
$\displaystyle=\mathrm{Tr}(V_{1}^{T}\left(\mathbb{I}_{2N}-\frac{1}{2}YY^{T}\right)V_{2})$
and (44) $\displaystyle\hat{g}_{Y}(W_{1},W_{2})$
$\displaystyle=\mathrm{Tr}(W_{1}^{T}\left(\mathbb{I}_{2N}-\frac{1}{2}\mathbb{J}_{2N}^{T}Y(Y^{T}Y)^{-1}Y^{T}\mathbb{J}_{2N}\right)W_{2}(Y^{T}Y)^{-1}).$
(45)

For the first of these equations see e.g. [8] and for the second see e.g. [5].
See definition 1 for the Riemannian gradient.

Here $\langle\cdot,\cdot\rangle$ is the dual pairing between the tangent space
and its dual. For a matrix manifold $\mathcal{M}$ and a function
$L:\mathcal{M}\to\mathbb{R}$ we have:
$\langle{}dL,V\rangle=\mathrm{Tr}((\nabla_{Y}L)^{T}V)$.

It can now easily be verified that equation (30) shows the correct Riemannian
gradients. For the first case:

$\displaystyle g_{Y}(\mathrm{grad}_{Y}L,V)$
$\displaystyle=\mathrm{Tr}((\nabla_{Y}L)^{T}V-(Y^{T}\nabla{}L)Y^{T}V-\frac{1}{2}(\nabla_{Y}L)^{T}YY^{T}V+\frac{1}{2}(Y^{T}\nabla{}L)Y^{T}YY^{T}V)$
(46)
$\displaystyle=\mathrm{Tr}((\nabla_{Y}L)^{T}V-(Y^{T}\nabla{}L)Y^{T}V-\frac{1}{2}(\nabla_{Y}L)^{T}YY^{T}V+\frac{1}{2}(Y^{T}\nabla{}L)Y^{T}V)$
(47)
$\displaystyle=\mathrm{Tr}((\nabla_{Y}L)^{T}V-\frac{1}{2}(Y^{T}\nabla{}L)Y^{T}V-\frac{1}{2}(\nabla_{Y}L)^{T}YY^{T}V)$
(48) $\displaystyle=Tr((\nabla_{Y}L)^{T}V).$ (49)

Here we have used that $Y^{T}V=-V^{T}Y$ and $\mathrm{Tr}(AB)=\mathrm{Tr}(BA)$
for arbitrary (compatible) matrices $A$ and $B$ as well as
$\mathrm{Tr}(A)=\mathrm{Tr}(A^{T})$ for $A$ a square matrix. To sum up:

$\displaystyle-\mathrm{Tr}((Y^{T}\nabla_{y}L)Y^{T}V)-\mathrm{Tr}((\nabla_{Y}L)^{T}YY^{T}V)$
$\displaystyle=-\mathrm{Tr}((\nabla_{Y}L)Y^{T}VY^{T})+\mathrm{Tr}((\nabla_{Y}L)^{T}YV^{T}Y)$
(50)
$\displaystyle=\mathrm{Tr}((\nabla_{Y}L)Y^{T}VY^{T})+\mathrm{Tr}(YV^{T}Y(\nabla_{Y}L)^{T})=0.$
(51)

This concludes the proof that the first line in equation (30) shows the
correct Riemannian gradient for $St(n,N)$. With a similar calculation we can
also do this for $Sp(2n,2N)$ (also see [5]).

## References

  * Absil et al. [2008] P-A Absil, Robert Mahony, and Rodolphe Sepulchre.  _Optimization algorithms on matrix manifolds_.  Princeton University Press, 2008.
  * Antoulas et al. [2001] Athanasios C Antoulas, Danny C Sorensen, and Serkan Gugercin.  A survey of model reduction methods for large-scale systems.  _Contemp.Math._ , 280:193–220, 2001.
  * Arnold [1978] VI Arnold.  Mathematical methods of classical mechanics.  _Ann Arbor_ , 1001:48109, 1978.
  * Bank et al. [2020] Dor Bank, Noam Koenigstein, and Raja Giryes.  Autoencoders.  _arXiv preprint arXiv:2003.05991_ , 2020.
  * Bendokat and Zimmermann [2021] Thomas Bendokat and Ralf Zimmermann.  The real symplectic stiefel and grassmann manifolds: metrics, geodesics and applications.  _arXiv preprint arXiv:2108.12447_ , 2021.
  * Bishop and Goldberg [1980] Richard L Bishop and Samuel I Goldberg.  _Tensor Analysis on Manifolds_.  Courier Corporation, 1980.
  * Blickhan [2023] Tobias Blickhan.  A registration method for reduced basis problems using linear optimal transport.  _arXiv preprint arXiv:2304.14884_ , 2023.
  * Brantner [2023] Benedikt Brantner.  Generalizing adam to manifolds for efficiently training transformers.  _arXiv preprint arXiv:2305.16901_ , 2023.
  * Buchfink et al. [2023] Patrick Buchfink, Silke Glas, and Bernard Haasdonk.  Symplectic model reduction of hamiltonian systems on nonlinear manifolds and approximation with weakly symplectic autoencoder.  _SIAM Journal on Scientific Computing_ , 45(2):A289–A311, 2023.  doi: 10.1137/21M1466657.  URL https://doi.org/10.1137/21M1466657.
  * Edelman et al. [1998] Alan Edelman, Tomás A Arias, and Steven T Smith.  The geometry of algorithms with orthogonality constraints.  _SIAM journal on Matrix Analysis and Applications_ , 20(2):303–353, 1998.
  * Fresca et al. [2021] Stefania Fresca, Luca Dede, and Andrea Manzoni.  A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized pdes.  _Journal of Scientific Computing_ , 87(2):1–36, 2021.
  * Gao et al. [2021] Bin Gao, Nguyen Thanh Son, P-A Absil, and Tatjana Stykel.  Riemannian optimization on the symplectic stiefel manifold.  _SIAM Journal on Optimization_ , 31(2):1546–1575, 2021.
  * Gao et al. [2022] Bin Gao, Nguyen Thanh Son, and Tatjana Stykel.  Optimization on the symplectic stiefel manifold: Sr decomposition-based retraction and applications.  _arXiv preprint arXiv:2211.09481_ , 2022.
  * Goodfellow et al. [2016] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  _Deep learning_.  MIT press, 2016.
  * Greif and Urban [2019] Constantin Greif and Karsten Urban.  Decay of the kolmogorov n-width for wave problems.  _Applied Mathematics Letters_ , 96:216–222, 2019.
  * Gromov [1985] Mikhael Gromov.  Pseudo holomorphic curves in symplectic manifolds.  _Inventiones mathematicae_ , 82(2):307–347, 1985\.
  * Hairer et al. [2006] Ernst Hairer, Christian Lubich, and Gerhard Wanner.  _Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations_.  Springer, 2006.
  * Hornik et al. [1989] Kurt Hornik, Maxwell Stinchcombe, and Halbert White.  Multilayer feedforward networks are universal approximators.  _Neural networks_ , 2(5):359–366, 1989.
  * Innes et al. [2019] Mike Innes, Alan Edelman, Keno Fischer, Chris Rackauckas, Elliot Saba, Viral B Shah, and Will Tebbutt.  A differentiable programming system to bridge machine learning and scientific computing.  _arXiv preprint arXiv:1907.07587_ , 2019.
  * Jin et al. [2020] Pengzhan Jin, Zhen Zhang, Aiqing Zhu, Yifa Tang, and George Em Karniadakis.  Sympnets: Intrinsic structure-preserving symplectic networks for identifying hamiltonian systems.  _Neural Networks_ , 132:166–179, 2020.
  * Jin et al. [2022] Pengzhan Jin, Zhangli Lin, and Bo Xiao.  Optimal unit triangular factorization of symplectic matrices.  _Linear Algebra and its Applications_ , 2022.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kong et al. [2022] Lingkai Kong, Yuqing Wang, and Molei Tao.  Momentum stiefel optimizer, with applications to suitably-orthogonal attention, and optimal transport.  _arXiv preprint arXiv:2205.14173_ , 2022.
  * Kraus [2013] Michael Kraus.  Variational integrators in plasma physics.  _arXiv preprint arXiv:1307.5665_ , 2013.
  * Lassila et al. [2014] Toni Lassila, Andrea Manzoni, Alfio Quarteroni, and Gianluigi Rozza.  Model order reduction in fluid dynamics: challenges and perspectives.  _Reduced Order Methods for modeling and computational reduction_ , pages 235–273, 2014.
  * Lee and Carlberg [2020] Kookjin Lee and Kevin T Carlberg.  Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.  _Journal of Computational Physics_ , 404:108973, 2020.
  * Leimkuhler and Reich [2004] Benedict Leimkuhler and Sebastian Reich.  _Simulating hamiltonian dynamics_.  Cambridge university press, Cambridge, 2004.
  * McDuff and Salamon [2017] Dusa McDuff and Dietmar Salamon.  _Introduction to symplectic topology_ , volume 27.  Oxford University Press, 2017.
  * Moses et al. [2021] William S Moses, Valentin Churavy, Ludger Paehler, Jan Hückelheim, Sri Hari Krishna Narayanan, Michel Schanen, and Johannes Doerfert.  Reverse-mode automatic differentiation and optimization of gpu kernels via enzyme.  In _Proceedings of the international conference for high performance computing, networking, storage and analysis_ , pages 1–16, 2021.
  * Peng and Mohseni [2016] Liqian Peng and Kamran Mohseni.  Symplectic model reduction of hamiltonian systems.  _SIAM Journal on Scientific Computing_ , 38(1):A1–A27, 2016.
  * Raissi et al. [2019] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational physics_ , 378:686–707, 2019\.
  * Sanz-Serna and Calvo [2018] Jesus-Maria Sanz-Serna and Mari-Paz Calvo.  _Numerical hamiltonian problems_ , volume 7.  Dover Publications, Mineola, NY, 2018.
  * Tyranowski and Kraus [2022] Tomasz M Tyranowski and Michael Kraus.  Symplectic model reduction methods for the vlasov equation.  _Contributions to Plasma Physics_ , 2022.
  * Xu [2003] Hongguo Xu.  An svd-like matrix decomposition and its applications.  _Linear Algebra and its Applications_ , 368:1–24, 2003\.  ISSN 0024-3795.  doi: https://doi.org/10.1016/S0024-3795(03)00370-7.  URL https://www.sciencedirect.com/science/article/pii/S0024379503003707.