Data-driven identification of parametric partial differential equations

Samuel Rudy, Alessandro Alla, Steven L. Brunton, J. Nathan Kutz

Introduction

The extraction of physical laws from experimental data, often in the form of differential and partial differential equations, may be critical to science and engineering applications where governing equations are unknown. Time-series data collected from experiments, or as the macroscopic aggregate of small scale behavior, often obeys unknown governing equations that are parametrized by time-evolving parameters $\mu(t)$ . In some cases, it may be possible to derive physical laws from first principals using data as well as some knowledge of the system, but there are many cases where this is elusive such as the large scale networked dynamical systems of the power grid and the brain, or chemical kinetics of, for example, the Belousov-Zhabotinsky reaction. Recently, there has been a substantial research effort towards automating the process of data-driven model discovery in order to identify interpretable expressions for the dynamics in the form of ordinary and partial differential equations (ODEs and PDEs):

Figure 1 demonstrates two prototypical parametric dependencies: (a) a PDE model whose constant parameters change at fixed points in time and (b) a PDE model that depends continuously on the parameter $\mu(t)$ . Our proposed model discovery method provides a principled approach to efficiently handle these two parametric cases, thus advancing the field of model discovery by disambiguating the dynamics from parametric time dependence. Even if the governing equations are known, parametric dependencies in PDEs complicate numerical simulation schemes and challenge one’s ability to produce accurate future state predictions. Characterizing parametric dependence is also an critical task for model reduction in both time-independent and time-dependent PDEs . Thus the ability to explicitly extract the parametric dependence of a spatio-temporal system is necessary for accurate, quantitative characterization of PDEs.

More broadly, system identification using machine learning methods has emerged as a viable alternative to expert knowledge and first principles derivations. It is important to separate the field of system identification into two distinct categories: (i) methods that accurately reflect observed dynamics using black box functions (e.g. neural networks), and (ii) methods that recover closed form and interpretable expressions for the dynamics in the form of ordinary and partial differential equations (ODEs and PDEs). This duality reflects the two cultures narrative of machine learning and classical statistics made popular by Leo Breiman . On one hand, the research may assume a specific model for the data with known mechanism, while on the other the research is interested in algorithmic models that, while not necessarily reflecting the true mechanism, are accurate in prediction. While several recent works have made progress from both viewpoints, we focus on the former. The terms in a differential equation often have physical interpretations and motivation, e.g. diffusion and advection are ubiquitous in many physical systems and are characterized by prototypical expressions in a PDE. For systems where first principal derivations prove intractable, we may gain insights into the underlying physics of the system based on the terms in the identified PDE. Further, we view the process of extracting closed form equations as being more generalizable than fitting black box models to a specific dataset. Specifically, altering initial conditions will not change governing equations, but may break a machine learned black box solver.

Research towards the automated inference of dynamical systems from data is not new . Methods for extracting linear systems from time series data include the eigensystem realization algorithm (ERA) and Dynamic Mode Decomposition (DMD) . Identification of nonlinear systems has, until very recently, relied on black box methods. These include NARMAX , neural networks , equation free methods , and Laplacian spectral analysis . There has also been considerable recent work towards data-driven approximation of the Koopman operator via extensions of DMD , diffusion maps , delay-coordinates and neural networks .

The use of genetic algorithms for nonlinear system identification allowed for the derivation of physical laws in the form of ordinary differential equations. Genetic algorithms are highly effective in learning complex functional forms but are slower computationally than simple regression. Sparsity promoting methods have been used previously in dynamical systems , and sparse regression has been leveraged to identify parsimonious ordinary differential equations from a large library of allowed functional forms . Much work building on the sparse regression framework has followed and includes inferring rational functions , the use of information criteria for model validation , constrained regression for conservation laws , model discovery with highly corrupt data , the learning of bifurcation parameters , stochastic dynamics , weak forms of the observed dynamics , and regression with small amounts of data . In contrast to sparse regression, a neural network based approach was proposed to identify ordinary differential equations, sacrificing some interpretability for a richer class of allowed functional forms .

Sparse regression based methods for PDEs were first used in . These methods were demonstrated on a large class of PDEs and have the benefit of being highly interpretable, but struggle with numerical differentiation of noisy data. In Rudy et al. the noise was addressed by testing with only a small degree of noise (large SNR), while in Schaeffer et al. noise was added to the time derivative after it was computed from clean data. Alternatively, Gaussian processes were used to determine linear PDEs and nonlinear PDEs known up to a set of coefficients . Using Gaussian process regression requires less data than sparse regression and naturally manages noise, but the method is only applicable to PDEs with a known structure. Reference makes a substantial contribution by using neural networks to accurately learn partial differential equations with non constant coefficients. A neural network is constructed that mimics a forward Euler timestepping scheme and the accuracy of potential models is evaluated based on their future state prediction accuracy. While seemingly more robust than sparse regression, the method in does not penalize extraneous terms in the learned PDE and thus falls short of producing optimally parsimonious models. Furthermore, only tests the method on a nonlinear problem using a relatively strong Ansatz. Neural networks were also used in and to solve and to estimate parameters in partial differential equations with known terms to a high degree of accuracy. However, similar to , it is assumed that the PDE is known up to a set of coefficients. A more sophisticated neural network approach was used in to learn dynamics of systems with unknown terms. However, the approach in does not give closed form representations of the dynamics and the resulting neural network model therefore does not give insights into the underlying physics.

In this work, we present a sparse regression framework for identifying PDEs with non-constant coefficients, something that none of the previous PDE discovery methods are equipped to do. Specifically, we allow for variation in the value of a coefficient across time or space, but maintain that the active terms in the PDE must be consistent. This is an important innovation in practice as the parameters of physical systems often vary during the measurement process, so that the parametric dependencies be disambiguated from the PDE itself. Our method extends the sparse regression frameworks first proposed for PDE discovery by using group sparsity and results in a more parsimonious and interpretable model than neural networks. We are still limited by the accuracy of numerical differentiation and by the library terms in the sparse regression. Numerical differentiation using neural networks as shown in appears promising as a method for obtaining more accurate time derivatives from noisy data. The limitation based on terms included in the library seems more permanent. Any closed-form model expression must be representable with a finite set of building blocks. Here we only use monomials in our data and its derivatives, since these are the common terms seen in physics, but there is no limitation to the terms included in the library.

Methods

The parametric discovery method relies on several foundational mathematical tools. In the following subsections, we will discuss the identification of constant coefficient equations as well as regression methods for group sparsity. Finally, we combine these ideas to show how one may identify parametric PDEs and suggest a methodology for model selection that balances accuracy and the number of active terms in the PDE.

Several recent methods have been proposed for the identification of constant coefficient partial differential equations from data. In this work we expand on the sparse regression framework, PDE-FIND, used in . We will briefly elaborate on the method and refer the reader to the original paper for details. The PDE-FIND algorithm provides a principled technique for discovering the underlying PDE from spatial time series measurements alone using a library of candidate functions for the PDE and sparse regression. For the identification of constant coefficient PDEs, we have a dataset ${\bf U}$ , which is a discretization of a function $u(x,t)$ that we assume satisfies the PDE of the form given in (1):

We assume that the nonlinear expression $N(\cdot)$ may be expanded as a sum of simple monomial basis functions $N_{j}$ of $u$ and its derivatives. Note that this sum is not unique and that we can include extra basis functions by simply setting the corresponding $\xi_{j}$ to be zero. In PDE-FIND, which constructs an overcomplete library of many possible monomial basis functions and regresses to find $\xi$ , sparsity is used to ensure that basis functions that do not appear in the PDE are set to zero in the sum.

which is a large, overdetermined linear system of equations ${\bf A}{\bf x}={\bf b}$ . Note that here we have shown a problem where derivatives up to third order are multiplied by powers of $u$ up to cubic order, but one could include arbitrarily many library functions. Solving for $\xi$ and ensuring sparsity gives the PDE. PDE-FIND has been shown to accurately identify several partial differential equations from data alone. The sparsity constraint is a regularizer for the linear regression .

2 Group Sparsity

In a typical sparse regression, we seek a sparse solution to the linear system of equations ${\bf A}{\bf x}={\bf b}$ . Accuracy of the predictor, $\|$ Ax-b $\|$ , is balanced against the number of nonzero coefficients in ${\bf x}$ . Thus the sparse regularization enforces a solution ${\bf x}$ with many zeros (which is the variable $\boldsymbol{\xi}$ in (3)). In this paper, we use the notion of group sparsity to find time series representing each parameter in the PDE, rather than single values. We group collections of terms in ${\bf x}$ together and seek solutions to ${\bf A}{\bf x}={\bf b}$ that minimize the number of groups with nonzero coefficients.

One well studied method for solving regression problems with group sparsity is the group LASSO (GLASSO)

Here $\mathcal{G}$ is a collection of groups, each of which contains a subset of the indices enumerating the columns of ${\bf A}$ and coefficients in ${\bf x}$ . Note that the second term in the GLASSO corresponds to a convex relaxation of the number of groups containing a nonzero value.

The concept of group sparsity has been used in several previous methods for identifying dynamics given by ordinary differential equations . As it will be shown in the following, we find that the GLASSO performs poorly in the case of identifying PDEs. We instead use a sequential thresholding method based on ridge regression, similar to the method used in , but adapted for group sparsity. A sequential thresholding method was also used in for group sparsity but for ordinary and not partial differential equations. Our method, which we call Sequential Grouped Threshold Ridge Regression (SGTR), is summarized in algorithm 1.

Throughout the training, $\mathcal{G}$ tracks the groups that have nonzero coefficients, and it is paired down as we threshold coefficients with sufficiently small relevance, as measured by $f$ . We use the 2-norm of the coefficients in each group for $f$ but one could also consider arbitrary functions. In particular, for problems where the coefficients within each group have a natural ordering, as they do in our case as time series or spatial functions, one could consider smoothness or other properties of the functions. In practice, we normalize each column of ${\bf A}$ and ${\bf b}$ so that differences in scale between the groups do not affect the result of the algorithm. For the GLASSO we always perform an unregularized least squares regression on the nonzero coefficients after the sparsity pattern has been discovered to debias the coefficients. We found SGTR to outperform the GLASSO for the problem of correctly identifying the active terms in parametric PDEs.

3 Data-driven identification of parametric partial differential equations

In the identification of parametric PDEs, we consider equations of the form

Note that this equation is similar to (2) but has time-varying parametric dependence. To capture spatial variation in the coefficients, we simply replace $\xi(t)$ with $\xi(x)$ . The PDE is assumed to contain a small number of active terms, each with a time varying coefficient $\xi(t)$ . We seek to solve two problems: (i) determine which coefficients are nonzero and (ii) find the values of the coefficients for each $\xi_{j}$ at each timestep or spatial location for which we have data.

For time dependent problems, we construct a separate regression for each timestep, allowing for variation in the PDE between timesteps. Similar to the PDE-FIND method, we construct a library of candidate functions for the PDE using monomials in our data and its derivatives so that

where the set of $m$ equations is given by

Our goal is to solve the set of equations given by (7) with the constraint that each $\xi^{(j)}$ is sparse and that they all share the same sparsity pattern. That is, we want a fixed set of active terms in the PDE. To do this, we consider the set of equations as a single linear system and use group sparsity. Expressing the system of equations for the parametric equation as a single linear system we get the block diagonal structure given by

We solve (8) using SGTR with columns of the block diagonal library matrix grouped by their corresponding term in the PDE. Thus for $m$ timesteps and $d$ candidate functions in the library, groups are defined as $\mathcal{G}=\{{j+d\cdot i:i=1,\dots,m}:j=1,\ldots,d\}$ . This ensures a sparse solution to the PDE while also allowing arbitrary time series for each variable. To obtain the correct level of sparsity, we separate 20% of the data from each timestep to use as a validation set and search over the parameter $\lambda$ in the SGTR algorithm using cross validation to find the optimal $\lambda$ . For problems with spatial, rather than temporal, variation in the coefficients, we simply group by spatial rather than time coordinate. A similar block diagonal structure is obtained but with $n$ blocks of size $m\times d$ rather than $m$ blocks of size $n\times d$ . The groups are defined by $\mathcal{G}=\{{j+d\cdot i:i=1,\dots,n}:j=1,\ldots,d\}$ .

Since we are evaluating the relevance of groups based on their norm, it is important to consider differences in the scale of the candidate functions. For example, if $u\sim\mathcal{O}(10^{-2})$ then a cubic function will be $\mathcal{O}(10^{-6})$ and relatively large coefficients multiplying this data may not have a large effect on the dynamics, but will not be removed by the hard threshold due to its size. To remedy this, we normalize each candidate function represented in $\mathbf{\Theta}$ as well as each $u_{t}^{(j)}$ to have unit length prior to the group thresholding algorithm and then correct for the normalization after we have discovered the correct sparsity pattern.

4 Model Selection

We check 50 evenly spaced values of $\lambda$ between $10^{-5}\lambda_{\mbox{max}}$ and $\lambda_{\mbox{max}}$ on a logarithmic scale.

For SGTR, we search over the range of tolerances between $\epsilon_{min}$ and $\epsilon_{max}$ defined as

To select the optimal model generated via each method, we evaluate the models using the AIC-inspired loss function

where $k$ is the number of nonzero coefficients in the identified PDE, $\|\xi\|_{0}/m$ , and $N$ is the number of rows in $\mathbf{\Theta}$ , which is equal to the size of our original dataset $u$ .

Equation (11) is closely related to the Akaike Information Criterion (AIC) . Typically, the mean square error of a linear model is used to evaluate goodness of fit, but in our case there is error in computing the time derivative $u_{t}$ , so we assume that any linear model which perfectly fits the data is overfit. We have added $\epsilon=10^{-5}$ to the mean square error of each model as a floor in order to avoid overfitting. Without this addition, our algorithm selects insufficiently parsimonious representations of the dynamics.

Figure 2 illustrates the loss function from equation (11) evaluated on models derived from 50 values of $\epsilon$ and $\lambda$ using SGTR and GLASSO respectively. Initially, a low penalty in each algorithm yields a model that is overfit to the data given our sparsity criteria. For an intermediate value of $\epsilon$ or $\lambda$ , a more parsimonious but still predictive model is obtained. For sufficiently high values, the model is too sparse and is no longer predictive.

Computational Results of Parametric PDE Discovery

To test the parametric discovery of PDEs, we consider a solution of Burgers’ equation with a sinusoidally oscillating coefficient $a(t)$ for the nonlinear advection term

where a small amount of diffusion is added to regularize the evolution dynamics.

The time dependent Burgers’ equation was solved numerically using a spectral method on the interval $ $with periodic boundary conditions and$ t\in $with$ n=256 $grid points and$ m=256 $time steps. We search for parsimonious representations of the dynamics by including powers of$ u $up to cubic order, which can be multiplied by derivatives of$ u$ up to fourth order. For the noise-free dataset we use the discrete Fourier transform for computing derivatives. For the noisy dataset, we use polynomial interpolation to smooth the derivatives .

The resulting time series for the identified nonzero coefficients are shown in Fig. 4. SGTR correctly identified the active terms in the PDE for both the noise-free and noisy datasets, whereas GLASSO fails in both cases to produce the correct PDE and its parametric dependencies.

2 Navier-Stokes: Flow around a cylinder

We consider the fluid flow around a circular cylinder by simulating the Navier-Stokes vorticity equation

Data is generated using the Immersed Boundary Projection Method (IBPM) with $n_{x}=449$ and $n_{y}=199$ spatial points in $x$ and $y$ respectively, and 1000 timesteps with $dt=0.02$ . The Reynolds number is adjusted half way through the simulation from $\nu=100$ initially to $\nu=75$ . This is representative of the fluid velocity exhibiting a sudden decrease midway through the data collection. Our library of candidate functions is constructed using up to second order derivatives of the vorticity and multiplying by up to quadratic functions of the data. To keep the size of the machine learning problem tractable, we subsample 1000 random spatial location from the wake of the cylinder to construct our library at every tenth timestep . For the noise-free dataset, far fewer points are needed to accurately identify the dynamics. We suspect that with a more careful treatment of the numerical differentiation in the case of noisy data, such as that used in , the same would be true for the dataset with artificial noise, however such work is not the focus of this paper. The identified time series for the Navier-Stokes equation are shown in Fig. 6. SGTR and GLASSO both correctly identify the active terms in the PDE.

3 Spatially Dependent Advection-Diffusion Equation

The advection-diffusion equation is a simple model for the transport of a physical quantity in a velocity field with diffusion. Here, we adapt the equation to have a spatially dependent velocity

which models transport through a spatially varying vector field due to $c=c(x)$ . The PDE is solved on a periodic domain $[-L,L]$ with $L=5$ , $\epsilon=0.1$ , and $c(x)=-1.5+\cos(2\pi x/L)$ using a spectral method with $n=256$ and $m=256$ . The library consists of powers of $u$ up to cubic, multiplied by derivatives of $u$ up to fourth order. Results for the advection-diffusion equation are shown in Fig. 8. In the noise-free and noisy datasets, both SGTR and GLASSO Correctly identify the active terms in the PDE.

4 Spatially Dependent Kuramoto-Sivashinsky Equation

We now test the method on a Kuramoto-Sivashinsky equation with spatially varying coefficients

We use a periodic domain $[-L,L]$ with $L=20$ and coefficients $a(x)=1+\sin(2\pi x/L)/4$ , $b(x)=-1+e^{-(x-2)^{2}/5}/4$ and $c(x)=-1-e^{-(x+2)^{2}/5}/4$ .

The equation is solved numerically to $t=200$ using $n=512$ grid points and $m=1024$ timesteps. The second half of the data set is used, so as to only consider the region where the dynamics exhibited spatio-temporal chaos, resulting in a dataset containing 512 snapshots of 512 gridpoints. Since the Kuramoto Sivashinsky equation involves a fourth order derivative, it is very difficult to correctly identify it with noisy data since it is exceptionally difficult to accurately compute the fourth derivative. Indeed, our method fails to correctly identify the active terms when 1% noise is added. With 0.01% noise the correct terms were identified but with substantial error in coefficient value. We suspect that this shortcoming could at least be partially remedied by a more careful treatment of the numerical differentiation such as in . The results of our parametric identification are shown in Fig. 10.

Discussion

We have presented a method for identifying governing laws for physical systems which exhibit either spatially or temporally dependent behavior. The method builds on a growing body of work in the applied mathematics and machine learning community that seeks to automate the process of discovering physical laws. To the best of our knowledge, our method is the first approach for deriving parsimonious PDE expressions of spatio-temporal system in the case of non-constant coefficients. Specifically, we can disambiguate between the governing PDE evolution and its parametric dependencies. In all examples, the SGTR algorithm outperformed the GLASSO in correctly identifying active terms in the PDE. Errors from the latter were generally in the form of extra terms added with small coefficient values throughout the time series. It may seem reasonable to threshold these time series after the discovery algorithm, but doing so assumes that the terms importance in the PDE is directly related to its magnitude, an assumption which we do not make given the normalization prior to sparse regression.

In this work we have split the data into distinct timesteps or spatial locations in order to find PDE models for each subset of the data, resulting in coefficients that can vary in space or time. However, with a sufficiently fine grid, it seems feasible that one could bin the data by areas localized in space and time to determine a coefficient varying in both space and time with some loss of resolution. This same result may be achievable in a more stable manner by introducing a sparsity term to the work in .

As is the case with other sparse regression methods for identifying dynamical systems, this method is constrained by the ability of the user to accurately differentiate data. For ordinary differential equations, this may be circumvented by looking at the weak form of the dynamics , but doing so for PDEs seems difficult since there are derivatives that need to be evaluated with respect to multiple variables. We find the automatic differentiation approach used in promising and suspect that the inclusion of neural network based differentiation could radically improve the ability of our method to identify dynamics from noisy data. With sufficient knowledge of data it may also be possible to obtain better estimates through tuning the polynomial based differentiation .

Automating the identification of closed form physical laws from data will hopefully boost scientific progress in areas where deriving the same laws from first principals proves intractable. There are several limitations to many methods proposed in the field thus far. In particular, current methods have generally studied equations of the form $u_{t}=N(u,x,t)$ but many equations in physics are not in this class. Indeed, if measuring a system with parametric dependencies, then past methods are be unable to disambiguate between the evolution dynamics and its parametric dependencies $\mu(t)$ , thus greatly limiting model discovery. There is also a trade off between methods that are able to derive parsimonious representations, though which are limited to a finite set of library elements, and those that use black box models to represent larger classes of possible functions. The researcher may also find difficulties in attempting to infer dynamics from the wrong set of measurements. For example, one could not derive the Schrödinger by only looking at measurements of intensity. While not addressing these issues, this work makes a step towards generalizing the class of equations which may be accurately identified via machine learning methods.

Code: https://github.com/snagcliffs/parametric-discovery