A Comprehensive Review of Latent Space Dynamics Identification Algorithms for Intrusive and Non-Intrusive Reduced-Order-Modeling

In this chapter, we consider physical phenomenons described by governing PDEs with the following form:

In Equation (1), the solution $\mathbf{u}$ may represent either a scalar or vector field over the time-space domain $[0,t_{\text{max}}]\times\Omega$. The spatial domain $\Omega$ can be of any dimension, and the differential operator $\mathbf{f}$ might include a mix of linear and/or non-linear combinations of spatial derivatives and source terms. The governing equations, along with their initial and boundary conditions, are influenced by a set of parameters, represented by a vector $\boldsymbol{\mu}$. The parameter space is designated as $\mathcal{D}\subseteq\mathbb{R}^{n_{\mu}}$ and can have any dimension (for instance, a 2D-parameter space with $n_{\mu}=2$ is considered in subsequent sections). For a specific parameter vector $\boldsymbol{\mu}^{(i)}\in\mathcal{D}$, it’s assumed that we can access the corresponding discretized solution of Equation (1), denoted as $\mathbf{U}^{(i)}$. Here, $\mathbf{U}^{(i)}$ is a matrix composed of sequential snapshots $\mathbf{u}_{n}^{(i)}$ at each time step $n$, arranged as $\mathbf{U}^{(i)}=[\mathbf{u}_{0}^{(i)},\dots,\mathbf{u}_{N_{t}}^{(i)}]^{\top}\in\mathbb{R}^{(N_{t}+1)\times N_{u}}$. This solution is acquired either through a full-order model solver or an experiment. In this chapter, all the FOM solutions are compiled into a $3^{\text{rd}}$ order tensor dataset $\mathbf{U}\in\mathbb{R}^{N_{\mu}\times(N_{t}+1)\times N_{u}}$, where $N_{t}$ and $N_{u}$ are the number of time steps and degrees of freedom, respectively, and $N_{\mu}$ is the number of available FOM solutions.

An auto-encoder [34, 28] is a type of neural network specifically tailored for compressing large datasets by diminishing their dimensions through nonlinear transformation. It comprises two interconnected neural networks: the encoder, labeled as $\phi_{e}$ and parameterized by $\boldsymbol{\theta}_{\text{enc}}$, and the decoder, labeled as $\phi_{d}$ and parameterized by $\boldsymbol{\theta}_{\text{dec}}$. The encoder processes an input data snapshot $\mathbf{u}_{n}^{(i)}\in\mathbb{R}^{N_{u}}$ and generates a compressed version $\mathbf{z}_{n}^{(i)}\in\mathbb{R}^{N_{z}}$ in a latent space. Here, $N_{z}$ denotes the latent dimension, i.e., the number of variables in the latent space. This number is chosen based on design preferences, typically ensuring $N_{z}<\!\!<N_{u}$. Similar to the matrix $\mathbf{U}^{(i)}$, we concatenate the latent representations at each time step into a matrix $\mathbf{Z}^{(i)}=[\mathbf{z}_{0}^{(i)},\dots,\mathbf{z}_{N_{t}}^{(i)}]^{\top}\in\mathbb{R}^{(N_{t}+1)\times N_{z}}$. The latent variables for each time step and for each parameter $\boldsymbol{\mu}^{(i)}$ are compiled into a third-order tensor $\mathbf{Z}\in\mathbb{R}^{N_{\mu}\times(N_{t}+1)\times N_{z}}$, similar in structure to the tensor $\mathbf{U}$. The decoder takes each $\mathbf{z}_{n}^{(i)}$ as input and generates a reconstructed version of $\mathbf{u}_{n}^{(i)}$, denoted as $\hat{\mathbf{u}}_{n}^{(i)}$.

$\boldsymbol{\theta}_{\text{enc}}$ and $\boldsymbol{\theta}_{\text{dec}}$ are learned using a numerical optimization algorithm that minimizes the $L_{2}$ norm of the difference between the set of input solutions $\mathbf{U}$ and the set of reconstructed solutions $\hat{\mathbf{U}}$. The reconstruction loss is defined as:

It should be noted that the auto-encoder acts as a non-linear data projection. Alternatively, a linear projection method, such as POD, may be employed. This option alleviates the cost of training an auto-encoder, but it may result in lower accuracy, especially for advection-dominated problems. A more detailed comparison between auto-encoders and POD is outlined in LaSDI [25] and WLaSDI [72].

The encoder’s role is to compress high-dimensional physical data, such as solutions of PDEs that span both space and time, into a more compact array of discrete and abstract latent variables, defined over the time dimension. As a result, the latent space can effectively be viewed as a dynamical system, but now governed by ordinary differential equations (ODEs). This feature is a key element in LaSDI algorithms. For each discrete time step, the behavior of the latent variables within this latent space can be represented by a specific equation form. At each time step, the dynamics of the latent variables in the latent space can be described by an equation of the following form:

where $\psi_{DI}$ represents a Dynamics Identification (DI) function governing the latent space dynamics, which can be defined as a system of ODEs and identified using the Sparse Identification of Nonlinear Dynamics (SINDy) method [11]. SINDy involves creating a library, $\boldsymbol{\Theta}(\mathbf{Z}^{(i)})\in\mathbb{R}^{(N_{t}+1)\times N_{l}}$, which is composed of $N_{l}$ linear and nonlinear potential terms that could be part of the ODEs. This method is based on the premise that the time derivatives of $\mathbf{Z}$, denoted as $\dot{\mathbf{Z}}$, can be represented as a linear combination of these selected terms. Consequently, the equation that describes the system’s behavior on the right-hand side can be estimated in the following manner:

where $\mathbf{\Xi}^{(i)}\in\mathbb{R}^{N_{z}\times N_{l}}$ denotes an ODE coefficient matrix associated with $\boldsymbol{\mu}^{(i)}$. The selection of terms in $\boldsymbol{\Theta}(\cdot)$ is a design choice. To capture the latent space dynamics more accurately, it may be desirable to include a broad variety of terms, although this may lead to sets of ODE more challenging and longer to solve numerically. $\boldsymbol{\Theta}(\cdot)$ can be written explicitly as:

where $\mathbf{b}_{q}$ are basis functions that output arbitrary linear and/or non-linear transformations of each terms in $\mathbf{z}_{n}^{(i)}=[z_{n,1}^{(i)},\dots,z_{n,N_{z}}^{(i)}]$ to build the SINDy library. $z^{(i)}_{n,j}$ represents the $j^{\text{th}}$ latent variable for parameter $\boldsymbol{\mu}^{(i)}$ at time step $n$. Sparse linear regressions are performed between $\boldsymbol{\Theta}(\mathbf{Z}^{(i)})$ and $dz^{(i)}_{:,j}/dt$ for each $j\in[\![1,N_{z}]\!]$. The resulting coefficients associated with each SINDy term are stored in a vector $\boldsymbol{\xi}_{j}^{(i)}=[\xi_{j,1}^{(i)},\dots,\xi_{j,N_{l}}^{(i)}]\in\mathbb{R}^{N_{l}}$. The system of SINDy regressions can be expressed as:

Each set of coefficients $\boldsymbol{\xi}_{j}^{(i)}$ is compiled into a coefficient matrix $\mathbf{\Xi}^{(i)}=[\boldsymbol{\xi}_{1}^{(i)\top},\dots,\boldsymbol{\xi}_{N_{z}}^{(i)\top}]^{\top}\in\mathbb{R}^{N_{z}\times N_{l}}$. Since a unique set of ODEs controls the dynamics within the latent space for each parameter vector $\boldsymbol{\mu}^{(i)}\in\mathcal{D}$, several SINDy regressions are conducted simultaneously to identify the respective sets of ODE coefficients $\mathbf{\Xi}^{(i)}$, where $i\in[\![1,N_{\mu}]\!]$. The ensemble of ODE coefficient matrices $\mathbf{\Xi}=[\mathbf{\Xi}^{(1)},\dots,\mathbf{\Xi}^{(N_{\mu})}]\in\mathbb{R}^{N_{\mu}\times N_{z}\times N_{l}}$, corresponding to each ODE system, is learned by minimizing the mean-squared-error loss specific to SINDy.

The time derivatives $\dot{z}_{n,j}^{(i)}$ may be computed using a finite difference [8] or, alternatively, using the chain rule [25, 31]:

The aforementioned approach is a local DI method where each parameter has its own coefficient matrix. During testing, the local SINDy of training parameters can be exploited to estimate the SINDy associated with the testing parameter. More details will be presented in Section 2.7. Alternatively, a global DI approach can also be adopted in which only one coefficient matrix is employed to model the latent space dynamics for all parameters.

The coefficient matrix $\mathbf{\Xi}^{(i)}=[\boldsymbol{\xi}_{1}^{(i)\top},\dots,\boldsymbol{\xi}_{N_{z}}^{(i)\top}]^{\top}$ referenced in equation (5) can alternatively be determined using weak-form equation learning methods. This approach replaces the need to directly approximate pointwise derivatives from data (in the presence of noise, this is a particularly challenging task). Instead, the variance-reduction nature of the weak-form facilitates a more robust and accurate system recovery. In this section, we present WLaSDI [72], a direct extension of LaSDI [25] that leverages the weak-form equation learning technique.

We begin by transforming equation (7) into the weak-form through multiplication by an absolutely continuous test function $\phi(t):\mathbb{R}\to\mathbb{R}$ and integrating it over the time domain:

Integration by parts yields:

Choosing $\phi$ to be compactly supported in $(t_{a},t_{b})$, i.e., $\phi(t_{b})=\phi(t_{a})=0$, we have:

We use the trapezoidal rule to discretize the integrals in equation (10) numerically and assume that observations of the system’s state have a uniform time step of $\Delta t$.^1††1It is possible to use a non-uniform time interval, but this necessitates modifications to the test function matrices. Let $\{\phi_{k}\}_{m=1}^{N_{k}}$ be the set of compactly supported test functions placed uniformly along the time domain, we define the test function matrices:

From there, we have

for $z_{:,j}^{(i)}=[z_{0,j}^{(i)},\cdots,z_{N_{t},j}^{(i)}]^{\top}$ and $j\in[\![1,N_{z}]\!]$. By solving the least square problem, we identify the coefficient vector $\boldsymbol{\xi}_{j}^{(i)\top}$ subject to:

The weak-SINDy loss can then be defined as:

The test functions employed in weak-form equation learning include high order polynomial or $C^{\infty}$ bump functions. For further details about how to choose these test functions, please refer to [10, 52, 53].

tLaSDI [56] utilizes a neural network-based model to embed the General Equation for Non-Equilibrium Reversible-Irreversible Coupling (GENERIC) formalism [29, 55] into the latent space dynamics. The GENERIC formalism is a mathematical framework that covers both conservative and dissipative systems, offering a general description of beyond-equilibrium thermodynamic systems. This formalism models dynamical systems through four key functions – the scalar functions $E$ and $S$ that representing the system’s total energy and entropy, respectively, and two matrix-valued functions, $L$ and $M$, referred to as the Poisson and friction matrices, respectively.

The degeneracy and structural conditions ensure the first and second laws of thermodynamics, i.e., the properties of energy conservation and non-decreasing entropy, $\frac{d}{dt}E(\mathbf{z})=0$ and $\frac{d}{dt}S(\mathbf{z})\geq 0$.

In LaSDI [25], the auto-encoder and SINDy are trained separately. Hower, for enhanced accuracy and robustness when identifying the latent space dynamics, in gLaSDI [31], the auto-encoder and SINDy are trained simultaneously. gLaSDI also introduced an additional loss term: the mean-squared-error of the velocity, computed through the chain rule:

In GPLaSDI [8], a penalty is also added to the SINDy coefficients. This enforce smaller values of the ODE coefficients, which leads to better conditioned ODEs. The LaSDI training loss, weighted with hyperparameters $\beta_{1}$, $\beta_{2}$, $\beta_{3}$ and $\beta_{4}$ is:

where $\mathcal{L}_{\text{SINDy}}$ may be either the vanilla SINDy loss (equation (8)) or the weak SINDy loss (equation (14)). In tLaSDI, the latent space dynamics is not learned through SINDy, so the training loss has a slightly different (yet equivalent) formalism:

where $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ and $\lambda_{4}$ are hyperparameters (analog to $\beta_{1}$, $\beta_{2}$, $\beta_{3}$ and $\beta_{4}$ in equation (20)). The first and second terms of $\mathcal{L}_{\text{MOD}}$ are analog to $\mathcal{L}_{\text{SINDy}}$ and $\mathcal{L}_{\text{VEL}}$ in the gLaSDI training loss, respectively. Furthermore, the Jacobian loss $\mathcal{L}_{\text{JAC}}$ is defined as

The loss components are based on the error estimates of ROM approximation provided in [56].

The auto-encoder learns a latent space dynamical representation of each training parameter. In the local Dynamics Identification (DI) approach, a set of SINDy’s or weak-SINDy’s learns the ODEs governing the latent space dynamics of training parameters. However, the set of ODE coefficients associated with a new parameter $\boldsymbol{\mu}^{(*)}$ (distinct from the set of training parameters $\boldsymbol{\mu}^{(i)}\in\mathcal{D}$) remains unknown. It can be estimated by finding an interpolant $\psi:\boldsymbol{\mu}^{(*)}\mapsto\mathbf{\Xi}^{(*)}$ of the learned ODE coefficients, given $\boldsymbol{\mu}^{(i)}\in\mathcal{D}$ and $\mathbf{\Xi}^{(i)}$, $i\in[\![1,N_{\mu}]\!]$. Several interpolants have been proposed in the literature, and we briefly cover three of them here, RBF interpolation, as used in LaSDI [25], $k-$NN, as used in gLaSDI [31], and Gaussian process regression (GP), as introduced in GPLaSDI [8].

Once the auto-encoder has been properly trained and the latent space governing ODEs have been learned (section 2.2, 2.3, 2.4 and 2.6), ROM prediction for an arbitrary test parameter $\boldsymbol{\mu}^{(*)}$ is straightforward and can be made through the following steps:

1. 1.

   Estimate $\mathbf{\Xi}^{(*)}$ by interpolating the set of learned ODE coefficients using one of the interpolator introduced in section 2.7. If using GP regression, in-lieu of $\mathbf{\Xi}^{(*)}$, a sample from the predictive distribution may be used: $\mathbf{\Xi}^{(d)}\sim\mathcal{N}(\mathbf{\Xi}^{(*)}|m^{(*)},s^{(*)2})$, $d\in[\![1,N_{s}]\!]$, with $N_{s}$ an arbitrary number of samples.

2. 2.

   Build and solve the set of ODEs associated with $\boldsymbol{\mu}^{(*)}$ (equation (5)) using a standard numerical integrator. The initial condition $\mathbf{z}_{0}^{(*)}$ can be computed using the encoder ($\mathbf{z}_{0}^{(*)}=\phi_{e}(\mathbf{u}_{0}^{(*)}\,|\,\boldsymbol{\theta}_{\text{enc}})$), where $\mathbf{u}_{0}^{(*)}$ is a function of $\boldsymbol{\mu}^{(*)}$. The approximate latent-space solution is noted $\tilde{\mathbf{Z}}^{(*)}$.

3. 3.

   Make a forward pass through the decoder to map the ROM prediction back into the full-order space. $\tilde{\mathbf{U}}^{(*)}=\phi_{d}(\tilde{\mathbf{Z}}^{(*)}\,|\,\boldsymbol{\theta}_{\text{dec}})$

Note that if using GP regression to interpolate the latent space governing equations, step 1, 2 and 3 may be repeated for a finite number of samples. The ROM prediction can then be taken as the average prediction over all the samples ($\tilde{\mathbf{U}}^{(*)}\equiv\mathbb{E}[\tilde{\mathbf{U}}^{(*)}]$). Similarly, the standard deviation $\mathbb{V}[\tilde{\mathbf{U}}^{(*)}]^{1/2}$ can also be computed to estimate prediction uncertainty. In tLaSDI, the prediction do not involve interpolation methods as GFINN is a global DI model independent of parameters. The tLaSDI prediction for a test parameter $\boldsymbol{\mu}^{(*)}$ can be conducted through the following analog steps:

1. 1.

   Evaluate the network parameters of the encoder and decoder using corresponding hypernetworks, i.e., $\boldsymbol{\theta}_{\text{enc}}^{(*)}:=\phi_{e}^{\text{hyp}}(\boldsymbol{\mu}^{(*)}|\boldsymbol{\theta}_{\text{enc}}^{\text{hyp}})$ and $\boldsymbol{\theta}_{\text{dec}}^{(*)}:=\phi_{d}^{\text{hyp}}(\boldsymbol{\mu}^{(*)}|\boldsymbol{\theta}_{\text{dec}}^{\text{hyp}})$.

2. 2.

   Solve the latent space dynamics (equation (15)) represented by GFINNs using a numerical integrator. The initial condition is given by $\mathbf{z}_{0}^{(*)}=\phi_{e}(\mathbf{u}_{0}^{(*)}\,|\,\boldsymbol{\theta}_{\text{enc}}^{(*)})$. The approximate solution is denoted as $\tilde{\mathbf{Z}}^{(*)}$.

3. 3.

   Use the decoder to return the prediction in the full-order space, i.e., $\tilde{\mathbf{U}}^{(*)}=\phi_{d}(\tilde{\mathbf{Z}}^{(*)}\,|\,\boldsymbol{\theta}_{\text{dec}}^{(*)})$.

For a given test parameter, the ROM prediction error may be estimated using the time-averaged PDE residual of the FOM governing equation:

$\mathbf{r}$ refers to the discretized residual of the governing PDE (equation (1)). $N_{ts}$ is the number of time steps used to estimate the residual. When $\mathbf{r}$ is expensive to evaluate, it may be desirable to evaluate the residual only at a handfull of time steps, and we may thus consider $N_{ts}\ll N_{t}$.

In gLaSDI [31], the training process is initialized using a limited number ($N_{\mu}$) of training FOM simulations. The training is then paused every $N_{up}$ epochs, and the (current) LaSDI model is used to make predictions at a finite number of test parameters $\boldsymbol{\mu}^{(*)}\in\mathcal{D}^{h}$ ($\mathcal{D}^{h}$ is a discretized representation of the parameter space). Each prediction $\tilde{\mathbf{U}}^{(*)}$ is plugged into equation (29) to compute the residual error, and the next sampling parameter can be selected as the one associated with the ROM solution yielding the largest residual:

A high-fidelity simulation can then be run for parameter $\boldsymbol{\mu}^{(N_{\mu}+1)}$, and the resulting solution is added to the training dataset before resuming the training process. This operation may be repeated until the budget for running FOM simulations has been exhausted, or the ROM predictions have become sufficiently accurate. The residual-based greedy sampling is also employed in tLaSDI [56] following the algorithm proposed in [31].

As described in section 2.8, when using GP interpolation, the uncertainty in the latent space dynamics may be propagated through the decoder, and the ROM prediction variance can be computed for any time step $n$:

$\tilde{\mathbf{z}}_{n}^{(*,d)}$ represent the latent space dynamics solved for a sampled ODE coefficient matrix $\mathbf{\Xi}^{(d)}$ (equation (28)). The variance at each time step can be computed and the variance across time and space is written $\mathbb{V}[\tilde{\mathbf{U}}^{(*)}]=[\mathbb{V}[\tilde{\mathbf{u}}_{0}^{(*)}],\dots,\mathbb{V}[\tilde{\mathbf{u}}_{N_{t}}^{(*)}]]$.

An uncertainty based strategy for picking the next sampling parameter $\boldsymbol{\mu}^{(N_{\mu}+1)}$ can now easily be implemented. The general approach is very similar to the residual-based greedy sampling strategy presented in section 3.1, but instead of picking the parameter yielding the largest residual error, we pick the parameter yielding the largest ROM variance (at any point in time and space):

Note that in order to properly quantify the ROM variance ($\mathbb{V}[\tilde{\mathbf{U}}^{(*)}]$), for each $\boldsymbol{\mu}^{(*)}\in\mathcal{D}^{h}$, the GP predictive distribution needs to be sampled $N_{s}$ times, and the corresponding ODE needs to be solved and fed to the decoder. Therefore, this approach requires more ROM prediction than the residual-based greedy sampling approach. However, this approach does not require to evaluate the PDE residual at any time. Thus, GPLaSDI is non-intrusive and may be particularly suitable when the residual is difficult to implement, expensive to evaluate, or simply unknown (e.g. legacy codes for high-fidelity simulations).

In this first application, we use a simple benchmark problem to compare the different LaSDI algorithms introduced in the previous two sections. We consider the following inviscid 1D Burgers equation:

The initial condition is parameterized by $\boldsymbol{\mu}=\{a,w\}\in\mathcal{D}$, and the parameter space is defined as $\mathcal{D}=[0.7,0.9]\times[0.9,1.1]$:

The parameter space is discretized into a square grid $\mathcal{D}^{h}$ with a stepping of $\Delta a=\Delta w=0.01$, resulting in a total of 441 grid points ($21$ values in each dimension). The high-fidelity solver employs a classic backward Euler time integration, and a finite difference discretization in space ($\Delta x=6\cdot 10^{-3}$ and $\Delta t=10^{-3}$).

In this second example, we consider the non-linear heat equation in 2D introduced in gLaSDI [31]:

The diffusion coefficients are chosen as $\kappa=0.5$ and $\alpha=0.01$. The initial condition is parameterized by $\boldsymbol{\mu}=\{a,w\}\in\mathcal{D}$, and the parameter space is defined as $\mathcal{D}=[1,1.4]\times[4,4.3]$: