# SNS: A Solution-based Nonlinear Subspace method for time-dependent model order reduction

Several reduced order models have been successfully developed for nonlinear dynamical systems. To achieve a considerable speed-up, a hyper-reduction step is needed to reduce the computational complexity due to nonlinear terms. Many hyper-reduction techniques require the construction of nonlinear term basis, which introduces a computationally expensive offline phase. A novel way of constructing nonlinear term basis within the hyper-reduction process is introduced. In contrast to the traditional hyper-reduction techniques where the collection of nonlinear term snapshots is required, the SNS method avoids collecting the nonlinear term snapshots. Instead, it uses the solution snapshots that are used for building a solution basis, which enables avoiding an extra data compression of nonlinear term snapshots. As a result, the SNS method provides a more efficient offline strategy than the traditional model order reduction techniques, such as the DEIM, GNAT, and ST-GNAT methods. The SNS method is theoretically justified by the conforming subspace condition and the subspace inclusion relation. It is especially useful for ST-GNAT that has shown promising results, such as a good accuracy with a considerable online speed-up for hyperbolic problems in a recent paper, because ST-GNAT involves an expensive offline cost related to collecting nonlinear term snapshots. Error analysis shows that the oblique projection error bound of the SNS method depends on the condition number of the volume matrix generated from a discretization of a specific numerical scheme. Numerical results support that the accuracy of the solution from the SNS method is comparable to the traditional methods and a considerable speed-up (i.e., a factor of two to a hundred) is achieved in the offline phase.

Comments

There are no comments yet.

## Authors

• 20 publications
• 2 publications
• 6 publications
09/11/2018

### SNS: A Solution-based Nonlinear Subspace method for time-dependent nonlinear model order reduction

Several reduced order models have been successfully developed for nonlin...
01/14/2021

### An EIM-degradation free reduced basis method via over collocation and residual hyper reduction-based error estimation

The need for multiple interactive, real-time simulations using different...
04/01/2019

### Finite strain homogenization using a reduced basis and efficient sampling

The computational homogenization of hyperelastic solids in the geometric...
11/13/2020

### Efficient nonlinear manifold reduced order model

Traditional linear subspace reduced order models (LS-ROMs) are able to a...
03/13/2020

### Model Reduction of Time-Dependent Hyperbolic Equations using Collocated Residual Minimisation and Shifted Snapshots

We develop a non-linear approximation for solution manifolds of parametr...
03/02/2021

### Projection-tree reduced order modeling for fast N-body computations

This work presents a data-driven reduced-order modeling framework to acc...
04/26/2022

### A Reduced Order Model for Joint Assemblies by Hyper-Reduction and Model-Driven Sampling

The dynamic behavior of jointed assemblies exhibiting friction nonlinear...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Time-dependent nonlinear problems arise in many important disciplines such as engineering, science, and technologies. They are numerically solved if it is not possible to solve them analytically. Depending on the complexity and size of the governing equations and problem domains, the problems can be computationally expensive to solve. It may take a long time to run one forward simulation even with high performance computing. For example, a simulation of the powder bed fusion additive manufacturing procedure shown in [27] takes a week to finish with cores. Other computationally expensive simulations include the 3D shocked spherical Helium bubble simulation appeared in [4] and the inertial confinement fusion implosion dynamics simulations appeared in [1]. The computationally expensive simulations are not desirable in the context of parameter study, design optimization, uncertainty quantification, and inverse problems where several forward simulations are needed. A Reduced Order Model (ROM) can be useful in this context to accelerate the computationally expensive forward simulations with good enough approximate solutions.

We consider projection-based ROMs for nonlinear dynamical systems. Such ROMs include the Empirical Interpolation method (EIM)

[6, 22], the Discrete Empirical Interpolation Method (DEIM) [12, 15]

and the Gauss–Newton with Approximated Tensors

[10, 11], the Best Point Interpolation Method (BPIM) [31]

, the Missing Point Estimation (MPE)

[5], and Cubature Methods [3, 19, 20, 24]. There a hyper-reduction (The term first used in the context of ROMs by Ryckelynck in [33]

) is necessary to efficiently reduce the complexity of nonlinear terms for a considerable speed-up compared to a corresponding full order model. The DEIM and GNAT approaches take discrete nonlinear term snapshots to build a nonlinear term basis. Then, they select a subset of each nonlinear term basis vector to either interpolate or data-fit in a least-squares sense. In this way, they reduce the computational complexity of updating nonlinear terms in an iterative solver for nonlinear problems. The EIM, BPIM, and MPE approaches take the similar hyper-reduction to the ones in DEIM and GNAT except for the fact that they build nonlinear term basis functions in a continuous space. Whether they work on a discrete or continuous space, they follow the framework of reconstructing “gappy” data, first introduced in the context of reduced order models by

[18]. The cubature methods, in contrast, takes a different approach. They approximate nonlinear integral as a finite sum of positive scalar weights multiplied by the integrand evaluated at sampled elements. The cubature methods developed in [3, 19, 20] do not require building nonlinear term basis. They solve the Non-Negative Least-Squares (NNLS) problem directly with the nonlinear term snapshots. Recently, Hernandez, et al., in [24] developed the Empirical Cubature Method (ECM) approach that builds nonlinear term basis to solve a smaller NNLS problem. The requirement of building nonlinear term basis in hyper-reduction results in computational cost and storage in addition to solution basis construction. The cost is significant if the corresponding Full Order Models (FOMs) are large-scale. The large-scale problem requires large additional storage for nonlinear term snapshots and large-scale compression techniques to build a basis. In particular, the cost of the hyper-reduction in the recently-developed space–time ROM (i.e., ST-GNAT) [13] is even more significant than aforementioned spatial ROMs (e.g., EIM, DEIM, GNAT, BPIM, MPE, and ECM). The best nonlinear term snapshots of the ST-GNAT method are obtained from the corresponding space–time ROMs without hyper-reduction (i.e., ST-LSPG) [13], which is not practical for a large-scale problem due to the cumbersome size.111

The full size, implicitly handled by the ST-LSPG approach due to the nonlinear term and corresponding Jacobian updates, is the FOM spatial degrees of freedom multiplied by the number of time steps.

Therefore, a novel and efficient way of constructing nonlinear term basis needs to be developed.

This paper shows a practical way of avoiding nonlinear term snapshots for the construction of nonlinear term basis. The idea comes from a simple fact that the nonlinear terms are related with solution snapshots through underlying time integrators. In fact, many time integrators approximate time derivative terms as a linear combination of the solution snapshots. It implies that the nonlinear term snapshots belong to the subspace spanned by the solution snapshots. Furthermore, a subspace needed for nonlinear terms in a hyper-reduction is determined by the range space of the solution basis matrix possibly multiplied by a nonsingular matrix (e.g., the volume matrix). Therefore, the solution snapshots can be used to construct nonlinear term basis. This leads to our proposed method, the Solution-based Nonlinear Subspace (SNS) method, that provides two savings for constructing nonlinear term basis because of

1. no additional collection of nonlinear term snapshots (i.e., storage saving).

2. no additional compression of snapshots (e.g., no additional singular value decomposition of nonlinear term snapshots, implying computational cost saving).

The first saving is especially big for GNAT and ST-GNAT becuase they involve expensive collection procedure for their best performace.

### 1.1 Organization of the paper

We start our discussion by describing the time-continuous representation of the FOM in Section 2. Section 2 also describes the time-discrete representation of the FOM with one-step Euler time integrators. The subspace inclusion relation between the subspaces spanned by the solution and nonlinear term snapshots is described for the Euler time integrators. Several projection-based ROMs (i.e., the DEIM, GNAT, and ST-GNAT models) are considered in Section 3 for the SNS method to be applied. The SNS method is described in Section 4 and applied to those ROMs. Section 4 also introduces the conforming subspace condition to justify the SNS method. Section 5 shows an error analysis for the SNS method regarding the oblique projection error bound. It analyzes the effect of the nonsingular matrix that is used to form a nonlinear term basis. Section 6 reports numerical results that support benefits of the SNS method and Section 7 concludes the paper with summary and future work. Although the SNS method is mainly illustrated with the Euler time integrators throughout the paper, it is applicable to other time integrators. Appendix A considers several other time integrators and the subspace inclusion relation for each time integrator. The following time integrators are included: the Adams–Bashforth, Adams–Moulton, BDF, and midpoint Runge–Kutta methods.

## 2 Full order model

A parameterized nonlinear dynamical system is considered, characterized by a system of nonlinear ordinary differential equations (ODEs), which can be considered as a resultant system from semi-discretization of Partial Differential Equations (PDEs) in space domains

 (2.1) Mdxdt=f(x,t;μ),x(0;μ)=x0(μ),

where denotes a nonsingular matrix, denotes time with the final time , and denotes the time-dependent, parameterized state implicitly defined as the solution to problem (2.1) with . Further, with denotes the scaled velocity of , which we assume to be nonlinear in at least its first argument. The initial state is denoted by , and denotes the parameters with parameter domain .

A uniform time discretization is assumed throughout the paper, characterized by time step and time instances for with , , and . To avoid notational clutter, we introduce the following time discretization-related notations: , , , and .

Two different types of time discretization methods are considered: explicit and implicit time integrators. As an illustration purpose, we mainly consider the forward Euler time integrator for an explicit scheme and the backward Euler time integrator for an implicit scheme. Several other time integrators are shown in Appendix A.

The explicit Forward Euler (FE) method numerically solves Eq. (2.1), by time-marching with the following update:

 (2.2) Mxn−Mxn−1=Δtfn−1.

Eq. (2.2) implies the following subspace inclusion:

 (2.3) span{fn−1}⊆span{Mxn−1,Mxn}.

By induction, we conclude the following subspace inclusion relation:

 (2.4) span{f0,…,fNt−1}⊆span{Mx0,…,MxNt},

which shows that the span of nonlinear term snapshots is included in the span of -scaled solution snapshots. The residual function with the forward Euler time integrator is defined as

 (2.5) rnFE(xn;xn−1,μ):=M(xn−xn−1)−Δtfn−1

The implicit Backward Euler (BE) method numerically solves Eq. (2.1), by solving the following nonlinear system of equations for at -th time step:

 (2.6) Mxn−Mxn−1=Δtfn.

Eq. (2.6) implies the following subspace inclusion:

 (2.7) span{fn}⊆span{Mxn−1,Mxn}.

By induction, we conclude the following subspace inclusion relation:

 (2.8) span{f1,…,fNt}⊆span{Mx0,…,MxNt},

which shows that the span of nonlinear term snapshots is included in the span of -scaled solution snapshots. The residual function with the backward Euler time integrator is defined as

 (2.9) rnBE(xn;xn−1,μ):=M(xn−xn−1)−Δtfn.

## 3 Projection-based reduced order models

Projection-based reduced order models are considered for nonlinear dynamical systems. Especially, the ones that require building a nonlinear term basis are of our interest: the DEIM, GNAT, and ST-GNAT approaches.

### 3.1 Deim

The DEIM approach applies spatial projection using a subspace with . Using this subspace, it approximates the solution as (i.e., in a trial subspace) or equivalently

 (3.1) ~x=x0+Φ^x

where denotes a basis matrix and with denotes the generalized coordinates. Replacing with in Eq. (2.1) leads to the following system of equations with reduced number of unknowns:

 (3.2) MΦd^xdt=f(x0+Φ^x,t;μ).

For constructing , Proper Orthogonal Decomposition (POD) is commonly used. POD [7] obtains

from a truncated Singular Value Decomposition (SVD) approximation to a FOM solution snapshot matrix. It is related to principal component analysis in statistical analysis

[26] and Karhunen–Loève expansion [29] in stochastic analysis. POD forms a solution snapshot matrix, , where is a solution state at th time step with parameter for and . Then, POD computes its thin SVD:

 (3.3) X=UΣVT,

where and are orthogonal matrices and is a diagonal matrix with singular values on its diagonals. Then POD chooses the leading columns of to set (i.e., in MATLAB notation). The POD basis minimizes over all with orthonormal columns, where denotes the Frobenius norm of a matrix , defined as with being an -th element of . Since the objective function does not change if is post-multiplied by an arbitrary orthogonal matrix, the POD procedure seeks the optimal –-dimensional subspace that captures the snapshots in the least-squares sense. For more details on POD, we refer to [25, 28].

Note that Eq. (3.2) has more equations than unknowns (i.e., an overdetermined system). It is likely that there is no solution satisfying Eq. (3.2). In order to close the system, the Galerkin projection solves the following reduced system with :

 (3.4) ΦTMΦd^xdt=ΦTf(x0+Φ^x,t;μ).

Applying a time integrator to Eq. (3.4) leads to a fully discretized reduced system, denoted as the reduced OE. Note that the reduced OE has unknowns and equations. If an implicit time integrator is applied, a Newton–type method can be applied to solve for unknown generalized coordinates each time step. If an explicit time integrator is applied, time marching updates will solve the system. However, we cannot expect any speed-up because the size of the nonlinear term and its Jacobian, which need to be updated for every Newton step, scales with the FOM size. Thus, to overcome this issue, the DEIM approach applies a hyper-reduction technique. That is, it projects onto a subspace and approximates as

 (3.5) f≈Φf^f,

where , , denotes the nonlinear term basis matrix and denotes the generalized coordinates of the nonlinear term. The generalized coordinates, , can be determined by the following interpolation:

 (3.6) ZTf=ZTΦf^f,

where is the sampling matrix and is the

th column of the identity matrix

. Therefore, Eq. (3.5) becomes

 (3.7) f≈PDEIMf,

where is the DEIM oblique projection matrix. The DEIM approach does not construct the sampling matrix . Instead, it maintains the sampling indices and corresponding rows of and . This enables DEIM to achieve a speed-up when it is applied to nonlinear problems.

The original DEIM paper [12] constructs the nonlinear term basis by applying another POD on the nonlinear term snapshots222the nonlinear term snapshots are with the backward Euler time integrator and with the forward Euler time integrator obtained from the FOM simulation at every time step. This implies that DEIM requires two separate SVDs and storage for two different snapshots (i.e., solution state and nonlinear term snapshots). Section 4 discusses how to avoid the collection of nonlinear term snapshots and an extra SVD without losing the quality of the hyper-reduction.

The sampling indices (i.e., ) can be found either by a row pivoted LU decomposition [12] or the strong column pivoted rank-revealing QR (sRRQR) decomposition [15, 16]. Depending on the algorithm of selecting the sampling indices, the DEIM projection error bound is determined. For example, the row pivoted LU decomposition in [12] results in the following error bound:

 (3.8) ∥f−Pf∥2≤κ∥(INs−ΦfΦTf)f∥2,

where is the condition number of and it is bounded by

 (3.9) κ≤(1+√2Ns)nf−1∥ϕf,1∥−1∞.

On the other hand, the sRRQR factorization in [16] reveals tighter bound than (3.9):

 (3.10) κ≤√1+η2nf(Ns−nf)

where is a tuning parameter in the sRRQR factorization (i.e., in Algorithm 4 of [23]).

### 3.2 Gnat

In contrast to DEIM, the GNAT method takes the Least-Squares Petrov–Galerkin (LSPG) approach. The LSPG method projects a fully discretized solution space onto a trial subspace. That is, it discretizes Eq. (2.1) in time domain and replaces with for in residual functions defined in Section 2 and Appendix A. Here, we consider only implicit time integrators because the LSPG projection is equivalent to the Galerkin projection when an explicit time integrator is used as shown in Section 5.1 in [9]. The residual functions for implicit time integrators are defined in (2.9), (A.8), and (A.12) for various time integrators. For example, the residual function with the backward Euler time integrator333Although the backward Euler time integrator is used extensively in the paper as an illustration purpose, many other time integrators introduced in Appendix A can be applied to all the ROM methods dealt in the paper in a straight forward way. after the trial subspace projection becomes

 (3.11) ~rnBE(^xn;^xn−1,μ):=rnBE(x0+Φ^xn;x0+Φ^xn−1,μ)=MΦ(^xn−^xn−1)−Δtf(x0+Φ^xn,t;μ).

The basis matrix can be found by the POD as in the DEIM approach. Note that Eq. (3.11) is an over-determined system. To close the system and solve for the unknown generalized coordinates, , the LSPG ROM takes the square norm of the residual vector function and minimize it at every time step:

 (3.12) ^xn=argmin^v∈Rns12∥∥rnBE(^v;^xn−1,μ)∥∥22.

The Gauss–Newton method with the starting point is applied to solve the minimization problem (3.12) in GNAT. However, as in the DEIM approach, a hyper-reduction is required for a speed-up due to the presence of the nonlinear residual vector function. The GNAT method approximates the nonlinear residual term with gappy POD [18], whose procedure is similar to DEIM, as

 (3.13) r≈Φr^r,

where , , denotes the residual basis matrix and denotes the generalized coordinates of the nonlinear residual term. In contrast to DEIM, the GNAT method solves the following least-squares problem to obtain the generalized coordinates at n-th time step:

 (3.14) ^rn=argmin^v∈Rnr12∥∥ZT(r−Φr^v)∥∥22.

where , , is the sampling matrix and is the th column of the identity matrix . The solution to Eq. (3.14) is given as

 (3.15) ^rn=(ZTΦr)†ZTr,

where the Moore–Penrose inverse of a matrix with full column rank is defined as . Therefore, Eq. (3.13) becomes

 (3.16) r≈PGNATr,

where is the GNAT oblique projection matrix. Note that it has a similar structure to . The GNAT projection matrix has a pseudo-inverse instead of the inverse becuase it allows . The GNAT method does not construct the sampling matrix . Instead, it maintains the sampling indices and corresponding rows of and . This enables GNAT to achieve a speed-up when it is applied to nonlinear problems.

The sampling indices (i.e., ) can be determined by Algorithm 3 of [11] for computational fluid dynamics problems and Algorithm 5 of [10] for other problems. These two algorithms take greedy procedure to minimize the error in the gappy reconstruction of the POD basis vectors . The major difference between these sampling algorithms and the ones for the DEIM method is that these algorithms for the GNAT method allows oversampling (i.e., ), resulting in solving least-squares problems in the greedy procedure. These selection algorithms can be viewed as the extension of Algorithm 1 in [12] (i.e., a row pivoted LU decomposition) to the oversampling case. The nonlinaer residual term projection error associated with these sampling algorithms is presented in Appendix D of [11]. That is,

 (3.17) ∥r−PGNATr∥2≤∥R−1∥2∥r−ΦrΦTrr∥2

where is a triangle matrix from QR factorization of (i.e., ).

We present another sampling selection algorithm that was not considered in any GNAT papers (e.g., [10, 11]). It is to use the sRRQR factorization developed originally in [23] and further utilized in the W-DEIM method of [16] for the case of . That is, applying Algorithm 4 of [23] to produces an index selection operator whose projection error satisfies

 (3.18) ∥r−PGNATr∥2≤√1+η2nr(Ns−nr)∥r−ΦrΦTrr∥2.

This error bound can be obtained by setting the identity matrix as a weight matrix in Theorem 4.8 of [16].

Finally, the GNAT method minimizes the following least-squares problem at every time step, for example, with the backward Euler time integrator:

 (3.19) ^xn=argmin^v∈Rns12∥∥ (ZTΦr)†ZTrnBE(^v;^xn−1,μ)∥∥22.

The GNAT method applies another POD to a nonlinear residual term snapshots to construct . The original GNAT paper [11] collects residual snapshots at every Newton iteration from the LSPG simulations444We denote the LSPG simulation to be the ROM simulation without hyper-reduction. That is, LSPG solves Eq. (3.12) without any hyper-reduction if the backward Euler time integrator is used. for several reasons:

1. The GNAT method takes LSPG as a reference model (i.e., denoted as Model II in [11]). Its ultimate goal is to achieve the accuracy of Model II.

2. The residual snapshots taken every time step (i.e., at the end of Newton iterations at every time step) of the FOM are small in magnitude.

3. The residual snapshots taken from every Newton step gives information about the path that the Newton iterations in Model II take. GNAT tries to mimic the Newton path that LSPG takes.

Some disadvantages of the original GNAT approach include:

1. The GNAT method requires more storage than DEIM to store residual snapshots from every Newton iteration (cf., The DEIM approach stores only one nonlinear term snapshot per each time step).

2. The GNAT method requires more expensive SVD for nonlinear residual basis construction than the DEIM approach because the number of nonlinear residual snapshots in the GNAT method is larger than the number of nonlinear term snapshots in DEIM.

3. The GNAT method requires the simulations of Model II which are computationally expensive. For a parametric global ROM, it is computationally expensive because it requires as many Model II simulations as there are training points in a given parameter domain.

Section 4 discusses how to avoid the collection of nonlinear term snapshots and the extra SVD without losing the quality of the hyper-reduction.

### 3.3 Space–time GNAT

The ST-GNAT method takes the space–time LSPG (ST-LSPG) approach. Given a time integrator, one can rewrite a fully discretized system of equations to Eq. (2.1) in a space–time form. For example, if the backward Euler time integrator555 Here we only consider the backward Euler time integrator for the simplicity of illustration. However, the extension to other linear multistep implicit time integrators is straight forward. is used for time domain discretization, then the corresponding space–time formulation to Eq. (2.6) becomes

 (3.20) ABE¯x=Δt¯f+¯q0,

where

 (3.21)

Note that denotes the parameterized space–time state implicitly defined as the solution to the problem (3.20) with and . Further, denotes the space–time nonlinear term that is nonlinear in at least its first argument with . The space–time residual function corresponding to Eq. (3.20) is defined as

 (3.22) ¯rBE(¯x;μ):=ABE¯x(μ)−Δt¯f(¯x;μ)−¯q0(μ).=0

To reduce both the spatial and temporal dimensions of the full-order model, the ST-LSPG method enforces the approximated numerical solution to reside in an affine ‘space–time trial subspace’

 (3.23) ~y∈ST⊆RNsNt,

where with and whose elements are all one. Here, denotes the Kronecker product of two matrices and the product of two matrices and is defined as

 (3.24) A⊗B=⎡⎢ ⎢⎣a11B⋯a1JB⋮⋱⋮aI1B⋯aIJB⎤⎥ ⎥⎦.

Further, denotes a space–time basis vector. Thus, the ST-LSPG method approximates the numerical solution as

 (3.25) ¯x(μ)≈~y(μ)=1Nt⊗x0(μ)+nst∑i=1¯ϕi^yi(μ)

where , denotes the generalized coordinate of the ST-LSPG solution.

A space–time residual vector function can now be defined with the generalized coordinates as an argument. Replacing in Eq. (3.22) with gives the residual vector function

 (3.26) ¯rBE(^y;μ):=ABE¯Φ^y−Δt¯f(¯x0+¯Φ^y;μ)−¯q0(μ)+ABE¯x0(μ),

where and denotes a space–time basis matrix that is defined as and denotes the generalized coordinate vector that is defined as

 (3.27) ^y:=[^y1(μ)⋯^ynst(μ)]T.

Note that vanishes. Eq. (3.26) becomes

 (3.28) ¯rBE(^y;μ):=ABE¯Φ^y−Δt¯f(¯x0+¯Φ^y;μ).

Reduced space–time residual vector functions for other time integrators can be defined similarly. We denote as a reduced spate–time residual vector function with a generic time integrator.

The space–time basis matrix can be obtained by tensor product of spatial and temporal basis vectors. The spatial basis vectors can be obtained via POD as in the DEIM and GNAT approaches. The temporal basis vectors can be obtained via the following three tensor decompositions described in [13]:

• Fixed temporal subspace via T-HOSVD

• Fixed temporal subspace via ST-HOSVD

• Tailored temporal subspace via ST-HOSVD

The ST-HOSVD method is a more efficient version of T-HOSVD. Thus, we will not consider T-HOSVD. The tailored temporal subspace via ST-HOSVD has a LL1 form that has appeared, for example, in [34, 14]. Therefore, we will refer to it as the LL1 decomposition.

Because of the reduction in spatial and temporal dimension, the space–time residual vector cannot achieve zero most likely. Thus, the ST-LSPG method minimizes the square norm of and computes the ST-LSPG solution:

 (3.29) ^y(μ)=argmin^v∈Rnst∥¯r(^v;μ)∥22.

The ST-LSPG ROM solves Eq. (3.29) without any hyper-reduction. As in the DEIM and GNAT approaches, a hyper-reduction is required for a considerable speed-up due to the presence of the space–time nonlinear residual vector. Therefore, the ST-GNAT method approximates the space–time nonlinear residual terms with gappy POD [18], which in turn requires construction of a space–time residual basis. Similar to the GNAT method, the ST-GNAT method approximates the space–time nonlinear residual term as

 (3.30) ¯r≈¯Φr^r,

where , , denotes the space–time residual basis matrix and denotes the generalized coordinates of the nonlinear residual term. The ST-GNAT solves the following space–time least-squares problem to obtain the generalized coordinates, :

 (3.31) ^r=argmin^v∈Rnr12∥∥¯ZT(¯r−¯Φr^v)∥∥22.

where , , is the sampling matrix and is the th column of the identity matrix . The solution to Eq. (3.31) is given by

 (3.32) ^r=(¯ZT¯Φr)†¯ZT¯r.

Therefore, Eq. (3.30) becomes

 (3.33) ¯r≈PST-GNAT¯r,

where is the ST-GNAT oblique projection matrix. Note that has the same structure as . The ST-GNAT method does not construct the sampling matrix . Instead, it maintains the sampling indices and corresponding rows of and . This enables the ST-GNAT to achieve a speed-up when it is applied to nonlinear problems.

Section 5.3 in [13] discusses three different options to determine the sampling indices (i.e., ). However, all these three options are simple variations of Algorithm 3 in [11] and Algorithm 5 in [10]. They all minimize the error in the gappy reconstruction of the POD basis vectors for nonlinear space and time residuals. Therefore, the space–time nonlinear residual projection error due to (3.33) is similar to the ones in the GNAT method. That is,

 (3.34) ∥¯r−¯Φr(¯ZT¯Φr)†¯ZT¯r∥2≤∥R−1∥2∥¯r−¯Φr¯ΦTr¯r∥2,

where is a triangle matrix from QR factorization of (i.e., ). On the other hand, one can also apply the sRRQR factorization in Algorithm 4 with a tuning parameter of [23] to to obtain that is associated with a tighter error bound for the projection error:

 (3.35) ∥¯r−¯Φr(¯ZT¯Φr)†¯ZT¯r∥2≤√1+η2nr(Ns−nr)∥¯r−¯Φr¯ΦTr¯r∥2.

This error bound can be obtained by setting the identity matrix as a weight matrix in Theorem 4.8 of [16].

Finally, the ST-GNAT solves the following least-squares problem at every time step, for example, with the backward Euler time integrator:

 (3.36) ^y(μ)=argmin^v∈Rns12∥∥ (¯ZT¯Φr)†¯ZT¯rBE(^v;μ)∥∥22.

The original ST-GNAT paper introduces three different ways of collecting space–time residual snapshots, that are in turn used for the space–time residual basis construction (see Section 5.2 in [13]). Below is a list of the approaches introduced in [13] and explains advantages and disadvantages of each:

1. ST-LSPG ROM training iterations. This approach takes the space–time residual snapshot from every Newton iteration of the ST-LSPG simulations at training points in . This case leads to the number of residual snapshots, , where is the number of Newton iterations taken to solve the ST-LSPG simulation for the training point . This approach is not realistic for a large-scale problem because it requires training simulations of the computationally expensive ST-LSPG ROM, where denotes the cardinality of the set . Furthermore, it requires the extra SVD on the residual snapshots, which is not necessary for our proposed method.

2. Projection of FOM training solutions. This approach takes the following steps:

1. take FOM state solution at every Newton iteration.

2. re-arrange them in the space–time form (i.e., in Eq. (3.21)).666

Note that these are FOM solutions from time marching algorithms in which each time step results in different number of Newton iterations if implicit time integrators are used. Some time steps take a smaller number of Newton iterations than other time steps. However, in order to re-arrange each Newton iterate state in the space–time form, we must have the same number of Newton iterations at each time step. Therefore, for the time steps that have converged with a smaller number of Newton iterations than other time steps, we pad the solution state vectors of the Newton iterations beyond the convergence with the converged solution. This only applies to an implicit time integrator because an explicit time integrator does not require any Newton solve.

3. project them onto the space–time subspace, (i.e., ).

4. compute the corresponding space–time residual (e.g., in Eq. (3.22) in the case of the backward Euler time integrator).

5. use those residuals as residual snapshots.

This approach simply requires projections and evaluations of the space–time residual. However, it requires the extra SVD on the residual snapshots, which is not necessary for our proposed method.

3. Random samples.

This approach generates a random space–time solution state samples (e.g., via Latin hypercube sampling or random sampling from uniform distribution) and computes the corresponding space–time residual (e.g.,

in Eq. (3.22) in the case of the backward Euler time integrator). This approach simply requires random sample generations and evaluations of the space–time residual. However, random samples are hardly correlated with actual data. Therefore, it is likely to generate poor space–time residual subspace. Furthermore, it requires the extra SVD on the residual snapshots, which is not necessary for our proposed method.

## 4 Solution-based Nonlinear Subspace (SNS) method

Finally, we state our proposed method that avoids collecting nonlinear term snapshots and additional POD for the DEIM and GNAT approaches or additional tensor decomposition for the ST-GNAT method. We propose to use solution snapshots to construct nonlinear term basis in the DEIM, GNAT, and ST-GNAT approaches. A justification for using the solution snapshots comes from the subspace inclusion relation between the subspace spanned by the solution snapshots and the subspace spanned by the nonlinear term snapshots as shown in Eqs. (2.4) and (2.8) with the forward and backward Euler time integrators.777Subspace inclusion relations for other time integrators are shown in Appendix A.

### 4.1 Deim-Sns

We are going back to Eq. (3.2) and replace the nonlinear term with the approximation in Eq. (3.7):

 (4.1) MΦd^xdt=Φf(ZTΦf)−1ZTf(x0+Φ^x,t;μ).

Eq. (4.1) is an over-determined system, so it is likely that it will not have a solution. However, if there is a solution, then necessary conditions for Eq. (4.1) to have a non-trivial solution (i.e., ) are and

 (4.2) Ran(MΦ)∩Ran(Φf)≠{0}.

The second condition says that the intersection of and needs to be non-trivial if there is a non-trivial solution to Eq. (4.1). Typically, we build first, using the POD approach explained in Section 3.1 and is set by . Therefore, the intersection of and can be controlled by the choice of we made. The larger the intersection of those two range spaces are, the more likely it is that there is a solution to Eq. (4.1). Given , the largest subspace intersection it can be is , i.e.,

 (4.3) Ran(MΦ)∩Ran(Φf)=Ran(MΦ).

We call this condition as the conforming subspace condition. The conforming subspace condition leads to two obvious choices for :

• The first choice is to ensure . If , then the range space of the left and right-hand sides of Eq. (4.1) are the same. This leads Eq. (4.1) to become

 (4.4) MΦd^xdt=MΦ(ZTMΦ)−1ZTf(x0+Φ^x,t;μ).

Because Eq. (4.4) is an over-determined system and unlikely to have a solution, applying the Galerkin projection to Eq. (4.4) becomes:

 (4.5) ΦTMΦd^xdt=ΦTMΦ(ZTMΦ)−1ZTf(x0+Φ^x,t;μ).

Assuming is invertible, Eq. (4.5) becomes:

 (4.6) d^xdt=(ZTMΦ)−1ZTf(x0+Φ^x,t;μ).

For the special case of being an identity matrix, Eq. (4.6) becomes:

 (4.7) d^xdt=(ZTΦ)−1ZTf(x0+Φ^x,t;μ).
• The second choice is to ensure . This can be achieved by taking an extended solution basis, with and . Then we set , which leads to . The obvious choice for is to take a larger truncation of the left singular matrix from SVD of the solution snapshot matrix than . Note that this particular choice of results in the first columns of being the same as (i.e., where with and ). By setting , Eq. (4.1) becomes

 (4.8) MΦd^xdt=MΦe(ZTMΦe)−1ZTf(x0+Φ^x,t;μ).

Because Eq. (4.8) is unlikely to have a solution, applying the Galerkin projection to Eq. (4.8) becomes:

 (4.9) ΦTMΦd^xdt =ΦTMΦe(ZTMΦe)−1ZTf(x0+Φ^x,t;μ) (4.10) =ΦTM[ΦΦE](ZTMΦe)−1ZTf(x0+Φ^x,t;μ).

Assuming is invertible, Eq. (4.9) becomes:

 (4.11) d^xdt=[Ins(ΦTMΦ)−1(ΦTMΦE)](ZTMΦe)−1ZTf(x0+Φ^x,t;μ).

For the special case of being an identity matrix, Eq. (4.11) becomes

 (4.12) d^xdt=[Ins0](ZTΦe)−1ZTf(x0+Φ^x,t;μ).

The DEIM-SNS approach solves either Eq. (4.6) or Eq. (4.11) depending on the choice of above. Applying a time integrator to Eq. (4.6) or Eq. (4.12) leads to a reduced OE, whose solution, , can be lifted to find the full order approximate solution via .

Additionally, the subspace inclusion relations888 Eq. (2.4) of the forward Euler time integrator, Eq. (A.3) of the Adams–Moulton time integrator, and Eq. (A.16) of the midpoint Runge-Kutta method show the subspace inclusion relations for explicit time integrators. show that the subspace spanned by the solution snapshots include the subspace spanned by the nonlinear term snapshots. This fact further motivates the use of solution snapshots to build a nonlinear term basis. Indeed, numerical experiments show that the solution accuracy obtained by the DEIM-SNS approach is comparable to the one obtained by the traditional DEIM approach. For example, see Figs. 3, and 4 in Section 6.1.1.

The obvious advantage of the DEIM-SNS approach over DEIM is that no additional SVD or eigenvalue decomposition is required, which can save the computational cost of the offline phase.

###### Remark 4.1

Although the numerical experiments show that the two choices of above give promising results, the error analysis in Section 5 shows that the nonlinear term projection error bound increases by the condition number of (see Theorem 5). This is mainly because the orthogonality of is lost when or . This issue is resolved by orthogonalizing (e.g., apply economic QR factorization ) before using it as nonlinear term basis. Then apply, for example, Algorithm 4 of [23] to the transpose of orthogonalized one (e.g., ) to generate a sampling matrix . This procedure eliminates the condition number of in the error bound because (3.10) is valid.

###### Remark 4.2

Inspired by the weighted inner product space introduced for DEIM in [16], another oblique DEIM-SNS projection is possible. With the weight matrix , the selection operator , and the basis according to Section 4.2 of [16], the weighted oblique DEIM-SNS projection can be defined as:

 (4.13) PSNS =^U(STW^U)†STW (4.14) =MΦ(ZTΦ)†ZTM−1.

### 4.2 Gnat-Sns

The GNAT method needs to build a nonlinear residual term basis, , as described in Section 3.2. The nonlinear residual term is nothing more than linear combinations of the nonlinear term and the time derivative approximation as in Eq. (3.11). Thus, the subspace spanned by the nonlinear term residual snapshots are included in the subspace spanned by the solution snapshots. This motivates the use of the solution snapshots for the construction of a nonlinear residual term basis. Therefore, the same type of the nonlinear term basis in the DEIM-SNS approach can be used to construct the nonlinear residual term basis in the GNAT-SNS method:

• The first choice is to set . Thus, for example, the GNAT-SNS method solves the following least-squares problem with the backward Euler time integrator:

 (4.15) ^xn=argmin^v∈Rns12∥∥MΦ(ZTMΦ)†ZT(MΦ(^v−^xn−1)−Δtf(x0+Φ^xn,t;μ))∥∥22,

which can be manipulated to

 (4.16) ^xn=argmin^v∈Rns12∥∥MΦ(^v−^xn−1)−ΔtMΦ(ZTMΦ)†ZTf(x0+Φ^xn,t;μ)∥∥22.

Note that the terms in norm in Eq. (4.16) is very similar to the discretized DEIM residual before Galerkin projection (i.e., applying the backward Euler time integrator to Eq. (4.4) gives the terms in norm in Eq. (4.16)). They are only different by the fact that one is an inverse and the other one is the Moore–Penrose inverse. In fact, if , then Eq. (4.16) is equivalent to applying the backward Euler time integrator to Eq. (4.4) and minimize the norm of the corresponding residual.

For the special case of being an identity matrix, Eq. (4.16) becomes

 (4.17) ^xn=argmin^v∈Rns12∥∥^v−^