Hidden Physics Models: Machine Learning of Nonlinear Partial Differential Equations

Introduction

There are more than a trillion sensors in the world today and according to some estimates there will be about 50 trillion cameras worldwide within the next five years, all collecting data either sporadically or around the clock. However, in scientific experiments, quality and error-free data is not easy to obtain – e.g., for system dynamics characterized by bifurcations and instabilities, hysteresis, and often irreversible responses. Admittedly, as in all everyday applications, in scientific experiments too, the volume of data has increased substantially compared to even a decade ago but analyzing big data is expensive and time-consuming. Data-driven methods, which have been enabled in the past decade by the availability of sensors, data storage, and computational resources, are taking center stage across many disciplines of science. We now have highly scalable solutions for problems in object detection and recognition, machine translation, text-to-speech conversion, recommender systems, and information retrieval. All of these solutions attain state-of-the-art performance when trained with large amounts of data. However, purely data driven approaches for machine learning present difficulties when the data is scarce relative to the complexity of the system. Hence, the ability to learn in a sample-efficient manner is a necessity in these data-limited domains. Less well understood is how to leverage the underlying physical laws and/or governing equations to extract patterns from small data generated from highly complex systems. In this work, we propose a modeling framework that enables blending conservation laws, physical principles, and/or phenomenological behaviors expressed by partial differential equations with the datasets available in many fields of engineering, science, and technology. This paper should be considered a direct continuation of a preceding one in which we addressed the problem of inferring solutions of time dependent and nonlinear partial differential equations using noisy observations. Here, a similar methodology is employed to deal with the problem of learning, system identification, or data-driven discovery of partial differential equations .

Problem Setup

Let us consider parametrized and nonlinear partial differential equations of the general form

Here, $h^{n}(x)=h(t^{n},x)$ is the hidden state of the system at time $t^{n}$. Approximating the nonlinear operator on the left-hand-side of equation (2) by a linear one we obtain

involved in the Burgers’ equation can be approximated by the linear operator

where $h^{n-1}(x)$ is the state of the system at the previous time $t^{n-1}$.

The Basic Model

Similar to Raissi et al. , we build upon the analytical property of Gaussian processes that the output of a linear system whose input is Gaussian distributed is again Gaussian. Specifically, we proceed by placing a Gaussian processGaussian processes (see ) provide a flexible prior distribution over functions and enjoy analytical tractability. They can be viewed as a prior on one-layer feed-forward Bayesian neural networks with an infinite number of hidden units . Gaussian processes are among a class of methods known as kernel machines (see ) and are analogous to regularization approaches (see ). prior over the latent function $h^{n}(x)$; i.e.,

Here, $\theta$ denotes the hyper-parameters of the covariance function $k$. Without loss of generality, all Gaussian process priors used in this work are assumed to have a squared exponentialFrom a theoretical point of view, each kernel (i.e., covariance function) gives rise to a Reproducing Kernel Hilbert Space (RKHS) that defines a class of functions that can be represented by this kernel. In particular, the squared exponential covariance function implies smooth approximations. For a more systematic treatment of the kernel-selection problem we would like to refer the readers to . Furthermore, more complex function classes can be accommodated by employing nonlinear warping of the input space to capture discontinuities . covariance function, i.e.,

where $\theta=\left(\gamma,w_{1},\cdots,w_{D}\right)$ are the hyper-parameters and $x$ is a $D$-dimensional vector. The Gaussian process prior assumption (4) along with equation (3) enable us to capture the entire structure of the operator $\mathcal{L}_{x}^{\lambda}$ in the resulting multi-output Gaussian process

It is worth highlighting that the parameters $\lambda$ of the operators $\mathcal{L}^{\lambda}_{x}$ and $\mathcal{N}^{\lambda}_{x}$ turn into hyper-parameters of the resulting covariance functions. The specific forms of the kernelsIt should be noted that for all examples studied in this work the kernels are generated at the push of a button using Wolfram Mathematica, a mathematical symbolic computation program.

are direct functions of equation (3) as well as the prior assumption (4); i.e.,

We call the multi-output Gaussian process (5) a hidden physics model, because its matrix of covariance functions explicitly encodes the underlying laws of physics expressed by equations (1) and (3).

Learning

Given the noisy data $\{\bm{x}^{n-1},\bm{h}^{n-1}\}$ and $\{\bm{x}^{n},\bm{h}^{n}\}$ on the latent solution at times $t^{n-1}$ and $t^{n}$, respectively, the hyper-parameters $\theta$ of the covariance functions and more importantly the parameters $\lambda$ of the operators $\mathcal{L}_{x}^{\lambda}$ and $\mathcal{N}_{x}^{\lambda}$ can be learned by employing a Quasi-Newton optimizer L-BFGS to minimize the negative log marginal likelihood

where $\bm{h}=\begin{bmatrix}\bm{h}^{n}\\ \bm{h}^{n-1}\end{bmatrix}$, $p(\bm{h}|\theta,\lambda,\sigma^{2})=\mathcal{N}\left(\bm{0},\bm{K}\right)$, and $\bm{K}$ is given by

Here, $N$ is the total number of data points in $\bm{h}$. Moreover, $\sigma^{2}$ is included to capture the noise in the data and is also learned by minimizing the negative log marginal likelihood. The implicit underlying assumption is that $\bm{h}^{n}=h^{n}(\bm{x}^{n})+\bm{\epsilon}^{n}$ and $\bm{h}^{n-1}=h^{n-1}(\bm{x}^{n-1})+\bm{\epsilon}^{n-1}$ with $\bm{\epsilon}^{n}\sim\mathcal{N}(0,\sigma^{2}I)$ and $\bm{\epsilon}^{n-1}\sim\mathcal{N}(0,\sigma^{2}I)$ being independent. The negative log marginal likelihood (6) does not simply favor the models that fit the training data best. In fact, it induces an automatic trade-off between data-fit and model complexity. Specifically, minimizing the term $\bm{h}^{T}\bm{K}^{-1}\bm{h}$ in equation (6) targets fitting the training data, while the log-determinant term $\log|\bm{K}|$ penalizes model complexity. This regularization mechanism automatically meets the Occam’s razor principle which encourages simplicity in explanations. The aforementioned regularization mechanism of the negative log marginal likelihood (6) effectively guards against overfitting and enables learning the unknown model parameters from very fewRegularization is important even in data abundant regimes as witnessed by the recently growing literature on discovering ordinary and partial differential equations from data using sparse regression techniques . noisy observations. However, there is no theoretical guarantee that the negative log marginal likelihood does not suffer from multiple local minima. Our practical experience so far with the negative log marginal likelihood seems to indicate that local minima are not a devastating problem, but certainly they do exist. Moreover, it should be highlighted that, although not pursued here, a fully Bayesian and more robust estimate of the linear operator parameters $\lambda$ can be obtained by assigning priors on $\{\theta,\lambda,\sigma^{2}\}$. However, this would require more costly sampling procedures such as Markov Chain Monte Carlo (see , chapter 5) to train the model. Furthermore, the most computationally intensive part of learning using the negative log marginal likelihood (6) is associated with inverting dense covariance matrices $\bm{K}$. This scales cubically with the number $N$ of training data in $\bm{h}$. While it has been effectively addressed by the recent works of , this cubic scaling is still a well-known limitation of Gaussian process regression.

Results

The proposed framework provides a general treatment of time-dependent and nonlinear partial differential equations, which can be of fundamentally different nature. This generality will be demonstrated by applying the algorithm to a dataset originally proposed in , where sparse regression techniques are used to discover partial differential equations from time series measurements in the spatial domain. This dataset covers a wide range of canonical problems spanning a number of scientific domains including the Navier-Stokes, Schrödinger, and Kuramoto-Sivashinsky equations. Moreover, all data and codes used in this manuscript are publicly available on GitHub at https://github.com/maziarraissi/HPM.

Burgers’ equation arises in various areas of applied mathematics, including fluid mechanics, nonlinear acoustics, gas dynamics, and traffic flow . It is a fundamental partial differential equation and can be derived from the Navier-Stokes equations for the velocity field by dropping the pressure gradient term. Burgers’ equation, despite its relation to the much more complicated Navier-Stokes equations, does not exhibit turbulent behavior. However, for small values of the viscosity parameters, Burgers’ equation can lead to shock formation that is notoriously hard to resolve by classical numerical methods. In one space dimension the equation reads as

with $(\lambda_{1},\lambda_{2})$ being the unknown parameters. The original data-set proposed in contains 101 time snapshots of a solution to the Burgers’ equation with a Gaussian initial condition, propagating into a traveling wave. The snapshots are $\Delta t=0.1$ apart. The spatial discretization of each snapshot involves a uniform grid with 256 cells. As depicted in figure 1 using only two of these snapshots (randomly selected) with 71 and 69 data points, respectively, the algorithm is capable of identifying the correct parameter values up to a relatively good accuracy. It should be noted that we are using only $140=71+69$ data points out of a total of $25856=101\times 256$ in the original data set. This surprising performance is achieved at the cost of explicitly encoding the underlying physical laws expressed by the Burgers’ equation in the covariance functions of the hidden physics model (5). For a systematic study of the performance of the method, let us carry out the same experiment as the one illustrated in figure 1 for every pair of consecutive snapshots in the original dataset. We are still using the same number of data points (i.e., 71 and 69) for each pair of snapshots, albeit in different locations. The resulting statistics for the learned parameter values are reported in table 1. As is clearly demonstrated in this table, more noise in the data leads to less confidence in the estimated values for the parameters. Moreover, let us recall the main assumption of this work that the gap $\Delta t$ between the pair of snapshots should be small enough so that we can employ the backward Euler scheme (see equation (2)). To test the importance of this assumption, let us use the exact same setup as the one explained in figure 1, but increase $\Delta t$. The results are reported in table 2. Therefore, the most important facts about the proposed methodology are that more data, less noise, and a smaller gap $\Delta t$ between the two snapshots enhance the performance of the algorithm.

2 The KdV Equation

As a mathematical model of waves on shallow water surfaces one could consider the Korteweg-de Vries (KdV) equation. This equation can also be viewed as Burgers’ equation with an added dispersive term. The KdV equation has several connections to physical problems. It describes the evolution of long one-dimensional waves in many physical settings. Such physical settings include shallow-water waves with weakly non-linear restoring forces, long internal waves in a density-stratified ocean, ion acoustic waves in a plasma, and acoustic waves on a crystal lattice. Moreover, the KdV equation is the governing equation of the string in the Fermi-Pasta-Ulam problem in the continuum limit. The KdV equation reads as

with $(\lambda_{1},\lambda_{2})$ being the unknown parameters. The original dataset proposed in contains a two soliton solution to the KdV equation with 512 spatial points and 201 time-steps. The snapshots are $\Delta t=0.1$ apart. As depicted in figure 2 using only two of these snapshots (randomly selected) with 111 and 109 data points, respectively, the algorithm is capable of identifying the correct parameter values up to a relatively good accuracy. In particular, we are using $220=111+109$ out of a total of $102912=201\times 512$ data points in the original data set. This level of efficiency is a direct consequence of equation (5) where the covariance functions explicitly encode the underlying physical laws expressed by the KdV equation. As a sensitivity analysis of the reported results, let us perform the same experiment as the one illustrated in figure 2 for every pair of consecutive snapshots in the original dataset. We are still using the same number of data points (i.e., 111 and 109) for each pair of snapshots, albeit in different locations. The resulting statistics for the learned parameter values are reported in table 3. As is clearly demonstrated in this table, more noise in the data leads to less confidence in the estimated values for the parameters. Moreover, to test the sensitivity of the results with respect to the gap between the two time snapshots, let us use the exact same setup as the one explained in figure 2, but increase $\Delta t$. The results are reported in table 4. These results verify the most important facts about the proposed methodology that more data, less noise, and a smaller gap $\Delta t$ between the two snapshots enhance the performance of the algorithm.

3 Kuramoto-Sivashinsky Equation

The Kuramoto-Sivashinsky equation has similarities with Burgers’ equation. However, because of the presence of both second and fourth order spatial derivatives, its behavior is far more complicated and interesting. The Kuramoto-Sivashinsky is a canonical model of a pattern forming system with spatio-temporal chaotic behavior. The sign of the second derivative term is such that it acts as an energy source and thus has a destabilizing effect. The nonlinear term, however, transfers energy from low to high wave numbers where the stabilizing fourth derivative term dominates. The first derivation of this equation was by Kuramoto in the study of reaction-diffusion equations modeling the Belousov-Zabotinskii reaction. The equation was also developed by Sivashinsky in higher space dimensions in modeling small thermal diffusive instabilities in laminar flame fronts and in small perturbations from a reference Poiseuille flow of a film layer on an inclined plane. In one space dimension it has also been used as a model for the problem of Bénard convection in an elongated box, and it may be used to describe long waves on the interface between two viscous fluids and unstable drift waves in plasmas. In one space dimension the Kuramoto-Sivashinsky equation reads as

where $(\lambda_{1},\lambda_{2},\lambda_{3})$ are the unknown parameters. The original dataset proposed in contains a direct numerical solution of the Kuramoto-Sivashinsky equation with 1024 spatial points and 251 time-steps. The snapshots are $\Delta t=0.4$ apart. As depicted in figure 3 using only two of these snapshots (randomly selected) with 301 and 299 data points, respectively, the algorithm is capable of identifying the correct parameter values up to a relatively good accuracy. In particular, we are using $600=301+299$ out of a total of $257024=251\times 1024$ data points in the original data set. This is possible because of equation (5) where the covariance functions explicitly encode the underlying physical laws expressed by the Kuramoto-Sivashinsky equation. For a sensitivity analysis of the reported results, let us perform the same experiment as the one illustrated in figure 3 for every pair of consecutive snapshots in the original dataset. We are still using the same number of data points (i.e., 301 and 299) for each pair of snapshots, albeit in different locations. The resulting statistics for the learned parameter values are reported in table 5. As shown in this table, more noise in the data leads to less confidence in the estimated parameter values. Moreover, to test the sensitivity of the results with respect to the gap between the two time snapshots, let us use the exact same setup as the one explained in figure 3, but increase $\Delta t$. The results are reported in table 6. These results indicate that more data, less noise, and a smaller gap $\Delta t$ between the two snapshots enhance the performance of the algorithm.

4 Nonlinear Schrödinger Equation

The one-dimensional nonlinear Schrödinger equation is a classical field equation that is used to study nonlinear wave propagation in optical fibers and/or waveguides, Bose-Einstein condensates, and plasma waves. In optics, the nonlinear term arises from the intensity dependent index of refraction of a given material. Similarly, the nonlinear term for Bose-Einstein condensates is a result of the mean-field interactions of an interacting, N-body system. The nonlinear Schrödinger equation is given by

where $(\lambda_{1},\lambda_{2})$ are the unknown parameters. Let $u$ denote the real part of $h$ and $v$ the imaginary part. Then, the nonlinear Schrödinger equation can be equivalently written as

Employing the backward Euler time stepping scheme, we obtain

The above equations can be approximated by

which involves only linear operations. Here, $u^{n-1}(x)$ and $v^{n-1}(x)$ are the real and imaginary parts of the state of the system at the previous time step, respectively. We proceed by placing two independent Gaussian processes on $u^{n}(x)$ and $v^{n}(x)$; i.e.,

Here, $\theta_{u}$ and $\theta_{v}$ are the hyper-parameters of the kernels $k_{u}$ and $k_{v}$, respectively. The prior assumptions (14) along with equations (13) enable us to encode the underlying laws of physics expressed by the nonlinear Schrödinger equation in the resulting hidden physics model

The specific forms of the covariance functions involved in model (15) is a direct function of the prior assumptions (14) as well as equations (13). The hyper-parameters $\theta_{u}$ and $\theta_{v}$ along with the parameters $\lambda_{1}$ and $\lambda_{2}$ are learned by minimizing the negative log marginal likelihood as outlined in section 4. The original data-set proposed in contains 501 time snapshots of a solution to the nonlinear Schrödinger equation with a Gaussian initial condition. The snapshots are $\Delta t=0.0063$ apart. The spatial discretization of each snapshot involves a uniform grid with 512 elements. As depicted in figure 4 using only two of these snapshots (randomly selected) with 49 and 51 data points, respectively, the algorithm is capable of identifying the correct parameter values up to a relatively good accuracy. It should be noted that we are using only $100=49+51$ data points out of a total of $256512=501\times 512$ in the original data set. Such a performance is achieved at the cost of explicitly encoding the underlying physical laws expressed by the nonlinear Schrödinger equation in the covariance functions of the hidden physics model (15). For a systematic study of the performance of the method, let us carry out the same experiment as the one illustrated in figure 4 for every pair of consecutive snapshots in the original dataset. We are still using the same number of data points (i.e., 49 and 51) for each pair of snapshots. The resulting statistics for the learned parameter values are reported in table 7. As is clearly demonstrated in this table, more noise in the data leads to less confidence in the estimated values for the parameters. Moreover, let us recall the main assumption of this work that the gap $\Delta t$ between the pair of snapshots should be small enough so that we can employ the backward Euler scheme (see equation (12)). To test the importance of this assumption, let us use the exact same setup as the one explained in figure 4, but increase $\Delta t$. The results are reported in table 8. Therefore, the most important facts about the proposed methodology are that more data, less noise, and a smaller gap $\Delta t$ between the two snapshots enhance the performance of the algorithm.

5 Navier-Stokes Equations

Navier-Stokes equations describe the physics of many phenomena of scientific and engineering interest. They may be used to model the weather, ocean currents, water flow in a pipe and air flow around a wing. The Navier-Stokes equations in their full and simplified forms help with the design of aircraft and cars, the study of blood flow, the design of power stations, the analysis of the dispersion of pollutants, and many other applications. Let us consider the Navier-Stokes equations in two dimensionsIt is straightforward to generalize the proposed framework to the Navier-Stokes equations in three dimensions (3D). (2D) given explicitly by

where $u(t,x,y)$ denotes the $x$-component of the velocity field, $v(t,x,y)$ the $y$-component, and $p(t,x,y)$ the pressure. Here, $\lambda=(\lambda_{1},\lambda_{2})$ are the unknown parameters. Solutions to the Navier-Stokes equations are searched in the set of divergence-free functions; i.e.,

This extra equation is the continuity equation for incompressible fluids that describes the conservation of mass of the fluid. Applying the backward Euler time stepping scheme to the Navier-Stokes equations (16) we obtain

where $u^{n}(x,y)=u(t^{n},x,y)$ and $v^{n}(x,y)=v(t^{n},x,y)$. We make the assumption that

for some latent function $\psi^{n}(x,y)$. Under this assumption, the continuity equation (17) will be automatically satisfied. We proceed by placing a Gaussian process prior on

where $\theta$ are the hyper-parameters of the kernel $k((x,y),(x^{\prime},y^{\prime});\theta)$. This will result in the following multi-output Gaussian process

By construction (see equation (19)), any samples generated from this multi-output Gaussian process will satisfy the continuity equation (17). Moreover, independent from $\psi^{n}(x,y)$, we will place a Gaussian process prior on $p^{n}(x,y)$; i.e.,

We linearize the backward Euler time stepping scheme by employing the states $u^{n-1}(x,y)$ and $v^{n-1}(x,y)$ of the system at the previous time step and writing

The above equations (23) can be rewritten as

by defining the linear operator $\mathcal{L}_{(x,y)}^{\lambda}$ to be given by

This will allow us to obtain the following hidden physics model encoding the structure of the Navier-Stokes equations and the backward Euler time stepping scheme in its kernels; i.e.,

The lower triangular portion of the matrix of covariance functions (26) is not shown due to symmetry. The hyper-parameters $\theta$ and $\theta_{p}$ along with the parameters $\lambda=(\lambda_{1},\lambda_{2})$ are learned by minimizing the negative log marginal likelihood as outlined in section 4. As for the data, following the exact same instructions as the ones provided in and , we simulate the Navier-Stokes equations describing the two-dimensional fluid flow past a circular cylinder at Reynolds number 100 using the Immersed Boundary Projection Method . This approach utilizes a multi-domain scheme with four nested domains, each successive grid being twice as large as the previous one. Length and time are nondimensionalized so that the cylinder has unit diameter and the flow has unit velocity. Data is collected on the finest domain with dimensions $9\times 4$ at a grid resolution of $449\times 199$. The flow solver uses a 3rd-order Runge Kutta integration scheme with a time step of t = 0.02, which has been verified to yield well-resolved and converged flow fields. After simulations converge to steady periodic vortex shedding, flow snapshots are saved every $\Delta t=0.02$. As depicted in figure 5 using only two snapshots of the velocityIt is worth emphasizing that we are not making use of any data on the pressure or vorticity fields. In practice, unlike velocity (e.g., Particle Image Velocimetry (PIV) data), obtaining direct measurements of the pressure or vorticity fields are more demanding if not impossible. Our method circumvents the need for having data on the pressure simply because of the prior assumption (21) where any samples generated from this multi-output Gaussian process satisfy the continuity equation (17). field with 251 and 249 data points, respectively, the algorithm is capable of identifying the correct parameter values up to a relatively good accuracy. It should be noted that we are using only two snapshots with a total of $500=251+249$ data points. This surprising performance is achieved at the cost of explicitly encoding the underlying physical laws expressed by the Navier-Stokes equations in the covariance functions of the hidden physics model (26). For a sensitivity analysis of the reported results, let us perform the same experiment as the one illustrated in figure 5 for 501 pairs of consecutive snapshots. We are still using the same number of data points (i.e., 251 and 249) for each pair of snapshots. The resulting statistics for the learned parameter values are reported in table 9. As is clearly demonstrated in this table, more noise in the data leads to less confidence in the estimated values for the parameters. Moreover, to test the sensitivity of the results with respect to the gap between two time snapshots, let us use the exact same setup as the one explained in figure 5, but increase $\Delta t$. The results are reported in table 10. These results verify the most important facts about the proposed methodology that more data, less noise, and a smaller gap $\Delta t$ between the two snapshots enhance the performance of the algorithm. In particular, the results reported in table 10 indicate that to obtain more accurate estimates of the Reynolds number $1/\lambda_{2}$ one needs to utilize a smaller gap $\Delta t$ between the pair of snapshots. To verify the validity of this conjecture let us decrease the gap $\Delta t$ between the pair of time snapshots while employing the exact same setup as the one explained in figure 5. The results are reported in table 11. As is clearly demonstrated in this table, a smaller $\Delta t$ leads to more accurate estimates of the Reynolds number $1/\lambda_{2}$ in the absence of noise in the data. However, a smaller $\Delta t$ seems to make the algorithm more susceptible to noise in the data.

6 Fractional Equations

Let us consider the one dimensional fractional equation

where $(\lambda_{1},\lambda_{2})$ are the unknown parameters. In particular, $\lambda_{2}$ is the fractional order of the operator $\mathcal{D}^{\lambda_{2}}_{-\infty,x}$ that is defined in the Riemann-Liouville sense . Fractional operators often arise in modeling anomalous diffusion processes and other non-local interactions. Integer values such as $\lambda_{2}=1$ and $\lambda_{2}=2$ can model classical advection and diffusion phenomena, respectively. However, under the fractional calculus setting, $\lambda_{2}$ can assume real values and thus continuously interpolate between inherently different model behaviors. The proposed framework allows $\lambda_{2}$ to be directly inferred from noisy data, and opens the path to a flexible formalism for model discovery and calibration. Applying the backward Euler time stepping scheme to equation (27) we obtain

Here, $u^{n}(x)=u(t^{n},x)$ is the hidden state of the system at time $t^{n}$. We make the prior assumption that

The prior assumption (29) along with the backward Euler scheme (28) allow us to obtain the following hidden physics model corresponding to the fractional equation (27); i.e.,

The only technicality induced by fractional operators has to do with deriving the kernels $k^{n,n-1}$, $k^{n-1,n}$, and $k^{n-1,n-1}$. Here, $k^{n,n-1}(x,x^{\prime};\theta,\lambda_{1},\lambda_{2})$ was obtained by taking the inverse Fourier transform of

where $\widehat{k}(w,w^{\prime};\theta)$ is the Fourier transform of the kernel $k(x,x^{\prime};\theta)$. Similarly, one can obtain $k^{n-1,n}$ and $k^{n-1,n-1}$. The hyper-parameters $\theta$ along with the parameters $\lambda_{1}$ and $\lambda_{2}$ are learned by minimizing the negative log marginal likelihood as outlined in section 4. We use the hidden physics model (30) to identify the long celebrated relation between Brownian motion and the diffusion equation . The Fokker-Planck equation for a Brownian motion with $x(t+\Delta t)\sim\mathcal{N}(x(t),dt)$, associated with a particle’s position, is $u_{t}=0.5u_{xx}$. We simulated a Brownian motion at evenly spaced time points and generated two histograms of the particle’s displacement. These two histograms are $\Delta t=0.01$ apart. As depicted in figure 6 using only two histograms with 100 bins for each one, the algorithm is capable of identifying the correct fractional order and parameter values up to a relatively good accuracy. Moreover, let us now consider the one dimensional fractional equation

where $\alpha$ is the unknown parameter and $(-\nabla^{\alpha}_{x})$ is the fractional Laplacian operator . The fractional Laplacian is the operator with symbol $|w|^{\alpha}$. In other words, the Fourier transform of $(-\nabla^{\alpha}_{x})u(x)$ is given by $|w|^{\alpha}\widehat{u}(w)$. The fractional Laplacian operator can also be defined as the generator of $\alpha$-stableStable distributions are a rich class of probability distributions that allow skewness and heavy tails. Stable distributions have been proposed as a model for many types of physical and economic systems. In particular, it is argued that some observed quantities are the sum of many small terms – the price of a stock, the noise in a communication system, etc. – and hence a stable model should be used to describe such systems. Lévy processes. Motivated by this observation, we simulated an $\alpha$-stable Lévy process and employed the hidden physics model resulting from equation (31) to identify the fractional order $\alpha$. As depicted in figure 7 using only two histograms with 100 bins for each one, the algorithm is capable of identifying the correct fractional order up to a relatively good accuracy.

Summary and Discussion

We have introduced a structured learning machine which is explicitly informed by the underlying physics that possibly generated the observed data. Exploiting this structure is critical for constructing data-efficient learning algorithms that can effectively distill information in the data-scarce scenarios appearing routinely when we study complex physical systems. We applied the proposed framework to the problem of identifying general parametric nonlinear partial differential equations from noisy data. This generality was demonstrated using various benchmark problems with different attributes. This work should be considered a direct follow up on in which a similar methodology was employed to infer solutions to time-dependent and nonlinear partial differential equations, and effectively quantify and propagate uncertainty due to noisy initial or boundary data. The ideas introduced in these two papers provide a natural platform for learning from noisy data and computing under uncertainty. Perhaps the most pressing limitation of this work in its present form stems from the cubic scaling with respect to the total number of training data points. However, ideas such as recursive Kalman updates , variational inference , and parametric Gaussian processes can be used to address this limitation.

Moreover, the examples studied in the current work were inspired by the pioneering work recently presented in . The authors of followed a sparse regression approach and a full set of spatio-temporal time series measurements consisting of thousands of data points. In contrast, here we used much smaller datasets with only hundreds of points and two snapshots of the systems. However, unlike the work in , here we did not use a dictionary of all possible terms involved in the partial differential equation. We could possibly include such a dictionary in our formulation but that would make our kernel evaluations more expensive. Moreover, in some systems, e.g., in an advection-diffusion-reaction system we know most of the terms of the equation, i.e., advection and diffusion but typically the reaction term is unknown. In this case, we would seek to obtain the parameters in front of the advection-diffusion and discover the functional form of the reaction term along with any parameters using the methodology outline in this paper. In comparison to , our method does not require numerical differentiation as the kernels are obtained analytically. Moreover, we do not require a regular lattice as in and can work with scattered data. An additional advantage of our approach is that it can estimate parameters appearing anywhere in the formulation of the partial differential equation while the method of is only suitable for parameters appearing as coefficients. For example, they cannot estimate the fractional order in the last example we presented in our paper or the parameters of partial differential equations (e.g., the sine-Gordon equation) involving a term like $\sin(\lambda u(x))$ with $\lambda$ being the parameter. Also, the treatment of the noise is somewhat complex in the method of as it involves some sort of filtering via e.g., singular value decomposition whereas our method can filter arbitrarily noisy data automatically via the Gaussian process prior assumptions. We believe that both methods can be used in different contexts effectively and we anticipate that this is only the beginning of a new way of thinking and formulating new and possibly simpler equations, e.g., by employing fractional operators that are naturally captured in our framework.

Acknowledgements

This work received support by the DARPA EQUiPS grant N66001-15-2-4055, the MURI/ARO grant W911NF-15-1-0562, and the AFOSR grant FA9550-17-1-0013. All data and codes used in this manuscript are publicly available on GitHub at https://github.com/maziarraissi/HPM.

References