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  takes a week to finish with cores. Other computationally expensive simulations include the 3D shocked spherical Helium bubble simulation appeared in  and the inertial confinement fusion implosion dynamics simulations appeared in . 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, 18], the Discrete Empirical Interpolation Method (DEIM) [11, 14]
and the Gauss–Newton with Approximated Tensors[9, 10], the Best Point Interpolation Method (BPIM) 
, the Missing Point Estimation (MPE), and Cubature Methods [3, farhat2014dimensional, 17, 19]. There a hyper-reduction (The term first used in the context of ROMs by Ryckelynck in 
) is necessary to efficiently reduce the complexity of nonlinear terms for a considerable speedup 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. 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, farhat2014dimensional, 17] 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  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)  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) , 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. 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
no additional collection of nonlinear term snapshots (i.e., storage saving).
no additional compression of snapshots (e.g., no additional singular value decomposition of nonlinear term snapshots, implying computational cost saving).
1.1 Organization of the paper
Notations are introduced in Section 2. We start our discussion by describing the time-continuous representation of the FOM in Section 3. Section 3 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 4 for the SNS method to be applied. The SNS method is described in Section 5 and applied to those ROMs. Section 5 also introduces the conforming subspace condition to justify the SNS method. 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.
The cardinality of a subspace, , is denoted as . The Kronecker product of two matrices and is denoted by and is defined as
3 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
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 (3.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. (3.1), by time-marching with the following update:
Eq. (3.2) implies the following subspace inclusion:
By induction, we conclude the following subspace inclusion relation:
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
The implicit Backward Euler (BE) method numerically solves Eq. (3.1), by solving the following nonlinear system of equations for at -th time step:
Eq. (3.6) implies the following subspace inclusion:
By induction, we conclude the following subspace inclusion relation:
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
4 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.
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
where denotes a basis matrix and with denotes the generalized coordinates. Replacing with in Eq. (3.1) leads to the following system of equations with reduced number of unknowns:
For constructing , Proper Orthogonal Decomposition (POD) is commonly used. POD  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 and Karhunen–Loève expansion  in stochastic analysis. POD forms a solution snapshot matrix, and computes its thin SVD:
where and are orthonormal 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. 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 [20, 23].
Note that Eq. (4.2) has more equations than unknowns (i.e., an overdetermined system). It is likely that there is no solution satisfying Eq. (4.2). In order to close the system, the Galerkin projection solves the following reduced system:
Applying a time integrator to Eq. (4.4) leads to a 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 speedup 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
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:
where is the sampling matrix and is the
th column of the identity matrix. Therefore, Eq. (4.5) becomes
The DEIM approach does not construct the sampling matrix .
Instead, it maintains the sampling indices
222The sampling indices can be found
either by a row pivoted LU decomposition  or a
column pivoted QR decomposition
or a column pivoted QR decomposition and corresponding rows of and . This enables DEIM to achieve a speedup when it is applied to nonlinear problems.
The original DEIM paper  constructs the nonlinear term basis by applying another POD on the nonlinear term snapshots333the 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 5 discusses how to avoid the collection of nonlinear term snapshots and an extra SVD without losing the quality of the hyper-reduction.
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. (3.1) in time domain and replaces with for in residual functions defined in Section 3 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 . The residual functions for implicit time integrators are defined in (3.9), (A.8), and (A.12) for various time integrators. For example, the residual function with the backward Euler time integrator444Although 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
The basis matrix can be found by the POD as in the DEIM approach. Note that Eq. (4.8) 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:
The Gauss–Newton method with the starting point is applied to solve the minimization problem (4.9) in GNAT. However, as in the DEIM approach, a hyper-reduction is required for a speedup due to the presence of the nonlinear residual vector function. The GNAT method approximates the nonlinear residual term with gappy POD , whose procedure is similar to DEIM. The GNAT method approximates the nonlinear residual term as
where , , denotes the residual basis matrix and denotes the generalized coordinates of the nonlinear residual term. In constrast to DEIM, the GNAT method solves the following least-squares problem to obtain the generalized coordinates at n-th time step:
where , , is the sampling matrix and is the th column of the identity matrix . The solution to Eq. (4.11) is given as
where the Moore–Penrose inverse of a matrix is defined as . Therefore, Eq. (4.10) becomes
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 speedup when it is applied to nonlinear problems.
Finally, the GNAT method minimizes the following least-squares problem at every time step, for example, with the backward Euler time integrator:
The GNAT method applies another POD to a nonlinear residual term snapshots to construct . The original GNAT paper  collects residual snapshots at every Newton iteration from the LSPG simulations555We denote the LSPG simulation to be the ROM simulation without hyper-reduction. That is, LSPG solves Eq. (4.9) without any hyper-reduction if the backward Euler time integrator is used. for several reasons:
The GNAT method takes LSPG as a reference model (i.e., denoted as Model II in ). Its ultimate goal is to achieve the accuracy of Model II.
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.
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:
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).
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.
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 5 discusses how to avoid the collection of nonlinear term snapshots and the extra SVD without losing the quality of the hyper-reduction.
4.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. (3.1) in a space–time form. For example, if the backward Euler time integrator666 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. (3.6) becomes
Note that denotes the parameterized space–time state implicitly defined as the solution to the problem (4.15) 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. (4.15) is defined as
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’
where with and whose elements are all one. Further, denotes a space–time basis vector. Thus, the ST-LSPG method approximates the numerical solution as
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. (4.17) with gives the residual vector function
where and denotes a space–time basis matrix that is defined as and denotes the generalized coordinate vector that is defined as
Note that vanishes. Eq. (4.20) becomes
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 :
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 [30, 13]. 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:
The ST-LSPG ROM solves Eq. (4.23) without any hyper-reduction. As in the DEIM and GNAT approaches, a hyper-reduction is required for a considerable speedup 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 , 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
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, :
where , , is the sampling matrix and is the th column of the identity matrix . The solution to Eq. (4.25) is given by
Therefore, Eq. (4.24) becomes
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 speedup when it is applied to nonlinear problems.
Finally, the ST-GNAT solves the following least-squares problem at every time step, for example, with the backward Euler time integrator:
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 ). Below is a list of the approaches introduced in  and explains advantages and disadvantages of each:
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. Furthermore, it requires the extra SVD on the residual snapshots, which is not necessary for our proposed method.
Projection of FOM training solutions. This approach takes the following steps:
take FOM state solution at every Newton iteration.
re-arrange them in the space–time form (i.e., in Eq. (4.16)).777
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.
project them onto the space–time subspace, (i.e., ).
compute the corresponding space–time residual (e.g., in Eq. (4.17) in the case of the backward Euler time integrator).
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.
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. (4.17) 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.
5 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. (3.4) and (3.8) with the forward and backward Euler time integrators.888Other subspace inclusion relations are shown in Appendix A.
Eq. (5.1) is an over-determined system, so it is likely that it will not have a solution. However, if there is a solution, then a necessary condition for Eq. (5.1) to have a non-trivial solution is that
This condition says that the intersection of and needs to be non-trivial if there is a non-trivial solution to Eq. (5.1). Typically, we build first, using the POD approach explained in Section 4.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. (5.1). Given , the largest subspace intersection it can be is , i.e.,
We call this condition as the conforming subspace condition. The conforming subspace condition leads to two obvious choices for :
Assuming is invertible, Eq. (5.5) becomes:
For the special case of being an identity matrix, Eq. (5.6) becomes:
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. (5.1) becomes
(5.8) (5.9) (5.10)
Assuming is invertible, Eq. (5.9) becomes:
For the special case of being an identity matrix, Eq. (5.11) becomes
The DEIM-SNS approach solves either Eq. (5.6) or Eq. (5.11) depending on the choice of above. Applying a time integrator to Eq. (5.6) or Eq. (5.12) leads to a reduced OE, whose solution, , can be lifted to find the full order approximate solution via .
Additionally, the subspace inclusion relations999Eq. (3.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, 4, 5, and 6 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.
The GNAT method needs to build a nonlinear residual term basis, , as described in Section 4.2. The nonlinear residual term is nothing more than linear combinations of the nonlinear term and the time derivative approximation as in Eq. (4.8). 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:
which can be manipulated to
Note that the terms in norm in Eq. (5.14) is very similar to the discretized DEIM residual before Galerkin projection (i.e., applying the backward Euler time integrator to Eq. (5.4) gives the terms in norm in Eq. (5.14)). 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. (5.14) is equivalent to applying the backward Euler time integrator to Eq. (5.4) and minimize the norm of the corresponding residual.
For the special case of being an identity matrix, Eq. (5.14) becomes
The GNAT-SNS method solves either Eq. (5.14) or Eq. (5.16) depending on the choice of above. For the special case of being an identity matrix, the GNAT-SNS method solves either Eq. (5.15) or Eq. (5.17) depending on the choice of . The reduced solution can be lifted to find the full order approximate solution via .
The ST-GNAT method needs to build a space–time nonlinear residual term basis, , as described in Section 4.3. We are going back to Eq. (4.22) to find the conforming subspace condition for ST-GNAT-SNS. In order to increase the chance of making the space–time residual function defined in Eq. (4.22) zero, the following conforming subspace condition can be made:
Therefore, we propose the following bases of the space–time nonlinear residual term:
The first choice is to set . Thus, for example, the ST-GNAT-SNS method solves the following least-squares problem with the backward Euler time integrator:
which can be rewritten as
This is what the ST-GNAT-SNS method solves if .
The second choice is to set where with and . The obvious choice for is to take a larger truncation of the factor matrices from the tensor decomposition (e.g., ST-HOSVD and LL1) of the solution snapshot tensor than . Note that this particular choice of results in the first columns of being the same as (i.e., where with ). In this case, the ST-GNAT-SNS method solves the following least-squares problem:
In addition to the choices above, we propose the following two choices for the special case of being an identity matrix:
The first choice is to set . Thus, for example, the ST-GNAT-SNS method solves the following least-squares problem with the backward Euler time integrator:
The second choice is to set . In this case, the ST-GNAT-SNS method solves the following least-squares problem: