{'fraction_non_alphanumeric': 0.02672834746851709, 'fraction_numerical': 0.005718324338216397, 'mean_word_length': 4.894527551600076, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Latent-space time evolution of non-intrusive reduced-order models using Gaussian process emulation

Non-intrusive reduced-order models (ROMs) have recently generated considerable interest for constructing computationally efficient counterparts of nonlinear dynamical systems emerging from various domain sciences. They provide a low-dimensional emulation framework for systems that may be intrinsically high-dimensional. This is accomplished by utilizing a construction algorithm that is purely data-driven. It is no surprise, therefore, that the algorithmic advances of machine learning have led to non-intrusive ROMs with greater accuracy and computational gains. However, in bypassing the utilization of an equation-based evolution, it is often seen that the interpretability of the ROM framework suffers. This becomes more problematic when black-box deep learning methods are used which are notorious for lacking robustness outside the physical regime of the observed data. In this article, we propose the use of a novel latent space interpolation algorithm based on Gaussian process regression. Notably, this reduced-order evolution of the system is parameterized by control parameters to allow for interpolation in space. The use of this procedure also allows for a continuous interpretation of time which allows for temporal interpolation. The latter aspect provides information, with quantified uncertainty, about full-state evolution at a finer resolution than that utilized for training the ROMs. We assess the viability of this algorithm for an advection-dominated system given by the inviscid shallow water equations.


Introduction
Once this basis (or latent representation) is identified, we need a cost-effective strategy for accurate nonlinear dynamical system evolution to reproduce the full-order spatiotemporal complexity of the problem in the reduced basis. For example, linear reduced-basis construction allows for the use of intrusive methods (which project the governing equations onto the reduced-basis), as seen in [15,34], which use a Galerkin projection; or [6,49,11], which use the Petrov-Galerkin method (see [5] for the comparison of these two methods). Such extensions are not straightforward for autoencoder based latent space constructions, since projecting to a basis space is infeasible. We note, though, that the use of convolutional autoencoder architectures has been demonstrated for the Petrov-Galerkin method, where there is no requirement to project the governing equations onto a trial space [24]. As mentioned previously, neural networks are also commonly used for the temporal evolution of these latent space representations. Mainstream ML has generated a large body of time-series forecasting literature that lends itself readily to latent-space evolution of dynamical systems. Recently, long short-term memory architectures (LSTMs) have become popular for the non-intrusive characterization of dynamical systems [45,1,32,33,46] as they allow for the construction of non-i.i.d, directional relationships between the states of a dynamical system at different times before a future forecast is performed. Other methods that have been utilized for such tasks include the temporal convolutional network [50], non-autoregressive variants for the LSTM [29] and system-identification within latent space [7]. In most of these developments, the evolution framework (in time) is deterministic in nature. These frameworks highlight a crucial drawback when coupled with the black-box nature of purely data-driven ROMs -it is imperative to embed some notion of uncertainty into the time evolution of ROMs in the latent space. Past work has shown that LSTMs suffer from stability and interpretability issues when utilized for the reduced-order modeling of advection-dominated systems [29] and this study proposes an alternative to provide interpretable forecasts with quantified uncertainty. To that end, we propose the use of a Gaussian process regression (GPR) framework for evolving the low-dimensional representation of data obtained from a PDE evolution. In addition to providing forecasts at the temporal resolution of the training data, the use of the our framework allows for interpolation in time. This is because time is interpreted as a continuous variable. When coupled with quantified uncertainty, the framework allows a ROM user to interrogate the behavior of the emulated system at a finer temporal resolution, with implications for data sources that are temporally sparse.
To summarize, the major contributions of this article are: • We introduce a technique to construct parametric non-intrusive ROMs with quantified uncertainty during latent-space time evolution of dynamical systems.
• We detail the use of Gaussian processes (GPs) that are customized for the time-evolution of parametric nonlinear dynamical systems.
• We demonstrate the ability of the proposed time-evolution algorithm for systems compressed by linear reduced-basis methods such as POD, as well as nonlinear compression frameworks such as variational and convolutional autoencoders.
• We test the ability of a ROM constructed from coarse training data to interpolate on a finer temporal resolution with quantified uncertainty. In the following, we shall introduce our experimental setup which relies on the inviscid shallow water equations in Section 2, our various compression mechanisms in Section 3, the UQ embedded emulation strategy in Section 4, followed by results and conclusions in Sections 5 and 6 respectively.

The shallow water equations
The inviscid shallow water equations belong to a prototypical system of equations for geophysical flows. In particular, the shallow water equations admit solutions where advection dominates dissipation and pose challenges for conventional ROMs [47]. These governing equations are given by where, η corresponds to the total fluid column height, and (u, v) is the fluid's horizontal flow velocity, averaged across the vertical column, g is acceleration due to gravity, and ρ is the fluid density, typically set to 1.0. Here, t, x and y are the independent variables of time and the spatial coordinates of the two-dimensional system. Equation 1 captures the law of mass conservation, whereas Equations 2 and 3 denote the conservation of momentum. The initial conditions of the problem are given by ρηu(x, y, t = 0) = 0 (5) ρηv(x, y, t = 0) = 0, i.e., a Gaussian perturbation at a particular location on the grid [x,ȳ] ≡ w. We solve the system of equations until t = 0.5 with a time-step of 0.001 seconds on a square two-dimensional grid with 64 collocation points to completely capture the advection and gradual decay of this perturbation. Note that these numbers may vary according to the forecasting and fidelity requirements of a particular problem and perturbation. The initial and boundary conditions for this particular shallow-water equation experiment represent a tightly-controlled traveling wave problem that is translation invariant. Different realizations of the initial condition lead to translationally shifted trajectories. We also note the presence of mirror symmetries with respect to x =x and y =ȳ coupled with a rotational symmetry of π/2 radians about the origin. However, our motivation for a first assessment of our emulators on this system stems from the well-known fact that POD and Galerkin-projection based methods are severely limited in their ability to forecast on these simple traveling-wave systems [22,31] and require special treatment with intrinsic knowledge of the flow dynamics. This is in addition to the fact that equation-based models are impossible to construct because of the absence of information from the other variables of the PDE. We seek to build predictive models solely from observations of ρη conditioned on w mimicking a real-world scenario where complete observations of all relevant variables (in this case, velocities) are unavailable.

Proper orthogonal decomposition
In this section, we review the POD technique for the construction of a reduced basis [20,2]. The interested reader may also find an excellent explanation of POD and its relationship with other dimension-reduction techniques in [42]. The POD procedure is tasked with identifying a space which approximates snapshots optimally with respect to the L 2 -norm. The process of ϑ generation commences with the collection of snapshots in the snapshot matrix where N s is the number of snapshots, andq i h : T × P → R N h corresponds to an individual snapshot in time of the discrete solution domain with the mean value removed, i.e., with q h : P → R N h being the time-averaged solution field. Our POD bases can then be extracted efficiently through the method of snapshots where we solve the eigenvalue problem on the correlation matrix C = S S ∈ R Ns×Ns . Then where Λ = diag {λ 1 , λ 2 , · · · , λ Ns } ∈ R Ns×Ns is the diagonal matrix of eigenvalues and W ∈ R Ns×Ns is the eigenvector matrix. Our POD basis matrix can then be obtained by In practice a reduced basis ψ ∈ R N h ×Nr is built by choosing the first N r columns of ϑ for the purpose of efficient ROMs, where N r N s . This reduced basis spans a space given by The coefficients of this reduced basis (which capture the underlying temporal effects) may be extracted as The POD approximation of our solution is then obtained viâ whereq i h : T × P → R N h corresponds to the POD approximation toq i h . The optimal nature of reconstruction may be understood by defining the relative projection error which exhibits that with increasing retention of POD bases, increasing reconstruction accuracy may be obtained. We remark that for dimension d > 1, the solution variables may be stacked to obtain this set of bases that are utilized for the reduction of each PDE within the coupled system. Another approach may be to obtain reduced bases for each dependent variable within the coupled system and evolve each PDE on a different manifold. Each dependent variable is projected onto bases constructed from its snapshots alone. This affects the computation of the nonlinear term for computing the updates for each dimension in q. In practice, this operation manifests itself in the concatenation of reduced bases to obtain one linear operation for reconstruction of all field quantities.

Convolutional autoencoders
Autoencoders are neural networks that learn a new representation of the input data, usually with lower dimensionality. The initial layers, called the encoder, map the input x ∈ R m to a new representation z ∈ R k . The remaining layers, called the decoder, map z back to R m with the goal of reconstructing x. The objective is to minimize the reconstruction error. Autoencoders are unsupervised; the data x is given, but the representation z must be learned.
More specifically, we use convolutional autoencoders (CAEs) that have convolutional layers. In a convolutional layer, instead of learning a matrix that connects all m neurons of layer's input to all n neurons of the layer's output, we learn a set of filters. Each filter f i is convolved with patches of the layer's input. Suppose a one-dimensional convolutional layer has filters of length m fi . Then each of the layer's output neurons corresponding to filter f i is connected to a patch of m fi of the layer's input neurons. In particular, a one-dimensional convolution of filter f and patch p is defined as f * p = j f j p j (for neural networks, convolutions are usually technically implemented as cross-correlations). Then, for a typical one-dimensional convolutional layer, the layer's output neuron y ij = ϕ(f i * p j + B i ), where ϕ is an activation function and B i are the entries of a bias term. As j increases, patches are shifted by stride s. For example, a one-dimensional convolutional layer with a filter f 0 of length m f0 = 3 and stride s = 1 could be defined so that y 0j involves the convolution of f 0 and inputs j − 1, j, and j + 1. To calculate the convolution, it is common to add zeros around the inputs to a layer, which is called zero padding. In the decoder, we use deconvolutional layers to return to the original dimension. These layers upsample with nearest-neighbor interpolation.
Two-dimensional convolutions are defined similarly, but each filter and each patch are two-dimensional. A twodimensional convolution sums over both dimensions, and patches are shifted both ways. For a typical two-dimensional convolutional layer, the output neuron y ijk = ϕ(f i * p jk +B i ). Input data can also have a "channel" dimension, such as red, green and blue for images. The convolutional operator sums over channel dimensions, but each patch contains all of the channels. The filters remain the same size as patches, so they can have different weights for different channels.
In our study, we use solely one channel for the spatial magnitude of ρη.
It is common to follow a convolutional layer with a pooling layer, which outputs a sub-sampled version of the input.
In this paper, we specifically use max-pooling layers. Each output of a max-pooling layer is connected to a patch of the input, and it returns the maximum value in the patch.

Variational autoencoders
As opposed to CAEs, the variational autoencoder (VAE) [18] takes a Bayesian approach to latent space modeling. We utilize a convolutional VAE architecture, where outer convolutional, pooling and upscaling layers are identical to CAE. The difference only arises in the bottleneck, where the encoder final layers are fully connected layers that project the input x ∈ R n onto the latent-variable z space. Thus, the encoder effectively is a function q(z|x). The fully connected part of the decoder network p(x |z) generates newly sampled x from the latent space. The output of this undergoes upsampling and convolutions (similar to that of the CAE) to reconstruct the image. The VAE also constrains the latentspace variables to follow a normal distribution z ∼ N (µ, σ 2 ) via the Kullback-Leibler divergence (KL divergence) term D KL (q(z|x i )||p(z|x i )). The architecture of VAE (shown in Figure 3) is similar to that of the conventional CAE shown in Figure 2, except at the bottleneck, where the encoder network outputs the mean µ and variance σ 2 .
The inference itself is undertaken using variational inference (VI) by modeling the true distribution q(z|x) using a simple Gaussian distribution, and then minimizing the difference via the KL divergence as an addition to the loss The KL-divergence loss is applied such that the distribution on z is as close to the normal distribution as possible.
With the inclusion of additional complexity to the loss function and the bottleneck, VAEs have more parameters to tune than CAEs. However, this also gives a significant control over the latent-space distribution. Depending on the data, dimensionality reduction using CAEs may not allow for a straightforward interpolation due to the presence of discontinuous clusters in the latent-space representation. On the other hand, VAEs constrain the latent-space representation to follow a pre-specified distribution, hence by design, facilitate easier interpolation.

Gaussian process regression
Dimensionality reduction performed using POD, CAE or VAE techniques results in a representation of the original data, along with a model (the POD bases or decoders) to reconstruct the data for any point in the latent space. We then only require a temporal interpolation scheme fitted on the representation space, so that a time evolution of the dynamical system can be reconstructed, as schematically shown in Fig 1. In our approach, we deploy the use of GPs [48] as our interpolation algorithm. While GPs perform Bayesian regression tasks with a high level of interpretabilty, they are generally computationally expensive for large data sizes. We achieve a considerable reduction of computational cost by fitting in the space of reduced dimensions. In addition, the GPR may also be restricted to a less noisy space compared with the original dataset of the dynamical evolution.
A GP is an accumulation of random variables, every finite collection of which follows a multivariate Gaussian distribution. It can be perceived as a generalization of a multivariate Gaussian distribution with infinite space. GPs are a popular choice for non-linear regression due to their flexibility and ease of use. In addition, one of their main advantages, is that they incorporate a principled way of measuring the uncertainty information, since they provide predictions in distributional form. For the purpose of this paper, a GPR model is used to fit the reduced space from the data-compression algorithms of Section 3. Subsequently, the mean prediction, which corresponds to the maximum a posteriori (MAP) estimate, is used for the reconstructions. We use the GPflow library for the experiments [10].
A GP can be completely specified by its second-order statistics. Thus, a mean function m(x) equal to zero can be assumed, and a positive definite covariance function (kernel) k(x, x ), which can be perceived as a measure of similarity between x and x , is the only requirement to specify the GP.
For the experiments in Section 5 we initially experimented with changepoint kernels [25]. The intuition came from the data (see e.g. Figure 11), where two typical behaviours are commonly observed; a steep increase or decrease of the latent dimension values for early times, and a subsequent, smoother change in direction, that eventually leads to stabilization. At first we examined a changepoint kernel only for the time feature, which was then added to a regular kernel that accounted for the other variables. The results were discouraging, which we attribute to the fact that this kernel structure leads to loss of correlational information between time and the other variables. Subsequently, we examined changepoint kernels that accounted for all parameters in both their sub-kernels. Even though this type of kernel was successful in producing acceptable results, and also detecting adequately the position of the changepoint, the output was only a slight improvement when compared to a standard kernel. Furthermore, the computational time was substantially larger, which led us to abandon the pursuit of a changepoint kernel. Instead, we settled on a Matérn 3/2 kernel, due to its versatility, flexibility and smoothness, and more specifically its automatic relevance determination (ARD) extension [3], which incorporates a separate parameter for each input variable and gave a significant improvement in our results.
For a GPR model, we considered a GP f and noisy training observations y of n datapoints x, derived from the true values f (x) with additive i.i.d. Gaussian noise with variance σ 2 n . In mathematical form, that is: where k(·, ·) is the kernel. We obtain the complete specification of the GP, by maximizing the marginal likelihood, which we can acquire by integrating the product of the Gaussian likelihood and the GP prior over f : For testing input x and testing output f , we derive the joint marginal likelihood: where I is the identity matrix.
Finally, by conditioning the joint distribution on the training data and the testing inputs, we derive the predictive distribution During the reconstruction phase, we focus on the predictions that correspond tof .

Experiments
In this section, we outline several experiments designed to assess the various compression frameworks and how they interface with latent-space emulation using our aforementioned GPs. A first series of assessments is solely targeted at assessing the fidelity of reconstruction (i.e., which framework offers the most efficient compression). Following this, we interface latent-space representations of our compressed fields with GP emulators to obtain low-dimensional surrogates with embedded uncertainty quantification. Finally, we outline the ability for the GP emulators to predict the dynamics' evolution at finer temporal resolutions.
For the purpose of training the compression frameworks and the GP emulators, we generate 100 forward-model solves which are obtained by a Latin hypercube sampling of different initial conditions w between −0.5 and 0.5 for each dimensionx andȳ. For each of these 100 simulations, 10 evenly-spaced snapshots in time are obtained to construct our total data set (i.e., we have 1000 flow-fields for ρη of resolution 64×64). We remind the reader that the equationbased simulations require the solution of a system of PDEs (for ρη, ρηu and ρηv), but our emulators will be built from ρη information alone. We split the 100 simulations into 80 for training, 10 for validation and 10 for testing. The validation data set is primarily utilized for early stopping criteria in the deep neural network autoencoders. Note that the training and validation data are combined into one data set for training GP emulators. All statistical and qualitative assessments in the following will be shown for the testing data alone.

Reconstruction
We begin by assessing the ability of our different compression frameworks, i.e., the POD, CAE and VAE. This comparison is obtained by training multiple encoders with varying degrees of freedom (DOF) in the latent space. The results of these experiments can be seen in Table 1. Here, we have chosen 2, 4, 8, 16, 30 and 40 DOF for all of our compression frameworks and have compared the fidelity of the reconstruction. We use metrics given by the coefficient of determination (R 2 ), mean squared error (MSE), and mean absolute error (MAE) to compare the true and reconstructed fields for testing data sets. Our analysis of the metrics indicate that the CAE is able to reach optimal reconstruction accuracy faster than both the POD and VAE. Both the VAE and CAE are seen to possess an advantage over the POD, due to their ability to find a non-linear low-dimensional manifold. POD, instead, obtains a linear affine subspace of the high-dimensional system. Interestingly, with increasing DOF in latent space (∼40) the POD method is seen to outperform the VAE. We also note that the CAE and VAE frameworks obtain their peak accuracy at around 8 DOF in the latent space and proceed to saturate in accuracy with greater latent space dimensions. We remark that this aligns with the lack of any guarantees on convergence with increasing DOF for these nonlinear compression frameworks. POD, instead, is guaranteed to converge with increasing latent space dimensions.
Following these quantitative assessments, we assess the reconstruction fidelity of the different frameworks by comparing contours from the different methods with varying DOF in latent space. Figure 4 shows the performance for our three compression frameworks for four DOF in the latent space. At this coarse resolution (in latent space), the linear compression of POD is inadequate at capturing the coherent features in the solution field upon reconstruction. This is due to the advective nature of the dynamics of this data set and the associated high Kolmogorov width. In contrast, both CAE and VAE are able to identify coherent structures in the flow field after being reconstructed from a latent space. Note, however, that the CAE is able to identify the crests and troughs in the flow field in a more accurate manner. We observe similar results from an eight-dimensional latent space where the CAE and VAE are still seen to outperform POD (although improvements in the latter may be observed). For both cases (i.e., four and eight DOF), the CAE and VAE are seen to struggle with reconstructing the dissipating coherent features later in the evolution of the system. For completeness, we also show a result for a forty-dimensional latent space in Figure 6, where POD can be seen to capture the spatio-temporal trends of the true solution appropriately. Note that improvements in the CAE and VAE are marginal with the POD outperforming both these frameworks later in the evolution of system (at t = 0.09).

Latent space forecasts
Next we test the ability of our trained GP emulators to forecast the evolution of systems in their latent space representations. For this assessment, we choose trained encoders (with 4,8 and 40 latent space dimensions) to obtain training and validation data for fitting our previously introduced GPs. Once trained, the GPs are tasked with predicting the evolution of the latent state in reduced space for a set of test initial conditions. Figure 7 shows the latent space evolution of a testing simulation over time for the three different compression methodologies. This result utilizes solely four latent dimensions. It is readily apparent that CAE compression leads to a smooth evolution of the system in latent space. In contrast, POD and VAE methods display significant oscillations. The GP is able to capture the behavior of the system evolution and also provides confidence intervals which are based on two standard deviations around the mean. With the exception of a few instances in time, the confidence intervals are able to envelope the true evolution of the system. We see similar results in Figure 8 where eight latent space DOF are obtained for each of our compression frameworks. While the POD is seen to provide oscillatory system evolution (as before), the CAE and VAE are smooth. In either case, the constructed GP is able to recover the evolution well.

Reconstruction from latent space forecasts
We proceed by assessing the fidelity of the reconstructions from latent space using the GP forecasts of the previous subsection. Qualitative comparisons for the reconstruction of a test simulation are shown in Figure 9 for a fourdimensional latent space at three different times. The CAE and VAE compression is seen to outperform the linear reduced-basis constructed by POD. This aligns with past studies where nonlinear compression methods have outperformed the POD. The CAE is seen to be more accurate than the VAE due to its deterministic formulation during compression. We observe the CAE and VAE to outperform the POD for the eight-dimensional latent space as well, as shown in Figure 10. Qualitatively, the CAE is seen to be the best compression technique of all the methods. Table 2 shows different metrics to establish these conclusions quantitatively. Note that all these metrics are evaluated on the reconstructed data in physical space.

Temporal super-resolution
It must be highlighted that ML-based time-series forecasting methods are generally formulated in a discrete fashion where the temporal resolution of the training data determines the resolution of the ROM deployment. Also, most studies of ML-based forecasting in time assume a regular sampling of the state in time. In practice, due to simulation or experimental limitations, state information may be available sparsely in time and at irregular intervals. Therefore,   obtained in the discrete training and testing data. This is made possible due to the continuous function approximation property of the GPR.
To assess this capability, we re-generate our testing data which is now sampled 10 times more finely in the temporal dimension. We utilize our pretrained CAE (solely trained on the coarsely sampled training data), to compress this system evolution to an eight-DOF latent space. Following this, our previously trained GP is tasked with sampling at the intermediate points in time, which correspond to this finely sampled testing data. Note that the GP is trained only with the coarse data set as well. Thus, this assessment represents interpolation in both space and time. The results for the GP forecast in this assessment are shown in Figure 11. A good agreement can be observed between the true latent space trajectory and the GP interpolated counterpart. We also direct our attention to the forecast behavior for latent dimensions 2, 6 and 7 where uncertainties are seen to oscillate corresponding to the coarsely sampled training points. We qualitatively assess the accuracy of the temporal interpolation as shown in Figure 12, where the reconstruction from the true latent space trajectory and the GP interpolation are compared against the truth. A good agreement is seen. Naturally, the Root Mean Square Error (RMSE) and MAE for this finely sampled test data set is seen to be higher that the case for interpolating solely on the coarser sample locations (see Table 2 for comparison), but this adds a useful utility to this ROM strategy.

Conclusions
The development of parameteric non-intrusive reduced-order models for advective PDE systems has great applications for cost reductions in large numerical simulation campaigns across multiple domain sciences. This article addresses Figure 6: Reconstruction from true latent space evolutions for a forty-dimensional latent space for time t = 0.03 (top), t = 0.06 (middle), t = 0.09 (bottom) for a testing initial condition. Qualitatively, for the same data set and similar architectures, the CAE outperforms the VAE and is comparable to the POD encoding. This is because of the greater DOF in the latent space.  their limitations associated with reduced interpretability by proposing the use of GPRs, conditioned on time and system control parameters that provide quantification of uncertainty. This is particularly useful when nonlinear compression techniques, such as autoencoders, are used for efficient DOF reduction. In addition to the ability to interpolate in initial condition space, we also investigate the ability of the proposed framework to interpolate in time. This addresses the fact that the sampling of a dynamical system for training these ROMs may not match the emulation temporal resolution requirement. We also remark that the modular nature of the compression and time evolution allows for the use of conventional reduced-basis methods such as the POD for dynamics which are intrinsically low-dimensional.
Our results indicate that the proposed model-order reduction technique is successful at dealing with advective dynamics through assessments on the inviscid shallow-water equations. We establish this by testing on unseen initial conditions for our system evolution where a low-dimensional evolution successfully replicates high-dimensional results when coupled with a convolutional or variational autoencoder. The non-intrusive nature of our framework also allows for construction of emulators from remote-sensing or experimental data. This is of value when the underlying governing PDEs are not known a priori. Extensions to the present study shall investigate the integration of a feedback loop to sample points in control parameter space with the knowledge of prediction uncertainty. Through this, we aim to establish continually-learning model-order reduction frameworks for advective problems spanning large physical regimes.     Figure 12: True fields (left), CAE reconstructions (middle) and GP-interpolated reconstructions (right) for a testing simulation that has been sampled 10 times more finely. Note that the shown instances in time are not the usual sampling points utilized for our CAE or GP training.