{'fraction_non_alphanumeric': 0.04461477204857655, 'fraction_numerical': 0.015508660163249054, 'mean_word_length': 4.419246952206279, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Projection-based model reduction of dynamical systems using space-time subspace and machine learning

This paper considers the creation of parametric surrogate models for applications in science and engineering where the goal is to predict high-dimensional spatiotemporal output quantities of interest, such as pressure, temperature and displacement fields. The proposed methodology develops a low-dimensional parametrization of these quantities of interest using space-time bases combining with machine learning methods. In particular, the space-time solutions are sought in a low-dimensional space-time linear trial subspace that can be obtained by computing tensor decompositions of usual state-snapshots data. The mapping between the input parameters and the basis expansion coefficients (or generalized coordinates) is approximated using four different machine learning techniques: multivariate polynomial regression, k-nearest-neighbors, random forest and neural network. The relative costs and effectiveness of the four machine learning techniques are explored through three engineering problems: steady heat conduction, unsteady heat conduction and unsteady advective-diffusive-reactive system. Numerical results demonstrate that the proposed method performs well in terms of both accuracy and computational cost, and highlight the important point that the amount of model training data available in an engineering setting is often much less than it is in other machine learning applications, making it essential to incorporate knowledge from physical models. In addition, simpler machine learning techniques are seen to perform better than more elaborate ones.


Introduction
Many tasks in computational science and engineering are many query in nature, as they require the repeated simulations of a parameterized large-scale computational model. Model reduction has become a popular approach to make such tasks tractable. Most of such techniques first perform an "offline" training stage that simulates the computational model for multiple input-parameter instances to build a low-dimensional subspace or manifold. Then during an "online" deployed stage, these techniques reduce the dimensionality and complexity of the original computational model at arbitrary unseen input-parameter instances by performing a projection process of the original computational model onto that low-dimensional subspace or manifold. The projection framework can be combined with various ways of representing parametric dependence, where the input parameters might describe, for example, material properties, system geometry, system configuration, initial conditions and/or boundary conditions. The output generally consists of quantities of interest in the form of whole fields or some extractions/functions of the field such as pressure, temperature or displacement fields.
For initial boundary value problems, many model reduction approaches based on time and frequency domains have been proposed. In particular, moment-matching methods are commonly used for frequency domain-based model reduction, such as the Lanczos algorithm [1], the Arnoldi algorithm [2], and the Krylov bases) [24,37,38,39,40,41] and very few works use only parameter-dependent generalized coordinates (i.e., space-time bases) [42,43]. In comparison to these and previous works, the list below summarizes our contributions to the current state-of-the-art detailed in this paper: • Our proposed method can be considered as non-intrusive model reduction (or data-driven model reduction) and PDE-agnostic which means that one does not need to know (and solve) any underlying reduced partial differential equation (PDE) online. The only requirement is that full-order model snapshots are already available (to create space-time bases and train surrogate models).
• We compare various machine learning techniques to approximate the mapping between input parameters and output space-time bases expansion coefficients: multivariate polynomial regression, k-nearest neighbors, random forest and neural network.
• In the offline stage, we collect usual (space-time) state snapshots to build space-time bases and train the machine learning surrogate models. In the online stage, for an arbitrary input parameter we express the space-time solution as linear combinations of those space-time bases and compute the coefficients using the trained surrogate models.
• The space-time basis expansion coefficients only depend on input parameters and not on time.
• The proposed method is verified on three benchmark problems: parameterized steady heat conduction PDE, parameterized unsteady heat conduction PDE and parameterized advective-diffusive-reactive problem. Numerical results show that the online computation of machine learning surrogate models is extremely fast and the relative solution error is less than 1% on average.
• We show that while our proposed approach is less accurate than the least-square Petrov-Galerkin (LSPG) model reduction (without hyper-reduction) [44,45,46], it is significantly faster.
The paper is structured as follows. Section 2 characterizes the full-order models for both steady and unsteady problems. Section 3 details the proposed projection-based model reduction technique that address both steady (Section 3.1) and unsteady problems (Section 3.2). Section 4 describes four machine learning methods used to construct the mapping between input parameters and the output basis expansion coefficients. Section 5 reports the numerical experiments, and Section 6 concludes the paper.

Full-order model
We consider the (high-fidelity) full-order model (FOM) of a dynamical system represented as a parameterized system of ordinary differential equations (ODEs) r(x; µ) := f (x, t; µ) −ẋ = 0, x(0; µ) = x 0 (µ), (2.1) where t ∈ [0, T ] denotes time with final time T ∈ R + , and x : [0, T ] × D → R Ns denotes the time-dependent, parameterized state as the solution of problem (2.1) for the parameters µ ∈ D, and r : R Ns × D → R Ns is time-continuous residual (trivially zero) at time t. Here, D ⊆ R nµ is the parameter space, x 0 : D → R Ns is the parameterized initial condition, and f : R Ns × [0, T ] × D → R Ns is the velocity. Differentiation of a variable x with respect to time is denoted byẋ. Dynamical systems such as these may arise from the semidiscretization of a time-dependent PDE model. We refer to Eq. (2.1) as the FOM ODE. We note that whenẋ vanishes, the FOM model reduces to a steady parameterized system of nonlinear algebraic equations. In such a case, all time-derivatives in the following treatment disappear.
Numerically solving Eq. (2.1) requires a time-discretization method. In this work we use a linear multistep method. A linear k-step method applied to Eq. (2.1) leads to the algebraic system r n (x n ; µ) = 0, n = 1, . . . , N t , (2.2) where the time-discrete residual r n : R Ns × D → R Ns is defined as r n : (ξ; ν) → α 0 ξ − ∆tβ 0 f (ξ, t n ; ν) + k j=1 α j x n−j − ∆t k j=1 β j f (x n−j , t n−j ; ν). (2.3) Here, x k is the numerical approximation to x(k∆t; µ), ∆t ∈ R + is the time step, and coefficients α j and β j , j = 0, . . . , k with k j=0 α j = 0 define a particular linear multistep scheme. We obtain an implicit timediscretization scheme if β 0 = 0. In this paper we use a uniform time step ∆t and a fixed k for each time instance. We refer to Eq. (2.2) as the FOM O∆E.
The time-continuous solution field is a function of both space and time. Thus, the full solution lives in the tensor product manifold of both space and time, i.e., where H is the set of functions from {t n } Nt n=0 to R. We refer to R Ns ⊗ H as the 'FOM trial space'. Let g be an isomorphic function that 'unrolls' time from H to R Nt such that where p ∈ N is an arbitrary dimension. This yields an equivalent (discrete) representation of the FOM space as g(y(·; µ)) = y( Let us also define a vectorization function which converts the matrix solution to a single column vector 3 Projection-based linear model reduction We now consider applying projection-based model reduction for the full order models presented in Section 2: Section 3.1 addresses steady problems and Section 3.2 addresses unsteady problems.

Proper orthogonal decomposition for steady problems
We begin by prescribing the subset of the original state space R Ns on which the ROM techniques will seek approximate solutions. Assume that we have constructed reduced bases Φ = [φ 1 , . . . , φ p ] ∈ R Ns×p with p ≤ N s , where R m×n denotes the non-compact Stiefel manifold: the set of full-column-rank m × n realvalued matrices. (We shall describe the popular snapshot-POD approach for constructing such bases later in this section.) We then approximate the solution in the associated p-dimensional trial subspace with wherex(µ) = [x 1 (µ), . . . ,x p (µ)] T ∈ R p is the generalized/reduced coordinates (or coefficients), x ref (µ) is a fixed "particular solution" or also referred in some literature as "static correction". (We refer to [24] for a specific way to compute x ref (µ).) Remark 1. Given the reduced bases Φ ∈ R Ns×p , we approximate the input-output mapping µ ∈ R nµ → x(µ) ∈ R p using different machine learning regression models, such as polynomials, k-nearest neighbors, random forests, and dense neural networks that will be described in Section 4.
Basis construction This section describes how the reduced bases Φ can be constructed assuming that a full-system state-snapshot matrix X : ⊆ D has been precomputed during an "offline" training stage. Algorithm 1 lists the widely-used proper orthogonal decomposition algorithm that we employ to construct all spatial reduced bases in this work. Throughout, υ ∈ [0, 1] is the "energy criterion" used to determine the applied truncation.

Space-time projection-based model reduction for unsteady problems
As in the previous section, we begin by formulating the reduced basis solution, followed by the construction of the space-time basis. Let Φ = [φ 1 , . . . , φ p ] ∈ R Ns×ns represent the spatial basis matrix and Ψ = [ψ 1 , . . . , ψ k ] ∈ R Nt×nt be the temporal basis matrix, then each spatiotemporal basis is contained in the span of the tensor product space π := Φ ⊗ Ψ = {φ i ⊗ ψ j } ns,nt i,j=1 ∈ R NsNt×nst where n st = n s n t . While this tensor product construction of space-time basis is most straightforward, in this paper we have chosen a tailored set of spatiotemporal bases, which is slightly different, but without loss of generality, can be represented in a similar tensor product form.
Following the idea in [43], we seek the approximated numerical solutionx(µ) to reside in an affine space-time trial subspacex wherex i (µ) ∈ R, 1 ≤ i ≤ n st are the generalized coordinates (or coefficients), π i ∈ R NsNt are the spatiotemporal bases in Φ ⊗ Ψ, and x ref (µ) is a reference solution with possible µ dependence, commonly referred to as the "static correction" or centered solution (over both space and time). In the case where the reference solution is the initial condition, then Here, the key idea of (3.2) is that the time dependence of the approximated solution is moved from the generalized coordinates to the basis vectors; this enables fewer generalized coordinates to be computed in order to characterize the complete space-time solution. We emphasize again that this is quite new as most of works in literature proposed to use parameter-and time-dependent generalized coordinates (i.e., spatial bases) [24,37,38,39,40,41]; and very few works used only parameter-dependent generalized coordinates (i.e., space-time bases) [42,43].
. ,x nst (µ)] T ∈ R nst using different machine learning regression models, such as polynomials, k-nearest neighbors, random forests, and dense neural networks that will be described in Section 4.

Space-time basis construction
In order to compute the space-time basis, we evaluate the FOM (2.1) at the set of training parameter instances D train := {µ 1 train , . . . , µ ntrain train } defined by, e.g., uniform sampling, Latin-hypercube sampling, greedy, etc., using a linear multistep method (2.2). This gives us 'snapshots' of the solution space x(t j ; µ k train ) ∈ R Ns , j ∈ N(N t ), k ∈ N(n train ) in both space and time; here N(N ) := {1, . . . , N }. This data collection procedure is commonly referred to as the 'offline' stage in model reduction.
The space-time trial subspace is computed by applying tensor-decomposition techniques to the three-way 'state tensor' X ∈ R Ns×Nt×ntrain with elements We refer readers to [43] for a complete research of such space-time bases. We expect our resulting ROM be accurate if the solution exhibits separable behavior in space and time.
In order to construct the basis from this three-way tensor, we introduce two different unfoldings of X . These unfoldings or reshapings allow us to utilize standard matrix decomposition techniques such as PCA (Principal Component Analysis) with simplicity. The mode-1 and mode-2 unfoldings of X are defined as In short, X (1) removes time as a variable by treating time and the parameter space the same, and X (2) removes the spatial dimension as a variable by treating space and parameter samples the same. In the modelreduction literature, the matrix X (1) is typically referred to as the 'global snapshot matrix'; we refer to it in this work as the 'spatial snapshot matrix' as its columns comprise snapshots of the spatial solution over time and parameter variation. Similarly, we refer to X (2) as the 'temporal snapshot matrix' as its columns comprise snapshots of the time-evolution of the solution over variation in space and parameter. Now we can begin computing the spatiotemporal basis vectors.
Alternatively, we propose using the "tailored temporal bases" defined in [43] due to its advantageous capabilities in capturing the time evolution of each individual spatial basis vector. In particular, we use ST-HOSVD (sequentially truncated high-order SVD) [47,48]. The idea can be summarized as follows. First, for each spatial basis vector φ i we project each space-time snapshot X(µ m train ), m ∈ N(n train ) onto φ i , i ∈ N(n s ). This computation can be performed by a single matrix-vector multiplication to create X(φ i ) (2) = X (2) φ i ∈ R Nt×ntrain . Then, in order to obtain the n i t -dimensional tailored temporal bases corresponding to φ i , we simply apply an SVD, as in Algorithm 1, giving us a collection of tailored temporal basis vectors Ψ i ∈ R Nt×n i t , where n i t may be different for each spatial basis vector φ i . We then repeat this for all n s spatial basis vectors to get a total of ns i=1 n i t temporal basis vectors. Finally, the Kronecker product is implemented corresponding to each spatial basis vector to create resulting spatiotemporal bases π = {φ i ⊗ Ψ i } ∈ R NsNt×nst , i ∈ N(n s ). Detailed implementations are listed in Algorithm 2.
There are several notes for this approach described as follow. First, the maximum dimension of each temporal basis is limited by the number of training-parameter instances, i.e., n i t ≤ n train , i ∈ N(n s ). Second, this approach is less computationally expensive than applying directly SVD on X (2) above since we shall perform n s SVD algorithms on X(φ i ) (2) ∈ R Nt×ntrain , while X (2) ∈ R Nt×Nsntrain and typically n s N s . Third, Ref. [43] has shown that comparing with other temporal basis approaches, this approach is very efficient since it provides similar accuracy levels but using far fewer degrees of freedom (in the context of GNAT hyper-reduction technique) as each temporal basis vector is tailored to its associated spatial basis vector. Finally this was the approach employed to construct temporal bases in previous works based on forecasting [49,50]. Input: Spatial snapshot matrix X (1) ∈ R Ns×Ntm , energy criterion υ ∈ [0, 1] Output: Spatiotemporal basis matrix π ∈ R NsNt×p n i t 1: Compute (thin) singular value decomposition (SVD) for X (1) Note that here we generalize so that the number n i t can be adapted for each spatial vector, but for simplicity, we will assume a single n t for a total of n s n t spatiotemporal basis vectors.

Machine learning surrogate models
Four different machine learning methods were used to approximate the mappings µ →x(µ) in Remark 1 and Remark 2, respectively. This section describes the setup and brief overview of each machine learning method. Specific implementation choices for each modeling approach are given in Section 5.

Learning problem setup
In general, a surrogate model is constructed for the mapping α(µ) : D → A , e.g., α ≡x and A ⊆ R p (or A ⊆ R nst ), respectively. The outputs α(µ) are the coefficients defined in Eq. (3.1) or Eq. (3.2) and the inputs are system parameters. In general, let µ ∈ D ⊆ R nµ and A ⊆ R r .
Consider the case that we have m snapshots, where each snapshot corresponds to a different input parameter. (Note that for steady problem, a snapshot is a spatial FOM solution while for unsteady problem, a snapshot is a space-time FOM solution. Snapshots correspond to different input parameters.) The input parameters corresponding to each snapshot are collected in the matrix D ∈ R m×nµ ; while the corresponding basis coefficients for each snapshot are collected in the matrix A ∈ R m×r , respectively. Subsequent sections provide an overview of the four methods used to make predictions. For each of these methods, input and output data D and A are divided into training and testing sets denoted by (D train , A train ) and (D test , A test ), where m = n train + n test respectively. The goal is to learn the mapping α : D → A from the training data (D train , A train ), and test the performance of the ML methods on (D test , A test ).
With the available snapshots data x(µ) and constructed reduced basis Φ (or π, respectively), we can obtain the "best" linear ROM approximations by projecting all FOM solutions onto the reduced bases and then compute associated "best" relative error as: wherex best (µ) ∈ R r are the "best" generalized coordinates (in linear ROM sense). The training and testing data is then defined as:

Polynomial regression
The first model considered is multivariate polynomial regression (MPR), which approximates the outputs (or functions) α(µ) using a multivariate polynomial expansion and least squares regression. In this paper, we use polynomial functions up to fifth order to define the model. The model for a coefficient α l (µ) can be written as where b i , i = 0, 1, . . . , n µ ; c ij , i, j ∈ N(n µ ) and d ijk , i, j, k ∈ N(n µ ) are the coefficients of the approximation model, and µ i , µ j , µ k are components of µ. Note that in this type of expansion, known as a full tensor product grid, the number of basis elements scales like s nµ where s is the polynomial degree. In the offline/training stage, these unknown coefficients are determined via least squares regression that minimizes the mean squared error between the predicted and actual training outputs. Once these coefficients are determined, in the online/evaluation stage one can quickly compute α(µ test ) for any given µ test . (Note that in all numerical experiments in Section 5, the number of training samples are always larger than number of unknown coefficients assuring that the least square problem is well-defined.) The largest computational complexity arises from solving the least squares systems, i.e. inverting the D train T D train covariance term. The computational complexity depends on n train , n µ and the degree of the regression polynomials. As shown in Section 5, the polynomial regression model is much cheaper to train and to evaluate than other more complicated ML models (such as neural network or random forest regressors). 2

k-nearest neighbor regression
The k-nearest-neighbors (kNN) model can be interpreted as a localized form of multivariate regression: the approximation is computed over a local set of k training samples that are the closest neighbors of the evaluation point µ in the input space. Due to the local nature of kNN approximations, the resulting models do not provide a global interpretation of the underlying relationship between inputs and outputs. Therefore kNN is an unstructured but flexible scheme, balancing the simplicity of polynomial regression with the flexibility of localized models.
In this paper, we consider a standard implementation of kNN that approximates α(µ) by averaging the outputs associated with the k nearest neighbors. In the offline stage, a k-d tree partitions the data along each input dimension and can be constructed in O(n µ n train ) time without explicitly calculating any n µdimensional distances. This also allows fast nearest neighbor searches that have an average computational complexity of O(log(n train )) [51,52]. In the online stage, the predicted output of the kNN model is the weighted average of the training outputs at each of the neighbors; this is essentially a localized linear regression model.

Random forest regression
The third model considered is random forest regression, which is an ensemble of decision trees, generally trained via the bagging method (or pasting method) [53]. A decision tree partitions the input domain into many regions and makes local average estimates of the output. This is performed in a top down fashion where the data are recursively partitioned at each node using a greedy algorithm to group similar data together. Each node is assigned a certain split criteria, which aims to make a split that creates similar valued subsets of data. Nodes in the tree continue to split the data into left and right child nodes using a series of logical statements. This process continues until certain stopping criteria are satisfied and nodes cease to split, becoming leaf nodes. After the data are partitioned into multiple regions, a piecewise constant regression model is built. Within each region (i.e., leaf node), a constant function is fit as the average of each value in the region.
Many efficient methods for constructing decision trees are available [54]. Once a tree has beed learned, predictions can be made rapidly. Given an input, only simple logic statements are processed to reach a leaf node where an average value is computed to obtain the output. Partitioning the training data along each input dimension has a computational complexity of O(n µ n train log(n train )) to build the tree. If a decision tree is approximately balanced, this allows fast predictions on the order of O(log(n train )). One of the disadvantages of a random forrest regression model is that it is not globally differentiable, as opposed to a polynomial model which is infinitely differentiable.

Neural network regression
Input Layer Hidden Layer Hidden Layer Output Layer The last model considered is a fully connected, feed forward neural network, or multi-layer perceptron model. Neural network models estimate the output α(µ) at an input µ using weights and biases that are adapted to the data during training. With the use of adaptivity, these models are structured but can achieve great flexibility. A disadvantage of using neural networks is that interpretation of the resulting model may be difficult.
A neural network consists of sequential layers of nodes that define a mapping between an input and an output. The basic structure of a neural network is shown in Figure 1. The first layer of nodes is called the input layer, followed by intermediate hidden layers, and ended with an output layer. The input layer directly receives the input to the model: the number of nodes in the input layer is equal to the number of inputs. Next, in each hidden (i.e., not an input or output) layer l, there are n l nodes. The ith node in layer l is denoted as η i l for i ∈ {1, . . . , n l }. Each node takes an input, evaluates an activation function g, and produces an intermediate output o i l , which is used as an input to nodes in the next layer l + 1. Every node in a given layer is connected to every node in the following layer (fully connected) and information is passed forward through the network (feed forward). Each of these connections are given a weight w i l and a bias b i l , which defines the importance and the effect of a particular input to the output [55]. The activation function at each node is a function of h i l = w i l o i l−1 + b i l and thus depends on the input to the node, the weight assigned to the connection to that node and the bias. The output for node η i l is determined by evaluating the activation The last layer is the output layer which produces the final output: the number of nodes in the last layer is equal to the number of output values.
To train a neural network, training inputs D train and corresponding training outputs A train are provided to the model. The training process iterates forward and backward over the network, adjusting weights and biases to minimize the mean squared error between predicted and actual training outputs [56]. Each pair of forward and backward passes is referred to as an epoch. The computational complexity of training a neural network is approximately linear with the size of the network (nodes and layers), which in turn depends upon the dimension of the input and output, and with the number of epochs, assuming using a smooth activation function [57].

Numerical experiments
This section reports numerical experiments that assess the performance of the proposed space-time machine learning ROM (ST-ML-ROM) method on three benchmark problems: a steady heat PDE, an unsteady heat PDE and an advective-diffusive-reactive problem. In particular, the first numerical problem is to check our methods against established results [24], the second problem is to extend to unsteady problems, and the third problem is to check feasibility when there are many spatial/temporal scales and modes are depend nonlinearly on input parameters. We focus on the accuracy and speed of four regression methods (polynomial, kNN, random forest and neural network) and compare the following models: • FOM. This model corresponds to the full-order model, i.e., the solution satisfying Eq. (2.1).
• ST-ML-ROM. This model corresponds to data-driven space-time ROMs, i.e., surrogate models constructed via using each ML technique described in sections 4.2-4.5.
For steady and unsteady problems, the accuracy of ROM solutionsx(µ) are measured by computing the 2 -norm of relative errors: Intel Core i9 with 16 GB RAM of memory. Reported timings comprise the average time over specified numbers of repeated runs for offline and online stages of each machine learning method used. We also note that the implementation of all ML methods in this paper are pretty straightforward with standard settings and there is no special treatments such as regularization or any ad hoc tuning.  We first consider the model example introduced in Refs. [59,13,46]. This is a parametric nonlinear 2D heat problem which consists of computing u(x, µ) with x ≡ (x 1 , x 2 ) ∈ Ω = [0, 1] 2 and µ ≡ (µ 1 , µ 2 ) ∈ D = [0.01, 10] 2 and homogeneous Dirichlet boundary condition on Γ ≡ ∂Ω satisfying

Full-order model
This model can be interpreted as a 2D stationary diffusion problem with a nonlinear interior heat source. The resulting solution exhibits a strongly nonlinear dependence on the parameters µ. Note that we use this steady problem to verify all machine learning methods developed earlier (such as the work [24]). For spatial discretization, we apply the finite-element method using a mesh characterized by 3600 (60x60) bilinear quadrilateral (Q1) elements. Figure 2a depicts the mesh, while Table 1 reports the corresponding discretization parameters. Figure 2b plots the FOM reference solution on the mesh with an input parameter µ = (10, 10). Applying these finite-element discretizations to Eq. (5.3) leads to a parameterized system of nonlinear algebraic equations of the form (2.1) withẋ = 0.

Model reduction by snapshot-POD
We use Latin Hypercube Sampling with a uniform distribution [60] to generate a total of 300 sample points for the input parameters µ in D that is denoted as D. These 300 sample points are split into 200 training samples (D train ) and 100 testing samples (D test ), where D train is used for training purpose and D test is used for testing purpose. Figure 3a depicts Latin Hypercube Sampling with a uniform distribution, while Figure 3b plots the training and testing samples sets (D train and D test ), respectively. We also note that To generate the reduced bases required for the reduced-order models, we solve the FOM (2.1) for all µ ∈ D train , i.e., X train = [x(µ 1 train ), . . . , x(µ ntrain train )], then feed X train to Algorithm 1 to obtain reduced bases. Figure 4 shows the results of Algorithm 1: Figure 4a, 4b plot the singular values and corresponding cumulative energy for spatial SVD. Figure 4b implies that we only need first 5 POD bases to capture more than 99.99% energy of the system, hence we choose p = 5 for all computations afterward. Note that we also solve the FOM (2.1) for all µ ∈ D test , i.e., X test = [x(µ 1 test ), . . . , x(µ ntest test )] to compute the relative solution error of all proposed ML-ROM techniques later.
Next, we can obtain the "best" linear ROM approximations by projecting all FOM solutions onto the POD bases and then compute associated "best" relative errors following the computation in Eq. (4.1). Figure 5 plots "best" relative solution errors with µ ∈ D.    We now use the four machine learning techniques described in Section 4 to approximate the mapping µ →x best (µ) by µ →x(µ) in the offline (training) stage; then compute the approximated ROM solution Eq. (3.1) and associated relative error Eq. (5.1) in the online (prediction) stage. Note that D train is used to build surrogate models and D test is used to assess the performance (accuracy and speed) of these surrogate models as described in Section 4.1.

ML-ROM: one online computation
We compare the four methods for fixed values of their parameters, and for a single randomly selected online point µ test = (9.0350, 9.2472) ∈ D test . In particular, we use 2 nd order polynomial regressor, kNN regressor with 3 neighbors, random forest regressor with configuration (max features=2, max leaf nodes=16, n estimators=271) and neural network with configuration (n hidden=2, n neurons=50, learning rate=10 −3 ,    Table 2 at point µ test = (9.0350, 9.2472) ∈ D test . activation=tanh), respectively 3 . Note that the configurations and hyper-parameters of these regressors are chosen randomly to check their performance, we will perform a more rigorous parameter study in the next Section 5.1.4. Table 2 reports the chosen hyper-parameters and associated performance of the methods, while    Table 2 at point µ test = (9.0350, 9.2472) ∈ D test .  Figure 7 visualize the ML-ROM solutions and solution errors of all methods, respectively. The results in Table 2 and Figures 6, 7 verify that the proposed methods yield accurate solutions as anticipated [24].

ML-ROM: parameter study
Because assessing a given method's performance for a single instance of its parameters does not lend insight into the model's complete error-cost performance tradeoff, this section subjects each of the proposed ROMs to a parameter study wherein each model parameter is varied within relatively wide ranges. In this spirit, the degree of polynomial regressor will be varied in the range 1 ≤ r ≤ 5. For kNN regressor, the number of neighbors will be varied in the range 1 ≤ k ≤ 19. For random forest regressor, we first perform a cross-validation randomized search [61] with hyper-parameter limits specified in Table 3. Table 4 reports hyper-parameter values associated with nine best configurations of random forest regressors resulting from the tanh sigmoid relu tanh relu sigmoid relu relu tanh randomized search. For neural network regressor, we also perform first a cross-validation randomized search [61] with hyper-parameter ranges specified in Table 5. Table 6 reports hyper-parameter values associated with nine best configurations of neural network regressors after the randomized search. We then investigate the accuracy of these configurations by computing their relative solution errors for all µ ∈ D test . The computational results are shown in Figure 8. In particular, the relative errors are presented in the form of data boxes where on each box the median, 25 th , 75 th percentiles, most extreme data points and outliers are all shown together. For reference, the "best" relative error defined in (4.1) is also added to Figure 8 as the lowermost data box. Figure 8 shows that the proposed methods yield accurate solutions and reasonable error levels for engineering problems.  We also examine the computational time of the considered configurations. To quantify the computation time, the offline and online stages of each configuration of the four methods are repeated many times and the reported time is averaged from that of these repeated runs. Table 7 lists the number of repeated runs mentioned previously. We emphasize that the reported time only measures the time to compute µ →x(µ), and does not includex(µ) →x(µ) as the stepx(µ) →x(µ) is the same for all regressors. Figure 9a and 9b report the training time and prediction time of some configurations of four regressors. In general, we observe that neural network takes much more time to train and predict than the other three techniques.
From the results in Figures 8 and 9b, we construct a Pareto front to characterize the error-cost for each method. This Pareto front is characterized by the method hyper-parameters that minimize the competing objectives of relative error and computational time. Figure 9c reports the Pareto fronts of the four methods considered. From Figure 9c, we observe a very similar trend to that in the work [24] in which simple ML techniques (such as polynomial, kNN and random forest regressors) outperform complex ML ones (such as   neural network) in terms of accuracy and computational cost (for this specific example). For example, with the same accuracy level neural network is about two order of magnitude slower than polynomial and kNN regressors, and one order of magnitude slower than the random forest regressor. The conclusion is somewhat similar to that in [24]: in engineering problems where the training data is often expensive to obtain and sparse, simple ML techniques can outperform the complex ones because the limitation of training data is limiting the performance of the complicated ML techniques.

Accuracy comparison with LSPG-ROM
We finally provide an accuracy comparison of the proposed methods with the well-known LSPG-ROM method [44,11,45]. The LSPG-ROM method is implemented using our in-house code [46] with standard settings. Figure 10 shows all relative solution errors of proposed methods and the LSPG-ROM method, together with the "best" relative errors (4.1). (In other words, Figure 10 comprises Figure 8 adding relative errors of LSPG-ROM method.) From Figure 10, we observe that LSPG-ROM provides very accurate solutions that are very close to the "best" possible linear ROM solutions; and the LSPG-ROM solutions are much more accurate than that of the proposed methods. However, note that the computational cost of LSPG-ROM  Figure 10: Steady heat equation, relative solution error with µ ∈ D test . The lowermost box is for "best" relative error defined in 4.1, the middle box is for LSPG-ROM (green color) [44], and all upper boxes for ML-ROM using four machine learning techniques.
(without hyper-reduction) is much higher than that of our proposed methods 4 .

Full-order model
We now consider the model example introduced in Refs. [59,13]. This is a parametric nonlinear 2D unsteady heat problem which consists of computing u(x, t, µ) with x ≡ (x 1 , x 2 ) ∈ Ω = [0, 1] 2 , t ∈ [0, T ] with T = 2 and µ ≡ (µ 1 , µ 2 ) ∈ D = [0.01, 10] 2 and homogeneous Dirichlet boundary condition on Γ ≡ ∂Ω, such that: where ρ = 10 4 is a constant, b = sin(2πt) is the control input. This model can be interpreted as a 2D dynamically diffusion problem with a nonlinear interior heat source density. It presents nonlinear dependence of the solution with respect the input parameter µ.
For spatial discretization, we use a finite element mesh with 3600 (60×60) bilinear quadrilateral (Q1) elements, which results in a system of parameterized ODEs of the form (2.1). The same Table 1 also reports the spatial discretization parameters used for this model. For time discretization, we employ the first order backward Euler scheme which is a linear multistep method characterized by k(t n ) = 1, α n 0 = β n 0 = 1, We use a uniform time step ∆t = 8 × 10 −3 , leading to N t = 250 time instances. Figure 11a depicts the mesh, while Figure 11b plots the FOM reference solution at time instance t = 1.5s with µ = (3, 3), respectively.   Similar to section 5.1.2, we use Latin Hypercube Sampling with a uniform distribution [60] to generate a total of 300 sample points for the input parameter µ in D that is denoted as D. These 300 sample points are split into 200 training samples (D train ) and 100 testing samples (D test ), where D train is used for training purposes and D test is used for testing purposes. Figure 12a depicts Latin Hypercube Sampling with a uniform distribution, while Figure 12b plots the training and testing samples sets (D train and D test ), respectively. We also note that D = D train ∪ D test and D train ∩ D test = ∅.

Model reduction with space-time bases
To generate the space-time reduced bases required for the reduced-order models, we solve the FOM (2.2) for all µ ∈ D train , i.e., X train = [x(t n , µ 1 train ), . . . , x(t n , µ ntrain train )], n ∈ N(N t ) then feed X train to Algorithm 2 to obtain the space-time bases. Figure 13 shows the results after performing this step. In particular, Figure 13a, 13b present singular values and cumulative energy for spatial SVD, while Figure 13c and 13d plot singular values and cumulative energy for temporal SVD associated with Φ s 2 . From Figures 13b and  13d, we choose n s = 3 spatial bases, n t = 2 temporal bases, hence n st = n s n t = 6 space-time reduced bases for all computations afterward. Again, we also solve the FOM (2.2) for all µ ∈ D test , i.e., X test = [x(t n , µ 1 test ), . . . , x(t n , µ ntest test )], n ∈ N(N t ) to compute the relative solution error of all proposed ML-ROM techniques. We finally obtain the "best" linear ROM approximations by projecting all space-time FOM solutions onto the space-time bases and then compute associated "best" relative errors following Eq. (4.1).  Figure 14: Unsteady heat equation, "best" relative error for µ ∈ D total with n s = 3, n t = 2, p = n s n t = 6. Figure 14 plots "best" relative solution errors with µ ∈ D using n s , n t and n st chosen above.
We now use the four machine learning techniques described in Section 4 to approximate the mapping µ →x(µ) in the offline (training) stage; then compute the approximated ROM solution Eq. (3.2) and associated relative error Eq. (5.2) in the online (prediction) stage. Note again that D train is used to build surrogate models and D test is used to assess the performance (accuracy and speed) of these surrogate models. We now compare the four methods for fixed values of their parameters, and for a single randomly selected online point µ test = (9.6289, 7.2334) ∈ D test . In particular, we use 2 nd order polynomial regressor, kNN regressor with 3 neighbors, random forest regressor with configuration (max features=2, max leaf nodes=44, n estimators=162) and neural network with configuration (n hidden=4, n neurons=75, learning rate=3.74× 10 −4 , activation=sigmoid), respectively. Note that the configurations and hyper-parameters of these regressors are chosen randomly to check their performance, we will perform a more rigorous parameter study in the next Section 5.2.4. Table 8 reports the chosen hyper-parameters and associated performance of the methods, while Figure 15 and Figure 16 visualize the ST-ML-ROM solutions and solution errors of all methods at time instance t = 1.6s, respectively. The results on Table 8 and Figures 15, 16 verify that the proposed methods yield very accurate solutions.  Table 8 at time t = 1.6s, µ test = (9.6289, 7.2334) ∈ D test .

ST-ML-ROM: parameter study and comparison with space-time LSPG-ROM
We again compare the performance of the ROM methods across a wide variation of all method parameters. In particular, the degree of polynomial regressor will be varied in the range 1 ≤ r ≤ 5. For kNN regressor, the number of neighbors will be varied in the range 1 ≤ k ≤ 19. For random forest and neural network regressors, we employ the same approach as described in Section 5.1.4 to find their several optimal configurations.  Table 8 at time t = 1.6s, µ test = (9.6289, 7.2334) ∈ D test .
Since the input-output mapping of this problem (R 2 → R 6 ) is very similar to that of the previous problem in Section 5.1 (R 2 → R 5 ), we use the same randomized search ranges for hyper-parameters of random forest and neural network regressors as specified in Tables 3 and 5, respectively. Tables 9 and 10 report hyperparameter values associated with nine best configurations of random forest and neural network regressors resulting from the corresponding cross-validation randomized searches.
We then investigate the accuracy of these configurations by computing their relative solution errors together with the relative solution errors of the space-time LSPG-ROM method [43] and "best" relative error (4.1) for all µ ∈ D test . Figure 17 reports all these relative solution errors. We observe a similarity between Figure 17 and Figure 10 that is the space-time LSPG-ROM provides very accurate solutions that are close to the "best" possible linear ROM solutions; and the space-time LSPG-ROM solutions are slightly more accurate than our proposed methods. Again, we also note that the computational cost of space-time LSPG-ROM (without hyper-reduction) is much higher than that of our proposed methods (see Section 5.1.5).
The computational time of the considered configurations are also examined. We apply the same procedure described in Section 5.1.4 to measure computational times; the same Table 7 lists the number of repeated offline and online runs. Figure 18a and Figure 18b report the training time and Pareto front for online predictions, respectively. Again, Figure 18b shows a very similar trend to Figure 9c in which simple ML techniques (such as polynomial, kNN and random forest regressors) outperform complex ML ones (neural network) in terms of accuracy and computational cost for this problem.

Results with 75% training data
We investigate the effect of the amount of training data on computational results. In particular, within training data set D train , we pick a random nonstratified partition [62] that has 75% size of D train to train all ML models, then evaluate the accuracy of these models similarly to Section 5.2.4 above. Figure 19 reports the relative solution errors together with the "best" relative error (4.1) for all µ ∈ D test . Comparing Figure 19 Figure 17: Unsteady heat equation, relative solution error with µ ∈ D test . The lowermost box is for "best" relative error defined in Eq. (4.1), the middle box is for space-time LSPG-ROM (green color) [43], and all upper boxes for ST-ML-ROM using four machine learning techniques.
and Figure 17 shows that all ML algorithms still work well and provide accurate solutions even when using only 75% training data. (It is generally obvious that using much less training data will have adverse effect on the performance of ML algorithms.)

Model reduction with space-time bases
We again use Latin Hypercube Sampling with uniform distribution to generate a total of 400 sample points for the input parameter µ in D that is denoted as D. These 400 sample points are split into 320 training samples (D train ) and 80 testing samples (D test ), where D train is used for training purpose and D test is used for testing purpose. Within D train , we use only first 150 points to build space-time bases due to implementation To generate the space-time reduced bases required for the reduced-order models, we solve the FOM (2.2) for all µ ∈ D train , i.e., X train = [x(t n , µ 1 train ), . . . , x(t n , µ ntrain train )], n ∈ N(N t ) then feed X train to Algorithm 2 to obtain the space-time bases. Figure 22 shows the results after performing this step: Figure 22a, 22b present singular values and cumulative energy for spatial SVD, while Figure 22c and 22d plot singular values and cumulative energy for temporal SVD associated with Φ s 2 . Based on Figures 22b, 22d and some trial-error experiments, we choose n s = 20 spatial bases, n t = 3 temporal bases, hence n st = n s n t = 60 space-time reduced bases for all computations afterward.
Again, we also solve the FOM (2.2) for all µ ∈ D test , i.e., X test = [x(t n , µ 1 test ), . . . , x(t n , µ ntest test )], n ∈ N(N t ) to compute the relative solution error of all proposed ST-ML-ROM techniques. We finally obtain the "best" linear ROM approximations by projecting all space-time FOM solutions onto the space-time bases and then compute associated "best" relative errors following Eq. (4.1). Figure 23 plots "best" relative solution errors with µ ∈ D using n s , n t and n st chosen above. We then build the surrogate models (using four machine learning techniques described in Section 4) for the mapping µ →x(µ). The surrogate models are built in the offline (training) stage using data in D train ; and their performance (accuracy and speed) are evaluated in the online (prediction) stage using data in D test , respectively.

ST-ML-ROM: one online computation
We now compare the four methods for fixed values of their parameters, and for a single randomly selected online point µ test = (5.9315 × 10 12 , 3.5416) ∈ D test . In particular, we use 2 nd order polynomial regressor, kNN regressor with 3 neighbors, random forest regressor with configuration (max leaf nodes=116, n estimators=81) and neural network with configuration (n hidden=10, n neurons=19, learning rate=1.27×  10 −3 , activation=sigmoid), respectively. Table 11 reports the chosen hyper-parameters and associated performance of the methods, while Figure 24 visualizes the ST-ML-ROM solutions and solution errors of the 2nd order polynomial regressor. The results on Table 11 and Figure 24 verify that our proposed methods  Figure 23: Reacting flow, "best" relative error for µ ∈ D with n s = 20, n t = 3, p = n s n t = 60.
yield very accurate solutions.   We again compare the performance of the ROM methods across a wide variation of all method parameters. In particular, the degree of polynomial regressor will be varied in the range 1 ≤ r ≤ 5. For kNN regressor, the number of neighbors will be varied in the range 1 ≤ k ≤ 19. For random forest and neural network regressors,   Table 11.

ST-ML-ROM: parameter study
we employ the same approach as described in Section 5.1.4 to find their several optimal configurations. In this spirit, the hyper-parameter limits of the cross-validation randomized searches for random forest and neural network regressors are listed on Table 12 and 14; while Table 13 and 15 report the hyper-parameter  values associated with their nine and seven best configurations after the randomized searches, respectively. We then investigate the accuracy of these configurations by computing their relative solution errors together with the "best" relative error Eq. (4.1) for all µ ∈ D test . Figure 25 reports all these relative solution errors. We observe from Figure 25 again that the proposed methods yield accurate solutions and reasonable error levels for engineering settings.
Finally, the computational time of the considered configurations are also investigated. We apply the same   Table 7 lists the number of repeated offline and online runs for each configuration. Figure 26a and Figure 26b report the training time and Pareto front for online predictions, respectively. Again, Figure 26b shows a very similar trend to Figure 18b and 9c in which simple ML techniques (such as polynomial, kNN and random forest regressors) outperform complex ML ones (neural network) in terms of accuracy and computational cost; especially random forest is the best choice for this complicated problem.

Conclusion
Numerical experiments in this paper have demonstrated that the use of space-time subspace bases combined with machine learning methods provide a tractable and effective way to map between input parameters to a time-dependent high-dimensional output quantity of interest. The method can be applied straightforwardly to applications that exhibit a fast decaying Kolmogorov n-width 6 (see other approaches, e.g., [63,64,65] to handle problems exhibiting slowly decaying Kolmogorov n-width). The numerical experiments in this paper have also highlighted the important point that the availability of data in an engineering setting is typically much less than it is in other machine learning applications. This is because engineering data are often generated from expensive physics-based codes or from sensors embedded on a physical system. In either case, it is reasonable to expect hundreds or thousands of training data points, but not millions as is the case in many non-engineering machine learning applications. The strategy to choose appropriate machine learning methods is thus affected by this fact. In particular, we have observed that in this setting simpler polynomial regression or random forest models may be better choices (than complicated neural networks)even when the underlying physical phenomena exhibit strongly nonlinear behavior. Another possible reason explaining why polynomial regression, kNN and random forest provide more accurate solutions than neural network does (in this engineering setting) is that these 3 methods build one surrogate model for each mapping µ → α i (µ), i ∈ N(r) (hence totally r models to approximate whole mapping), while neural network build one unique surrogate model for the whole mapping µ → α(µ).