{'fraction_non_alphanumeric': 0.041762898612141115, 'fraction_numerical': 0.01339124907299502, 'mean_word_length': 4.369894185914211, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
S-OPT: A Points Selection Algorithm for Hyper-Reduction in Reduced Order Models

While projection-based reduced order models can reduce the dimension of full order solutions, the resulting reduced models may still contain terms that scale with the full order dimension. Hyper-reduction techniques are sampling-based methods that further reduce this computational complexity by approximating such terms with a much smaller dimension. The goal of this work is to introduce a points selection algorithm developed by Shin and Xiu [SIAM J. Sci. Comput., 38 (2016), pp. A385--A411], as a hyper-reduction method. The selection algorithm is originally proposed as a stochastic collocation method for uncertainty quantification. Since the algorithm aims at maximizing a quantity S that measures both the column orthogonality and the determinant, we refer to the algorithm as S-OPT. Numerical examples are provided to demonstrate the performance of S-OPT and to compare its performance with an over-sampled Discrete Empirical Interpolation (DEIM) algorithm. We found that using the S-OPT algorithm is shown to predict the full order solutions with higher accuracy for a given number of indices.

1. Introduction. Physical simulation is the key to the developments of science, engineering and technology. Various physical processes are mathematically modeled by time-dependent nonlinear partial differential equations (PDEs). Since analytical solutions to such problems are not available in general, one has to resort to numerical methods to effectively approximate the solution. State-of-the-art numerical methods have been proven successful in obtaining accurate approximations of the groundtruth observations in various application problems. However, subject to the complexity and the scale of the problem domain, the computational cost of such numerical methods could be prohibitively high. Even with high-performance computing, a single forward simulation could take a very long time. Yet, multiple forward simulations are typically required in some real world decision-making applications such as design optimization [60,23,22,62], optimal control [19,13], uncertainty quantification [57,5], and inverse problems [31,5], which make such problems computationally intractable.
Constructing a reduced order model (ROM) is a popular and powerful computational technique to obtain sufficiently accurate numerical solutions with considerable speed-up compared to the corresponding full order model (FOM). Various model reduction schemes have been proposed. Many of them seek to extract an intrinsic solution manifold using a condensed solution representation. Depending on how these representations are constructed, two major approaches exist-linear subspace reduced order models (LS-ROM) and nonlinear manifold reduced order models (NM-ROM). In either case, the governing equations are projected onto the solution manifold as part of the reduction strategy, and therefore these approaches are referred to as projectionbased reduced order models (PROMs). This differs from other ROM techniques such as interpolation or data fitting.
In LS-ROM, the reduced basis vectors are obtained through, for example, proper orthogonal decomposition (POD). The number of degrees of freedom is then reduced by substituting the ROM solution representation into the (semi-)discretized governing equation. This approach takes advantage of both the known governing equations and the solution data generated from the corresponding FOM simulations to form LS-ROM. Example applications include, but are not limited to, nonlinear diffusion equations [38,29], Burgers equation and the Euler equations in small-scale [17,18,8] and large-scale, convection-diffusion equations [48,43], the Navier-Stokes equations [63,6], the compressible Euler equations in a moving Lagrangian frame [21,12], rocket nozzle shape design [2], flutter avoidance wing shape optimization [15], topology optimization of wind turbine blades [20], lattice structure design [47], porous media flow/reservoir simulations [32,40,65,59], computational electro-cardiology [64], inverse problems [30], shallow water equations [66,58], Boltzmann transport problems [16], computing electromyography [49], spatio-temporal dynamics of a predator-prey system [24], acoustic wave-driven microfluidic biochips [3], and the Schrödinger equation [11]. However, in advection-dominated problems, the intrinsic solution space cannot be approximated by subspaces with a small dimension, i.e., the solution space with slowly decaying Kolmogorov n-width. As an alternative to LS-ROM, we can replace the linear subspace solution representation with a nonlinear manifold. This type of ROM is known as NM-ROM. A neural network-based reduced order model is developed in [45] and extended to preserve the conserved quantities in the physical conservation laws [44]. Recently, Kim, et al., [42,41] have achieved a considerable speed-up with NM-ROMs via autoencoder.
The main advantage of a ROM (either LS-ROM or NM-ROM) is to reduce the computational cost by using a low-dimensional structure for representation of state variables. However, in nonlinear systems of PDEs, the performance of ROM is degraded due to the bottleneck issue of lifting to FOM size. That is, the nonlinear terms need to be evaluated in every time step as the state variables evolve in the time marching process. Since such evaluation scales with the FOM size, we cannot expect any speed-up without special treatment even if reduced representation is used to approximate the state variables. To overcome this issue, a hyper-reduction technique [53] is used to efficiently evaluate the nonlinear source terms by approximation. We note that a majority of the aforementioned literature achieved a true speed-up by applying a hyper-reduction technique. The key idea is to approximate the nonlinear terms using a small number of basis vectors, while keeping the number of evaluations of nonlinear terms as small as possible. To this end, the hyper-reduction technique requires one to strategically select a set of indices that leads to accurate approximations to the nonlinear terms.
One of the most well-known selection algorithms is the Discrete Empirical Interpolation Method (DEIM) [10], which aims at minimizing the operator norm from the error estimate (Lemma 3.2 of [10]). DEIM is implemented through a greedy algorithm that sequentially selects one index at a time with respect to a certain criterion. Carlberg, et al., [9] and [7] extend this idea to allow oversampling. Q-DEIM is introduced in [26] as a new framework for constructing the DEIM-related operator via the QR factorization with column pivoting. The stability and oversampling of DEIM is also investigated in [51].
The goal of this paper is to introduce the S-OPT sampling method [55] as a hyper-reduction technique in ROMs, and compare its performance with those by the DEIM family algorithms. S-OPT was first developed by Shin and Xiu in [55] as a points selection algorithm for least-squares based stochastic-collocation methods in Uncertainty Quantification (UQ) [56,50,35,34]. In the context of UQ, the goal is to find the best subset of points that yields the most accurate least-squares solution, which is closely aligned with the goal of the hyper-reduction. S-OPT aims to find a set of indices (rows or points) that maximizes the S quantity [1] (to be introduced in Section 3.1.2). The S quantity measures both the mutual column orthogonality and the determinant. The S-OPT algorithm is fundamentally different from the DEIM algorithm. The core principle of DEIM lies at maximizing the smallest singular value (spectral norm) of the underlying projection matrix, while S-OPT seeks to maximize both the product of all the singular values (determinant) and the column orthogonality of the underlying projection matrix. We employ the S-OPT algorithm as an index selection operator for hyper-reduction.
1.1. Paper organization. The FOM is described in Section 2.1 to introduce some background information and notation. We then describe the PROM formulation in Section 2.2, which leads to the description of the hyper-reduction procedure in Section 3. The sampling algorithms are part of the hyper-reduction procedure, so they are described in Section 3.1. Specifically, a DEIM algorithm with oversampling is outlined in Section 3.1.1, and the S-OPT algorithm is described in Section 3.1.2. The performance and comparison of the two algorithms are presented in Section 4 using four examples: 1D Burgers problem, 2D laminar viscous flow around airfoil, and two hydrodynamics examples, i.e., a 2D Gresho vortex problem and a 3D Sedov blast problem. The paper is concluded with summary and discussion in Section 5.
2. Problem Formulation. We start by defining the FOM and some notation used throughout the paper. We use two different PROMs in the example problems. First, we outline the LS-ROM formulation, followed by NM-ROM formulation for a general ordinary differential equation (ODE). This section precedes the main contribution of this paper, which is a sampling algorithm used with both of the PROM formulations.
For the rest of the paper, · is understood as either the standard Euclidean norm or the spectral matrix norm.
2.1. Full Order Model. Consider a system of nonlinear ODEs resulting from the semidiscretization of a system of PDEs in the space domain where t ∈ [0, T ] denotes time, u(t; µ) ∈ R N denotes the state vector of dimension N, u 0 (µ) ∈ R N denotes the initial condition, µ ∈ D denotes a vector of parameters defining the operating point of interest within the parameter domain D ⊆ R Nµ , A(µ) ∈ R N×N denotes a nonsingular matrix, and f : R N × R × D → R N is a nonlinear function and boundary conditions. The parameter and time dependence of A and variables are dropped in the rest of the paper for notational simplicity, and are implied. Furthermore, the dot notation denotes the derivative with respect to time.
The FOM system in Eq. (2.1) can be written in a residual form as follows: The time derivative term above can be approximated by various time integration schemes. Suppose the temporal domain [0, T ] is partitioned by {t n } Nt n=0 , where N t is the number of subintervals, t n denotes a discrete moment in time with t 0 = 0, t Nt = T , and t n−1 < t n for n ∈ {1, . . . , N t }. Throughout the paper, we use a superscript n to denote the time-discrete counterpart of a function evaluated at t = t n . For the numerical experiments, we use the implicit Backward Euler (BE), method and the second order explicit Runge-Kutta average (RK2-average) method to numerically solve Eq. (2.1), but other numerical time integration schemes are also applicable. For example, the BE method solves for u n at the n-th time step in Eq. (2.3): where f n = f (u n , t n , µ) and t n is the n-th time. The residual function with the BE time integrator is then defined by Although we continue the discussion using the BE time integrator and its residual as an example of the demonstration, the residuals of other types of time integrators can replace r n BE in a similar fashion. For example, we refer to [21] for the residual of the RK2-average method.

Projection-Based Reduced Order Model.
A PROM formulation relies on the concept that full state solutions can be represented in lower-dimensional manifold. As such, a PROM projects the governing equations to a manifold, resulting in lower-dimensional equations. If the manifolds for the solution field and the equations are the same, then the projection is called Galerkin. On the other hand, if the manifolds for the solution field and the equations are different, then the projection is called Petrov-Galerkin. Both the Galerkin and Petrov-Galerkin projection methods are considered. We also consider a PROM with a linear subspace solution representation (LS-ROM), as well as a model with nonlinear manifold solution representation (NM-ROM).

Linear Subspace Reduced Order Model.
A LS-ROM reduces the spatial dimension by approximating the full solution using a subspace W := span{φ i ∈ R N : i = 1, . . . , k} with dim(W) = k N, also called a trial subspace. The approximation u of the full solution is where u ref ∈ R N is a reference state, Φ ∈ R N×k denotes a basis matrix whose i-th column is φ i , and y ∈ R k is a vector of unknown generalized coordinates. The reduced subspace, Φ, is commonly found using POD [4], which is related to principal component analysis (PCA) in statistical analysis and the Karhunen-Loève expansion in stochastic analysis [39,46]. While we use POD to define the basis in this paper, the basis can be built using other options, such as Fourier modes. In POD, a set of basis functions are built by performing a singular value decomposition (SVD) over a solution snapshot matrix. These snapshots are based on solutions of the FOM, either steady-state solutions for multiple parameters, or time-varying solutions. Additionally, the reference state, u ref , can be found by taking the average of the collected snapshots.
Substitutingũ of Eq. . Note that r(y, v, t, µ) is linear in v. By projecting the system of equations onto a test subspace with basis matrix Ψ ∈ R N×k (which can be different than Φ) whose i-th column is ψ i , that is, we solve for v, which gives the governing equation for y Equation (2.7) corresponds to Galerkin projection when the test subspace is the same as the trial subspace, i.e. Ψ = Φ. When the test subspace differs, we have Petrov-Galerkin projection. However, Petrov-Galerkin projection is generally applied after discretizing in time, which leads us to describe the nonlinear Least-Squares Petrov-Galerkin (LSPG) projection procedure. A LSPG ROM substitutes u n = u ref + Φy n for u into Eq. (2.4) and minimizes the residual at each time instance. Using the BE time discretization as an example, we have where the LS-ROM backward Euler reduced residual is defined as The necessary first-order optimality condition for Eq. (2.8) is Here J f (·,tn,µ) ∈ R N×N is the Jacobian of f (·, t n , µ). This shows that Eq. (2.8) corresponds to a Petrov-Galerkin projection of the FOM equations with a trial subspace Φ and a test subspace Ψ = J n Φ, hence the name LSPG projection. A note on notation: throughout this paper, the hat (r) is used for reduced variables that lie in a smaller subspace such as R k , and the tilde (ũ) is used for variables in R N that approximate their non-accented counterparts.

Nonlinear Manifold Reduced Order
Model. The LS-ROM relies on a linear subspace (Φ) for the solution manifold. In this section, we outline a ROM with nonlinear solution representation. As a generalization of the linear subspace representation in (2.5), NM-ROM seeks approximation of the full state solutions in a trial manifold by where u ∈ R N denotes an approximation of the full solution, u ref ∈ R N denotes a reference state, g : R k → R N denotes a nonlinear function, and y ∈ R k denotes a vector of unknown latent variables. The nonlinear function is found by neural network training as a decoder in an autoencoder architecture. In this work, we adopt the shallow masked autoencoder as in [42], where the nonlinear function g is a scaled decoder with one single hidden layer and a sparsity mask in the output layer. The trainable parameters, i.e. the weight and bias in the decoder and the encoder networks, are optimized against the mismatch between the original training data and the corresponding autoencoder output. The trained decoder and its Jacobian, J g ∈ R N×k , is used to formulate the system of equations. More precisely, the counterpart of (2.6) is The NM-ROM Galerkin projection is then given bẏ Similarly, the NM-ROM LSPG projection is given by where the NM-ROM backward Euler reduced residual is defined as A note on time integrators: although the BE method is used for the discussion of LS-ROM and NM-ROM, other numerical time integrators can be used in a similar fashion.

Hyper-Reduction and Sampling.
Due to the term f still being nonlinear in the reduced subspace, the requirement to compute the Jacobian, and to use the Jacobian in a matrix-matrix multiplication for LSPG, the reduced, low-dimensional equations may still not be computationally more efficient than the FOM. However, the ROM framework allows for a further approximation of the equations using a hyper-reduction strategy.
Specifically, when Galerkin projection is utilized, the (Ψ AΦ) −1 Ψ part of Eq. (2.7) is time-independent, and can be computed offline for each parameter µ. This leaves the nonlinear function f ( u, t, µ) to be calculated at every iteration. However, f still scales with the full dimension. Furthermore, for Petrov-Galerkin projection, the projection matrix Ψ must be calculated at every iteration because the Jacobian relies on the current state. In the case of NM-ROM, both projection methods rely on the Jacobian. The full dimension dependency in the Jacobian and the nonlinear function limits the speed-up performance of the PROM.
There are two ways hyper-reduction can be applied. Firstly, the nonlinear term is approximated so that it scales with the reduced dimension. This is applicable when the projection matrix can be computed offline, so approximating only the nonlinear term results in sufficient computational savings. Secondly, the residual is approximated instead of independently working with the nonlinear term. This is especially favorable in cases where the nonlinear term may be difficult to access in the FOM solver, such as for residual minimization solvers. Either way, the aforementioned terms are approximated by using only a carefully constructed index subset and omitting gaps for ignored entries, so this procedure is called gappy tensor approximation, as first introduced in [28].
Hyper-reduction on f . The nonlinear term f ∈ R N can be approximated in a least-squares sense from a gappy form using a reduced basis Φ f ∈ R N×n f built for the nonlinear term and the set of sampled indices, Z = {i 1 , . . . , i n }. Given a gappy tensor f ∈ R n f , the approximation f of f is: , · · · , e in ] ∈ R N×n is a sampling matrix which contains n f ≤ n N columns of the N × N identity matrix, and e i is the standard basis in R N . In other words, Z samples f at only n indices. Choosing which indices to select is the subject of Section 3.1. The (minimum norm) solution to Eq.
where the superscript † denotes the Moore-Penrose pseudoinverse. A hyper-reduction of the nonlinear term f is then given as which is termed as "gappy POD" hyper-reduction. We remark that, The procedure that combines Galerkin projection and gappy POD for the nonlinear term is known as the discrete empirical interpolation method (DEIM). In the original DEIM paper [10], Φ f was computed by collecting snapshots of the nonlinear term and computing the POD. However, it has also been shown that solution snapshots can be used through the subspace relation in [18], i.e., Φ f = AΦ.
One can consider a simple approximation f = ZZ f to f , which only requires Z. This is termed "collocation." Even in this case, however, Φ f may still be used to build Z during the offline phase, yet is not used during the online phase of the ROM. This is justified by the assumption that f is well represented by the subspace spanned by Φ f .
Hyper-reduction on the residual term. Given v, t n , y n−1 and µ, the nonlinear term r := r n BE (v; y n−1 , µ) (defined in Eq. (2.9)) is similarly approximated by a gappy tensor r ∈ R nr using a reduced basis Φ r ∈ R N×nr built for the residual and a set of sampled indices, Z = {i 1 , . . . , i n }: which is the gappy-POD approach to hyper-reduction. Similarly as before, one can also use the collocation approximation r = ZZ r. The procedure that combines LSPG projection and gappy POD hyper-reduction for the residual term is called the Gauss-Newton with approximated tensors (GNAT) procedure [7]. Error analysis of hyper-reduction approximation. In practice, the sampling matrix Z is not built; rather, the selected indices are maintained along with the corresponding rows of Φ f , Φ r , f and r. Regardless of which approach is used, the goal is to approximate a vector b of size N (either b = f or r) using a predefined basis . Hence, we are faced with the following optimization problem: Find the optimal sampling matrix Z * ∈ R N ×n such that In Theorem 3.1, we quantify and estimate the oblique projection error b −b(Z) , which becomes the theoretical basis of many existing sampling methods.
Proof. Let a * =ã(I). Observe that the quantities A M , Q M , y M , S m ,β M ,β m stated in Theorem 3.1 of [55] are M , Q, b, Z , a * , a Z , respectively, in the notation of the current paper. It then follows from Theorem 3.1 of [55] which gives (3.6). Furthermore, it follows from (Z) = (Z Q) † Z proj M ⊥ b and the above equality that Since Q(Z Q) † Z is a projection operator, the proof is completed by observing Q = Z = 1.
Remark 3.2. Theorem 3.1 shows that the optimization problem (3.5) is equivalent to as proj M ⊥ b is independent of the sampling matrix Z. Yet, since the optimal sampling matrix Z * requires the full knowledge of b, the true optimality is hardly achieved in practice. 3.1. Sampling Algorithms. One of the most popular algorithms for the construction of the sampling matrix Z is based on the error bound of (3.7). These algorithms aim at constructing Z that makes (Z Q) † as small as possible. This may be viewed as E-optimality [52] in the optimal design community, which maximizes the smallest nonzero eigenvalue of (Z Q) Z Q. The commonly used algorithms in the PROM community, such as DEIM [10], Q-DEIM [26], follow this principle.
In the UQ community, a similar problem has been addressed in the context of least squares based stochastic collocation [55,56,50,35,34]. In particular, [55] proposed a method based on the error equality of (3.6). Since the true optimum is not available (as also mentioned in Remark 3.2), [55] developed the so-called S-optimality [1] that maximizes both the column orthogonality of Z Q and the determinant (Z Q) Z Q.
Maximizing the column orthogonality of Z Q can increase the numerical stability, which has been an issue for E-optimality based sampling methods, such as DEIM. In this paper, we refer to this method as "S-OPT".
The goal of this paper is to introduce the S-OPT sampling method to the PROM community as a hyper-reduction method, and compare its performance with those by the DEIM family algorithms.
In what follows, the following notations are used. Let Q be an orthogonal matrix obtained from a QR factorization of M if the columns of M are not orthonormal, and let Q = M otherwise. Given a set of indices Z = {i 1 , . . . , i } and a vector φ ∈ R N , the ith component of φ is denoted by φ(i).
3.1.1. Oversampled DEIM Algorithm. The oversampled DEIM algorithm used in this paper differs from the original DEIM algorithm [10] and more closely follows the sampling method from the LSPG paper (Algorithm 5 in [7]). Unlike the original DEIM, this allows for oversampling, i.e., selecting more samples than the number of POD modes, i.e., n > p. We call the algorithm "oversampled DEIM" in this paper, while it could be fairly named a greedy or gappy-POD algorithm, among other options. This algorithm has also been modified for space-time LSPG in [17].
By closely following notation of [61], we present the pseudo-code for the oversampled DEIM in Algorithm 3.1. Starting with an appropriate basis Q = [φ 1 , · · · , φ p ] as an input, as well as the desired number of samples n, the algorithm selects indices based on a greedy method. The first selection is simply the index of the largest absolute value entry of φ 1 . The algorithm then loops over the rest of the columns of Q, retaining all prior columns within Q 1:j = [φ 1 , · · · , φ j ]. Each column of Q is approximated using a gappy POD reconstruction (line 8), and the index at which there is the largest error is selected and included in Z. This loop continues until the sampled set contains the desired number of unique indices. If necessary, the set Z of total indices (sampled and required neighboring nodes) is also populated according to the numerical scheme.

S-OPT: A Quasi-Optimal Points Selection
Algorithm. This subsection introduces the sampling method (S-OPT) proposed in [55]. The underlying principle of S-OPT is to find a sampling matrix Z that makes (Z) as small as possible. According to Theorem 3.1, (Z) should satisfy If the sampling matrix Z is constructed to make Z Q preserve the same column orthogonality as Q, the right-hand side of (3.9) becomes zero, leading to (Z) = 0. However, this will only be the case when Z is the identity matrix with n = N and Construct Z = [e i1 , · · · , e i ] from Z 8: if |Z| = n then return 12: end for 13: end for in general, (Z) will be nonzero. With the goal of having a small (Z), the S-OPT [55] seeks to (a) maximize the column orthogonality of Z Q so that the right-hand side of (3.9) is minimized, and (b) maximize the determinant of (Z Q) Z Q so that the nonzero solution of (3.9) is small. A quantity denoted by S (the precise definition is given in (3.10)) that measures the mutual column orthogonality and the determinant was developed in [55]. The S-OPT aims at finding a sampling matrix Z that maximizes the quantity S(Z Q).
Since the true optimality Z * of (3.8) is not achievable in practice, the selection method based on the S-optimality (3.10) was termed "quasi-optimal" in [55] and "near-optimal/S-optimality" in [56]. Here, we simply refer to it as "S-OPT." Let A be an N × p matrix and α i be its ith column. Assuming α i = 0 for all i = 1, . . . , p, S(A) is defined to be It was shown in [55] that S(A) = 1 if and only if the columns of A are mutually orthonormal. Hence, maximizing S enforces both mutual column orthogonality and a larger determinant.
In the context of hyper-reduction, the S-OPT seeks the optimal sampling matrix Z S-OPT that maximizes S, i.e., Solving the above optimization problem requires one to compute multiple evaluations of S. The evaluation of S, however, can be expensive as it requires the computation of determinants. However, [55] presented an efficient way of evaluating S without computing determinants, based on a greedy algorithm and the Sherman-Morrison formula. The pseudo-code for the S-OPT index selection is presented in Algorithm 3.2 for the reader's convenience. We refer to the original paper [55] for full algorithmic details.
Algorithm 3.2 S-OPT selection algorithm [55]. Input: Desired number of sampled indices n ∈ {p, . . . , N }, and an orthonormal basis Z ← Z ∪ {i j } 6: end for 4. Main results. In this section we show four example problems using both sampling algorithms. They are a 1D Burgers equation problem, an aerodynamic laminar flow around airfoil, and two Lagrangian hydrodynamic problems.
The 1D Burgers equation compares the use of the sampling methods for reduced order models built with linear or nonlinear subspaces; that is, LS-ROM and NM-ROM. Meanwhile, the laminar airfoil problem builds a parametric LS-ROM using shape parameters. The airfoil problem also explores how truncating the POD basis may affect the results when using either sampling algorithm. The final two problems are the 2D Gresho vortex and the 3D Sedov blast Lagrangian hydrodynamics problems.
In these examples, we show that the sampling methods result in different error bounds, and we further examine the performance of the two sampling algorithms for varying the number of sampled indices.
Furthermore, the examples in this paper use different combinations of projection and hyper-reduction methods. Before presenting the problem descriptions and results, we present the specific methods used for each of the four examples here, for comparison and completeness purposes.
While the examples in Sections 4.1 and 4.2 both apply LSPG projection, the example in Section 4.1 applies gappy POD hyper-reduction (or the GNAT procedure), and the airfoil example in Section 4.2 applies collocation hyper-reduction for the residual term. Finally, the two examples in Section 4.3 use the DEIM approach; i.e., a Galerkin projection with a gappy POD approximation of the nonlinear term.
Applying LSPG with gappy POD hyper-reduction, as done in Section 4.1, gives which aims at minimizing the approximation (3.4) of r n BE . Applying LSPG with collocation hyper-reduction, as done in Section 4.2, gives In Equations (4.1) and (4.2), the fully discretized residual r n BE is used as an example, and a residual defined using a different discretization can be used in a similar fashion.
And lastly, applying Galerkin projection with gappy POD hyper-reduction of the nonlinear term, as done in Section 4.3, giveṡ In the case of Galerkin projection, remember that the projection matrix Ψ is taken to be the solution basis Φ. The front matter (Ψ AΦ) −1 Ψ Φ f (Z Φ f ) † can be precomputed once, and the nonlinear term Z f can be evaluated at only the selected indices. Equations (4.1, 4.2, 4.3) no longer scale with the full dimension N, but rather with the number of selected indices n and the reduced dimension k.
Finally, it is important to note that computing the residual at the selected nodes may depend on the solution at neighboring nodes. For example, for a finite volume discretization involving inviscid fluid dynamics (Euler) equations (e.g., Section 4.2) as well as the finite element method (e.g., Section 4.3), the solution at the neighboring nodes is required to compute the flux, and then the residual at the selected node can be computed using the flux. For this reason it is necessary to maintain the indices of all required neighboring nodes, in addition to the selected nodes.

1D Burgers Equation.
As the first example, we consider the 1D inviscid Burgers equation with the periodic boundary condition u(t, 0) = u(t, 2) for t ∈ [0, 0.5] and the initial condition u(0, x) = 1 + 1 2 (sin(2πx − π 2 ) + 1). We decompose the spatial domain into a uniform mesh with mesh size ∆x = 0.002, which consists of 1001 grid points x i = i∆x, i ∈ {0, 1, 2, . . . , 1000}, at which the discrete solution is defined by u i (t) = u(t, x i ). The semi-discretization of the Burgers equation is given by: where u = (u 1 , u 2 , . . . , u 1000 ) ∈ R 1000 is the coefficient vector, denotes the entrywise product, and To obtain a fully discrete scheme, we similarly decompose the temporal domain into N t = 500 subintervals to form a uniform mesh with mesh size ∆t = 0.001, and use a backward difference to numerically approximate the temporal derivative. Figure 1 shows the initial condition and the final-time solution for 1D Burgers equation. It can be seen that a shock is eventually developed at x = 1.5. In [42], a linear-subspace reduced order model (LS-ROM) and a nonlinear-manifold reduced order model (NM-ROM) are developed for (4.5), where DEIM is used for hyper-reduction in both approaches. In what follows, we compare the numerical results of reduced order models using the sampling algorithms in Section 3.1, i.e. the oversampled DEIM and S-OPT. The maximum-in-time L 2 relative error for the ROM solution u is measured against the corresponding FOM solution u, which is defined as: All the 1D Burgers equation simulations in this subsection use Lassen in Livermore Computing Center 1 , on Intel Power9 CPUs with 256 GB memory and NVIDIA V100 GPUs, peak TFLOPS of 23,047.20, and peak single CPU memory bandwidth of 170 GB/s.
In Figure 2, we depict the maximum-in-time relative discrete L 2 error of the solution against the number of sampling indices for both sampling algorithms. The number of sampling indices n s takes value within n f ≤ n s ≤ n f + 30, where n f = 30 in LS-ROM and n f = 47 in NM-ROM. In LS-ROM, despite only a minor difference in the error, the DEIM algorithm yields more oscillation in the error with respect to the number of sampling indices. In NM-ROM, the oscillation is even more severe in DEIM. In other words, the S-OPT algorithm yields more stable results in the Burgers equation with increasing number of sampling indices. Figure 3 and Figure 4 illustrate the selected nodes by the sampling algorithms in LS-ROM and NM-ROM respectively. In both cases, both S-OPT and DEIM select nodes around the expended shock wave, where the nonlinearity occurs, in priority. As the oversampling number increases, S-OPT tends to sample nodes in a more widespread manner, while DEIM still tends to densely select nodes close to the shock, as shown in Figure 1. 4.2. 2D Laminar Airfoil. The next example considers a 2D laminar airfoil using the steady compressible Navier-Stokes equations with a low Reynolds number (i.e., no turbulence). For the steady problem, we define the following residual definition, which is used in the construction of a LS-ROM: where F c is the convective flux, and F v is the viscous flux. For this 2D case, the conservative variables are given by u = {ρ, ρv, ρE} with the velocity vector defined   as v = {v 1 , v 2 } ∈ R 2 , ρ denoting the fluid density, and E denoting the total energy per unit mass. Furthermore, F c and F v are defined as: where p is the pressure field, q is the heat flux, and τ is the viscous stress tensor defined as: with µ as viscosity. The airfoil surface is set as an adiabatic wall (zero heat flux) and the far field is set to the free stream conditions. The FOM is solved using a pseudo-time stepping method, specifically an implicit, local time stepping method, to march the solution forward from the initial condition (freestream fluid state) to the steady state solution. The local time stepping technique allows for quicker convergence by allowing the time step to vary between elements. This approach is typical for steady aerodynamics. We use converged steady solutions as simulation data to construct the reduced basis.
The full order model domain contains 8741 mesh points with a NACA0012 airfoil surface in the center. We use shape parameters for a parameterized ROM study. In this case, the airfoil shape is slightly altered using three Hicks-Henne bumps [37] on the upper surface and three on the bottom surface. Hicks-Henne bumps are smooth bump formulations, and each shape parameter defines the amplitude of the Hicks-Henne bump. The six-dimensional parameter space is sampled using 73 Latin hypercube samples, resulting in 73 total snapshots. Both the FOM and ROM equations are solved using the open-source CFD code SU2 [27], aided by the libROM code [14] from Lawrence Livermore National Lab for the snapshot collection and POD computation.
The laminar airfoil problem is a good academic problem for quick testing, because it does not require a large amount of computational resources. The 73 snapshots are collected on the Sherlock cluster operated by the Stanford Research Computing Center. Then the reduced order model simulations are run serially on a macOS laptop with 2.4 GHz Intel Core i9 processor and 16 GB memory.
The POD basis has a dimension k = 73, but is then truncated to 20, 10, and 5 POD modes. Figure 5 shows the first four modes, or singular vectors, for the density field from the POD computation. The variations of the flow field due to the airfoil shape differences are captured by the POD modes, so the flow solution for any set of shape parameters can be approximated using a combination of these modes. Both sampling methods are tested and the results are shown in Figure 6. In all cases, the goal of the ROM is to predict the baseline NACA0012 airfoil shape, which is not included in the 73 snapshots.
It is also important to note that for finite-volume CFD models there are multiple equations at each node, so the input basis for the sampling algorithms needs to be condensed so that the size of the first dimension and the number of nodes are the same. We do this by taking the norm of the basis values at each node, resulting in one basis value per node. This is in contrast to Section 4.3, where we also show a different scenario of building a separate basis for a different field, and therefore, merging the sampling indices from different fields.
The measurements of error shown in Figure 6 are defined as follows for the L 2 error and the maximum relative error: where the truth values are obtained from the full order model solution of the prediction case (the NACA0012 baseline shape) using SU2.  In Figure 6 we see that as the number of sampled indices decreases, the S-OPT algorithm outperforms the DEIM algorithm as shown by the error measures. At the most extreme level of hyper-reduction, using only 437 samples, the ROM built with S-OPT algorithm and 5 POD modes performs better than the ROM built with DEIM and 20 POD modes.
However, there is a cost associated with the performance improvement. The ROMs built with the S-OPT algorithm generally take longer to converge. Figure 7 shows one reason the S-OPT algorithm ROMs take longer. Since this is a viscous problem, two levels of node-neighbors are required to compute the residual at each selected node. The S-OPT algorithm tends to choose nodes that are more spread out in the domain, while the DEIM algorithm chooses nodes near the airfoil surface where most of the changes from solution to solution take place. Since the set of neighboring nodes is greater for the S-OPT ROM, the total dimension for the hyper-reduced ROM is larger and each iteration takes longer. It is important to note that this dimension still scales with the reduced dimension, k.

Lagrangian hydrodynamics.
In the next two examples, we consider advectiondominated problems arising in compressible gas dynamics. The Euler equation is used to model the high-speed flow and shock wave propagation in a complex multimaterial setting, and numerically solved in a moving Lagrangian frame [36], assuming no external body force is exerted: Here, d dt = ∂ ∂t +v ·∇ is the material derivative, ρ denotes the density of the fluid, x and v denote the position and the velocity of the particles in a deformable medium Ω(t) in the Eulerian coordinates, σ denotes the deformation stress tensor, and e denotes the internal energy per unit mass. In gas dynamics, the stress tensor is isotropic, and we write σ = −pI +σ a , where p denotes the thermodynamic pressure, and σ a denotes the artificial viscosity stress. The thermodynamic pressure is described by the equation of state, and can be expressed as a function of the density and the internal energy.
With the assumption of polytropic ideal gas with an adiabatic index γ > 1, which yields the equation of state (4.13) p = (γ − 1)ρe.
The system is prescribed with an initial condition and a boundary condition v · n = g, where n is the outward normal unit vector on the domain boundary. Using a high-order curvilinear finite element method (FEM) for Lagrangian hydrodynamics [25] for semi-discretization of the Euler equation results in the differential system: (4.14) where v, e, x : [0, T ] → R denotes the FEM coefficient vector functions for velocity v, internal energy density e and position x. In order to obtain a fully discretized system of equations, we apply RK2-average scheme as the time integrator, a modification of the midpoint Runge-Kutta second-order scheme. The RK2-average scheme is written as x n+ 1 2 = x n + (∆t n /2)v n+ 1 2 , x n+1 = x n + ∆t nv where the state w n = (v n ; e n ; x n ) ∈ R N is used to compute the updates F n 1 = (F (w n )) · 1, F n tv = (F (w n )) · v n+ 1 2 , (4. 18) in the first stage. Similarly, w n+ 1 2 = (v n+ 1 2 ; e n+ 1 2 ; x n+ 1 2 ) ∈ R N is used to compute the updates withv n+ 1 2 = (v n + v n+1 )/2 in the second stage. Since explicit Runge-Kutta methods are used, we need to control the time step size in order to maintain the stability of the fully discrete schemes. We follow the automatic time step control algorithm described in Section 7.3 of [25].
We use a linear-subspace reduced order model technique for (4.14) as developed in [21]. The solution nonlinear subspace (SNS) procedure in [18] is used to build the nonlinear term bases Φ F1 and Φ Ftv on the right-hand-side of the momentum conservation equation and energy conservation equation separately, and the sampling indices for hyper-reduction of each of nonlinear terms are used to construct the sampling matrices Z F1 and Z Ftv , respectively. In the following benchmark experiments, we compare the numerical results of reduced order model using the sampling algorithms discussed in Section 3.1, i.e., the oversampled DEIM and S-OPT. With the automatic time step control algorithm, it is very likely that the temporal discretization used in the hyper-reduced system is different from the full order model even with the same problem setting. To this end, we denote by N t the number of time steps in the fully discrete hyper-reduced system, to differentiate it from the notation N t for the full order model. The L 2 relative error for each ROM field is measured against the corresponding FOM solution at the final time T , which is defined as: All the Lagrangian hydrodynamics simulations in this subsection use Quartz in Livermore Computing Center 2 , on Intel Xeon CPUs with 128 GB memory, peak TFLOPS of 3251.4, and peak single CPU memory bandwidth of 77 GB/s. 4.3.1. 2D Gresho vortex. The Gresho vortex problem is a two-dimensional standard benchmark test for the incompressible inviscid Navier-Stokes equations [33]. In this problem, we consider a manufactured smooth solution from extending the steady state Gresho vortex solution to the compressible Euler equations. For the detailed description of the set-up of the Gresho vortex problem, we refer the readers to Sections 6.1.1 and 6.2 of [21]. The final time is taken as T = 0.1. Figure 8 shows the initial condition and the final-time solution for 2D Gresho vortex. It can be seen that the vortex is rotating. We first investigate how the sampling algorithms affect the projection error of some sampled snapshots of the right hand side of the energy conservation equation. Following the notations in Theorem 3.1, with b = F n tv , M = Φ Ftv , and Z = Z Ftv , we examine the oblique projection error F n tv −F n tv (Z Ftv ) and use the orthogonal projection error proj ⊥ Φ F tv F n tv as reference. In Figure 9, we illustrate the effects of the choice of sampling algorithm with varying number of sampling indices on the oblique projection error in some nonlinear snapshot samples, which is a crucial component of error bounding in PROM for nonlinear problems.
In this test, the dimension of the nonlinear term subspaces, i.e. the number of columns in the nonlinear term bases, of the nonlinear term in momentum conservation equation and energy conservation equation is 32 and 72, respectively. In each of the nonlinear term evaluations, the number of sampling indices is taken as the product of the nonlinear term basis dimension and the oversampling ratio, which takes a value between 2 and 15. The oblique projection error in the S-OPT has a much faster and smoother decay to the orthogonal projection error than the oversampled DEIM selection for all the selected snapshot samples.
In Figure 10, we depict the final-time L 2 error of the reduced order model solution against the number of sampling indices for both sampling algorithms. When the oversampled DEIM algorithm with the oversampling ratio less than 5 is used, the reduced order model is unstable and is not able to yield meaningful approximation. Moreover, with the oversampling ratio between 5 and 11, the final-time error for DEIM is oscillatory and significantly larger than the S-OPT for all solution components in Lagrangian hydrodynamics. This suggests the advantage of both accuracy and stability by using the S-OPT over the oversampled DEIM selection. Figure 11 depicts the sample mesh, which consists of all the elements containing a sampling node, of both sampling algorithms at an oversampling ratio of 3. Similar to our previous observation, S-OPT tends to sample nodes in a more widespread manner. In this example, DEIM densely selects nodes close to the rotating vortex, as shown in Figure 8. 4.3.2. 3D Sedov blast. The Sedov blast problem is a three-dimensional standard shock hydrodynamic benchmark test [54], where we consider a delta source of internal energy initially deposited at the origin of a three-dimensional cube. The final  Again, we begin with investigating how the sampling algorithms affect the projection error of some sampled snapshots of the right hand side of the energy conservation equation by examining the oblique projection error F n tv −F n tv (Z Ftv ) , with the orthogonal projection error proj ⊥ Φ F tv F n tv as reference. In Figure 13, we illustrate the effects of the choice of sampling algorithm with varying number of sampling indices on the oblique projection error in some nonlinear snapshot samples, which, again, is a crucial component of error bound in PROM for nonlinear problems. In this test, the dimensions of the nonlinear term subspaces, i.e. the number of columns in the nonlinear term bases, for the momentum conservation equation and energy conservation equation are 53 and 13, respectively. In each of the nonlinear term evaluations, the number of sampling indices is taken as the product of the nonlinear term basis dimension and the oversampling ratio, which takes value between 2 and 15. Similar to the Gresho vortex problem, while the oblique projection error in both sampling algo-rithms asymptotically decays to the orthogonal projection error, it can be observed that the decay is much faster with the S-OPT selection than the oversampled DEIM selection for all the selected snapshot samples. In Figure 14, we depict the final-time L 2 error of the solution against the number of sampling indices for both sampling algorithms. Again, it can be seen that for all solution components in Lagrangian hydrodynamics, the S-OPT selection gives much more stable decay in the solution error than the oversampled DEIM selection.

5.
Conclusions. This work proposes the use of the S-OPT algorithm for selecting indices for hyper-reduction of projection-based reduced order models. The algorithm chooses indices while trying to keep the POD basis orthogonal to enhance the numerical stability, while other selection methods, such as DEIM, do not.
As shown in the results, the indices chosen by the two algorithms tend to lie in the regions experiencing the most change from snapshot to snapshot. However, for the S-OPT algorithm those indices also tend to be more spread out.
The resulting error for the ROM is smaller when using the S-OPT algorithm, especially when using a small number of total selected indices. It is expected that the relative error of the ROM will increase when selecting fewer indices, but the error when using the S-OPT algorithm appears to increase more smoothly, whereas selecting fewer indices using the oversampled DEIM algorithm has a greater adverse effect on the ROM performance. A topic of future research is comparing S-OPT to other selection algorithms, and investigating whether or not there is an indication of how many indices to select for a certain error bound.