LaSDI: Parametric Latent Space Dynamics Identification

Physical simulations have influenced developments in engineering, technology, and science more rapidly than ever before. The widespread application of simulations as digital twins is one recent example. Advances in optimization and control theory have also increased simulation predictability and practicality. Physical simulations can reproduce responses in a variety of systems, from quantum mechanics to astrophysics, with considerable detail and accuracy, and they provide information that often cannot be obtained from experiments or in-situ measurements due to associated high risks and measurement difficulties. However, high-fidelity forward physical simulations are computationally expensive and, thus, make intractable any decision-making applications, such as inverse problems, design optimization, optimal controls, uncertainty quantification, and parameter studies, where many forward simulations are required.

To compensate for the computational expense issue, researchers develop several surrogate models to accelerate the physical simulations with high accuracy. One particular type is projection-based reduced order model (pROM), in which the full state fields are approximated by applying linear or nonlinear compression techniques. A popular linear compression technique include proper orthogonal decomposition (POD) [1], reduced basis method [2], and balanced truncation method [3], while a popular nonlinear compression technique is auto-encoder (AE) [4, 5, 6, 7]. The linear compression-based pROM is successfully applied to many different problems, such as Lagrangian hydrodynamics [8, 9], Burgers equations [10, 11, 12], nonlinear heat conduction problem [13], aero-elastic wing design problem [14], Navier–Stokes equation [15], computational fluid dynamics simulation for aircraft [16], convection–diffusion equation [17], Boltzman transport problem [18], topology optimization [19], the shape optimization of a simplified nozzle inlet model and the design optimization of a chemical reaction [20], and lattice-type structure design [21]. These projection-based ROMs are further categorized by its intrusiveness. The intrusive ROMs plug the reduced solution representation into the underlying physics law, governing equations, and its numerical discretization methods, such as the finite element, finite volume, and finite difference methods. Therefore, they are physics-constrained data-driven approaches and require less data than the methods that only require data to achieve the same level of accuracy. However, one needs to understand the underlying numerical methods of solving the high-fidelity simulation to implement the intrusive ROMs. Furthermore, the intrusive ROMs are only applicable when the source code for a high-fidelity physics solver is available, which is not the case for certain applications. On the other hand, the non-intrusive ROMs do not require the access to the source code of a high-fidelity physics solver. Therefore, we will focus on developing the non-intrusive ROMs, which use only data to approximate the full state fields.

Among many non-intrusive methods, various interpolation techniques are used to build a nonlinear map that predicts new outputs for new inputs. The interpolation techniques include, but not limited to, Gaussian processes [22, 23], radial basis functions [24, 25], Kriging [26, 27], and convolutional neural networks [28, 29]. Among them, neural networks have been the most popular framework because of their rich representation capability supported by the universal approximation theorem. Such surrogate models have been applied to various physical simulations, including, but not limited to, particle simulation [30], nanophotonic particle design [31], porous media flow [32, 33, 34, 35], storm prediction [36], fluid dynamics [37], hydrology [38, 39], bioinformatics [40], high-energy physics [41], turbulence modeling [42, 43, 44, 45], uncertainty propagation in a stochastic elliptic partial differential equation [46], bioreactor with unknown reaction rate [47], barotropic climate models [48], and deep Koopman dynamical models [49]. These methods have the potential to be extended to experimental video data as done in [50], where intrinsic dimension is sought through auto-encoder training. However, these methods lack the interpretability due to the black-box nature caused by its complex underlying structure of interpolators, e.g., neural networks. Furthermore, they do not aim to predict the state field solution in general. Therefore, if the target output quantity of interest changes, the model needs to go through the expensive training phase again, which shows the lack of generalizability.

To improve the interpretability and generalizability, we focus on developing a latent space dynamics learning algorithm, where the whole state field data is compressed into a reduced space and the dynamics in the reduced space is learned. Two different types of compression are possible, i.e., linear and nonlinear compression. The popular linear compression can be realized by the singular value decomposition (SVD), justified by the proper orthogonal decomposition (POD) framework. The popular nonlinear compression can be accomplished by the auto-encoder (AE), where the encoder and decoder are designed with neural networks. After the compression, the data size is reduced tremendously. Moreover, the dynamics of the data within the reduced space is often much simpler than the dynamics of the full space. For example, Figure 1 shows the reduced space dynamics corresponding to 2D radial advection problems. This motivates our current work to identify the simplified and reduced but latent dynamics with a system identification regression technique.

There are many existing latent space learning algorithms. DeepFluids [51] uses the auto-encoder for nonlinear compression and applies the latent space time integrator to propagate the solution in the reduced space. The authors in [34] and [52] use both linear and nonlinear compressions and apply some interpolation techniques, such as artificial neural networks (ANNs) and radial basis function (RBF) interpolations within the reduced space to predict the solution for new parameter value. Xie, et al., in [53] uses POD as linear compression technique and apply linear multi-step neural network to predict and propagate the latent space dynamical solutions. Hoang, et al., in [54] compresses both space and time solution space with POD, which was first introduced in [11]. Then they use several interpolation techniques, such as polynomial regression, $k$-nearest neighbor regression, random forest regression, and neural network regression, within the reduced space formed by the space–time POD. However, all the aforementioned methods use a complex form as a latent space dynamics model, whose structure is not well understood and lacking the interpretability. On the other hand, Champion, et al., in [55] use an auto-encoder for the compression and the sparse identification of nonlinear dynamics (SINDy) to identify ordinary differential equations (ODEs) for the latent space dynamics, which is simpler than neural networks, improving the interpretability. However, they fail to generalize the method to achieve a robust parametric model. It is partly because one system of ODEs is not enough to cover all the variations within the latent space due to parameter change.

Therefore, we propose to identify latent space dynamics with a set of local system of ODEs that is tailored to improve the accuracy on a local area of the parameter space. This implies that each local model has a trust region, which covers a sub-region of the whole parameter space. If the set of the trust region covers the whole parameter space, then our model is guaranteed to achieve a certain accuracy level everywhere, achieving truly parametric model based on latent space learning algorithm. At the same time, the framework enables a faster calculation than the corresponding high-fidelity model. We call this framework, LaSDI, which stands for Latent Space Dynamics Identification.

The procedure of LaSDI is summarized by four distinct steps below and its schematics are depicted in Figure 2:

1. 1.

   Data Generation: Generate full order model (FOM) simulation data, i.e., parametric time dependent solution data

2. 2.

   Compression: Apply compression on the simulation data, through either singular value decomposition or auto-encoder to form latent space data.

3. 3.

   Dynamics Identification: Identify the governing equations that best matches the latent space data in a least-squares sense.

4. 4.

   Prediction: Use the identified governing equation to predict the latent space solution for a new parameter point. In this paper, we assume that the parameter affects only the initial condition. The predicted latent space solution is reconstructed to the full state by decompression.

In fact, the four-step procedure above is prevalent in developing data-driven reduced-order models to accelerate physical simulations. The traditional data-driven projection-based reduced order models mentioned earlier follow the four-step procedure if the dynamics identification model is derived from the existing physical governing equations. The dynamic mode decomposition (DMD) approaches [56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69] also follow the four-step procedure, where linear compression via SVD is used for Step 2 and linear ODEs are generally used for Step 3. Additionally, the operator inference (OpInf) approaches [70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82], which utilize linear compression and fit the latent space dynamics data into polynomial models (mostly quadratic models), also follow the four-step procedure above. As a matter of fact, the proposed LaSDI in this paper is similar to both DMD and OpInf approaches. The major two differences include (i) in Step 2, LaSDI allows for nonlinear compression via an auto-encoder, while DMD and OpInf use a linear compression via the SVD; and (ii) in Step 3, DMD and OpInf only considers polynomial terms in the candidate library, while LaSDI allows general nonlinear terms. Therefore, LaSDI can be viewed as a generalization of both DMD and OpInf.

The main contributions of this paper include

• •

  A generalized data-driven physical simulation framework, i.e., LaSDI, is introduced.

• •

  A novel local latent space learning algorithm that demonstrates robust parametric prediction, is introduced.

• •

  Three novel dynamics identification algorithms are introduced, i.e., global, local, and interpolated ones.

• •

  Both good accuracy and speed-up are demonstrated with several numerical experiments.

The paper is organized in the following way: Section 3 technically describes LaSDI, where each four distinct steps are explained in details. Section 3.1 describes how to generate and store the simulation data. Section 3.2 discusses two different types of compression techniques. Section 3.3 elaborates the procedure of the dynamics identification and subsequently describes three different types. Section 3.4 elucidates the prediction step. We demonstrate the performance of LaSDIs for four different numerical experiments in Section 4, where global, local, and interpolated dynamics identification algorithms are compared and analyzed. Finally, we summarize and discuss the implications, limitations, and potential future directions for LaSDI in Section 5.

We formally state a parameterized dynamical system, characterized by the following time dependent ordinary differential equations (ODEs):

where $t\in[0,t_{f}]$ denotes time with the final time $t_{f}\in\mathbb{R}_{+}$, and $\boldsymbol{u}(t;\boldsymbol{\mu})$ denotes the time-dependent, parameterized state implicitly defined as the solution to System (1) with $\boldsymbol{u}:[0,t_{f}]\times\mathcal{D}\rightarrow\mathbb{R}^{N_{s}}$. Further, $\boldsymbol{f}:\mathbb{R}^{N_{s}}\times[0,t_{f}]\rightarrow\mathbb{R}^{N_{s}}$ with $(\boldsymbol{w},\tau)\rightarrow\boldsymbol{f}(\boldsymbol{w},\tau)$ denotes the scaled velocity of $\boldsymbol{u}$, which can be either linear or nonlinear in its first argument. The initial state is denoted by $\boldsymbol{u}_{0}:\mathcal{D}\rightarrow\mathbb{R}^{N_{s}}$, and $\boldsymbol{\mu}\in\mathcal{D}$ denotes the parameters with parameter domain $\mathcal{D}\subseteq\mathbb{R}^{n_{\mu}}$. We assume that the parameter affects only the initial condition. System (1) can be considered as semi-discretized version of a system of partial differential equations (PDEs), whose spatial domain is denoted as $\Omega\in\mathbb{R}^{d}$, $d\in\mathbb{N}(3)$, where $\mathbb{N}(N):=\{1,\ldots,N\}$. The spatial discretization can be done through many different numerical methods, such as the finite difference, finite element, and finite volume methods.

Many different time integrators are available to approximate the time derivative term, $d\boldsymbol{u}/dt$, e.g., explicit and implicit time integrators. A uniform time discretization is assumed throughout the paper, characterized by time step $\Delta t\in\mathbb{R}_{+}$ and time instances $t_{n}=t_{n-1}+\Delta t$ for $n\in\mathbb{N}({N_{t}})$ with $t_{0}=0$, ${N_{t}}\in\mathbb{N}$. To avoid notational clutter, we introduce the following time discretization-related notations: $\boldsymbol{u}_{n}:=\boldsymbol{u}(t^{n};\boldsymbol{\mu})$, $\tilde{\boldsymbol{u}}_{n}:=\tilde{\boldsymbol{u}}(t^{n};\boldsymbol{\mu})$, $\hat{\boldsymbol{u}}_{n}:=\hat{\boldsymbol{u}}(t^{n};\boldsymbol{\mu})$, and $\boldsymbol{f}_{n}:=\boldsymbol{f}(\boldsymbol{u}(t^{n};\boldsymbol{\mu}),t^{n};\boldsymbol{\mu})$.

Although many advanced numerical methods are available to solve System (1), as the problem size increases, i.e., a large $N_{s}$, and the computational domain $\Omega$ gets geometrically complicated, the overall solution time becomes impractically slow, e.g., taking a day or week for one forward simulation. The proposed method, i.e., LaSDI, accelerates the computationally expensive simulations.

This section mathematically describes LaSDI. As shown in Introduction 1, LaSDI consists of four steps, i.e., data generation, compression, identification, and prediction. Each step is described in the subsequent sections.

We demonstrate the performance of LaSDI for four different numerical examples, i.e., 1D and 2D Burgers equations, heat conduction, and radial advection problems. The governing equations, initial conditions, boundary conditions, parameter space, and domain are specified in Table 1.

For both 1D and 2D Burgers problems, a uniform space discretization (i.e., $dx=6/1000$ for 1D and $dx=dy=1/10$ for 2D Burgers problems) is used for the full order model solve. The first order spatial derivative term is approximated by the backward difference scheme, while the diffusion terms are approximated by the central difference scheme. This generates a semi-discretized PDE characterized by a system of nonlinear ordinary differential equations specified in Eq. (1). We integrate each of these systems over a uniform time-step (i.e., $\Delta t=1/1000$ for 1D and $\Delta t=2/1500$ for 2D Burgers problems) using the implicit backward Euler time integrator, i.e., $\boldsymbol{u}_{n}-\boldsymbol{u}_{n-1}=\Delta t\boldsymbol{f}_{n}$.

For the full order model of the nonlinear heat conduction problem, we use the second order finite element to discretize the spatial domain. The spatial discretization starts with $64\times 64$ uniform squares, which are divided into triangular elements, resulting in 8192 triangular elements in total. A uniform time step of $\Delta t=.01$ is used with the SDIRK3 implicit L-stable time integrator (i.e., see Eq. (229) of [87]). For the temperature field dependent conductivity coefficient, we linearize the problem by using the temperature field $\boldsymbol{u}$ from the previous time step.¹¹1The source code of the full order model for the nonlinear heat conduction problem can be found at https://github.com/mfem/mfem/blob/master/examples/ex16.cpp

For the full order model of the radial advection problem, we use the third order discontinuous finite element to discretize the spatial domain. The spatial discretization is dictated by a square-mesh with periodic boundary conditions. We use $24\times 24$ grid of square finite elements. The finite element data is interpolated to generate a $64\times 64$ uniform grid across the spatial domain. We use a uniform time step of $\Delta t=.0025$ with the RK4 explicit time integrator.²²2The source code of the full order model for the radial advection problem can be found at https://github.com/mfem/mfem/blob/master/examples/ex9.cpp

These full order models are used to generate training data. Table 2 shows testing parameter set, $\mathcal{T}$, as well as all the training samples, $\mathcal{S}$, for four problems. Note that the testing points include the training points. The accuracy is measured by the maximum relative error, $e:\mathcal{D}\rightarrow\mathbb{R}_{+}$, at each testing parameter point, which is defined as

where $\boldsymbol{\mu}^{*}$ is a testing parameter point.

The computational cost is measured in terms of the wallclock time. Specifically, timing is obtained by performing calculations on a IBM Power9 @ 3.50 GHz and DDR4 Memory @ 1866 MT/s. The auto-encoders are trained on a NVIDIA V100 (Volta) GPU with 3168 NVIDIA CUDA Cores and 64 GB GDDR5 GPU Memory using PyTorch [55] which is the open source machine learning frame work.

We have introduced a framework for latent space dynamics identification (LaSDI) that builds efficient parametric models by exploiting both local and global regression techniques. Two types of local dynamics identification models are introduced: (i) region-based and (ii) point-wise interpolated DIs. Latent spaces are generated in two different ways as well, i.e., using linear and nonlinear compression techniques. The compression techniques, such as proper orthogonal decomposition and auto-encoder, reduce large-scale simulation data to small-scale data whose size is manageable for regression techniques. The framework is general enough to be applicable to any physical simulations. We have applied various LaSDI models to four different problems, i.e., 1D and 2D Burgers, heat conduction, and radial advection problems.

We observe that local versus global DI models need to be determined case by case. For example, a local DI model outperforms global DI for LaSDI-NM for 1D Burgers problem with nine training points (see the case of Degree 0 dynamical system in Figure 8) and 2D Burgers problem with 25 training points (see Figure 13) and all the cases for radial advection problems (see Figure 18), while a global DI or local DI with larger number of points outperform a local DI with smaller number of points for LaSDI-NM for 1D Burgers problem with 25 training points (see Figure 10).

The nonlinear compression techniques, such as auto-encoder, outperforms the linear compression technique for problems with slowly decaying Kolmogorov’s width, while the linear compression techniques, such as proper orthogonal decomposition, outperforms the nonlinear compression technique for problems with fast decaying Kolmogorov’s width. The nonlinear compression was able to bring the projection error down with small latent space dimension even for the problems with slowly decaying Komogorov’s width. However, the computational cost of the nonlinear manifold compression is much higher than the linear compression, so we recommend using linear compression if the problem shows fast singular value decay.

There are several future directions to consider in the context of latent space dynamics identification framework. First, we have only considered predefined uniformly distributed training parameters, which is fine with a small size parameter space dimension of one or two. However, the number of uniformly distributed training points increases exponentially as the dimension of the parameter space increases, which makes the generation of the simulation data extremely challenging. Therefore, adaptive and sparse parameter sampling for LaSDI needs to be developed, which will lead to an optimal sampling. Second, we have only considered the parameters that affects an initial condition. However, a parameter of interest may not alter the initial condition, e.g., material properties. The LaSDI framework introduced in this paper cannot handle this case. Therefore, parameterization other than the initial condition needs to be developed. Third, we have extensively used a regression technique to identify dynamics in a latent space. However, this is not the only option. For example, any system identification techniques or operator learning algorithms are applicable. Fourth, the architecture of the neural network for LaSDI-NM may not have been optimally chosen. The more thorough neural architecture search would lead to the better performance of LaSDI-NM models. Finally, because LaSDI is a non-intrusive method, it can be applied to data sets that might not be generated by simulations. These methods, along with adding regularization to the training process can help discover more robust latent-space dynamical systems for potentially noisy and imperfect data as done in [50].

With these future directions, the potential for LaSDI is to offer fast and accurate data-driven simulation capability. We expect to apply LaSDI to various important applications, such as climate science, manufacturing, and fusion energy.

This work was performed at Lawrence Livermore National Laboratory and partially funded by two LDRDs (21-FS-042 and 21-SI-006). Lawrence Livermore National Laboratory is operated by Lawrence Livermore National Security, LLC, for the U.S. Department of Energy, National Nuclear Security Administration under Contract DE-AC52-07NA27344 and LLNL-JRNL-831849.

This research was supported in part by an appointment with the National Science Foundation (NSF) Mathematical Sciences Graduate Internship (MSGI) Program sponsored by the NSF Division of Mathematical Sciences. This program is administered by the Oak Ridge Institute for Science and Education (ORISE) through an interagency agreement between the U.S. Department of Energy (DOE) and NSF. ORISE is managed for DOE by ORAU. All opinions expressed in this paper are the author’s and do not necessarily reflect the policies and views of NSF, ORAU/ORISE, or DOE.