Robust and optimal sparse regression for nonlinear PDE models

TL;DR

This paper presents a robust sparse regression method using weak formulation to accurately identify nonlinear PDE models from noisy data, outperforming existing techniques especially under varying noise conditions.

Contribution

It introduces an optimized sparse regression approach with a scaling relation linking model accuracy to data resolution and noise levels, specifically for nonlinear PDEs.

Findings

Achieves significantly better accuracy than existing methods.

Derives a scaling relation between model accuracy, data resolution, and noise.

Demonstrates effectiveness on the Kuramoto-Sivashinsky equation.

Abstract

This paper investigates how models of spatiotemporal dynamics in the form of nonlinear partial differential equations can be identified directly from noisy data using a combination of sparse regression and weak formulation. Using the 4th-order Kuramoto-Sivashinsky equation for illustration, we show how this approach can be optimized in the limits of low and high noise, achieving accuracy that is orders of magnitude better than what existing techniques allow. In particular, we derive the scaling relation between the accuracy of the model, the parameters of the weak formulation, and the properties of the data, such as its spatial and temporal resolution and the level of noise.

Equations

$$\sum _{n=1}^N{c}_n{\mathbf{f}}_n(\mathbf{x},t,\mathbf{u},{\partial }_t\mathbf{u},{\partial }_t^2\mathbf{u},\nabla \mathbf{u},{\nabla }^2\mathbf{u},\cdots )=0$$

$$Q\mathbf{c}=0$$

$$\mathbf{c}=\underset{\Vert \mathbf{c}\Vert =1}{\mathrm{argmin}}\Vert Q^{T}Q\mathbf{c}\Vert .$$

$$c_1{\partial }_{t}u+c_{2}u{\partial }_{x}u+c_3{\partial }_x^{2}u+c_4{\partial }_x^{4}u=0$$

$$q_1^{jk}={\int }_{{\Omega }_k}{w}_j{\partial }_{t}ud\Omega ,q_2^{jk}={\int }_{{\Omega }_k}{w}_{j}u{\partial }_{x}ud\Omega ,$$

$$q_3^{jk}={\int }_{{\Omega }_k}{w}_j{\partial }_x^{2}ud\Omega ,q_4^{jk}={\int }_{{\Omega }_k}{w}_j{\partial }_x^{4}ud\Omega ,$$

$$q_1^{jk}=-{\int }_{{\Omega }_k}u{\partial }_t{w}_{j}d\Omega ,q_2^{jk}=-\frac{1}{2}{\int }_{{\Omega }_k}{u}^2{\partial }_x{w}_{j}d\Omega ,$$

$$q_3^{jk}={\int }_{{\Omega }_k}u{\partial }_x^2{w}_{j}d\Omega ,q_4^{jk}={\int }_{{\Omega }_k}u{\partial }_x^4{w}_{j}d\Omega$$

$$c_4(\mathbf{x},t)=\sum _p{c}_p^{\prime }{g}_p(\mathbf{x},t)$$

$${\int }_{{\Omega }_k}{w}_j{c}_4{\partial }_x^{4}ud\Omega =\sum _p{c}_p^{\prime }{q}_p^{jk},$$

$$q_p^{jk}={\int }_{{\Omega }_k}u{\partial }_x^4(g_p{w}_j)d\Omega .$$

$${\Omega }_k=\{(x,t):|x-x_k|\le H_x,|t-t_k|\le H_t\}$$

$$w_j=({\underline{x}}^2-1{)}^{\alpha }({\underline{t}}^2-1{)}^{\beta }{e}^{\pm il\pi \underline{x}}{e}^{\pm im\pi \underline{t}},$$

$$F_n^{lm}={\int }_{-1}^{1}d\underline{t}{\int }_{-1}^{1}d\underline{x}{f}_n^{\alpha \beta }(\underline{x},\underline{t})e^{\pm il\pi \underline{x}}{e}^{\pm im\pi \underline{t}},$$

$$f_n^{\alpha \beta }(\underline{x},\underline{t})=f_n(u,\underline{x},\underline{t})({\underline{x}}^2-1{)}^{\alpha }({\underline{t}}^2-1{)}^{\beta },$$

$$\Delta c_n=\left|\frac{c_n-{\bar{c}}_n}{{\bar{c}}_n}\right|,$$

$${\epsilon }_n\sim \sigma u_s\sqrt{\frac{{\delta }_x{\delta }_t}{F_x{F}_t}}.$$

$$I={\int }_0^{L}g(x)dx,$$

$${\epsilon }_d\sim \{\begin{array}{ll}{h}^{\mu +2},& \mu \mathrm{even}\\ h^{\mu +1},& \mu \mathrm{odd}\end{array}$$

$${\epsilon }_d\sim h^{\alpha -2}$$

$${\epsilon }_d\sim h^{\alpha }$$

$$Q=\bar{Q}+\epsilon \hat{Q}$$

$$\Vert {\bar{\mathbf{q}}}_n\Vert \propto \{\begin{array}{ll}{F}_x{F}_t,& F_x\ll {\ell }_x,F_t\ll {\ell }_t\\ e^{-{\lambda }_n{F}_x}{\left(\frac{{\ell }_t}{F_t}\right)}^{{\xi }_n}{\ell }_x{\ell }_t,& F_x\gg {\ell }_x,F_t\gg {\ell }_t\end{array}$$

$$\mathbf{c}=\bar{\mathbf{c}}-\epsilon {\bar{Q}}^{+}\hat{Q}\bar{\mathbf{c}},$$

$$\Delta c_n\propto \epsilon \frac{F_x{F}_t{K}^{-1/2}}{\psi (F_x,F_t)},$$

$$\psi (F_x,F_t)=e^{-\lambda F_x}{\left(\frac{{\ell }_t}{F_t}\right)}^{\xi }{\ell }_x{\ell }_t$$

$$\Delta c_n\propto \sigma \sqrt{\frac{{\delta }_x{\delta }_t}{K{F}_x{F}_t}}\frac{F_x{F}_t}{\psi (F_x,F_t)}.$$

$$\Delta c_n\propto \epsilon N_d^{-1/2}$$

$$N_d=\frac{L_x{L}_t}{F_x{F}_t}$$

Full text

data-driven discovery, machine learning, sparse regression, partial differential equations

In recent years, data-driven discovery of mathematical models of spatially extended systems described by nonlinear PDEs has emerged as a promising alternative to more traditional modeling approaches. Existing approaches to model discovery such as sparse regression have several major weaknesses, however. Most notably, they break down for data with high levels of noise and have to be tuned empirically to produce meaningful results, making them ill-suited for analyzing experimental data. We show how these weaknesses can be addressed using a weak formulation of the model PDE. The weak formulation has substantial freedom that makes it extremely powerful and flexible, but the question arises of how this freedom can be used to robustly obtain the most accurate model. This question is addressed here for the first time.

I Introduction

Partial differential equations (PDEs) provide a natural description for the temporal evolution of spatially extended systems in various fields of science and engineering. Historically and practically important examples include wave equations arising in many areas of physics, the Schrödinger equation in quantum mechanics, the Navier-Stokes equations in fluid dynamics, and reaction-diffusion equations used to model physical, chemical, or biological systems. In the past, models of such systems were almost always constructed from first principles or using a suitable empirical approach. However, in recent years, a data-driven paradigm for learning the dynamics has emerged, which leverages the modern prevalence of data and computational power to create models when the underlying governing laws have eluded first-principles derivation.

Many indirect methods for learning the dynamics that do not require a PDE have been proposed. Notable examples include equation-free modeling Kevrekidis et al. (2003), artificial neural networks Hsu et al. (1997); Raissi (2018); Pathak et al. (2018), dynamic mode decomposition Tu et al. (2013) and Koopman operator approaches Mezić (2013), balanced truncation Rowley (2005), and resolvent-based analysis McKeon and Sharma (2010). While these techniques can provide an economical approximate description of the dynamics, this is done at the cost of losing the mathematical structure that affords physical intuition or interpretability. Symbolic regression, which was originally used to derive nonlinear ordinary differential equations describing low-dimensional systems Bongard and Lipson (2007); Schmidt and Lipson (2009), offers an enticing alternative by allowing construction of exact models and discovery of conservation laws. The genetic algorithms used in these earlier studies are however computationally expensive, preventing application of this approach to high-dimensional systems. Thus, the recent emergence of a sparse regression approach for model discovery Xu and Khanmohamadi (2008); Khanmohamadi and Xu (2009); Brunton, Proctor, and Kutz (2016); Rudy et al. (2017); Reinbold and Grigoriev (2019) has made a significant impact. Applied to spatially extended systems, this approach allows data-driven discovery of governing equations in the form of PDEs by evaluating a library of candidate terms containing partial derivatives at a large number of points and using a regularized regression procedure to compute the coefficients of each term and select a parsimonious model.

Sparse regression has proven computationally efficient and capable of reconstructing numerous canonical PDEs Rudy et al. (2017); Li et al. (2019), but it faces serious difficulties when used for analysis of experimental data. One complication is that the proper choice of parsimonious model is often unclear. In many implementations, it relies on a manual Pareto analysis to balance model accuracy and complexity Brunton, Proctor, and Kutz (2016) or on an automatic but complex thresholding procedure (e.g. sequential threshold ridge regression Rudy et al. (2017)) that tends to be sensitive to the choice of parameters. More importantly, existing sparse regression methods often suffer from low accuracy even in the absence of noise and completely break down at noise levels characteristic of realistic applications. This is because they inherently require explicit numerical evaluation of partial derivatives of the data, which is a notably ill-conditioned problem.

In this paper, we present a weak formulation of the sparse regression problem that eliminates this fundamental issue. We also suggest a simple thresholding procedure that can always identify the correct form of the governing PDE even in the presence of extremely high noise. Finally, we explore how this extremely flexible and robust approach can be optimized and tuned to the properties of the underlying data set to maximize accuracy. This paper has the following structure. Section II describes our approach and the system used to test it. Results are presented and interpreted in Section III, and conclusions are discussed in Section IV.

II Methods

We consider the problem of using the data ${\bf u}({\bf x},t)$ to identify a parsimonious mathematical model in the form of a PDE

[TABLE]

where each term in the sum is a function of ${\bf u}$ and its partial derivatives in space and time with constant coefficients $c_{n}$. In most applications, the form of the basis functions $\mathbf{f}_{n}$ can be restricted based on physical considerations, such as symmetries, conservation laws, etc. Bär, Hegger, and Kantz (1999); Reinbold and Grigoriev (2019). Typically, $\mathbf{f}_{n}$ are taken to be products of powers of independent variables (${\bf x}$, $t$) and dependent variables (${\bf u}$ and its various derivatives), although the form can be arbitrary in theory. Our goal is to determine the constants $c_{n}$ for the terms that should be present in the model while eliminating the dynamically insignificant and thus likely spurious terms. Sparse regression aims to convert the PDE (1) to a tractable (and ideally, robust) linear algebra problem. Conventionally this is done by evaluating all of the terms in the PDE at a random collection of points $({\bf x}_{k},t_{k})$ using finite differences Vallette, Jacobs, and Gollub (1997); Bär, Hegger, and Kantz (1999), spectral methods Xu and Khanmohamadi (2008); Khanmohamadi and Xu (2009), or polynomial approximation Rudy et al. (2017); Reinbold and Grigoriev (2019). All of these approaches are extremely sensitive to noise, especially when high-order derivatives are present. We will instead pursue a weak formulation of the problem that can be obtained by multiplying (1) by a weight $\mathbf{w}_{j}(\mathbf{x},t)$ and then integrating the result over a domain $\Omega_{k}$. Repeating the process for $K$ distinct combinations of weight functions and integration domains yields the linear system

[TABLE]

where $\mathbf{c}=[c_{1},\dots,c_{N}]^{T}$ and $Q=[\mathbf{q}_{1},\dots,\mathbf{q}_{N}]$ is a “library” matrix, with each column $\mathbf{q}_{n}\in\mathbb{R}^{K}$ consisting of the integrals of the function $\mathbf{f}_{n}$ with all $K$ combinations of weights $\mathbf{w}_{j}$ and domains $\Omega_{k}$.

Note that there is an extra degree of freedom in (2) corresponding to the normalization of $\mathbf{c}$. Conventionally this is dealt with by assuming that ${\bf f}_{1}=\partial_{t}{\bf u}$, setting $c_{1}=1$, and solving the overdetermined system that corresponds to the choice $K\gg N$ using least squares or some regularized version of it Rudy et al. (2017). This is however not always a valid assumption: it is usually unknown a priori whether any given temporal derivative should be included in the PDE at all, whereas in this case a particular term is forced into the model. Moreover, even if this term should be present in the model, the regression effectively assumes that the time derivative was computed without error, which reduces the practical accuracy of the procedure.

We will therefore not make the assumption that the model has the form of an evolution equation and consider the linear problem (2) in its most general form. The normalization of $\mathbf{c}$ can be fixed by adding an extra row with arbitrary nonzero elements to $Q$, after which the resulting equation (2) can be solved by ordinary least squares. A more elegant solution pursued in the present study is to instead compute $\mathbf{c}$ as the right singular vector of $Q$ corresponding to the smallest singular value. Note that this corresponds to the solution of a constrained least squares problem for $Q^{T}Q\mathbf{c}=0$:

[TABLE]

Once a suitable solution has been obtained by further constraining the problem, the resulting parsimonious model can be rewritten in the form of an evolution equation by solving for a term such as $\partial_{t}\mathbf{u}$ (or $\partial_{t}^{2}\mathbf{u}$ for a wave equation).

To obtain a parsimonious model, we employ an iterative procedure to eliminate unnecessary terms from (1). At each step $i$, singular value decomposition is used to obtain the solution ${\bf c}^{i}$ given the matrix $Q^{i}$, and the residual $\eta^{i}=\|Q^{i}{\bf c}^{i}\|$ is computed. We then find the term with the smallest $\|c^{i}_{n}\mathbf{q}_{n}\|/\|\mathbf{q}_{n}\|$ and construct $Q^{i+1}$ by eliminating the column ${\bf q}_{n}$ from $Q^{i}$. The corresponding term is eliminated from the model if $\eta^{i+1}<\gamma\eta^{i}$, where $\gamma>1$ is some fixed constant (we use $\gamma=1.4$ in the present study). The iteration terminates at step $i$ if $\eta^{i+1}>\gamma\eta^{i}$, yielding a parsimonious model. We find that this method compares favorably to alternatives such as sequential threshold ridge regression Rudy et al. (2017) as it robustly eliminates spurious terms without requiring extremely careful choice of parameters. Moreover, the sparsification parameter has a simple interpretation: $\gamma-1$ is the maximum acceptable relative increase in the residual resulting from discarding a single library term.

We illustrate the advantages of our approach by applying it to the Kuramoto-Sivashinsky equation Kuramoto (1978); Sivashinsky (1985)

[TABLE]

which has posed a significant challenge in past studies of sparse regression Xu and Khanmohamadi (2008); Rudy et al. (2017) because it contains a fourth-order partial derivative that is difficult to evaluate numerically with adequate accuracy. Here $c_{1}=\cdots=c_{4}=1$ are all constants, although our approach can easily be extended even to the case when these coefficients are functions of time and/or space, as discussed below. Since this is a scalar equation in one spatial and one temporal dimension, we use scalar weight functions $w^{j}(x,t)$. If we denote the terms in the model (4) by $f_{1},\cdots,f_{4}$, then

[TABLE]

where $d\Omega=dx\,dt$. The key feature of the weak formulation is that it can almost always be used to completely eliminate, or at least reduce the order of, the derivatives acting on the noisy data by integrating by parts. In our particular case,

[TABLE]

under the assumption that $w_{j}$ and its first three partial derivatives with respect to $x$ vanish on the boundary $\partial\Omega_{k}$. In our implementation, we use the composite trapezoidal rule to evaluate the integrals numerically.

Note that although this particular PDE features constant coefficients, terms with variable coefficients can be treated in a similar manner. For instance, suppose that the coefficient of the term $\partial_{x}^{4}u$ is a function of ${\bf x}$ and $t$ that can be expanded in some (finite) basis as

[TABLE]

with some constants $c^{\prime}_{p}$. Then

[TABLE]

where

[TABLE]

Sparse regression for a model including such a term would then simply require expanding the library $Q$ to include additional columns ${\bf q}_{p}$ with entries $q_{p}^{jk}$. In this case as well, no derivatives of the noisy $u$ are used in finding the elements of $Q$.

Although in principle integration domains of any shape can be used, here we will only consider rectangular domains of a fixed size

[TABLE]

where the centers $(x_{k},t_{k})$ of the rectangles $\Omega_{k}$ are chosen randomly. Similarly, there are many possible choices for the weight functions satisfying the boundary conditions on $\partial\Omega_{k}$; we focus on functions of the form

[TABLE]

where $\underline{x}=(x-x_{k})/H_{x}$, $\underline{t}=(t-t_{k})/H_{t}$ are nondimensionalized independent variables and $\alpha\geq 4$, $\beta\geq 1$, $l\geq 0$, and $m\geq 0$ are integers. Note that there are four weight functions (corresponding to the four different choices of the signs in the exponentials) for each pair of nonzero $l$ and $m$. The integrals $q_{n}^{jk}$ are all of the form

[TABLE]

where

[TABLE]

so $F_{n}^{lm}$ are the coefficients of the two-dimensional Fourier series for $f_{n}^{\alpha\beta}(\underline{x},\underline{t})$. Although $f_{n}(u,\underline{x},\underline{t})$ is not periodic on $\Omega_{k}$, the functions $f_{n}^{\alpha\beta}(\underline{x},\underline{t})$ are. Moreover, $f_{n}^{\alpha\beta}(\underline{x},\underline{t})$ has at least $\alpha-1$ continuous derivatives in $\underline{x}$ and ${\beta-1}$ continuous derivatives in $\underline{t}$, so the Fourier coefficients decay according to $F_{n}^{lm}\sim l^{-\alpha}m^{-\beta}$. The powers $\alpha$ and $\beta$ therefore control the width of the Fourier spectrum of the entries $q_{n}^{jk}$ in the library $Q$, while the choice of $l$ and $m$ allows us to tune the frequencies of the weights to the spectral properties of the data. The convergence rate of Fourier series turns out to control the accuracy with which the integrals are evaluated using data that are available only on a discrete grid. For simplicity, we will assume that the same weight functions are integrated on every domain. It is possible to use either weight functions involving only a single pair of frequencies (e.g., $l$ and $m$) or a range of frequencies in space and/or time.

To test our sparse regression approach, we computed a solution of the Kuramoto-Sivashinsky equation, using the integrator described in Ref. Rudy et al., 2017 to generate data on a physical domain with dimensions $L_{x}=32\pi$ and $L_{t}=500$. The numerical integration generated data with spatial resolution $\Delta x=0.0491$ using a time step $\Delta t=0.005$, which was then downsampled to a lower spatial resolution $\delta_{x}$ and temporal resolution $\delta_{t}$. Unless noted otherwise, the results presented below are for $\delta_{x}=0.1964$ and $\delta_{t}=1$. For reference, the solution has a correlation length $\ell_{x}\sim 1.67\approx 8.5\delta_{x}$ and correlation time $\ell_{t}\sim 8=8\delta_{t}$. To test the effects of noise, Gaussian noise with standard deviation $\sigma s_{u}$ was added to the data for various choices of $\sigma$, where $s_{u}\approx 1.3$ is the sample standard deviation of $u$ on the whole domain.

To test the ability of the algorithm to eliminate spurious terms, in addition to the terms present in the Kuramoto-Sivashinsky equation (4), we also included terms $\partial_{x}u$, $\partial_{x}^{3}u$, $u$, $u^{2}$, $u^{3}$, and $1$ (which represents a hypothetical forcing) in our library. The corresponding integrals were rewritten using integration by parts to remove derivatives acting on $u$, as described previously. In the next section, we quantify the performance of our sparse regression approach using two key metrics: how well the algorithm can discriminate between the essential and spurious terms and how accurately it can determine the coefficients of the essential terms. Since the data were generated using a known model, we know which terms are essential (those contained in the PDE (4)). If the reference model is unavailable, ensemble regression Reinbold and Grigoriev (2019) may be used instead to help distinguish essential terms from spurious ones.

III Results

As discussed previously, the elements of the library matrix $Q$ are given by the Fourier coefficients of the different terms included in the generalized model (windowed by the envelope $E^{\alpha\beta}(\underline{x},\underline{t})=(\underline{x}^{2}-1)^{\alpha}(\underline{t}^{2}-1)^{\beta}$ on each domain $\Omega_{k}$); hence knowledge of the Fourier spectrum of the data is crucial for an optimal choice of the size of the integration domains $\Omega_{k}$ and the weight functions $w_{j}$. The power spectrum (or, more precisely, the absolute value of the Fourier coefficients) of the noiseless data on the entire physical domain is shown in Figure 1. In space, the spectrum is sharply peaked around a wave number $\kappa\approx 0.625$. At high wave numbers, the spectrum decays exponentially, $P\propto e^{-\kappa/\bar{\kappa}}$ where $\bar{\kappa}\approx 0.3$. In time, the spectrum is peaked at zero frequency $\omega$ and decays as a power law, $P\propto\omega^{-\chi}$ with $\chi\approx 2.5$.

Having characterized the data, we turn to the investigation of how the performance of our algorithm depends on the choice of various parameters. Since the number of parameters is quite large, instead of exploring the entire parameter space, we focus on the dependence on one or two parameters at a time, with the remaining parameters staying fixed. Specifically, the noise level $\sigma$ is fixed to 3% and we use the following near-optimal parameters in the sparse regression. The dimensions of the integration domain are $F_{x}=2H_{x}=14.73$ and $F_{t}=2H_{t}=75$. This choice corresponds to an equal number of grid points in both directions, $F_{x}/\delta^{*}_{x}=F_{t}/\delta^{*}_{t}=75$. Unless noted otherwise, we use a single set of weights with $\alpha=\beta=8$, $l=1$, $m=2$, and the sparsification parameter is $\gamma=1.4$. We generally use every combination of 4 weight functions over 50 integration domains, so that the total number of library rows is $K=200$. To characterize the stochastic effects, for each set of parameters, we used an ensemble of $M=100$ trials featuring different random distributions of the integration domains and realizations of noise.

First, we tested the ability of the method to reconstruct the correct form of the PDE (4) for various values of $\gamma$ with all other parameters fixed at their near-optimal values. Our iterative regression procedure proved very robust for a fairly wide range of values of $\gamma$. In particular, at a noise level of 30%, it performed perfectly for $1.1\leq\gamma\leq 2$, with the reconstructed model containing no missing or spurious terms in all of the trials. For the highest noise level considered here (100%), we found perfect performance for $1.2\leq\gamma\leq 1.5$. In some fraction of the trials, spurious terms appeared at lower $\gamma$ and missing terms at higher $\gamma$, as shown in Figure 2. For reference, without the benefit of the weak formulation, sparse regression failed Rudy et al. (2017) to correctly reconstruct the lambda-omega reaction-diffusion system, which is only second-order, for noise level as low as 1%.

The accuracy of regression (i.e., model identification) was quantified by computing the relative error in each parameter of the Kuramoto-Sivashinsky equation

[TABLE]

where $\bar{c}_{n}$ and $c_{n}$ are the true and estimated values of the model parameters, respectively, for $n=2,3,4$. We normalize the estimated parameters so that $c_{1}=1$. In all of the following figures, we plot the estimated mean value of $\Delta c_{n}$ with 95% error bars, where all of the parameters in the regression procedure are held at their near-optimal values stated previously unless noted otherwise.

In particular, Fig. 3 shows the accuracy of regression as a function of $\sigma$ for two different choices of data resolution ($\delta_{x}$ and $\delta_{t}$). It is worth noting that the average relative error is $\Delta c_{n}\sim 10^{-10}$ for all of the parameters for noiseless data with the higher of the two resolutions. However, even for $1\%$ noise, $\Delta c_{n}\sim 2\times 10^{-4}$, which is more than three orders of magnitude smaller than what had been achieved in previous studies Rudy et al. (2017). The results are very similar for all three parameters; as this is generally the case, in subsequent figures, we only show the generally largest error $\Delta c_{4}$, which corresponds to the term $\partial_{x}^{4}u$ involving the highest-order derivative. We find two distinct regimes. At higher noise levels, the error in evaluating the library matrix entries is due primarily to the averaged effect of noise. Applying the central limit theorem, we find that the relative error scales as

[TABLE]

At low noise levels, the parameter accuracy is controlled by numerical error, which has two different sources. The first source is a numerical error in the data itself, which is due to the finite accuracy of the integrator that “solves” the Kuramoto-Sivashinsky equation. This source dominates for smaller $\delta_{x}$ and $\delta_{t}$. For experimental data, this source would correspond to systematic error. For larger $\delta_{x}$ and $\delta_{t}$, the parameter inaccuracy is mainly due to the error in computing the library matrix entries based on data that are available on a discrete grid. Suppose we want to use numerical quadratures to evaluate an integral

[TABLE]

where $g(x)\in C^{m}$ (i.e., has $m$ continuous derivatives) and $g^{(i)}(0)=g^{(i)}(L)$ for all $0\leq i<m$. Then, for the composite trapezoidal rule on a grid with spacing $h$, the relative error associated with the discretization can be estimated using exact Euler-Maclaurin formulasTrefethen and Weideman (2014) and is found to scale as $h^{m+2}|g^{(m+2)}|$ for $m$ even (or $h^{m+1}|g^{(m+1)}|$ for $m$ odd), where a characteristic value of the derivative on the interval $[0,L]$ is used.

Generalizing this result to two dimensions (and assuming $\delta_{x}\ll\min(\ell_{x},F_{x})$, $\delta_{t}\ll\min(\ell_{t},F_{t})$), we find an estimate of the relative discretization error for an element of the library matrix $Q$ that involves a temporal derivative of order $\nu_{t}$ and/or spatial derivative of order $\nu_{x}$:

[TABLE]

where $h=\delta_{t}/\ell_{t}\approx\delta_{x}/\ell_{x}$ and $\mu=\min(\alpha-\nu_{x},\beta-\nu_{t})$. It is easy to check that, due to the conditions on $\alpha$ and $\beta$, we always have $\mu\geq 0$, as it should be for the trapezoidal rule. The Kuramoto-Sivashinsky equation features terms that all involve derivatives, with the lowest order being one and the highest being four; hence, for even $\alpha=\beta\geq 4$, the exponent $\mu$ ranges between $\alpha-2$ and $\alpha+1$. Therefore the scaling

[TABLE]

dominates for lower $h$, while the scaling

[TABLE]

dominates for higher $h$.

The error $\Delta c_{n}$ can be found using perturbation theory. Let $\bar{Q}$ be the library matrix evaluated using a continuous noiseless solution so that $\bar{Q}\bar{\mathbf{c}}=0$ exactly (we assume that $\bar{Q}$ corresponds to the parsimonious model). In the presence of measurement noise and/or discretization error, the error in evaluating each entry $q_{n}^{jk}$ of the library matrix is proportional to $\varepsilon=\max(\varepsilon_{d},\varepsilon_{n})$, so

[TABLE]

for some matrix $\hat{Q}$ whose entries are distributed as white Gaussian noise. Note that the entries of $\hat{Q}$ are $O(F_{x}F_{t})$. The entries of $\bar{Q}$ have a more complicated scaling that is determined by the Fourier spectrum of the data (i.e., exponential in space, power law in time). Specifically, we find (cf. Fig. 4)

[TABLE]

where $\lambda_{n}=O(\ell_{x}^{-1})$ and $\xi_{n}=O(1)$ are some positive constants. To leading order in $\varepsilon$, the least squares solution to (2) is given by

[TABLE]

where $\bar{Q}^{+}$ is the Moore-Penrose pseudoinverse of $\bar{Q}$. Since the elements of $\hat{Q}$ can be considered uncorrelated, we have for $F_{x}\gg\ell_{x}$ and $F_{t}\gg\ell_{t}$

[TABLE]

where the numerator and denominator describe the scaling of the entries of $\hat{Q}$ and $\bar{Q}$, respectively. Following from (21),

[TABLE]

with some positive constants $\lambda=O(\ell_{x}^{-1})$ and $\xi=O(1)$. For low $\sigma$, we have $\varepsilon=\varepsilon_{d}$ and therefore $\Delta c_{n}$ is independent of $\sigma$. For high $\sigma$, we have $\varepsilon=\varepsilon_{n}$, so combining (23) and (15) we find

[TABLE]

The predicted scaling of $\Delta c_{n}$ with $\sigma$ in both regimes is consistent with the results shown in Fig. 3. In particular, we find that the effect of changing the resolution of the data is quite minor at high $\sigma$, where $\Delta c_{n}\propto h$ according to (25). At low $\sigma$, the effect is much stronger: for $\alpha=\beta=8$, we have $\Delta c_{n}\propto h^{6}$ according to (17). The dependence of the scaling in (17) on $\alpha$ and $\beta$ is further confirmed by Fig. 5, which shows results for noiseless data. In the $\alpha=\beta=4$ case, we observe the scaling law $\Delta c_{n}\propto h^{2}$ corresponding to (18) in the entire range of $h$ we examined. When $\alpha=\beta=6$, the parameter error scales according to $\Delta c_{n}\propto h^{4}$ for small $h$ and $\Delta c_{n}\propto h^{6}$ for large $h$, which correspond to the limiting cases (18) and (19), respectively. We should also note that for $h$ as large as $1/4$, the accuracy remains very good. Thus, the method is suitable for fairly sparse data.

As illustrated in Fig. 6, we also observe the scaling for $\Delta c_{n}$ with $K$ predicted by (25). This scaling is expected to break down when the total area of the integration domains exceeds the area of the physical domain due to the loss of statistical independence between the data on different integration domains, leading to an increased linear dependence of the rows of the library matrix $Q$. We can expect the error to asymptote to

[TABLE]

for $K\gg N_{d}$, where

[TABLE]

is the area ratio. For the reference set of parameters, saturation did not occur over the range of $K$ we tested. To more easily observe the saturation effect, we set $l=m=0$, so that only one weight function is used and the number of integration domains equals $K$ (rather than $K/4$ for nonzero $l$ and $m$). Furthermore, we reduce the size of the physical domain to $L_{x}=16\pi$ and $L_{t}=250$, so that $N_{d}\approx 11$ is relatively small. As Fig. 6 illustrates, for large $K$, the parameter accuracy indeed asymptotes to a constant.

The scaling described by (26) can also be observed in the dependence of $\Delta c_{n}$ on the size of the physical domain (and hence $N_{d}$) with all other parameters fixed. This dependence is quite important, since it determines how much data needs to be collected to identify the model with meaningful precision. As Fig. 7 illustrates, choosing the physical domain to be just double the size of the (optimal) integration domain in both directions (which corresponds to $N_{d}=4$) already yields a rather acceptable accuracy when only one weight function is used. When $l$ and $m$ are nonzero, accurate reconstruction is possible even if $N_{d}$ is only slightly greater than $1$.

Next, we consider how the error in the estimated coefficients depends on the choice of the integration domain size. Figure 8 shows the dependence of the error $\Delta c_{4}$ on the size of the integration domains for two different choices of the weight functions. In panel (a), $F_{t}$ is fixed to 75 and $F_{x}$ is taken to vary, and in panel (b), $F_{x}$ is fixed to 14.73 with $F_{t}$ varying. In both cases, we find that there is an optimal domain size with $F_{x}\approx 14.73$ and $F_{t}\approx 75$; moreover, the optimal values remain approximately the same even if we vary the size of the other dimension or the choice of weight functions. For small $F_{x}$ and/or $F_{t}$, the error is large because (a) the integration domain is too small to effectively average out the influence of noise and (b) the numerical quadrature error is large (both $\epsilon_{n}$ and $\epsilon_{d}$ increase as $F_{x}$ and/or $F_{t}$ decrease). For large $F_{x}$ and $F_{t}$, we enter the regime described by (23), which predicts that the error should grow exponentially in $F_{x}$ and as a power of $F_{t}$. Indeed, this is exactly what we observe in Fig. 8. Based on (21), it appears that the optimal choice of $F_{x}$ and $F_{t}$ corresponds to the crossover between these two regimes, i.e., $F_{x}\propto\ell_{x}$ and $F_{t}\propto\ell_{t}$. Our numerical results suggest that the optimal choice corresponds to $F_{x}/\ell_{x}\approx F_{t}/\ell_{t}\approx 8$.

Finally, let us address the optimal choice of frequencies appearing in the weight functions (11). Figure 9 shows the effect of varying either $l$ or $m$ with all other parameters fixed at their reference values. Specifically, we plot $\Delta c_{4}$ versus $\kappa_{l}=2\pi l/F_{x}$ and $\omega_{m}=2\pi m/F_{t}$. (Note that when $l$ or $m$ is 0, the number of distinct weight functions is halved, so we correspondingly double the number of integration domains to keep the number of rows in the library constant.) One could assume that the optimal values would be given by the dominant frequencies of the original data (as we discussed previously, the dominant wave number is $\kappa\approx 0.625$ and the dominant temporal frequency is $\omega=0$). According to Fig. 1, windowing the data broadens the peaks but leaves both dominant frequencies roughly the same: $\kappa\approx 0.8$ ($l=2$) and $\omega=0$ ($m=0$). Unfortunately, it turns out that we cannot use the spectra to exactly predict the optimal frequencies, which are $\kappa\approx 0.4$ ($l=1$) and $\omega\approx 0.2$ ($m=2$ or $3$). However, choosing the frequencies based on the spectra still produces reasonably good accuracy (within a factor of $4$ or so of the optimal result).

These results suggest that using weight functions with a combination of different frequencies may be more robust and/or accurate. To test this hypothesis, we considered the case in which weight functions with a range of frequencies in space or time were included, with the total number of library rows fixed at $200$. However, this approach yielded a decrease in the accuracy, as the broader choice of weight functions did not compensate for a decrease in the number of integration domains. This suggests that the optimal strategy is to use a large number of integration domains while keeping the frequencies of the weight functions fixed.

IV Conclusions

We have introduced a robust and flexible approach to data-driven discovery of models in the form of nonlinear PDEs. The approach uses a weak formulation, coupled with a novel sparse regression procedure, to obtain a parsimonious description. We have demonstrated its capability to identify PDEs, even with high-order derivatives, from extremely noisy data with unprecedented accuracy. For instance, with 1% noise, we were able to reduce the error in estimating the parameters of the 4th-order Kuramoto-Sivashinsky equation from 50% Rudy et al. (2017) to just $2\times 10^{-4}$. Furthermore, whereas correct identification of the functional form of the underlying PDE has been far from guaranteed at any noise level using past approaches, our algorithm was able to reconstruct the Kuramoto-Sivashinsky equation accurately in $100\%$ of cases from data with a signal-to-noise ratio of 100%.

This impressive performance is achieved by shifting the partial derivatives from the data onto a known smooth weight function using integration by parts, thus avoiding the large errors incurred by repeated numerical differentiation. Our method also proved to be well-adapted to sparse data, maintaining errors of less than $0.1\%$ for a grid resolution only $4$ times finer than the correlation length/time. Such reliability and high accuracy in the presence of noisy or sparse data is indispensable for analysis of experimental data. Notably, even in the absence of noise, our results compare very favorably with those of previous studies Xu and Khanmohamadi (2008); Rudy et al. (2017) because the discretization error of the algorithm can be made extremely small: for the Kuramoto-Sivashinsky equation, the relative error in all parameters can easily be reduced to $10^{-10}$. It is also important to mention that the computational cost of our algorithm is comparable to that of existing sparse regression methods.

We also derived the scaling laws that describe the accuracy of the regression as a function of the parameters used in the algorithm and the properties of the data. These scaling laws can be used to fully exploit the flexibility of the weak formulation approach by tuning its various paramters. In particular, the size of the input used by the regression can be controlled by choosing both the number of different integration domains and the number of different weight functions. We have shown that the number of integration domains plays a much more important role than the number of weight functions: the best results can be obtained by using a set of weight functions with a fixed shape (frequency and envelope) and a large number of integration domains. Furthermore, we have determined the optimal shape of the weights and the optimal size of the integration domains. The latter turned out to be determined by the correlation length and time describing the data (with the size roughly an order of magnitude larger than these characteristic scales). We have also shown that, although the error can be reduced further by using data on ever-larger physical domains, satisfactory results can be obtained for physical domains that are just a factor of two larger than the optimal integration domain in each dimension.

Acknowledgements.

This material is based upon work supported by the National Science Foundation under Grant No. CMMI-1725587. DG gratefully acknowledges the support of the Letson Undergraduate Research Scholarship.

V References