Sampling low-dimensional Markovian dynamics for pre-asymptotically recovering reduced models from data with operator inference

Introduction

Reduced models have become a ubiquitous tool to make tractable computations that require large numbers of model evaluations in, e.g., uncertainty quantification, optimization, and inverse problems. Traditional model reduction derives reduced models from high-dimensional (full) models of systems that typically are given in the form of partial differential equations (PDEs) and their corresponding discretized operators. The properties of reduced models have been extensively studied by the model reduction community and even rigorous error estimation has been established for certain classes of problems . The aim of data-driven model reduction methods is to learn reduced models from data alone and so to extend the scope of model reduction to settings where the governing equations and the corresponding discrete operators of the high-dimensional systems are unavailable; however, the models learned from data alone typically are only approximations of the reduced models obtained with traditional model reduction and thus establishing the same rigor for the learned models as for reduced models is challenging. In contrast, this work presents an approach to learn low-dimensional models from data that exactly match the reduced models that are obtained with traditional model reduction as if the governing equations and discrete operators of the high-dimensional systems were available. This guarantee of exactly recovering reduced models from data holds pre-asymptotically in the number of data points and for a wide class of high-dimensional systems with polynomial nonlinear terms under certain conditions. Thus, models learned with the proposed approach are the reduced models of traditional model reduction and therefore directly inherit their well-studied properties.

There is a large body of literature on learning dynamical-system models from data. We review only the works that are most relevant for the proposed approach. First, there is system identification that originated in the systems and control community . The Loewner approach was introduced by Antoulas and collaborators and has been extended from linear time-invariant systems to parametrized , bilinear and quadratic-bilinear systems . Under certain conditions, the models learned with the Loewner approach are the reduced models that are obtained with interpolatory model reduction; however, Loewner models are learned from frequency-response data rather than from time-domain data. The work builds on Loewner to learn reduced models of linear time-invariant systems from time-domain data; however, learning from time-domain data can introduce errors and so the learned models can differ from the corresponding Loewner models derived from frequency-response data. Second, there is dynamic mode decomposition that best-fits linear operators to state trajectories with respect to the $L_{2}$ norm. Methods based on the Koopman operator have been developed as one path to extending dynamic mode decomposition to nonlinear dynamical systems . Third, there are methods that learn parsimonious models by exploiting sparsity in the high-dimensional systems, e.g., the work by Schaeffer and collaborators and the work by Kutz, Brunton, and collaborators . The learned models typically are either continuous in the sense that terms of PDEs are learned from a dictionary or high-dimensional models are learned that inherit sparsity from, e.g., finite-element discretizations of the governing equations of the systems of interest. In contrast, we aim to learn low-dimensional models that help to reduce computational costs in applications that require many model evaluations .

Instead of aiming to find models that best-fit data, we aim to exactly recover reduced models from data so that our models inherit the reduced models’ well-studied properties. Our approach is based on operator inference , which has been derived from and is a data-driven model reduction approach that learns approximations of reduced models from state trajectories. In , operator inference has been introduced for systems with polynomial nonlinear terms and in operator inference is combined with the transform & learn approach to obtain models of systems with more general nonlinear terms. Operator inference projects trajectories of systems of interest onto low-dimensional subspaces of the high-dimensional state spaces and then fits operators to the projected trajectories via least-squares regression. However, as is known from, e.g., the Mori-Zwanzig formalism from statistical physics , the projected trajectories correspond to non-Markovian dynamics in the low-dimensional subspaces even though the high-dimensional trajectories and the corresponding high-dimensional systems are Markovian. The non-Markovian dynamics are related to the closure error in model reduction . To account for the non-Markovian dynamics, methods have been proposed that learn non-Markovian terms and that use time-delay and other embeddings ; however, since we aim to exactly recover the Markovian reduced models that are obtained with traditional model reduction, neither of these remedies are applicable in our situation. Instead, we propose a data sampling scheme that iterates between time stepping the high-dimensional systems and projecting onto low-dimensional subspaces to generate trajectories that correspond to low-dimensional Markovian dynamics. We then show that, under certain conditions, applying operator inference to these re-projected trajectories gives the same operators that are obtained with traditional model reduction methods. The result is a pre-asymptotic guarantee to exactly recover reduced models from finite amounts of data for a wide class of systems with polynomial nonlinear terms. Our numerical results demonstrate these theoretical results in practice by learning low-dimensional models that match the reduced models from traditional model reduction up to numerical errors.

Section 2 discusses preliminaries on dynamical systems, traditional model reduction, operator inference, and formulates the problem. Section 3 introduces data sampling with re-projection to obtain trajectories that correspond to low-dimensional Markovian dynamics and provides an analysis that shows that operators fitted to these re-projected trajectories are the operators obtained with traditional model reduction. The overall computational approach is presented in Algorithm 2 in Section 4 and numerical results are given in Section 5. Conclusions are drawn in Section 6.

Preliminaries

The focus of this work is on dynamical systems with polynomial nonlinear terms, which we introduce in Section 2.1 together with traditional model reduction for these systems in Section 2.2. A building block of our approach is operator inference for learning reduced models from data, which we discuss in Section 2.3. The problem we aim to address is formulated in Section 2.4.

2 Model reduction of systems with polynomial nonlinear terms

For $j=1,\dots,m$, the reduced operators are constructed via, e.g., Galerkin projection

3 Operator inference

Operator inference proceeds in three steps. First, state trajectories $\bm{X}(\bm{\mu}_{1}),\dots,\bm{X}(\bm{\mu}_{m})$ and $\bm{Y}(\bm{\mu}_{1}),\dots,\bm{Y}(\bm{\mu}_{m})$ are obtained by querying the system (1) at parameters $\bm{\mu}_{1},\dots,\bm{\mu}_{m}\in\mathcal{D}$ to derive a reduced space spanned by the columns of $\bm{V}_{n}=[\bm{v}_{1},\dots,\bm{v}_{n}]$. Many of the basis construction techniques developed in traditional model reduction can be applied; see references given in Section 2.2. In the following, we will use POD to construct $\bm{V}_{n}$ as described in Section 2.2. The second step of operator inference is to project the trajectories onto the reduced space $\mathcal{V}_{n}$ spanned by the columns of $\bm{V}_{n}$ and so to obtain the projected trajectories

In the third step of operator inference, the operators

3.2 Data matrix

It will be convenient to write (7) for each $j=1,\dots,m$ as

4 Problem formulation

By fitting operators to projected trajectories with operator inference as described in Section 2.3 and in , the closure error (11) is introduced into the learned operators, which means that the learned operators can fail to approximate the dynamics of the intrusive reduced model.

Sampling Markovian dynamics via re-projection

In this section, we focus on learning reduced models corresponding to a single parameter $\bm{\mu}_{j}$, which then is subsequently repeated for all parameter $j=1,\dots,m$. To ease exposition, we drop the dependence on $\bm{\mu}_{j}$ in this section.

which gives with an inductive argument that

2 Data sampling with re-projection to avoid non-Markovian dynamics

We now describe our sampling scheme with re-projection. Consider an initial condition $\bm{x}_{0}\in\mathcal{V}_{n}$ and set $\bar{\bm{x}}_{0}=\bm{V}_{n}^{T}\bm{x}_{0}$. Note that $\bm{V}_{n}\bar{\bm{x}}_{0}=\bm{x}_{0}$ because $\bm{x}_{0}\in\mathcal{V}_{n}$. Our scheme proceeds iteratively, see Figure 2. In the first iteration, system (1) is queried at initial condition $\bm{V}_{n}\bar{\bm{x}}_{0}$ and input $\bm{u}_{0}$ to obtain

Then, the re-projected state $\bar{\bm{x}}_{1}=\bm{V}_{n}^{T}\bm{x}_{\text{tmp}}$ is computed by projecting $\bm{x}_{\text{tmp}}$ onto $\mathcal{V}_{n}$. In the second iteration, system (1) is queried for a single time step at the initial condition $\bm{V}_{n}\bar{\bm{x}}_{1}$ and input $\bm{u}_{1}$ to obtain $\bm{f}(\bm{V}_{n}\bar{\bm{x}}_{1},\bm{u}_{1})$ and to compute $\bar{\bm{x}}_{2}$ via projection $\bar{\bm{x}}_{2}=\bm{V}_{n}^{T}\bm{f}(\bm{V}_{n}\bar{\bm{x}}_{1},\bm{u}_{1})$. This process is repeated to generate the re-projected states $\bar{\bm{x}}_{0},\bar{\bm{x}}_{1},\dots,\bar{\bm{x}}_{K}$ and to collect them into the re-projected trajectories $\bar{\bm{X}}=[\bar{\bm{x}}_{0},\bar{\bm{x}}_{1},\dots,\bar{\bm{x}}_{K-1}]$ and $\bar{\bm{Y}}=[\bar{\bm{x}}_{1},\dots,\bar{\bm{x}}_{K}]$.

Algorithm 1 summarizes our data sampling scheme with re-projection. The inputs to Algorithm 1 are the high-dimensional system $\bm{f}$, a basis matrix $\bm{V}_{n}$, an initial condition $\bm{x}_{0}\in\mathcal{V}_{n}$, a parameter $\bm{\mu}\in\mathcal{D}$, and inputs $\bm{u}_{0},\dots,\bm{u}_{K-1}$. Line 2 projects the initial condition $\bm{x}_{0}$ to obtain $\bar{\bm{x}}_{0}$. The for loop on line 3 iterates over the time steps $k=0,\dots,K-1$ and generates the re-projected state $\bar{\bm{x}}_{k+1}$ by querying the high-dimensional system for a single time step in line 4. The re-projected trajectories $\bar{\bm{X}}$ and $\bar{\bm{Y}}$ are returned in line 7.

3 Exact recovery of reduced models from re-projected trajectories

derived from the re-projected trajectory $\bar{\bm{X}}$, cf. the data matrix $\breve{\bm{D}}$ derived from the projected trajectory $\breve{\bm{X}}$ defined in (10). If $\bar{\bm{D}}$ has full rank, then the least-squares problem

Computational procedure and practical aspects

This section summarizes the overall computational procedure of operator inference with re-projected trajectories in Algorithm 2 and discusses practical aspects as well as limitations of the approach.

The computational costs of Algorithm 2 are typically dominated by querying the high-dimensional system. The costs of assembling the data matrix on line 8 and the costs of solving the corresponding least-squares problem on line 9 typically are negligible. In the for loop in line 2, the high-dimensional system is time stepped to generate the trajectories for constructing the POD basis matrix, which is similar to traditional, intrusive model reduction. The for loop in line 6 requires time stepping the high-dimensional systems once more to sample the re-projected trajectories with Algorithm 1. Thus, the computational costs of learning a reduced model with operator inference with re-projection is twice as high as the costs of constructing a model with operator inference without re-projection. Note, however, that it is unnecessary to sample re-projected trajectories of length $K$. Sampling shorter re-projected trajectories can significantly reduce the computational costs of operator inference with re-projection.

2 Practical aspects and condition of least-squares problem

We make three remarks of practical aspects of operator inference with re-projection. First, Corollary 1 states that operator inference from re-projected trajectories gives the intrusive reduced models if condition (15) is satisfied and if the data matrix $\bar{\bm{D}}$ defined in (16) has full rank. It is straightforward to numerically verify these two conditions in practice and so to determine if Corollary 1 applies and if the intrusive reduced model is obtained up to numerical errors.

Second, to sample the re-projected trajectories with Algorithm 1, it is necessary to have available the high-dimensional system in the sense that it can be time stepped for a single time step with initial condition $\bar{\bm{x}}_{k}$ for $k=0,\dots,K-1$. This is in contrast to operator inference without re-projection, which is applicable even if only the trajectories $\bm{X}(\bm{\mu}_{1}),\dots,\bm{X}(\bm{\mu}_{m})$ and the corresponding inputs $\bm{U}(\bm{\mu}_{1}),\dots,\bm{U}(\bm{\mu}_{m})$ are available and the high-dimensional system cannot be queried. However, note that it is unnecessary to time step the high-dimensional system at arbitrary initial conditions. The re-projected states are close to the states of the high-dimensional system if the space $\mathcal{V}_{n}$ is sufficiently rich, which typically is a necessary requirement for the success of model reduction in any case.

and then use (19) and $\bar{\bm{Y}}(\bm{\mu}_{j})$ obtained from $\bar{\bm{Y}}_{1}(\bm{\mu}_{j}),\dots,\bar{\bm{Y}}_{m^{\prime}}(\bm{\mu}_{j})$ in the least-squares problem (17) to learn a model. This is a similar process as suggested in .

Numerical results

The numerical results in this section demonstrate that the proposed data sampling strategy with re-projection leads to low-dimensional models that match reduced models derived with traditional model reduction methods up to numerical errors in practice. The toy example introduced in the problem formulation in Section 2.4 is revisited in Section 5.1. Section 5.2 derives models for the viscous Burgers’ equation and Section 5.3 for the Chafee-Infante equation. Both of these examples are one dimensional in the spatial domain. Section 5.4 demonstrates learning models from re-projected trajectories on a diffusion-reaction equation with two spatial dimensions.

We revisit the toy example introduced in Section 2.4. Let $\bar{\bm{X}}$ be the re-projected trajectory obtained with Algorithm 1. Following the least-squares problem (17) described in Corollary 1, we learn a model from the re-projected trajectory $\bar{\bm{X}}$ and time step the learned model to obtain the trajectory $\hat{\bm{X}}$, which is plotted in Figure 3a. The trajectory of the model learned from the re-projected trajectory closely follows the trajectory of the intrusive reduced model, which is in stark contrast to the model learned from the trajectory $\breve{\bm{X}}$ without re-projection. Thus, the results in Figure 3a are in agreement with Corollary 1.

Now consider the data matrix $\bar{\bm{D}}$ defined in (16). Figure 3b shows the condition number of $\bar{\bm{D}}^{T}\bar{\bm{D}}$ for dimensions $n\in\{2,4,6\}$ and various numbers of time steps $K$. In this example, the condition number grows with the dimension $n$. This means that even though condition (15) together with a full-rank data matrix are sufficient to recover the intrusive reduced model, numerical errors are introduced into the learned operators because of the potentially high condition number of $\bar{\bm{D}}^{T}\bar{\bm{D}}$; cf. Section 4.2. Figure 3c demonstrates that the difference

2 Burgers’ equation

A similar setup as in is used for demonstrating the proposed approach on the viscous Burgers’ equation.

2.2 Results

is plotted in Figure 4c. The models learned from re-projected trajectories achieve similar errors as the intrusive reduced models, in contrast to models learned from trajectories without re-projection. Similar observations can be made for $\bar{n}=15$ as shown in Figure 4b for training parameters and training inputs and in Figure 4d for test parameters and test inputs.

between the trajectories of the intrusive reduced models and the trajectories computed with the learned models. Thus, $\bm{Z}(\mu_{i}^{\text{test}})$ in (24) is either the trajectory obtain with $\hat{\bm{f}}(\cdot,\cdot;\mu_{i}^{\text{test}})$ or with $\breve{\bm{f}}(\cdot,\cdot;\mu_{i}^{\text{test}})$ for $i=1,\dots,m_{\text{test}}$. The difference (24) is plotted in Figure 5. The models learned from re-projected trajectories achieve a difference to the intrusive reduced model of less than $10^{-10}$, whereas the models learned from trajectories without re-projection are up to 8 orders of magnitude worse in terms of difference (24) and diverge in most cases (missing values in the plots).

3 Chafee-Infante equation

A similar setup as in is used in this section.

with the spatial coordinate $\xi\in\Omega$ and time $t\in[0,T]$. Note that we consider a parameter-independent version of the Chafee-Infante equation. The boundary conditions are

for $K=4\times 10^{5}$ and $N=128$ and where the input matrix $\bm{B}$ corresponds to the discretization of the boundary conditions.

3.2 Results

4 Diffusion-reaction equation

The setup of the diffusion-reaction equation in this section follows the example in .

for $K=10^{4}$. The dimension $N$ of the state $\bm{x}_{k}$ at time step $k$ is $N=64^{2}=4096$. Plots of $\bm{x}_{K}(\mu)$ for $\mu=1.0625$ and $\mu=1.4375$ are given in Figure 8.

To construct a reduced space, we select $m=10$ equidistant parameters $\mu_{1},\dots,\mu_{m}\in\mathcal{D}$ and set the inputs to be realizations of the random variables uniformly distributed in $$. From these trajectories, the basis matrix$\bm{V}_{\bar{n}}$with$\bar{n}=10$columns is computed with POD. Then, re-projected trajectories are sampled up to time$t=5$(instead of end time$T=100$). For each$\mu_{i}$, 10 re-projected trajectories with different random inputs are derived, and concatenated together as described in Section 4.2. The concatenation of trajectories ensures that the data matrix$\bar{\bm{D}}$has full rank in this example. Models are learned with operator inference from the re-projected trajectories to obtain$\hat{\bm{f}}(\cdot,\cdot;\mu_{1}),\dots,\hat{\bm{f}}(\cdot,\cdot;\mu_{m})$. The same process is repeated for the trajectories without re-projection to obtain the models$\breve{\bm{f}}(\cdot,\cdot;\mu_{1}),\dots,\breve{\bm{f}}(\cdot,\cdot;\mu_{m})$. The rest of the setup is the same as in Section 5.2. Test parameters are 7 equidistantly chosen parameters in$\mathcal{D}$. Test inputs are realizations of random variables with uniform distribution in$$.

4.2 Results

Figure 9a shows the error (22) for the training parameters and training inputs. The model fitted to trajectories without re-projection numerically diverged to NaNs during time stepping for all dimensions $n>2$. The model fitted to re-projected trajectories closely matches the behavior of the intrusive reduced model as expected from the analysis presented in Corollary 1. The same observations can be made for the error (23) with the test parameters and test inputs.

Conclusions

The presented approach exactly recovers reduced models from data under certain conditions. This result holds pre-asymptotically in the number of data points and the dimension of the reduced space as long as the corresponding data matrix is full rank. The optimization problem underlying operator inference with re-projected trajectories is convex and can be solved with standard numerical linear algebra packages. Numerical experiments demonstrate that reduced models are learned up to numerical errors in practice for a wide class of systems with polynomial nonlinear terms.

Acknowledgments

The author would like to thank Elizabeth Qian, Nihar Sawant, and Karen Willcox for many helpful discussions. This work was partially supported by US Department of Energy, Office of Advanced Scientific Computing Research, Applied Mathematics Program (Program Manager Dr. Steven Lee), DOE Award DESC0019334. The numerical experiments were computed with support through the NYU IT High Performance Computing resources, services, and staff expertise.

References