Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders

TL;DR

This paper introduces a novel nonlinear manifold projection framework for dynamical systems using deep convolutional autoencoders, overcoming limitations of traditional linear subspace methods especially for advection-dominated problems.

Contribution

It proposes a new approach combining manifold projection with deep autoencoders, providing theoretical analysis and error bounds, and demonstrating superior performance on benchmark problems.

Findings

Outperforms linear subspace ROMs on advection-dominated problems

Provides a posteriori error bounds for the proposed methods

Uses deep convolutional autoencoders to efficiently compute nonlinear manifolds

Abstract

Nearly all model-reduction techniques project the governing equations onto a linear subspace of the original state space. Such subspaces are typically computed using methods such as balanced truncation, rational interpolation, the reduced-basis method, and (balanced) POD. Unfortunately, restricting the state to evolve in a linear subspace imposes a fundamental limitation to the accuracy of the resulting reduced-order model (ROM). In particular, linear-subspace ROMs can be expected to produce low-dimensional models with high accuracy only if the problem admits a fast decaying Kolmogorov $n$-width (e.g., diffusion-dominated problems). Unfortunately, many problems of interest exhibit a slowly decaying Kolmogorov $n$-width (e.g., advection-dominated problems). To address this, we propose a novel framework for projecting dynamical systems onto nonlinear manifolds using minimum-residual formulations at the time-continuous and time-discrete levels; the former leads to manifold Galerkin projection, while the latter leads to manifold least-squares Petrov--Galerkin (LSPG) projection. We perform analyses that provide insight into the relationship between these proposed approaches and classical linear-subspace reduced-order models; we also derive a posteriori discrete-time error bounds for the proposed approaches. In addition, we propose a computationally practical approach for computing the nonlinear manifold, which is based on convolutional autoencoders from deep learning. Finally, we demonstrate the ability of the method to significantly outperform even the optimal linear-subspace ROM on benchmark advection-dominated problems, thereby demonstrating the method's ability to overcome the intrinsic $n$-width limitations of linear subspaces.

Figures

Tables

Table 1. Table 1 : Parameters of the autoencoder architecture described in Figure 1 applied to both benchmark problems. The encoder consists of n L = 5 subscript 𝑛 L 5 n_{\text{L}}={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}5}} layers with n conv = 4 subscript 𝑛 conv 4 n_{\text{conv}}=4 convolutional layers, followed by n full = 1 subscript 𝑛 full 1 n_{\text{full}}={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1}} fully-connected layers. The decoder consists of n full = 1 subscript 𝑛 full 1 n_{\text{full}}={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1}} fully-connected layers, followed by n conv = 4 subscript 𝑛 conv 4 n_{\text{conv}}=4 transposed-convolution layers. For the parameterized 1D Burgers’ equation, we employ 1D convolutional kernel filters with the kernel length 25 at every convolutional layer and transposed- convolutional layer, and apply length- 2 stride ( s = 2 𝑠 2 s={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2}} ) encoder layer 1 and the decoder layer 5 , and length- 4 stride ( s = 4 𝑠 4 s={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}4}} ) for other convolutional layers and transposed-convolutional layers. For the chemically reacting flow, we employ 5 × 5 5 5 5\times 5 convolutional kernel filters at every convolutional layer and transposed-convolutional layer, and apply length-2 stride ( s = 2 𝑠 2 s=2 ). For the definition of stride, we refer to B . For zero-padding, we use half padding [ 26 , 30 ] . With these settings, p n conv = p ¯ n full = 128 subscript 𝑝 subscript 𝑛 conv subscript ¯ 𝑝 subscript 𝑛 full 128 p_{n_{\text{conv}}}=\bar{p}_{n_{\text{full}}}=128 for the 1D Burgers’ equation and p n conv = p ¯ n full = 512 subscript 𝑝 subscript 𝑛 conv subscript ¯ 𝑝 subscript 𝑛 full 512 p_{n_{\text{conv}}}=\bar{p}_{n_{\text{full}}}=512 for the chemically reacting flow.

Encoder network Convolution layers Layer Number of filters 1 8 2 16 3 32 4 64 Fully-connected layers Layer Input dimension Output dimension 5 $p_{n_{\text{conv}}}$ $p$

Decoder network Fully-connected layers Layer Input dimension Output dimension 1 $p$ ${\overline{p}}_{n_{\text{full}}}$ Transposed-convolutional layers Layer Number of filters 2 64 3 32 4 16 5 1

Table 2. Table 2 : 1D Burgers’ equation. Relative errors of the ROM solutions (Eq. ( 7.1 )) at online-parameter instances 𝝁 test 1 superscript subscript 𝝁 test 1 \bm{{\mu}}_{\text{test}}^{1} and 𝝁 test 2 superscript subscript 𝝁 test 2 \bm{{\mu}}_{\text{test}}^{2} . In this experiment, the offline training is performed both with and without the early stopping strategy.

	${𝝁}_{\text{test}}^1=(4.3,0.021)$		${𝝁}_{\text{test}}^2=(5.15,0.0285)$
Early stopping	with	without	with	without
Manifold LSPG	$7.12\times {10}^{-4}$	$6.78\times {10}^{-4}$	$8.05\times {10}^{-4}$	$4.26\times {10}^{-4}$
Manifold Galerkin	$2.62\times {10}^{-2}$	$1.15\times {10}^{-2}$	$2.77\times {10}^{-2}$	$1.20\times {10}^{-2}$
Optimal projection error	$4.20\times {10}^{-4}$	$2.87\times {10}^{-4}$	$3.98\times {10}^{-4}$	$2.81\times {10}^{-4}$

Equations

$$d_n(\mathcal{M}):=\underset{{\mathcal{S}}_n}{\inf }\underset{f\in \mathcal{M}}{\sup }\underset{g\in {\mathcal{S}}_n}{\inf }\Vert f-g\Vert ,$$

$$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},t;\mathbf{\mu }),\mathbf{x}(0;\mathbf{\mu })={\mathbf{x}}^0(\mathbf{\mu }),$$

$${\mathbf{r}}^n({\mathbf{x}}^n;\mathbf{\mu })=0,n=1,\dots ,N_t,$$

$${\mathbf{r}}^n:(\mathbf{\xi };\mathbf{\nu })\mapsto {\alpha }_0\mathbf{\xi }-\Delta t{\beta }_0\mathbf{f}(\mathbf{\xi },t^n;\mathbf{\nu })+\sum _{j=1}^k{\alpha }_j{\mathbf{x}}^{n-j}-\Delta t\sum _{j=1}^k{\beta }_j\mathbf{f}({\mathbf{x}}^{n-j},t^{n-j};\mathbf{\nu }).$$

$$\tilde{\mathbf{x}}(t;\mathbf{\mu })={\mathbf{x}}_{\text{ref}}(\mathbf{\mu })+\mathbf{g}(\hat{\mathbf{x}}(t;\mathbf{\mu })),$$

$$\dot{\tilde{\mathbf{x}}}(t;\mathbf{\mu })=\mathbf{J}(\hat{\mathbf{x}}(t;\mathbf{\mu }))\dot{\hat{\mathbf{x}}}(t;\mathbf{\mu }),$$

$$T_{\hat{\mathbf{x}}}\mathcal{S}\equiv \text{Ran}(\mathbf{J}(\hat{\mathbf{x}})),$$

$${\mathbf{x}}_{\text{ref}}(\mathbf{\mu })={\mathbf{x}}^0(\mathbf{\mu })-\mathbf{g}({\hat{\mathbf{x}}}^0(\mathbf{\mu })).$$

$$\tilde{\mathbf{x}}(t;\mathbf{\mu })={\mathbf{x}}^0(\mathbf{\mu })+\mathbf{\Phi }\hat{\mathbf{x}}(t;\mathbf{\mu })\text{and}\dot{\tilde{\mathbf{x}}}(t;\mathbf{\mu })=\mathbf{\Phi }\dot{\hat{\mathbf{x}}}(t;\mathbf{\mu }),$$

$$\mathbf{r}:(\mathbf{v},\mathbf{\xi },\tau ;\mathbf{\nu })\mapsto \mathbf{v}-\mathbf{f}(\mathbf{\xi },\tau ;\mathbf{\nu })$$

$$\mathbf{r}\left(\dot{\mathbf{x}},\mathbf{x},t;\mathbf{\mu }\right)=0,\mathbf{x}(0;\mathbf{\mu })={\mathbf{x}}^0(\mathbf{\mu }).$$

$$\dot{\hat{\mathbf{x}}}(t;\mathbf{\mu })=\underset{\hat{\mathbf{v}}\in {\mathbb{R}}^p}{\mathrm{argmin}}{\Vert \mathbf{r}\left(\mathbf{J}(\hat{\mathbf{x}}(t;\mathbf{\mu }))\hat{\mathbf{v}},{\mathbf{x}}_{\text{ref}}(\mathbf{\mu })+\mathbf{g}(\hat{\mathbf{x}}(t;\mathbf{\mu })),t;\mathbf{\mu }\right)\Vert }_2^2$$

$$\dot{\hat{\mathbf{x}}}(t;\mathbf{\mu })=\mathbf{J}(\hat{\mathbf{x}}(t;\mathbf{\mu }){)}^{+}\mathbf{f}\left({\mathbf{x}}_{\text{ref}}(\mathbf{\mu })+\mathbf{g}(\hat{\mathbf{x}}(t;\mathbf{\mu })),t;\mathbf{\mu }\right),\hat{\mathbf{x}}(0;\mathbf{\mu })={\hat{\mathbf{x}}}^0(\mathbf{\mu }),$$

$${\hat{\mathbf{r}}}_G^n({\hat{\mathbf{x}}}^n;\mathbf{\mu })=0,$$

$$\begin{array}{rl}{\hat{\mathbf{r}}}_G^n:(\hat{\mathbf{\xi }};\mathbf{\nu })\mapsto & {\alpha }_0\hat{\mathbf{\xi }}-\Delta t{\beta }_0\mathbf{J}(\hat{\mathbf{\xi }}{)}^{+}\mathbf{f}({\mathbf{x}}_{\text{ref}}(\mathbf{\nu })+\mathbf{g}(\hat{\mathbf{\xi }}),t^n;\mathbf{\nu })+\\ & \sum _{j=1}^k{\alpha }_j{\hat{\mathbf{x}}}^{n-j}-\Delta t\sum _{j=1}^k{\beta }_j\mathbf{J}({\hat{\mathbf{x}}}^{n-j}{)}^{+}\mathbf{f}({\mathbf{x}}_{\text{ref}}(\mathbf{\nu })+\mathbf{g}({\hat{\mathbf{x}}}^{n-j}),t^{n-j};\mathbf{\nu }).\end{array}$$

$$\dot{\tilde{\mathbf{x}}}(t;\mathbf{\mu })=\underset{\mathbf{v}\in T_{\hat{\mathbf{x}}(t;\mathbf{\mu })}\mathcal{S}}{\mathrm{argmin}}{\Vert \mathbf{v}-\mathbf{f}\left({\mathbf{x}}_{\text{ref}}(\mathbf{\mu })+\mathbf{g}(\hat{\mathbf{x}}(t;\mathbf{\mu })),t;\mathbf{\mu }\right)\Vert }_2^2,$$

$$\dot{\hat{\mathbf{x}}}(t;\mathbf{\mu })=\underset{\hat{\mathbf{v}}\in {\mathbb{R}}^p}{\mathrm{argmin}}{\Vert \mathbf{r}\left(\mathbf{J}(\hat{\mathbf{x}}(t;\mathbf{\mu }))\hat{\mathbf{v}},{\mathbf{x}}_{\text{ref}}(\mathbf{\mu })+\mathbf{g}(\hat{\mathbf{x}}(t;\mathbf{\mu })),t;\mathbf{\mu }\right)\Vert }_{\mathbf{\Theta }}^2,$$

$$\dot{\hat{\mathbf{x}}}(t;\mathbf{\mu })=[{\mathbf{\Theta }}^{1/2}\mathbf{J}(\hat{\mathbf{x}}(t;\mathbf{\mu })){]}^{+}{\mathbf{\Theta }}^{1/2}\mathbf{f}\left({\mathbf{x}}_{\text{ref}}(\mathbf{\mu })+\mathbf{g}(\hat{\mathbf{x}}(t;\mathbf{\mu })),t;\mathbf{\mu }\right),\hat{\mathbf{x}}(0;\mathbf{\mu })={\hat{\mathbf{x}}}^0(\mathbf{\mu }),$$

$${\hat{\mathbf{r}}}_G^n:(\hat{\mathbf{\xi }};\mathbf{\nu })\mapsto$$

$$\sum _{j=1}^k{\alpha }_j{\hat{\mathbf{x}}}^{n-j}-\Delta t\sum _{j=1}^k{\beta }_j[{\mathbf{\Theta }}^{1/2}\mathbf{J}({\hat{\mathbf{x}}}^{n-j}){]}^{+}{\mathbf{\Theta }}^{1/2}\mathbf{f}({\mathbf{x}}_{\text{ref}}(\mathbf{\nu })+\mathbf{g}({\hat{\mathbf{x}}}^{n-j}),t^{n-j};\mathbf{\nu }).$$

$$\dot{\hat{\mathbf{x}}}(t;\mathbf{\mu })=\bar{\mathbf{g}}\left(\mathbf{f}\left({\mathbf{x}}_{\text{ref}}(\mathbf{\mu })+\mathbf{g}(\hat{\mathbf{x}}(t;\mathbf{\mu })),t;\mathbf{\mu }\right)\right),\hat{\mathbf{x}}(0;\mathbf{\mu })=\bar{\mathbf{g}}({\mathbf{x}}^0(\mathbf{\mu })).$$

$${\hat{\mathbf{x}}}^n(\mathbf{\mu })=\underset{\hat{\mathbf{\xi }}\in {\mathbb{R}}^p}{\mathrm{argmin}}{\Vert {\mathbf{r}}^n\left({\mathbf{x}}_{\text{ref}}(\mathbf{\mu })+\mathbf{g}(\hat{\mathbf{\xi }});\mathbf{\mu })\right)\Vert }_2^2,$$

$${\mathbf{\Psi }}^n({\hat{\mathbf{x}}}^n;\mathbf{\mu }{)}^T{\mathbf{r}}^n\left({\mathbf{x}}_{\text{ref}}(\mathbf{\mu })+\mathbf{g}({\hat{\mathbf{x}}}^n);\mathbf{\mu }\right)=0,$$

$${\mathbf{\Psi }}^n:(\hat{\mathbf{\xi }};\mathbf{\nu })$$

$${\hat{\mathbf{x}}}^n(\mathbf{\mu })=\underset{\hat{\mathbf{\xi }}\in {\mathbb{R}}^p}{\mathrm{argmin}}{\Vert {\mathbf{r}}^n\left({\mathbf{x}}_{\text{ref}}(\mathbf{\mu })+\mathbf{g}(\hat{\mathbf{\xi }});\mathbf{\mu })\right)\Vert }_{\mathbf{\Theta }}^2,$$

$${\mathbf{\Psi }}^n:(\hat{\mathbf{\xi }};\mathbf{\nu })$$

$$\mathbf{f}(\mathbf{x};\mathbf{\mu })=0,$$

$$\hat{\mathbf{x}}(\mathbf{\mu })=\underset{\hat{\mathbf{\xi }}\in {\mathbb{R}}^p}{\mathrm{argmin}}{\Vert \mathbf{f}\left({\mathbf{x}}_{\text{ref}}(\mathbf{\mu })+\mathbf{g}(\hat{\mathbf{\xi }});\mathbf{\mu })\right)\Vert }_2^2,$$

$$\mathbf{\Psi }(\hat{\mathbf{x}}(\mathbf{\mu });\mathbf{\mu }{)}^T\mathbf{f}\left({\mathbf{x}}_{\text{ref}}(\mathbf{\mu })+\mathbf{g}(\hat{\mathbf{x}}(\mathbf{\mu }));\mathbf{\mu }\right)=0,$$

$$\mathbf{\Psi }:(\hat{\mathbf{\xi }};\mathbf{\nu })$$

$$\frac{\partial {\hat{\mathbf{r}}}_G^n}{\partial {\hat{\xi }}_i}:(\hat{\mathbf{w}};\mathbf{\nu })\mapsto {\alpha }_0{\hat{\mathbf{e}}}_i-\Delta t{\beta }_0\left(\frac{\partial {\mathbf{J}}^{+}}{\partial {\hat{\xi }}_i}\left(\hat{\mathbf{w}}\right)\mathbf{f}({\mathbf{x}}_{\text{ref}}(\mathbf{\nu })+\mathbf{g}(\hat{\mathbf{w}}),t^n;\mathbf{\nu })+\mathbf{J}(\hat{\mathbf{w}}{)}^{+}\frac{\partial \mathbf{f}}{\partial \mathbf{\xi }}({\mathbf{x}}_{\text{ref}}(\mathbf{\nu })+\mathbf{g}(\hat{\mathbf{w}}),t^n;\mathbf{\nu })\mathbf{J}(\hat{\mathbf{w}}){\hat{\mathbf{e}}}_i\right)$$

Full text

url]sandia.gov/ ktcarlb

1 Introduction

Physics-based modeling and simulation has become indispensable across many applications in engineering and science, ranging from aircraft design to monitoring national critical infrastructure. However, as simulation is playing an increasingly important role in scientific discovery, decision making, and design, greater demands are being placed on model fidelity. Achieving high fidelity often necessitates including fine spatiotemporal resolution in computational models of the system of interest; this can lead to very large-scale models whose simulations consume months on thousands of computing cores. This computational burden precludes the integration of such high-fidelity models in important scenarios that are real time or many query in nature, as these scenarios require the (parameterized) computational model to be simulated very rapidly (e.g., model predictive control) or thousands of times (e.g., uncertainty propagation).

Projection-based reduced-order models (ROMs) provide one approach for overcoming this computational burden. These techniques comprise two stages: an offline stage and an online stage. During the offline stage, these methods perform computationally expensive training tasks (e.g., simulating the high-fidelity model at several points in the parameter space) to compute a representative low-dimensional ‘trial’ subspace for the system state. Next, during the inexpensive online stage, these methods rapidly compute approximate solutions for different points in the parameter space via projection: they compute solutions in the low-dimensional trial subspace by enforcing the high-fidelity-model residual to be orthogonal to a low-dimensional test subspace of the same dimension.

As suggested above, nearly all projection-based model-reduction approaches employ linear trial subspaces. This includes the reduced-basis technique [62, 67] and proper orthogonal decomposition (POD) [37, 15] for parameterized stationary problems; balanced truncation [54], rational interpolation [4, 32], and Craig–Bampton model reduction [20] for linear time invariant (LTI) systems; and Galerkin projection [37], least-squares Petrov–Galerkin projection [13], and other Petrov–Galerkin projections [78] with (balanced) POD [37, 44, 78, 66] for nonlinear dynamical systems.

The Kolmogorov $n$-width [61] provides one way to quantify the optimal linear trial subspace; it is defined as

[TABLE]

where the first infimum is taken over all $n$-dimensional subspaces of the state space, and $\mathcal{M}$ denotes the manifold of solutions over all time and parameters. Assuming the dynamical system has a unique trajectory for each parameter instance, the intrinsic solution-manifold dimensionality is (at most) equal to the number of parameters plus one (time). For problems that exhibit a fast decaying Kolmogorov $n$-width (e.g., diffusion-dominated problems), employing a linear trial subspace is theoretically justifiable [3, 6] and has enjoyed many successful demonstrations. Unfortunately, many computational problems exhibit slowly decaying Kolmogorov $n$-width (e.g., advection-dominated problems). In such cases, the use of low-dimensional linear trial subspaces often produces inaccurate results; the ROM dimensionality must be significantly increased to yield acceptable accuracy [58]. Indeed, the Kolmogorov $n$-width with $n$ equal to the intrinsic solution-manifold dimensionality is often quite large for such problems.

Several approaches have been pursued to address this $n$-width limitation of linear trial subspaces. One set of approaches transforms the trial basis to improve its approximation properties for advection-dominated problems. Such methods include separating transport dynamics via ‘freezing’ [57], applying a coordinate transformation to the trial basis [40, 56, 10], shifting the POD basis [63], transforming the physical domain of the snapshots [77, 76], constructing the trial basis on a Lagrangian formulation of the governing equations [53], and using Lax pairs of the Schrödinger operator to construct a time-evolving trial basis [27]. Other approaches pursue the use of multiple linear subspaces instead of employing a single global linear trial subspace; these local subspaces can be tailored to different regions of the time domain [24, 22], physical domain [71], or state space [2, 60]. Ref. [11] aims to overcome the limitations of using a linear trial subspace by providing online-adaptive $h$-refinement mechanism that constructs a hierarchical sequence of linear subspaces that converges to the original state space. However, all of these methods attempt to construct, manually transform, or refine a linear basis to be locally accurate; they do not consider nonlinear trial manifolds of more general structure. Further, many of these approaches rely on substantial additional knowledge about the problem, such as the particular advection phenomena governing basis shifting.

This work aims to address the fundamental $n$-width deficiency of linear trial subspaces. However, in contrast to the above methods, we pursue an approach that is both more general (i.e., it should not be limited to piecewise linear manifolds or mode transformations) and only requires the same snapshot data as typical POD-based approaches (e.g., it should require no special knowledge about any particular advection phenomena). To accomplish this objective, we propose an approach that (1) performs optimal projection of dynamical systems onto arbitrary nonlinear trial manifolds (during the online stage), and (2) computes this nonlinear trial manifold from snapshot data alone (during the offline stage).

For the first component, we perform optimal projection onto arbitrary (continuously differentiable) nonlinear trial manifolds by applying minimum-residual formulations at the time-continuous (ODE) and time-discrete (O$\Delta$E) levels. The time-continuous formulation leads to manifold Galerkin projection, which can be interpreted as performing orthogonal projection of the velocity onto the tangent space of the trial manifold. The time-discrete formulation leads to manifold least-squares Petrov–Galerkin (LSPG) projection, which can also be straightforwardly extended to stationary (i.e., steady-state) problems. We also perform analyses that illustrate the relationship between these manifold ROMs and classical linear-subspace ROMs. Manifold Galerkin and manifold LSPG projection require the trial manifold to be characterized as a (generally nonlinear) mapping from the low-dimensional reduced state to the high-dimensional state; the mapping from the high-dimensional state to the low-dimensional state is not required.

The second component aims to compute a nonlinear trial manifold from snapshot data alone. Many machine-learning methods exist to perform nonlinear dimensionality reduction. However, many of these methods do not provide the required mapping from the low-dimensional embedding to the high-dimensional input; examples include Isomap [73], locally linear embedding (LLE) [65], Hessian eigenmaps [23], spectral embedding [5], and t-SNE [49]. Methods that do provide this required mapping include self-organizing maps [43], generative topographic mapping [7], kernel principal component analysis (PCA) [70], Gaussian process latent variable model [45], diffeomorphic dimensionality reduction [75], and autoencoders [35]. In principle, manifold Galerkin and manifold LSPG projection could be applied with manifolds constructed by any of the methods in the latter category. However, this study restricts focus to autoencoders—more specifically deep convolutional autoencoders—due to their expressiveness and scalability, as well as the availability of high-performance software tools for their construction.

Autoencoders (also known as auto-associators [21]) comprise a specific type of feedforward neural network that aim to learn the identity mapping: they attempt to copy the input to an accurate approximation of itself. Learning the identity mapping is not a particularly useful task unless, however, it associates with a dimensionality-reduction procedure comprising data compression and subsequent recovery. This is precisely what autoencoders accomplish by employing a neural-network architecture consisting of two parts: an encoder that provides a nonlinear mapping from the high-dimensional input to a low-dimensional embedding, and a decoder that provides a nonlinear mapping from the low-dimensional embedding to an approximation of the high-dimensional input. Convolutional autoencoders are a specific type of autoencoder that employ convolutional layers, which have been shown to be effective for extracting representative features in images [51]. Inspired by the analogy between images and spatially distributed dynamical-system states (e.g., when the dynamical system corresponds to the spatial discretization of a partial-differential-equations model), we propose a specific deep convolutional autoencoder architecture tailored to dynamical systems with states that are spatially distributed. Critically, training this autoencoder requires only the same snapshot data as POD; no additional problem-specific information is needed.

In summary, new contributions of this work include:

Manifold Galerkin (Section 3.2) and manifold LSPG (Section 3.3) projection techniques, which project the dynamical-system model onto arbitrary continuously-differentiable manifolds. We equip these methods with

1. (a)

the ability to exactly satisfy the initial condition (Remark 3.1), and 2. (b)

quasi-Newton solvers (Section 3.4) to solve the system of algebraic equations arising from implicit time integration. 2. 2.

Analysis (Section 4), which includes

1. (a)

demonstrating that employing an affine trial manifold recovers classical linear-subspace Galerkin and LSPG projection (Proposition 4.1), 2. (b)

sufficient conditions for commutativity of time discretization and manifold Galerkin projection (Theorem 4.1), 3. (c)

conditions under which manifold Galerkin and manifold LSPG projection are equivalent (Theorem 4.2), and 4. (d)

a posteriori discrete-time error bounds for both the manifold Galerkin and manifold LSPG projection methods (Theorem 4.3). 3. 3.

A novel convolutional autoencoder architecture tailored to spatially distributed dynamical-system states (Section 5.2) with accompanying offline training algorithm that requires only the same snapshot data as POD (Section 6). 4. 4.

Numerical experiments on advection-dominated benchmark problems (Section 7). These experiments illustrate the ability of the method to outperform even the projection of the solution onto the optimal linear subspace; further, the proposed method is close to achieving the optimal performance of any nonlinear-manifold method. This demonstrates the method’s ability to overcome the intrinsic $n$-width limitations of linear trial subspaces.

We note that the methodology is applicable to both linear and nonlinear dynamical systems.

To the best of our knowledge, Refs. [33, 41] comprise the only attempts to incorporate an autoencoder within a projection-based ROM. These methods seek solutions in the nonlinear trial manifold provided by an autoencoder; however, these methods reduce the number of equations by applying the encoder to the velocity. Unfortunately, as we discuss in Remark 3.5, this approach is kinematically inconsistent, as the velocity resides in the tangent space to the manifold, not the manifold itself. Thus, encoding the velocity can produce significant approximation errors. Instead, the proposed manifold Galerkin and manifold LSPG projection methods produce approximations that associate with minimum-residual formulations and adhere to the kinematics imposed by the trial manifold.

Relatedly, Ref. [31] proposes a general framework for projection of dynamical systems onto nonlinear manifolds. However, the proposed method constructs a piecewise linear trial manifold by generating local linear subspaces and concatenating those subspaces. Then, the method projects the residual of the governing equations onto a nonlinear test manifold that is also piecewise linear; this is referred to as a ‘piecewise linear projection function’. Thus, the approach is limited to piecewise-linear manifolds, and the resulting approximation does not associate with any optimality property.

We also note that autoencoders have been applied to various non-intrusive model-reduction methods that are purely data driven in nature and are not based on a projection process. Examples include Ref. [52], which constructs a nonlinear model of wall turbulence using an autoencoder; Ref. [29], which applies an autoencoder to compress the state, followed by a recurrent neural network (RNN) [68] to learn the dynamics; Refs. [72, 59, 48, 55], which apply autoencoders to learn approximate invariant subspaces of the Koopman operator; and Ref. [17], which applies hierarchical dimensionality reduction comprising autoencoders and PCA followed by dynamics learning to recover missing CFD data.

The paper is organized as follows. Section 2 describes the full-order model, which corresponds to a parameterized system of (linear or nonlinear) ordinary differential equations. Section 3 describes model reduction on nonlinear manifolds, including the mathematical characterization of the nonlinear trial manifold (Section 3.1), manifold Galerkin projection (Section 3.2), manifold LSPG projection (Section 3.3), and associated quasi-Newton methods to solve the system of algebraic equations arising at each time instance in the case of implicit time integration (Section 3.4). Section 4 provides the aforementioned analysis results. Section 5 describes a practical approach for constructing the nonlinear trial manifold using deep convolutional autoencoders, including a brief description of autoencoders (Section 5.1), the proposed autoencoder architecture applicable to spatially distributed states (Section 5.2), and the way in which the proposed autoencoder can be used to satisfy the initial condition (Section 5.3). When the manifold Galerkin and manifold LSPG ROMs employ this choice of decoder, we refer to them as Deep Galerkin and Deep LSPG ROMs, respectively. Section 6 describes offline training, which entails snapshot-based data collection (Section 6.1), data standardization (Section 6.2), and autoencoder training (Section 6.3); Algorithm 1 summarizes the offline training stage. Section 7 assesses the performance of the proposed Deep Galerkin and Deep LSPG ROMs compared to (linear-subspace) POD–Galerkin and POD–LSPG ROMs on two advection-dominated benchmark problems. Finally, Section 8 concludes the paper.

2 Full-order model

This work considers the full-order model (FOM) to correspond to a dynamical system expressed as a parameterized system of ordinary differential equations (ODEs)

[TABLE]

where $t\in[0,T]$ denotes time with final time $T\in\mathbb{R}_{+}$, and $\bm{x}:[0,T]\times\mathcal{D}\rightarrow\mathbb{R}^{N}$ denotes the time-dependent, parameterized state implicitly defined as the solution to problem (2.1) given the parameters $\bm{{\mu}}\in\mathcal{D}$. Here, $\mathcal{D}\subseteq\mathbb{R}^{{n_{\mu}}}$ denotes the parameter space, $\bm{x}^{0}:\mathcal{D}\rightarrow\mathbb{R}^{N}$ denotes the parameterized initial condition, and $\bm{f}:\mathbb{R}^{N}\times[0,T]\times\mathcal{D}\rightarrow\mathbb{R}^{N}$ denotes the velocity, which may be linear or nonlinear in its first argument. We denote differentiation of a variable $x$ with respect to time by $\dot{x}$. Such dynamical systems may arise from the semidiscretization of a partial-differential-equations (PDE) model, for example. We refer to Eq. (2.1) as the FOM ODE.

Numerically solving the FOM ODE (2.1) requires application of a time-discretization method. For simplicity, this work restricts attention to linear multistep methods; see Ref. [12] for analysis of both (linear-subspace) Galerkin and LSPG reduced-order models with Runge–Kutta schemes. A linear $k$-step method applied to numerically solve the FOM ODE (2.1) leads to solving the system of algebraic equations

[TABLE]

where the time-discrete residual $\bm{r}^{n}:\mathbb{R}^{N}\times\mathcal{D}\rightarrow\mathbb{R}^{N}$ is defined as

[TABLE]

Here, $\Delta t\in\mathbb{R}_{+}$ denotes the time step, $\bm{x}^{k}$ denotes the numerical approximation to $\bm{x}(k\Delta t;\bm{{\mu}})$, and the coefficients $\alpha_{j}$ and $\beta_{j}$, $j=0,\ldots,k$ with $\sum_{j=0}^{k}\alpha_{j}=0$ define a particular linear multistep scheme. These methods are implicit if $\beta_{0}\neq 0$. For notational simplicity, we have assumed a uniform time step $\Delta t$ and a fixed number of steps $k$ for each time instance; the extensions to nonuniform grids and a non-constant value of $k$ are straightforward. We refer to Eq. (2.2) as the FOM O$\Delta$E.

Remark 2.1** **(Intrinsic solution-manifold

dimensionality).

Assuming the initial value problem (2.2) has a unique solution for each parameter instance $\bm{{\mu}}\in\mathcal{D}$, the intrinsic dimensionality of the solution manifold $\{\bm{x}(t;\bm{{\mu}})\;|\;t\in[0,T],\ \bm{{\mu}}\in\mathcal{D}\}$ is (at most) $p^{\star}={n_{\mu}}+1$, as the mapping $(t;\bm{{\mu}})\mapsto\bm{x}$ is unique in this case. This provides a practical lower bound on the dimension of a nonlinear trial manifold for exactly representing the dynamical-system state.

3 Model reduction on nonlinear manifolds

This section proposes two classes of residual-minimizing ROMs on nonlinear manifolds. The first minimizes the (time-continuous) FOM ODE residual and is analogous to classical Galerkin projection, while the second minimizes the (time-discrete) FOM O$\Delta$E residual and is analogous to least-squares Petrov–Galerkin (LSPG) projection [13, 16, 12]. Section 3.1 introduces the notion of a nonlinear trial manifold, Section 3.2 describes the manifold Galerkin ROM resulting from time-continuous residual minimization, and Section 3.3 describes the manifold LSPG ROM resulting from time-discrete residual minimization.

3.1 Nonlinear trial manifold

We begin by prescribing the subset of the original state space $\mathbb{R}^{N}$ on which the ROM techniques will seek approximate solutions. Rather than introducing a classical affine trial subspace for this purpose, we instead introduce a continuously differentiable nonlinear trial manifold. In particular, we seek approximate solutions $\tilde{\bm{x}}\approx\bm{x}$ of the form

[TABLE]

where \tilde{\bm{x}}:\mathbb{R}_{+}\times\mathcal{D}\rightarrow{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{x}_{\text{ref}}(\bm{{\mu}})+}}\mathcal{S} and $\mathcal{S}\vcentcolon=\{\bm{g}(\hat{\bm{\xi}})\,|\,\hat{\bm{\xi}}\in\mathbb{R}^{p}\}$ denotes the nonlinear trial manifold from the extrinsic view. Here, $\bm{x}_{\text{ref}}:\mathcal{D}\rightarrow\mathbb{R}^{N}$ denotes a parameterized reference state and $\bm{g}:\hat{\bm{\xi}}\mapsto\bm{g}(\hat{\bm{\xi}})$ with $\bm{g}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{N}$ and $p\leq N$ denotes the parameterization function, which we refer to in this work as the decoder, which comprises a nonlinear mapping from low-dimensional generalized coordinates $\hat{\bm{x}}:\mathbb{R}_{+}\times\mathcal{D}\rightarrow\mathbb{R}^{p}$ to the high-dimensional state approximation. From application of the chain rule, the approximated velocity is then

[TABLE]

where $\bm{J}:\hat{\bm{\xi}}\mapsto\frac{d\bm{g}}{d\hat{\bm{\xi}}}(\hat{\bm{\xi}})$ with $\bm{J}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{N\times p}$ denotes the Jacobian of the decoder. This Jacobian defines the tangent space to the manifold

[TABLE]

where $\text{Ran}(\bm{A})$ denotes the range of matrix $\bm{A}$. From Remark 2.1, we observe that the nonlinear trial manifold dimension must be greater than or equal to the intrinsic solution-manifold dimension, i.e., $p\geq p^{\star}$, to be able to exactly represent the state.

Remark 3.1** (Initial-condition satisfaction).**

Satisfaction of the initial conditions requires the initial generalized coordinates $\hat{\bm{x}}(0;\bm{{\mu}})=\hat{\bm{x}}^{0}(\bm{{\mu}})$ to satisfy $\bm{g}(\hat{\bm{x}}^{0}(\bm{{\mu}}))=\bm{x}^{0}(\bm{{\mu}})$. This can be achieved for any choice of $\hat{\bm{x}}^{0}(\bm{{\mu}})$ provided the reference state is set to

[TABLE]

However, doing so must ensure that the decoder can accurately represent deviations from this reference state (3.4). In the context of the proposed autoencoder-based trial manifold, Section 5.3 describes a strategy for computing the initial generalized coordinates $\hat{\bm{x}}^{0}(\bm{{\mu}})$ and defining training data for manifold construction such that the decoder accomplishes this.

Remark 3.2** (Linear trial subspace).**

The classical linear (or more precisely affine) trial subspace corresponds to a specific case of a nonlinear trial manifold; this occurs when the decoder is linear. In this case, the decoder can be expressed as $\bm{g}:\hat{\bm{\xi}}\mapsto\bm{\Phi}\hat{\bm{\xi}}$, where $\bm{\Phi}\in\mathbb{R}_{\star}^{N\times p}$ denotes trial-basis matrix and $\mathbb{R}_{\star}^{n\times m}$ denotes the set of full-column-rank $n\times m$ matrices (the non-compact Stiefel manifold). In this case, the approximated state and velocity can be expressed as

[TABLE]

respectively, such that the trial manifold is affine, i.e., \mathcal{S}={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\text{Ran}(\bm{\Phi})}}, and the decoder Jacobian is the constant matrix $\bm{J}(\hat{\bm{x}}(t;\bm{{\mu}}))=\bm{J}=\bm{\Phi}$. Note that this approach can enforce the initial condition by setting the initial generalized coordinates $\hat{\bm{x}}^{0}(\bm{{\mu}})$ to zero and subsequently setting $\bm{x}_{\text{ref}}(\bm{{\mu}})=\bm{x}^{0}(\bm{{\mu}})$. Common choices for computing the reduced-basis matrix $\bm{\Phi}$ when the velocity is linear in its first argument include balanced truncation [54], rational interpolation [4, 32], the reduced-basis method [62, 67], Rayleigh–Ritz eigenvectors [20]; common choices when the velocity is nonlinear include POD [36] and balanced POD [44, 78, 66].

3.2 Manifold Galerkin ROM: time-continuous residual minimization

We now derive the ROM corresponding to time-continuous residual minimization. To do so, define the FOM ODE residual as

[TABLE]

with $\bm{r}:\mathbb{R}^{N}\times\mathbb{R}^{N}\times\mathbb{R}_{+}\times\mathcal{D}\rightarrow\mathbb{R}^{N}$. It can be easily verified that the FOM ODE (2.1) is equivalent to

[TABLE]

To derive the time-continuous residual-minimizing ROM, we substitute the approximated state $\bm{x}\leftarrow\tilde{\bm{x}}$ defined in Eq. (3.1) and velocity $\dot{\bm{x}}\leftarrow\dot{\tilde{\bm{x}}}$ defined in Eq. (3.2) into the FOM ODE (3.6) and minimize the $\ell^{2}$-norm of the residual, i.e.,

[TABLE]

with initial condition $\hat{\bm{x}}(0;\bm{{\mu}})=\hat{\bm{x}}^{0}(\bm{{\mu}})$. If the Jacobian $\bm{J}(\hat{\bm{x}}(t;\bm{{\mu}}))$ has full column rank, then Problem (3.7) is convex and has a unique solution that leads to the manifold Galerkin ODE

[TABLE]

where the superscript $+$ denotes the Moore–Penrose pseudoinverse.

We note that although the state evolves on the nonlinear manifold $\mathcal{S}$, the generalized coordinates evolve in the Euclidean space $\mathbb{R}^{p}$, which facilitates time integration of the manifold Galerkin ODE (3.8). Indeed, this can be accomplished using any time integrator. If a linear multistep scheme is employed, then the manifold Galerkin O$\Delta$E corresponds to

[TABLE]

where the manifold Galerkin O$\Delta$E residual is

[TABLE]

We note that the manifold Galerkin O$\Delta$E is nonlinear if either the velocity is nonlinear in its first argument or if the trial manifold is nonlinear.

Remark 3.3** **(Manifold Galerkin projection as orthogonal projection onto the

tangent space).

Eq. (3.7) can be written equivalently as

[TABLE]

where the tangent space is defined in Eq. (3.3). The solution to this minimization problem is $\dot{\tilde{\bm{x}}}(t;\bm{{\mu}})=\bm{J}(\hat{\bm{x}}(t;\bm{{\mu}}))\bm{J}(\hat{\bm{x}}(t;\bm{{\mu}}))^{+}\bm{f}\left(\bm{x}_{\text{ref}}(\bm{{\mu}})+\bm{g}(\hat{\bm{x}}(t;\bm{{\mu}})),t;\bm{{\mu}}\right)$. Thus, Manifold Galerkin projection can be interpreted as computing the orthogonal (i.e., $\ell^{2}$-optimal) projection of the velocity onto the tangent space $T_{\hat{\bm{x}}(t;\bm{{\mu}})}\mathcal{S}$ of the trial manifold.

Remark 3.4** (Weighted $\ell^{2}$-norms for manifold Galerkin projection).**

We note that the minimization problem (3.7) could be posed in other norms. For example, if we consider residual minimization in the weighted $\ell^{2}$-norm $\|\cdot\|_{\bm{\Theta}}$ satisfying $\|\bm{x}\|_{\bm{\Theta}}\equiv(\bm{x},\bm{x})_{\bm{\Theta}}$ with $(\bm{x},\bm{y})_{\bm{\Theta}}\equiv\sqrt{\bm{x}^{T}{\bm{\Theta}}\bm{y}}$ and ${\bm{\Theta}}\equiv\bm{\Theta}^{T/2}\bm{\Theta}^{1/2}$, then Problem (3.7) becomes

[TABLE]

leading to the manifold Galerkin ODE

[TABLE]

and manifold Galerkin O$\Delta$E residual

[TABLE]

For notational simplicity, this work restricts attention to the $\ell^{2}$-norm; future work will consider such weighted norms, e.g., to enable hyper-reduction.

Remark 3.5** (Alternative Galerkin projection).**

Refs. [33, 41] provide an alternative way to perform Galerkin projection on nonlinear manifolds constructed from autoencoders. In contrast to the proposed manifold Galerkin ROM (and manifold LSPG ROM introduced in Section 3.3 below), these methods additionally require the existence of an encoder $\bar{\bm{g}}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{p}$ that satisfies $\bm{g}\circ\bar{\bm{g}}:\bm{x}\mapsto\tilde{\bm{x}}$ with $\tilde{\bm{x}}\approx\bm{x}$, which is trained on state-snapshot data. These methods then form the low-dimensional system of ODEs

[TABLE]

One problem with this approach is that it implicitly assumes that the velocity $\bm{f}$ can be represented on the manifold used to represent the state. This is kinematically inconsistent, as the velocity resides in the tangent space to this manifold as derived in Eq. (3.2), i.e., $\dot{\tilde{\bm{x}}}(t;\bm{{\mu}})\in T_{\hat{\bm{x}}(t;\bm{{\mu}})}\mathcal{S}$; the tangent space and the manifold coincide if and only if the trial manifold associates with a linear trial subspace. Thus, encoding the velocity is likely to produce a poor approximation of the reduced-state velocity $\dot{\hat{\bm{x}}}(t;\bm{{\mu}})$. Instead, the proposed manifold Galerkin ROM performs orthogonal projection of the velocity onto the tangent space (see Remark 3.3). We also note the approaches proposed in Refs. [33, 41] compute the initial state by encoding the initial condition, which is not guaranteed to satisfy the initial conditions as discussed in Remark 3.1.

3.3 Manifold LSPG ROM: time-discrete residual minimization

Analogously to the previous section, we now derive the ROM corresponding to time-discrete residual minimization. To do so, we simply substitute the approximated state $\bm{x}\leftarrow\tilde{\bm{x}}$ defined in Eq. (3.1) into the FOM O$\Delta$E (2.2) and minimize the $\ell^{2}$-norm of the residual, i.e.,

[TABLE]

which is solved sequentially for $n=1,\ldots,{N_{t}}$ with initial condition $\hat{\bm{x}}(0;\bm{{\mu}})=\hat{\bm{x}}^{0}(\bm{{\mu}})$. Necessary optimality conditions for Problem (3.17) correspond to stationarity of the objective function, i.e., the solution $\hat{\bm{x}}^{n}$ satisfies the manifold LSPG O$\Delta$E

[TABLE]

where the test basis matrix $\bm{\Psi}^{n}:\mathbb{R}^{p}\times\mathcal{D}\rightarrow\mathbb{R}^{N\times p}$ is defined as

[TABLE]

As with the manifold Galerkin O$\Delta$E, the manifold LSPG O$\Delta$E is nonlinear if either the velocity is nonlinear in its first argument or if the trial manifold is nonlinear.

Remark 3.6** (Weighted $\ell^{2}$-norms for manifold LSPG projection).**

As discussed in Remark 3.4, the minimization problem in (3.17) could be posed in other norms. If the weighted $\ell^{2}$-norm $\|\cdot\|_{\bm{\Theta}}$ is employed, then Problem (3.17) becomes

[TABLE]

and the test basis matrix in the manifold LSPG O$\Delta$E(3.18) becomes

[TABLE]

As with manifold Galerkin projection, this work restricts attention to the $\ell^{2}$-norm for notational simplicity.

Remark 3.7** (Application to stationary problems).**

Manifold LSPG projection can also be applied to stationary (i.e., steady-state) problems. Stationary problems are characterized by computing $\bm{x}(\bm{{\mu}})$ as the implicit solution to Eq. (2.1) with $\dot{\bm{x}}=\bm{0}$, i.e.,

[TABLE]

where the (time-independent) velocity is equivalent to the stationary-problem residual. Substituting the approximated state $\bm{x}(\bm{{\mu}})\leftarrow\tilde{\bm{x}}(\bm{{\mu}})=\bm{x}_{\text{ref}}(\bm{{\mu}})+\bm{g}(\hat{\bm{x}}(\bm{{\mu}}))$ into Eq. (3.21) and subsequently minimizing the $\ell^{2}$-norm of the residual yields

[TABLE]

which has the same nonlinear-least-squares form as manifold LSPG projection applied to the dynamical system model (3.17). In this case, necessary optimality conditions for Problem (3.22) correspond to stationarity of the objective function, i.e., the solution $\hat{\bm{x}}(\bm{{\mu}})$ satisfies the system of algebraic equations

[TABLE]

where the test basis matrix $\bm{\Psi}:\mathbb{R}^{p}\times\mathcal{D}\rightarrow\mathbb{R}^{N\times p}$ is defined as

[TABLE]

Note that we do not consider a manifold Galerkin projection applied to stationary problems, as the manifold Galerkin ROM was derived by minimizing the time-continuous residual, and a Galerkin-like projection of the form $\bm{J}(\hat{\bm{x}}(\bm{{\mu}}))^{{T}}\bm{f}\left(\bm{x}_{\text{ref}}(\bm{{\mu}})+\bm{g}(\hat{\bm{x}}(\bm{{\mu}}));\bm{{\mu}}\right)=\bm{0}$ does not generally associate with any optimality property.

3.4 Quasi-Newton methods for implicit integrators

When an implicit time integrator is employed such that $\beta_{0}\neq 0$ and the trial manifold is nonlinear, then solving the manifold Galerkin O$\Delta$E (3.10) and manifold LSPG O$\Delta$E (3.18) using Newton’s method is challenging, as the residual Jacobians involve high-order derivatives of the decoder. For this purpose, we propose quasi-Newton methods that approximate these Jacobians while retaining convergence of the nonlinear solver to the solution of the O$\Delta$Es under certain conditions.

Manifold Galerkin. We first consider the manifold Galerkin case. If an implicit integrator is employed (i.e., $\beta_{0}\neq 0$), then the manifold Galerkin O$\Delta$E (3.9) corresponds to a system of algebraic equations. Applying Newton’s method to solve this system can be challenging for nonlinear decoders, as the Jacobian of the residual requires computing a third-order tensor of second derivatives in this case. Indeed, the $i$th column of the Jacobian is

[TABLE]

for $i=1,\ldots,p$, where $\hat{\bm{e}}_{i}\in\{0,1\}^{p}$ denotes the $i$th canonical unit vector. The gradient of the pseudo-inverse of the decoder Jacobian, which appears in the second term of the right-hand side, can be computed from the gradients of this Jacobian as

[TABLE]

The primary difficulty of computing this term is the requirement of computing second derivatives of the form ${\partial\bm{J}}/{\partial\hat{{\xi}}_{i}}={\partial^{2}\bm{g}}/{\partial\hat{\bm{\xi}}\partial\hat{{\xi}}_{i}}$. For this reason, we propose to approximate the residual Jacobian as $\tilde{\bm{J}}^{n}_{\hat{r}_{G}}\approx{\partial\hat{\bm{r}}_{G}}/{\partial\hat{\bm{\xi}}}$ by neglecting this term such that

[TABLE]

The resulting quasi-Newton method to solve the manifold Galerkin O$\Delta$E (3.9) takes the form

[TABLE]

for $k=0,\ldots,K$. Here, ${\hat{\bm{x}}}^{n(0)}$ is the initial guess (often taken to be $\hat{\bm{x}}^{n-1}$ or a data-driven extrapolation [18]), and ${\alpha}^{n(k)}\in\mathbb{R}$ is the step length. Convergence of the quasi-Newton iterations (3.24) to a root of $\hat{\bm{r}}_{G}(\cdot;\bm{{\mu}})$ is ensured if

[TABLE]

and ${\alpha}^{n(k)}\in\mathbb{R}$, $k=0,\ldots,K$ is computed to satisfy the strong Wolfe conditions, as this ensures that the search direction ${\bm{p}}^{n(k)}$ is a descent direction of the function $\|\hat{\bm{r}}_{G}(\hat{\bm{\xi}};\bm{\nu})\|_{2}^{2}$ evaluated at $(\hat{\bm{\xi}};\bm{\nu})=({\hat{\bm{x}}}^{n(k)};\bm{{\mu}})$ for $k=0,\ldots,K$ [50].

Manifold LSPG. We now consider the manifold LSPG case. As in previous works that considered applying LSPG on affine trial subspaces [13, 16, 12], we propose to solve the nonlinear least-squares problem (3.17) using the Gauss–Newton method, which leads to the iterations

[TABLE]

for $k=0,\ldots,K$ with ${\hat{\bm{x}}}^{n(0)}$ a provided initial guess, and ${\alpha}^{n(k)}\in\mathbb{R}$ is a step length chosen to satisfy the strong Wolfe conditions, which ensure global convergence to a local minimum of the manifold LSPG objective function in (3.17) that satisfies the stationarity conditions associated with the LSPG O$\Delta$E (3.18).

To show the connection between the Gauss–Newton method and the quasi-Newton approach proposed for the manifold Galerkin method, we define the discrete residual associated with manifold LSPG projection as

[TABLE]

such that the manifold LSPG O$\Delta$E (3.18) can be expressed simply as the system of algebraic equations

[TABLE]

Solving Eq. (3.27) using Newton’s method requires computing the residual Jacobian, whose $i$th column is

[TABLE]

Thus, it is clear that the Gauss–Newton iterations (3.25) are equivalent to employing a quasi-Newton method to solve the manifold LSPG O$\Delta$E (3.27) with an approximated residual Jacobian $\tilde{\bm{J}}^{n}_{\hat{r}_{L}}\approx{\partial\hat{\bm{r}}_{L}^{n}}/{\partial\hat{\bm{\xi}}}$ defined as $\tilde{\bm{J}}^{n}_{\hat{r}_{L}}:(\hat{\bm{w}};\bm{\nu})\mapsto\bm{\Psi}^{n}(\hat{\bm{w}};\bm{\nu})^{{T}}\bm{\Psi}^{n}(\hat{\bm{w}};\bm{\nu})$, which is obtained from neglecting the term containing second derivatives of the decoder, namely ${\partial\bm{J}}/{\partial\hat{{\xi}}_{i}}={\partial^{2}\bm{g}}/{\partial\hat{\bm{\xi}}\partial\hat{{\xi}}_{i}}$.

4 Analysis

We now perform analysis of the proposed manifold Galerkin and manifold LSPG projection methods. For notational simplicity, this section omits the dependence of operators on parameters $\bm{{\mu}}$.

Proposition 4.1** **(Affine trial manifold recovers classical Galerkin and

LSPG).

If the trial manifold is affine, then manifold LSPG projection is equivalent to classical linear-subspace LSPG projection. If additionally the decoder mapping associates with an orthgonal matrix, then manifold Galerkin projection is equivalent to classical linear-subspace Galerkin projection.

Proof.

If the trial manifold is affine, then the decoder can be expressed as the linear mapping

[TABLE]

with $\bm{\Phi}\in\mathbb{R}^{N\times p}$.

Manifold Galerkin. Substituting the linear decoder (4.1) into the manifold Galerkin ODE (3.8) yields

[TABLE]

Similarly, substituting the linear decoder (4.1) into the manifold Galerkin O$\Delta$E Eq. (3.9) yields the same system, but with the residual defined as

[TABLE]

If additionally the trial-basis matrix is orthogonal, i.e., then $\bm{\Phi}^{+}=\bm{\Phi}^{{T}}$ and the ODE (4.2) is equivalent to the classical Galerkin ODE

[TABLE]

while the O$\Delta$E residual (4.3) is equivalent to the classical Galerkin O$\Delta$E residual

[TABLE]

Manifold LSPG. Substituting the linear decoder (4.1) into the manifold LSPG O$\Delta$E (3.18) yields the classical LSPG O$\Delta$E

[TABLE]

with classical LSPG test basis

[TABLE]

∎

The manifold Galerkin O$\Delta$E (3.9)–(3.10) was derived using a project-then-discretize approach, as we formulated the residual-minimization problem (3.7) at the time-continuous level, and subsequently discretized the equivalent ODE (3.8) in time. We now derive conditions under which an equivalent model could have been derived using a discretize-then-project approach, i.e., by substituting the approximated state $\bm{x}\leftarrow\tilde{\bm{x}}$ defined in Eq. (3.1) into the FOM O$\Delta$E (2.2) and subsequently premultiplying the overdetermined system by a test basis.

Theorem 4.1** (Commutativity of time discretization and manifold Galerkin projection).**

For manifold Galerkin projection, time discretization and projection are commutative if either (1) the trial manifold corresponds to an affine subspace, or (2) the nonlinear trial manifold is twice continuously differentiable; $\|\hat{\bm{x}}^{n-j}-\hat{\bm{x}}^{n}\|=O(\Delta t)$ for all n\in{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\{1,\ldots,{N_{t}}\}}}, j\in{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\{1,\ldots,k\}}}; and the limit $\Delta t\rightarrow 0$ is taken.

Proof.

The discrete residual characterizing the manifold Galerkin O$\Delta$E (3.9), which was derived using a project-then-discretize approach, is defined as $\hat{\bm{r}}_{G}^{n}$ in Eq. (3.10).

Substituting the approximated state $\bm{x}\leftarrow\tilde{\bm{x}}$ defined in Eq. (3.1) into the FOM O$\Delta$E (2.2) yields $\bm{r}^{n}(\bm{x}_{\text{ref}}+\bm{g}(\hat{\bm{x}}^{n}))=\bm{0}$. This is an overdetermined system of equations and thus may not have a solution. As such, we perform projection by enforcing orthgonality of this residual to a test basis $\bm{\Psi}_{G}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{N\times p}$ such that the discretize-then-project manifold Galerkin O$\Delta$E corresponds to $\hat{\bm{r}}_{G,\text{disc}}^{n}(\hat{\bm{x}}^{n})=\bm{0}$ with

[TABLE]

The discretize-then-project manifold Galerkin ROM is equivalent to the proposed manifold Galerkin ROM if and only if there exists a full-rank matrix $\bm{A}^{n}_{G}:\mathbb{R}^{p}\rightarrow\mathbb{R}_{\star}^{p\times p}$ such that

[TABLE]

Case 1. We now consider the general (non-asymptotic) case. Eq. (4.8) holds for any sequence of approximate solutions $\hat{\bm{x}}^{n}$, $n=1,\ldots,{N_{t}}$ if and only if each term in this expansion matches, i.e.,

[TABLE]

where we have used $\sum_{j=0}^{k}\alpha_{j}=0$. We first consider the first and third conditions of (4.9). Because the inverse of a nonlinear operator—if it exists—must also be nonlinear, a necessary condition for these two requirements is that the trial manifold corresponds is affine such that the decoder satisfies (4.1). Substituting Eq. (4.1) into the first and third of the conditions of (4.9) yields

[TABLE]

for $j=1,\ldots,k$. A necessary and sufficient conditions for these to hold is

[TABLE]

We now consider the second and fourth conditions of (4.9). Substituting Eqs. (4.1) and (4.11) into these conditions yields

[TABLE]

for $j=1,\ldots,k$. Noting that $\bm{\Phi}^{+}=[\bm{\Phi}^{{T}}\bm{\Phi}]^{-1}\bm{\Phi}^{{T}}$, a necessary and sufficient condition for (4.12) to hold is $\bm{\Phi}^{+}=[\bm{\Psi}_{G}(\hat{\bm{\xi}})^{{T}}\bm{\Phi}]^{-1}\bm{\Psi}_{G}(\hat{\bm{\xi}})^{{T}}.$ Because $\bm{\Phi}$ has no dependence on the generalized state, the test basis must also be independent of the generalized state, which yields

[TABLE]

Thus, without using any asymptotic arguments, the proposed manifold Galerkin ROM could have been derived using a discretize-then-project approach if and only if the test basis $\bm{\Psi}_{G}$ applied to the overdetermined FOM O$\Delta$E acting on the trial manifold satisfies Eq. (4.13). Examples of test bases that satisfy this condition include $\bm{\Psi}_{G}=[\bm{\Phi}^{+}]^{{T}}$ and $\bm{\Psi}_{G}=\gamma\bm{\Phi}$ for some nonzero scalar $\gamma\in\mathbb{R}$.

Case 2. We now consider asymptotic arguments. Assuming the nonlinear trial manifold is twice continuously differentiable, and $\|\hat{\bm{x}}^{n-j}-\hat{\bm{x}}^{n}\|=O(\Delta t)$, then for $\hat{\bm{\xi}}$ in a neighborhood of $\hat{\bm{x}}^{n}$ such that $\|\hat{\bm{\xi}}-\hat{\bm{x}}^{n}\|=O(\Delta t)$, we have

[TABLE]

Substituting these expressions into the definition of the discretize-then-project manifold Galerkin residual $\hat{\bm{r}}_{G,\text{disc}}^{n}$ and enforcing equivalence of each term in the matching conditions (4.8) and taking the limit $\Delta t\rightarrow 0$ yields

[TABLE]

where we have used $\sum_{j=0}^{k}\alpha_{j}=0$. A necessary and sufficient condition for the first and third conditions to hold is

[TABLE]

Then, a necessary and sufficient condition for the second and fourth conditions to hold is

[TABLE]

Examples of test bases that satisfy Eq. (4.17) include $\bm{\Psi}_{G}(\hat{\bm{\xi}})=[\bm{J}(\hat{\bm{x}}^{n})^{+}]^{{T}}$ and $\bm{\Psi}_{G}=\gamma\bm{J}(\hat{\bm{x}}^{n})$ for some nonzero scalar $\gamma\in\mathbb{R}$. ∎

Theorem 4.1 shows that manifold Galerkin can be derived using a discretize-then-project approach for nonlinear trial manifolds. This shows that previous analysis [12, Theorem 3.4] related to commutativity of (residual-minimizing) Galerkin projection and time discretization extends to the case of nonlinear trial manifolds.

We now show that the limiting equivalence result reported in [12, Section 5] between Galerkin and LSPG projection for linear subspaces also extends to nonlinear trial manifolds.

Theorem 4.2** (Limiting equivalence).**

Manifold Galerkin and manifold LSPG are equivalent if either (1) the trial manifold corresponds to an affine subspace and either an explicit scheme is employed or the limit $\Delta t\rightarrow 0$ is taken, or (2) the nonlinear manifold is twice continuously differentiable; $\|\hat{\bm{x}}^{n-j}-\hat{\bm{x}}^{n}\|=O(\Delta t)$ for all n\in{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\{1,\ldots,{N_{t}}\}}}, j\in{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\{1,\ldots,k\}}}; and the limit $\Delta t\rightarrow 0$ is taken.

Proof.

The proof follows similar steps to those applied in the proof of Theorem 4.1. The discrete residual $\hat{\bm{r}}_{L}^{n}$ defined in (3.27), which characterizes the manifold LSPG O$\Delta$E (3.18), can be written as

[TABLE]

where we have used $\sum_{i=0}^{k}\alpha_{i}=0$. Manifold LSPG projection is equivalent to manifold Galerkin projection if and only if the O$\Delta$Es characterizing the two methods are equivalent, which is equivalent to requiring the existence of a full-rank matrix $\bm{A}^{n}:\mathbb{R}^{p}\rightarrow\mathbb{R}_{\star}^{p\times p}$ such that

[TABLE]

Case 1. We now consider the general (non-asymptotic) case. Eq. (4.19) holds for any sequence of approximate solutions $\hat{\bm{x}}^{n}$, $n=1,\ldots,{N_{t}}$ if and only if each term in this expansion matches, i.e.,

[TABLE]

for $j=1,\ldots,k$. Comparing Conditions (4.20) with (4.9), we see that the requirements are the same as in Theorem 4.1, but with $\bm{\Psi}^{n}(\hat{\bm{\xi}};\bm{\nu})$ replacing $\bm{\Psi}_{G}$ and $\bm{A}^{n}$ replacing $\bm{A}^{n}_{G}$. Thus, we arrive at the same result: necessary and sufficient conditions for equivalence are that the trial manifold corresponds to an affine subspace such that the decoder satisfies (4.1) and the test basis $\bm{\Psi}^{n}(\hat{\bm{\xi}})$ satisfies

[TABLE]

Note that necessary conditions for Eq. (4.21) to hold are $\text{Ran}(\bm{\Psi}^{n}(\hat{\bm{\xi}}))=\text{Ran}(\bm{\Phi})$. Examples of test basis that satisfy this condition include $\bm{\Psi}^{n}(\hat{\bm{\xi}})=[\bm{\Phi}^{+}]^{{T}}$ and $\bm{\Psi}^{n}(\hat{\bm{\xi}})=\gamma\bm{\Phi}$ for some nonzero scalar $\gamma\in\mathbb{R}$. However, unlike in Theorem 4.1, we are not free to choose the test basis; rather, it is defined as

[TABLE]

For the necessary condition $\text{Ran}(\bm{\Psi}^{n}(\hat{\bm{\xi}}))=\text{Ran}(\bm{\Phi})$ to hold, the second term in Eq. (4.22) must be zero. This occurs if and only if either (1) an explicit scheme is employed such that $\beta_{0}=0$, or (2) the limit $\Delta t\rightarrow 0$ is taken. In both cases, $\bm{\Psi}^{n}(\hat{\bm{\xi}})=\gamma\bm{\Phi}$ holds with $\gamma=\alpha_{0}$.

Case 2. We now consider asymptotic arguments. Assuming the nonlinear trial manifold is twice continuously differentiable, and $\|\hat{\bm{x}}^{n-j}-\hat{\bm{x}}^{n}\|=O(\Delta t)$, then for $\hat{\bm{\xi}}$ in a neighborhood of $\hat{\bm{x}}^{n}$ such that $\|\hat{\bm{\xi}}-\hat{\bm{x}}^{n}\|=O(\Delta t)$, we obtain the expressions (4.14). Substituting these expressions into the definition of the $\hat{\bm{r}}_{L}^{n}$ in Eq. (4.18), enforcing equivalence of each term in the matching conditions (4.19), and taking the limit $\Delta t\rightarrow 0$ yields

[TABLE]

where we have used $\sum_{j=0}^{k}\alpha_{j}=0$. A necessary and sufficient condition for the first and third conditions to hold is

[TABLE]

Then, a necessary and sufficient condition for the second and fourth conditions to hold is

[TABLE]

As in Case 1 above, a necessary condition for (4.25) to hold is that $\text{Ran}(\bm{\Psi}^{n}(\hat{\bm{\xi}}))=\text{Ran}(\bm{J}(\hat{\bm{x}}^{n}))$. Due to the definition of the test basis, this requires the second term on the right-hand side of Eq. (3.19) to be zero. This is already satisfied by the stated assumption that the limit $\Delta t\rightarrow 0$ is taken, in which case $\bm{\Psi}^{n}(\hat{\bm{\xi}})=\gamma\bm{J}(\hat{\bm{\xi}})$ holds with $\gamma=\alpha_{0}$. It can be easily verified that this test basis satisfies the condition (4.25). ∎

We now derive a posteriori local discrete-time error bounds for both the manifold Galerkin and manifold LSPG projection methods in the context of linear multistep methods, and demonstrate that manifold LSPG projection sequentially minimizes the error bound in time. For notational simplicity, we suppress the second argument of the velocity $\bm{f}$, as the time superscript appearing in the second argument always matches that of the first argument.

Theorem 4.3** (Error bound).**

If the velocity $\bm{f}$ is Lipschitz continuous, i.e., there exists a constant $\kappa>0$ such that $\|\bm{f}(\bm{x})-\bm{f}(\bm{y})\|_{2}\leq\kappa\|\bm{x}-\bm{y}\|_{2}$ for all $\bm{x},\,\bm{y}\in\mathbb{R}^{N}$, and the time step is sufficiently small such that $\Delta t<{\left|\alpha_{0}\right|}/{\left|\beta_{0}\right|\kappa}$, then

[TABLE]

where $\hat{\bm{x}}_{G}^{n}$ denotes the manifold Galerkin solution satisfying (3.9) and $\hat{\bm{x}}_{P}^{n}$ denotes the manifold LSPG solution satisfying (3.17). Here, $h\vcentcolon=\left|\alpha_{0}\right|-\left|\beta_{0}\right|\kappa\Delta t$ and $\gamma_{j}\vcentcolon=\left(\left|\alpha_{j}\right|+\left|\beta_{j}\right|\kappa\Delta t\right)/h$.

Proof.

We begin by defining linear multistep residual operators associated with the FOM sequence of solutions $\{\bm{x}^{j}\}_{j=1}^{n}$ and an (arbitrary) sequence of approximate solutions $\{{\tilde{\bm{x}}}^{j}\}_{j=1}^{n}$, i.e.,

[TABLE]

We note that ${\tilde{\bm{x}}}^{j}=\bm{x}_{\text{ref}}+\bm{g}(\hat{\bm{x}}_{G}^{j})$, $j=1,\ldots,n$ in the case of manifold Galerkin projection, and ${\tilde{\bm{x}}}^{j}=\bm{x}_{\text{ref}}+\bm{g}(\hat{\bm{x}}_{P}^{j})$, $j=1,\ldots,n$ in the case of manifold LSPG projection. Subtracting $\tilde{\bm{r}}^{n}\left({\tilde{\bm{x}}}^{n}\right)$ from $\bm{r}_{\star}^{n}\left(\bm{x}^{n}\right)$ and noting from Eq. (2.2) that $\bm{r}_{\star}^{n}\left(\bm{x}^{n}\right)=\bm{0}$ yields

[TABLE]

Rearranging this expression gives

[TABLE]

We proceed by bounding $\|\text{(I)}\|_{2}$ from below, and $\|\text{(II)}\|_{2}$ from above. To bound $\|\text{(I)}\|_{2}$ from below, we apply the reverse triangle to obtain

[TABLE]

We now use Lipschitz continuity of the velocity $\bm{f}$ to bound $\|\text{(I.a)}\|_{2}$ from above as

[TABLE]

If time-step-restriction condition $\Delta t<{\left|\alpha_{0}\right|}/{\left|\beta_{0}\right|\kappa}$ holds, then we can combine inequalities (4.31) and (4.32) as

[TABLE]

To bound $\|\text{(II)}\|_{2}$ in Eq. (4.30) from above, we apply the triangle inequality and employ Lipschitz continuity of the velocity $\bm{f}$, which gives

[TABLE]

Combining Eq. (4.30) with inequalities (4.33) and (4.34) gives

[TABLE]

The manifold-Galerkin error bound (4.26) results from substituting ${\tilde{\bm{x}}}^{j}=\bm{x}_{\text{ref}}+\bm{g}(\hat{\bm{x}}_{G}^{j})$, $j=1,\ldots,n$ in (4.35), while the manifold-LSPG error bound (4.27) results from substituting ${\tilde{\bm{x}}}^{j}=\bm{x}_{\text{ref}}+\bm{g}(\hat{\bm{x}}_{P}^{j})$, $j=1,\ldots,n$ in (4.35) and noting that the manifold LSPG solution $\bm{x}_{\text{ref}}+\bm{g}(\hat{\bm{x}}_{P}^{j})$ satisifies the time-discrete residual-minimization property (3.17). ∎

This result illustrates that the time-discrete residual minimization property of manifold LSPG projection allows its approximation to sequentially minimize the error bound, as the first term on the right-hand side of bound (4.27) corresponds to the new contribution to the error bound at time instance $t^{n}$, while the second term on the right-hand side comprises the recursive term in the bound.

5 Nonlinear trial manifold based on deep convolutional autoencoders

This section describes the approach we propose for constructing the decoder $\bm{g}:\mathbb{R}^{p}\rightarrow\mathbb{R}^{N}$ that defines the nonlinear trial manifold. As described in the introduction, any nonlinear-manifold learning method equipped with a continuously differentiable mapping from the generalized coordinates (i.e., the latent state) to an approximation of the state (i.e., the data) is compatible with the manifold Galerkin and manifold LSPG methods proposed in Section 3. Here, we pursue deep convolutional autoencoders for this purpose, as they are very expressive and scalable; there is also high-performance software available for their construction. When the nonlinear trial manifold is constructed using the proposed deep convolutional autoencoder, we refer to the Manifold Galerkin ROM as the Deep Galerkin ROM and the Manifold LSPG ROM as the Deep LSPG ROM.111This is analogous to the naming convention employed for linear-subspace ROMs, as linear-subspace Galerkin projection with POD is referred to as a POD–Galerkin ROM.

Section 5.1 describes the mathematical structure of autoencoders, which provides a neural-network model for the decoder $\bm{g}(\cdot;\bm{\theta}):\mathbb{R}^{p}\rightarrow\mathbb{R}^{N}$, where $\bm{\theta}$ denotes the neural-network parameters to be computed during training; Section 5.2 describes the proposed autoencoder architecture, which is tailored to spatially distributed dynamical-system states; Section 5.3 describes how we propose to satisfy the initial condition using the proposed autoencoder in line with Remark 3.1.

5.1 Autoencoder

An autoencoder is a type of feedforward neural network that aims to learn the identity mapping, i.e., $\bm{h}:\bm{x}\mapsto\tilde{\bm{x}}$ with $\tilde{\bm{x}}\approx\bm{x}$ and $\bm{h}:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{x}}$. Autoencoders achieve this using an architecture composed of two parts: the encoder $\bm{h}_{\text{enc}}:\bm{x}\mapsto\hat{\bm{x}}$ with $\bm{h}_{\text{enc}}:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}^{n_{\hat{x}}}$ and $n_{\hat{x}}\ll n_{x}$, which maps a high-dimensional vector $\bm{x}$ to a low-dimensional code $\hat{\bm{x}}$; and the decoder $\bm{h}_{\text{dec}}:\hat{\bm{x}}\mapsto\tilde{\bm{x}}$ with $\bm{h}_{\text{dec}}:\mathbb{R}^{n_{\hat{x}}}\rightarrow\mathbb{R}^{n_{x}}$, which maps the code $\hat{\bm{x}}$ to an approximation of the original high-dimensional vector $\tilde{\bm{x}}$. Thus, the resulting autoencoder takes the form

[TABLE]

If $\bm{h}(\bm{x})\approx\bm{x}$ over a data set $\bm{x}\in\{\bm{x}^{(i)}\}_{i}$, then the low-dimensional codes $\bm{h}_{\text{enc}}(\bm{x}^{(i)})$ contain sufficient information to recover accurate approximations of the data $\tilde{\bm{x}}^{(i)}\approx\bm{x}^{(i)}$ via the application of the decoder $\bm{h}_{\text{dec}}$. In this work, we propose to employ the autoencoder decoder $\bm{h}_{\text{dec}}$ for the decoder $\bm{g}$ used to define the nonlinear trial manifold introduced in Section 3.1.

In feedforward networks, each network layer typically corresponds to a vector or tensor, whose values are computed by applying an affine transformation to the previous layer followed by a nonlinear activation function. An encoder with $n_{\text{L}}$ layers takes the form

[TABLE]

where $\bm{h}_{i}(\cdot;\bm{\Theta}_{i}):\mathbb{R}^{p_{i-1}}\rightarrow\mathbb{R}^{p_{i}}$, $i=1,\ldots,n_{\text{L}}$ denotes the function applied at layer $i$ of the neural network; $\bm{\Theta}_{i}$, $i=1,\ldots,n_{\text{L}}$ denotes the weights and the biases employed at layer $i$ with $\bm{\theta}_{\text{enc}}\equiv(\bm{\Theta}_{1},\ldots,\bm{\Theta}_{n_{\text{L}}})$; and $p_{i}$, $i=1,\ldots,n_{\text{L}}$ denotes the dimensionality of the output at layer $i$. The input has dimension $p_{0}=n_{x}$ and the final layer produces a code with dimension $p_{n_{\text{L}}}=n_{\hat{x}}$. The nonlinear activation function is applied in layers 1 to $n_{\text{L}}$ to a function of the weights, biases, and the outputs from the previous layer such that

[TABLE]

where $\bm{\phi}_{i}$ is an element-wise nonlinear activation function. For fully-connected layers as in the traditional multilayer perceptron (MLP), $h_{i}(\bm{\Theta}_{i},\bm{x})=\bm{\Theta}_{i}[1,\bm{x}^{{T}}]^{{T}}$ with $\bm{\Theta}_{i}\in\mathbb{R}^{p_{i}\times(p_{i-1}+1)}$ a real-valued matrix. For convolutional layers, $h_{i}$ corresponds to a convolution operator with $\bm{\Theta}_{i}$ providing the convolutional-filter weights.

A decoder with $\bar{n}_{\text{L}}$ layers also corresponds to a feedforward network; it takes the form

[TABLE]

with $\bar{\bm{h}}_{i}(\cdot;\bar{\bm{\Theta}}_{i}):\mathbb{R}^{\bar{p}_{i-1}}\rightarrow\mathbb{R}^{\bar{p}_{i}}$, $i=1,\ldots,\bar{n}_{\text{L}}$, and $\bm{\theta}_{\text{dec}}\equiv(\bar{\bm{\Theta}}_{1},\ldots,\bar{\bm{\Theta}}_{\bar{n}_{\text{L}}})$. The dimension of the input is $\bar{p}_{0}=n_{\hat{x}}$ and the dimension of the final layer (i.e., the output layer) is $\bar{p}_{\bar{n}_{\text{L}}}=n_{x}$. Again, the nonlinear activation is applied to the output of the previous layer such that

[TABLE]

As in the encoder, $\bar{h}_{i}(\bar{\bm{\Theta}}_{i},\bm{x})=\bar{\bm{\Theta}}_{i}[1,\bm{x}^{{T}}]^{{T}}$ with $\bar{\bm{\Theta}}_{i}\in\mathbb{R}^{\bar{p}_{i}\times(\bar{p}_{i-1}+1)}$ a real-valued matrix for fully-connected layers, while $\bar{h}_{i}$ corresponds to a convolution operator with $\bar{h}_{i}$ providing transposed convolutional-filter weights for convolutional layers.

Because MLP autoencoders are fully connected, the number of parameters (corresponding to the edge weights and biases) can be extremely large when the number of inputs $n_{x}$ is large; in such scenarios, these models typically require a large amount of training data. As this work aims to enable model reduction for large-scale dynamical systems, directly applying an MLP autoencoder to the state such that $n_{x}=N$ is not practical in many scenarios. To address this, alternative neural-network architectures have been devised that make use of parameter sharing to reduce the total number of parameters in the model, and thus the amount of data needed for training. In the context of dynamical systems characterized by a state that can be represented as spatially distributed data, we propose to apply convolutional autoencoders. Such methods are applicable to multi-channel spatially distributed input data and employ parameter sharing such that they can be trained with less data. Further, such models tend to generalize well to unseen test data [46, 47] because they exploit three key properties of natural signals: local connectivity, parameter sharing, and equivariance to translation [30, 46].

5.2 Deep convolutional autoencoders for spatially distributed

states

Many dynamical systems are characterized by a state that can be represented as spatially distributed data, e.g., spatially discretized partial-differential-equations models. In such cases, there is a restriction operator $\bm{R}$ that maps the state to a tensor representing spatially distributed data, i.e.,

[TABLE]

where $n_{i}$ denotes the number of discrete points in spatial dimension $i$; $d\in\{1,2,3\}$ denotes the spatial dimension, and $n_{\text{chan}}$ denotes the number of channels. For images, typically $n_{\text{chan}}=3$ (i.e., red, green, and blue). For dynamical system models the number of channels $n_{\text{chan}}$ is equal to the number of state variables defined at a given spatial location; for example, $n_{\text{chan}}$ is equal to the number of conserved variables when the dynamical system arises from the spatial discretization of a conservation law. We also write the associated prolongation operator, which aims to provide the inverse mapping such that

[TABLE]

For coarse discretizations of the spatial domain, it is possible to ensure $\bm{R}\circ\bm{P}(\cdot)$ corresponds to the identity map; for fine discretizations, it is possible to ensure $\bm{P}\circ\bm{R}(\cdot)$ corresponds to the identity map; for cases where underlying grid provides an isomorphic representation of the state, it is possible to achieve both [17].

After reformatting the state from a vector to a tensor by applying the restriction operator $\bm{R}$, we apply an invertible affine scaling operator $\bm{S}:\mathbb{R}^{n_{1}\times\cdots\times n_{d}\times n_{\text{chan}}}\rightarrow\mathbb{R}^{n_{1}\times\cdots\times n_{d}\times n_{\text{chan}}}$ for data-standardization purposes. Section 6.2 defines the specific elements of this operator, which are computed from the training data.

Figure 1 depicts the architecture of the proposed deep convolutional autoencoder. Before the first convolutional layer, the autoencoder applies the restriction operator $\bm{R}$ and scaling operator $\bm{S}$, while after the last convolutional layer, the autoencoder applies the inverse scaling operator $\bm{S}^{-1}$ and prolongation operator $\bm{P}$. We note that the restriction and prolongation operators are (deterministic) operations that simply reshape the state into an appropriate format for the autoencoder, while the scaling operator is defined explicitly from the range of the training data; thus, these quantities are not subject to training. We consider the combination of convolutional layers (gray boxes) and fully-connected layers (blue rectangles). B provides a description of convolutional layers, including hyperparameters and parameters that are subject to optimization. The encoder network is composed of restriction and scaling, followed by a sequence of $n_{\text{conv}}$ convolutional layers and $n_{\text{full}}$ fully-connected layers. The decoder network is composed of a sequence of $n_{\text{full}}$ fully-connected layers and $n_{\text{conv}}$ transposed-convolutional layers with no nonlinear activation in the last layer, followed by inverse scaling and prolongation. As a result, the dimension of the input to the first convolutional layer is $p_{0}=n_{\text{chan}}\prod_{i=1}^{d}n_{i}$. The dimension $p_{n_{\text{conv}}}$ of the output of encoder layer $n_{\text{conv}}$ (i.e., the final convolutional layer) is determined by the hyperparameters defining the kernels used in the convolutional layers (i.e., depth, stride, zero-padding). Similarly, the dimension $\bar{p}_{n_{\text{full}}}$ of the output of decoder layer $n_{\text{full}}$ (i.e., the final fully connected layer) is determined by the hyperparameters defining the kernels in the subsequent convolutional layers. The dimension of the output of the final decoder layer is $\bar{p}_{n_{\text{L}}}=p_{0}=n_{\text{chan}}\prod_{i=1}^{d}n_{i}$.

A particular instance of the network architecture can be defined by specifying the number of convolutional layers $n_{\text{conv}}$, the number of fully-connected layers $n_{\text{full}}$, the number of units in each layer, and the types of nonlinear activations. Such architecture-related parameters are typically considered hyperparameters, as they are not subject to optimization during training.

We propose to set the decoder for the proposed manifold ROMs to the decoder arising from the deep convolutional autoencoder architecture described in Figure 1, i.e., $\bm{g}=\bm{h}_{\text{dec}}$. When the Manifold Galerkin and Manifold LSPG ROMs employ this choice of decoder, we refer to them as Deep Galerkin and Deep LSPG ROMs, respectively.

5.3 Initial condition satisfaction

As discussed in Remark 3.1, the initial generalized coordinates $\hat{\bm{x}}^{0}(\bm{{\mu}})$ can be ensured to satisfy the initial conditions by employing a reference state defined by Eq. (3.4); however, this implies that the decoder must be able to accurately represent deviations from this reference state. To accomplish this for the proposed autoencoder, we propose (1) to train the autoencoder with snapshot data centered on the initial condition along with the zero vector (as described in Section 6.1), and (2) to set the initial generalized coordinates to an encoding of zero, i.e., $\hat{\bm{x}}^{0}(\bm{{\mu}})=\bm{h}_{\text{enc}}(\bm{0})$ for all $\bm{{\mu}}\in\mathcal{D}$. Then, setting the reference state according to Eq. (3.4) leads to a reference state of $\bm{x}_{\text{ref}}(\bm{{\mu}})=\bm{x}^{0}(\bm{{\mu}})-\bm{g}(\bm{h}_{\text{enc}}(\bm{0}))$. Further, if $\bm{h}_{\text{dec}}(\bm{h}_{\text{enc}}(\bm{0}))\approx\bm{0}$, as is encouraged by including the zero vector in training, then the reference state comprises a perturbation of the initial condition, and the decoder $\bm{h}_{\text{dec}}$ must only represent deviations from this perturbation. This is consistent with the training of the autoencoder, as the snapshots have been centered on the initial condition.

6 Offline training

This section describes the offline training process used to train the deep convolutional autoencoder proposed in Section 5.2. The approach employs precisely the same snapshot data used by POD. Section 6.1 describes the (snapshot-based) data collection procedure, which is identical to that employed by POD. Section 6.2 describes how the data are scaled to improve numerical stability of autoencoder training. Section 6.3 summarizes the gradient-based-optimization approach employed for training the autoencoder (i.e., computing the optimal parameters $\bm{\theta}^{\star}$) given the scaled training data. Algorithm 1 provides a summary of offline training.

6.1 Snapshot-based data collection

The first step of offline training is snapshot-based data collection. This requires solving the FOM O$\Delta$E (2.2) for training-parameter instances $\bm{{\mu}}\in\mathcal{D}_{\text{train}}\equiv\{\bm{{\mu}}_{\text{train}}^{i}\}_{i=1}^{{n_{\text{train}}}}\subset\mathcal{D}$ and assembling the snapshot matrix

[TABLE]

with $n_{\text{snap}}\vcentcolon={N_{t}}{n_{\text{train}}}$ and $\bm{W}(\bm{{\mu}})\vcentcolon=[\bm{x}^{1}(\bm{{\mu}})-\bm{x}^{0}(\bm{{\mu}})\ \cdots\ \bm{x}^{{N_{t}}}(\bm{{\mu}})-\bm{x}^{0}(\bm{{\mu}})].$

Remark 6.1** (Proper orthogonal decomposition).**

POD employs the snapshot matrix $\bm{W}$ to compute a trial basis matrix $\bm{\Phi}$ used to define an affine trial subspace $\bm{x}_{\text{ref}}(\bm{{\mu}})+\text{Ran}(\bm{\Phi})$. To do so, POD computes the singular value decomposition (SVD) and sets the trial basis matrix to be equal to the first $p$ left singular vectors, i.e.,

[TABLE]

The resulting POD trial basis matrix $\bm{\Phi}$ satisfies the minimization problem

[TABLE]

where $\bm{w}^{i}$ denotes the $i$th column of $\bm{W}$ and $V_{k}(\mathbb{R}^{n})$ and denotes the set of orthogonal $k\times n$ matrices (the compact Stiefel manifold). The solution is unique up to a rotation. This technique is equivalent (up to the data-centering process) to principal component analysis [38], an unsupervised machine-learning method for linear dimensionality reduction.

6.2 Data standardization

As described in Section 5.2, the first layer of the proposed autoencoder applies a restriction operator $\bm{R}$, which reformats the input vector into a tensor compatible with convolutional layers, followed by an affine scaling operator $\bm{S}$; the last layer applies the inverse of this scaling operator $\bm{S}^{-1}$ and subsequently applies the prolongation operator $\bm{P}$ to reformat the data into a vector. We now define the scaling operator from the training data to ensure that all elements of the training data lie between zero and one. This scaling improves numerical stability of the gradient-based optimization for training [39, 69], which will be described in Section 6.3. We adopt a standard scaling procedure also followed, e.g., by Ref. [29]. Namely, defining the restriction of the $i$th snapshot as ${\mathcal{W}}^{i}\vcentcolon=\bm{R}(\bm{w}^{i})\in\mathbb{R}^{n_{1}\times\cdots\times n_{d}\times n_{\text{chan}}}$, we set the elements of the scaling operator $\bm{S}$ to

[TABLE]

with

[TABLE]

6.3 Autoencoder training

Once the autoencoder architecture has been defined (including the restriction, prolongation and scaling operators), it takes the form $\bm{h}(\bm{x};\bm{\theta})=\bm{h}_{\text{dec}}(;\bm{\theta}_{\text{dec}})\circ\bm{h}_{\text{enc}}(\bm{x};\bm{\theta}_{\text{enc}})$, where the undefined parameters $\bm{\theta}$ correspond to convolutional-filter weights and weights and biases for the fully connected layers. We compute optimal values of these parameters $\bm{\theta}^{\star}$ using a standard approach from deep learning: stochastic gradient descent (SGD) with minibatching and early stopping [8]. A provides additional details, where Algorithm 2 provides the training algorithm.

7 Numerical experiments

This section assesses the performance of the proposed Deep Galerkin and Deep LSPG ROMs, which employ nonlinear trial manifolds, compared to POD–Galerkin and POD–LSPG ROMs, which employ affine trial subspaces. We consider two advection-dominated benchmark problems: 1D Burgers’ equation and a chemically reacting flow. We employ the numerical PDE tools and ROM functionality provided by pyMORTestbed [79], and we construct the autoencoder using TensorFlow [1].

For both benchmark problems, the Deep Galerkin and Deep LSPG ROMs employ a 10-layer convolutional autoencoder corresponding to the architecture depicted in Figure 1. The encoder $\bm{h}_{\text{enc}}$ consists of n_{\text{L}}={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}5}} layers with $n_{\text{conv}}=4$ convolutional layers, followed by n_{\text{full}}={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1}} fully-connected layer. The decoder $\bm{h}_{\text{dec}}$ consists of n_{\text{full}}={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1}} fully-connected layer, followed by $n_{\text{conv}}=4$ transposed-convolution layers. The latent code of the autoencoder is of dimension $p$ (i.e., $p_{n_{\text{L}}}=\bar{p}_{0}=p$), which will vary during the experiments to define different reduced-state dimensions. Table 1 specifies attributes of the kernels used in convolutional and transposed-convolutional layers.

For the nonlinear activation functions $\bm{\phi}_{i}$, $i=1,\ldots,n_{\text{L}}$ and $\bar{\bm{\phi}}_{i}$, $i=1,\ldots,\bar{n}_{\text{L}}-1$, we use exponential linear units (ELU) [19], which is defined as

[TABLE]

and an identity activation function in the output at layer $\bar{n}_{\text{L}}$ in decoder (as is common practice). This choice of activation ensures that the resulting decoder $\bm{g}$ is continuously differentiable everywhere, and is twice continuously differentiable almost everywhere. We employ the $\ell^{2}$-loss function defined in Eq. (A.2) in the minimization problem (A.1), which is equivalent to the loss function minimized by POD (see Remark 6.1). We apply the Adam optimizer [42], which is compatible with the stochastic gradient descent (SGD) algorithm reported in Algorithm 2; here, the adaptive learning rate strategy computes rates for different parameters using estimates of first and second moments of the gradients.

Using the same snapshots as that to train the autoencoder, we also compute a POD basis $\bm{\Phi}$ following the steps discussed in Remark 6.1.

We compare the performance of four ROMs: 1) POD–Galerkin: linear-subspace Galerkin projection (Eq. (4.4)) with the POD basis defining $\bm{\Phi}$, 2) POD-LSPG: linear-subspace LSPG projection (Eq. (4.6)) with the POD basis defining $\bm{\Phi}$, 3) Deep Galerkin: manifold Galerkin projection (Eq. (3.8)) with the deep convolutional decoder defining $\bm{g}$, and 4) Deep LSPG: manifold LSPG projection (Eq. (3.18)) with the deep convolutional decoder defining $\bm{g}$. All ROMs enforce the initial condition by setting the reference state according to Remark 3.1; this implies $\bm{x}_{\text{ref}}(\bm{{\mu}})=\bm{x}^{0}(\bm{{\mu}})$ for the linear-subspace ROMs.

To solve the O$\Delta$Es arising at each time instance, we apply Newton’s method for POD–Galerkin, the Gauss–Newton method for POD–LSPG, and the quasi-Newton methods proposed in Section 3.4 for the Deep Galerkin and Deep LSPG. We terminate the (quasi)-Newton iterations when the residual norm drops below $10^{-6}$ of its initial guess at that time instance; the initial guess corresponds to the solution at the previous time instance.

To assess the ROM accuracy, we compute the relative $\ell^{2}$-norm of the state error

[TABLE]

We also include the projection error of the solution

[TABLE]

onto both (1) the POD basis, in which case $\bm{\Phi}$ is the POD basis employed by POD–Galerkin and POD–LSPG, and (2) the optimal basis, in which case $\bm{\Phi}=\bm{\Phi}^{\star}$, which consists of the first $p$ left singular vectors of the snapshot matrix collected at the online-parameter instance $\bm{W}(\bm{{\mu}})$. We refer to the former metric as the POD projection error, as it provides a lower bound for the POD–Galerkin and POD–LSPG relative errors; we refer to the latter metric as the optimal projection error, as it provides an $\ell^{2}$-norm counterpart to the Kolmogorov $p$-width. Finally, we compute the projection of the FOM solution on the trial manifold

[TABLE]

and compute its relative error by employing $\tilde{\bm{x}}^{n}\leftarrow\tilde{\bm{x}}_{\star}^{n}$ in Eq. (7.1); we refer to this as the manifold projection error.

7.1 1D Burgers’ equation

We first consider a parameterized inviscid Burgers’ equation [64], as it comprises a very simple benchmark problem for which linear subspaces are ill suited due to its slowly decaying Kolmogorov $n$-width. The governing system of partial differential equations (with initial and boundary conditions) is

[TABLE]

where the flux is $f(w)=0.5w^{2}$ and there are ${n_{\mu}}=2$ parameters; thus, the intrinsic solution-manifold dimension is $p^{\star}=3$ (see Remark 2.1). We set the parameter domain to $\mathcal{D}=[4.25,5.5]\times[0.015,0.03]$ and the final time to $T=35$. We apply Godunov’s scheme with 256 control volumes to spatially discretize Eq. (7.4), which results in a system of parameterized ODEs of the form (2.1) with $N=256$ spatial degrees of freedom and initial condition $\bm{x}^{0}(\bm{{\mu}})=\bm{x}^{0}=\bm{1}$. For time discretization, we use the backward-Euler scheme, which corresponds to a linear multistep scheme with $k=1$, $\alpha_{0}=\beta_{0}=1$, $\alpha_{1}=-1$, and $\beta_{1}=0$ in Eq. (2.3). We consider a uniform time step $\Delta t=0.07$, resulting in ${N_{t}}=500$ time instances.

For offline training, we set the training-parameter instances to \mathcal{D}_{\text{train}}=\{(4.25+(1.25/9)i,\ 0.015+(0.015/7)j{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0})}}\}_{i=0,\ldots,9;\ j=0,\ldots 7}, resulting in ${n_{\text{train}}}=80$ training-parameter instances. The restriction operator $\bm{R}$ and prolongation operator $\bm{P}$ correspond to reshaping operators (without interpolation); the restriction operator reshapes the state vector into a tensor corresponding to the finite-volume grid such that $d=1$, $n_{1}=256$, and $n_{\text{chan}}=1$ in definitions (5.1) and (5.2). Then, we apply Algorithm 1 with inputs specified above and the following SGD hyperparameters: the fraction of snapshots to use for validation $\omega=0.1$; Adam optimizer learning-rate strategy with an initial uniform learning rate $\eta=10^{-4}$; initial parameters $\bm{\theta}^{(0)}$ computed via Xavier initialization [28] for weights and zero for biases; the number of minibatches determined by a fixed batch size of ${m}^{\mathcal{I}(i)}=20$, $i=1,\ldots,n_{\text{batch}}$; a maximum number of epochs $n_{\text{epoch}}=1000$; and early-stopping enforced if the loss on the validation set fails to decrease over 100 epochs. For the online stage, we consider two online-parameter instances $\bm{{\mu}}_{\text{test}}^{1}=\left(4.3,0.021\right)$, and $\bm{{\mu}}_{\text{test}}^{2}=\left(5.15,0.0285\right)$, which are not included in $\mathcal{D}_{\text{train}}$.

Figure 2 reports solutions at four different time indices computed by using FOM O$\Delta$E and the four considered ROMs. All ROMs employ the same reduced dimension of $p=10$ in Figure 2 (left) and $p=20$ in Figure 2 (right). These results demonstrate that nonlinear-manifold ROMs Deep Galerkin and Deep LSPG produce extremely accurate solutions, while the linear-subspace ROMs—constructed using the same training data—exhibit significant errors. This is due to the fundamental ill-suitedness of linear trial subspaces to advection-dominated problems.

Figure 3 reports the convergence of the relative error as a function of reduced dimension $p$. These results illustrate the promise of employing nonlinear trial manifolds. First, we note that employing a linear trial subspace immediately introduces significant errors: the projection error onto the optimal basis of dimension $p=p^{\star}=3$ is over 10%; the optimal nonlinear trial manifold (corresponding to the solution manifold) of the same dimension yields zero error. Even with a reduced dimension $p=50$, the relative projection error onto the optimal basis has not yet reached 0.1%. Second, we note that with a reduced dimension of only $p=5=p^{\star}+2$, the Deep LSPG ROM realizes relative errors near 0.1%, while linear-subspace ROMs (including the projection error with the optimal basis) exhibit relative errors near 10%. Thus, the proposed convolutional autoencoder is very close to achieving the optimal performance of any nonlinear trial manifold; the dimension it requires to realize sub-0.1% errors is only two larger than the intrinsic solution-manifold dimension $p^{\star}=3$. We also observe that the POD–Galerkin and POD–LSPG ROMs are also nearly able to achieve optimal performance for linear-subspace ROMs, as their relative errors are close to the projection error onto the optimal basis; unfortunately, this error remains quite large relative to the Deep Galerkin and Deep LSPG ROMs. This highlights that the proposed nonlinear-manifold ROMs are able to overcome the fundamental limitations of linear-subspace ROMs on problems exhibiting a slowly decaying Kolmogorov $n$-width. We also observe that the POD basis is very close to the optimal basis, which implies that the training data are sufficient to accurately represent the online solution.

Additionally, we note that Deep LSPG outperforms Deep Galerkin for this problem, likely due to the fact that the residual-minimization problem is defined over a finite time step rather than time-instantaneously; similar results have been shown in the case of linear trial subspaces, e.g., in Refs. [12].

Finally, to illustrate the dependence of the proposed methods on the amount of training data, we vary the number of training-parameter instances in the set ${n_{\text{train}}}\in\{5,10,20,40,80,120\}$, and we use the resulting $n_{\text{snap}}={N_{t}}{n_{\text{train}}}$ snapshots to train both the autoencoder for the Deep Galerkin and Deep LSPG ROMs and to compute the POD basis for the POD–Galerkin and POD–LSPG ROMs. We fix the reduced dimension to $p=5$. For offline training, we define the maximum number of epochs and the early-stopping strategy in a manner that the Adam optimizer employs the same number of gradient computations to train each model. Figure 4 reports the relative error of the four ROMs for the diferent number of training-parameter instances. This figure demonstrates that the proposed Deep Galerkin and Deep LSPG ROMs yield accurate results—around 2% relative error for Deep Galerkin and less than 1% relative error for Deep LSPG—with only ${n_{\text{train}}}=5$ parameter instances. This highlights that—for this example—the proposed methods do not require an excessive amount of training data relative to standard POD–based ROMs. This lack of demand for a large amount of data is likely due to the use of convolutional layers in the autoencoder, which effectively employ parameter sharing to reduce significantly the total number of parameters in the autoencoder, and thus the amount of data needed for training.

7.2 Chemically reacting flow

We now consider a model of the reaction of a premixed $\text{H}_{2}$-air flame at constant uniform pressure [9]. The evolution of the flame is modeled by the nonlinear convection–diffusion–reaction equation

[TABLE]

where $\nabla$ denotes the gradient with respect to physical space, $\kappa$ denotes the molecular diffusivity, $\bm{v}$ denotes the velocity field, and

[TABLE]

denotes the thermo-chemical composition vector consisting of the temperature $w_{T}(\vec{x},t;\bm{{\mu}})$ and the mass fractions of chemical species $\text{H}_{2},\text{O}_{2}$, and $\text{H}_{2}\text{O}$, i.e., $w_{i}(\vec{x},t;\bm{{\mu}})$ for $i\in\{\text{H}_{2},\text{O}_{2},\text{H}_{2}\text{O}\}$. The nonlinear reaction source term $\bm{q}(\bm{w}(\vec{x},t;\bm{{\mu}});\bm{{\mu}})\equiv[q_{T}(\bm{w};\bm{{\mu}}),q_{\text{H}_{2}}(\bm{w};\bm{{\mu}}),q_{\text{O}_{2}}(\bm{w};\bm{{\mu}}),q_{\text{H}_{2}\text{O}}(\bm{w};\bm{{\mu}})]^{{T}}$ is of Arrhenius type and is defined as

[TABLE]

where $(\nu_{\text{H}_{2}},\nu_{\text{O}_{2}},\nu_{\text{H}_{2}\text{O}})=(2,1,-2)$ denote stoichiometric coefficients, $(W_{\text{H}_{2}},W_{\text{O}_{2}},W_{\text{H}_{2}\text{O}})=(2.016,31.9,18)$ denote molecular weights with units g$\cdot$mol*-1*, $\rho=1.39\times 10^{-3}$ g $\cdot$ cm*-3* denotes the density mixture, $R=8.314$ J $\cdot$ mol${}^{-1}\cdot$ K*-1* denotes the universal gas constant, and $Q=9800$K denotes the heat of the reaction. The ${n_{\mu}}=2$ parameters correspond to $\bm{{\mu}}=(A,E)$, which are the pre-exponential factor $A$ and the activation energy $E$; we set the corresponding parameter domain to $\mathcal{D}=[2.3375\times 10^{12},6.5\times 10^{12}]\times[5.625\times 10^{3},9\times 10^{3}]$. Thus, the intrinsic solution-manifold dimension is $p^{\star}=3$ (see Remark 2.1). We set the molecular diffusivity to $\kappa=2\;$cm${}^{2}\cdot$s*-1*, and the velocity field set to be constant and divergence-free with $\bm{v}=[50\;\text{cm}\cdot\text{s}^{-1},\;0]^{T}$. We set the final time to $T=0.06\;\text{s}$.

Figure 5 reports the geometry of the spatial domain. On the inflow boundary $\Gamma_{2}$, we impose Dirichlet boundary conditions $w_{\text{H}_{2}}=0.0282$, $w_{\text{O}_{2}}=0.2259$ and $w_{\text{H}_{2}\text{O}}=0$ for the chemical-species mass fractions and $w_{T}=950$K for the temperature. On boundaries $\Gamma_{1}$ and $\Gamma_{3}$, we impose homogeneous Dirichlet boundary conditions for the chemical-species mass fractions, and we set the temperature $w_{T}=300\;\text{K}$. On $\Gamma_{4},\Gamma_{5},$ and $\Gamma_{6}$, we impose homogeneous Neumann conditions on the temperature and mass fractions. We consider a uniform initial condition corresponding to a domain that is empty of chemical species such that $w_{\text{H}_{2}}=w_{\text{O}_{2}}=w_{\text{H}_{2}\text{O}}=0$ and we set the temperature to $w_{T}=300\;\text{K}$.

We employ a finite-difference method with 65 grid points in the horizontal direction and 32 grid points in the vertical direction to spatially discretize Eq. (7.5), which results in a system of parameterized ODEs of the form (2.1) with $N=8192$ degrees of freedom.

For time discretization, we employ the second order backward difference scheme (BDF2), which corresponds to a linear multistep scheme with $k=2$, $\alpha_{0}=1$, $\alpha_{1}=-\frac{4}{3}$, $\alpha_{2}=\frac{1}{3}$, $\beta_{0}=\frac{2}{3}$, and $\beta_{1}=\beta_{2}=0$ in Eq. (2.3). We consider a uniform time step $\Delta t=10^{-4}$, resulting in ${N_{t}}=600$ time instances.

For offline training, we set the training-parameter instances to \mathcal{D}_{\text{train}}=\{(2.3375\times 10^{12}+(3.2725\times 10^{12}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}/7)i}},\ 5.625\times 10^{3}+(3.375\times 10^{3}/7)j{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0})}}\}_{i=0,\ldots,7;\ j=0,\ldots 7}, resulting in ${n_{\text{train}}}=64$ training-parameter instances. The restriction operator $\bm{R}$ and prolongation operator $\bm{P}$ correspond to reshaping operators (without interpolation); the restriction operator reshapes the state vector into a tensor corresponding to the finite-difference grid such that $d=2$, $n_{1}=64$, $n_{2}=32$, and $n_{\text{chan}}=4$ in definitions (5.1) and (5.2). We emphasize that each of the 4 unknown variables is considered a different channel for the input data. Then, we apply Algorithm 1 with inputs specified above and the following SGD hyperparameters: the fraction of snapshots to use for validation $\omega=0.1$; Adam optimizer learning-rate strategy with an initial uniform learning rate $\eta=10^{-4}$; initial parameters parameters $\bm{\theta}^{(0)}$ computed via He initialization [34] for weights and zero for biases; number of minibatches determined by a fixed batch size of ${m}^{\mathcal{I}(i)}=20$, $i=1,\ldots,n_{\text{batch}}$; a maximum number of epochs $n_{\text{epoch}}=5000$; and early-stopping enforced if the loss on the validation set fails to decrease over 500 epochs. For the online stage, we consider parameter instances $\bm{{\mu}}_{\text{test}}^{1}=(2.5\times 10^{12},\,5.85\times 10^{3})$ and $\bm{{\mu}}_{\text{test}}^{2}=(3.2\times 10^{12},\,7.25\times 10^{3})$ that are not included in $\mathcal{D}_{\text{train}}$.

Figure 6 reports the FOM solution for online-parameter instance $\bm{{\mu}}_{\text{test}}^{1}$ and final time $t=T=0.06\;\text{s}$. Figures 7 and 8 report relative errors of the temperature solutions and the $\text{H}_{2}$ mass fraction solutions at the final time $t=T$ computed using all considered ROMs with the same reduced dimension $p=3$. Again, we observe that the proposed nonlinear-manifold ROMs produce significantly lower errors as compared with the linear-subspace ROMs. Figure 9 reports the convergence of the relative error as a function of reduced dimension $p$ along with the projection errors Eq. 7.2 onto the POD basis and the optimal basis. This figure again shows that the proposed manifold ROMs significantly outperform the linear-subspace ROMs, as well as the projection error onto both the POD and optimal bases. First, we observe that employing a linear trial subspace again introduces significant errors, as the projection error onto the optimal basis of dimension $p=p^{\star}=3$ is around 5%; the optimal nonlinear trial manifold (corresponding to the solution manifold) of the same dimension yields zero error. Second, we note that for a reduced dimension of only $p=p^{\star}=3$, both the Deep Galerkin and Deep LSPG ROMs yield relative errors of less than 0.1%, while the projection errors onto the POD and optimal bases exceed 5%, and the linear-subspace ROMs yield relative errors in excess of 40%. Thus, the proposed convolutional autoencoder nearly achieves optimal performance, as the dimension it requires to realize sub-0.1% errors is exactly equal to the intrinsic solution-manifold dimension $p^{\star}=3$.

We also observe that there is a significant gap between the performance of the linear-subspace ROMs and the POD projection error; this gap is attributable to the closure problem. We also observe that—in contrast with the previous example—there is a non-trivial gap between the POD projection error and the optimal-basis projection error, which suggests that the online solution is less well represented by the training data than in the previous case.

Additionally, we observe that for a reduced dimension of $p=10$, the projection error onto the optimal basis begins to become smaller than the relative error of the proposed Deep Galerkin and Deep LSPG ROMs; however this occurs for (already small) errors less than 0.1%. The error saturation of the Deep Galerkin and Deep LSPG ROMs is likely due to the fact that the online solution cannot be perfectly represented using the training data, as illustrated by the gap between the projection errors associated with POD and the optimal basis.

We vary the number of training-parameter instances in the set ${n_{\text{train}}}\in\{4,8,16,32,64\}$, and we use the resulting $n_{\text{snap}}={N_{t}}{n_{\text{train}}}$ snapshots to train both the autoencoder for the Deep Galerkin and Deep LSPG ROMs and to compute the POD basis for the POD–Galerkin and POD–LSPG ROMs. We fix the reduced dimension to $p=5$. Figure 10 reports the relative error of the four ROMs for these different amounts of training data. This figure demonstrates that the proposed Deep Galerkin and Deep LSPG ROMs yield accurate results—less than 1% relative error for both Deep Galerkin and Deep LSPG—with only ${n_{\text{train}}}=4$ parameter instances. These results again show that the proposed methods do not appear to require an excessive amount of training data relative to standard POD–based ROMs.

8 Conclusion

This work has proposed novel manifold Galerkin  and manifold LSPG projection techniques, which project dynamical-system models onto arbitrary continuously-differentiable nonlinear manifolds. We demonstrated how these methods can exactly satisfy the initial condition, and provided quasi-Newton solvers for implicit time integration.

We performed analyses that demonstrated that employing an affine trial manifold recovers classical linear-subspace Galerkin and LSPG projection. We also derived sufficient conditions for commutativity of time discretization and manifold Galerkin projection, as well as conditions under which manifold Galerkin and manifold LSPG projection are equivalent. In addition, we derived a posteriori time-discrete error bounds for the proposed methods.

We also proposed a practical strategy for computing a representative low-dimensional nonlinear trial manifold that employs a specific convolutional autoencoder tailored to spatially distributed dynamical-system states. When the Manifold Galerkin and Manifold LSPG ROMs employ this choice of decoder, we refer to them as Deep Galerkin and Deep LSPG ROMs, respectively.

Finally, numerical experiments demonstrated the ability of the method to produce substantially lower errors for low-dimensional reduced states than even the projection onto the optimal basis. Indeed, the proposed methodology is nearly able to achieve optimal performance for a nonlinear-manifold method; in the first case, the reduced dimension required to achieve sub-0.1% errors is only two larger than the intrinsic solution-manifold dimension; in the second case, the reduced dimension required to achieve such errors is exactly equal to the intrinsic solution-manifold dimension.

We note that one drawback of the method in its current form (relative to classical linear-subspace methods) is that it incurs a costlier training process, as training a deep convolutional autoencoder is significantly more computationally expensive than simply computing the left singular vectors of a snapshot matrix. In addition, the proposed method is characterized by significantly more hyperparameters—which relate to the autoencoder architecture—than classical methods.

Future work involves integrating the proposed manifold Galerkin and manifold LSPG projection methods within a full hyper-reduction framework to realize computational cost savings; this investigation will leverage sparse norms as discussed in Remarks 3.4 and 3.6, and will also consider specific instances of the proposed autoencoder architecture that enable computationally efficient hyper-reduction. Additional future work includes appending structure-preserving constraints to the minimum-residual formulations [14], implementing the proposed techniques in a production-level simulation code and demonstrating the methods on truly large-scale dynamical-system models.

Acknowledgments

The authors gratefully acknowledge Matthew Zahr for graciously providing the pyMORTestbed code that was modified to obtain the numerical results, as well as Jeremy Morton for useful discussions on the application of convolutional neural networks to simulation data. The authors also thank Patrick Blonigan and Eric Parish for providing useful feedback. This work was sponsored by Sandia’s Advanced Simulation and Computing (ASC) Verification and Validation (V&V) Project/Task #103723/05.30.02. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525.

Appendix A Stochastic gradient descent for autoencoder training

This section briefly summarizes stochastic gradient descent (SGD) with minibatching and early stopping, which is used to compute the parameters $\bm{\theta}^{\star}$ of the autoencoder.

We begin by randomly shuffling the snapshot matrix $\bm{W}$ (defined in Eq. (6.1)) into neural-network training snapshots $\bm{W}_{\text{train}}\in\mathbb{R}^{N\times m}$ and validation snapshots $\bm{W}_{\text{val}}\in\mathbb{R}^{N\times\bar{m}}$, i.e.,

[TABLE]

where $\bm{P}\in\{0,1\}^{n_{\text{snap}}\times n_{\text{snap}}}$ denotes a random permutation matrix and $m+\bar{m}=n_{\text{snap}}$, typically with $\bar{m}\approx 0.1m$.

Given the training data $\bm{W}_{\text{train}}$, we compute the optimal parameters $\bm{\theta}^{\star}$ by (approximately) solving the optimization problem

[TABLE]

which minimizes the empirical risk over the training data. Here, $\hat{p}_{\text{data}}$ denotes the empirical distribution associated with the training data $\bm{W}_{\text{train}}$ and the loss function $\mathcal{L}$ provides a measure of discrepancy between training snapshot $\bm{x}^{i}_{\text{train}}$ and its reconstruction $\bm{h}(\bm{x}^{i}_{\text{train}};\bm{\theta})$; for example, the $\ell^{2}$-loss function is often employed such that

[TABLE]

in which case the loss function is equivalent to that employed for POD (compare with Problem (6.3)). Note that autoencoder training is categorized as an unsupervised learning (or semi-supervised learning) problem, as there is no target response variable other than recovery of the original input data.

We apply SGD to (approximately) solve optimization problem (A.1), which leads to parameter updates at the $i$th iteration of the form

[TABLE]

Here, $\eta\in\mathbb{R}$ denotes the step length or learning rate and $\widetilde{\nabla\bm{J}}^{(i)}(\approx\nabla\bm{J})$ denotes a gradient approximation at optimization iteration $i$. The gradient approximation corresponds to the sample mean of the gradient over a minibatch $\bm{W}_{\text{mini},\mathcal{I}(i)}\in\mathbb{R}^{N\times{m}^{\mathcal{I}(i)}}$, where

[TABLE]

where $\mathcal{I}:i\mapsto((i-1)\mod n_{\text{batch}})+1$ provides the mapping from optimization iteration $i$ to batch index, and $\bm{P}_{\text{train}}\in\{0,1\}^{N\times m}$ denotes a random permutation matrix. That is, the gradient approximation at optimization iteration $i$ is

[TABLE]

For small batch sizes, this gradient approximation not only significantly reduces the per-iteration cost of the optimization algorithm, it can also improve generalization performance, and—for early iterations and for minibatch size 1—leads to the same sublinear rate of convergence of the expected risk as the empirical risk [8]. In practice, each gradient contribution $\nabla_{\bm{\theta}}\mathcal{L}(\bm{x}_{\text{mini},\mathcal{I}(i)}^{j},\bm{\theta})$ is computed from the chain rule via automatic differentiation, which is referred to in deep learning as backpropagation [68].

Some optimization methods employ adaptive learning rates that are tailored for each parameter, which comprises a modification of the parameter update (A.3) to

[TABLE]

where $\bm{\eta}^{(i)}$ denotes a vector of learning rates and $\odot$ denotes the (element-wise) Hadamard product. Examples include AdaGrad [25], RMSProp [74], and Adam.

Although we have presented the specific case of non-adaptive minibatches (which we employ in the numerical experiments), it is also possible to employ adaptive batch sizes; often, the batch size increases with iteration count to produce a lower-variance gradient estimate as a local minimum is approached [8].

Rather than terminating iterations when a local minimum of the objective function $\bm{J}$ is reached, we instead employ early stopping, which is a widely used form of regularization that has been shown to improve generalization performance in many applications [30, 8]. Here, we terminate iterations when the loss function on the validation snapshots

[TABLE]

does not decrease for a certain number of epochs, where an epoch is equivalent to $n_{\text{batch}}$ iterations, i.e., a single pass through the training data. Early stopping effectively treats the number of optimization iterations as a hyperparameter, as allowing a large number of iterations can be associated with a higher-capacity model.

Algorithm 2 describes autoencoder training using SGD with minibatching and early stopping. Note that the adaptive learning rate strategy, initial parameters $\bm{\theta}^{(0)}$, number of minibatches $n_{\text{batch}}$, maximum number of epochs $n_{\text{epoch}}$, and early-stopping criterion are all SGD hyperparameters that comprise inputs to the algorithm. Line 6 of Algorithm 2 applies the adaptive learning rate strategy. For example, AdaGrad scales the learning rate to be inversely proportional to the square root of the sum of all the historical squared values of the gradient [25]; Adam updates the learning rate based on estimates of the first moment (the mean) and the second moment (the uncentered variance) of the gradients [42]. Alternatively, the learning rate can be kept to a single constant for all parameters over the entire training procedure, which simplifies the parameter update procedure in Line (7) of Algorithm 2 to the typical SGD update (A.3). See Ref. [8] for a review of training deep neural networks.

Appendix B Basics of convolutional layers

This section provides notations and basic operations performed in convolutional layers. See [26] for more details of the convolution arithmetic for deep learning.

Units in convolutional layers are organized as feature maps ${\mathcal{H}}$, where each unit is connected to the local patches of the feature maps of the previous layer through a discrete convolution defined by a set of kernels, which is denoted by filter bank ${\mathcal{W}}$, followed by a nonlinear activation and a pooling opertion. The feature map at layer $l$ can be considered as a 3-dimensional tensor ${\mathcal{H}}^{l}\in\mathbb{R}^{n_{\text{chan}}^{l}\times n_{1}^{l}\times n_{2}^{l}}$ with element ${\mathcal{H}}^{l}_{i,j,k}$ representing a unit within channel $i$ at row $j$ and column $k$, and the filter banks at layer $l$ can be considered as a 4-dimensional tensor ${\mathcal{W}}^{l}\in\mathbb{R}^{n_{\text{filter}}^{l}\times n_{\text{chan}}^{l-1}\times k_{1}^{l}\times k_{2}^{l}}$ with element ${\mathcal{W}}^{l}_{i,j,m,n}$ connecting between a unit in channel $i$ of the output and a unit in channel $j$ of the input, with an offset of $m$ rows and $n$ columns between the output unit and the input unit. The number of filters in the filter bank is denoted by $n_{\text{filter}}$, and the kernel length is characterized by $k_{1}$ and $k_{2}$. Convolving a feature map ${\mathcal{H}}^{l-1}$ with a filter bank ${\mathcal{W}}^{l}$ can be written as

[TABLE]

for all valid $i$, where $\mathcal{B}^{l}$ is a tensor indicating bias. Here, $s$ denote the stride, which determines downsampling rate of each convolution; only every $s$ elements is sampled in each direction in the output. By having $s>1$, the dimension of the next feature map can be reduced by factor of $s$ in each direction. The filter banks ${\mathcal{W}}$ and the biases ${\mathcal{B}}$ are learnable parameters, whereas the kernel length [$k_{1}$, $k_{2}$] and the number of filters $n_{\text{filter}}$, and the stride $s$ are the hyperparameters. After the nonlinearity, a pooling function is typically applied to the output, extracting a statistical summary of the neighboring units of the output at certain locations.

Appendix C Additional numerical experiments: extrapolation and interpolation in

time

We now perform additional investigations that assess the ability of the proposed Deep Galerkin and Deep LSPG ROMs to both extrapolate (C.1) and interpolate (C.2) in time. We also assess the effect of early stopping on the performance of these methods (C.3). All numerical experiments in this section are performed on the 1D Burgers’ equation described in Section 7.1 with the same setup except when otherwise specified.

C.1 Time extrapolation

To assess time extrapolation, we collect snapshots associated with the first $n_{t}(\leq{N_{t}}=500)$ time instances for each training-parameter instance and use the resulting $n_{\text{snap}}=n_{t}{n_{\text{train}}}$ snapshots to train the autoencoder for the Deep Galerkin and Deep LSPG ROMs and to compute the POD basis for the POD–Galerkin and POD–LSPG ROMs. Figure 11 reports the relative error for each of these ROMs of dimension $p=5$ at the online-parameter instances $\bm{{\mu}}_{\text{test}}^{1}$ and $\bm{{\mu}}_{\text{test}}^{2}$ for the number of collected snapshots varying in the set $n_{t}\in\{200,300,400,500\}$. These results show that the proposed Deep Galerkin and Deep LSPG ROMs yield more accurate results than the POD–Galerkin and POD–LSPG ROMs for $n_{t}\geq 300$, with Deep LSPG yielding relative errors around 0.1% for $n_{t}\geq 400$. Thus, we conclude that—for this example—the proposed Deep Galerkin and Deep LSPG generally yield superior time-extrapolation results than the POD–Galerkin and POD–LSPG ROMs.

C.2 Time interpolation

To assess time interpolation, we now consider collecting $n_{t}$ snapshots that are equispaced in the time interval at each training-parameter instance, with the final snapshot collected at the final time instance $t^{N_{t}}$. As in the case of time extrapolation, we then use the resulting $n_{\text{snap}}=n_{t}{n_{\text{train}}}$ snapshots to train both the autoencoder for the Deep Galerkin and Deep LSPG ROMs and to compute the POD basis for the POD–Galerkin and POD–LSPG ROMs. Figure 12 reports the relative error for each of these ROMs of dimension $p=5$ at the online-parameter instances $\bm{{\mu}}_{\text{test}}^{1}$ and $\bm{{\mu}}_{\text{test}}^{2}$ for the number of collected snapshots varying in the set $n_{t}\in\{10,20,50,100,140,250,500\}$. First, these results show that the POD–Galerkin and POD–LSPG ROMs are almost entirely insensitive to the number of snapshots collected within the time interval for $n_{t}\geq 10$. Similarly, the proposed Deep Galerkin and Deep LSPG ROMs are insensitive to the number of snapshots for $n_{t}\geq 100$. For $n_{t}\in\{10,20,50\}$, the Deep Galerkin and Deep LSPG ROMs yield larger errors than for $n_{t}\geq 100$; however, these errors are still smaller than the errors produced by the POD–Galerkin and POD–LSPG ROMs. Thus, we conclude that—for this example—the proposed methods do not exhibit accuracy degradation when only 20% of the snapshots (collected equally spaced in time) are used to train the autoencoder.

C.3 Early-stopping study

We now assess the effect of early stopping—a common approach in training deep neural networks for improving generalization performance—on the performance of the proposed Deep Galerkin and Deep LSPG methods. To do so, we set the reduced dimension to $p=30$ and train the autoencoder both with and without the early stopping strategy. To train without early stopping, we set the maximum number of epochs to $n_{\text{epoch}}=1000$, save the network parameters at each epoch, and select the parameters (out of all $n_{\text{epoch}}=1000$ candidates) that yield the smallest value of the loss function $\mathcal{L}$ (defined in Eq. A.2) on the validation set. Table 2 reports the relative errors obtained by the proposed ROMs—as well as the optimal projection error onto the manifold—at the online-parameter instances. In this case, the early-stopping strategy terminated iterations after 506 epochs, while the strategy outlined above selected the network parameters arising at the 944th epoch. These results show that adopting the ‘no early stopping’ strategy can slightly improve performance over the early-stopping stratgy. We hypothesize this to be the case because—for this example—the training data are quite informative of the prediction task; thus, attempting to prevent overfitting essentially amounts to underfitting.

Appendix D Computational costs

Let us consider the $i$th layer of the autoencoder,

[TABLE]

Fully-connected layers: If the $i$th layer is a fully-connected layer, then $\bm{\Theta}_{i}\in\mathbb{R}^{p_{i}\times(p_{i-1}+1)}$ is a real-valued matrix and the operation $\bm{h}_{i}$ comprises a matrix–vector product $h_{i}(\bm{\Theta}_{i},\bm{x})=\bm{\Theta}_{i}[1,\bm{x}^{{T}}]^{{T}}$, followed by an element-wise nonlinear activation $\bm{\phi}_{i}$. Thus, evaluating $\bm{\phi}_{i}(h_{i}(\bm{\Theta}_{i},\bm{x}))$ incurs $2p_{i}(p_{i-1}+1)+c_{\bm{\phi}}p_{i}$ floating-point operations (flops), where $c_{\bm{\phi}}$ denotes the number of flops incurred by applying the nonlinear activation to a scalar.

Convolutional layers: If the $i$th layer is a (transposed-)convolutional layer, the operation $\bm{h}_{i}$ comprises discrete convolutions, followed by an element-wise nonlinear activation $\bm{\phi}_{i}$. The input and the output of the $i$th convolutional layer are 3-way tensors ${\mathcal{H}}^{i-1}\in\mathbb{R}^{n_{\text{chan}}^{i-1}\times n_{1}^{i-1}\times n_{2}^{i-1}}$, and ${\mathcal{H}}^{i}\in\mathbb{R}^{n_{\text{chan}}^{i}\times n_{1}^{i}\times n_{2}^{i}}$ (see B). Each element of the output tensor is computed by applying convolutional filters to a certain local region at every channel of the input tensor, where a single application of the discrete convolution with filter size $k_{1}^{i}\times k_{2}^{i}$ incurs $2k_{1}^{i}k_{2}^{i}n_{\text{chan}}^{i-1}$ flops. That is, computing $p_{i}(=n_{\text{chan}}^{i}n_{1}^{i}n_{2}^{i})$ elements with convolutional filters of size $k_{1}^{i}\times k_{2}^{i}$ requires $2p_{i}k_{1}^{i}k_{2}^{i}n_{\text{chan}}^{i-1}$ flops. Thus, evaluating $\bm{\phi}_{i}(h_{i}(\bm{\Theta}_{i},\bm{x}))$ incurs $2p_{i}k_{1}^{i}k_{2}^{i}n_{\text{chan}}^{i-1}+c_{\bm{\phi}}p_{i}$ flops and, again, $c_{\bm{\phi}}$ denotes the number of flops incurred by applying the nonlinear activation to a scalar.

Restrictor, prolongator, scaling and inverse scaling operator: In this study, the restrictor and the prolongator are restricted to only change the shapes of the snapshot matrix into/from a tensor representing spatially distributed data and, thus, they do not incur floating-point operations. The scaling and the inverse scaling operators are applied to each FOM solution snapshot and incur $2N$ flops.

Figure 13 reports the estimated computational costs of the decoder with the latent code of dimension $p=5$ for varying choices of the hyperparameters. On the left, the figure shows the computational costs required at each layer of the decoder for varying convolutional filter sizes $k=\{25,16,9,5,3\}$ with the number of the input channels and the number of convolutional filters specified as $\{64,32,16,8,1\}$, where the first element is the number of the input channels and the last four elements are the number of convolutional filters, and with strides $\{4,4,4,2\}$. On the right, the figure shows the total computational costs for varying convolutional filter sizes $k=\{25,16,9,5,3\}$ and different numbers of input channels and convolutional filters $\{\{64,32,16,8,1\},\{32,16,8,4,1\},\{16,8,4,2,1\},\{8,4,2,1,1\},\{4,2,1,1,1\},\{2,1,1,1,1\}\}$ and with strides $\{4,4,4,2\}$. The gray plane is the computational cost for the linear decoder (e.g., POD) with the reduced dimension $p=85$, which produces approximated solutions with the accuracy comparable to the accuracy of the manifold ROM solutions with the autoencoder configured as described in Table 1 with the reduced dimension $p=5$. The figure also reports the computational cost for the linear decoder (e.g., POD) with the latent dimension $p=5$ (the blue plane).