Symplectic model reduction of Hamiltonian systems using data-driven quadratic manifolds

ucsd]Harsh Sharmacor1
[cor1]Corresponding author
affil3]Hongliang Mu
affil4]Patrick Buchfink
affil2]Rudy Geelen
affil3]Silke Glas

ucsd]Boris Kramer

[ucsd]Department of Mechanical and Aerospace Engineering, University of California San Diego, CA, United States

[affil3] Department of Applied Mathematics, University of Twente, Enschede, The Netherlands

[affil4] Institute of Applied Analysis and Numerical Simulation, University of Stuttgart, Stuttgart, Germany
[affil2] Oden Institute for Computational Engineering and Sciences, University of Texas at Austin, TX, United States

This work presents two novel approaches for the symplectic model reduction of high-dimensional Hamiltonian systems using data-driven quadratic manifolds. Classical symplectic model reduction approaches employ linear symplectic subspaces for representing the high-dimensional system states in a reduced-dimensional coordinate system. While these approximations respect the symplectic nature of Hamiltonian systems, linear basis approximations can suffer from slowly decaying Kolmogorov $N$-width, especially in wave-type problems, which then requires a large basis size. We propose two different model reduction methods based on recently developed quadratic manifolds, each presenting its own advantages and limitations. The addition of quadratic terms to the state approximation, which sits at the heart of the proposed methodologies, enables us to better represent intrinsic low-dimensionality in the problem at hand. Both approaches are effective for issuing predictions in settings well outside the range of their training data while providing more accurate solutions than the linear symplectic reduced-order models.

Symplectic model reductionHamiltonian systems Data-driven modeling Quadratic manifolds Scientific machine learning

§ INTRODUCTION

Computational modeling, simulation, and control of dynamical systems characterized by Hamiltonian mechanics are essential for many science and engineering applications such as robotics, quantum mechanics, solid state physics, and climate modeling, see [37, 39]. Hamiltonian partial differential equations (PDEs) demonstrate complex dynamic behavior while possessing underlying mathematical structures in the form of symmetries, Casimirs, symplecticity, and energy conservation. The structure-preserving spatial discretization of Hamiltonian PDEs with finite element methods, spectral methods, or finite difference methods leads to finite-dimensional Hamiltonian systems with large state-space dimensions.

For a high-dimensional Hamiltonian system, it is desirable to construct a reduced-order model (ROM) that preserves the underlying geometric structure and in particular, the stability of the full-order model (FOM). To this end, the field of structure-preserving model reduction has developed efficient tools for reducing the computational complexity of high-dimensional FOMs while preserving the underlying geometric structure. Symplectic model reduction of Hamiltonian systems was introduced in [50], where the Galerkin projection-based ROM was modified so that the ROM retains the underlying symplectic structure. In another research direction, the work in [38] presented a greedy approach for symplectic model reduction of parametric Hamiltonian systems. A variety of symplectic model reduction approaches have been developed along these two research directions recently, e.g., basis generation for symplectic model reduction of Hamiltonian systems [6, 11, 8, 9], modified projection techniques [21, 29, 44], and data-driven operator inference [58, 56, 25]. However, all of the aforementioned methods for symplectic model reduction project onto linear subspaces which often requires simulating high-dimensional ROMs to obtain accurate approximate solutions.

Classical projection-based model reduction methods encounter substantial difficulties whenever the solution manifold has so-called Kolmogorov $N$-widths that decay slowly with increasing $N$ [35]. The Kolmogorov $N$-width quantifies the degree of accuracy by which a set can be approximated using linear subspaces of dimension $N$. If the $N$-widths decay slowly with increasing $N$, the associated ROMs can provide accurate approximate solutions only for large values of $N$, which in turn defeats the original purpose of achieving computational speedups. This poses substantial limitations to the spectrum of possible applications. Many physics-based systems do not exhibit a global low-rank structure and are characterized by slowly decaying Kolmogorov $N$-widths, such as those featuring coherent structures that propagate in space-time [24, 43]. Strategies aimed at overcoming this bottleneck typically fall into one of three categories: (1) dictionary approaches (local reduced basis methods, divide-and-conquer, [3, 18, 48, 14, 13]), (2) time-dependent basis or other transformations [30, 46, 51, 41, 42, 59, 34, 29, 54], (3) and the use of nonlinear approximation techniques (deep learning and other data-driven manifold techniques) [15, 16, 17, 22, 52, 47, 2, 5], which are sometimes referred as model reduction on manifolds. In the present work, we focus on the latter class of methods. Following the pioneering work from [36], projection-based reduced-order modeling approaches based on convolutional autoencoders have seen a surge of activity in recent years [15, 16, 33, 32, 2]. With various levels of success, these methods can produce ROMs of lower dimensionality—compared to classical approaches based on linear subspaces—for the same level of accuracy. In general, however, nonlinear model reduction methods are not tailored to preserve the underlying geometric structure and are thus ill-suited for Hamiltonian systems.

An exception in that regard is the work [10], which uses a nonlinear approximation while preserving the structure of Hamiltonian systems.
There, the symplectic manifold Galerkin (SMG) projection is introduced which generalizes the structure-preserving model reduction with linear subspaces from [50] to symplectic (nonlinear) approximation mappings.
The SMG projection requires a symplectic approximation mapping. While the SMG projection works for all symplectic approximation mappings, the numerical example in [10] however only considers a weakly symplectic mapping. Moreover, the mapping is based on a deep convolutional autoencoder which leads to a ROM that lacks interpretability. Designing autoencoders requires careful considerations towards reproducible research including the choice of network size and depth, activation function, optimizer, and weight initialization. Each physical problem usually requires its own hyperparameter tuning. Furthermore, different latent space sizes might lead to different optimal hyperparameter configurations, rendering the training process tedious and often unstable. Solutions to this problem have been proposed in the form of automated parameter choices, but more principled techniques have yet to be established.

In a similar direction, data-driven quadratic manifolds have emerged to build ROMs for transport-dominated problems in [4, 19]. Combining quadratic manifolds with Galerkin projection leads to a highly scalable approach for deriving efficient, yet interpretable, ROMs. However, the application of these approaches to Hamiltonian systems, in their standard formulation, leads to a violation of the underlying geometric structure. This makes the associated ROMs prone to displaying nonphysical behavior (see Section <ref>). The main goal of this work is to build on the recent successes of data-driven quadratic manifolds to construct structure-preserving state approximations, characterized by quadratic nonlinearities, for symplectic model reduction. We focus particularly on high-dimensional dynamical-system models that arise from the semi-discretization of Hamiltonian PDEs. The key contributions of this paper are:

* We propose the quadratic manifold cotangent lift (QMCL) and the resulting SMG-QMCL-ROM in Section <ref> that use a quadratic approximation to build a symplectic approximation mapping while using the SMG as projection. The resulting ROM is a Hamiltonian system of reduced dimension. Thus, the SMG-QMCL is strictly preserving the structure of Hamiltonian systems throughout the reduction.

* We propose the blockwise-quadratic (BQ) approximation mapping and the resulting Galerkin-BQ-ROM in Section <ref> that increases the accuracy of linear-subspace Hamiltonian ROMs by introducing quadratic terms to the state approximation. This approach augments the linear-subspace Hamiltonian ROM with interpretable and physics-preserving higher-order terms to derive an approximately Hamiltonian ROM with improved computational efficiency compared to the SMG-QMCL-ROM.

* Numerical experiments on the challenging cases of a thin moving pulse and a two-dimensional nonlinear wave equation illustrate that the proposed ROMs of the same dimensionality yield more accurate results than the symplectic linear subspace ROMs. In a parameter extrapolation study for the nonlinear example, the proposed methods yield accurate predictions of the high-dimensional state even for parameter values that are outside the range of the training parameters. In a time extrapolation study, the proposed ROMs provide long-time stable and accurate predictions outside the training time interval.

This paper is structured as follows. In Section <ref> we recall the main concepts and definitions of Hamiltonian systems and symplectic model reduction on linear subspaces. We also describe the data-driven construction of quadratic manifolds for nonlinear dimension reduction. Section <ref> presents the proposed SMG-QMCL-ROM and Galerkin-BQ-ROM for deriving ROMs of high-dimensional, Hamiltonian systems. In Section <ref> we present various numerical experiments for the parametric one-dimensional linear wave equation and two-dimensional nonlinear wave equation. Finally, Section <ref> summarizes the contributions and suggests a number of interesting avenues for future research.

§ BACKGROUND

Section <ref> provides an introduction to Hamiltonian PDE models, followed by details about their structure-preserving space discretization. In Section <ref> we review the basic principles of symplectic model reduction for Hamiltonian systems using linear symplectic subspaces. The data-driven learning of quadratic manifolds directly from high-dimensional data is recapitulated in Section <ref>.

§.§ Hamiltonian systems

In the following we consider parametric Hamiltonian systems. Let $\paramSet \subset \real^{n_{\bmu}}$ denote a parameter set for which we denote the $n_{\bmu}$-dimensional elements $\bmu \in \paramSet$ as the parameter vector. We consider infinite-dimensional Hamiltonian systems described by evolutionary PDEs of the form

\begin{equation}
\frac{\partial y(x,t;\bmu)}{\partial t} = {\Ra{\mathcal{S}}} \frac{\delta \mathcal{H}}{\delta y}(y(x,t;\bmu);\bmu),
\label{eq:HPDE}
\end{equation}

where $x$ is the spatial coordinate, $t$ is time, $\mathcal{S}$ is a skew-symmetric operator, and $\delta \mathcal{H}/\delta y$ is the variational derivative[
The variational derivative of $\mathcal{H}$ is defined through
$\frac{\text{d}}{\text{d}\varepsilon}\mathcal{H}[y+ \varepsilon v; \bmu]|_{\varepsilon=0}=\left\langle \frac{\delta \mathcal{H}}{\delta y}(y; \bmu),v\right\rangle$ where $v$ is an arbitrary function.
] of the Hamiltonian energy functional $\mathcal{H}$ with respect to the state variable $y$, and $\bmu \in \paramSet$ is a fixed but arbitrary parameter vector.

Hamiltonian PDEs possess important geometric properties in the form of symplecticity and conservation laws. These geometric properties are intimately related to the ability of the space-discretized FOMs to reproduce the long-time behavior of the solutions of Hamiltonian PDEs [7, 37]. Therefore, the underlying symplectic structure and the conservative nature of Hamiltonian PDEs should be respected in their spatial discretization. Space-discretized Hamiltonian FOMs are finite-dimensional systems that are derived from the Hamiltonian PDE via structure-preserving semi-discretization. Contrary to the traditional approach in which the governing PDEs are discretized directly, Hamiltonian FOMs are obtained through a discretization of the space-time continuous Hamiltonian functional $\mathcal{H}$ and the skew-symmetric operator $\mathcal{S}$ in space [7, 12]. The resulting Hamiltonian FOM is given by

\begin{equation}\label{eq:noncan}
\dot{\by}(t;\bmu)
={\Ra{\mathbf{S}}} \nabla_{\by}H(\by(t;\bmu);\bmu), \qquad \by(0;\bmu)=\by_0(\bmu),
\end{equation}

where $\by(t;\bmu) \in \real^{2n}$ is the high-dimensional state vector, $\by_0:\mathcal{P}\mapsto \real^{2n}$ is the parametric initial condition, ${\Ra{\mathbf{S}=-\mathbf{S}^\top}}$ is the skew-symmetric matrix approximation to ${\Ra{\mathcal{S}}}$ in (<ref>), and $H(\cdot; \bmu) \in \cC^1(\real^{2n})$ is the space-discretized parametric Hamiltonian function. In the above equation, and for the remainder of this paper, the dot notation denotes differentiation with respect to time.

We focus on high-dimensional Hamiltonian systems with canonical symplectic structure ${\Ra{\mathbf S}}=\Jn=\small \begin{pmatrix}
\bzero & \In \\
-\In & \bzero
\end{pmatrix}$, where $\In$ is the $n$-dimensional identity matrix. For canonical Hamiltonian systems, the state vector can be partitioned as $\by(t;\bmu)=(\bq(t;\bmu)^\top,\bp(t;\bmu)^\top)^\top \in \real^{2n}$ where $\bq(t;\bmu) \in \real^n$ is the generalized position vector, and $\bp(t;\bmu) \in \real^n$ is the generalized momentum vector. The governing equation for the canonical Hamiltonian FOM equals

\begin{equation}
\dot{\by}(t;\bmu)=\begin{pmatrix}
\dot{\bq}(t;\bmu) \\
\dot{\bp}(t;\bmu)
\end{pmatrix}=\Jn \nabla_{ \by }H(\bq(t;\bmu),\bp(t;\bmu);\bmu), \qquad \by(0;\bmu)=\by_0(\bmu).
\label{eq:HFOM}
\end{equation}

The flow map[For a dynamical system $\dot{\by}(t)=\mathbf{f}(\by(t))$  described by a differentiable vector field $\bf$ with initial value $\by(0) = \by_0$, the flow map $\varphi_t(\by_0) := \by(t)$ describes the solution as a function of the initial value.] $\varphi_t(\by(0;\bmu);\bmu)=\by(t;\bmu)$ for Hamiltonian FOMs (<ref>) preserves a skew-symmetric, bilinear form known as the canonical symplectic form $\omega = \sum_{i=1}^n \text{d}p_i \wedge \text{d}q_i$ and conserves the Hamiltonian function $\dH$, i.e., $H(\bq(0;\bmu),\bp(0;\bmu);\bmu)=H(\bq(t;\bmu),\bp(t;\bmu);\bmu)$ for all $t$.

The field of geometric numerical integration [27, 37, 57] has developed a wide variety of time integrators for Hamiltonian systems that are designed to respect some of the geometric features, such as energy, momentum, or the symplectic form. Structure-preserving integrators for finite-dimensional Hamiltonian systems can be divided into two categories: (i) symplectic methods [53, 60, 40] that preserve the canonical symplectic form $\omega$ and (ii) energy-preserving integrators [23, 12, 55] that conserve the space-discretized Hamiltonian $\dH$ exactly.

§.§ Symplectic model reduction using linear subspaces

The main goal in symplectic model reduction of high-dimensional Hamiltonian systems is the preservation of symplectic structure during the projection step. The methodologies described in Section <ref> may be viewed as nonlinear extensions of the symplectic model reduction method from [50], which covered linear subspaces and is briefly reviewed here.

To approximate the high-dimensional system state $\by(t;\bmu) \in \real^{2n}$ we seek an approximation mapping $\bqgamma: \real^{2r}\mapsto \real^{2n}$, with $r \ll n$, as follows:

\begin{equation}
\state(t; \bmu) \approx \bqgamma ( \stateRed ( t ; \bmu ) ) \in \real^{2n},
\end{equation}

where $\stateRed(t; \bmu)\in\real^{2r}$ denotes the reduced state vector of dimension $2r$ at time $t$ and parameter vector $\bmu$. We continue by defining the linear symplectic lift (LSL) as

\begin{equation}
\bqgamma_{\text{LSL}} ( \stateRed ( t ; \bmu )) := \by_{\text{ref}} + \bV\stateRed(t; \bmu),
\label{eq:linear_approximation}
\end{equation}

where $\mathbf{y}_\text{ref}\in \real^{2n}$ is a reference state that is used for centering the training data and the basis matrix $\bV \in \real^{2n \times 2r}$ is a symplectic matrix, that is a matrix that satisfies

\begin{equation}
\bV^\top\Jn\bV=\Jr.
\end{equation}

The set of all $2n \times 2r$ symplectic matrices is referred to as the symplectic Stiefel manifold and denoted with $\Sp(2r,\real^{2n})$. The symplectic inverse $\Vt \in \real^{2r \times 2n}$ of a symplectic matrix $\bV$ is defined by

\begin{equation}\label{eq:sympl_inv}
\Vt:=\Jr^\top\bV^\top\Jn
\qquad\text{such that}\qquad
\Vt \bV = \bI_{2r} \in \real^{2r \times 2r},
\end{equation}

which allows one to define the symplectic Galerkin projection with

Assuming $\by(t;\bmu)-\by_{\text{ref}} \in \colspan(\bV)$ for all $t$, the time evolution of the reduced state $\dot{\stateRed}(t;\bmu)$ is derived from the symplectic projection of the FOM (<ref>) on the symplectic subspace via $\Vt$. This results in the ROM equations

\begin{equation} \label{eq:Hr_linear_subs}
\begin{aligned}
\dot{\stateRed}(t;\bmu)&=\Vt\Jn \nabla_{ \by }\dH(\bqgamma_{\text{LSL}} ( \stateRed ( t ; \bmu )) ) \\
&=\Jr\bV^\top\nabla_{\by}\dH(\bqgamma_{\text{LSL}} ( \stateRed ( t ; \bmu )) ) \\
&=\Jr\nabla_{\stateRed}\dH(\bqgamma_{\text{LSL}} ( \stateRed ( t ; \bmu ))),
\end{aligned}
\end{equation}

where we used $\Vt \bV = \bI_{2r}$ on the left-hand side, $\bV^+\Jn=\Jr\bV^\top$ in the second step, and the chain rule
$\nabla_{\stateRed} \dH(\bqgamma_{\text{LSL}} ( \stateRed ( t ; \bmu )))
=\bV^\top \nabla_{\by} \dH(\bqgamma_{\text{LSL}} ( \stateRed ( t ; \bmu )))$ for $\stateRed(t;\bmu) \in \real^{2r}$ in the last step. Even if $\by(t;\bmu)-\by_{\text{ref}} \notin \colspan(\bV)$ for some $t$, the rightmost expression in (<ref>) is still well-defined.

By defining the reduced Hamiltonian $\widetilde{H}: \real^{2r}\mapsto \real$ as $\widetilde H(\stateRed(t;\bmu);\bmu):=\dH (\bqgamma_{\text{LSL}} ( \stateRed ( t ; \bmu )))$, we can rewrite (<ref>) as a canonical Hamiltonian system. Thus, the symplectic Galerkin projection of a $2n$-dimensional Hamiltonian system (<ref>) is given by a $2r$-dimensional LSL-ROM with

\begin{equation}
\dot{\stateRed}(t;\bmu)=\Jr\nabla_{\stateRed}\widetilde{H}(\stateRed(t;\bmu);\bmu).
\label{eq:Hr_linear}
\end{equation}

The basis matrix can be computed using one of many snapshot-based basis generation techniques. Such techniques collect solutions of the Hamiltonian FOM (the snapshots) for different time instances and parameter vectors. The proper symplectic decomposition [50] is a snapshot-based basis generation method to find a symplectic basis matrix $\bV \in \real^{2n\times 2r}$. The basis is required to minimize the projection error of the symplectic projection in the mean over all snapshots which is characterized by the optimization problem

\begin{equation}
\min_{\substack{\bV\in \real^{2n\times 2r} \\ \text{s.t.} \ \bV^\top\Jn\bV=\Jr}} \left\|
\bY_{\bmu} - \bV\Vt \bY_{\bmu}
\right\|_F,
\label{eq:PSD}
\end{equation}

where the centered snapshot data matrix $\bY_{\bmu}$ is defined as

\begin{equation}\label{eq:cotangent_snapshot}
\bY_{\bmu} := \left(
\bY_{\bmu_1}, \dots, \bY_{\bmu_M} \right)
\in \mathbb{R}^{2n \times \ns},
\end{equation}

which stacks $\ns := MK$ centered snapshots in its columns for $M$ parameter instances $\bmu_1,\dots,\bmu_M$, each with $K$ time steps with

\begin{equation}
\bY_{\bmu_j} :=
\begin{pmatrix}
| &  & | \\
(\state(t_1; \bmu_j)-\by_{\text{ref}}) &  \dots & (\state(t_{K}; \bmu_j)-\by_{\text{ref}}) \\
| &  & | \\
\end{pmatrix} \in \mathbb{R}^{2n \times K}.
\end{equation}

For canonical Hamiltonian systems, the centered snapshot data matrix $\bY_{\bmu}$ can be partitioned as $\bY_{\bmu}=(\bQ_{\bmu}^\top,\bP_{\bmu}^\top)^\top$ with the (centered) generalized position and momentum data matrices $\bQ_{\bmu} \in \mathbb{R}^{n \times \ns}$ and $\bP_{\bmu} \in \mathbb{R}^{n \times \ns}$, respectively. We denote the columns of $\bY_{\bmu}$, $\bQ_{\bmu}$,  and $\bP_{\bmu}$ as $\by_j$, $\bq_j$, and $\bp_j$, respectively, with $j=1, \dots, \ns$ in the following.

Cotangent lift algorithm [50]
[1]
Centered snapshot data matrix $\bY_{\bmu}=(\bQ_{\bmu}^\top,\bP_{\bmu}^\top)^\top \in \real^{2n \times \ns}$ arranged as in (<ref>) and reduced dimension $r$
Linear symplectic basis matrix $\bV \in \real^{2n \times 2r}$
Assemble the extended snapshot matrix $\bYe:=(\bQ_{\bmu},\bP_{\bmu})\in \real^{n \times 2\ns}$

Compute the SVD of $\bYe$
$\mathbf{\Phi} \in \mathbb{R}^{n \times r}$ $\leftarrow$ the $r$ leading left-singular vectors of $\bYe$
$\bV= \begin{pmatrix}
\mathbf{\Phi} & \bzero \\
\bzero & \mathbf{\Phi}
\end{pmatrix} \in \real^{2n \times 2r}$ with $\Vt=\bV^\top$ Assemble linear symplectic basis matrix

Since no general solution is known for (<ref>), we seek to find solutions in subsets of $\Sp(2r,\real^{2n})$ by imposing further assumptions on $\bV$. In this paper the cotangent lift algorithm[The cotangent lift algorithm [50] is an SVD-based algorithm for constructing a symplectic basis matrix $\bV$. Unlike the cotangent lift of a smooth map defined for general nonlinear diffeomorphisms [1], the cotangent lift algorithm in [50] only considers linear transformations.], summarized in Algorithm <ref>, is considered for computing a linear symplectic basis

\begin{align}\label{eq:basis_cotan_lift}
\bV = \begin{pmatrix}
\bPhi & \bzero\\
\bzero & \bPhi
\end{pmatrix} \in \real^{2n\times 2r}
\end{align}

of block-diagonal structure. For problems where $\bQ_{\bmu}$ and $\bP_{\bmu}$ have very disparate scales, one can scale the extended snapshot data matrix in Algorithm <ref>, i.e., $\bYe:=(\bQ_{\bmu},\gamma\bP_{\bmu})\in \real^{n \times 2\ns}$, before computing $\bPhi$ via SVD of the extended snapshot matrix, see [50] for more details. Other symplectic basis generation techniques consider greedy basis generation [38, 11], non-orthogonal, symplectic basis generation via the SVD-like decomposition [8], iterative basis generation via optimization on manifolds [6], or optimal symplectic basis generation under certain assumptions on the Hamiltonian system [9].

§.§ Data-driven learning of quadratic state approximations

Many datasets in problems involving high-dimensional state spaces are amenable to dimension reduction using linear techniques such as the POD. However, driving the projection error to an acceptable value can require a large reduced dimension $r$, rendering the associated ROMs inefficient. We therefore follow [19] in introducing nonlinear state approximations with quadratic dependence on the reduced state vector as

\begin{equation}
\state(t; \bmu) \approx \bgammaQ ( \stateRed ( t ; \bmu )) := \by_{\text{ref}}+ \bV \stateRed(t; \bmu) + \overline{\bV} (\stateRed(t; \bmu) \otimes \stateRed(t; \bmu)),
\label{eq:quadratic_approximation}
\end{equation}

where $\bV\in \real^{2n \times 2r}$ is a basis matrix that has $\mathbf{v}_j$ as its $j$th column; $\overline{\bV} \in \real^{2n \times r(2r+1)}$ is a matrix whose columns are populated by the vectors $\overline{\mathbf{v}}_j$; and $\otimes$ denotes the Kronecker product without the redundant terms[For a column vector $\bx=[x_1,x_2, \dots, x_m]^\top \in \real^m$, the column-wise Kronecker product without the redundant terms is defined by
$ \bx \otimes \bx := \begin{bmatrix}
x_{1}^{2} & x_{1}x_{2} & \dots & x_1x_m  & x_{2}^{2} & x_2x_3 & \dots & x_2x_m & \dots & x_m^2\end{bmatrix}^{\top} \in \real^{m(m+1)/2}.$]. Learning the quadratic state approximation (<ref>) from the centered snapshot set $\{ \state_j \}_{j=1}^{\ns} \subset \real^{2n}$ amounts to a representation learning problem in which the matrices $\bV$, $\overline{\bV}$, and the reduced-order state representations $\{ \stateRed_j \}_{j=1}^{\ns} \subset \real^{2r}$ are determined numerically through the solution of the optimization problem

\begin{equation}
\argmin_{\bV, \overline{\bV}, \{ \stateRed_j \}_{j=1}^{\ns} } \left( J(\bV,\overline{\bV},\{ \stateRed_j \}_{j=1}^{\ns}) + \gamma \left\| \overline{\bV} \right\|_F^2 \right),
\label{eq:regularized}
\end{equation}

\begin{equation}
J(\bV,\overline{\bV},\{ \stateRed_j \}_{j=1}^{\ns}) = \sum_{j=1}^{\ns} \left\| \mathbf{y}_j  - \bV \stateRed_j - \overline{\bV} ( \stateRed_j  \otimes \stateRed_j ) \right\|_2^2.
\label{eq:optim_problem2}
\end{equation}

A scalar regularization factor $\gamma \geq 0$ is used to avoid the overfitting of data.

We learn state approximations of the form (<ref>) in a two-step fashion. First, the columns of the basis matrix $\bV$ are chosen as the left-singular vectors corresponding to the $2r$ largest singular values of the centered snapshot data matrix $\bY_{\bmu}$. The representation of the data in the reduced coordinate system can then be computed via the projection $\stateRed_j = \bV^\top\mathbf{y}_j$ for every centered snapshot $\mathbf{y}_j$ with $j=1,\dots,\ns$. This leaves only the matrix operator $\overline{\bV}$ to be inferred from the data. Given the basis matrix $\mathbf{V}$ and snapshot set $\{ \stateRed_j \}_{j=1}^{\ns}$, the optimization problem (<ref>) simplifies to

\begin{equation}
\argmin_{\overline{\bV} } \left( J(\bV,\overline{\bV},\{ \stateRed_j \}_{j=1}^{\ns}) + \gamma \left\| \overline{\bV} \right\|_F^2 \right),
\label{eq:regularized2}
\end{equation}

which is a linear least-squares problem with the explicit solution

\begin{equation}\label{eq:def_Vbar}
\overline{\mathbf{V}} = (\mathbf{I} - \mathbf{V}\mathbf{V}^\top) \bY_{\text{\bmu}} \mathbf{W}^\top ( \mathbf{W} \mathbf{W}^\top + \gamma \mathbf{I})^{-1} \in \mathbb{R}^{2n \times r(2r+1)},
\end{equation}

\begin{equation}
\mathbf{W} :=
\begin{pmatrix}
| &  & | \\
\stateRed_1 \otimes \stateRed_1 & \dots & \stateRed_{\ns} \otimes \stateRed_{\ns}  \\
| &  & |
\end{pmatrix} \in \mathbb{R}^{r(2r+1) \times \ns}.
\label{eq:quadratic_data}
\end{equation}

By construction of $\Vbar$ in (<ref>), each column of $\overline{\mathbf{V}}$ is in the column space of the orthogonal complement of $\mathbf{V}$ in $\real^{2n}$ so that the orthogonality condition $\bV^\top\overline{\mathbf{V}}=\bzero$ holds, see [19] for more details. The inferred basis matrix then approximates $\state(t; \bmu) \approx \bgammaQ ( \stateRed ( t ; \bmu ))$ in a least-squares sense. Equation (<ref>) amounts to a regularized linear least-squares problem that can be solved for each row of $\overline{\bV}$ for a total of $2n$ such problems, each inferring $r(2r+1)$ unique entries.

§ SYMPLECTIC MODEL REDUCTION OF HAMILTONIAN
SYSTEMS USING QUADRATIC MANIFOLDS

We now present two different approaches that use, at their core, data-driven quadratic manifolds for building projection-based ROMs of Hamiltonian systems. Section <ref> motivates the need for preserving the symplectic structure while deriving ROMs when  using quadratic manifold constructions. Section <ref> formalizes the problem considered in this paper. In Section <ref>, we employ quadratic manifolds to develop a novel nonlinear mapping that is guaranteed to be symplectic and then use this mapping to derive a Hamiltonian ROM of (<ref>) via the SMG projection. In Section <ref>, we use data-driven quadratic manifolds to derive approximately Hamiltonian ROMs that are computationally more efficient than the Hamiltonian ROMs derived in Section <ref>. In this alternate approach, we augment the linear symplectic approximation mapping (based on the cotangent lift) with terms of the quadratic Kronecker
product to derive a quadratic approximation mapping and then combine it with a Galerkin projection to derive an approximately Hamiltonian reduced model of (<ref>).

§.§ A motivating example

We motivate the need for symplectic model reduction using quadratic manifolds by demonstrating how the conventional approach for deriving ROMs using quadratic manifolds [19] violates the
underlying Hamiltonian structure. To demonstrate this we consider ROMs for a nonparametric linear wave equation with periodic boundary conditions. The space-time continuous Hamiltonian for the linear wave equation, see definition (<ref>), is discretized using $n=2048$ grid points. This leads to snapshot vectors with $2n=4096$ entries. The scalar parameter value in this numerical experiment is fixed to $\mu=0.5$. Since we aim to reproduce the dynamics of only a single trajectory, we omit the dependence on the scalar parameter $\mu$ in our derivation of the governing ROM equations in this subsection. The corresponding FOM equations for the linear wave equation are

\begin{equation}
\dot{\by}(t)=\begin{pmatrix}
\dot{\bq}(t) \\
\dot{\bp}(t)
\end{pmatrix}=\bA
\begin{pmatrix}
\bq(t) \\
\bp(t)
\end{pmatrix}, \quad \text{with } \bA := \begin{pmatrix}
\bzero & \In \\
0.25\mathbf{D}_{\text{fd}} & \bzero
\end{pmatrix},
\end{equation}

where $\mathbf{D}_{\text{fd}}$ denotes the finite difference approximation of the spatial derivative $\partial_{xx}$. The Hamiltonian FOM is integrated using a symplectic integrator based on the implicit midpoint rule with a fixed time step of $\Delta t=10^{-3}$. For this numerical example, we center the training data about the initial condition (<ref>), that is $\by_{\text{ref}}=\by(0)$. Based on $\ns=K=4000$ FOM data snapshots in the range $t \in [0,4]$, we compute a POD basis matrix $\bV$ by stacking the leading $2r$ left-singular vectors of the centered snapshot data matrix. We then obtain $\Vbar$ by numerically solving the optimization problem (<ref>) via (<ref>). Employing the quadratic manifold approximation of the full state $\state(t)$, together with a Galerkin projection step, yields the following ROM equations

\begin{equation}
\dot{\stateRed}(t)
=\bV^\top \bA\by_{\text{ref}}+\bV^\top \bA\bV \stateRed(t) + \bV^\top \bA\Vbar \left( \stateRed(t) \otimes \stateRed(t) \right).
\label{eq:motivation-rom}
\end{equation}

We then solve (<ref>) using the same fixed time step as the FOM, after which we carry out a reconstruction of the reduced-state data in the original state space using (<ref>).

[subfigure]oneside,margin=1.8cm,0 cm

State error (training data)

FOM energy error

One-dimensional linear wave equation. Even though plot (a) shows low relative state error (<ref>) for some of the quadratic ROMs in the training data regime, the corresponding energy error (<ref>) behavior in plot (b) reveals that the conventional model reduction using quadratic manifolds violates the underlying Hamiltonian structure. The black dashed line indicates end of training time interval.

Figure <ref> shows the state approximation error in the training data regime for ROMs of different reduced dimensions. We observe that the relative state error (see (<ref>) for a precise definition) in the training data regime decreases from $2r=12$ to $2r=20$ and then increases from $2r=20$ to $2r=36$. For $2r=8$ and $2r=40$, the ROM solutions become unstable in the training data regime, which is the reason there are no markers for these points in Figure <ref>. For $2r=32$ and $2r=36$, we observe similar instabilities at approximately $t=8$ (not shown here). Despite some of the quadratic ROMs (<ref>) exhibiting a relative state error below $10^{-1}$ in the training data regime, they do not conserve the space-discretized Hamiltonian $\dH$. The Hamiltonian error plots (see (<ref>) for a precise definition) in Figure <ref> show that, as expected, the FOM energy is not preserved by standard quadratic ROMs of different sizes. This behavior can be attributed to the reduced quadratic operators not preserving the Hamiltonian structure. In the remainder of this section, we propose two novel approaches that employ data-driven quadratic manifolds to derive structure-preserving ROMs of Hamiltonian systems. The numerical experiments in Section <ref> show that the proposed methods lead to accurate and stable ROMs with predictive capabilities.

§.§ Problem formulation

We consider a canonical high-dimensional parametric Hamiltonian system

\begin{equation}
\dot{\by}(t;\bmu)=\Jn\nabla_{\by}\dH(\by(t;\bmu);\bmu), \qquad \by(0;\bmu)=\by_0(\bmu).
\label{eq:FOM_general}
\end{equation}

Specifically, we focus on parametric Hamiltonian systems with FOM Hamiltonians of the form

\begin{equation}\label{eq:quad_sep_Ham}
\dH(\bq,\bp; \bmu)
= \frac{1}{2} \rT\bp \bH_p(\bmu) \bp
+ \frac{1}{2} \rT\bq \bH_q(\bmu) \bq,
\end{equation}

where $\bH_q(\bmu), \bH_p(\bmu) \in \real^{n \times n}$ are obtained by structure-preserving spatial discretization of parametric Hamiltonian PDEs. The corresponding FOM equations are

\begin{equation}
\dot{\by}(t; \bmu)
\dot{\bq}(t; \bmu) \\ \dot{\bp}(t; \bmu)
\end{pmatrix}
=\Jn\nabla_{\by}\dH(\by(t; \bmu); \bmu)
\bzero & \In \\ -\In & \bzero
\end{pmatrix}
\begin{pmatrix}
\bH_q(\bmu) & \bzero\\
\bzero & \bH_p(\bmu)
\end{pmatrix}
\begin{pmatrix}
\bq(t; \bmu) \\
\bp(t; \bmu)
\end{pmatrix}.
\label{eq:lin_sep_Ham_sys}
\end{equation}

The block-diagonal structure in the FOM equations is related to the separable nature of (<ref>) where the FOM Hamiltonian is additively separable with respect to generalized position vector $\bq$ and  generalized momentum vector $\bp$. Separable Hamiltonian FOMs of the form (<ref>) appear as models in many science and engineering applications such as the linear wave equation in physics, the linear elasticity equation in solid mechanics, and Maxwell's equations in electromagnetics.

We consider an intrusive model reduction setting in this work where the FOM operators $\bH_q(\bmu)$ and $\bH_p(\bmu)$ are assumed to be known.

In the scope of this work, a projection-based model reduction technique consists of (a) an approximation mapping $\bqgamma: \real^{2r} \mapsto \real^{2n}$ and (b) a projection of the residual

\begin{equation}
\br(t; \bmu) :=  \frac{\text{d}}{\text{d}t}\left(\bqgamma(\stateRed(t;\bmu))\right) -  \Jn \nabla_{ \by }\dH(\bqgamma(\stateRed(t;\bmu));\bmu ) \in \real^{2n},
\label{eq:mor_residual}
\end{equation}

where $\stateRed(t; \bmu) \in \real^{2r}$ denotes the reduced state vector.

In this paper, we propose two novel approaches that employ data-driven quadratic manifolds, as introduced in Section <ref>, to derive accurate and stable ROMs for Hamiltonian systems. The two distinct model reduction techniques proposed in Sections <ref> and <ref> are named by the approximation map and the projection that are used to define it, SMG-QMCL (with the SMG projection and the QMCL approximation map) and Galerkin-BQ (with a Galerkin projection and the BQ approximation map).

The reduced-order matrix operators in both the SMG-QMCL-ROM and the Galerkin-BQ-ROM may be precomputed in the offline part of the approach.

Since the parameter vectors $\bmu \in \cP$ are fixed (but arbitrary) for each FOM or ROM, we simplify the notation by omitting the explicit dependence on the parameter: the FOM state vector $\by(t;\bmu)$ at time $t$ and parameter $\bmu$ is therefore often denoted as $\by(t)$. Analogously, the generalized position and momentum vectors $\bq(t;\bmu), \bp(t;\bmu)$ are denoted as $\bq(t)$ and $\bp(t)$, respectively.

For noncanonical Hamiltonian systems with a constant and non-degenrate skew-symmetric matrix $\mathbf{S}$ in (<ref>), there exists a congruent transformation which can transform the noncanonical Hamiltonian system to the canonical form, see [50, 25] for more details about the model reduction of noncanonical Hamiltonian systems. Thus, the proposed approaches can be applied to noncanonical systems with constant and non-dengerate $\mathbf{S}$ after transforming the noncanonical FOM to a canonical form.

§.§ Quadratic manifold cotangent lift (QMCL) state approximation and the SMG-QMCL-ROM

We introduce a model reduction technique that strictly preserves the structure of the Hamiltonian system.
To this end, an approximation mapping is constructed—based on data-driven quadratic manifolds—which is guaranteed to be symplectic. We begin by considering a quadratic approximation mapping for the generalized position vector $\stateqRed \in \real^r$ in the form of

\begin{align}\label{eq:bqgammaQMCL}
\bqgammaQMCL(\stateqRed)
:= \qref + \Vq \stateqRed + \Vqbar (\stateqRed \otimes \stateqRed),
\end{align}

with a constant vector $\qref \in \real^{n}$ and matrices $\Vq \in \real^{n \times r}$, $\Vqbar \in \real^{n \times r(r+1)/2}$ with

\begin{align}\label{eq:bqgammaQMCL_emb_matrices}
\rT\Vq \Vq =& \Ir
&\rT\Vq \Vqbar = \bzero \in \real^{r \times r(r+1)/2}.
\end{align}

Based thereon, we introduce the quadratic manifold cotangent lift (QMCL) as the approximation mapping of the state $\stateRed = (\stateqRed^\top, \statepRed^\top)^\top \in \real^{2r}$ with

\begin{align}\label{eq:quad_mnf_co_tan}
\bgammaQMCL(\stateqRed, \statepRed)
:= \begin{pmatrix}
\bqgammaQMCL(\stateqRed)\\
\bpgammaQMCL(\stateqRed, \statepRed)
\end{pmatrix}
:= \begin{pmatrix}
\qref + \Vq \stateqRed + \Vqbar (\stateqRed \otimes \stateqRed)\\
\pref + \Big( \Vq + \Vqbar \kronFun(\stateqRed) \Big) \vRed(\stateqRed, \statepRed)
\end{pmatrix},
\end{align}

with $\vRed(\stateqRed, \statepRed)
:= \invb{
\Ir + \rTb{\kronFun(\stateqRed)} \rT\Vqbar \Vqbar \kronFun(\stateqRed)
\statepRed \in \real^r$ where $\kronFun(\stateqRed) := \eval[\stateqRed]{\bD (\cdot \otimes \cdot)} \in \real^{r(r+1)/2 \times r}$ denotes the derivative of the Kronecker product and a constant vector $\pref \in \real^{n}$. We show later in this subsection (in a more general framework) that the approximation mapping $\bgammaQMCL$ is (i) continuously differentiable and (ii) a symplectic map, i.e. the Jacobian $\eval[\stateRed]{\bD \bgammaQMCL} \in \real^{2n \times 2r}$ is a symplectic matrix

\begin{align}\label{eq:sympl_map}
\rT{\eval[\stateRed]{\bD \bgammaQMCL}} \Jn \eval[\stateRed]{\bD \bgammaQMCL} = \Jr
\qquad \text{for all } \stateRed \in \real^{2r}.
\end{align}

These two properties allow us to formulate the symplectic inverse (<ref>) of the Jacobian $\eval[\stateRed]{\bD \bgammaQMCL}$ which we denote as $\eval[{\stateRed(t)}]{\bD \bgammaQMCL}^+$. Thus, we can directly apply the techniques from [10] to use the SMG projection.

The SMG projection requires that the symplectic projection of the residual (<ref>) with the Jacobian of the approximation mapping vanishes, that is

\begin{equation}\label{eq:smg_proj_mor_residual}
\eval[{\stateRed(t)}]{\bD \bgammaQMCL}^+ \br(t) = \bzero \in \real^{2r} \quad  {\rm{for\ all\ }} t.
\end{equation}

The ROM obtained via the SMG projection (<ref>) of a canonical Hamiltonian system with $\bgammaQMCL$ (<ref>) is again a canonical Hamiltonian system.

We apply the proof from <cit.> to the QMCL approximation mapping proposed in the present work: Substituting the residual (<ref>) into the definition of the SMG projection (<ref>) leads to

\begin{align*}
\dot{\stateRed}(t)
=&\eval[{\stateRed(t)}]{\bD \bgammaQMCL}^+ \Jn \nabla_{ \by } \dH(\bgammaQMCL(\stateRed(t)))\\
=&\Jr \eval[{\stateRed(t)}]{\bD \bgammaQMCL}^\top \nabla_{\by} \dH(\bgammaQMCL(\stateRed(t)))\\
\end{align*}

where we use that for all $\stateRed \in \real^{2r}$ it holds that
$\eval[{\stateRed}]{\bD \bgammaQMCL}^+
\eval[{\stateRed}]{\bD \bgammaQMCL}
= \bI_{2r}$ on the left-hand side,
$\eval[{\stateRed}]{\bD \bgammaQMCL}^+ \Jn
= \Jr \eval[{\stateRed}]{\bD \bgammaQMCL}^\top$
in the second step, and the chain rule for the gradient of the Hamiltonian
$\nabla_{\stateRed} \dH(\bgammaQMCL(\stateRed))
= \eval[{\stateRed}]{\bD \bgammaQMCL}^\top \nabla_{\by} \dH(\bgammaQMCL(\stateRed))$ in the last step.
With the definition of the reduced Hamiltonian $\dHRed := \dH \circ \bgammaQMCL \in \cC^1(\real^{2r})$, we see that this is indeed a canonical Hamiltonian system

\begin{align}\label{eq:jac_based_nonlin_Ham}
\dot{\stateRed}(t) = \Jr \nabla_{\stateRed}\dHRed(\stateRed(t)) \in \real^{2r},
\end{align}

of reduced dimension $2r$.

For the special case of a linear Hamiltonian system of the form (<ref>), the ROM obtained using the SMG-QMCL approach can be written as

\begin{align}\label{eq:jac_based_lin_Ham}
\dot{\stateRed}(t)
=& \Jr \rT{ \eval[\stateRed(t)]{\bD \bgammaQMCL} } \bH \bgammaQMCL(\stateRed(t)),&
\bH
\bH_q & \bzero\\
\bzero & \bH_p
\end{pmatrix}
\in \real^{2n \times 2n}.
\end{align}

To write the specific operations out more explicitly, we introduce the abbreviations

\begin{align}
\VV
\Vq                 & \bzero\\
\bzero & \Vq
\end{pmatrix}
\in \real^{2n \times 2r},&
\VbarVbar
&:= \begin{pmatrix}
\Vqbar                 & \bzero\\
\bzero & \Vqbar
\end{pmatrix}
\in \real^{2n \times r(r+1)},&
\yref
&:= \begin{pmatrix}
\qref\\
\pref
\end{pmatrix}
\in \real^{2n},
\end{align}

and for $\stateRed = (\stateqRed^\top, \statepRed^\top)^\top \in \real^{2r}$ the respective nonlinear, but low-dimensional, coefficient functions

\begin{align*}
\stateV(\stateqRed, \statepRed)
&:= \begin{pmatrix}
\stateqRed\\
\vRed(\stateqRed, \statepRed)\\
\end{pmatrix}
\in \real^{2r},&
\stateVbar(\stateqRed, \statepRed)
&:= \begin{pmatrix}
(\stateqRed \otimes \stateqRed)\\
\kronFun(\stateqRed)
\vRed(\stateqRed, \statepRed)
\end{pmatrix}
\in \real^{r(r+1)},
\end{align*}

such that the mapping $\bgammaQMCL$ and its derivative are

\begin{align*}
\bgammaQMCL(\stateRed)
\yref
+ \VV \stateV(\stateRed)
+ \VbarVbar \stateVbar(\stateRed)
\eval[\stateRed]{\bD\bgammaQMCL}
=& \VV \eval[\stateRed]{\bD\stateV}
+ \VbarVbar \eval[\stateRed]{\bD\stateVbar}.
\end{align*}

The SMG-QMCL-ROM for the linear Hamiltonian system (<ref>) can be written as

\begin{align}\label{eq:jac_based_lin_Ham_quad_gamma}
\begin{split}
\dot{\stateRed}(t)
=&\; \Jr
\rT{\eval[\stateRed(t)]{\bD\stateV}}
\left(
\bHVyref
\bHVV\;
\stateV(\stateRed(t))
\bHVVbar\;
\stateVbar(\stateRed(t))
\right)
\\
& + \Jr
\rT{\eval[\stateRed(t)]{\bD\stateVbar}}
\left(
\bHVbaryref
\rT{\bHVVbar}\;
\stateV(\stateRed(t))
\bHVbarVbar\;
\stateVbar(\stateRed(t))
\right),
\end{split}
\end{align}

with different projections of $\bH$ as

\begin{gather}
\bHVyref
:= \rT\VV \bH \yref \in \real^{2r},\quad
\bHVV
:= \rT\VV \bH \VV \in \real^{2r \times 2r},\quad
\bHVVbar
:= \rT\VV \bH \VbarVbar \in\real^{2r \times r(r+1)},\nonumber\\
\label{eq:jac_based_proj_ops}
\bHVbaryref
:= \rT\VbarVbar \bH \yref \in\real^{r(r+1)},\quad
\bHVbarVbar
:= \rT\VbarVbar \bH \VbarVbar \in\real^{r(r+1) \times r(r+1)}.
\end{gather}

Note that all multiplications in (<ref>) are independent of the full dimension $n$ and only depend on the reduced dimension $r$. Thus, this model is offline–online separable if $\bH$ is parameter-separable as then the reduced operators (<ref>) can be computed offline following standard techniques (see e.g. <cit.>). In the online stage, the SMG-QMCL approach requires $\bigO(r^5)$ multiplication operations to evaluate the right-hand side of the nonlinear ROM and $\bigO(r^5)$ multiplication operations to evaluate the Jacobian of this right-hand side.

The derived equations hold for all $\Vq \in \real^{n \times r}$ and $\Vbar \in \real^{n \times r(r+1)/2}$ which fulfill (<ref>). In our approach, we compute $\Vq=\bPhi$ through the cotangent lift algorithm from Algorithm <ref> and the matrix $\Vqbar$ for the quadratic part is computed as in <Ref> by only using the snapshot matrix $\bQ_{\bmu}$. This approach minimizes the reconstruction errors in $\bq$ based on $\bQ_{\bmu}$ snapshot data, which leads to a linear least-squares problem for computing $\Vqbar$. Alternatively, one could think about construction $\Vqbar$ based on reconstruction errors for $\bq$ and $\bp$. However, since $\bpgammaQMCL$ depends on the Jacobian of $\bqgammaQMCL$, this would lead to challanging nonlinear least-squares problem. The SMG-QMCL approach is summarized in Algorithm <ref>.

Since the SMG-QMCL approach is based on the SMG projection,
we can directly transfer results on the stability in the sense of Lyapunov:

If the Hamiltonian FOM (<ref>) has Lyapunov-stable states and these states are included in the image of the symplectic approximation mapping $\bgammaQMCL$ (<ref>), then these states are also Lyapunov-stable states of the SMG-QMCL-ROM (<ref>).

It remains to show that the proposed nonlinear approximation map $\bgammaQMCL$ is (i) continuously differentiable and (ii) a symplectic map. The proof is based on a more general approximation mapping, which we refer to as the manifold cotangent lift. We begin by introducing immersions and the Moore–Penrose inverse:

An immersion $\bf \in \cC^1(\real^k, \real^m)$ is a mapping for which the Jacobian $\eval[\bx]{\bD \bf} \in \real^{m \times k}$ is of rank $k$ at every point $\bx \in \real^m$.

For an immersion $\bf \in \cC^1(\real^k, \real^m)$ with $k \leq m$,
the Moore–Penrose inverse of the Jacobian

\begin{equation}\label{eq:mpinv}
\mpinv{ \eval[\bx]{\bD\bf} }
:= \invb{
\rT{\eval[\bx]{\bD\bf}}
\eval[\bx]{\bD\bf}
\rT{\eval[\bx]{\bD\bf}}
\end{equation}

is well-defined for all $\bx \in \real^k$.
If $\bf \in \cC^2(\real^k, \real^m)$,
we have that $\mpinv{ \eval[(\cdot)]{\bD\bf} } \in \cC^1(\real^k, \real^{m \times k})$.

Since $\bf \in \cC^1(\real^k, \real^m)$ is an immersion with $k \leq m$,
the Jacobian $\eval[\bx]{\bD \bf} \in \real^{m \times k}$ is of full column rank $k$ for all $\bx \in \real^k$.
Thus, the product $\rTb{\eval[\bx]{\bD \bf}} \eval[\bx]{\bD \bf}$ is of full rank and thus (<ref>) is well-defined.
For an immersion $\bf \in \cC^2(\real^k, \real^m)$,
we know that $\eval[(\cdot)]{\bD \bf} \in \cC^1(\real^k, \real^{m \times k})$ with constant rank $k$ which is enough to show $\mpinv{ \eval[(\cdot)]{\bD\bf} } \in \cC^1(\real^k, \real^{m \times k})$ (see <cit.>).

We propose a nonlinear generalization of the linear symplectic lift employed in the cotangent lift algorithm, which we refer to as the manifold cotangent lift.

For a given immersion $\bqgammaMCL \in \cC^2(\real^r, \real^n)$,
we define the manifold cotangent lift (MCL) embedding as

\begin{align}\label{eq:map_mnf_co_tan}
\bgammaMCL(\stateqRed, \statepRed)
:=& \begin{pmatrix}
\bqgammaMCL(\stateqRed)\\
\bpgammaMCL(\stateqRed, \statepRed)
\end{pmatrix}&
&\bpgammaMCL(\stateqRed, \statepRed)
:= \pref + \rTb{\mpinv{ \eval[\stateqRed]{\bD\bqgammaMCL} }} \statepRed,
\end{align}

for some $\pref \in \real^{r}$.

We show that an MCL embedding is (i) continuously differentiable and (ii) a symplectic map.

The MCL embedding $\bgammaMCL$ from (<ref>) is continuously differentiable.

We argue about the differentiability of the
two components of $\bgammaMCL$ (<ref>) in three steps.
First, we know that the first component $\bqgammaMCL \in \cC^2(\real^{r}, \real^n)$ is continuously differentiable by assumption.
Second, $\bpgammaMCL$ is smooth in $\statepRed$ since the function is linear in $\statepRed$.
Finally, since $\bqgammaMCL \in \cC^2(\real^r, \real^n)$ is an immersion, <Ref> applies and we know that the Moore–Penrose pseudo inverse is continuously differentiable.
Then, we know that $\bpgammaMCL$ is continuously differentiable in $\stateqRed$ by composition of continuously differentiable functions.

The MCL embedding $\bgammaMCL$ from (<ref>) is a symplectic map (<ref>).

For all $\stateRed = (\stateqRed^\top, \statepRed^\top)^\top \in \real^{2r}$, the derivative of the nonlinear mapping $\bgammaMCL$ can be written as

\begin{align*}
\eval[\stateRed]{ \bD \bgammaMCL }
= \begin{pmatrix}
\eval[\stateqRed]{\bD \bqgammaMCL} & \bzero_{n \times r}\\
\eval[\stateRed]{\Dq \bpgammaMCL}
& \rTb{\mpinv{ \eval[\stateqRed]{\bD\bqgammaMCL} }}
\end{pmatrix}
\in \real^{2n \times 2r},
\end{align*}

where $\bzero_{n \times r} \in \real^{n \times r}$ is the matrix of all zeros
and $\Dq(\cdot)$ denotes the partial derivative w.r.t. $\stateqRed$.
Using the above expression for $\eval[\stateRed]{ \bD \bgammaMCL }$, we can write the left-hand side of (<ref>) as

\begin{align}\label{eq:proof_sympl_intermed_term}
\rT{\eval[\stateRed]{ \bD \bgammaMCL }} \Jn \eval[\stateRed]{ \bD \bgammaMCL }
= \begin{pmatrix}
2 \skewPart \left(
\rT{\eval[\stateqRed]{\bD \bqgammaMCL}}
\eval[\stateRed]{\Dq \bpgammaMCL}
\right)
& \rT{\eval[\stateqRed]{\bD \bqgammaMCL}}
\rTb{\mpinv{ \eval[\stateqRed]{\bD\bqgammaMCL} }}\\
-\mpinv{ \eval[\stateqRed]{\bD\bqgammaMCL} }
\eval[\stateqRed]{\bD \bqgammaMCL}
& \bzero_{r \times r}
\end{pmatrix},
\end{align}

where $\skewPart(\bA) := \frac{1}{2} (\bA - \rT\bA)$ denotes the skew-symmetric part of $\bA \in \real^{r \times r}$.
To prove the symplecticity, we need to show that (<ref>) is equal to $\Jr$. We observe that the off-diagonal terms already match due to

= \mpinv{ \eval[\stateqRed]{\bD\bqgammaMCL} } \eval[\stateqRed]{\bD \bqgammaMCL}
= \rTb{\mpinv{ \eval[\stateqRed]{\bD\bqgammaMCL} } \eval[\stateqRed]{\bD \bqgammaMCL}}
= \rT{\eval[\stateqRed]{\bD \bqgammaMCL}} \rTb{\mpinv{ \eval[\stateqRed]{\bD\bqgammaMCL} }}.$$

Thus, it remains to show that the term with the skew-symmetric part in (<ref>) equals zero for which we use that

\begin{align}\label{eq:proof_sympl_statepRed_term}
\statepRed
= \bI_{r} \statepRed
= \rT{\eval[\stateqRed]{\bD \bqgammaMCL}}
\rTb{\mpinv{ \eval[\stateqRed]{\bD\bqgammaMCL} }} \statepRed
= \rT{\eval[\stateqRed]{\bD \bqgammaMCL}} \bpgammaMCL(\stateqRed, \statepRed).
\end{align}

For a better presentation, we reformulate this equation to index notation and use the notation
$(\bqgammaMCL)_{i,j} = \left( \eval[\stateqRed]{\bD \bqgammaMCL} \right)_{ij}$ for the Jacobian with respect to $\stateqRed$.
Then, we can write (<ref>) in index notation for $1\leq i \leq r$ as

\begin{align*}
= \sum_{k=1}^r (\bqgammaMCL)_{k,i}
\left( \bpgammaMCL \right)_k.
\end{align*}

Deriving this equation w.r.t. $\stateqRed$ yields by using the product rule for $1\leq i,j \leq r$

\begin{align*}
= \sum_{k=1}^r (\bqgammaMCL)_{k,ij}
\left( \bpgammaMCL \right)_k
+ \underbrace{
\sum_{k=1}^r
(\bqgammaMCL)_{k,i} \left(
\bpgammaMCL
\right)_{k,j}
= \left(
\rTb{\eval[\stateqRed]{\bD \bqgammaMCL}}
\eval[\stateRed]{\Dq \bpgammaMCL}
\right)_{ij}
\end{align*}

where we identify in the underbrace the term which we want to make an assertion about for (<ref>).
We move this term to the left-hand side and obtain for $1 \leq i,j \leq r$

\begin{align}\label{eq:proof_sympl_sym_term}
\left(
\rT{\eval[\stateqRed]{\bD \bqgammaMCL}}
\eval[\stateRed]{\Dq \bpgammaMCL}
\right)_{ij}
= - \sum_{k=1}^r (\bqgammaMCL)_{k,ij} \left( \bpgammaMCL \right)_k.
\end{align}

Since $\bqgammaMCL \in \cC^2(\real^{2r}, \real^{2n})$, we know by the Lemma of Schwarz that the Hessian of each component $(\bqgammaMCL)_i$ is symmetric, i.e. in the presented index notation $(\bqgammaMCL)_{k,ij} = (\bqgammaMCL)_{k,ji}$.
Thus we know that the terms in (<ref>) are symmetric in $i$ and $j$ and the skew-symmetric part vanishes.
This concludes the proof that (<ref>) equals $\Jr$ and thus $\bgammaMCL$ is a symplectic map.

We can show that  the cotangent lift and the QMCL can be interpreted as an MCL embedding. For this purpose, we have to show (a) that these methods can be formulated with (<ref>) via a specific choice for $\bqgammaMCL$ and (b) that this specific $\bqgammaMCL$ is indeed an immersion.

For a matrix $\Vq \in \real^{n \times r}$ with orthonormal columns, $\rT\Vq \Vq = \bI_r$,
the choice $\bqgammaMCL(\stateqRed) ={\Ra{\qref +}} \Vq \stateqRed$ in the MCL embedding (<ref>) recovers the cotangent lift with $\Vq = \bPhi$.

The Jacobian $\bD\bqgammaMCL \equiv \Vq$ is of full column rank by the assumption of orthonormal columns in $\Vq$ and thus this linear choice for $\bqgammaMCL$ is an immersion.
Since $\bD\bqgammaMCL \equiv \Vq$ with $\rT\Vq \Vq = \bI_{r}$, we know that $\bpgammaMCL(\stateqRed, \statepRed) ={\Ra{ \pref +}} \invb{\rT\Vq \Vq} \Vq \statepRed ={\Ra{\pref +}}  \Vq \statepRed$.
Thus, the approximation mapping from (<ref>) equals the approximation mapping in the linear symplectic lift (<ref>) with $\bV$ from (<ref>).

For the quadratic function $\bqgammaMCL \equiv \bqgammaQMCL$ from (<ref>), the MCL recovers the QMCL from the beginning of this section.

With the requirements (<ref>) on $\Vq$ and $\Vqbar$, we know

\begin{align*}
\rT{\eval[\stateqRed]{\bD\bqgammaMCL}} \eval[\stateqRed]{\bD\bqgammaMCL}
= \Ir + \rTb{\kronFun(\stateqRed)} \rT\Vqbar \Vqbar \kronFun(\stateqRed),
\end{align*}

and thus $\bpgammaMCL(\stateqRed, \statepRed) = \bpgammaQMCL(\stateqRed, \statepRed)$.
Again, using (<ref>), we can show that $\bqgammaQMCL$ is an immersion since

\begin{align*}
&\eval[\stateqRed]{\bD\bqgammaMCL} = \Vq + \Vqbar \kronFun(\stateqRed)&
&\text{and thus}&
&\rT\Vq \eval[\stateqRed]{\bD\bqgammaMCL} = \Ir.
\end{align*}

Then, with a rank argument $\eval[\stateqRed]{\bD\bqgammaMCL}$ is of full column rank $r$ for all $\stateqRed \in \real^r$.

MCL embeddings offer a blueprint for generating (nonlinear) mappings which are guaranteed to be symplectic. This is a clear advantage over [10] which uses a weakly symplectic approach in the numerical experiments to determine nonlinear approximations of the state. Moreover, the proposed QMCL approximation mapping allows an offline–online separation for parameter-separable, linear Hamiltonian systems (which is another missing key piece in [10]) for symplectic model reduction using nonlinear approximation mappings.

§.§ Blockwise-quadratic (BQ) state approximation and the Galerkin-BQ-ROM

The SMG-QMCL approach in Section <ref> employs a state-dependent projection (<ref>) which leads to a nonlinear Hamiltonian ROM (<ref>) with computational complexity $\bigO(r^5)$ in the online stage for linear problems. In this section, we propose an alternative approach to derive approximately Hamiltonian ROMs that are computationally more efficient in online computations. This approach also retains the physical interpretation of the state variables at the reduced level by choosing the following quadratic approximation for $\stateRed = (\stateqRed^\top, \statepRed^\top)^\top \in \real^{2r}$

\begin{align}
\label{eq:gammqp}
\bgammaBQ(\stateRed) :=& \begin{pmatrix}
\bqgammaBQ(\stateqRed)\\
\bpgammaBQ(\statepRed)
\end{pmatrix} :=
\begin{pmatrix}
\bq_{\text{ref}} + \bV_q \stateqRed + \Vbar_q ( \stateqRed \otimes \stateqRed) \\
\bp_{\text{ref}} + \bV_p \statepRed + \Vbar_p ( \statepRed \otimes \statepRed )
\end{pmatrix},
\end{align}

which we refer to as blockwise-quadratic (BQ) approximation mapping.
We then choose $\bV_q=\bV_p=\bPhi$ in accord with cotangent lift algorithm <ref>. The matrices $\Vbar_q$ and $\Vbar_p$ are obtained from centered snapshots $\{ \state_j \}_{j=1}^{\ns}$ through the numerical solution of the following pair of optimization problems:

\begin{equation}
\argmin_{\overline{\bV}_q} \left( J_q(\bV_q,\overline{\bV}_q,\{\stateqRed_j\}^{n_s}_{j=1}) + \gamma_q \left\| \overline{\bV}_q \right\|_F^2 \right);
\quad
\argmin_{\overline{\bV}_p} \left( J_p(\bV_p,\overline{\bV}_p,\{\statepRed_j\}^{n_s}_{j=1}) + \gamma_p \left\| \overline{\bV}_p \right\|_F^2 \right),
\label{eq:optim_problem}
\end{equation}

\begin{equation}
\begin{aligned}
J_q(\bV_q,\overline{\bV}_q,\{\stateqRed_j\}^{n_s}_{j=1}) &:= \sum_{j=1}^{\ns} \left\| \bq_j - \bPhi\stateqRed_j -\Vbar_q ( \stateqRed_j  \otimes \stateqRed_j ) \right\|_2^2, \\
J_p(\bV_p,\overline{\bV}_p,\{\statepRed_j\}^{n_s}_{j=1}) &:= \sum_{j=1}^{\ns} \left\| \bp_j - \bPhi\statepRed_j -\Vbar_p ( \statepRed_j  \otimes \statepRed_j ) \right\|_2^2,
\end{aligned}
\end{equation}

with the reduced snapshots $\stateqRed_j = \rT\bPhi \stateq_j$ and $\statepRed_j = \rT\bPhi \statep_j$ from $\state_j = (\stateq_j^\top,\statep_j^\top)^\top \in \real^{2n}$ for $j = 1, \dots, \ns$ and the scalar regularization parameters $\gamma_q,\gamma_p \geq 0$.

Since the BQ state approximation is not a symplectic map, we can not use the SMG projection for deriving ROMs. Instead, we use a Galerkin projection that requires the residual (<ref>) to be orthogonal to the columns of $\bV=  \begin{pmatrix}
\bPhi & \bzero \\
\bzero & \bPhi
\end{pmatrix}$, that is $\bV^\top \br(t; \bmu) = \bzero$.
Fusing BQ approximation (<ref>) with the Galerkin projection, we obtain the Galerkin-BQ-ROM

\begin{equation}\label{pod_based_nonlinear_Ham}
\dot{\stateRed}(t) = \bV^\top\Jn\nabla_{\by} \dH(\bgammaBQ(\stateRed(t)))=\Jr\bV^\top\nabla_{\by}\dH(\bgammaBQ(\stateRed(t))).
\end{equation}

For linear Hamiltonian systems of the form (<ref>), the Galerkin-BQ-ROM simplifies to

\begin{equation}
\begin{aligned}
\dot{\stateqRed}(t)
&= \phantom{-}\bPhi^\top \bH_p \bp_{\text{ref}} + \bPhi^\top \bH_p \bPhi \statepRed(t) + \bPhi^\top\bH_p \Vbar_p \left(     \statepRed(t) \otimes \statepRed(t) \right), \\
\dot{\statepRed}(t)
&=- \bPhi^\top \bH_q \bq_{\text{ref}} - \bPhi^\top \bH_q \bPhi \stateqRed(t) - \bPhi^\top \bH_q \Vbar_q \left( \stateqRed(t) \otimes \stateqRed(t)
\right),
%    \label{eq:phatrom}
\end{aligned}
\end{equation}

where we used the orthogonality properties of the basis matrices $\bPhi^\top \Vbar_q  = \bPhi^\top \Vbar_p = \bzero$. This system of equations can be rewritten as

\begin{equation}\label{eq:pod_Ham_rom}
\begin{pmatrix}
\dot{\stateqRed}(t) \\
\dot{\statepRed}(t)
\end{pmatrix}
= \begin{pmatrix}
\phantom{-}\tbHrefp\\
\end{pmatrix}
\bzero & \Ir \\
-\Ir & \bzero
\end{pmatrix}
\begin{pmatrix}
\tbHq & \bzero \\
\bzero & \tbHp
\end{pmatrix}
\begin{pmatrix}
\stateqRed(t) \\
\statepRed(t)
\end{pmatrix}\\
+ \begin{pmatrix}
\bzero & \tbHpbar  \\
-\tbHqbar & \bzero
\end{pmatrix}  \begin{pmatrix}
\stateqRed(t)  \otimes \stateqRed(t) \\
\statepRed(t)  \otimes \statepRed(t)
\end{pmatrix},
\end{equation}

where the reduced operators are

\begin{align} \label{eq:pod_rom_operators}
\begin{aligned}
\tbHrefq &:= \bPhi^\top \bH_q \qref \in \real^{r},&
\tbHq &:=\bPhi^\top \bH_q \bPhi\in \real^{r\times r},&
\tbHqbar &:=\bPhi^\top \bH_q \Vbar_q\in \real^{r\times r(r+1)/2}, \\
\tbHrefp &:= \bPhi^\top \bH_p \pref\in \real^{r},&
\tbHp &:=\bPhi^\top \bH_p \bPhi\in \real^{r\times r},&
\tbHpbar &:=\bPhi^\top \bH_p \Vbar_p\in \real^{r\times r(r+1)/2}.
\end{aligned}
\end{align}

Approximation (<ref>) of the SMG-QMCL-ROM transforms the high-dimensional problem (<ref>) into a reduced nonlinear Hamiltonian system. In sharp contrast, the use of the Galerkin-BQ-ROM leads to a reduced linear Hamiltonian system with a quadratic perturbation term. Nonetheless, (<ref>) exploits knowledge about the canonical structure of the FOM, leading to an approximately Hamiltonian ROM (<ref>) that retains the intrinsic coupled structure of the FOM and has interpretable states. The numerical results in Section <ref> demonstrate that this is indeed a very accurate strategy. The Galerkin-BQ-ROM is summarized in Algorithm <ref>.

The Galerkin-BQ-ROM requires $\bigO(r^3)$ online multiplication operations to evaluate the right-hand side of the ROM and $\bigO(r^4)$ multiplication operations to evaluate the Jacobian of this right-hand side. The SMG-QMCL-ROM, on the other hand, requires $\bigO(r^5)$ multiplications for both of these tasks. In other words, by allowing an approximate Hamiltonian structure we achieve a computationally more efficient online phase compared to the SMG-QMCL-ROM from Section <ref>. Importantly, we do not compromise the improved accuracy from data-driven quadratic manifolds.

Offline phase for the SMG-QMCL-ROM
[1]
Centered snapshot data matrix $\bY_{\bmu}=(\bQ_{\bmu}^\top,\bP_{\bmu}^\top)^\top \in \real^{2n \times \ns}$ arranged as in (<ref>), reference state $\by_{\text{ref}}=(\bq_{\text{ref}}^\top,\bp_{\text{ref}}^\top)^\top$, Hamiltonian FOM operators $\bH_q$ and $\bH_p$ (<ref>), and reduced dimension $r$.
Reduced operators $\bHVyref$, $\bHVbaryref$, $\bHVV$, $\bHVbarV, \bHVbarVbar$ for Hamiltonian ROM (<ref>) and basis matrices $\Vq$ and $\Vqbar$ (<ref>).
$\Vq = \bPhi$ $\leftarrow$ Compute symplectic basis matrix using Algorithm <ref> Basis computation
$\widetilde{\bq}_j=\bPhi^\top\bq_j$ for $j=1,\dots,\ns$ Represent ${\Rc{\bq}}_j$ in reduced-dim. coordinate system
$\Vqbar$ $\leftarrow$ Solve linear least-squares problem (<ref>) Representation learning problem
$\bHVyref$, $\bHVbaryref$, $\bHVV, \bHVVbar, \bHVbarVbar$ $\leftarrow$ Compute matrix operators from (<ref>) ROM operators

Offline phase for the Galerkin-BQ-ROM
[1]
Centered snapshot data matrix $\bY_{\bmu}=(\bQ_{\bmu}^\top,\bP_{\bmu}^\top)^\top \in \real^{2n \times \ns}$ arranged as in (<ref>), reference state $\by_{\text{ref}}=(\bq_{\text{ref}}^\top,\bp_{\text{ref}}^\top)^\top$, Hamiltonian FOM operators $\bH_q$ and $\bH_p$ (<ref>), and reduced dimension $r$.
Reduced operators $\tbHrefq$, $\tbHrefp$, $\tbHq$, $\tbHp$, $\tbHqbar$, $\tbHpbar$ for Hamiltonian ROM (<ref>), and basis matrices $\bV_q$, $\bV_p$, $\Vbar_q$, and $\Vbar_p$ (<ref>).
$\bV_q=\bV_p=\bPhi$ $\leftarrow$ Compute symplectic basis matrix using Algorithm <ref> Basis computation
$\widetilde{\bq}_j=\bPhi^\top\bq_j$ for $j=1,\dots,\ns$ Represent ${\Rc{\bq}}_j$ in reduced-dim. coordinate system
$\widetilde{\bp}_j=\bPhi^\top\bp_j$ for $j=1,\dots,\ns$ Represent ${\Rc{\bp}}_j$ in reduced-dim. coordinate system
$\Vbar_q, \Vbar_p$ $\leftarrow$ Solve linear least-squares problems (<ref>) Representation learning problem
$\tbHrefq, \tbHrefp, \tbHq, \tbHp, \tbHqbar, \tbHpbar$ $\leftarrow$ Compute matrix operators from (<ref>) ROM operators

§ NUMERICAL RESULTS

In this section, the proposed model reduction methods are applied to two parametrized wave equations, which are prototypical of transport problems for which slowly-decaying Kolmogorov $N$-widths have been observed in similar settings in [24]. Section <ref> provides details about the numerical implementations and the reported error measures. In Section <ref>, we demonstrate that the proposed approaches yield more accurate ROMs than the linear symplectic ROMs for a model that generalizes in the parameter and extrapolates in time. Finally, in Section <ref>, we demonstrate the parameter extrapolation capabilities of the proposed approaches on a two-dimensional nonlinear wave equation.

§.§ Practical considerations & error measures

All numerical experiments in this paper are conducted using MATLAB version 2022a. For time integration we use the implicit midpoint rule for all FOM and ROM simulations. For time-continuous dynamical systems the corresponding time-marching equations are

\begin{equation*}
\frac{\by_{k+1}-\by_k}{\Delta t}=\bf\left(\frac{\by_k + \by_{k+1}}{2}\right),
\end{equation*}

where $\Delta t$ is the fixed time step. The implicit midpoint rule is a second-order symplectic integrator that exhibits bounded energy error for nonlinear Hamiltonian systems [27]. The Galerkin-BQ-ROM (<ref>) is solved at every time step using Newton's method. The SMG-QMCL-ROM (<ref>), on the other hand, uses a quasi-Newton scheme  that neglects the second-order derivatives of $\bgammaQMCL$ for an improvement of the online runtime, following the procedures from [36].

The regularization parameters $\gamma_q, \gamma_p$ in (<ref>) are chosen to ensure accurate and stable ROMs throughout the range of time integration. They are found by means of a two-dimensional grid search across a sufficiently wide range of parameter values. For both the linear and the nonlinear wave equation example, we use the same $\gamma_q$ value for the SMG-QMCL-ROMs and the Galerkin-BQ-ROMs which leads to the same projection error in $\bq$ for both approaches. Furthermore, all the numerical examples in this work use the same regularization parameter values for different basis sizes.

In the following, we introduce the error measures used in the remainder of the paper with a general approximation $\bqgamma$, which can either be $\bgammaLSL$ (<ref>), $\bgammaQMCL$ (<ref>), or $\bgammaBQ$ (<ref>). With $\bqgamma_q$, we denote the restriction of the mapping to the variable $\bq$. The average relative projection error in $\bq(t;\bmu)$ is computed as

\begin{equation}\label{eq:err_rel_proj_train}
{\Rb{\text{err}_{\text{proj,q}} =}}  \frac{1}{M}\sum_{i=1}^M \frac{\lVert \bQ(\bmu_i)- {\color{black}\bqgamma_q} (\bV^\top \bQ(\bmu_i)) \rVert^2_F}{\lVert \bQ(\bmu_i) \rVert^2_F},
\end{equation}

where $\bmu_1, \dots, \bmu_M \in \mathcal{P}$ denotes a set of parameters. This error is valid for both the training and the test parameters.

An equivalent metric ${\Rb{\text{err}_{\text{proj,p}}}}$ is defined for $\bp(t;\bmu)$. The average relative state error is computed as

\begin{equation}\label{eq:err_rel_sim_train}
{\Rb{\text{err}_{\text{sim}}=}}  \frac{1}{M}\sum_{j=1}^M \frac{\lVert \bY(\bmu_j)- {\color{black}\bqgamma} (\widetilde{\bY}(\bmu_j)) \rVert^2_F}{\lVert \bY(\bmu_j)\rVert^2_F},
\end{equation}

where $\widetilde{\bY}(\bmu_{j})$ is obtained from the ROM simulations and $\bqgamma(\widetilde{\bY}(\bmu_j))$ is the reconstruction of a trajectory at parameter value $\bmu_j$ in the original state space. We introduce two separate energy error measures for comparing the preservation of the Hamiltonian in the different ROMs. Since the FOM state trajectories for the linear wave equation preserve the energy exactly, the error in the Hamiltonian for the linear wave equation example in Section <ref> is computed as

\begin{equation}\label{eq:err_hamiltonian}
{\Rb{\Delta H_{\text{lin}}(t)=}}  \left| \dH(\bqgamma(\stateRed(t; \bmu)); \bmu)-\dH(\bqgamma(\stateRed(0; \bmu)); \bmu) \right|,
\end{equation}

where $\dH(\bqgamma(\stateRed(t; \bmu)); \bmu)$ is the FOM energy approximation. For nonlinear Hamiltonian FOMs, the FOM state trajectories do not preserve the energy exactly because the the implicit midpoint scheme only preserves invariances up to the quadratic order. As a result, the error in the Hamiltonian for the nonlinear wave equation example in Section <ref> is computed as

\begin{equation}\label{eq:err_hamiltonian_non}
{\Rb{\Delta H_{\text{nonlin}}(t)=}}  \left| \dH(\bqgamma(\stateRed(t; \bmu)); \bmu)-\dH(\state(t; \bmu); \bmu) \right|.
\end{equation}

§.§ Parametrized linear wave equation

We revisit the linear wave example from Section <ref>. This model problem is similar to the parametric linear wave example in [38]. Let $\Omega=\left(-0.5,0.5\right) \subset \real $ be the spatial domain and consider the parametrized one-dimensional wave equation

\begin{equation}
\frac{\partial^2 }{\partial t^2}\varphi(x,t;\mu)=\mu^2 \frac{\partial^2 }{\partial x^2}\varphi(x,t;\mu),
\label{eq:wave1d_param}
\end{equation}

with the state $\varphi(x,t;\mu)$ at the spatial coordinate $x \in \Omega$, time $t \in  (0,T]$, and the scalar parameter $\mu \in \mathcal {P}=[5/12,5/7]$. Homogenous Dirichlet boundary conditions

\begin{equation}
\varphi(-0.5,t;\mu)=\varphi(0.5,t;\mu)=0,
\end{equation}

are imposed for $t \in (0,T]$ and $\mu \in \mathcal{P}$. We consider a parametric initial condition of the form $\varphi_0(x;\mu):=h(s(x;\mu))$ based on the spline function

\begin{equation}\label{eq:parameteric_ic}
h(s(x;\mu)):= \begin{cases}
1-\frac{3}{2}s(x;\mu)^2 + \frac{3}{4}s(x;\mu)^3& 0\leq s\leq 1, \\
\frac{1}{4}(2-s(x;\mu))^3 & 1\leq s\leq 2, \\
0 & s>2
\end{cases}
\end{equation}

with $s(x;\mu):=\left(4\left|x + \frac{1}{2}-\frac{\mu}{2} \right|\right)/\mu$ for which the exact solution is given by $\varphi(x,t;\mu)=\varphi_0(x-\mu t;\mu)$.

§.§.§ Hamiltonian PDE formulation and FOM implementation details

We rewrite (<ref>) as an infinite-dimensional Hamiltonian system with $q(x,t;\mu):=\varphi(x,t;\mu)$ and $p(x,t;\mu):=\partial \varphi(x,t;\mu)/\partial t$. The associated Hamiltonian functional is

\begin{equation}
\mathcal{H}(q(x,t;\mu), p(x,t;\mu); \mu)
\left[
\frac{1}{2}p(x,t;\mu)^2
+ \frac{1}{2}\mu^2 \left(
\frac{\partial}{\partial x}q(x,t;\mu)
\right)^2
\right] \ \text{d}x,
\label{eq:functional}
\end{equation}

and the original PDE can be recast as a Hamiltonian PDE

\begin{equation}
\frac{\partial }{\partial t}q(x,t;\mu)
=\frac{\delta \cH}{\delta p}{\Rc{(q,p;\mu)}}
=p(x,t;\mu), \qquad
\frac{\partial }{\partial t}p(x,t;\mu)
=-\frac{\delta \cH}{\delta q}{\Rc{(q,p;\mu)}}
=\mu^2\frac{\partial^2 }{\partial x^2}q(x,t;\mu).
\end{equation}

We discretize the spatial domain $\Omega$ with $n=2048$ equally spaced grid points leading to a Hamiltonian FOM of dimension $2n=4096$. Using a finite difference scheme, we obtain the following space-discretized Hamiltonian

\begin{equation}
\dH( \bq(t;\mu), \bp(t;\mu))
=\Delta x\sum^n_{i=1} \left[
\frac{p_i(t;\mu)^2}{2}
+ \frac{\mu^2 (q_{i+1}(t;\mu)-q_i(t;\mu))^2}{4\Delta x^2} + \frac{\mu^2 (q_{i}(t;\mu)-q_{i-1}(t;\mu))^2}{4\Delta x^2}
\right] ,
\end{equation}

\begin{align*}
\bq(t;\mu)&=(q_1(t;\mu), \ldots, q_n(t;\mu) )^\top \in \real^n,\\
p_i(t;\mu)&:=\frac{\partial}{\partial t} \varphi(x_i,t:\mu),&
\bp(t;\mu)&=(p_1(t;\mu), \ldots, p_n(t;\mu) )^\top \in \real^n.
\end{align*}

The corresponding parametrized Hamiltonian FOM equals

\begin{equation}
\dot{\by}(t;\mu) =\begin{pmatrix}
\dot{\bq}(t;\mu)  \\
\dot{\bp}(t;\mu)
\end{pmatrix}=\Jn \nabla_{ \by }\dH(\by(t;\mu))
\bzero & \In \\
\mu^2\mathbf{D}_{\text{fd}} & \bzero
\end{pmatrix}
\begin{pmatrix}
\bq(t;\mu) \\
\bp(t;\mu)
\end{pmatrix}.
\label{eq:wave_FOM}
\end{equation}

§.§.§ Generalization to unseen parameters

Let $\mu_1=0.417$, $\mu_2=0.476$, $\mu_3=0.536$, $\mu_4=0.595$, $\mu_5=0.655$, $\mu_6=0.714 \in \mathcal P$ be $M=6$ parameters equidistantly distributed in $\mathcal{P}$. In this study, a training dataset is built by integrating the Hamiltonian FOM (<ref>) for each training parameter with the implicit midpoint method until final time ${\Rc{T}}=1$. We use a fixed time step of $\Delta t=2.5 \times 10^{-4}$. For this study, we do not shift the trajectory snapshot data, i.e., $\bq_{\text{ref}}=\bp_{\text{ref}}=\bzero$. From these six trajectories, we construct the nonlinear approximation functions $\bgammaQMCL$ (<ref>) and the corresponding SMG-QMCL-ROM, as well as $\bgammaBQ$ (<ref>) and the corresponding Galerkin-BQ-ROM. For this study, we found $\gamma_q=\gamma_p=10^2$ to be a robust choice for SMG-QMCL and Galerkin-BQ ROMs of dimension $4\leq2r\leq 20$. We provide more details on this selection in Section <ref>. We consider $M_{\text{test}}=2$ test parameters $\mu_{\text{test},1}=0.51 \in \mathcal{P}$ and $\mu_{\text{test},2}=0.625\in \mathcal{P}$ to evaluate how these data-driven ROMs generalize for parameter values outside the training dataset.

[subfigure]oneside,margin=1.8cm,0 cm

Relative projection error in $\bq$

Relative projection error in $\bp$

One-dimensional linear wave equation (generalization to unseen parameters). Data-driven approximations based on quadratic manifolds yield lower relative projection error (<ref>) than the linear symplectic subspaces for $\bq$ and $\bp$ variables. The regularization factors are chosen to be $\gamma_q=\gamma_p=10^{2}$.

[subfigure]oneside,margin=1.8cm,0 cm

Training parameters

Test parameters

One-dimensional linear wave equation (generalization to unseen parameters). Galerkin-BQ-ROMs and SMG-QMCL-ROMs achieve lower state error (<ref>) than the linear symplectic ROMs for both training and test parameters.

In Figure <ref>, we compare the relative projection error from (<ref>) for the training data for different values of the reduced dimension. For both $\bq$ and $\bp$ the quadratic manifold approximations yield higher accuracy compared to the linear symplectic subspaces for all reduced dimension values. For this study, we use the same quadratic approximations for $\bq$ in both the SMG-QMCL and the Galerkin-BQ approach, and therefore the relative projection error for $\bq$ in Figure <ref> is the same. The projection error comparison for $\bp$ is shown in Figure <ref> where we observe that the Galerkin-BQ yields marginally lower projection error than the SMG-QMCL.

The comparison of the relative state error (<ref>) between nonlinear SMG-QMCL-ROM and quadratic Galerkin-BQ-ROM is shown in Figure <ref> and Figure <ref> for the training and the test parameters, respectively. Compared to the linear symplectic ROMs, we observe that SMG-QMCL-ROM and Galerkin-BQ-ROM obtain lower relative state error in both training and testing data regimes. The plots in Figure <ref> and Figure <ref> show that the linear symplectic subspace approach produces the least accurate ROMs whereas the Galerkin-BQ-ROM admits the highest accuracy in both training and test data regimes.

The energy error (<ref>) plots in Figures <ref> and <ref> show that the approximately Hamiltonian Galerkin-BQ-ROMs exhibit bounded error in the Hamiltonian for both $\mu_{\text{test},1}=0.51$ and $\mu_{\text{test},2}=0.625$. The Hamiltonian SMG-QMCL-ROMs demonstrate a substantially lower Hamiltonian error due to the usage of the exactly symplectic mapping $\bgammaQMCL$ and the structure-preserving SMG projection. For both approaches, we observe a decrease in the Hamiltonian error when we increase the ROM size from $2r=16$ to $2r=20$.

We compare the approximate numerical solution for $\mu_{\text{test},1}=0.51$ of the linear symplectic ROM, the SMG-QMCL-ROM and the Galerkin-BQ-ROM of size $2r=16$ in Figure <ref>. Even though the FOM solution snapshots for $\mu_{\text{test},1}=0.51$ are not included in the training data, the ROMs based on an quadratic approximation, SMG-QMCL-ROM and Galerkin-BQ-ROM, capture the correct wave shape at $t=1$. The approximate solutions of the linear symplectic ROM, on the other hand, suffer from spurious oscillations, which grow in magnitude as the solution evolves with time.

§.§.§ Sensitivity to regularization prameters

We study the effect of the regularization parameters $\gamma_q$ and $\gamma_p$ on the accuracy of the proposed ROMs in Figure <ref>. There, we fix the reduced dimension to $2r=16$, set $\gamma_p=\gamma_q$, and vary the regularization parameter $\gamma_q$ from $10^{-3}$ to $10^3$. In Figure <ref>, we observe that the relative projection error in $\bq$ for both approaches remain approximately constant up to $\gamma_q=10^{1}$. However, as $\gamma_q$ increases further, the relative projection error increases in magnitude.  For the comparison in Figure <ref>, we observe that the projection error in $\bp$ increases marginally for the SMG-QMCL approach whereas the Galerkin-BQ approach demonstrates nearly constant projection error. In Figure <ref>, we observe that both the SMG-QMCL and the Galerkin-BQ approach yield ROMs with approximately same accuracy up to $\gamma_q=1$. For $\gamma_q>1$, we observe a marginal decrease in the accuracy with increasing $\gamma_q$ for both approaches. Overall, the relative state error values for SMG-QMCL and Galerkin-BQ ROMs of different dimensions remain largely unaffected by the regularization parameter in the investigated range and hence we choose $\gamma_q$ and $\gamma_p$ values that lead to accurate and stable ROMs over a range of reduced dimensions. The results from this sensitivity study illustrate that the proposed approaches provide higher accuracy than the linear symplectic subspace approach over a wide range of regularization parameters. Therefore, the regularization parameters used for the numerical results in Section <ref> are quite robust and can be used to derive accurate and stable ROMs without extensive fine-tuning.

[subfigure]oneside,margin=1.8cm,0 cm

One-dimensional linear wave equation (generalization to unseen parameters). The Hamiltonian error (<ref>) for the Galerkin-BQ-ROMs remains below $10^{-2}$ for $\mu_{\text{test},1}=0.51$ and $\mu_{\text{test},2}=0.625$ whereas the Hamiltonian error for the SMG-QMCL-ROMs remains below $10^{-6}$ for both test parameter values.

One-dimensional linear wave equation (generalization to unseen parameters). Plots show the numerical approximation of the solution of (<ref>) for $\mu_{\text{test}}=0.51$ using low-dimensional ($2r=16$) LSL-ROM, Galerkin-BQ-ROM and SMG-QMCL-ROM at different $t$ values. The Galerkin-BQ-ROMs and SMG-QMCL-ROMs capture the correct wave shape at $t=1$ whereas the LSL-ROM $2r=16$ exhibits spurious oscillations from $t=0.25$ onwards.

Relative projection error in $\bq$

Relative projection error in $\bp$

Relative state error (training)

One-dimensional linear wave equation (sensitivity to regularization parameters). The accuracy of the proposed ROMs with reduced dimension $2r=16$ does not change with the regularization parameter up to $\gamma_q=10^1$. For $\gamma_q>10^{1}$, we observe a marginal decrease in the accuracy for both the SMG-QMCL-ROM and the Galerkin-BQ-ROM.

§.§.§ Time extrapolation study

A key motivation for symplectic model reduction using data-driven quadratic manifolds is to obtain low-dimensional ROMs that can capture the periodic behavior in an accurate and stable manner. In this study, we fix the scalar parameter in (<ref>) to $\mu=0.5$ and build a training dataset of snapshots over one cycle by integrating the Hamiltonian FOM with the implicit midpoint method until time $t=4$ with a fixed time step $\Delta t=10^{-3}$. Note that we use the same training dataset with the same fixed time step $\Delta t=10^{-3}$ as in the motivational example in Section <ref>. For this time extrapolation study, we found $\gamma_q=\gamma_p=10^6$ to be a robust choice for the proposed approaches. For test data, we consider the FOM solution snapshots over the next nine cycles, i.e., from $t=4$ to $t=40$. For this time extrapolation study, we center the trajectory snapshot data about the initial condition with $\bq_{\text{ref}}=\bq(0)$ and $\bp_{\text{ref}}=\bp(0)$.

[subfigure]oneside,margin=1.8cm,0 cm

State error (training data)

FOM energy error

One-dimensional linear wave equation (time extrapolation study).  Plot (a) shows that LSL-ROMs, Galerkin-BQ-ROMs, and SMG-QMCL-ROMs achieve similar accuracy in the training data regime. Plot(b) compares the energy error behavior for ROMs of different sizes. The nonlinear SMG-QMCL-ROMs exhibit bounded error in the Hamiltonian (<ref>) outside the training data regime. In contrast, the Hamiltonian errors for the quadratic Galerkin-BQ-ROMs slowly grow with time.  The dashed black line indicates end of training time interval at $t=4$.

We compare the accuracy and the energy error performance of the proposed ROMs in Figure <ref>. In Figure <ref>, we observe that both the SMG-QMCL and the Galerkin-BQ approach yield ROMs with similar accuracy in the training data regime. Compared with Figure <ref>, we observe that the proposed approaches yield stable ROMs that demonstrate a decrease in state error with increasing reduced dimension. The energy error plots in Figure <ref> demonstrate stability outside the training data regime for both the SMG-QMCL-ROMs and the Galerkin-BQ-ROMs of different reduced dimensions. For both approaches, we observe that the error in the Hamiltonian decreases marginally with an increase in the reduced dimension $2r$. All three SMG-QMCL-ROMs exhibit bounded energy error and the error in the Hamiltonian for all three ROMs remain below an absolute value of $10^{-7}$ outside the training time interval. Despite the Hamiltonian error for the Galerkin-BQ-ROM remaining below an absolute value of $10^{-1}$ over the entire observation interval $t \in [0,40]$, we do observe a slow growth because of its approximately symplectic structure.

The space-time evolution of the FOM solution field is compared with the space-time evolution of the approximate solutions obtained using ROMs of dimension $2r=24$ in Figure <ref>. While all three ROMs approximate the FOM solution accurately in the training data regime, we observe that the linear symplectic ROM exhibits spurious oscillations that gradually become more pronounced as we march forward in time in Figure <ref>. In Figure <ref> and Figure <ref>, we observe that both proposed approaches provide stable and accurate approximate solutions at time $t=40$, which is $900\%$ outside the training time interval. Unlike the conventional quadratic manifolds approach in Section <ref> which can lead to unstable ROMs, the long-time stable and accurate predictions of future cycles in this time extrapolation study highlights the utility of the proposed methods for high-dimensional systems with periodic behavior.

FOM solution

One-dimensional linear wave equation (time extrapolation study). Plots show the numerical approximation of the
FOM solution using low-dimensional ($2r=24$) LSL-ROM, SMG-QMCL-ROM, and Galerkin-BQ-ROM for $t \in [0,40]$. The dashed black line indicates the end of the training time interval at $t=4$. Plot (b) shows that the LSL-ROM yields inaccurate solutions after $t=8$ whereas plot (c) and plot (d) show that SMG-QMCL-ROM and Galerkin-BQ-ROM provide accurate approximate solutions even at $t=40$.

§.§ Parametrized two-dimensional nonlinear wave equation

For a second numerical example, we consider the two-dimensional nonlinear wave equation with parametric nonlinearity. This parametric model problem is similar to the two-dimensional nonlinear wave equation example in [31]. Let $\Omega=(-10,10) \times (-10,10) \subset \real^2$ be the spatial domain and consider the parametrized two-dimensional nonlinear wave equation

\begin{equation}\label{eq:non_pde}
\frac{\partial^2 \varphi}{\partial t^2}(x,y,t;\mu)=\frac{\partial^2 \varphi}{\partial x^2} (x,y,t;\mu)+ \frac{\partial^2 \varphi}{\partial y^2}(x,y,t;\mu) - \mu\varphi(x,y,t;\mu)^3,
\end{equation}

with the state $\varphi(x,y,t;\mu)$ at $(x,y)\in \Omega$, time $t\in(0,T]$, and the scalar parameter $\mu \in \mathcal P=[0.1,3]$. Periodic boundary conditions

\begin{equation}
\varphi(-10,y,t;\mu)=\varphi(10,y,t;\mu), \qquad \varphi(x,-10,t;\mu)=\varphi(x,10,t;\mu),
\end{equation}

are imposed for $t \in (0,T]$ and $\mu \in \mathcal P$.

§.§.§ Hamiltonian PDE formulation and FOM implementation details

We rewrite (<ref>) as an infinite-dimensional canonical Hamiltonian system with $q(x,y,t;\mu):=\varphi(x,y,t;\mu)$ and $p(x,y,t;\mu):=\partial \varphi(x,y,t;\mu)/\partial t$. The associated space-time continuous Hamiltonian functional is

\begin{multline}\label{eq:fham_non}
\mathcal H(q(x,y,t;\mu),p(x,y,t;\mu);\mu)=\frac{1}{2}\int_{\Omega} \left[ p(x,y,t;\mu)^2 + \left(\frac{\partial}{\partial x}q(x,y,t;\mu)\right)^2 +  \left(\frac{\partial}{\partial y}q(x,y,t;\mu)\right)^2 \right] \mathrm{d}x\mathrm{d}y \\
+ \int_{\Omega} \frac{\mu}{4}q(x,y,t;\mu)^4  \mathrm{d}x\mathrm{d}y
\end{multline}

and the original nonlinear PDE can be recast as a canonical Hamiltonian PDE

\begin{align*}
\frac{\partial }{\partial t}q(x,y,t;\mu)&=\frac{\textcolor{black}{\delta} \mathcal H}{\textcolor{black}{\delta} p}(q,p;\mu)=p(x,y,t;\mu),\\
\frac{\partial }{\partial t}p(x,y,t;\mu)&=-\frac{\textcolor{black}{\delta} \mathcal H}{\textcolor{black}{\delta} q}(q,p;\mu)=\frac{\partial ^2}{\partial x^2}q(x,y,t;\mu)+\frac{\partial ^2}{\partial y^2}q(x,y,t;\mu)-\mu q(x,y,t;\mu)^3.
\end{align*}

We discretize the two-dimensional spatial domain $\Omega$ with $n_x=n_y=100$ equally spaced grid points in both spatial directions leading to system states of dimension $2n$ with $n=n_xn_y=10,000$. Using a finite difference scheme in both spatial directions, we obtain the following space-discretized Hamiltonian

\begin{multline}\label{eq:dham_non}
H(\bq(t;\mu),\bp(t;\mu);\mu)=\Delta x \Delta y \sum_{i,j=1}^{n_x,n_y} \left[ p_{i,j}(t;\mu)^2 + \frac{ (q_{i+1,j}(t;\mu)-q_{i,j}(t;\mu))^2}{4\Delta x^2} + \frac{ (q_{i,j}(t;\mu)-q_{i-1,j}(t;\mu))^2}{4\Delta x^2}\right] \\
+\Delta x \Delta y  \sum_{i,j=1}^{n_x,n_y} \left[ \frac{ (q_{i,j+1}(t;\mu)-q_{i,j}(t;\mu))^2}{4\Delta y^2} + \frac{ (q_{i,j}(t;\mu)-q_{i,j-1}(t;\mu))^2}{4\Delta y^2} +  \frac{\mu}{4}q_{i,j}(t;\mu)^4 \right],
\end{multline}

\begin{align*}
\bq(t;\mu)&=(q_{1,1}(t;\mu), \ldots, q_{n_x,n_y}(t;\mu) )^\top \in \real^n,\\
p_{i,j}(t;\mu)&:=\frac{\partial}{\partial t} \varphi(x_i,y_j,t;\mu),&
\bp(t;\mu)&=(p_{1,1}(t;\mu), \ldots, p_{n_x,n_y}(t;\mu) )^\top \in \real^n.
\end{align*}

The corresponding parametrized nonlinear Hamiltonian FOM is

\begin{equation}
\dot{\by}(t;\mu) =\begin{pmatrix}
\dot{\bq}(t;\mu)  \\
\dot{\bp}(t;\mu)
\end{pmatrix}=\Jn \nabla_{ \by }\dH(\by(t;\mu))
\bzero & \In \\
\mathbf{D}_{\text{fd},2d} & \bzero
\end{pmatrix}
\begin{pmatrix}
\bq(t;\mu) \\
\bp(t;\mu)
\end{pmatrix} - \mu\begin{pmatrix}
\bzero \\
\bq(t;\mu)^3
\end{pmatrix} ,
\label{eq:wave_2d_FOM}
\end{equation}

where $\mathbf{D}_{\text{fd},2d}$ denotes the finite difference approximation in the two-dimensional setting and the vector $ \bq(t;\mu)^3\in \real^{n}$ contains as components the entry-wise cubic exponential of the generalized state vector $\bq(t;\mu)$ .

[subfigure]oneside,margin=1.8cm,0 cm

Relative projection error in $\bq$

Relative projection error in $\bp$

Two-dimensional nonlinear wave equation. Both the QMCL and the BQ state approximation yield a substantially lower relative projection error (<ref>) than the linear symplectic subspaces for $\bq$ and $\bp$ variables. The regularization factors are chosen to be $\gamma_q=\gamma_p=10^{-1}$.

[subfigure]oneside,margin=1.8cm,0 cm

Training parameters

Test parameters

Two-dimensional nonlinear wave equation. The proposed SMG-QMCL and Galerkin-BQ ROMs  achieve lower state error (<ref>) than the LSL-ROMs for both training and test parameters.

§.§.§ Parameter extrapolation study

Let $\mu_1,\dots,\mu_{10} \in \mathcal P_{\text{train}}$ be $M=10$ parameters equidistantly distributed (including the endpoints) in $\mathcal P_{\text{train}}=[0.1,1]$. For this parameter extrapolation study, we build a training dataset by integrating the nonlinear Hamiltonian FOM (<ref>) for each training parameter with the implicit midpoint method until final time $T=8$. We use a fixed time step of $\Delta t=0.1$. In this study, we do not shift the trajectory snapshot data, i.e., $\bq_{\text{ref}}=\bp_{\text{ref}}=\bzero$. From this training dataset, we first construct data-driven nonlinear approximation functions, i.e., $\bgammaQMCL$ (<ref>) and $\bgammaBQ$ (<ref>). We then derive the corresponding nonlinear SMG-QMCL-ROM (<ref>) and nonlinear Galerkin-BQ-ROM (<ref>) for different reduced dimensions. For this parameter extrapolation study, we found $\gamma_q=\gamma_p=10^{-1}$ to be a robust choice for both approaches. We also consider $M_{\text{test}}=4$ test parameters $\mu_{\text{test},1}=1.25$, $\mu_{\text{test},2}=1.5$, $\mu_{\text{test},3}=2$, and $\mu_{\text{test},4}=3$ to evaluate how the proposed nonlinear data-driven ROMs generalize for parameter values outside the training dataset.

In Figure <ref>, we compare the relative projection error in $\bq$ and $\bp$ for the training data consisting of ten trajectories. For both $\bq$ and $\bp$, the nonlinear approximations based on data-driven quadratic manifolds yield a substantially higher accuracy compared to the linear symplectic subspaces for all reduced dimensions. Similar to the linear wave equation example, we use the same quadratic approximations for $\bq$ in both the SMG-QMCL and the Galerkin-BQ approach, and therefore the relative projection error for $\bq$ in Figure <ref> is the same. In Figure <ref>, we observe that the BQ approximation mapping yields a lower projection error than the QMCL approximation for $\bp$.

In Figure <ref>, we observe that the proposed SMG-QMCL and Galerkin-BQ approaches provide ROMs with lower relative state error than the linear symplectic ROMs for both training and test parameters. The comparison plots for the training parameters are shown in Figure <ref> where we observe that the proposed approaches achieve substantially lower state error than the symplectic ROMs based on linear symplectic subspaces. The comparison for the test parameters in Figure <ref> shows that the ROMs based on SMG-QMCL and Galerkin-BQ approaches achieve two to four times lower relative state error than the ROMs based on the linear symplectic subspace approach.

The Hamiltonian error (<ref>) plots in Figure <ref> show that the approximate solutions obtained using SMG-QMCL and Galerkin-BQ ROMs of different reduced dimensions  accurately approximate the FOM Hamiltonian for both $\mu_{\text{test},1}=1.25$ and $\mu_{\text{test},4}=3$. For $\mu_{\text{test},1}=1.25$, the error in the Hamiltonian in Figure <ref> for the proposed approaches remain below $10^{-2}$. In Figure <ref>, we observe that the energy error for the Hamiltonian ROMs based on the SMG-QMCL approach level off at approximately $5 \times 10^{-3}$. The approximately Hamiltonian ROMs based on the Galerkin-BQ approach, on the other hand, exhibit oscillatory error behavior with the error in the Hamiltonian slowly growing with time.

Finally, we study the accuracy of the proposed ROMs over the two-dimensional computational domain by comparing the pointwise error in $\bq$ between the FOM solution and the reconstructed solution $\bGamma_{\bq}(\widetilde{\bQ}(\mu))$ of the ROM solution $\widetilde{\bQ}(\mu)$ for two test parameter values. The time-evolution of the pointwise error in $\bq$ for the nonlinear Hamiltonian ROM of size $2r=48$ based on the linear symplectic subspace is compared with the corresponding errors for the SMG-QMCL-ROM and the Galerkin-BQ-ROM in Figure <ref> and Figure <ref> for $\mu_{\text{test},1}=1.25$ and $\mu_{\text{test},4}=3$, respectively. For $\mu_{\text{test},1}=1.25$, we observe that the nonlinear ROMs based on SMG-QMCL and Galerkin-BQ approaches yield reconstructed solution fields for $\bq$ with a maximal pointwise error of approximately $10^{-3}$ at all three time instances. In contrast, we observe that the maximal pointwise error in $\bq$ for the linear symplectic subspace approach increases from approximately $5 \times 10^{-3}$ to $2 \times 10^{-2}$ as we march forward in time. For $\mu_{\text{test},4}=3$, we observe that both the SMG-QMCL and the Galerkin-BQ ROM exhibit a maximal pointwise error of approximately $10^{-2}$ whereas the  linear symplectic subspace approach exhibits a maximal pointwise error of approximately $10^{-1}$ which is an order of magnitude higher than the error for proposed approaches. These results demonstrate the ability of both the SMG-QMCL and the Galerkin-BQ approach to provide accurate numerical solutions for multi-dimensional nonlinear problems, even for parameter values outside the range of training parameters.

[subfigure]oneside,margin=1.8cm,0 cm

Two-dimensional nonlinear wave equation. The error in the Hamiltonian (<ref>) for SMG-QMCL-ROMs and Galerkin-BQ-ROMs remains below $10^{-2}$ for $\mu_{\text{test},1}=0.125$  whereas both approaches exhibit bounded error in the Hamiltonian below $10^{-1}$ for $\mu_{\text{test},4}=3$.

Two-dimensional nonlinear wave equation. Plots compare the pointwise error in $\bq$ between the FOM solution and the reconstructed solution $\bGamma_{\bq}(\widetilde{\bQ}(\mu))$ of the ROM solution $\widetilde{\bQ}(\mu)$ using low-dimensional ($2r=48$) LSL-ROM, SMG-QMCL-ROM, and Galerkin-BQ-ROM at selected time instances $t\in\{ 6,7,8\}$ for $\mu_{\text{test,1}}=1.25$. The pointwise error for the linear symplectic ROM increases as we march forward in time whereas the SMG-QMCL ROM and Galerkin-BQ-ROM yield accurate approximate solutions even at $t=8$.

Two-dimensional nonlinear wave equation. Plots compare the pointwise error in $\bq$ between the FOM solution and the reconstructed solution $\bGamma_{\bq}(\widetilde{\bQ}(\mu))$ of the ROM solution $\widetilde{\bQ}(\mu)$ using low-dimensional ($2r=48$) LSL-ROM, SMG-QMCL-ROM, and Galerkin-BQ-ROM at selected time instances $t\in\{ 6,7,8\}$ for $\mu_{\text{test,4}}=3$. The proposed approaches yield a lower pointwise error than the LSL-ROM at all three time instances.

§ CONCLUSIONS & FUTURE WORK

We have presented two projection-based model reduction approaches that use data-driven quadratic manifolds to derive accurate structure-preserving reduced-order models of canonical Hamiltonian systems. The SMG-QMCL approach derives Hamiltonian ROMs using the quadratic manifold cotangent lift mapping which is based on a quadratic manifold approximation for the generalized position vector. More generally, we proposed the manifold cotangent lift approximation mapping which is a symplectic map and can thus be used with the SMG projection to build a symplectic SMG-MCL-ROM. Note, that this construction does not require the approximation to be a quadratic mapping, as in the QMCL, so the SMG-MCL-ROM could be used to formulate symplectic MOR with other nonlinear approximation mappings like e.g., higher-order polynomials or autoencoders in future work. The Galerkin-BQ approach, on the other hand, augments the linear symplectic subspace with quadratic Kronecker product terms to derive approximately Hamiltonian ROMs. Both approaches are offline–online separable for linear Hamiltonian systems which means that the ROMs of both approaches are independent of the FOM dimension. The Galerkin-BQ-ROM approach is particularly attractive from a computational efficiency standpoint due to the ROMs admitting a convenient linear-quadratic model structure with computational complexity on the order of $\bigO(r^4)$ compared to $\bigO(r^5)$ for the SMG-QMCL-ROM. On the other hand, the SMG-QMCL-ROM is favorable in terms of long-term energy conservation. Both of these novel approaches together constitute a first step towards the model reduction of dynamical systems on nonlinear manifolds using interpretable (e.g., polynomial) manifold constructions that ensure that the approximate solution satisfies key physical properties as dictated by the original high-dimensional problem.

The numerical experiments with the parametrized linear and nonlinear wave equations demonstrate the advantages of using data-driven quadratic manifold approximations for structure-preserving model reduction of transport-dominated problems. The numerical results also show that the proposed approaches produce stable ROMs with higher accuracy and better predictive capabilities than their linear counterparts. Notably, the proposed methods yield accurate predictions of the high-dimensional state even for parameter values that were not included in the training data.

This work has opened a number of avenues for future work. Improving the computational efficiency of the proposed nonlinear ROMs based on data-driven quadratic manifolds is necessary. To reduce the online computational cost, the proposed approaches could be combined with the very recently developed structure-preserving DEIM [45] to improve the computational efficiency of the nonlinear SMG-QMCL and Galerkin-BQ ROMs. To further reduce the computational cost, alternative formulations for learning manifold approximations in high-dimensional state spaces could be exploited to formulate more suitable symplectic mappings. In another direction, the proposed methodologies could be extended to noncanonical Hamiltonian systems with state-dependent and degenerate forms by building on the recent linear subspace approaches [28]. Finally, we seek to broaden the application of the presented framework by learning Hamiltonian dynamics directly from time-domain simulation data in a non-intrusive fashion using operator inference methods [49, 58, 25].

§ ACKNOWLEDGMENTS
H.S. and B.K. were in part financially supported by the U.S. Office of Naval Research under award number N00014-22-1-2624 and the Ministry of Trade, Industry and Energy (MOTIE) and the Korea Institute for Advancement of Technology (KIAT) through the International Cooperative R&D program (No. P0019804, Digital twin based intelligent unmanned facility inspection solutions). R.G. has been supported in part by the U.S. Department of Energy AEOLUS MMICC center under award DE-SC0019303, program manager W. Spotz, and by the AFOSR MURI on physics-based machine learning, award FA9550-21-1-0084, program manager F. Fahroo.
P.B. is funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC 2075 - 390740016.
P.B. acknowledges the support by the Stuttgart Center for Simulation Science (SimTech).

[1]
R. Abraham and J. E. Marsden.
Foundations of mechanics.
Number 364. American Mathematical Soc., 2008.

[2]
S. E. Ahmed, O. San, A. Rasheed, and T. Iliescu.
Nonlinear proper orthogonal decomposition for convection-dominated
Physics of Fluids, 33(12):121702, 2021.

[3]
D. Amsallem, M. J. Zahr, and C. Farhat.
Nonlinear model order reduction based on local reduced-order bases.
International Journal for Numerical Methods in Engineering,
92(10):891–916, 2012.

[4]
J. Barnett and C. Farhat.
Quadratic approximation manifold for mitigating the Kolmogorov
barrier in nonlinear projection-based model order reduction.
Journal of Computational Physics, 464:111348, 2022.

[5]
J. L. Barnett, C. Farhat, and Y. Maday.
Neural-network-augmented projection-based model order reduction for
mitigating the Kolmogorov barrier to reducibility of CFD models.
arXiv preprint arXiv:2212.08939, 2022.

[6]
T. Bendokat and R. Zimmermann.
Geometric optimization for structure-preserving model reduction of
Hamiltonian systems.
In Proceedings of the 10th Vienna International Conference on
Mathematical Modelling, 2022.

[7]
T. J. Bridges and S. Reich.
Numerical methods for Hamiltonian PDEs.
Journal of Physics A: Mathematical and General, 39(19):5287,

[8]
P. Buchfink, A. Bhatt, and B. Haasdonk.
Symplectic model order reduction with non-orthonormal bases.
Mathematical and Computational Applications, 24(2):43, 2019.

[9]
P. Buchfink, S. Glas, and B. Haasdonk.
Optimal bases for symplectic model order reduction of canonizable
linear Hamiltonian systems.
IFAC-PapersOnLine, 55(20):463–468, 2022.
10th Vienna International Conference on Mathematical Modelling
MATHMOD 2022.

[10]
P. Buchfink, S. Glas, and B. 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.

[11]
P. Buchfink, B. Haasdonk, and S. Rave.
Proceedings of the Conference Algoritmy 2020, chapter
PSD-Greedy Basis Generation for Structure-Preserving Model Order Reduction of
Hamiltonian Systems, pages 151–160.
Vydavateľstvo SPEKTRUM, 2020.

[12]
E. Celledoni, V. Grimm, R. I. McLachlan, D. McLaren, D. O'Neale, B. Owren, and
G. Quispel.
Preserving energy resp. dissipation in numerical PDEs using the
“Average Vector Field” method.
Journal of Computational Physics, 231(20):6770–6789, 2012.

[13]
T. Daniel, F. Casenave, N. Akkari, A. Ketata, and D. Ryckelynck.
Physics-informed cluster analysis and a priori efficiency criterion
for the construction of local reduced-order bases.
Journal of Computational Physics, 458:111120, 2022.

[14]
T. Daniel, F. Casenave, N. Akkari, and D. Ryckelynck.
Model order reduction assisted by deep neural networks (ROM-net).
Advanced Modeling and Simulation in Engineering Sciences,
7(1):1–27, 2020.

[15]
S. Fresca, L. Dede, and A. 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.

[16]
S. Fresca and A. Manzoni.
POD-DL-ROM: Enhancing deep learning-based reduced order models for
nonlinear parametrized PDEs by proper orthogonal decomposition.
Computer Methods in Applied Mechanics and Engineering,
388:114181, 2022.

[17]
H. Gao, J.-X. Wang, and M. J. Zahr.
Non-intrusive model reduction of large-scale, nonlinear dynamical
systems using deep learning.
Physica D: Nonlinear Phenomena, 412:132614, 2020.

[18]
R. Geelen and K. Willcox.
Localized non-intrusive reduced-order modeling in the operator
inference framework.
Philosophical Transactions of the Royal Society A, 380(2229),

[19]
R. Geelen, S. Wright, and K. Willcox.
Operator inference for non-intrusive model reduction with quadratic
Computer Methods in Applied Mechanics and Engineering,
403:115717, 2023.

[20]
G. H. Golub and V. Pereyra.
The differentiation of pseudo-inverses and nonlinear least squares
problems whose variables separate.
SIAM Journal on Numerical Analysis, 10(2):413–432, 1973.

[21]
Y. Gong, Q. Wang, and Z. Wang.
Structure-preserving Galerkin POD reduced-order modeling of
Hamiltonian systems.
Computer Methods in Applied Mechanics and Engineering,
315:780–798, 2017.

[22]
F. J. Gonzalez and M. Balajewicz.
Deep convolutional recurrent autoencoders for learning
low-dimensional feature dynamics of fluid systems.
arXiv preprint arXiv:1808.01346, 2018.

[23]
O. Gonzalez.
On the stability of symplectic and energy-momentum algorithms for
non-linear Hamiltonian systems with symmetry.
Computer Methods in Applied Mechanics and Engineering,
134(3-4):197–222, 1996.

[24]
C. Greif and K. Urban.
Decay of the Kolmogorov ${N}$-width for wave problems.
Applied Mathematics Letters, 96:216–222, 2019.

[25]
A. Gruber and I. Tezaur.
Canonical and noncanonical Hamiltonian operator inference.
arXiv preprint arXiv:2304.06262, 2023.

[26]
B. Haasdonk and M. Ohlberger.
Efficient reduced models and a posteriori error estimation for
parametrized dynamical systems by offline/online decomposition.
Mathematical and Computer Modelling of Dynamical Systems,
17(2):145–161, 2011.

[27]
E. Hairer, C. Lubich, and G. Wanner.
Geometric Numerical Integration: Structure-Preserving Algorithms
for Ordinary Differential Equations, volume 31.
Springer Science & Business Media, 2006.

[28]
J. Hesthaven and C. Pagliantini.
Structure-preserving reduced basis methods for Poisson systems.
Mathematics of Computation, 90(330):1701–1740, 2021.

[29]
J. S. Hesthaven, C. Pagliantini, and N. Ripamonti.
Rank-adaptive structure-preserving model order reduction of
Hamiltonian systems.
ESAIM: Mathematical Modelling and Numerical Analysis,
56(2):617–650, 2022.

[30]
O. Issan and B. Kramer.
Predicting solar wind streams from the inner-heliosphere to earth via
shifted operator inference.
Journal of Computational Physics, 473:111689, 2023.

[31]
C. Jiang, Y. Wang, and W. Cai.
A linearly implicit energy-preserving exponential integrator for the
nonlinear Klein-Gordon equation.
Journal of Computational Physics, 419:109690, 2020.

[32]
T. Kadeethum, F. Ballarin, Y. Choi, D. O’Malley, H. Yoon, and N. Bouklas.
Non-intrusive reduced order modeling of natural convection in porous
media using convolutional autoencoders: Comparison with linear subspace
Advances in Water Resources, 160:104098, 2022.

[33]
Y. Kim, Y. Choi, D. Widemann, and T. Zohdi.
A fast and accurate physics-informed neural network reduced order
model with shallow masked autoencoder.
Journal of Computational Physics, 451:110841, 2022.

[34]
O. Koch and C. Lubich.
Dynamical low-rank approximation.
SIAM Journal on Matrix Analysis and Applications,
29(2):434–454, 2007.

[35]
A. Kolmogoroff.
Uber die beste Annaherung von Funktionen einer gegebenen
Annals of Mathematics, pages 107–110, 1936.

[36]
K. Lee and K. T. Carlberg.
Model reduction of dynamical systems on nonlinear manifolds using
deep convolutional autoencoders.
Journal of Computational Physics, 404:108973, 2020.

[37]
B. Leimkuhler and S. Reich.
Simulating Hamiltonian Dynamics.
Cambridge Monographs on Applied and Computational Mathematics.
Cambridge University Press, 2005.

[38]
B. Maboudi Afkham and J. S. Hesthaven.
Structure preserving model reduction of parametric Hamiltonian
SIAM Journal on Scientific Computing, 39(6):A2616–A2644, 2017.

[39]
J. E. Marsden and T. S. Ratiu.
Introduction to Mechanics and Symmetry: A Basic Exposition of
Classical Mechanical Systems, volume 17.
Springer Science & Business Media, 2013.

[40]
R. McLachlan.
Symplectic integration of Hamiltonian wave equations.
Numerische Mathematik, 66:465–492, 1993.

[41]
M. Nonino, F. Ballarin, G. Rozza, and Y. Maday.
Overcoming slowly decaying Kolmogorov n-width by transport maps:
Application to model order reduction of fluid dynamics and fluid–structure
interaction problems.
arXiv preprint arXiv:1911.06598, 2019.

[42]
M. Ohlberger and S. Rave.
Nonlinear reduced basis approximation of parameterized evolution
equations via the method of freezing.
Comptes Rendus Mathematique, 351(23-24):901–906, 2013.

[43]
M. Ohlberger and S. Rave.
Reduced basis methods: Success, limitations and future challenges.
In Proceedings of Algoritmy, pages 1–12, 2016.

[44]
C. Pagliantini.
Dynamical reduced basis methods for Hamiltonian systems.
Numerische Mathematik, 148:409–448, 2021.

[45]
C. Pagliantini and F. Vismara.
Gradient-preserving hyper-reduction of nonlinear dynamical systems
via discrete empirical interpolation.
arXiv preprint arXiv:2206.01792, 2022.

[46]
D. Papapicco, N. Demo, M. Girfoglio, G. Stabile, and G. Rozza.
The neural network shifted-proper orthogonal decomposition: A machine
learning approach for non-linear reduction of hyperbolic equations.
Computer Methods in Applied Mechanics and Engineering,
392:114687, 2022.

[47]
S. Pawar, S. Rahman, H. Vaddireddy, O. San, A. Rasheed, and P. Vedula.
A deep learning enabler for nonintrusive reduced order modeling of
fluid flows.
Physics of Fluids, 31(8):085101, 2019.

[48]
B. Peherstorfer, D. Butnaru, K. Willcox, and H.-J. Bungartz.
Localized discrete empirical interpolation method.
SIAM Journal on Scientific Computing, 36(1):A168–A192, 2014.

[49]
B. Peherstorfer and K. Willcox.
Data-driven operator inference for nonintrusive projection-based
model reduction.
Computer Methods in Applied Mechanics and Engineering,
306:196–215, 2016.

[50]
L. Peng and K. Mohseni.
Symplectic model reduction of Hamiltonian systems.
SIAM Journal on Scientific Computing, 38(1):A1–A27, 2016.

[51]
J. Reiss, P. Schulze, J. Sesterhenn, and V. Mehrmann.
The shifted proper orthogonal decomposition: A mode decomposition for
multiple transport phenomena.
SIAM Journal on Scientific Computing, 40(3):A1322–A1344, 2018.

[52]
O. San, R. Maulik, and M. Ahmed.
An artificial neural network framework for reduced order modeling of
transient flows.
Communications in Nonlinear Science and Numerical Simulation,
77:271–287, 2019.

[53]
J. M. Sanz-Serna.
Symplectic integrators for Hamiltonian problems: An overview.
Acta numerica, 1:243–286, 1992.

[54]
T. P. Sapsis and P. F. Lermusiaux.
Dynamically orthogonal field equations for continuous stochastic
dynamical systems.
Physica D: Nonlinear Phenomena, 238(23-24):2347–2360, 2009.

[55]
H. Sharma, J. Borggaard, M. Patil, and C. Woolsey.
Performance assessment of energy-preserving, adaptive time-step
variational integrators.
Communications in Nonlinear Science and Numerical Simulation,
114:106646, 2022.

[56]
H. Sharma and B. Kramer.
Preserving Lagrangian structure in data-driven reduced-order
modeling of large-scale dynamical systems.
arXiv preprint arXiv:2203.06361, 2022.

[57]
H. Sharma, M. Patil, and C. Woolsey.
A review of structure-preserving numerical methods for engineering
Computer Methods in Applied Mechanics and Engineering,
366:113067, 2020.

[58]
H. Sharma, Z. Wang, and B. Kramer.
Hamiltonian operator inference: Physics-preserving learning of
reduced-order models for canonical Hamiltonian systems.
Physica D: Nonlinear Phenomena, 431:133122, 2022.

[59]
T. Taddei, S. Perotto, and A. Quarteroni.
Reduced basis techniques for nonlinear conservation laws.
ESAIM: Mathematical Modelling and Numerical Analysis,
49(3):787–814, 2015.

[60]
H. Yoshida.
Construction of higher order symplectic integrators.
Physics Letters A, 150(5-7):262–268, 1990.