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 dif-ficulties 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 [1] 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 [2].
Let us consider parametrized and nonlinear partial differential equations of the general form
where h(t, x) denotes the latent (hidden) solution, is a nonlinear operator parametrized by
, and Ω is a subset of
. As an example, the one dimensional Burgers’ equation corresponds to the case where
and
). Here, the subscripts denote partial differentiation in either time or space. Given noisy measurements of the system, one is typically interested in the solution of two distinct problems. The first problem is that of inference or filtering and smoothing, which states: given fixed model parameters
what can be said about the unknown hidden state h(t, x) of the system? This question is the topic of a preceding paper [1] of the authors in which we introduce the concept of numerical Gaussian processes and address the problem of inferring solutions of time dependent and nonlinear partial differential equations using noisy observations. The second problem is that of learning, system identification, or data driven discovery of partial differ-ential equations [2] stating: what are the parameters
that best describe the observed data? Here we assume that all we observe are two snapshots
of the system at times
, respectively, which are ∆
apart. The main assumption is that ∆t is small enough so that we can apply the backward Euler time stepping scheme
equation (1) and obtain the discretized equation
Here, ) is the hidden state of the system at time
proximating the nonlinear operator on the left-hand-side of equation (2) by a linear one we obtain
For instance, the nonlinear operator
involved in the Burgers’ equation can be approximated by the linear operator
where ) is the state of the system at the previous time
Similar to Raissi et al. [3, 4], 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 processprior over the latent function
Here, 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 exponential
covariance function, i.e.,
where ) 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
in the resulting multi-output Gaussian process
It is worth highlighting that the parameters of the operators
turn into hyper-parameters of the resulting covariance functions. The specific
forms of the kernels
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).
Given the noisy data on the latent solution at times
, respectively, the hyper-parameters
of the covariance functions and more importantly the parameters
of the operators
can be learned by employing a Quasi-Newton optimizer L-BFGS [22] to minimize the negative log marginal likelihood [5]
Here, N is the total number of data points in h. Moreover, 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
) 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
in equation (6) targets fitting the training data, while the log-determinant term log |K| penalizes model complexity. This regularization mechanism automatically meets the Occam’s razor principle [23] 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 few
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 [25] and more robust estimate of the linear operator parameters
can be obtained by assigning priors on
. However, this would require more costly sampling procedures such as Markov Chain Monte Carlo (see [5], 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 K. This scales cubically with the number N of training data in h. While it has been effectively addressed by the recent works of [26, 27, 28], this cubic scaling is still a well-known limitation of Gaussian process regression.
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 [2], 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¨odinger, and Kuramoto-Sivashinsky equations. Moreover, all data and codes used in this manuscript are publicly available on GitHub at https://github.com/maziarraissi/HPM.
5.1. Burgers’ Equation
Burgers’ equation arises in various areas of applied mathematics, including fluid mechanics, nonlinear acoustics, gas dynamics, and traffic flow [29]. 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 () being the unknown parameters. The original data-set proposed in [2] contains 101 time snapshots of a solution to the Burgers’ equation with a Gaussian initial condition, propagating into a traveling wave. The snapshots are ∆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 =
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 ∆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 ∆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 ∆t between the two snapshots enhance the performance of the algorithm.
Figure 1: A solution to the Burgers’ equation is depicted in the top panel. The two white vertical lines in this panel specify the locations of the two randomly selected snapshots. These two snapshots are ∆t = 0.1 apart and are plotted in the middle panel. The red crosses denote the locations of the training data points. The correct partial differential equation along with the identified ones are reported in the lower panel.
Table 1: Resulting statistics for the learned parameter values.
Table 2: Effect of increasing the gap ∆t between the pair of snapshots.
5.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 [30] in the continuum limit. The KdV equation reads as
with () being the unknown parameters. The original dataset proposed in [2] contains a two soliton solution to the KdV equation with 512 spatial points and 201 time-steps. The snapshots are ∆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
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 ∆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 ∆t between the two snapshots enhance the performance of the algorithm.
Figure 2: The KdV equation: A solution to the KdV equation is depicted in the top panel. The two white vertical lines in this panel specify the locations of the two randomly selected snapshots. These two snapshots are ∆t = 0.1 apart and are plotted in the middle panel. The red crosses denote the locations of the training data points. The correct partial differential equation along with the identified ones are reported in the lower panel.
Table 3: The KdV equation: Resulting statistics for the learned parameter values.
Table 4: The KdV equation: Effect of increasing the gap ∆t between the pair of snapshots.
5.3. Kuramoto-Sivashinsky Equation
The Kuramoto-Sivashinsky equation [31, 32, 33] 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´enard 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 () are the unknown parameters. The original dataset proposed in [2] contains a direct numerical solution of the Kuramoto-Sivashinsky equation with 1024 spatial points and 251 time-steps. The snapshots are ∆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
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 con-fidence 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 ∆t. The results are reported in table 6. These results indicate that more data, less noise, and a smaller gap ∆t between the two snapshots enhance the performance of the algorithm.
Figure 3: Kuramoto-Sivashinsky equation: A solution to the Kuramoto-Sivashinsky equation is depicted in the top panel. The two white vertical lines in this panel specify the locations of the two randomly selected snapshots. These two snapshots are ∆t = 0.4 apart and are plotted in the middle panel. The red crosses denote the locations of the training data points. The correct partial differential equation along with the identified ones are reported in the lower panel.
Table 5: Kuramoto-Sivashinsky equation: Resulting statistics for the learned parameter values.
Table 6: Kuramoto-Sivashinsky equation: Effect of increasing the gap ∆t between the pair of snapshots.
5.4. Nonlinear Schr¨odinger Equation
The one-dimensional nonlinear Schr¨odinger 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¨odinger equation is given by
where () are the unknown parameters. Let u denote the real part of h and v the imaginary part. Then, the nonlinear Schr¨odinger 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, ) 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
Here, are the hyper-parameters of the kernels
tively. The prior assumptions (14) along with equations (13) enable us to encode the underlying laws of physics expressed by the nonlinear Schr¨odinger 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 along with the parameters
are learned by minimizing the negative log marginal likelihood as outlined in section 4. The original data-set proposed in [2] contains 501 time snapshots of a solution to the nonlinear Schr¨odinger equation with a Gaussian initial condition. The snapshots are ∆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
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¨odinger 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 ∆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 ∆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 ∆t between the two snapshots enhance the performance of the algorithm.
Table 7: Resulting statistics for the learned parameter values.
Table 8: Effect of increasing the gap ∆t between the pair of snapshots.
Figure 4: A solution to the nonlinear Schr¨odinger equation is depicted in the top two panels. The two black vertical lines in these two panels specify the locations of the two randomly selected snapshots. These two snapshots are ∆t = 0.0063 apart and are plotted in the two middle panels. The red crosses denote the locations of the training data points. The correct partial differential equation along with the identified ones are reported in the lower panel. Here, u is the real part of h and v is the imaginary part.
5.5. Navier-Stokes Equations
Navier-Stokes equations describe the physics of many phenomena of sci-entific 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 NavierStokes 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 dimensions(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, ) 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 ). We make the assumption that
for some latent function Under this assumption, the continuity equation (17) will be automatically satisfied. We proceed by placing a Gaussian process prior on
where are the hyper-parameters of the kernel
). This will result in the following multi-output Gaussian process
where
By construction (see equation (19)), any samples generated from this multi-output Gaussian process will satisfy the continuity equation (17). Moreover, independent from ), we will place a Gaussian process prior on
i.e.,
We linearize the backward Euler time stepping scheme by employing the states ) of the system at the previous time step and writing
The above equations (23) can be rewritten as
by defining the linear operator 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
where
and
The lower triangular portion of the matrix of covariance functions (26) is not shown due to symmetry. The hyper-parameters along with the parameters
) 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 [34] and [2], we simulate the NavierStokes equations describing the two-dimensional fluid flow past a circular cylinder at Reynolds number 100 using the Immersed Boundary Projection Method [35, 36]. 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
4 at a grid resolution of 449
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 ∆t = 0.02. As depicted in figure 5 using only two snapshots of the velocity
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 ∆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 ∆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
to utilize a smaller gap ∆t between the pair of snapshots. To verify the validity of this conjecture let us decrease the gap ∆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 ∆t leads to more accurate estimates of the Reynolds number 1
in the absence of noise in the data. However, a smaller ∆t seems to make the algorithm more susceptible to noise in the data.
Figure 5: Navier-Stokes equations: A single snapshot of the vorticity field of a solution to the Navier-Stokes equations for the fluid flow past a cylinder is depicted in the top panel. The black box in this panel specifies the sampling region. Two snapshots of the velocity field being ∆t = 0.02 apart are plotted in the two middle panels. The black crosses denote the locations of the training data points. The correct partial differential equation along with the identified ones are reported in the lower panel. Here, u denotes the x-component of the velocity field, v the y-component, p the pressure, and w the vorticity field.
Table 9: Navier-Stokes equations: Resulting statistics for the learned parameter values.
Table 10: Navier-Stokes equations: Effect of increasing the gap ∆t between the pair of snapshots.
Table 11: Navier-Stokes equations: Effect of decreasing the gap ∆t between the pair of snapshots.
5.6. Fractional Equations
Let us consider the one dimensional fractional equation
where () are the unknown parameters. In particular,
is the fractional order of the operator
that is defined in the Riemann-Liouville sense [37]. Fractional operators often arise in modeling anomalous diffusion processes and other non-local interactions. Integer values such as
However, under the fractional calculus setting,
can assume real values and thus continuously interpolate between inherently different model behaviors. The proposed framework allows
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, ) is the hidden state of the system at time
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 obtained by taking the inverse Fourier transform [37] of
where ) is the Fourier transform of the kernel
). Similarly, one can obtain
. The hyper-parameters
along with the parameters
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 diffu-sion equation [2]. The Fokker-Planck equation for a Brownian motion with
), associated with a particle’s position, is
We simulated a Brownian motion at evenly spaced time points and generated two histograms of the particle’s displacement. These two histograms are ∆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 is the unknown parameter and (
) is the fractional Laplacian operator [37]. The fractional Laplacian is the operator with symbol
In other words, the Fourier transform of (
) is given by
The fractional Laplacian operator can also be defined as the generator of
L´evy processes. Motivated by this observation, we simulated an
stable L´evy process [39, 40] and employed the hidden physics model resulting from equation (31) to identify the fractional order
. 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.
Figure 6: Fractional Equation – Brownian Motion: A single realization of a Brownian motion is depicted in the top panel. Two histograms of the particle’s displacement, being ∆t = 0.01 apart, are plotted in the middle panel. The correct partial differential equation along with the identified ones are reported in the lower panel.
Figure 7: A single realization of an
L´evy process is depicted in the top panel. Two histograms of the particle’s displacement, being ∆t = 0.01 apart, are plotted in the middle panel. The correct partial differential equation along with the identified ones are reported in the lower panel.
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 [1] in which a similar methodology was employed to infer solutions to time-dependent and non-linear 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 [41], variational inference [27], and parametric Gaussian processes [28] can be used to address this limitation.
Moreover, the examples studied in the current work were inspired by the pioneering work recently presented in [2]. The authors of [2] 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 [2], 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 [2], our method does not require numerical differentiation as the kernels are obtained analytically. Moreover, we do not require a regular lattice as in [2] 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 [2] 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(being the parameter. Also, the treatment of the noise is somewhat complex in the method of [2] 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.
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.
[1] M. Raissi, P. Perdikaris, G. E. Karniadakis, Numerical gaussian pro- cesses for time-dependent and non-linear partial differential equations, arXiv preprint arXiv:1703.10230 (2017).
[2] S. H. Rudy, S. L. Brunton, J. L. Proctor, J. N. Kutz, Data-driven discovery of partial differential equations, Science Advances 3 (2017).
[3] M. Raissi, P. Perdikaris, G. E. Karniadakis, Inferring solutions of dif- ferential equations using noisy multi-fidelity data, Journal of Computational Physics 335 (2017) 736 – 746.
[4] M. Raissi, P. Perdikaris, G. E. Karniadakis, Machine learning of linear differential equations using gaussian processes, Journal of Computational Physics 348 (2017) 683 – 693.
[5] C. E. Rasmussen, C. K. Williams, Gaussian processes for machine learn- ing, volume 1, MIT press Cambridge, 2006.
[6] K. P. Murphy, Machine learning: a probabilistic perspective, MIT press, 2012.
[7] R. M. Neal, Bayesian learning for neural networks, volume 118, Springer Science & Business Media, 2012.
[8] V. Vapnik, The nature of statistical learning theory, Springer Science & Business Media, 2013.
[9] B. Sch¨olkopf, A. J. Smola, Learning with kernels: support vector ma- chines, regularization, optimization, and beyond, MIT press, 2002.
[10] M. E. Tipping, Sparse Bayesian learning and the relevance vector ma- chine, The journal of machine learning research 1 (2001) 211–244.
[11] A. Tikhonov, Solution of incorrectly formulated problems and the reg- ularization method, in: Soviet Math. Dokl., volume 5, pp. 1035–1038.
[12] A. N. Tikhonov, V. Y. Arsenin, Solutions of Ill-posed problems, W.H. Winston, 1977.
[13] T. Poggio, F. Girosi, Networks for approximation and learning, Pro- ceedings of the IEEE 78 (1990) 1481–1497.
[14] N. Aronszajn, Theory of reproducing kernels, Transactions of the Amer- ican Mathematical Society 68 (1950) 337–404.
[15] S. Saitoh, Theory of reproducing kernels and its applications, volume 189, Longman, 1988.
[16] A. Berlinet, C. Thomas-Agnan, Reproducing kernel Hilbert spaces in probability and statistics, Springer Science & Business Media, 2011.
[17] D. Duvenaud, J. R. Lloyd, R. Grosse, J. B. Tenenbaum, Z. Ghahramani, Structure discovery in nonparametric regression through compositional kernel search, arXiv preprint arXiv:1302.4922 (2013).
[18] R. Grosse, R. R. Salakhutdinov, W. T. Freeman, J. B. Tenenbaum, Exploiting compositionality to explore a large space of model structures, arXiv preprint arXiv:1210.4856 (2012).
[19] G. Malkomes, C. Schaff, R. Garnett, Bayesian optimization for auto- mated model selection, in: Advances in Neural Information Processing Systems, pp. 2900–2908.
[20] R. Calandra, J. Peters, C. E. Rasmussen, M. P. Deisenroth, Manifold Gaussian processes for regression, in: Neural Networks (IJCNN), 2016 International Joint Conference on, IEEE, pp. 3338–3345.
[21] M. Raissi, G. Karniadakis, Deep multi-fidelity Gaussian processes, arXiv preprint arXiv:1604.07484 (2016).
[22] D. C. Liu, J. Nocedal, On the limited memory BFGS method for large scale optimization, Mathematical programming 45 (1989) 503–528.
[23] C. E. Rasmussen, Z. Ghahramani, Occam’s razor, Advances in neural information processing systems (2001) 294–300.
[24] S. L. Brunton, J. L. Proctor, J. N. Kutz, Discovering governing equa- tions from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 113 (2016) 3932–3937.
[25] A. M. Stuart, Inverse problems: a Bayesian perspective, Acta Numerica 19 (2010) 451–559.
[26] E. Snelson, Z. Ghahramani, Sparse Gaussian processes using pseudo- inputs, in: Advances in neural information processing systems, pp. 1257–1264.
[27] J. Hensman, N. Fusi, N. D. Lawrence, Gaussian processes for big data, in: Proceedings of the Twenty-Ninth Conference on Uncertainty in Ar-tificial Intelligence, UAI 2013, Bellevue, WA, USA, August 11-15, 2013.
[28] M. Raissi, Parametric gaussian process regression for big data, arXiv preprint arXiv:1704.03144 (2017).
[29] C. Basdevant, M. Deville, P. Haldenwang, J. Lacroix, J. Ouazzani, R. Peyret, P. Orlandi, A. Patera, Spectral and finite difference solutions of the Burgers equation, Computers & fluids 14 (1986) 23–41.
[30] T. Dauxois, Fermi, Pasta, Ulam and a mysterious lady, arXiv preprint arXiv:0801.1590 (2008).
[31] J. M. Hyman, B. Nicolaenko, The Kuramoto-Sivashinsky equation: a bridge between pde’s and dynamical systems, Physica D: Nonlinear Phenomena 18 (1986) 113–126.
[32] B. I. Shraiman, Order, disorder, and phase turbulence, Physical review letters 57 (1986) 325.
[33] B. Nicolaenko, B. Scheurer, R. Temam, Some global dynamical prop- erties of the Kuramoto-Sivashinsky equations: nonlinear stability and attractors, Physica D: Nonlinear Phenomena 16 (1985) 155–183.
[34] J. N. Kutz, S. L. Brunton, B. W. Brunton, J. L. Proctor, Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems, volume 149, SIAM, 2016.
[35] K. Taira, T. Colonius, The immersed boundary method: a projection approach, Journal of Computational Physics 225 (2007) 2118–2137.
[36] T. Colonius, K. Taira, A fast immersed boundary method using a nullspace approach and multi-domain far-field boundary conditions, Computer Methods in Applied Mechanics and Engineering 197 (2008) 2131–2146.
[37] I. Podlubny, Fractional differential equations: an introduction to frac- tional derivatives, fractional differential equations, to methods of their solution and some of their applications, volume 198, Academic press, 1998.
[38] J. Nolan, Stable distributions: models for heavy-tailed data, Birkhauser New York, 2003.
[39] J. M. Chambers, C. L. Mallows, B. Stuck, A method for simulating stable random variables, Journal of the american statistical association 71 (1976) 340–344.
[40] A. Weron, R. Weron, Computer simulation of l´evy -stable variables and processes, ChaosThe Interplay Between Stochastic and Deterministic Behaviour (1995) 379–392.
[41] J. Hartikainen, S. S¨arkk¨a, Kalman filtering and smoothing solutions to temporal Gaussian process regression models, in: Machine Learning for Signal Processing (MLSP), 2010 IEEE International Workshop on, IEEE, pp. 379–384.