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

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 and proper orthogonal decomposition (POD) for parameterized stationary problems; balanced truncation , rational interpolation , and Craig–Bampton model reduction for linear time invariant (LTI) systems; and Galerkin projection , least-squares Petrov–Galerkin projection , and other Petrov–Galerkin projections with (balanced) POD for nonlinear dynamical systems.

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

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 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 . 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’ , applying a coordinate transformation to the trial basis , shifting the POD basis , transforming the physical domain of the snapshots , constructing the trial basis on a Lagrangian formulation of the governing equations , and using Lax pairs of the Schrödinger operator to construct a time-evolving trial basis . 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 , physical domain , or state space . Ref. 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 , locally linear embedding (LLE) , Hessian eigenmaps , spectral embedding , and t-SNE . Methods that do provide this required mapping include self-organizing maps , generative topographic mapping , kernel principal component analysis (PCA) , Gaussian process latent variable model , diffeomorphic dimensionality reduction , and autoencoders . 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 ) 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 . 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

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

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

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

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

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

a posteriori discrete-time error bounds for both the manifold Galerkin and manifold LSPG projection methods (Theorem 4.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).

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. 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. 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. , which constructs a nonlinear model of wall turbulence using an autoencoder; Ref. , which applies an autoencoder to compress the state, followed by a recurrent neural network (RNN) to learn the dynamics; Refs. , which apply autoencoders to learn approximate invariant subspaces of the Koopman operator; and Ref. , 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.

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)

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. 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

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.

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 . 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.

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.

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

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.

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 , rational interpolation , the reduced-basis method , Rayleigh–Ritz eigenvectors ; common choices when the velocity is nonlinear include POD and balanced POD .

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

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

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

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.

and manifold Galerkin O$\Delta$E residual

3 Manifold LSPG ROM: time-discrete residual minimization

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

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.

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

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.,

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

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.

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

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

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

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

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

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

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

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}}$.

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.

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

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

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

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

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

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.

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.,

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

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

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

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

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

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

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

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.

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.

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

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.,

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

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

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

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

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). ∎

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$.

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

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

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

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

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

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.

Nonlinear trial manifold based on deep convolutional autoencoders

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

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

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 because they exploit three key properties of natural signals: local connectivity, parameter sharing, and equivariance to translation .

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.,

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

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 .

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.

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.

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.

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

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}})].$

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.,

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

2 Data standardization

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 . A provides additional details, where Algorithm 2 provides the training algorithm.

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 , and we construct the autoencoder using TensorFlow .

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) , which is defined as

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.

We also include the projection error of the solution

We first consider a parameterized inviscid Burgers’ equation , 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

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 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. .

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.

2 Chemically reacting flow

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

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

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

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 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.

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 , 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.

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

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

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

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 . 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 .

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

where $\bm{\eta}^{(i)}$ denotes a vector of learning rates and $\odot$ denotes the (element-wise) Hadamard product. Examples include AdaGrad , RMSProp , 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 .

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 . Here, we terminate iterations when the loss function on the validation snapshots

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 ; Adam updates the learning rate based on estimates of the first moment (the mean) and the second moment (the uncentered variance) of the gradients . 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. 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 for more details of the convolution arithmetic for deep learning.

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.

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,

References