Space-time least-squares Petrov-Galerkin projection for nonlinear model reduction

Introduction

Reduced-order models (ROMs) of nonlinear dynamical systems are essential for enabling high-fidelity computational models to be used in many-query and real-time applications such as uncertainty quantification, design optimization, and control. Such ROMs reduce the spatial dimensionality of the dynamical system by performing a projection process on the governing system of nonlinear ordinary differential equations (ODEs). The resulting ROM is then resolved in time via numerical integration, typically with the same time integrator and time step employed for the high-fidelity model. Unfortunately, many applications require simulating the model over long time intervals, leading to high temporal dimensionality characterized by the number of time instances in the time discretization. For example, many applications in fluid dynamics require long-time simulations to compute adequate statistics such as power spectral densities; structural-dynamics applications can demand long-time integration when structures undergo significant deformations; long-time integration is necessary to assess stability of planetary orbits ; molecular dynamics simulations and condensed phase dynamics also require long-time integration.

As such, ROMs are often characterized by low spatial dimensionality, but high temporal dimensionality, which can limit realizable computational savings in practice. It also renders ROMs ineffective in applications that demand a low temporal dimension for computational tractability. For example, a high temporal dimension can render intrusive uncertainty quantification methods (e.g., stochastic Galerkin ) and simultaneous analysis and design (SAND) in PDE-constrained optimization computationally intractable, as the dimension of the system of equations arising in such applications scales with the temporal dimension of the problem. Further, rigorous error bounds for these ROMs typically grow exponentially in time , which renders certification challenging. This work aims to devise a model-reduction methodology that enables significant reduction in both the spatial and temporal dimensions of the dynamical system, while simultaneously producing error bounds that exhibit slower growth in time.

Several attempts have been made to address this temporal-complexity bottleneck in model reduction. First, several authors have demonstrated that larger stable time steps (and thus a smaller number of time instances) can be taken with a ROM relative to the high-fidelity model in the case of explicit time integration . However, this approach is not always feasible, as many problems (e.g., compressible fluid dynamics, chemical kinetics) exhibit stiff dynamics that require implicit time integration, where the time step is limited by accuracy rather than stability. Increasing the time step in these contexts could significantly degrade time-discretization accuracy.

Time-parallel methods (e.g., parareal , PITA , and MGRIT ) aim to reduce the (serial) wall time incurred by a fine temporal discretization. These approaches enable dynamical-system simulations to be parallelized in the temporal domain, and are well suited for reduced-order models, as spatial parallelism alone quickly saturates for such low-dimensional models. However, while time-parallel methods can reduce the wall-time of such simulations, they do not reduce temporal dimensionality; in fact, they increase the total computational cost of simulations (as measured in core–hours).

More recently, a ‘forecasting’ approach was proposed that employs time-domain data (i.e., snapshot-matrix right singular vectors) generated during the offline stage of model reduction to produce accurate forecasts of the solution during online ROM simulations via gappy proper orthogonal decomposition (POD) . These forecasts can be used (1) to generate accurate initial guesses for the Newton solver at each time step , or (2) as an accurate coarse propagator to accelerate convergence of time-parallel methods . While both approaches reduce the computational cost incurred by time integration (by reducing the total number of Newton iterations and the wall time, respectively), neither directly reduces the temporal dimension of the ROM.

Alternatively, space–time ROMs have been devised in the reduced basis , POD–Galerkin , and ODE-residual minimization contexts. These approaches successfully reduce the temporal dimension of the underlying model by performing projection with a low-dimensional space–time basis. In addition, these methods can remove spurious temporal modes (e.g., unstable growth, artificial dissipation) from the state space, which can in principle lead to more accurate long-time responses. Further, space–time reduced-basis ROMs are equipped with error bounds that are observed to grow linearly (rather than exponentially) in the final time. While these approaches are quite promising, they exhibit several drawbacks in terms of applicability to general large-scale nonlinear dynamical systems. First, the space–time reduced-basis approaches require a space–time finite-element discretization for the high-fidelity model. Such discretizations are uncommon, as most computational models used in practice are constructed via spatial discretization (e.g., with a finite difference, finite volume, or finite element method) followed by time integration (e.g., with a linear multistep or Runge–Kutta scheme). Second, these space–time ROM approaches (with the exception of a collocation-like approach proposed in Ref. ) provide no mechanism for complexity reduction (i.e., hyper-reduction), which precludes these techniques from reducing the computational complexity in the presence of general nonlinearities. Further, the above approaches (with the exception of Ref. ) compute only a single space–time basis vector per training simulation. This can severely limit the dimensionality (and accuracy) of the resulting space–time ROM in the case of large-scale nonlinear dynamical-system models, where the number of training simulations may be limited by computational-cost considerations.

To this end, we propose a novel space–time least-squares Petrov–Galerkin (ST-LSPG) method that combines advantages of the above space–time ROM approaches, as it: (1) reduces the spatial and temporal dimensions of the dynamical system; (2) is equipped with a priori error bounds that bound the solution error by the best space–time approximation error and whose stability constants exhibit subquadratic growth in time; (3) is applicable to general nonlinear dynamical-system models; (4) is equipped with hyper-reduction to reduce the complexity in the presence of general nonlinearities; and (5) can extract multiple space–time basis vectors from each training simulation via tensor decomposition. To realize these advantages, the approach adopts aspects of both the forecasting and space–time ROM approaches described above.

Specific contributions of this work include:

A novel ST-LSPG model-reduction method for parameterized nonlinear dynamical systems (Section 4), including choices of weighting matrices to enable hyper-reduction (Section 4.3).

Several strategies for computing the ‘ingredients’ characterizing the ST-LSPG method: the space–time trial subspace via tensor decomposition (Section 5.1), the space–time residual basis in the case of space–time GNAT (Section 5.2), the sampling matrix to enable space–time hyper-reduction (Section 5.3), and the initial guess for the Gauss–Newton solver used to compute the ST-LSPG solution (Section 5.4).

A posteriori error bounds that enable the error in the ST-LSPG solution to be bounded by the value of the objective function minimized by the method (Corollary 6).

Numerical experiments that demonstrate the ability of the method to produce significant computational-cost savings relative to existing spatial-projection-based nonlinear ROMs without sacrificing accuracy (Section 7).

Ref. , which also proposed a space–time residual-minimizing projection applicable to parameterized nonlinear dynamical systems, is perhaps the most closely related work to the proposed technique. Our work can be distinguished from that contribution in several ways. First, our approach applies residual minimization to the discretized ODE (i.e., O$\Delta$E) rather than the time-continuous ODE over all time. This facilitates deriving error bounds with respect to the (fully discrete) full-order-model solution (Section 6); Ref. did not provide error bounds for the residual-minimizing approximation (it provides an error bound only for the best linear-subspace approximation). Also, Ref. enforces the sum of generalized-coordinate values to equal one; our method does not require such a constraint, which can enable lower objective-function values. Further, our approach enables multiple space–time basis vectors to be extracted from each training simulation via tensor decomposition (Section 4.3); Ref. computes only a single space–time basis vector from each training simulation, which can severely limit the dimensionality of the space–time basis in practice. Further, our proposed approach provides several mechanisms for enabling hyper-reduction, i.e., complexity reduction of the low-dimensional model (Section 4.3); Ref. proposes one approach, which is analogous to the space–time collocation method described in Section 4.3.2.

We proceed by describing the time-continuous and time-discrete representations of the full-order model in Section 2, followed by a summary of the previously developed spatial-projection-based LSPG method in Section 3. Then, we present the ST-LSPG method in Section 4, followed by proposals for the ingredients characterizing the method in Section 5. Section 6 provides error analysis, Section 7 reports numerical experiments, and Section 8 concludes the paper.

Full-order model

We begin by deriving time-continuous (ODE) and time-discrete (O$\Delta$E) formulations of the full-order model (FOM).

We consider the FOM to be a parameterized nonlinear dynamical system characterized by a parameterized system of nonlinear ODEs

2 Time-discrete representation: linear multistep methods

Here, $k(t^{n})(\leq n)$ denotes the number of steps used by the linear multistep method at time instance $n$ and the residual is defined as

Least-squares Petrov–Galerkin method

This section describes the trial subspace (Section 3.1) and projection (Section 3.2) employed by the original spatial-projection-based LSPG method .

where $\mathscr{O}\in\mathscr{H}$ is defined as $\mathscr{O}:\{t^{n}\}_{n=0}^{{N_{t}}}\rightarrow 1$. Noting that $\boldsymbol{1}_{{N_{t}}}={{\boldsymbol{g}}}(\mathscr{O})$, where $\boldsymbol{1}_{p}$ denotes a $p$-vector of ones, we also have

2 Spatial projection

The LSPG ROM computes an approximate solution by sequentially minimizing the discrete residual arising at each time instance, i.e.,

Space–time least-squares Petrov–Galerkin method

We now derive the proposed space–time least-squares Petrov–Galerkin (ST-LSPG) projection method. We begin by specifying the space–time trial subspace in Section 4.1, followed by a description of the space–time least-squares Petrov–Galerkin projection process in Section 4.2. Then, Section 4.3 describes choices for the weighting matrix to enable hyper-reduction.

Restricting the ST-LSPG state to lie in the space–time trial subspace $\boldsymbol{x}^{0}(\boldsymbol{\mu})\otimes\mathscr{O}+\mathscr{ST}$ can enable spurious temporal modes (e.g., spurious time growth or dissipation) to be removed from the set of possible solutions. More precisely, if the subspace $\mathscr{ST}$ is computed from training data as will be described in Section 5.1, then this subspace will contain only temporal modes that have been observed during the training simulations.

2 Space–time least-squares Petrov–Galerkin projection

To derive the ST-LSPG projection, we begin by defining

and define the vectorized residual as a function of the full state

and as a function of the generalized coordinates

Necessary first-order optimality conditions for Problem (4.10) correspond to stationarity of the objective function, i.e., the solution $\hat{\boldsymbol{y}}(\boldsymbol{\mu})$ satisfies

denotes the elements of the test basis. We refer to this approach as a space–time least-squares Petrov–Galerkin (ST-LSPG) projection because Eq. (4.11) corresponds to a Petrov–Galerkin projection with test basis $\{\frac{\partial\boldsymbol{{\mathsf{r}}}}{\partial\hat{w}_{i}}(\cdot;\hat{\boldsymbol{y}};\boldsymbol{\mu})\}_{i\in{\rm Nature}{{n_{st}}}}$ and also satisfies necessary conditions for the nonlinear least-squares problem (4.9)–(4.10).

As Problem (4.9)–(4.10) is simply a nonlinear least-squares problem with $\bar{z}$ equations in ${n_{st}}(\leq\bar{z})$ unknowns, we can solve it with the Gauss–Newton method, which leads to the following sequence of iterates for $k=0,\ldots,k_{\text{max}}(\boldsymbol{\mu})-1$ given an initial guess $\hat{\boldsymbol{y}}^{(0)}$:

Using the present formalism, we can also derive a (discrete) space–time Galerkin projection by enforcing Galerkin orthogonality rather than the Petrov–Galerkin orthogonality in Eq. (4.11). In particular, this space–time Galerkin method computes the solution $\hat{\boldsymbol{y}}_{\text{G}}(\boldsymbol{\mu})$ satisfying

for some prescribed metric ${\boldsymbol{\Theta}}\in\text{SPSD}(N_{s}{N_{t}})$. However, because the Galerkin solution $\hat{\boldsymbol{y}}_{\text{G}}(\boldsymbol{\mu})$ does not associate with the solution to any optimization problem in general, we do not pursue this method further. We note that this approach is the discrete counterpart to the continuous Galerkin projection proposed in Refs. ; however, these contributions effectively employ ${\boldsymbol{\Theta}}=\boldsymbol{I}_{N_{s}{N_{t}}}$ and thus provide no mechanism for hyper-reduction as we do in Section 4.3.

3 Weighting matrix and hyper-reduction

The most obvious choice for the weighting matrix is $\bar{\boldsymbol{A}}=\boldsymbol{I}_{N_{s}{N_{t}}}$; this choice simply minimizes the sum of squares of elements of the residual over both space and time and leads to $\bar{z}=N_{s}{N_{t}}$ with $\bar{\boldsymbol{A}}^{n}=\boldsymbol{I}_{N_{s}}$, $n\in{\rm Nature}{{N_{t}}}$ in Problem (4.15). This is analogous to unweighted LSPG in the spatial-projection case, in which case $\boldsymbol{A}=\boldsymbol{I}_{N_{s}}$ in Problem (3.4)–(3.5).

However, because this approach requires evaluating all $N_{s}{N_{t}}$ of the space–time residual in order to compute the objective function, it precludes significant computational-cost savings. As was also pointed out in Ref. , this bottleneck is especially cumbersome for space–time ROM approaches. For this reason, alternative choices for the weighting matrix $\bar{\boldsymbol{A}}$ can be employed that lead to a ROM whose computational complexity is independent of the full spatiotemporal dimension $N_{s}{N_{t}}$. We now describe different choices for the weighting matrix $\bar{\boldsymbol{A}}$ that lead to such hyper-reduction.

3.2 Space–time collocation

We can extend collocation hyper-reduction to the ST-LSPG context by employing weighting matrix $\bar{\boldsymbol{A}}=\bar{\boldsymbol{Z}}$ with

where $\varphi:(i,j)\mapsto i+N_{s}(j-1)$, and ${\boldsymbol{e}}_{i}$ denotes the $i$th canonical unit vector. Here, $\mathfrak{st}:=\{(\mathcal{s}_{i},\mathcal{t}_{i})\}_{i\in{\rm Nature}{\bar{n}_{z}}}\subseteq{\rm Nature}{N_{s}}\times{\rm Nature}{{N_{t}}}$ denotes the set of space–time sample indices. Again, we require ${n_{st}}\leq\bar{n}_{z}(\leq N_{s}{N_{t}})$ to ensure nonsingular residual Jacobians.

Critically, note that applying $\bar{\boldsymbol{A}}=\bar{\boldsymbol{Z}}$ in Problem (4.9)–(4.10) leads to hyper-reduction, as evaluating the objective function requires evaluating only $\bar{n}_{z}<N_{s}{N_{t}}$ elements of the spatiotemporal residual. In practice, this implies that the residual will be evaluated only at a subset of time instances and spatial degrees of freedom. Further, this can also lead to a positive semidefinite metric ${\boldsymbol{\Theta}}=\bar{\boldsymbol{A}}^{T}\bar{\boldsymbol{A}}$, as $\mathop{\operator@font rank}\nolimits(\bar{\boldsymbol{A}}^{T}\bar{\boldsymbol{A}})=\bar{n}_{z}\leq N_{s}{N_{t}}$ in this case.

3.3 Space–time GNAT

where we have defined the gappy POD residual approximations as

Computing method ingredients

This section describes particular methods for constructing the ingredients required for the ST-LSPG method, namely the space–time trial subspace $\mathscr{ST}$; the sampling matrix $\bar{\boldsymbol{Z}}$ in the case of hyper-reduction; and the residual basis ${\bar{\boldsymbol{\Phi}}}_{r}$ in the case of ST-GNAT.

Previous work on space–time model reduction constructed the space–time trial subspace simply as the span of these snapshots, i.e.,

The resulting space–time trial subspace comprises the direct sum of Kronecker products of spatial and temporal subspaces, i.e.,

The mode-1 and mode-2 unfolding of $\mathcal{X}$ can be written as

respectively, where we have defined $\boldsymbol{X}(\boldsymbol{\mu}):={{\boldsymbol{g}}}(\boldsymbol{x}(\cdot;\boldsymbol{\mu})-\boldsymbol{x}^{0}(\boldsymbol{\mu}))$. In the model-reduction literature, the matrix $\boldsymbol{X}_{(1)}$ is typically referred to as the ‘global snapshot matrix’; we refer to it in this work as the ‘spatial snapshot matrix’, as its columns comprise snapshots of the spatial solution over time and parameter variation. Similarly, we refer to $\boldsymbol{X}_{(2)}$ as the ‘temporal snapshot matrix’, as its columns comprise snapshots of the time-evolution of the solution over variation in space and parameter.

where ${n_{s}}\leq\min(N_{s},{N_{t}}n_{\text{train}})$ and $\boldsymbol{U}_{s}\equiv\left[\boldsymbol{u}_{s}^{1}\ \cdots\ \boldsymbol{u}_{s}^{{N_{t}}n_{\text{train}}}\right]$. The spatial subspace requires ${n_{s}}N_{s}$ storage. We now describe three approaches to computing the temporal subspaces $\mathscr{T}_{i}$, $i\in{\rm Nature}{{n_{s}}}$ from the state tensor $\mathcal{X}$.

1.2 Fixed temporal subspace via T-HOSVD

where $n_{t}\leq N_{s}n_{\text{train}}$, $\boldsymbol{\psi}_{j}:={\boldsymbol{h}}(\boldsymbol{\uppsi}_{j})$, and $\boldsymbol{U}_{t}\equiv\left[\boldsymbol{u}_{t}^{1}\ \cdots\ \boldsymbol{u}_{t}^{N_{s}n_{\text{train}}}\right]$. This approach is equivalent to applying the truncated higher-order SVD(T-HOSVD) to the state tensor $\mathcal{X}$, and it requires $n_{t}{N_{t}}$ storage for the temporal subspace. We note that this approach is similar to that proposed in Ref. in the context of space–time Galerkin projection performed at the time-continuous level.

1.3 Fixed temporal subspace via ST-HOSVD

1.4 Tailored temporal subspaces via ST-HOSVD

We can further tailor the temporal bases to capture the time evolution of each individual spatial basis vector. To achieve this using the ST-HOSVD, we compute the bases as

2 Space–time residual basis

We propose three methods for determining these training instances $\{\hat{\boldsymbol{y}}_{\text{res}}^{i},\boldsymbol{\mu}^{i}_{\text{res}}\}_{i\in{\rm Nature}{n_{\text{res}}}}$.

ST-LSPG ROM training iterations. This approach employs

where $\hat{\boldsymbol{y}}^{(k)}(\boldsymbol{\mu})$ corresponds to the ST-LSPG solution at the $k$th Gauss–Newton iteration (4.13)–(4.14) for some specified weighting matrix $\bar{\boldsymbol{A}}$ that does not rely on data (e.g., $\bar{\boldsymbol{A}}=\boldsymbol{I}_{N_{s}{N_{t}}}$), and $\mathcal{D}_{\text{res}}\subset\mathcal{D}$ denotes a set of training parameter instances that is in general different from $\mathcal{D}_{\text{train}}$. This case leads to $n_{\text{res}}=\sum_{\boldsymbol{\mu}\in\mathcal{D}_{\text{res}}}(k_{\text{max}}(\boldsymbol{\mu})+1)$ and requires $|\mathcal{D}_{\text{res}}|$ training simulations of the ST-LSPG ROM.

Projection of FOM training solutions. This approach employs

where $\hat{\boldsymbol{x}}(\boldsymbol{\mu})$ is defined as

Given the residual tensor $\mathcal{R}$, we can compute the associated space–time residual basis in a manner analogous to the approaches proposed in Section 5.1. That is, we can compute spatial residual bases as

with ${n_{r,s}}\leq{N_{t}}n_{\text{res}}$ and temporal residual bases either via the T-HOSVD

with $n_{r,t}\leq N_{s}n_{\text{res}}$, the ST-HOSVD

with $n_{r,t}\leq{n_{r,s}}n_{\text{res}}$, or the tailored ST-HOSVD

with $n_{r,t}^{i}\leq n_{\text{res}}$, where $\mathcal{R}(\boldsymbol{V}):=\mathcal{R}\times_{1}\boldsymbol{V}$, $\boldsymbol{\psi}_{r,j}:={\boldsymbol{h}}(\boldsymbol{\uppsi}_{r,j})$, and $\boldsymbol{\psi}^{i}_{r,j}:={\boldsymbol{h}}(\boldsymbol{\uppsi}_{r,j}^{i})$.

where $\boldsymbol{\pi}_{r,\mathcal{I}_{r}(i,j)}={\boldsymbol{h}}(\boldsymbol{\phi}_{r,i}\otimes\boldsymbol{\uppsi}_{r,j}^{i})$ and $\mathcal{I}_{r}:(i,j)\mapsto\sum_{k=1}^{i-1}n_{r,t}^{k}+j$ provides a mapping from the spatial-basis and temporal-basis indices to a space–time basis index for the residual.

3 Sampling matrix

We propose three approaches for computing the space–time sample set $\mathfrak{st}:=\{(\mathcal{s}_{i},\mathcal{t}_{i})\}_{i\in{\rm Nature}{\bar{n}_{z}}}$ that defines the residual-sampling matrix $\bar{\boldsymbol{Z}}$ in Eq. (4.17).

Greedy sampling of space–time indices. This approach selects space–time indices in a greedy manner by executing Algorithm 1, which is a space–time adaptation of the greedy method presented in Ref. that allows for oversampling to enable least-squares regression via gappy POD.

Sequential greedy sampling of spatial then temporal indices. This approach computes space–time sample indices as the Cartesian product of spatial and temporal samples, i.e., $\mathfrak{st}=\mathfrak{s}\times\mathfrak{t}$. First, the approach selects spatial indices $\mathfrak{s}$ by executing Algorithm 3 with inputs ${\bar{\boldsymbol{\Phi}}}_{r}$, the desired number of spatial samples $\bar{n}_{s}$, and $\mathfrak{t}={\rm Nature}{{N_{t}}}$ (i.e., full temporal sampling). Then, the method selects temporal indices $\mathfrak{t}$ by executing Algorithm 2 with inputs ${\bar{\boldsymbol{\Phi}}}_{r}$, the desired number of temporal samples $\bar{n}_{t}$, and $\mathfrak{s}$ computed from Algorithm 3.

Sequential greedy sampling of temporal then spatial indices. This approach also computes space–time sample indices as $\mathfrak{st}=\mathfrak{s}\times\mathfrak{t}$. First, the approach selects temporal indices $\mathfrak{t}$ by executing Algorithm 2 with inputs ${\bar{\boldsymbol{\Phi}}}_{r}$, the desired number of temporal samples $\bar{n}_{t}$, and $\mathfrak{s}={\rm Nature}{N_{s}}$ (i.e., full spatial sampling). Then, the method selects spatial indices $\mathfrak{s}$ by executing Algorithm 3 with inputs ${\bar{\boldsymbol{\Phi}}}_{r}$, the desired number of spatial samples $\bar{n}_{s}$, and $\mathfrak{t}$ computed from Algorithm 2.

We note that enforcing $\mathfrak{st}=\mathfrak{s}\times\mathfrak{t}$ as in approaches 2 and 3 above comes with a practical advantage. Namely, a single sample mesh —which is tasked with computing spatial samples associated with $\mathfrak{s}$—can be employed for all sampled time instances $\mathfrak{t}$.

4 Initial guess

One practical challenge of ST-LSPG relative to (spatial-projection-based) LSPG is devising an accurate initial guess $\hat{\boldsymbol{y}}^{(0)}$ for the Gauss–Newton iterations (4.13)–(4.14). In LSPG, the initial guess employed when solving Problem (3.4)–(3.5) at a given time instance $t^{n}$ using the Gauss–Newton method can be set to the solution from the previous time instance, i.e., $\hat{\boldsymbol{x}}(t^{n(0)};\boldsymbol{\mu})=\hat{\boldsymbol{x}}(t^{n-1};\boldsymbol{\mu})$. This choice typically leads to rapid convergence due to the fact that the state undergoes limited variation between time instances, particularly for small time steps $\Delta t^{n}$. Alternatively, accurate initial guesses based on polynomial extrapolation or forecasting can be employed to further improve convergence.

However, in ST-LSPG, deriving an accurate initial guess $\hat{\boldsymbol{y}}^{(0)}(\boldsymbol{\mu})$ is less straightforward. We propose computing $\hat{\boldsymbol{y}}^{(0)}(\boldsymbol{\mu})$ as an interpolant of the generalized coordinates $\hat{\boldsymbol{y}}(\boldsymbol{\mu})$ in the parameter space. That is, given the training parameter instances $\mathcal{D}_{\text{train}}\subset\mathcal{D}$ for which the FOM has been solved, we can compute the projection of these solutions onto the space–time trial subspace as $\{\hat{\boldsymbol{x}}(\boldsymbol{\mu})\}_{\boldsymbol{\mu}\in\mathcal{D}_{\text{train}}}$ with $\hat{\boldsymbol{x}}(\boldsymbol{\mu})$ defined in Eq. (5.17). Then, we can compute $\hat{\boldsymbol{y}}^{(0)}(\boldsymbol{\mu})$ via interpolation (or least-squares regression) in the parameter space $\mathcal{D}$ using data $\{\hat{\boldsymbol{x}}(\boldsymbol{\mu})\}_{\boldsymbol{\mu}\in\mathcal{D}_{\text{train}}}.$

Error analysis

respectively. We begin by stating assumptions that will be leveraged in subsequent analyses:

There exists a constant $L_{\boldsymbol{f}}>0$ such that

The time step $\Delta t$ is sufficiently small such that

where $\boldsymbol{I}=\boldsymbol{I}_{N_{s}}$ here and $\sigma_{\text{min}}(\boldsymbol{A})$ and $\sigma_{\text{max}}(\boldsymbol{A})$ denote the minimum and maximum singular values of the matrix $\boldsymbol{A}$, respectively.

Under Assumption A1, the linear multistep residual is also Lipschitz continuous, i.e.,

Defining $\bar{\boldsymbol{f}}:\boldsymbol{w}\mapsto{\boldsymbol{h}}(\boldsymbol{f}(\boldsymbol{w}(\cdot),\cdot))$ and noting that $\|\bar{\boldsymbol{f}}(\boldsymbol{w})\|_{2}^{2}=\sum_{n=1}^{{N_{t}}}\|\boldsymbol{f}(\boldsymbol{w}(t^{n}),t^{n})\|_{2}^{2}\leq L_{\boldsymbol{f}}^{2}\sum_{n=1}^{{N_{t}}}\|\boldsymbol{w}(t^{n})\|_{2}^{2}=L_{\boldsymbol{f}}^{2}\|\boldsymbol{w}\|_{2}^{2}$, we have

i.e., the Lipschitz constant of $\bar{\boldsymbol{f}}$ is identical to that of $\boldsymbol{f}$. Further noting that

where $\boldsymbol{b}(t^{n})=\alpha_{n}^{n}\boldsymbol{x}^{0}-\Delta t\beta_{n}^{n}\boldsymbol{f}(\boldsymbol{x}^{0})$, we have from the triangle inequality

Under Assumptions and A1 and A2, the linear multistep residual is also inverse Lipschitz continuous, i.e.,

Applying the reverse triangle inequality and employing Assumption A2 yields

which directly leads to the desired result. ∎

Under Assumptions A1, A2, and A3, the error in the ST-LSPG solution at any time instance can be bounded by the best approximation error as

Substituting in the definitions of $L_{\boldsymbol{r}}$ and $K_{\boldsymbol{r}}$ from Eqs. (6.3) and (6.8), respectively, yields the stated result. ∎

Under Assumptions A1, A2, and A3, the error in the ST-LSPG solution at any time instance can be bounded by the best approximation error as

We now provide simplified variants of these error bounds in the case of unweighted LSPG (Section 4.3.1) for which $\bar{\boldsymbol{A}}^{n}=\boldsymbol{I}_{N_{s}}$.

If $\bar{\boldsymbol{A}}=\boldsymbol{I}_{N_{s}{N_{t}}}$, then under Assumptions A1 and A2, the error in the ST-LSPG solution at any time instance can be bounded by the best approximation error as

where we define the Lebesgue constant for a given time integrator and time step $\Delta t$ as

Proofs follows trivially from Theorems 3 and 4 by substituting $\bar{\boldsymbol{A}}=\boldsymbol{I}_{N_{s}{N_{t}}}$ and noting that Assumption A3 is automatically satisfied for this choice of weighting matrix $\bar{\boldsymbol{A}}$, as $P=1$ in this case. ∎

Under Assumptions A1, A2, and A3, the error in the any approximation $\boldsymbol{w}\in\mathcal{ST}$ can be bounded by the computed residual norm as

Further, the ST-LSPG solution is the particular solution for which this error bound is minimized, i.e.,

By invoking Assumption A3, Lemma 2, and norm equivalence $\|\boldsymbol{w}\|_{2}\geq\|\boldsymbol{w}\|_{\infty}$, we can derive

which yields the first desired result. The second result follows from applying the optimality property of the ST-LSPG solution (6.1). ∎

Numerical experiments

This section compares the performance of the following methods:

FOM. This model corresponds to the full-order model, i.e., the solution satisfying Eq. (2.2).

LSPG ROM. This model corresponds to the unweighted LSPG ROM, i.e., the solution that satisfies Eq. (3.4) with $\boldsymbol{A}=\boldsymbol{I}_{N_{s}}$.

GNAT ROM. This model corresponds to the GNAT ROM, i.e., the solution that satisfies Eq. (3.4) with $\boldsymbol{A}=(\boldsymbol{Z}\boldsymbol{\Phi}_{r})^{+}\boldsymbol{Z}$. Algorithm 5 in Ref. is used to construct the sampling matrix $\boldsymbol{Z}$.

ST-LSPG-1 ROM. This model corresponds to the unweighted ST-LSPG ROM, i.e., the solution that satisfies Eq. (4.9) with $\bar{\boldsymbol{A}}=\boldsymbol{I}_{N_{s}{N_{t}}}$. The method is also characterized by the following:

Tailored temporal state subspaces computed according to Eqs. (5.12)–(5.13).

As described in Section 5.4, interpolation to compute the initial guess. For this, we employ interpolation using linear radial basis functions as described in Ref. .

ST-LSPG-2 ROM. This model is identical to the ST-LSPG-1 ROM except that it employs fixed temporal subspaces computed according to Eqs. (5.10)–(5.11).

ST-GNAT-1 ROM. This model corresponds to the ST-GNAT ROM, i.e., the solution that satisfies Eq. (4.9) with $\bar{\boldsymbol{A}}=(\bar{\boldsymbol{Z}}{\bar{\boldsymbol{\Phi}}}_{r})^{+}\bar{\boldsymbol{Z}}$. Otherwise, it is identical to the ST-LSPG-1 ROM with the additional following attributes:

Tailored temporal residual subspaces computed according to Eqs. (5.24)–(5.25).

Method 1 in Section 5.2 to generate space–time residual samples, where $\mathcal{D}_{\text{res}}=\mathcal{D}_{\text{train}}$.

Method 3 in Section 5.3 to construct the sampling matrix.

ST-GNAT-2 ROM. This model is identical to the ST-GNAT-1 ROM except that it employs a fixed temporal state subspace computed according to Eqs. (5.10)–(5.11).

and we measure its computational cost in terms of the wall time incurred by the ROM relative to that incurred by the FOM; the speedup is the reciprocal of the relative wall time. All timings are obtained by performing calculations on an Intel(R) Xeon(R) CPU E5-2670 @ 2.60 GHz, 31.4 GB RAM using the MORTestbed in MATLAB. All reported timings are averaged over five simulations.

We first consider the parameterized inviscid Burgers’ equation described in Ref. , which corresponds to the following initial boundary value problem for $x\in$ and $t\in[0,T]$ with $T=0.5$:

After applying Godunov’s scheme for spatial discretization with 100 control volumes, Eqs. (7.2) leads to a parameterized initial-value ODE problem consistent with Eq. (2.1) with $N_{s}=100$ spatial degrees of freedom. For time discretization, we employ the backward Euler scheme, which is a linear multistep method characterized by $k(t^{n})=1$, $\alpha_{0}^{n}=\beta_{0}^{n}=1$, $\alpha_{1}^{n}=-1$, $\beta_{1}^{n}=0$, $n\in{\rm Nature}{{N_{t}}}$. We employ a uniform time step of $\Delta t=2.5\times 10^{-4}$, leading to ${N_{t}}=2000$ time instances. For this problem, all ROMs employ a training set $\mathcal{D}_{\text{train}}=\{1.2,1.3,1.4,1.5\}\times\{0.02,0.025\}$ such that $n_{\text{train}}=8$ at which the FOM is solved.

We emphasize that (1) all assessed models employ the same time discretization; this includes the ST-LSPG and ST-GNAT ROMs, as the space–time residual is defined from this time discretization, and (2) we do not consider adaptive time-step selection. Future work will investigate the effect of different time integrators—including those that employ adaptive time-step selection—on the relative performance of the methods.

Figure 2 plots a selection of spatial and temporal modes computed using the three different techniques proposed in Section 5.1. Note that the fixed temporal modes are nearly identical, regardless of whether the T-HOSVD or ST-HOSVD is employed. Thus, because the ST-HOSVD is significantly less computationally expensive, we no longer consider the fixed modes computed with T-HOSVD, which was the approach considered in Ref. . On the other hand, the tailored temporal modes are significantly different from the fixed temporal modes. Further, they appear to be well suited for their respective spatial modes, as the temporal bases for higher-index spatial POD modes exhibit higher frequencies, which is consistent with previous studies (e.g., Ref. ).

1.2 Model predictions

We now compare the methods for fixed values of their parameters, and for two randomly selected online points $\boldsymbol{\mu}^{1}=(1.35,0.0229)\not\in\mathcal{D}_{\text{train}}$ and $\boldsymbol{\mu}^{2}=(1.45,0.0201)\not\in\mathcal{D}_{\text{train}}$. Table 1 reports the method parameter values and the associated performance of the methods. Figure 3 reports snapshots of the methods’ responses for $t\in\{0,0.1665,0.3332,0.5\}$.

First, note that all ROMs generate accurate responses for this particular combination of parameters, as the relative errors are less than 1% in all cases. Second, note that LSPG generates the most accurate responses, but fails to generate any speedup due to its lack of hyper-reduction. GNAT also fails to generate speedup in this case due to the relatively small spatial dimension of the FOM and the larger number of Newton iterations required for convergence relative to LSPG. The proposed ST-LSPG methods incur slightly larger errors than the LSPG method, but they do so with orders of magnitude fewer space–time degrees of freedom. This highlights the promise of performing projection in both space and time: the dimensionality of the problem can be significantly reduced while retaining high levels of accuracy. However, due to their lack of hyper-reduction, the ST-LSPG methods do not generate speedups. Finally, by employing hyper-reduction, the ST-GNAT methods generate very accurate predictions with significant speedups. We note that ST-LSPG-1 and ST-GNAT-1 exhibit better overall performance than ST-LSPG-2 and ST-GNAT-2, respectively; this suggests that employing tailored temporal subspaces enables similar accuracy to be achieved using far fewer degrees of freedom, as each temporal basis vector is tailored to its associated spatial basis vector.

1.3 Method-parameter study

This section compares the performance of the ROM methods across a variation of all method parameters. This study is essential to objectively compare the methods, as the particular method-parameter values employed in Section 7.1.2 did not necessarily yield optimal performance for a given method. For this reason, we subject each model to a parameter study wherein each model parameter is varied between specified limits; Table 2 reports the tested parameter values for each method. We consider all elements in the resulting set if they satisfy the following constraints: $1.5{n_{s}}\leq\bar{n}_{r}\leq n_{z}$ for GNAT and $1.5{n_{st}}\leq\bar{n}_{r}\leq\bar{n}_{s}\bar{n}_{t}$ for ST-GNAT-1 and ST-GNAT-2. From these results, we then construct a Pareto front for each method, which is characterized by the method parameters that minimize the competing objectives of relative error and relative wall time.

Figure 4 reports these Pareto fronts for the two online points, as well as an ‘overall’ Pareto front that selects the Pareto-optimal methods across all parameter variations. Table 3 reports values of the method parameters that yielded Pareto-optimal performance. The proposed ST-GNAT-1 method is Pareto-optimal for relative wall times less than one (i.e., faster than the FOM simulation). While the proposed ST-GNAT-2 method does produce speedups, it is dominated by ST-GNAT-1; this provides further evidence of the advantage of employing a tailored relative to a fixed temporal basis. We note that the worst-performing methods correspond to the ST-LSPG-1, and ST-LSPG-2 methods, as their lack of hyper-reduction leads to significant wall times that far exceed that of the FOM. Further, we note for a fixed error below a certain threshold, the ST-GNAT-1 method is nearly two orders of magnitude faster than the original GNAT method; this can be attributed to the fact that this approach reduces both the spatial and temporal complexities of the FOM. Finally, we note that because the spatial trial subspace employed by LSPG and GNAT has a (relatively large) spatiotemporal dimension of ${n_{s}}{N_{t}}$, while the space–time trial subspace employed by ST-LSPG and ST-LSPG has a (relatively small) spatiotemporal dimension of ${n_{st}}(\ll{n_{s}}{N_{t}})$, the LSPG and GNAT methods are able to generate smaller errors than the space–time methods. However, this is achieved at significant computational cost that exceeds that of the FOM in this case (i.e., relative wall times greater than one). Thus, for this problem, LSPG is Pareto-optimal and outperforms the space–time ROMs for relative errors less than $10^{-6}$, although this regime is not useful because it incurs relative wall times greater than one.

2 Quasi 1D Euler equation

We now consider a parameterized quasi-1D Euler equation associated with modeling inviscid compressible flow in a one-dimensional converging–diverging nozzle with a continuously varying cross-sectional area [33, Chapter 13]; Figure 5 depicts the problem geometry. The governing system of nonlinear PDEs is

Here, $\rho$ denotes density, $u$ denotes velocity, $p$ denotes pressure, $\epsilon$ denotes potential energy per unit mass, $e$ denotes total energy density, $\gamma$ denotes the specific heat ratio, and $A$ denotes the converging–diverging nozzle cross-sectional area. We employ a specific heat ratio of $\gamma=1.3$, a specific gas constant of $R=355.4$ $\text{m}^{2}/\text{s}^{2}/\text{K}$, a total temperature of $T_{t}=300$ K, and a total pressure of $p_{t}=10^{6}$ $\text{N}/\text{m}^{2}$. The cross-sectional area $A(x)$ is determined by a cubic spline interpolation over the points $(x,A(x))\in\{(0,0.2),(0.25,0.173),(0.5,0.17),(0.75,0.173),(1,0.2)\}$, which results in

We assume a perfect gas that obeys the ideal gas law (i.e., $p=\rho RT$). The initial flow field is created in several steps. First, the following isentropic relations are used to generate a zero pressure-gradient flow field:

Applying a finite-volume spatial discretization with 50 equally spaced control volumes and fully implicit boundary conditions leads to a parameterized system of nonlinear ODEs consistent with Eq. (2.1) with $N_{s}=150$ spatial degrees of freedom. The Roe flux difference vector splitting method is used to compute the flux at each intercell face [33, Chapter 9]. For time discretization, we again apply the backward Euler scheme and a uniform time step of $\Delta t=0.001$ s, leading to ${N_{t}}=600$.

For this problem, we use the following two parameters: the pressure factor $\mu_{1}=P_{\text{exit}}$ and the Mach number at the middle of the nozzle $\mu_{2}=M_{m}$. All ROMs employ a training set at which the FOM is solved of $\mathcal{D}_{\text{train}}=\{1.7+0.01i\}_{i=0}^{3}\times\{1.7,1.72\}$ such that $n_{\text{train}}=8$.

Figure 6 plots several spatial and temporal modes computed using the three different techniques proposed in Section 5.1. As with the Burgers equation, the ‘fixed’ temporal modes are nearly identical, regardless of whether the T-HOSVD or ST-HOSVD is employed, rendering the ST-HOSVD more appealing due to its reduced computational cost. In addition, the tailored temporal modes are significantly different, with the temporal basis exhibiting higher frequencies for higher-index spatial modes as expected.

2.2 Model predictions

We now compare the methods for fixed values of their parameters, and for two randomly selected online points $\boldsymbol{\mu}^{1}=(1.7125,1.71)\not\in\mathcal{D}_{\text{train}}$ and $\boldsymbol{\mu}^{2}=(1.7225,1.705)\not\in\mathcal{D}_{\text{train}}$. Table 4 reports the method parameter values and the associated performance of the methods. Figure 7 reports snapshots of the methods’ responses for $t\in\{0,T\}$.

Conclusions are similar to those derived from the Burgers’ equation results. First, note that all ROMs except for GNAT in the case of $\boldsymbol{\mu}^{1}$ generate accurate responses, as the relative errors are less than 3% in all cases. Second, as before, LSPG generates the most accurate responses, but fails to generate any speedup due to its lack of hyper-reduction. The proposed ST-LSPG methods incur sub-2% errors, but they do so with orders of magnitude fewer spatiotemporal degrees of freedom relative to the LSPG and GNAT methods, which highlights the promise of performing projection in both space and time. Again, as these methods do not employ hyper-reduction, they do not generate speedups. Finally, the ST-GNAT methods generate both accurate predictions with significant speedups. We note that ST-GNAT-1 performs better than ST-GNAT-2, providing further evidence of the ability of tailored bases to produce accurate responses with fewer degrees of freedom.

2.3 Method-parameter study

We again compare the performance of the ROM methods across a wide variation of all method parameters. Table 5 reports the tested parameter values for each method. We consider all elements in the resulting set if they satisfy constraints $1.5{n_{s}}\leq\bar{n}_{r}\leq n_{z}$ for GNAT and $1.5{n_{st}}\leq\bar{n}_{r}\leq\bar{n}_{s}\bar{n}_{t}$ for ST-GNAT-1 and ST-GNAT-2. From these results, we then construct a Pareto front for each method, which is characterized by the method parameters that minimize the competing objectives of relative error and relative wall time.

Figure 8 reports these Pareto fronts for the two online points, as well as an overall Pareto front that selects the Pareto-optimal methods across all parameter variations. Table 6 reports values of the method parameters that yielded Pareto-optimal performance. These results show that—as before—the proposed ST-GNAT-1 method is Pareto optimal for relative wall time less than 0.9 and relative errors less than 20%. While the proposed ST-GNAT-2 method produces speedups, it is again dominated by ST-GNAT-1, further highlighting the advantage of tailored versus fixed temporal bases. Again, the worst-performing methods correspond to the ST-LSPG-1, and ST-LSPG-2 methods, as their lack of hyper-reduction leads to significant wall times that far exceed that of the FOM. Further, we note that for a fixed error below a certain threshold, the ST-GNAT-1 method is over one order of magnitude faster than the original GNAT method; this can be attributed to the fact that this approach reduces both the spatial and temporal complexities of the FOM. Finally, we again note that LSPG is Pareto-optimal and outperforms the space–time ROMs for extremely small relative errors less than approximately $3\times 10^{-5}$ due to the higher spatiotemporal dimensionality of the spatial trial subspace; however, this regime is not useful for this problem, as it leads to LSPG models roughly as expensive as the FOM (i.e., relative wall times near one).

Conclusions

its ability to reduce both the spatial and temporal dimensions of the dynamical system (Remark 4.1),

a priori error bounds that bound the solution error by the best space–time approximation error and whose stability constants exhibit subquadratic growth in time (Section 6),

applicability to general nonlinear dynamical systems,

hyper-reduction that reduces the complexity in the presence of general nonlinearities (Section 4.3), and

its ability to extract multiple space–time basis vectors from each training simulation via tensor decomposition (Section 6).

In addition to introducing the novel ST-LSPG method, this work proposed specific approaches to computing the method’s ingredients: the space–time trial subspace (Section 5.1), the space–time residual basis in the case of ST–GNAT (Section 5.2), the sampling matrix in the case of hyper-reduction (Section 5.3), and the initial guess used in the Gauss–Newton method applied to solve the nonlinear least-squares problem (Section 5.4). Numerical experiments demonstrated the ability of the proposed method to generate orders-of-magnitude speedups over existing spatial-projection-based ROMs without sacrificing accuracy.

Future work entails implementing the method in parallel computational-mechanics codes, devising techniques to reduce the amount of storage required for the state and residual tensors, and assessing the effect of different time integrators—as well as adaptive time-stepping—on method performance.

Acknowledgments

The authors gratefully acknowledge Tamara Kolda and Grey Ballard for insightful discussions that led to the tensor-decomposition approaches for computing the space–time trial subspace. The authors also acknowledge Professors Benjamin Peherstorfer and Masayuki Yano for useful comments received at the Model Reduction for Parametrized Systems (MoRePaS) III workshop and the 2017 SIAM Conference on Computational Science and Engineering, respectively. The authors also gratefully acknowledge the helpful comments provided by the anonymous reviewers. This work was performed at Sandia National Laboratories and was supported by the LDRD program (project 190968). Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and 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-NA-0003525. Lawrence Livermore National Laboratory is operated by Lawrence Livermore National Security, LLC, for the U.S. Department of Energy, National Nuclear Security Administration under Contract DE-AC52-07NA27344.

References