Model order reduction with novel discrete empirical interpolation methods in space-time

This work proposes novel techniques for the efficient numerical simulation of parameterized, unsteady partial differential equations. Projection-based reduced order models (ROMs) such as the reduced basis method employ a (Petrov-)Galerkin projection onto a linear low-dimensional subspace. In unsteady applications, space-time reduced basis (ST-RB) methods have been developed to achieve a dimension reduction both in space and time, eliminating the computational burden of time marching schemes. However, nonaffine parameterizations dilute any computational speedup achievable by traditional ROMs. Computational efficiency can be recovered by linearizing the nonaffine operators via hyper-reduction, such as the empirical interpolation method in matrix form. In this work, we implement new hyper-reduction techniques explicitly tailored to deal with unsteady problems and embed them in a ST-RB framework. For each of the proposed methods, we develop a posteriori error bounds. We run numerical tests to compare the performance of the proposed ROMs against high-fidelity simulations, in which we combine the finite element method for space discretization on 3D geometries and the Backward Euler time integrator. In particular, we consider a heat equation and an unsteady Stokes equation. The numerical experiments demonstrate the accuracy and computational efficiency our methods retain with respect to the high-fidelity simulations.


INTRODUCTION
Conventional high-fidelity approximations of phenomena governed by partial differential equations (PDEs) generally involve fine space and time discretizations, requiring large computational resources.When the PDE is parameterized, the parameter-to-solution map is a high-dimensional manifold.Instead of performing the high fidelity (HF) simulation every time we want to get the solution for a given parameter, one can design a surrogate model for the parameter-to-solution map.Projection-based reduced order models (ROMs), such as the reduced basis reduced basis (RB) method, create surrogate models by approximating the solution manifold by a lower-dimensional vector space.These surrogates provide accurate solutions in a much shorter time and with much fewer computational resources.
Traditionally, these methods only aim to reduce the spatial complexity of the full order model (FOM) by employing (Petrov-)Galerkin projections on the ordinary differential equation (ODE) resulting from the spaceonly discretization, while using the same HF time marching scheme for the integration in time.Space-time ROMs can be applied to the fully discrete system of FOM equations to retain a reduction along the temporal and spatial dimensions.In particular, the space-time reduced basis (ST-RB) method initially builds a spatial and a temporal basis spanning a low-dimensional subspace, using either a greedy algorithm or a truncated proper orthogonal decomposition (TPOD) approach.Then, under a suitable norm, it minimizes the FOM residual on such subspace.
A comprehensive literature review of methods to reduce temporal complexity can be found in [12].The methods presented in [39,40,43,44] are characterized by a reduction of both spatial and temporal dimensions.However, they exhibit some relevant drawbacks, the most salient ones being the need for (uncommon) finite element (FE) discretizations of the time domain and the absence of hyper-reduction techniques to handle non-linearities efficiently.The authors of [11,22,37] successfully achieve a temporal complexity reduction by implementing ST-RB in the case of both 2D and 3D problems, demonstrating promising results both in terms of accuracy and computational efficiency.However, the applications discussed in these works are restricted to simple parameter-to-solution maps, i.e., affine dependence of operators on parameters.
In general, the performance of linear ROMs such as ST-RB depends on whether the FOM is reducible or not, which in turn is determined by two main factors: • There exists a low-dimensional spatio-temporal manifold that accurately approximates the manifold of the FOM solutions.• The implementation of the method admits an efficient splitting into an offline phase, (where all the expensive, parameter-independent operations are executed once and for all) and an online phase (where the inexpensive parameter-dependent computations are performed for every parameter of interest).The first condition is generally met when dealing with elliptic or dissipative time-dependent applications whose parameter-to-solution map is sufficiently regular [32,38].The second condition stems from the parametric complexity of the problem.It is satisfied when the parameterizations of the FOM are affine, i.e. they can be written as a linear combination of terms.When this does not occur, we can use a hyper-reduction technique to approximate the nonaffinities as a sum of a small number of affine quantities.
Hyper-reduction strategies are divided into two categories: collocation procedures [25,30,34], such as the Gauss-Newton approximated tensors (GNAT) method [8], which rely on sampling procedures applied to the nonaffine function; and function-reconstruction approaches, such as the empirical interpolation method (EIM) [4,18,20,26,28], its discrete variant (DEIM) [7,10,19,42], and the empirical interpolation method in matrix form (MDEIM) [5,9,36].The latter is a purely algebraic interpolation approach that was proposed to approximate parameterized bilinear operators.Although MDEIM has been successfully applied to unsteady problems [17,27], it was not originally conceived to reduce the temporal complexity of time-dependent bilinear operators.Additionally, the cost of its offline phase explodes when dealing with unsteady problems since an evaluation of the algebraic operator is required at every time step.
We list hereunder the contributions of our work: • We extend the implementation of MDEIM to compress bilinear operators in space and time.In particular, we approximate them as a linear combination of ST-RB basis vectors by coefficients satisfying an interpolation condition on a set of space-time indices chosen via a greedy algorithm.The same idea can be used to approximate the right hand side (RHS) of the problem.• We propose a functional alternative to the algebraic implementation of MDEIM to reduce the offline cost.Here, we first approximate the fields in the bilinear forms that depend nonlinearly on the temporal variable and the parameters (e.g.physical coefficients), similarly to what is done in [1,16,24,26].
Then, the bilinear operators obtained by substituting the nonaffine fields with their compressed versions are approximated with the standard algebraic MDEIM.The computational cost of this step is greatly reduced thanks to the preliminary compression of the fields.• We integrate the proposed hyper-reductions with a ST-RB state approximation, which results in a space-time MDEIM-RB (ST-MDEIM-RB) method.We derive a posteriori error bounds for each of the ST-MDEIM-RB approaches, which we show decay according to a user-defined tolerance that determines the accuracy of the methods, under suitable well-posedness and reducibility assumptions.• We investigate the numerical properties of the proposed ST-MDEIM-RB approaches to the heat equation with parameterized thermal diffusivity and boundary values and the unsteady Stokes equations with a paramaterized and time-dependent fluid viscosity and inflow rate.Both tests are conducted on 3D geometries.This article is organized as follows.In Sect.2, we provide the notation we employ throughout our work.In Sect.3, we discuss the nature and the implementation of the ST-RB approaches.In Sect.4, we review the current literature on MDEIM for steady problems, and we propose our novel MDEIM strategies for unsteady applications.In Sect.5, we present the ST-MDEIM-RB methods and prove a posteriori error estimates.The derivations in Sections 3-5 are shown using the heat equation as model problem.We consider the Stokes equations in Sec. 6.In Sect.7, we discuss the numerical results obtained when applying ST-MDEIM-RB to our test cases.Finally, in Sect.8 we draw some conclusions and possible extensions.

HYPERMATRICES AND BASIC NOTATION
In this work, we manipulate 3-dimensional arrays, which we call hypermatrices.Given a hypermatrix, we label its axes with specific subscripts (kept separate by commas).In particular, we use the subindices s, t and µ to indicate the spatial, temporal and parametric axes, respectively.For example, U s,t,µ ∈ R Ns×Nt×Nµ can be interpreted as an array of N µ HF (snapshot) matrices of size N s × N t , where N s , N t and N µ represent the number of HF spatial degrees of freedom (DOFs), HF time steps, and parameters.Our algorithms extensively employ the TPOD operation to compress hypermatrices in specific indices.We use the subscripts ŝ and t to indicate a reduction in the spatial and temporal dimensions, respectively, and we denote the number of reduced spatial and temporal DOFs as n s and n t .
We now define two basic operations on hypermatrices, i.e. permutation of axes and flattening.The former is a generalization of matrix transposition in multiple dimensions.For example, U t,s,µ ∈ R Nt×Ns×Nµ is the result of permuting space and time axes of U s,t,µ , and is marked by a different ordering of the subindices.The latter is equivalent to the operation of merging axes, e.g., U s,tµ ∈ R Ns×NtNµ is the result of flattening the time and parameter indices of U s,t,µ into a single index, ending up with a matrix.Since these two operations are isomorphisms, we do not need to distinguish between the hypermatrix before and after the permutation/flattening.We use the standard notation ∼ = to denote congruence by isometry of two objects.For example, U t,sµ ∼ = U s,t,µ , as they provide the same information while being indexed differently (due to flattening and permutation of indices).For the product of hypermatrices, we use the same notation as for matrices, e.g., U s,ŝ V ŝ ∈ R Ns .
The computation of a hypermatrix TPOD is a generalization of a matrix standard TPOD.We define the space-axis TPOD as: where I stands for an identity matrix (with axes determined by the subindices), || • || F indicates the Frobenius norm of a matrix, and n s is the rank of the truncation, defined as in (9).We are flattening the time and parameter axes and compressing the N t N µ space snapshots into n s space basis vectors.The result is a matrix Φ s,ŝ that collects the reduced space basis vectors in the columns.We interpret such a matrix as a basis spanning an n s -dimensional spatial subspace of R Ns on which we approximate the spatial evolution of the HF solution.
The matrix Φ ŝ,s is the transpose of Φ s,ŝ and is used to spatially compress the snapshots by projecting them onto range(Φ ŝ,s ).Using our notation, this corresponds to computing U ŝ,tµ = Φ ŝ,s U s,tµ , where • denotes a compressed quantity, i.e. it can be written as a linear combination of the RB vectors.Note that U s,tµ = Φ s,ŝ U ŝ,tµ is still in the reduced space, even though written in the HF basis.Different options are available to achieve a reduction in time.In this work, we choose the sequentially truncated high order singular value decomposition (ST-HOSVD) approach defined in [6].We first compute the space-axis TPOD, and then we compute the time-axis TPOD of the reduced space snapshots as: Here, we are compressing the n s N µ space-compressed and parameter snapshots into n t time basis vectors.The time-axis TPOD on the already compressed hypermatrix in space aims to reduce the computational cost of the operation.Similarly to the spatial basis, we interpret Φ t, t as a basis spanning an n t -dimensional temporal subspace of R Nt on which we approximate the temporal evolution of the HF solution.
We can now define the space-time RB as where ⊗ denotes the standard Kronecker product between matrices.We aim to approximate the evolution of the HF solution in space-time on the subspace spanned by the n st = n s n t columns of Φ st,ŝ t.We denote the space-time compression of U s,t as U ŝ, t ∼ = Φ ŝt ,st U st , whereas U s,t ∼ = Φ st,ŝ t U ŝt is the expansion of the compressed quantity on the HF DOFs.Lastly, A B (resp.A B, A B) denotes A ≤ cB (resp.A ≥ cB, A = cB ) for a constant c; and we use N(M ) to indicate {1, . . ., M }, with M ∈ N + .

SPACE-TIME REDUCED BASIS METHODS FOR A PARAMETERIZED HEAT EQUATION
We commence this section by introducing the FOM associated with a parameter-dependent heat equation.Consequently, we present in detail the offline and online phases of the ST-RB method applied to the FOM.
where Let us introduce the Hilbert spaces Let us consider a conforming quasi-uniform partition T h of Ω with meshsize h and a uniform partition of the time domain (0, T ) with time-step size δ = T /N t ; we define the sequence of time steps {t n } Nt n=0 with t n = nδ.For the spatial approximation of (4), we consider a conforming FE space V h ⊂ V on T h and V 0 h = V h ∩ V 0 Γ D .For the time discretization, we consider the Backward Euler (BE) scheme for simplicity in the exposition.To strongly impose non-homogeneous Dirichlet boundary conditions, we define a FE lifting L h (g µ ) of the Dirichlet value g µ .Let u µ h,0 be a FE interpolation of u µ 0 .The fully discrete problem at time value n ∈ N(N t ) reads: given where We can write this problem in algebraic form as: where U µ s,n ∈ R Ns is the vector of DOFs of u µ h,n .M s,s (which is independent of (t, µ) throughout this work), A µ s,s and L µ s are the mass matrix, the stiffness matrix and the RHS vector, respectively.We also introduce the symmetric, positive definite matrix X s,s = M s,s + A s,s , representing the V inner product on V h .(Note that A s,s is the discrete Laplacian with unit diffusivity.)We can also state the space-time algebraic system where K µ st,st is a bi-diagonal block-matrix with N t diagonal blocks δ −1 M s,s + A µ s,s (t i ) and lower diagonal blocks δ −1 M s,s .We can define a global spatio-temporal norm matrix X st,st , which is a block-diagonal matrix with N t blocks δX s,s ; the δ serves for this matrix to represent the L 2 (0, T ; V) product for the FE basis defined above.We also introduce the (X st,st , X −1 st,st ) space-time matrix norm: In the course of this work, we use the property where ||X −1 st,st || 2 acts as an equivalence constant between the (X st,st , X −1 st,st ) and 2 norms, and scales with δ −1 (h −d + h 2−d ).(In this work, the || • || 2 of a square matrix is to be understood as its spectral radius.)We lastly introduce the coercivity constant of K µ st,st which is positive by virtue of the coercivity of K µ st,st .When N st = N s N t is large, solving (7) for a variety of parameters becomes prohibitively expensive, even when exploiting parallel computations [23].In this scenario, devising a ROM such as the ST-RB method (which we address in the following subsections) that overcomes this issue becomes paramount.
3.2.Offline phase.We consider a set of parameters D off = {µ k } Nµ k=1 ⊂ D, and the FOM snapshots hypermatrix U s,t,µ , with U µ k s,t being the solution of (8) for µ k .By following the approaches in Sect.2, we can extract from U s,t,µ the spatial basis Φ s,ŝ , the temporal basis Φ t, t, and consequently the space-time basis Φ st,ŝ t.
We now provide the details regarding the choice of the TPOD ranks n s and n t .In the following discussions, we take the spatial basis as a point of reference to avoid redundancy.We start by recalling the following standard result of TPOD.We refer to [32] for additional details.
Proposition 1.Let {σ s,i } Nσ s i=1 be the N σs = min{N s , N t N µ } singular values (sorted by decreasing absolute value) of U s,tµ .Then: Prop. 1 provides a criterion for the selection of n s : given a user-defined tolerance ε ∈ R + representing a level of accuracy, we consider the minimum value of We essentially require the relative energy (i.e. the squared 2 norm) retained by the n s columns of Φ s,ŝ to be larger than 1 − ε 2 .We represent with TPOD s,ŝ (•) the map that, given the snapshot matrix U s,tµ , returns the RB Φ s,ŝ that holds (1) for n s in (9).We omit ε, since we will consider the same value in all appearances.Similarly, TPOD t, t (•) represents the time-axis TPOD that takes U t,ŝµ and computes Φ t, t that holds (2) up to ε, using a temporal version of the criterion (9).For the RB method to provide a good approximation, we must ensure that: • the problem we approximate is reducible (see Sect. 1).If not, the TPOD s,ŝ (•) will select n s ∼ N s , i.e. (close to) the entire column space of U s,tµ , thus preventing us from achieving a meaningful compression; • N µ is large enough and the choice of the sampling space is appropriate, otherwise the testing error will likely be (considerably) larger than the training error; • ε is low enough to attain a good enough accuracy, but not so low that TPOD s,ŝ (•) selects n s ∼ N s basis vectors.
Remark 1.Given a symmetric, positive definite matrix Y s,s , we can run a modified TPOD s,ŝ (•) so that the resulting spatial basis is orthogonal in the Y s,s -norm instead of the 2 norm.By definition, Y s,s admits a Cholesky decomposition Y s,s = H T s,s H s,s , with H s,s upper-triangular.We compute s Φ s,ŝ = TPOD s,ŝ (H s,s U s,tµ ), and then Φ s,ŝ is found as Φ s,ŝ = H −1 s,s s Φ s,ŝ .We use Y s,s = X s,s to compute the spatial basis in the numerical experiments in Sect.7. In time, 2 -orthogonality of Φ t, t is sufficient (i.e.Y t,t = I t,t ).
At this point, we have discussed the implementation details for constructing the spatial and temporal bases and, by combining them, the space-time basis.Now, we show that the ST-RB approximation is well-defined.Proposition 2. Let U st,µ = Φ st,ŝ tΦ ŝt ,st U st,µ be the space-time TPOD approximation in (3), and let U s,tµ = Φ s,ŝ Φ ŝ,s U s,tµ be the hypermatrix of space-only approximations.Moreover, let {σ t,i } Nσ t i=1 be the N σt = min{N t , n s N µ } singular values (sorted by decreasing absolute value) of U t,ŝµ .The following error bound holds: Proof.We split the total error as , where the first term accounts for the space reduction and the second for the time reduction.By virtue of Prop. 1 we have: Now we focus on the term U µ s,t − U µ s,t ; exploiting Prop. 1, the definition of TPOD t, t (•), the definition of Frobenius norm, and the property of invariance to cyclic permutations of the trace of a matrix we have: The statement follows from ( 11)-( 12).
We mention the following corollary of Prop. 2 and ( 9): Although this error bound is not as sharp as the one in (10), it has the advantage of relating the error to a user-defined quantity.
3.3.Online phase.In this phase, two main tasks are carried out.First, the projection of ( 8) onto the space-time reduced subspace.Second, solving the resulting (linear) reduced system of equations.In our work, we consider a Galerkin projection.We refer to [13] for a detailed review of the more general Petrov-Galerkin projections in an RB context.Adopting an algebraic notation, the ST-RB problem reads as We can also express (14) as K µ ŝt ,ŝ t U µ ŝt = L µ ŝt , where st , denote the Galerkin compression of the space-time left hand side (LHS) and RHS, respectively.These quantities are computed following the steps outlined in [37].

EMPIRICAL INTERPOLATION METHODS FOR STEADY AND UNSTEADY PARAMETERIZED PDES
We commence this Sect.by recalling the properties of MDEIM applied to a steady problem in Sect.4.1.Then, in Sects.4.2-4.4,we propose new MDEIM architectures for unsteady applications.4.1.Review of MDEIM for steady PDEs.Let us consider a parameterized, steady, nonsingular matrix A µ s,s ∈ R Ns×Ns .We introduce the vector of nonzero entries A µ s ∈ R Nz associated to A µ s,s (the two are congruent by isometry).MDEIM is a data-driven procedure hinging on: • running a TPOD on the snapshots hypermatrix to form a basis Φ a s,ŝ ∈ R Nz×n a s spanning a subspace for range(A s,µ ); • the iterative selection of a list of indices {I 1 , . . ., I n a s }, with I i ∈ N(N z ) ∀i, which are chosen to minimize an interpolation error over nested subspaces of range(A s,µ ), in the ∞ norm.We associate with this list the sampling matrix P a s,ŝ = e s (I 1 ), . . ., e s (I n a s ) ∈ {0, 1} Nz×n a s , where e s (I) ∈ R Nz is the I th vector of the N z -dimensional canonical basis.The implementation of MDEIM is summarized in Algorithm (1).
Algorithm 1 MDEIM algorithm Update P a s,ŝ = P a s,ŝ , e s (I k ) , where I k = arg max |r| 9: end for At the k th iteration, we compute the error between the k th basis vector and its projection on the first k − 1 columns of the MDEIM basis, evaluated at the k − 1 sampling points encoded in P a s,ŝ .Then, P a s,ŝ is updated by adding the index corresponding to the largest entry (in absolute value) of the residual.The online phase of MDEIM consists in computing the approximation where the vector of coefficients A ŝ is computed by enforcing the interpolation condition: An approximation A µ s,s ≈ A µ s,s on the full set of HF DOFs can be recovered by following the steps in [27].We remark that the procedure is well posed, since P a ŝ,s Φ a s,ŝ is nonsingular [10].
Remark 3. The MDEIM approximation A µ s,s of a symmetric, positive definite matrix A µ s,s preserves the symmetry, but not necessarily the positive definiteness.Alternative approaches that bypass this issue include finding the vector of coefficients A ŝ(µ) subject to a generalized linear constraint that forces A µ s,s to be coercive [9], or directly performing interpolation on the manifold of symmetric, positive definite matrices [8].Finally, we underline that A µ s,s is always non-singular, even in the case of non-symmetric input matrices [27].This observation alone supports the derivation of a posteriori bounds for the error that the approximation in (15) commits.
Remark 4. As shown in (16), computing A ŝ(µ) requires assembling the HF matrix A µ s,s at the indices encoded in P a s,ŝ during the online phase.This can be done efficiently by constructing a reduced integration domain, as stated in [27].Since, in a FE setting, assembling A µ s,s requires performing an element-wise numerical integration, we can devise an efficient FE assembling procedure that performs the numerical integration only over the elements whose contribution to the entries of P a ŝ,s A µ s is nonzero.

4.2.
Standard MDEIM for unsteady PDEs.In the framework of unsteady parameterized PDEs, the authors of [27] present an algorithm that progressively builds a global spatial reduced basis for the state variable and one for each of the nonaffinely parameterized FE matrices and vectors.However, even though their simultaneous system approximation and state-space reduction approach appears to be reasonable regarding memory footprint, its wall time is considerable.This claim becomes even stronger when the reduction occurs in space and time, as in our framework.For this reason, we propose an alternative procedure to separate the system approximation and state reduction (discussed in Sect.3) steps.For the former, we can devise a standard MDEIM (STD-MDEIM) that conceptually reuses the same steps of a MDEIM for steady problems.This method is outlined in Algorithm 2.
Algorithm 2 STD-MDEIM algorithm.The STD-MDEIM approximation of the nonaffine stiffness matrix reads, for any µ: This approach aims to find a spatial approximation of the nonaffine matrix at every time step, and as a result, it does not entail a reduction of the time complexity.Thanks to Prop. 1, we can bound the error of the STD-MDEIM approximation in (17) as: where χ a s = ||(P a ŝ,s Φ a s,ŝ ) −1 || F .This quantity represents the error committed by the MDEIM interpolation.Since STD-MDEIM requires computing an n a s -dimensional vector at N t time steps, it achieves a dimension reduction proportional to N z /n a s .In applications where N t is large, recovering a reduction in both space and time is desirable.The method outlined in Sect.4.3 is designed to achieve this property.
Remark 5. Unfortunately, the theoretical lower bound for χ a s is exponential with respect to the reduced dimension n a s .As a remedy, we can utilize alternative selection strategies of the MDEIM indices that enjoy error bounds that grow polynomially.We refer to [15,29] for more details.In any case, the exponential bound is observed to be quite pessimistic in practice; the MDEIM compression achieves far better results in most scenarios.We employ the standard selection process of the MDEIM indices for the numerical experiments in Sect.7.

4.3.
Space-Time MDEIM for unsteady PDEs.The STD-MDEIM procedure aims to find a suitable subspace for the efficient approximation of A µ s (t) at a given time instant t.We propose the space-time MDEIM (ST-MDEIM) approach to construct a reduced subspace for the space of the stiffness matrices across every time instant and parameter.The procedure is outlined in Algorithm (3).The ST-MDEIM approximation of the nonaffine stiffness matrix is the generalization of (17) at every time step, and reads for any µ as: ŝ,s ⊗ P a t,t identifies a space-time reduced integration domain.Thus, by extending the procedure outlined in Remark 4 to space-time, the ST-MDEIM coefficients can be efficiently computed with an online cost scaling with n a st , i.e. with a dimension reduction proportional to N z N t /n a st .By virtue of ( 13) we can bound the error of the space-time approximation in (19) as where || F , and A ŝ,tµ = Φ a ŝ,s A s,tµ .For ST-MDEIM to achieve the same accuracy as STD-MDEIM, we must ensure that the manifold of stiffness matrices is reducible in time, as well as in space.Additionally, depending on the nature of the problem at hand, it might be necessary to select a higher value of N µ to achieve a space-time compression at the specified tolerance.We further elaborate on this claim while discussing our numerical results in Sect.7.

4.4.
Functional MDEIM for unsteady PDEs.Even by employing the hyper-reducing techniques outlined in Remark 4, constructing the snapshots hypermatrix A s,t,µ used in Algorithms 2-3 requires a considerable number of operations.To limit the computational burden of this task, we can devise an approach leveraging a functional representation of A µ s,s (t).We recall that this quantity is obtained by integrating the bilinear form defined in (6) over the FE quadrature points for all combinations of trial and test functions.Instead of computing a TPOD directly on several snapshots of the matrix, the idea of this method is first to attain a compression of the nonlinearly (t, µ)-dependent field α µ (t) via TPOD, similarly to what is done in [1,16,24,26]; and then, to run MDEIM on the matrix snapshots that arise by plugging in the newly obtained (t, µ)-independent fields in the original form.This step serves to improve the compression achieved by the aforementioned methods, and when applicable limits the additional complexity introduced by geometrical parameterizations, thus overcoming the inefficiency that [27] points out as a deterrent from functional-based (e.g.EIM or DEIM) hyper-reducing techniques.Additionally, any problem specificity and intrusive changes to the FE solver which (D)EIM might suffer -as mentioned in [27] -can be bypassed by an appropriate implementation of the numerical algorithm (see Remark 6).Lastly, the functional approaches we propose can be easily extended in space-time, thus recovering the same dimensional reduction achieved by ST-MDEIM.
Let us consider the hypermatrix α s,t,µ ∈ R Nq×Nt×Nµ of evaluations of α µ (t) on the N q FE quadrature points, from which we extract the spatial basis Φ α s,ŝ = TPOD s,ŝ (α s,t,µ ) ∈ R Nq×n α s .We also introduce a discontinuous Galerkin space Q h whose nodes coincide with the N q quadrature points.Moreover, we introduce the functions φ α s,1 , . . ., φ α s,n α s ∈ Q h , obtained by interpolating the n α s basis vectors of Φ α s,ŝ on Q h .(We underline that φ α s,i ∼ = φ α s,i ∀i, where φ α s,i is the i th column of Φ α s,ŝ , and thus we make no distinction between the two.)Let us now consider the bilinear form obtained by substituting α µ k (t) with φ α s,i in the bilinear form defined in (6).We now introduce the FE matrix s A s,s (φ α s,i ) obtained by performing numerical integration on the form in (21), the corresponding vector of nonzero entries s A s (φ α s,i ) and the matrix that collects all these vectors The scope of this functional version of MDEIM (FUN-MDEIM) is the computation of s A s,ŝ itself, from which we can then retrieve a matrix approximation by following the steps outlined in Sect.4.2.Specifically, we compute s Φ a s and s P a s by applying STD-MDEIM on s A s,ŝ , and we approximate A µ s (t) for any (t, µ) as: The approximation above is derived without assembling a huge snapshots hypermatrix and running an expensive TPOD, as is done in STD-MDEIM.Instead, FUN-MDEIM proposes to compute a small hypermatrix α s,t,µ (N q N z ), compress it with a relatively inexpensive TPOD, and then run STD-MDEIM on the much smaller s A s,ŝ (n α s N t N µ ).Before formally deriving the error bounds for this method, we enunciate the following result.

Lemma 1. Let us consider s
A s,ŝ defined in (22).The following result holds: Proof.We recall that for any j ∈ N(N t ) and k ∈ N(N µ ), we can express A µ k s,s (t j ) as: , where T a : R Nq → R Ns×Ns is a continuous operator representing the FE assembling procedure.(T a depends on hyperparameters related to the FE assembly, such as the order of quadrature and mesh, not explicitly mentioned in the definition, as the steps in this proof are independent of them.)We can express V s and s V s as Exploiting the continuity of T a , we have: Now, we choose s V jk α µ k s (t j ), and thus by Prop. 1 and definition of TPOD Algorithm 4 STFUN-MDEIM algorithm.
where s Note that this bound is not necessarily less sharp than the one in (18), and indeed our results in Sect.7 suggest that STD-MDEIM and FUN-MDEIM exhibit the same order in ε.
We can extend the aforementioned concepts to the space-time case.We introduce the time basis Φ α t, t ∈ R Nt×n α t extracted via time-axis TPOD on α t,sµ , and its associated sampling matrix P α t, t.This methodspace-time functional MDEIM (STFUN-MDEIM) -aims to approximate the stiffness matrix for any µ as: The implementation of STFUN-MDEIM is reported in Algorithm 4 (FUN-MDEIM can be executed by computing only the spatial quantities in this algorithm).To derive the error bounds for STFUN-MDEIM, we enunciate the following space-time generalization of Lemma 1.

Lemma 2. Let s
A s,ŝ be defined as in Lemma 1, let α ŝ,tµ = Φ α ŝ,s α s,tµ and let Φ α t, t be the time basis extracted from α t,ŝµ .Then Proof.We define n α st = n α s n α t , and we associate with any given space-time index i ∈ N(n α st ) a pair of spatial, temporal indices (i s , i t ) ∈ N(n α s ) × N(n α t ).Let us initially consider the case V st = A µ 1 st ; proceding as in Lemma 1, we have that: The statement follows if we consider a generic , we follow a space-time version of (24) and we exploit (13).
Lemma 2 allows us to derive the following error bound for the approximation (26): where s || F .As we will show in Sect.7, FUN-MDEIM and STFUN-MDEIM achieve a cheap offline phase compared to their fully algebraic counterparts, thanks to the compression attained by compressing the data encoded in α s,t,µ .

EMBEDDING MDEIM IN ST-RB: THE ST-MDEIM-RB APPROACHES
In this section, we discuss the properties of the methods that we obtain by combining the MDEIM approximation of the nonaffine HF quantities with the ST-RB procedure.We can identify three main steps in their implementation: • we find the space and time bases Φ s,ŝ and Φ t, t spanning the reduced solution subspace, following the steps outlined in Sect.3; • we perform MDEIM of the nonaffine operators following one of the four approaches presented in Sect.4, retrieving an approximate affine decomposition; • thanks to the previous step, we can efficiently compute the space-time Galerkin projection of the problem at hand.
In the following subsections, we firstly present the problem formulation of a ST-MDEIM-RB method, and secondly we discuss a posteriori theoretical error bounds for this class of ROMs.We make use of the symbol •, which refers to a double compression: we approximate a quantity with MDEIM (first compression), then we reduce it with ST-RB (second compression).

Implementation of ST-MDEIM-RB.
The MDEIM strategies outlined in Sect. 4 can be leveraged to recover the following approximated affine decompositions: We denote with K µ st,st the quantity obtained by substituting the stiffness matrix K µ st,st with its MDEIM approximation A µ s,s .The statement of a ST-MDEIM-RB problem applied to our heat equation reads as: The error that a ST-MDEIM-RB method commits with respect to the FOM can be split into two contributions: the first one is due to the MDEIM approximation of the nonaffine FE quantities (system approximation) and the second one arises from the ST-RB compression (state approximation).In Sect.5.2, we perform an a posteriori error analysis for ST-MDEIM-RB methods.

Error analysis of ST-MDEIM-RB.
In this subsection, we derive a posteriori error estimates of the proposed ST-MDEIM-RB approaches under the spatio-temporal X st,st norm.For the sake of conciseness, we assume that both the Dirichlet and Neumann boundary conditions are homogeneous; however, the following results can be easily generalized to non-homogeneous cases.
Theorem 1.Let us consider the well-posed problem given by (8), and its well-posed ST-MDEIM-RB approximation in (28).The following approximation holds: , where U µ st is the approximate solution obtained with STD-MDEIM.
Proof.From the FOM (8), we have the identities which we conveniently split into the MDEIM error contribution E M st and the ST-RB error contribution E RB st .Left multiplying (30) by K µ st,st −1 , taking the X st,st norm and using the definition of β yields: corresponds to the last term in the statement (29).On the other hand, the MDEIM contribution can be bounded by recalling (18): From ( 31), the quantity can be bounded by the sum of the first two terms in the statement (29), thus concluding the proof.
When our problem features a reducible system, the error contribution associated with MDEIM (E M st in the proof) explicitly vanishes as the quality of the approximation improves.On the other hand, the quantity E RB st is a residual-based a posteriori ST-RB error estimate which also vanishes with ε when the state variable is reducible.This can be shown by using the following heuristics: when both the FOM and ST-RB formulations are well posed, the ST-RB solution will approximate the projection of U st on range(Φ st,ŝ t), which we have shown in (13) commits an approximation error of order ε.
As a corollary of Th. 1, we address the case of the remaining ST-MDEIM-RB methods.
Corollary 1.Let us consider the well-posed problem given by (8), and its well-posed approximation in (28).
The following approximation holds for ST-MDEIM: .
The following approximation holds for FUN-MDEIM: .
The following approximation holds for STFUN-MDEIM: Proof.The proof follows the same steps as the one for Thm. 1.The only difference comes when we bound the term , where instead of using (18) as we do in (31) we now employ: (20) for ST-MDEIM; (25) for FUN-MDEIM; and (27) for STFUN-MDEIM.

STOKES EQUATIONS
In this section, we briefly discuss the application of the ST-MDEIM-RB method to the unsteady, parameterized Stokes equations.We consider an inf-sup stable pair of FE spaces for velocity and pressure.By employing the BE time integrator, we obtain the following system of equations: We refer to [21,35] for a detailed exposition of the discretization of the Stokes problem.We represent with U µ s,n and P µ s,n the vectors of DOF values for velocity and pressure at t n , respectively.M s,s is the velocity mass matrix, A µ s,s (t n ) accounts for the viscous term, B s,s is the discrete divergence matrix and its transpose the discrete gradient; M s,s and B s,s are (t, µ)-independent in this work.
The ST-RB method for this problem features a Galerkin projection onto a cartesian product subspace, i.e. span(Φ u st,ŝ t × Φ p st,ŝ t), where Φ u st,ŝ t and Φ p st,ŝ t are the N u st × n u st -and N p st × n p st -dimensional ST-RB bases for velocity and pressure, respectively.To guarantee the inf-sup stability of the ST-RB framework, it is necessary to augment the reduced basis for the velocity with additional basis vectors -called supremizers -both in space and in time.For more details regarding supremizers, we refer to [3,14,31,33,37].In this work, we construct the supremizers following the ideas in [37].The ST-MDEIM-RB problem for the Stokes equation reads as: find We underline that the blocks B T ŝt ,ŝ t ∈ R n u st ×n p st and B ŝt ,ŝ t ∈ R n p st ×n u st are obtained without running a MDEIM procedure, given the independence of B st,st from (t, µ).Additionally, we recall that B T ŝt ,ŝ t = B T ŝt ,ŝ t (see [37] for more details).

NUMERICAL RESULTS
In this section, we investigate the numerical properties of the proposed ST-MDEIM-RB methods in the case of the heat and Stokes equations.In particular, we assess the following: • For the online phase, we measure (1) the errors between the HF solutions and the ones obtained with ST-MDEIM-RB, and (2) the computational speedup our algorithms achieve compared to the HF simulations.
As we wish to empirically verify the correctness of the error estimates derived in Sect.5.2, we make the measurements for several values of the TPOD tolerance.In particular, we consider ε ∈ {10 −i } 4 i=2 .• For the offline phase, we compare the wall time and memory footprint of the MDEIM approximation of the tangent matrix, for each of the proposed ST-MDEIM-RB strategies.Since the value of ε has a minimal impact on these quantities, we present the results obtained with ε = 10 −4 .We run the simulations on a fixed geometry, which is shown in Fig. 1.For both tests, we select a training set according to a uniform distribution on D, i.e.D off = {µ k } Nµ k=1 with N µ = 80, and we compute the FOM solutions for every µ ∈ D off .We use D off to construct the ST-RB approximation, and only the first 30 parameters of D off to obtain the MDEIM approximation, as empirical evidence suggests that a smaller number of snapshots suffices for this task.We evaluate the online performance of our algorithms on a testing set D on = {µ k } Non k=1 ⊂ D, with N on = 10, for which we also need to compute the corresponding HF solutions.To correctly estimate the performances, it is imperative to ensure that D off ∩ D on = ∅.The accuracy of the methods is evaluated according to the following measures: .
The norm matrix X u st,st is equal to X st,st and X p st,st ∈ R N p st ×N p st is the matrix representing the discrete L 2 (Ω)-2 norm in space-time (which we employ for the pressure in the Stokes test case).Concerning the computational efficiency of the methods, we compute the speedup defined as the ratio between the average wall time of a FOM  simulation and that of a ST-MDEIM-RB simulation.All the numerical tests are run on a local computer with 66Gb of RAM and an Intel Core i7 running at 3.40GHz.
Remark 6.For the generation of the snapshots, we employ Gridap [2,41], a library supporting a FE approximation of PDEs in the Julia programming language.Among the numerous useful features of Gridap, we mention that this library resorts to a cell-wise representation of fields while integrating the terms of the weak formulation associated with the problem at hand.This allows us to implement with ease our functional approaches, bypassing the implementation issues that [27] points out as deterrents from functional methods.Indeed, we can simply compress the nonaffinely parameterized fields in a cell-wise manner (i.e. on the quadrature points), and then define an approximated weak formulation where the parameterized fields are replaced with their compressed versions.The numerical solvers that Gridap makes use of can still be used in this scenario.
7.1.Heat equation.In this part, we solve the heat equation introduced in Sect.3, with the following parametric data: and f µ (x, t) = 1, u µ 0 (x) = 0.The Dirichlet boundary condition (BC) is imposed strongly, meaning that a lifting operator for the Dirichlet datum is explicitly computed and moved to the RHS, as shown in (5).The unknown temperature is measured in Celsius.
Tb. 2 collects the information related to the MDEIM offline phases.In particular, the functional methods improve by 4 -5 times the performance of their fully algebraic counterparts, both in terms of wall time, and memory allocation.This result reflects the cheaper construction of the snapshot matrices, as well as the reduced cost of the TPODs we run to extract the bases, as discussed in Sect.4.4.We also infer that the cost of computing the temporal basis of the tangent matrix is marginal, as we can tell from the negligible difference between the MDEIM methods in space-time, compared to their space-only MDEIM counterparts.Concerning the online phase, Tb. 3 presents the efficiency and accuracy measures the ST-MDEIM-RB methods achieve with respect to the HF simulations, for three different TPOD tolerances.We first observe that the error estimates decrease at the same rate as ε, which empirically validates the convergence estimates derived in Sect.5.2.We deduce that the use of a functional MDEIM approach does not negatively impact the online accuracy, thus providing empirical evidence of Lemmas 1-2.We thus infer that the functional MDEIM retains the same properties as the fully algebraic ones online, while greatly reducing the offline cost, confirming our prior expectations.On the other hand, employing a space-time reducing technique significantly increases the achieved speedup: in particular, ST-STRB and STFUN-STRB are roughly 12 to 20 times faster than their space-only counterparts, depending on the value of ε.This is entirely due to the compression of the temporal information these methods achieve.Even greater speedups are expected when considering a larger number of temporal DOFs than in the current setting (N t = 60).Such compression is also responsible for a marginal loss of accuracy compared to the spatial methods, as we can observe from the error estimates.An important takeaway from Tb. 2 is that it is necessary to select carefully the number of snapshots we employ to construct the MDEIM bases when employing ST-STRB or STFUN-STRB, if we desire to achieve an accuracy that decays with ε.As we can notice from the error measures in brackets (corresponding to the case in which only 20 snapshots are used, instead of 30), the temporal bases fail to correctly approximate the evolution in time of the system when ε = 10 −4 .We have investigated the source of the error, and we have discovered that there exists some µ * ∈ D on such that || L µ * st − L µ * st || 2 10χ l s ε||L s,tµ || F , i.e. a tenfold increase of a priori expectations, thus worsening the average relative error.We conclude that the significant online speedup that characterizes ST-STRB or STFUN-STRB comes at the (potential) cost of performing a more expensive offline phase, because a (potentially) higher number of snapshots is required to ensure an accurate temporal compression.In Fig. 2 we compare the temperature obtained by solving the FOM in (7) and f µ (x, t) = h µ (x, t) = 0, u µ 0 (x) = 0. Here, n Γ D is the outer normal vector to the inlet Γ D , where we impose a parabolic in space, and periodic in time inflow profile.The Dirichlet BC is once again imposed strongly, and we measure the unknowns in the cgs (centimeter-gram-second) unit system, i.e. we express the velocity in cm/s and the pressure in dyn/cm 2 , where dyn = g • cm/s 2 .As we show in Tb. 4, due to the large number of entries of the spatio-temporal tangent matrix, the offline phases of the fully algebraic MDEIM approaches are too expensive for our computational resources.Therefore, in this Sect.we are forced only to present the numerical results obtained with the functional methods; these approaches take approximately 1min to compute the MDEIM approximant of the LHS matrix (albeit with a considerable memory footprint).In Tb. 5 we collect the results relative to their online phase.The conclusions we draw are coherent with the ones of the heat equation: both methods achieve an accuracy that decays with ε, and STFUN-MDEIM is computationally far more efficient in terms of wall time (depending on the value of ε, it is 23 to 33 times faster than FUN-MDEIM).The difference in wall time between the two methods is more pronounced than in the heat equation because of the greater number of DOFs characterizing our Stokes equations.In general, we expect this difference in performance to increase with the the number of DOFs, given that the difference in dimension reduction (see Sect. 4.2-4.3)will tend to increase with the number of DOFs.The error estimates are smaller than the ones obtained in the heat equation because the state and system approximation errors are both more contained.The former is due to the shorter time interval, which results in a more accurate description of the solution manifold.The latter comes from the simpler parameterization of the problem, as the Neumann BC is homogeneous in this test case.In this test case, extracting the MDEIM bases from 20 snapshots instead of 30 has a marginal effect on the error estimates relative to STFUN-MDEIM (the performance worsens when considering 10-15 snapshots).Lastly, in Fig. 3 we compare the velocity and pressure obtained by solving the FOM in (32) and the ones computed by employing STFUN-MDEIM, at t = 0.075s, with µ = [5.58,1.20, 6.12], and employing ε = 10 −4 .for the approximation of the tangent matrix of the Stokes equations.The results relative to STD-MDEIM and FUN-MDEIM are not displayed since the memory that they require to perform the task exceeds the maximum memory limit that our computational resources allow.

CONCLUSIONS
In this work, we discuss the implementation of four ST-MDEIM-RB schemes, consisting of a MDEIMbased system approximation technique embedded in a ST-RB procedure.The first method reduces the spatial dimensionality of nonaffinely parameterized operators at every time step.The second approach achieves an efficient compression of the operator in both space and time.These two strategies proceed in a purely algebraic manner.On the other hand, the remaining two "functional" methods firstly compress the nonaffine fields evaluated on the quadrature points, thus obtaining a reduced version of the operators; then, they compress the reduced operators in space or space-time following the same steps as their algebraic counterparts.The numerical analysis of the ST-MDEIM-RB algorithms shows that the error decays with a user-defined tolerance that determines the dimension of the RB approximation.We assess the performance of the proposed methods on two test cases, namely a heat equation and a viscous flow governed by the Stokes equations on fixed 3D domains.In each of the tests we consider, the ST-MDEIM-RB algorithms converge as predicted by the theoretical analysis, and substantially shorten the wall time of the HF simulations.The methods that rely on a space-time compression are the most efficient ones.On the other hand, the functional approaches are characterized by a cheap offline phase, thanks to the compression of nonaffine fields.Therefore, we can identify a particular approach, i.e.STFUN-STRB, that achieves both cheaper offline and online phases.
An immediate extension of our work consists of considering nonlinear test cases, such as the Navier-Stokes equations.The proposed methods achieve good results on problems that are reducible when using linear ROMs.We plan to consider space-time nonlinear ROMs in the future, e.g. using neural network-based models with an encoder-decoder structure [12].

3. 1 .
Full order model.We consider an open, polyhedral Lipschitz domain Ω ⊂ R d with boundary ∂Ω, and a temporal domain [0, T ] ⊂ R + ∪ {0}.We also consider a parameter space D ⊂ R p .Given µ ∈ D, a generic heat equation defined on the space-time domain Ω × [0, T ] reads: are the Dirichlet and Neumann boundary data, respectively, u µ 0 : Ω → R is the initial condition, n is the outward unit normal vector to ∂Ω, and Γ D and Γ N are respectively the Dirichlet and Neumann boundaries (such that {Γ D , Γ N } forms a partition of ∂Ω).The superindex (•) µ indicates the dependence of a variable (•) upon µ.

FIGURE 1 .
FIGURE 1. Geometries used for the heat equation (left) and for the Stokes equations (right).We assign appropriate boundary tags in both cases: Γ D , Γ 0 D , and Γ N P D denote a non-homogeneous, a homogeneous, and a no-penetration Dirichlet boundary; Γ N , and Γ 0 N denote a non-homogeneous and homogeneous Neumann boundary.

FIGURE 3 .
FIGURE 3. From left to right: surface line integral convolution (SLIC) of the HF velocity; SLIC of the STFUN-MDEIM velocity; magnitude of the pointwise velocity error; HF pressure; STFUN-MDEIM pressure; and pointwise pressure error.The results are shown from a top view of the geometry longitudinal mid-section, at t = 0.075s, with µ = [5.58,1.20, 6.12] and employing ε = 10 −4 .

TABLE 1 .
From left to right: information related to the spatial domain, temporal domain, parameter space, and the FOM.For the latter, we indicate the FE polynomial basis, the spatiotemporal FOM DOFs, and the average wall time of a HF simulation.The quantities (L, W, H) denote the length, width and height of the geometry, whereas R is the radius of the three cylinders.
with the one computed by employing STFUN-STRB, at t = 0.15s, with µ = [7.37,9.58, 4.05], and employing ε = 10 −4 .(The solutions obtained with the other ST-MDEIM-RB procedures are omitted to avoid redundancy, given that they are graphically indistinguishable from the ones obtained with STFUN-STRB.)

TABLE 2 .
Wall time (in s) and memory footprint (in Gb) of running MDEIM with ε = 10 −4 for the approximation of the tangent matrix of the heat equation.

TABLE 3 .
Online computational speedup and accuracy achieved by the ST-MDEIM-RB approach compared to the HF simulations in the heat equation test case, for three different values of ε.The column relative to the accuracy of ST-STRB and STFUN-STRB displays in brackets the errors obtained when the number of snapshots used to compute the MDEIM bases is 20, instead of 30.The values are not displayed for STD-STRB and FUN-STRB, since changing the number of snapshots from 30 to 20 does not affect their error.The results are averaged over 10 different values of µ.7.2.Stokes equation.In this subsection, we solve the Stokes equations with the following parametric data: