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

TL;DR

This paper introduces a novel space-time least-squares Petrov-Galerkin projection method for nonlinear model reduction, simultaneously reducing spatial and temporal dimensions to improve efficiency and stability in dynamical system simulations.

Contribution

The paper presents a new space-time projection approach that combines tensor decompositions and residual minimization, enabling hyper-reduction and improved stability over traditional methods.

Findings

Achieves significant computational savings in fluid dynamics models.

Removes spurious temporal modes, enhancing model stability.

Provides error bounds with slower growth over time.

Abstract

This work proposes a space-time least-squares Petrov-Galerkin (ST-LSPG) projection method for model reduction of nonlinear dynamical systems. In contrast to typical nonlinear model-reduction methods that first apply (Petrov-)Galerkin projection in the spatial dimension and subsequently apply time integration to numerically resolve the resulting low-dimensional dynamical system, the proposed method applies projection in space and time simultaneously. To accomplish this, the method first introduces a low-dimensional space-time trial subspace, which can be obtained by computing tensor decompositions of state-snapshot data. The method then computes discrete-optimal approximations in this space-time trial subspace by minimizing the residual arising after time discretization over all space and time in a weighted $\ell^2$-norm. This norm can be defined to enable complexity reduction (i.e., hyper-reduction) in time, which leads to space-time collocation and space-time GNAT variants of the ST-LSPG method. Advantages of the approach relative to typical spatial-projection-based nonlinear model reduction methods such as Galerkin projection and least-squares Petrov-Galerkin projection include: (1) a reduction of both the spatial and temporal dimensions of the dynamical system, (2) the removal of spurious temporal modes (e.g., unstable growth) from the state space, and (3) error bounds that exhibit slower growth in time. Numerical examples performed on model problems in fluid dynamics demonstrate the ability of the method to generate orders-of-magnitude computational savings relative to spatial-projection-based reduced-order models without sacrificing accuracy.

Tables

Table 1. Table 1: Burgers’ equation . ROM method performance for fixed method parameters at randomly selected online points 𝝁 1 = ( 1.35 , 0.0229 ) ∉ 𝒟 train superscript 𝝁 1 1.35 0.0229 subscript 𝒟 train \boldsymbol{\mu}^{1}=(1.35,0.0229)\not\in\mathcal{D}_{\text{train}} and 𝝁 2 = ( 1.45 , 0.0201 ) ∉ 𝒟 train superscript 𝝁 2 1.45 0.0201 subscript 𝒟 train \boldsymbol{\mu}^{2}=(1.45,0.0201)\not\in\mathcal{D}_{\text{train}} .

method	LSPG	GNAT	ST-LSPG-1	ST-LSPG-2	ST-GNAT-1	ST-GNAT-2
$n_s$	15	15	15	15	15	15
$n_z$		55
$n_r$		55
$n_t^i$			2		2
$n_t$				20		20
${\overline{n}}_s$					30	30
${\overline{n}}_t$					120	120
$n_{r,s}$					100	100
$n_{r,t}^i$					3	10
spatiotemporal dimension	$3\times {10}^4$	$3\times {10}^4$	30	300	30	300
relative error for ${𝝁}^1$	0.00074	0.011	0.0025	0.0011	0.0058	0.0063
speedup for ${𝝁}^1$	0.82	0.34	0.34	0.079	7.19	1.78
relative error for ${𝝁}^2$	0.0012	0.017	0.0038	0.0040	0.0077	0.0082
speedup for ${𝝁}^2$	0.80	0.39	0.33	0.080	6.22	1.72

Table 2. Table 2: Burgers’ equation . Parameters used for the method-parameter study. The set of tested parameters comprises the Cartesian product of the specified parameter sets. We consider all elements in the resulting set if they satisfy the following guidelines: 1.5 ​ n s ≤ n ¯ r ≤ n z 1.5 subscript 𝑛 𝑠 subscript ¯ 𝑛 𝑟 subscript 𝑛 𝑧 1.5{n_{s}}\leq\bar{n}_{r}\leq n_{z} for GNAT and 1.5 ​ n s ​ t ≤ n ¯ r ≤ n ¯ s ​ n ¯ t 1.5 subscript 𝑛 𝑠 𝑡 subscript ¯ 𝑛 𝑟 subscript ¯ 𝑛 𝑠 subscript ¯ 𝑛 𝑡 1.5{n_{st}}\leq\bar{n}_{r}\leq\bar{n}_{s}\bar{n}_{t} for ST-GNAT-1 and ST-GNAT-2.

method	LSPG	GNAT	ST-LSPG-1	ST-LSPG-2	ST-GNAT-1	ST-GNAT-2
$n_s$	${\left\{10\times i\right\}}_{i=1}^5$	${\left\{10\times i\right\}}_{i=1}^5$	${\left\{10\times i\right\}}_{i=1}^5$	${\left\{10\times i\right\}}_{i=1}^5$	${\left\{10\times i\right\}}_{i=1}^5$	${\left\{10\times i\right\}}_{i=1}^5$
$n_z$		$\{20,30,40,60,80,90\}$
$n_r$		$\{20,30,40,60,80,90\}$
$n_t^i$			$\{3,4,5,6,8\}$		$\{3,4,5,6,8\}$
$n_t$				$\{5,10,20,30\}$		$\{5,10,20,30\}$
${\overline{n}}_s$					$\{30,40\}$	$\{30,40,60,70,80\}$
${\overline{n}}_t$					$\{60,120\}$	$\left\{120\right\}$
$n_{r,s}$					$\left\{100\right\}$	$\left\{100\right\}$
$n_{r,t}^i$					$\left\{3\right\}$	$\left\{10\right\}$

Table 3. Table 3: Burgers’ equation . Parameter values yielding Pareto-optimal performance for the ST-GNAT-1 method. Figure 4 provides labels.

label	$n_s$	$n_t^i$	${\overline{n}}_s$	${\overline{n}}_t$	$n_{r,s}$	$n_{r,t}^i$
1-(1)	10	3	30	60	100	3
1-(2)	10	4	30	60	100	3
1-(3)	30	3	40	120	100	3
2-(1)	20	3	30	60	100	10
2-(2)	30	3	40	120	100	10

Table 4. Table 4: Quasi-1D Euler equation . ROM method performance for fixed method parameters at randomly selected online points 𝝁 1 = ( 1.7125 , 1.71 ) ∉ 𝒟 train superscript 𝝁 1 1.7125 1.71 subscript 𝒟 train \boldsymbol{\mu}^{1}=(1.7125,1.71)\not\in\mathcal{D}_{\text{train}} and 𝝁 2 = ( 1.7225 , 1.705 ) ∉ 𝒟 train superscript 𝝁 2 1.7225 1.705 subscript 𝒟 train \boldsymbol{\mu}^{2}=(1.7225,1.705)\not\in\mathcal{D}_{\text{train}} .

method	LSPG	GNAT	ST-LSPG-1	ST-LSPG-2	ST-GNAT-1	ST-GNAT-2
$n_s$	50	50	50	50	50	50
$n_z$		145
$n_r$		145
$n_t^i$			3		3
$n_t$				30		30
${\overline{n}}_s$					120	140
${\overline{n}}_t$					20	100
$n_{r,s}$					150	150
$n_{r,t}^i$					10	10
spatiotemporal dimension	$3\times {10}^4$	$3\times {10}^4$	150	$1.5\times {10}^3$	150	$1.5\times {10}^3$
relative error for ${𝝁}^1$	$7.78\times {10}^{-6}$	0.55	0.012	$6.3\times {10}^{-4}$	0.0023	0.0048
speedup for ${𝝁}^1$	0.77	1.04	0.84	0.58	21.79	0.49
relative error for ${𝝁}^2$	$8.31\times {10}^{-6}$	0.026	0.0076	0.0021	0.0025	0.0040
speedup for ${𝝁}^2$	0.81	1.01	0.85	0.41	22.52	2.79

Table 5. Table 5: Quasi-1D Euler equation . Parameters used for the method-parameter study. The set of tested parameters comprises the Cartesian product of the specified parameter sets. We consider all elements in the resulting set if they satisfy the following guidelines: 1.5 ​ n s ≤ n ¯ r ≤ n z 1.5 subscript 𝑛 𝑠 subscript ¯ 𝑛 𝑟 subscript 𝑛 𝑧 1.5{n_{s}}\leq\bar{n}_{r}\leq n_{z} for GNAT and 1.5 ​ n s ​ t ≤ n ¯ r ≤ n ¯ s ​ n ¯ t 1.5 subscript 𝑛 𝑠 𝑡 subscript ¯ 𝑛 𝑟 subscript ¯ 𝑛 𝑠 subscript ¯ 𝑛 𝑡 1.5{n_{st}}\leq\bar{n}_{r}\leq\bar{n}_{s}\bar{n}_{t} for ST-GNAT-1 and ST-GNAT-2.

method	LSPG	GNAT	ST-LSPG-1	ST-LSPG-2	ST-GNAT-1	ST-GNAT-2
$n_s$	${\left\{10\times i\right\}}_{i=1}^6$	${\left\{10\times i\right\}}_{i=1}^6$	${\left\{10\times i\right\}}_{i=1}^6$	${\left\{10\times i\right\}}_{i=1}^6$	${\left\{10\times i\right\}}_{i=1}^6$	${\left\{10\times i\right\}}_{i=1}^6$
$n_z$		${\left\{10\times i\right\}}_{i=2}^{10}\cup \{120,145\}$
$n_r$		${\left\{10\times i\right\}}_{i=2}^{10}\cup \{120,145\}$
$n_t^i$			${\left\{i\right\}}_{i=3}^8$		${\left\{i\right\}}_{i=3}^8$
$n_t$				${\left\{10\times i\right\}}_{i=2}^5$		${\left\{10\times i\right\}}_{i=2}^5$
${\overline{n}}_s$					$\left\{120\right\}$	$\left\{140\right\}$
${\overline{n}}_t$					$\{10,20,30,40,60\}$	$\{30,50,100,150\}$
$n_{r,s}$					150	150
$n_{r,t}^i$					10	10

Table 6. Table 6: Quasi-1D Euler equation . Parameter values yielding Pareto-optimal performance for the ST-GNAT-1 method. Figure 8 provides labels.

label	$n_s$	$n_t^i$	${\overline{n}}_s$	${\overline{n}}_t$	$n_{r,s}$	$n_{r,t}^i$
1-(1)	10	3	120	10	150	10
1-(2)	60	3	120	10	150	10
1-(3)	30	4	120	10	150	10
1-(4)	50	8	120	10	150	10
2-(1)	30	3	120	10	150	10
2-(2)	50	4	120	10	150	10
2-(3)	60	7	120	10	150	10

Equations

$$\frac{d{\mathbf{x}}^{\star }}{dt}=\mathbf{f}({\mathbf{x}}^{\star },t;\mathbf{\mu }),{\mathbf{x}}^{\star }(0;\mathbf{\mu })={\mathbf{x}}^0(\mathbf{\mu }),$$

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

$$\begin{array}{rl}{\mathbf{r}}^n& :({\mathbf{w}}^0,\dots ,{\mathbf{w}}^{k(t^n)};\mathbf{\nu })\mapsto \sum _{j=0}^{k(t^n)}{\alpha }_j^n{\mathbf{w}}^j-\Delta t^n\sum _{j=0}^{k(t^n)}{\beta }_j^n\mathbf{f}({\mathbf{w}}^j,t^{n-j};\mathbf{\nu })\\ & :{\mathbb{R}}^{N_s}\otimes {\mathbb{R}}^{k(t^n)+1}\times \mathcal{D}\to {\mathbb{R}}^{N_s},\end{array}$$

$$\mathbf{x}(\cdot ;\mathbf{\mu })\in {\mathbb{R}}^{N_s}\otimes \mathcal{H}$$

$$\begin{array}{rl}\mathbf{g}& :\mathbf{u}\mapsto \left[\mathbf{u}(t^1)\cdots \mathbf{u}(t^{N_t})\right]\\ & :{\mathbb{R}}^p\otimes \mathcal{H}\to {\mathbb{R}}^p\otimes {\mathbb{R}}^{N_t},\end{array}$$

$$\mathbf{g}(\mathbf{x}(\cdot ;\mathbf{\mu }))=\left[\mathbf{x}(t^1;\mathbf{\mu })\cdots \mathbf{x}(t^{N_t};\mathbf{\mu })\right]\in {\mathbb{R}}^{N_s}\otimes {\mathbb{R}}^{N_t}.$$

$$\tilde{\mathbf{x}}(t^n;\mathbf{\mu })={\mathbf{x}}^0(\mathbf{\mu })+\sum _{i=1}^{n_s}{\mathbf{\varphi }}_i{\hat{x}}_i(t^n;\mathbf{\mu })$$

$$\tilde{\mathbf{x}}(\cdot ;\mathbf{\mu })\in {\mathbf{x}}^0(\mathbf{\mu })\otimes \mathcal{O}+\mathcal{S}\otimes \mathcal{H}\subseteq {\mathbb{R}}^{N_s}\otimes \mathcal{H},$$

$$\mathbf{g}(\tilde{\mathbf{x}}(\cdot ;\mathbf{\mu }))\in {\mathbf{x}}^0(\mathbf{\mu })\otimes 1_{N_t}+\mathcal{S}\otimes {\mathbb{R}}^{N_t}\subseteq {\mathbb{R}}^{N_s}\otimes {\mathbb{R}}^{N_t}.$$

$$\tilde{\mathbf{x}}(t^n;\mathbf{\mu })=\underset{\mathbf{v}\in {\mathbf{x}}^0(\mathbf{\mu })+\mathcal{S}}{\mathrm{argmin}}{\Vert \mathbf{A}{\mathbf{r}}^n(\mathbf{v},\tilde{\mathbf{x}}(t^{n-1};\mathbf{\mu }),\dots ,\tilde{\mathbf{x}}(t^{n-k(t^n)};\mathbf{\mu });\mathbf{\mu })\Vert }_2^2,n=1,\dots ,N_t,$$

$$\hat{\mathbf{x}}(t^n;\mathbf{\mu })=\underset{\hat{\mathbf{v}}\in {\mathbb{R}}^{n_s}}{\mathrm{argmin}}{\Vert \mathbf{A}{\mathbf{r}}^n({\mathbf{x}}^0(\mathbf{\mu })+\mathbf{\Phi }\hat{\mathbf{v}},\tilde{\mathbf{x}}(t^{n-1};\mathbf{\mu }),\dots ,\tilde{\mathbf{x}}(t^{n-k(t^n)};\mathbf{\mu });\mathbf{\mu })\Vert }_2^2,n=1,\dots ,N_t,$$

$$\tilde{\mathbf{y}}(\cdot ;\mathbf{\mu })\in \mathcal{ST}\subseteq {\mathbb{R}}^{N_s}\otimes \mathcal{H}$$

$$\begin{array}{rl}\tilde{\mathbf{y}}(t^n;\mathbf{\mu })& ={\mathbf{x}}^0(\mathbf{\mu })+\sum _{i=1}^{n_{st}}{\text{uppi}}_i(t^n){\hat{y}}_i(\mathbf{\mu })\end{array}$$

$$\mathbf{g}(\tilde{\mathbf{y}}(\cdot ;\mathbf{\mu }))\in \mathcal{S}\mathcal{T}\subseteq {\mathbb{R}}^{N_s}\otimes {\mathbb{R}}^{N_t}.$$

$$\begin{array}{rl}\mathbf{r}& :({\tau }^n;\mathbf{w};\mathbf{\nu })\mapsto {\mathbf{r}}^n(\mathbf{w}({\tau }^n),\dots ,\mathbf{w}({\tau }^{n-k({\tau }^n)});\mathbf{\nu })\\ & :\{t^n{\}}_{n=0}^{N_t}\times {\mathbb{R}}^{N_s}\otimes \mathcal{H}\times \mathcal{D}\to {\mathbb{R}}^{N_s}\end{array}$$

$$\begin{array}{rl}\mathsf{r}& :({\tau }^n;\hat{\mathbf{w}};\mathbf{\nu })\mapsto \mathbf{r}({\tau }^n;{\mathbf{x}}^0(\mathbf{\mu })+\sum _{i=1}^{n_{st}}{\text{uppi}}_i(\cdot ){\hat{w}}_i;\mathbf{\nu })\\ & :\{t^n{\}}_{n=0}^{N_t}\times {\mathbb{R}}^{n_{st}}\times \mathcal{D}\to {\mathbb{R}}^{N_s}.\end{array}$$

$$\begin{array}{rl}\mathbf{h}& :\mathbf{u}\mapsto \text{vec}(\mathbf{g}(\mathbf{u}))\\ & :{\mathbb{R}}^p\otimes \mathcal{H}\to {\mathbb{R}}^{p{N}_t},\end{array}$$

$$\begin{array}{rl}\bar{\mathbf{r}}& :(\mathbf{w};\mathbf{\mu })\mapsto \mathbf{h}(\mathbf{r}(\cdot ;\mathbf{w};\mathbf{\mu }))\\ & :{\mathbb{R}}^{N_s}\otimes \mathcal{H}\times \mathcal{D}\to {\mathbb{R}}^{N_s{N}_t}\end{array}$$

$$\begin{array}{rl}\bar{\mathsf{r}}& :(\hat{\mathbf{w}};\mathbf{\nu })\mapsto \mathbf{h}(\mathsf{r}(\cdot ;\hat{\mathbf{w}};\mathbf{\nu }))\\ & :{\mathbb{R}}^{n_{st}}\times \mathcal{D}\to {\mathbb{R}}^{N_s{N}_t}.\end{array}$$

$$\begin{array}{r}\tilde{\mathbf{y}}(\cdot ;\mathbf{\mu })=\underset{\mathbf{v}\in \mathcal{S}\mathcal{T}}{\mathrm{argmin}}{\Vert \mathbf{r}(\cdot ;\mathbf{v};\mathbf{\mu })\Vert }_{{\bar{\mathbf{A}}}^T\bar{\mathbf{A}}}^2=\underset{\mathbf{v}\in \mathcal{S}\mathcal{T}}{\mathrm{argmin}}{\Vert \bar{\mathbf{A}}\bar{\mathbf{r}}(\mathbf{v};\mathbf{\mu })\Vert }_2^2,\end{array}$$

$$\begin{array}{rl}\hat{\mathbf{y}}(\mathbf{\mu })& =\underset{\hat{\mathbf{v}}\in {\mathbb{R}}^{n_{st}}}{\mathrm{argmin}}{\Vert \mathsf{r}(\cdot ;\hat{\mathbf{v}};\mathbf{\mu })\Vert }_{{\bar{\mathbf{A}}}^T\bar{\mathbf{A}}}^2=\underset{\hat{\mathbf{v}}\in {\mathbb{R}}^{n_{st}}}{\mathrm{argmin}}{\Vert \bar{\mathbf{A}}\bar{\mathsf{r}}(\hat{\mathbf{v}};\mathbf{\mu })\Vert }_2^2,\end{array}$$

$$(\mathsf{r}(\cdot ;\hat{\mathbf{y}}(\mathbf{\mu });\mathbf{\mu }),\frac{\partial \mathsf{r}}{\partial {\hat{w}}_i}(\cdot ;\hat{\mathbf{y}}(\mathbf{\mu });\mathbf{\mu }){)}_{{\bar{\mathbf{A}}}^T\bar{\mathbf{A}}}=0,i\in \mathrm{Nature}{n}_{st},$$

$$\displaystyle\begin{split}\frac{\partial\boldsymbol{{\mathsf{r}}}}{\partial\hat{w}_{i}}(t^{n};\hat{\boldsymbol{y}};\boldsymbol{\mu})&=\frac{\partial\boldsymbol{r}}{\partial\boldsymbol{w}}\Bigl{(}t^{n};\boldsymbol{x}^{0}(\boldsymbol{\mu})+\sum_{\ell=1}^{{n_{st}}}\boldsymbol{\uppi}_{\ell}(\cdot)\hat{y}_{\ell};\boldsymbol{\mu}\Bigr{)}\boldsymbol{\uppi}_{i}(\cdot)\\ &=\sum_{j=0}^{k(t^{n})}\frac{\partial\boldsymbol{r}^{n}}{\partial\boldsymbol{w}^{j}}\Bigl{(}\boldsymbol{x}^{0}(\boldsymbol{\mu})+\sum_{\ell=1}^{{n_{st}}}\boldsymbol{\uppi}_{\ell}(t^{n})\hat{y}_{\ell}(\boldsymbol{\mu}),\ldots,\boldsymbol{x}^{0}(\boldsymbol{\mu})+\sum_{\ell=1}^{{n_{st}}}\boldsymbol{\uppi}_{\ell}(t^{n-k(t^{n})})\hat{y}_{\ell}(\boldsymbol{\mu});\boldsymbol{\mu}\Bigr{)}\boldsymbol{\uppi}_{i}(t^{n-j})\\ &=\sum_{j=0}^{k(t^{n})}\left[\alpha_{j}^{n}-\Delta t^{n}\beta_{j}^{n}\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{w}}\Bigl{(}\boldsymbol{x}^{0}(\boldsymbol{\mu})+\sum_{\ell=1}^{{n_{st}}}\boldsymbol{\uppi}_{\ell}(t^{n-j})\hat{y}_{\ell}(\boldsymbol{\mu}),t^{n-j};\boldsymbol{\mu}\Bigr{)}\right]\boldsymbol{\uppi}_{i}(t^{n-j})\end{split}$$

$${\left(\frac{\partial \mathsf{r}}{\partial \hat{\mathbf{w}}}(\cdot ;{\hat{\mathbf{y}}}^{(k)};\mathbf{\mu })\delta {\hat{\mathbf{y}}}^{(k)},\frac{\partial \mathsf{r}}{\partial \hat{\mathbf{w}}}(\cdot ;{\hat{\mathbf{y}}}^{(k)};\mathbf{\mu })\right)}_{{\bar{\mathbf{A}}}^T\bar{\mathbf{A}}}$$

$${\hat{\mathbf{y}}}^{(k+1)}$$

$$\tilde{\mathbf{y}}(\cdot ;\mathbf{\mu })=\underset{\mathbf{w}\in \mathcal{S}\mathcal{T}}{\mathrm{argmin}}\sum _{n=1}^{N_t}{\Vert {\bar{\mathbf{A}}}^n{\mathbf{r}}^n(\mathbf{w};\mathbf{\mu })\Vert }_2^2\text{and}\hat{\mathbf{y}}(\mathbf{\mu })=\underset{\hat{\mathbf{v}}\in {\mathbb{R}}^{n_{st}}}{\mathrm{argmin}}\sum _{n=1}^{N_t}{\Vert {\bar{\mathbf{A}}}^n\mathsf{r}(t^n;\hat{\mathbf{v}};\mathbf{\mu })\Vert }_2^2,$$

$$(\mathsf{r}(\cdot ;{\hat{\mathbf{y}}}_{\text{G}};\mathbf{\mu }),{\text{uppi}}_i(\cdot ){)}_{\mathbf{\Theta }}=0,i\in \mathrm{Nature}{n}_{st}$$

$$\bar{\mathbf{Z}}:={\left[{\mathbf{e}}_{\phi (s_1,t_1)}\cdots {\mathbf{e}}_{\phi (s_{{\bar{n}}_z},t_{{\bar{n}}_z})}\right]}^T\in \{0,1{\}}^{{\bar{n}}_z\times N_s{N}_t},$$

$$\tilde{\mathbf{y}}(\cdot ;\mathbf{\mu })=\underset{\mathbf{v}\in \mathcal{S}\mathcal{T}}{\mathrm{argmin}}{\Vert \tilde{\mathbf{r}}(\cdot ;\mathbf{v};\mathbf{\mu })\Vert }_2^2\text{and}\hat{\mathbf{y}}(\mathbf{\mu })=\underset{\hat{\mathbf{v}}\in {\mathbb{R}}^{n_{st}}}{\mathrm{argmin}}{\Vert \tilde{\mathsf{r}}(\cdot ;\hat{\mathbf{v}};\mathbf{\mu })\Vert }_2^2,$$

$$\begin{array}{rl}\mathbf{h}(\tilde{\mathbf{r}}(\cdot ;\mathbf{v};\mathbf{\mu }))& =\underset{\bar{\mathbf{v}}\in \text{range}({\bar{\mathbf{\Phi }}}_r)}{\mathrm{argmin}}\Vert \bar{\mathbf{Z}}\bar{\mathbf{v}}-\bar{\mathbf{Z}}\bar{\mathbf{r}}(\mathbf{v};\mathbf{\mu }){\Vert }_2^2={\bar{\mathbf{\Phi }}}_r(\bar{\mathbf{Z}}{\bar{\mathbf{\Phi }}}_r{)}^{+}\bar{\mathbf{Z}}\bar{\mathbf{r}}(\mathbf{v};\mathbf{\mu })\\ \mathbf{h}(\tilde{\mathsf{r}}(\cdot ;\hat{\mathbf{v}};\mathbf{\mu }))& =\underset{\bar{\mathbf{v}}\in \text{range}({\bar{\mathbf{\Phi }}}_r)}{\mathrm{argmin}}\Vert \bar{\mathbf{Z}}\bar{\mathbf{v}}-\bar{\mathbf{Z}}\bar{\mathsf{r}}(\hat{\mathbf{v}};\mathbf{\mu }){\Vert }_2^2={\bar{\mathbf{\Phi }}}_r(\bar{\mathbf{Z}}{\bar{\mathbf{\Phi }}}_r{)}^{+}\bar{\mathbf{Z}}\bar{\mathsf{r}}(\hat{\mathbf{v}};\mathbf{\mu }).\end{array}$$

$$\mathcal{ST}={\mathbf{x}}^0(\mathbf{\mu })\otimes \mathcal{O}+\mathrm{span}\{\mathbf{x}(\cdot ;{\mathbf{\mu }}_{\text{train}}^i)-{\mathbf{x}}^0({\mathbf{\mu }}_{\text{train}}^i){\}}_{i=1}^{n_{\text{train}}}\subseteq {\mathbb{R}}^{N_s}\otimes \mathcal{H}$$

Full text

1 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 [27]; molecular dynamics simulations [20] and condensed phase dynamics [42] 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 [3, 25]) and simultaneous analysis and design (SAND) in PDE-constrained optimization [26, 36, 15, 16, 37, 39, 4] 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 [38, 28, 34, 11], 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 [29, 32, 43]. 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 [31], PITA [23], and MGRIT [22]) 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) [21]. These forecasts can be used (1) to generate accurate initial guesses for the Newton solver at each time step [12], or (2) as an accurate coarse propagator to accelerate convergence of time-parallel methods [10]. 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 [45, 46, 50, 49], POD–Galerkin [48, 6], and ODE-residual minimization [17] 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 [45, 46, 50, 49] 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. [17]) 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. [6]) 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.

The original spatial-projection-based LSPG method [9, 11, 8] performed the following steps: (1) apply temporal discretization111When explicit time integration is used, LSPG projection is equivalent to Galerkin projection [8]. to the system of ODEs characterizing the high-fidelity model, (2) introduce a low-dimensional spatial trial subspace, and (3) compute the solution in the spatial subspace that minimizes (in a weighted $\ell^{2}$-norm) the discrete residual arising at each time step. This approach does not reduce temporal dimensionality, as the ROM and high-fidelity model time steps are typically the same. Instead, the proposed ST-LSPG method executes the following steps: (1) apply time integration to the system of nonlinear ODEs with an implicit or explicit linear multistep method, (2) introduce a low-dimensional space–time trial subspace, and (3) compute the solution in the space–time trial subspace that minimizes the discrete residual over all space and time in a weighted $\ell^{2}$-norm. This norm can be selected to enable hyper-reduction based on both collocation and gappy POD; this is in analogy to collocation [30, 1, 41] and gappy POD [1, 11] methods applied to spatial-projection-based ROMs.

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 priori error bounds that enable the error in the ST-LSPG solution to be bounded by the best approximation error as measured in the $\ell^{2}$-norm over all time (Theorem 3), and the $\ell^{\infty}$-norm over all time (Theorem 4). Critically, the stability constants for these bounds grow linearly and subquadratically in time, respectively (Remark 6.1).

• •

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. [17], 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. [17] did not provide error bounds for the residual-minimizing approximation (it provides an error bound only for the best linear-subspace approximation). Also, Ref. [17] 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. [17] 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. [17] 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.

2 Full-order model

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

2.1 Time-continuous representation

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

[TABLE]

where $t\in[0,T]$ denotes time with $T\in\mathbb{R}_{+}$ denoting the final time, and $\boldsymbol{x}^{\star}(t;\boldsymbol{\mu})$ denotes the time-dependent, parameterized state implicitly defined as the solution to problem (2.1) with $\boldsymbol{x}^{\star}:[0,T]\times\mathcal{D}\rightarrow\mathbb{R}^{N_{s}}$ and $\boldsymbol{x}^{\star}(\cdot;\boldsymbol{\mu})\in\mathbb{R}^{N_{s}}\otimes H$. Here, $H$ denotes the set of sufficiently smooth functions from $[0,T]$ to $\mathbb{R}$ (e.g., $H=H^{1}([0,T])$) under consideration and $\otimes$ denotes the tensor product. Further, $\boldsymbol{f}:\mathbb{R}^{N_{s}}\times[0,T]\times\mathcal{D}\rightarrow\mathbb{R}^{N_{s}}$ with $(\boldsymbol{w},\tau;\boldsymbol{\nu})\mapsto\boldsymbol{f}(\boldsymbol{w},\tau;\boldsymbol{\nu})$ denotes the velocity, which we assume to be nonlinear in at least its first argument, $\boldsymbol{x}^{0}:\mathcal{D}\rightarrow\mathbb{R}^{N_{s}}$ denotes the initial state, and $\boldsymbol{\mu}\in\mathcal{D}$ denotes the parameters with parameter domain $\mathcal{D}\subseteq\mathbb{R}^{n_{\mu}}$.

2.2 Time-discrete representation: linear multistep methods

We now introduce a (generally nonuniform) time discretization characterized by time step $\Delta t^{n}\in\mathbb{R}_{+}$ and time instances $t^{n}=t^{n-1}+\Delta t^{n}$, $n\in{\rm Nature}{{N_{t}}}$ with $t^{0}=0$, ${N_{t}}\in\mathbb{N}$, and ${\rm Nature}{N}:=\{1,\ldots,N\}$.222Note that by assuming a fixed time discretization, we do not allow for adaptive time-step selection; enabling the proposed ST-LSPG method to be applied in the context of adaptive time stepping is the subject of future research. Applying a linear $k$-step method333We consider only linear multistep methods for simplicity of presentation. One could develop the proposed ST-LSPG method for alternative time integrators in a rather straightforward manner; see, e.g., Ref. [8], which develops LSPG models for explicit, fully implicit, and diagonally implicit Runge–Kutta schemes. to numerically solve Eq. (2.1) using this discretization yields an O$\Delta$E, which is characterized by the following system of nonlinear algebraic equations to be solved for the numerical solution $\boldsymbol{x}(t^{n};\boldsymbol{\mu})\in\mathbb{R}^{N_{s}}$ at each time instance:

[TABLE]

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

[TABLE]

where coefficients $\alpha_{j}^{n},\beta_{j}^{n}\in\mathbb{R}$, $j=0,\ldots,k(t^{n})$ define a particular linear multistep scheme, $\alpha_{0}^{n}\neq 0$, and $\sum_{j=0}^{k(t^{n})}\alpha_{j}^{n}=0$ is necessary for consistency. We note that the numerical solution satisfies

[TABLE]

with $\mathscr{H}$ the set of functions from $\{t^{n}\}_{n=0}^{{N_{t}}}$ to $\mathbb{R}$. We refer to $\mathbb{R}^{N_{s}}\otimes\mathscr{H}$ as the ‘FOM trial space’ in which solutions are sought. An isomorphism exists between $\mathscr{H}$ and $\mathbb{R}^{{N_{t}}}$ provided by the (invertible) function ${\boldsymbol{g}}$ that ‘unrolls’ time according to the time discretization as

[TABLE]

where $p\in\mathbb{N}$ is an arbitrary dimension. This yields an equivalent (discrete) representation of the FOM space as

[TABLE]

3 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 [9, 11, 8].

3.1 Spatial trial subspace

The original LSPG method applies spatial projection using a subspace $\mathcal{S}:=\mathrm{span}\{\boldsymbol{\phi}_{i}\}_{i=1}^{n_{s}}\subseteq\mathbb{R}^{N_{s}}$ with $\dim(\mathcal{S})={n_{s}}\leq N_{s}$ (hopefully with ${n_{s}}\ll N_{s}$). Using this subspace, the LSPG method approximates the numerical solution at each time instance $t^{n}$, $n\in{\rm Nature}{{N_{t}}}$ as $\boldsymbol{x}(t^{n};\boldsymbol{\mu})\approx\tilde{\boldsymbol{x}}(t^{n};\boldsymbol{\mu})\in\boldsymbol{x}^{0}(\boldsymbol{\mu})+\mathcal{S}$ or equivalently

[TABLE]

where $\hat{x}_{i}(\cdot;\boldsymbol{\mu})\in\mathbb{R}\otimes\mathscr{H}$ with $\hat{x}_{i}(0;\boldsymbol{\mu})=0$, $i\in{\rm Nature}{{n_{s}}}$ denotes the generalized coordinates. This approach is equivalent to enforcing the approximated numerical solution $\tilde{\boldsymbol{x}}(t;\boldsymbol{\mu})$ with $\tilde{\boldsymbol{x}}:[0,T]\times\mathcal{D}\rightarrow\mathbb{R}^{N_{s}}$ to reside in an affine ‘spatial trial subspace’

[TABLE]

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

[TABLE]

Remark 3.1** (LSPG projection does not reduce the temporal dimension).**

Comparing (2.6) and (3.3) reveals that LSPG projection reduces spatiotemporal dimension of the problem from $\dim(\mathbb{R}^{N_{s}}\otimes\mathbb{R}^{{N_{t}}})=N_{s}{N_{t}}$ to $\dim(\boldsymbol{x}^{0}(\boldsymbol{\mu})\otimes\boldsymbol{1}_{{N_{t}}}+\mathcal{S}\otimes\mathbb{R}^{{N_{t}}})={n_{s}}{N_{t}}$. Thus, while the spatial dimension has been reduced from $N_{s}$ to ${n_{s}}$, the temporal dimension ${N_{t}}$ has not been reduced. This can preclude significant computational-cost savings when the original problem is characterized by a long time interval $T$, or if small time steps $\Delta t^{n}$ are needed for accuracy (e.g., for stiff dynamical systems).

3.2 Spatial projection

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

[TABLE]

or equivalently

[TABLE]

where $\hat{\boldsymbol{x}}(\cdot;\boldsymbol{\mu})\equiv[\hat{x}_{1}(\cdot;\boldsymbol{\mu})\ \cdots\ \hat{x}_{{n_{s}}}(\cdot;\boldsymbol{\mu})]^{T}\in\mathbb{R}^{{n_{s}}}\otimes\mathscr{H}$. Here, $\boldsymbol{A}\in\mathbb{R}^{z\times N_{s}}$ is a weighting matrix, where ${n_{s}}\leq z(\leq N_{s})$ is necessary for the residual Jacobian in the nonlinear least-squares problem (3.4)–(3.5) to be nonsingular. Examples of weighting matrices include $\boldsymbol{A}=\boldsymbol{I}_{N_{s}}$ in the case of unweighted LSPG, where $\boldsymbol{I}_{p}$ denotes the $p\times p$ identity matrix; however, this choice precludes computational-cost savings, as all $N_{s}$ elements of the residual vector $\boldsymbol{r}^{n}$ (as well as its Jacobian) must be computed during each iteration when solving nonlinear least-squares problem (3.4)–(3.5). To reduce the computational complexity, a ‘hyper-reduction’ approach is required that necessitates computing only a subset of residual elements. Particular weighting-matrix choices that lead to hyper-reduction include $\boldsymbol{A}=\boldsymbol{Z}\in\{0,1\}^{n_{z}\times N_{s}}$, where $n_{z}\leq N_{s}$ and $\boldsymbol{Z}$ comprises selected rows of the identity matrix $\boldsymbol{I}_{N_{s}}$ in the case of collocation [30]; and $\boldsymbol{A}=(\boldsymbol{Z}\boldsymbol{\Phi}_{r})^{+}\boldsymbol{Z}\in\mathbb{R}^{n_{r}\times N_{s}}$, where $\boldsymbol{\Phi}_{r}\in\mathbb{V}_{{n_{r}}}(\mathbb{R}^{{N_{s}}})$ denotes a basis for the residual and a superscript $+$ denotes the Moore–Penrose pseudoinverse in the case of GNAT [11]. Here, ${n_{s}}\leq n_{r}\leq n_{z}(\leq N_{s})$ is necessary for the residual Jacobian to be nonsingular, and $\mathbb{V}_{{k}}(\mathbb{R}^{{n}})$ denotes the Stiefel manifold: the set of orthogonal $k$-frames in $\mathbb{R}^{n}$.

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

4.1 Space–time trial subspace

To reduce both the spatial and temporal dimensions of the FOM, in analogy to Eq. (3.2) we enforce the approximated numerical solution $\tilde{\boldsymbol{y}}$ to reside in an affine ‘space–time trial subspace’

[TABLE]

where $\mathscr{ST}:=\boldsymbol{x}^{0}(\boldsymbol{\mu})\otimes\mathscr{O}+\mathrm{span}\{\boldsymbol{\uppi}_{i}\}_{i=1}^{n_{st}}\subseteq\mathbb{R}^{N_{s}}\otimes\mathscr{H}$ with $\dim(\mathscr{ST})={n_{st}}\ll N_{s}{N_{t}}$, and we enforce $\boldsymbol{\uppi}_{i}(0)={\bf 0}$, $i\in{\rm Nature}{{n_{st}}}$ such that $\tilde{\boldsymbol{y}}(0;\boldsymbol{\mu})=\boldsymbol{x}^{0}(\boldsymbol{\mu})$.

Thus, at a given time instance, the space–time LSPG method approximates the numerical solution at each time instance $n\in{\rm Nature}{{N_{t}}}$ as $\boldsymbol{x}(t^{n};\boldsymbol{\mu})\approx\tilde{\boldsymbol{y}}(t^{n};\boldsymbol{\mu})\in\mathscr{ST}$ or equivalently

[TABLE]

where $\hat{y}_{i}(\boldsymbol{\mu})\in\mathbb{R}$, $i\in{\rm Nature}{{n_{st}}}$ denotes the generalized coordinate of the ST-LSPG solution. Comparing Eqs. (3.1) and (4.2) reveals that the space–time trial subspace enables time dependence of the approximated solution to be moved from the generalized coordinates to the basis vectors; this enables fewer generalized coordinates to be computed in order to characterize the complete space–time solution.

Introducing $\boldsymbol{\pi}_{i}:={{\boldsymbol{g}}}(\boldsymbol{\uppi}_{i})\in\mathbb{R}^{N_{s}}\otimes\mathbb{R}^{{N_{t}}}$ and $\mathcal{ST}:=\boldsymbol{x}^{0}(\boldsymbol{\mu})\otimes\boldsymbol{1}_{{N_{t}}}+\mathrm{span}\{\boldsymbol{\pi}_{i}\}_{i=1}^{{n_{st}}}$, we can also write the space–time trial subspace in discrete form as

[TABLE]

Remark 4.1** (Space–time LSPG projection reduces the temporal dimension).**

Comparing (2.6) and (4.3) reveals that the proposed space–time LSPG trial subspace reduces the spatiotemporal dimension of the problem from $\dim(\mathbb{R}^{N_{s}}\otimes\mathbb{R}^{{N_{t}}})=N_{s}{N_{t}}$ to $\dim(\mathcal{ST})={n_{st}}$. Because the spatiotemporal dimension ${n_{st}}$ can be chosen to be independent of both the spatial and temporal dimensions $N_{s}$ and ${N_{t}}$, respectively, the proposed method can reduce both the spatial and temporal dimensions of the full-order model.

Remark 4.2** **(Space–time trial subspace can remove spurious temporal

modes).

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

4.2 Space–time least-squares Petrov–Galerkin

projection

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

[TABLE]

Note that the space–time residuals $\boldsymbol{r}(\cdot;\boldsymbol{w};\boldsymbol{\nu}),\boldsymbol{{\mathsf{r}}}(\cdot;\hat{\boldsymbol{w}};\boldsymbol{\nu})\in\mathbb{R}^{N_{s}}\otimes\mathscr{H}$ are defined from the O$\Delta$E residual $\boldsymbol{r}^{n}$; as such, they are defined directly from the integrator used to perform time discretization for the FOM. We now introduce a vectorization function

[TABLE]

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

[TABLE]

and as a function of the generalized coordinates

[TABLE]

We define the inner product $(\boldsymbol{u},\boldsymbol{v})_{\boldsymbol{\Theta}}:={\boldsymbol{h}}(\boldsymbol{v})^{T}{\boldsymbol{\Theta}}{\boldsymbol{h}}(\boldsymbol{u})$ and associated norm $\|\boldsymbol{u}\|_{\boldsymbol{\Theta}}=\sqrt{(\boldsymbol{u},\boldsymbol{u})_{\boldsymbol{\Theta}}}$ for $\boldsymbol{u},\boldsymbol{v}\in\mathbb{R}^{p}\otimes\mathscr{H}$ and ${\boldsymbol{\Theta}}\in\text{SPSD}(p{N_{t}})$, where $\text{SPSD}(p)$ denotes the set of $p\times p$ symmetric positive semidefinite matrices; we also define the inner product $(\boldsymbol{u},\boldsymbol{v})_{2}:={\boldsymbol{h}}(\boldsymbol{v})^{T}{\boldsymbol{h}}(\boldsymbol{u})$ and the associated norm $\|\boldsymbol{u}\|_{2}=\sqrt{(\boldsymbol{u},\boldsymbol{u})_{2}}$.444The metric will be rank deficient if we employ hyper-reduction as described in Section 4.3. Now, we propose computing the ST-LSPG solution by minimizing the residual in a weighted $\ell^{2}$-norm as

[TABLE]

where $\bar{\boldsymbol{A}}\in\mathbb{R}^{\bar{z}\times N_{s}{N_{t}}}$ is a space–time weighting matrix and ${n_{st}}\leq\bar{z}(\leq N_{s}{N_{t}})$ is necessary for the residual Jacobian in the nonlinear least-squares problem (4.9) to be nonsingular. We can also write Problem (4.9) in terms of the generalized coordinates as

[TABLE]

where $\hat{\boldsymbol{y}}\equiv[\hat{y}_{1}(\boldsymbol{\mu})\ \cdots\ \hat{y}_{{n_{st}}}(\boldsymbol{\mu})]^{T}\in\mathbb{R}^{{n_{st}}}$ and Eq. (4.2) relates $\hat{\boldsymbol{y}}$ to $\tilde{\boldsymbol{y}}$.

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

[TABLE]

where

[TABLE]

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)}$:

[TABLE]

where we set $\hat{\boldsymbol{y}}(\boldsymbol{\mu})=\hat{\boldsymbol{y}}^{(k_{\text{max}}(\boldsymbol{\mu}))}$ at convergence and $\alpha^{(k)}\in\mathbb{R}$ denotes a step length that can be computed to ensure global convergence (e.g., satisfy the strong Wolfe conditions). Again, we note that the condition ${n_{st}}\leq\bar{z}$ is necessary for singular values of the Jacobian $\bar{\boldsymbol{A}}\frac{\partial\boldsymbol{{\mathsf{r}}}}{\partial\hat{\boldsymbol{w}}}(\cdot;\hat{\boldsymbol{y}}^{(k)};\boldsymbol{\mu})\in\mathbb{R}^{\bar{z}\times{n_{st}}}$ to be uniformly bounded away from zero in the region of interest; this is one of the sufficient conditions required to prove convergence of the Gauss–Newton method (see, e.g., [35, Theorem 10.1]).

Remark 4.3** (Simplification for block-diagonal weighting matrices).**

If the weighting matrix is block-diagonal, i.e., $\bar{\boldsymbol{A}}=\mathop{\operator@font diag}\nolimits{\bar{\boldsymbol{A}}^{n}}$ with $\bar{\boldsymbol{A}}^{n}\in\mathbb{R}^{\bar{z}^{n}\times N_{s}}$, $\bar{z}^{n}\leq N_{s}$, and $\bar{z}=\sum_{n=1}^{N_{t}}\bar{z}^{n}$, then Problems (4.9) and (4.10) simplify to

[TABLE]

*respectively. *

Remark 4.4** (Space–time Galerkin projection).**

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

[TABLE]

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. [48, 6]; 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.

4.3 Weighting matrix and hyper-reduction

Hyper-reduction refers to reducing the computational complexity of nonlinear ROMs by introducing approximations of the nonlinear functions. This is typically achieved via collocation (wherein the nonlinear functions are simply sampled) [1, 41, 30] or function-reconstruction approaches (e.g., gappy POD [21], empirical interpolation [5, 13]), wherein the nonlinear function is approximated from a sample of its entries via interpolation or least-squares regression [1, 7, 13, 24, 19, 9, 11]. In the case of (spatial-projection-based) LSPG, it has been shown that hyper-reduction can be realized by particular choices of the weighting matrix [8]. We now propose several choices for the space–time weighting matrix $\bar{\boldsymbol{A}}$, some of which will lead to hyper-reduction for the ST-LSPG method.

All of these hyper-reduction methods will ensure that the computational cost of the ST-LSPG method is independent of both the spatial and temporal dimensions characterizing the FOM if the Jacobian of the space–time residual is sparse, i.e., the number of nonzeros in ${\boldsymbol{h}}(\frac{\partial r_{i}}{\partial\boldsymbol{w}}(t^{n};{\boldsymbol{v}};\boldsymbol{\mu}))\in\mathbb{R}^{N_{s}{N_{t}}}$, $i\in{\rm Nature}{N_{s}}$, $n\in{\rm Nature}{{N_{t}}}$ is ‘small’ and is independent of both the spatial and temporal dimensions $N_{s}$ and ${N_{t}}$, respectively. This is sometimes referred to as the $H$-independence condition [19] and is inherited directly from $H$-independence of the spatial residual $\boldsymbol{r}^{n}$, as $\mathrm{nnz}({\boldsymbol{h}}(\frac{\partial r_{i}}{\partial\boldsymbol{w}}(t^{n};{\boldsymbol{v}};\boldsymbol{\mu})))=\sum_{j=0}^{k(t^{n})}\mathrm{nnz}(\frac{\partial r_{i}^{n}}{\partial\boldsymbol{w}^{j}}({\boldsymbol{v}}^{n},\ldots,{\boldsymbol{v}}^{n-k(t^{n})};\boldsymbol{\mu}))$, where $\mathrm{nnz}$ denotes the number of nonzeros of its argument.

4.3.1 Unweighted LSPG

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. [6], 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.

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

[TABLE]

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.

4.3.3 Space–time GNAT

Similarly, we can extend gappy POD hyper-reduction to ST-LSPG by employing a weighting matrix $\bar{\boldsymbol{A}}=(\bar{\boldsymbol{Z}}{\bar{\boldsymbol{\Phi}}}_{r})^{+}\bar{\boldsymbol{Z}}\in\mathbb{R}^{\bar{n}_{r}\times N_{s}{N_{t}}}$, where ${\bar{\boldsymbol{\Phi}}}_{r}\in\mathbb{V}_{{\bar{n}_{r}}}(\mathbb{R}^{{N_{s}{N_{t}}}})$ denotes a basis (in matrix form) for the spatiotemporal residual. Again, we require ${n_{st}}\leq\bar{n}_{r}\leq\bar{n}_{z}(\leq N_{s}{N_{t}})$ for nonsingular residual Jacobians. When this weighting matrix is employed, Problem (4.9)–(4.10) is equivalent to minimizing the $\ell^{2}$-norm of the gappy POD-approximated residual, i.e., Problem (4.9)–(4.10) with $\bar{\boldsymbol{A}}=(\bar{\boldsymbol{Z}}{\bar{\boldsymbol{\Phi}}}_{r})^{+}\bar{\boldsymbol{Z}}\in\mathbb{R}^{\bar{n}_{r}\times N_{s}{N_{t}}}$ is equivalent to

[TABLE]

where we have defined the gappy POD residual approximations as

[TABLE]

As with space–time collocation, space–time gappy POD leads to hyper-reduction, as evaluating the objective function requires evaluating only $\bar{n}_{z}\leq N_{s}{N_{t}}$ elements of the spatiotemporal residual. It 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}_{r}\leq N_{s}{N_{t}}$ in this case. Due to its close relationship to the original GNAT method—in which case $\boldsymbol{A}=(\boldsymbol{Z}\boldsymbol{\Phi}_{r})^{+}\boldsymbol{Z}\in\mathbb{R}^{n_{r}\times N_{s}}$ in Problem (3.4)–(3.5), we refer to this approach as the space–time GNAT (ST-GNAT) method.

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

5.1 Space–time trial subspace

We first assume that a set of training parameter instances $\mathcal{D}_{\text{train}}:=\{\boldsymbol{\mu}^{1}_{\text{train}},\ldots,\boldsymbol{\mu}^{n_{\text{train}}}_{\text{train}}\}\subset\mathcal{D}$ has been defined (by, e.g., uniform sampling, Latin-hypercube sampling, greedy sampling) for which the FOM (2.1) has been solved numerically using a linear multistep method (2.2) to obtain state ‘snapshots’ $\boldsymbol{x}(\cdot;\boldsymbol{\mu}^{i}_{\text{train}})\in\mathbb{R}^{N_{s}}$, $i\in{\rm Nature}{n_{\text{train}}}$. This data-collection process is referred to as the ‘offline’ stage in model reduction.

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

[TABLE]

such that $\boldsymbol{\uppi}_{i}=\boldsymbol{x}(\cdot;\boldsymbol{\mu}^{i}_{\text{train}})-\boldsymbol{x}^{0}(\boldsymbol{\mu}^{i}_{\text{train}})\in\mathbb{R}^{N_{s}}\otimes\mathscr{H}$, $i\in{\rm Nature}{n_{\text{train}}}$ [45, 46, 50, 49, 17]. Unfortunately, because this approach extracts only a single space–time basis vector from each training simulation, it limits the dimension of the space–time ROM to ${n_{st}}=n_{\text{train}}$. In practical contexts where a single training simulation may incur significant computational costs, this can significantly limit the dimensionality (and resulting accuracy) of the space–time ROM. Further, storage costs can be significant in this case, as the basis requires ${n_{st}}N_{s}{N_{t}}$ storage.

To overcome these shortcomings, we propose to compute the space–time trial subspace by applying tensor-decomposition techniques to the three-way ‘state tensor’ $\mathcal{X}\in\mathbb{R}^{N_{s}\times{N_{t}}\times n_{\text{train}}}$ with elements

[TABLE]

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

[TABLE]

where $\mathcal{S}_{i}:=\mathrm{span}\{\boldsymbol{\phi}_{i}\}\subseteq\mathbb{R}^{N_{s}}$, $\mathscr{T}_{i}:=\mathrm{span}\{\boldsymbol{\uppsi}_{j}^{i}\}_{j=1}^{n_{t}^{i}}\subseteq\mathscr{H}$, and $n_{t}^{i}\leq{N_{t}}$ (hopefully with $n_{t}^{i}\ll{N_{t}}$) for $i\in{\rm Nature}{{n_{s}}}$. Here, $\boldsymbol{\uppsi}_{j}^{i}\in\mathscr{H}$ with $\boldsymbol{\uppsi}_{j}^{i}(0)=0$ denotes the $j$th temporal basis vector associated with the temporal behavior of $i$th spatial basis vector $\boldsymbol{\phi}_{i}$. Thus, in this case we have $\boldsymbol{\uppi}_{\mathcal{I}(i,j)}=\boldsymbol{\phi}_{i}\boldsymbol{\uppsi}_{j}^{i}\in\mathbb{R}^{N_{s}}\otimes\mathscr{H}$, $i\in{\rm Nature}{{n_{s}}}$, $j\in{\rm Nature}{n_{t}^{i}}$ and ${n_{st}}=\sum_{i=1}^{{n_{s}}}n_{t}^{i}$, where $\mathcal{I}:(i,j)\mapsto\sum_{k=1}^{i-1}n_{t}^{k}+j$ provides a mapping from the spatial-basis and temporal-basis indices to a space–time basis index. This approach enables a larger space–time ROM dimension, as ${n_{st}}=\sum_{i=1}^{{n_{s}}}n_{t}^{i}>n_{\text{train}}$ is possible; further, this approach requires only ${n_{s}}N_{s}+{n_{st}}{N_{t}}$ storage at most. We expect the resulting ST-LSPG ROM be accurate if the solution exhibits separable behavior in space and time.

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

[TABLE]

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.

5.1.1 Spatial subspaces

We propose to compute the spatial subspaces $\mathcal{S}_{i}:=\mathrm{span}\{\boldsymbol{\phi}_{i}\}$, $i\in{\rm Nature}{{n_{s}}}$ via proper orthogonal decomposition (POD) applied to the spatial snapshot matrix $\boldsymbol{X}_{(1)}$. In particular, we compute the spatial bases from the singular value decomposition (SVD) as

[TABLE]

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

5.1.2 Fixed temporal

subspace via T-HOSVD

The most straightforward approach is to compute a fixed temporal subspace $\mathscr{T}:=\mathrm{span}\{\boldsymbol{\uppsi}_{i}\}_{i=1}^{n_{t}}\subset\mathscr{H}$ such that $\mathscr{T}_{i}=\mathscr{T}$ and $n_{t}^{i}=n_{t}$, $i\in{\rm Nature}{{n_{s}}}$. We can compute such a subspace by applying POD to the temporal snapshot matrix $\boldsymbol{X}_{(2)}$. Specifically, we can compute the fixed temporal bases from the SVD of $\boldsymbol{X}_{(2)}$ as

[TABLE]

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) [18, 44] 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. [6] in the context of space–time Galerkin projection performed at the time-continuous level.

5.1.3 Fixed temporal

subspace via ST-HOSVD

Alternatively, we recall from Eq. (4.1) that each space–time basis vector is the Kronecker product of a spatial basis vector with a temporal basis vector. Thus, it is sensible to compute the temporal subspace according to the observed time evolution of the solution in the coordinates defined by the spatial subspace. Mathematically, we can achieve this by applying the sequentially truncated HOSVD (ST-HOSVD) [47, 2]. Rather than computing the SVD of $\boldsymbol{X}_{(2)}$, which is agnostic to dimensionality reduction in space, ST-HOSVD instead computes the SVD of the mode-2 unfolding of $\mathcal{X}(\boldsymbol{\Phi})$, where we have defined $\mathcal{X}(\boldsymbol{V}):=\mathcal{X}\times_{1}\boldsymbol{V}$ for $\boldsymbol{V}\in\mathbb{R}^{N_{s}\times p}$ such that $\boldsymbol{X}(\boldsymbol{V})_{(1)}=\boldsymbol{V}^{T}\boldsymbol{X}_{(1)}$. In particular, we have

[TABLE]

where $n_{t}\leq{n_{s}}n_{\text{train}}$ and $\boldsymbol{U}_{t}(\boldsymbol{\Phi})\equiv\left[\boldsymbol{u}_{t}^{1}(\boldsymbol{\Phi})\cdots\ \boldsymbol{u}_{t}^{{n_{s}}n_{\text{train}}}(\boldsymbol{\Phi})\right]$. In addition to enabling the temporal basis to be associated with the time evolution of the spatial basis $\boldsymbol{\Phi}$, this approach is less computationally expensive than applying the T-HOSVD, as $\boldsymbol{X}(\boldsymbol{\Phi})_{(2)}\in\mathbb{R}^{{N_{t}}\times{n_{s}}n_{\text{train}}}$, while $\boldsymbol{X}_{(2)}\in\mathbb{R}^{{N_{t}}\times N_{s}n_{\text{train}}}$ and typically ${n_{s}}\ll N_{s}$. This temporal subspace also requires $n_{t}{N_{t}}$ storage.

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

[TABLE]

where $\boldsymbol{\psi}^{i}_{j}:={\boldsymbol{h}}(\boldsymbol{\uppsi}_{j}^{i})$ and $\boldsymbol{U}_{t}(\boldsymbol{\phi}_{i}):=\left[\boldsymbol{u}_{t}^{1}(\boldsymbol{\phi}_{i})\ \cdots\ \boldsymbol{u}_{t}^{n_{\text{train}}}(\boldsymbol{\phi}_{i})\right]$. This approach generates a tailored temporal subspace $\mathscr{T}_{i}:=\mathrm{span}\{\boldsymbol{\uppsi}_{j}^{i}\}_{j=1}^{n_{t}^{i}}$ for each spatial subspace $\mathcal{S}_{i}:=\mathrm{span}\{\boldsymbol{\phi}_{i}\}$. Further, because $\boldsymbol{X}(\boldsymbol{\phi}_{i})_{(2)}\in\mathbb{R}^{{N_{t}}\times n_{\text{train}}}$, $i\in{\rm Nature}{{n_{s}}}$, the cost of computing the ${n_{s}}$ SVDs (5.12) is significantly less than computing the SVDs in either (5.8) or (5.10); this results from the quadratic dependence of the SVD cost on the number of columns in the matrix. However, the maximum dimension of each temporal basis is limited to the number of training-parameter instances, i.e., $n_{t}^{i}\leq n_{\text{train}}$, $i\in{\rm Nature}{{n_{s}}}$. This temporal subspace requires $\sum_{i=1}^{n_{s}}n_{t}^{i}{N_{t}}={n_{st}}{N_{t}}$ storage, which is larger than that required by the fixed temporal subspaces. We note that this was the approach employed to construct temporal bases in our previous work based on forecasting [12, 10].

5.2 Space–time residual basis

We propose to construct the space–time residual basis ${\bar{\boldsymbol{\Phi}}}_{r}\in\mathbb{V}_{{\bar{n}_{r}}}(\mathbb{R}^{{N_{s}{N_{t}}}})$ from training data comprising the space–time residual computed at a set of pairs of reduced coordinates and parameter instances, i.e., $\{\hat{\boldsymbol{y}}_{\text{res}}^{i},\boldsymbol{\mu}^{i}_{\text{res}}\}_{i\in{\rm Nature}{n_{\text{res}}}}$. In this case, the space–time residual ‘snapshots’ can be expressed in a residual tensor $\mathcal{R}\in\mathbb{R}^{N_{s}\times{N_{t}}\times n_{\text{res}}}$ with entries

[TABLE]

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

[TABLE]

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

Projection of FOM training solutions. This approach employs

[TABLE]

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

[TABLE]

where $\boldsymbol{\Pi}:=\begin{bmatrix}\boldsymbol{\pi}_{1}&\cdots&\boldsymbol{\pi}_{{n_{st}}}\end{bmatrix}\in\mathbb{R}^{N_{s}{N_{t}}\times{n_{st}}}$. This approach does not require any additional training simulations; it simply requires $n_{\text{res}}$ evaluations of the space–time residual. 3. 3.

Random samples. In this approach, the training set $\{\hat{\boldsymbol{y}}_{\text{res}}^{i},\boldsymbol{\mu}^{i}_{\text{res}}\}_{i\in{\rm Nature}{n_{\text{res}}}}$ comprises random samples (e.g., via Latin hypercube sampling) from $\mathcal{D}_{\hat{\boldsymbol{y}}}\times\mathcal{D}$ where $\mathcal{D}_{\hat{\boldsymbol{y}}}\subseteq\mathbb{R}^{{n_{st}}}$. This approach also requires only $n_{\text{res}}$ evaluations of the space–time residual.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

We compute the orthogonal residual basis ${\bar{\boldsymbol{\Phi}}}_{r}\in\mathbb{V}_{{\bar{n}_{r}}}(\mathbb{R}^{{N_{s}{N_{t}}}})$ from the QR factorization of the nonorthogonalized basis $\boldsymbol{\Pi}_{r}\equiv\begin{bmatrix}\boldsymbol{\pi}_{r,1}&\cdots&\boldsymbol{\pi}_{r,\bar{n}_{r}}\end{bmatrix}\in\mathbb{R}^{N_{s}{N_{t}}\times\bar{n}_{r}}$ as

[TABLE]

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.

5.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. [9, 11] that allows for oversampling to enable least-squares regression via gappy POD. 2. 2.

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. 3. 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 [11]—which is tasked with computing spatial samples associated with $\mathfrak{s}$—can be employed for all sampled time instances $\mathfrak{t}$.

5.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 [12] 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}}}.$

6 Error analysis

For simplicity, this section omits parameter dependence of all operators and assumes a uniform time step, i.e., $\Delta t^{n}=\Delta t$, $n\in{\rm Nature}{{N_{t}}}$. Thus, the FOM solution $\boldsymbol{x}$ and the ST-LSPG ROM solution $\tilde{\boldsymbol{y}}$ satisfy

[TABLE]

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

1. A1

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

[TABLE] 2. A2

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

[TABLE]

where

[TABLE]

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

The space–time weighting matrix $\bar{\boldsymbol{A}}$ is defined such that the residual in the weighted space–time norm is uniformly bounded below by the $\ell^{2}$-norm of the residual over all elements of the space–time trial subspace, i.e., there exists $P>0$ such that

[TABLE]

Lemma 1**.**

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

[TABLE]

with Lipschitz constant

[TABLE]

Proof.

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

[TABLE]

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

[TABLE]

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

[TABLE]

∎

Lemma 2**.**

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

[TABLE]

with inverse Lipschitz constant

[TABLE]

Proof.

Applying the reverse triangle inequality and employing Assumption A2 yields

[TABLE]

which directly leads to the desired result. ∎

Theorem 3** (a priori error bound with respect to $\ell^{2}$-optimal solution).**

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

[TABLE]

Proof.

We begin by defining the $\ell^{2}$-optimal solution $\boldsymbol{x}_{2}$, which satisfies

[TABLE]

Then, we can exploit the optimality properties of $\boldsymbol{x}_{2}$ and $\tilde{\boldsymbol{y}}$ from Eqs. (6.11) and (6.1), respectively; Lipschitz continuity (Lemma 1); residual-norm equivalence (Assumption A3); and inverse Lipschitz continuity of the residual (Lemma 2) to derive the following inequalities:

[TABLE]

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

Theorem 4** (a priori error bound with respect to $\ell^{\infty}$-optimal solution).**

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

[TABLE]

Proof.

We begin by defining the $\ell^{\infty}$-optimal solution $\boldsymbol{x}_{\infty}$, which satisfies

[TABLE]

where we have defined the $\ell^{\infty}$-norm as $\|\boldsymbol{w}\|_{\infty}:=\max_{n\in{\rm Nature}{{N_{t}}}}\|\boldsymbol{w}(t^{n})\|_{2}$. Then, we can exploit norm equivalence $\|\boldsymbol{w}\|_{\infty}\leq\|\boldsymbol{w}\|_{2}\leq\sqrt{n}\|\boldsymbol{w}\|_{\infty}$ for $\boldsymbol{w}\in\mathbb{R}^{n}$; Lipschitz continuity (Lemma 1); residual-norm equivalence (Assumption A3); and inverse Lipschitz continuity of the residual (Lemma 2); and the optimality properties of $\boldsymbol{x}_{\infty}$ and $\tilde{\boldsymbol{y}}$ from Eqs. (6.14) and (6.1), respectively, to derive the following inequalities:

[TABLE]

Noting that $\|\boldsymbol{x}-\tilde{\boldsymbol{y}}\|_{\infty}\geq\|\boldsymbol{x}(t^{n})-\tilde{\boldsymbol{y}}(t^{n})\|_{2}$, $\forall n\in{\rm Nature}{{N_{t}}}$ and substituting in the definitions of $L_{\boldsymbol{r}}$ and $K_{\boldsymbol{r}}$ from Eqs. (6.3) and (6.8), respectively, yields the stated result. ∎

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

Corollary 5** (Simplified a priori error bound).**

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

[TABLE]

[TABLE]

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

[TABLE]

Proof.

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

Remark 6.1** (Stability-constant growth).**

Figure 1 plots the dependence of the stability constants in the a priori error bounds (6.16) and (6.17) as a function of the final time $T$ for multiple linear multistep methods and fixed values of the time step and Lipschitz constant. Critically, note that the stability constant for the $\ell^{2}$-norm of the error in the ST-LSPG ROM solution grows only linearly in time, while the stability constant for the $\ell^{\infty}$-norm of the error in the ST-LSPG ROM solution exhibits polynomial growth in time with degree $3/2$. Further, this trend is valid for all assessed linear multistep schemes. This highlights one important feature of the proposed method: significantly slower time growth of the solution error in time relative to nonlinear model-reduction methods that perform only spatial projection, as such error bounds grow exponentially in time [38, 28, 34, 11]. This is similar to the slow time growth of the error bounds demonstrated in the context of the space–time reduced-basis method [45, 46, 50, 49].

We now provide computable a posteriori residual-based error bounds and show that the ST-LSPG solution minimizes this bound over all solutions in the space–time trial subspace.

Corollary 6** **(a posteriori error

bound).

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

[TABLE]

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

[TABLE]

Proof.

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

[TABLE]

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

7 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. [9] 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. [14].

• •

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

We assess the accuracy of any ROM solution $\tilde{\boldsymbol{x}}(\cdot;\boldsymbol{\mu})$ from its mean squared state-space error, i.e.,

[TABLE]

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 [51] in MATLAB. All reported timings are averaged over five simulations.

7.1 Parameterized Burgers’ equation

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

[TABLE]

where $w:[0,1]\times[0,T]\times\mathcal{D}\rightarrow\mathbb{R}$ is a conserved quantity and the $n_{\mu}=2$ parameters comprise the left boundary value and source-term coefficient with $\boldsymbol{\mu}\equiv(\mu_{1},\mu_{2})\in\mathcal{D}=[1.2,1.5]\times[0.02,0.025]$.

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.

7.1.1 Space–time bases

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. [6]. 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. [8]).

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

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

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

[TABLE]

where $T=0.6$ s and

[TABLE]

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

[TABLE]

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:

[TABLE]

where a subscript $m$ indicates the flow quantity at $x=0.5$ m, and $M$ denotes the Mach number. Then, a shock is placed at $x=0.85$ m of the flow field. We use the jump relations for a stationary shock and the perfect gas equation of state to derive the velocity across the shock $u_{2}$, which satisfies the quadratic equation

[TABLE]

Here, $m:=\rho_{2}u_{2}=\rho_{1}u_{1}$, $n:=\rho_{2}u_{2}^{2}+p_{2}=\rho_{1}u_{1}^{2}+p_{1}$, $h:=(e_{2}+p_{2})/\rho_{2}=(e_{1}+p_{1})/\rho_{1}$, and subscripts $1$ and $2$ denote a flow quantity to the left and to the right of the shock, respectively. We employ the solution $u_{2}$ to Eq. (7.7), which leads to a discontinuity (i.e., shock). Finally, the exit pressure is increased to a factor $P_{\text{exit}}$ of its original value in order to generate transient dynamics.

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

7.2.1 Space–time bases

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.

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

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

8 Conclusions

This work proposed a model-reduction method for nonlinear dynamical systems based on space–time least-squares Petrov–Galerkin projection. The method computes optimal approximations by minimizing the discrete space–time residual over all elements in a low-dimensional space–time trial subspace in a weighted $\ell^{2}$-norm. Advantages of the method include:

• •

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.