Numerical Gaussian Processes for Time-dependent and Non-linear Partial Differential Equations

Introduction

Data-driven methods are taking center stage across many disciplines of science, and machine learning techniques have achieved groundbreaking results across a diverse spectrum of pattern recognition tasks . Despite their disruptive implications, many of these methods are blind to any underlying laws of physics that may have shaped the distribution of the observed data. A natural question would then be how one can construct efficient learning machines that explicitly leverage such structured prior information? To answer this question we have to turn our attention to the immense collective knowledge originating from centuries of research in applied mathematics and mathematical physics. Modeling the physical world through the lens of mathematics typically translates into deriving conservation laws from first principles, which often take the form of systems of partial differential equations. In many practical settings, the solution of such systems is only accessible by means of numerical algorithms that provide sensible approximations to given quantities of interest. In this work, we aim to capitalize on the long-standing developments of classical methods in numerical analysis and revisit partial differential equations from a statistical inference viewpoint. The merits of this approach are twofold. First, it enables the construction of data-efficient learning machines that can encode physical conservation laws as structured prior information. Second, it allows the design of novel numerical algorithms that can seamlessly blend equations and noisy data, infer latent quantities of interest (e.g., the solution to a partial differential equation), and naturally quantify uncertainty in computations. This approach is aligned in spirit with the emerging field of probabilistic numerics , which roots all the way back to Poincaré’s courses on probability theory , and has been recently revived by the pioneering works of .

To illustrate the key ingredients of this study, let us start by considering linearNon-linear equations have to be studied on a case by case basis (see e.g., section 2.6). partial differential equations of the form

Here, $u^{n}(x)=u(t^{n},x)$. Building upon Raissi et al. , we place a Gaussian process prior on $u^{n-1}$, i.e.,

Here, $\theta$ denotes the hyper-parameters of the covariance function $k_{u,u}^{n-1,n-1}$. Gaussian process regression (see ) is a non-parametric Bayesian machine learning technique that provides a flexible prior distribution over functions, enjoys analytical tractability, and has a fully probabilistic work-flow that returns robust posterior variance estimates, which quantify uncertainty in a natural way. Moreover, Gaussian processes are among a class of methods known as kernel machines (see ) and are analogous to regularization approaches (see ). They can also be viewed as a prior on one-layer feed-forward Bayesian neural networks with an infinite number of hidden units . The Gaussian process prior assumption (3) along with the Euler scheme (2) will allow us to capture the entire structure of the differential operator $\mathcal{L}_{x}$ as well as the Euler time-stepping rule in the resulting multi-output Gaussian process

The specific forms of the kernels $k^{n,n}_{u,u}$ and $k^{n,n-1}_{u,u}$ are direct functions of the Euler scheme (2) as well as the prior assumption (3), and will be discussed in more detail later. The multi-output process (4) is an example of a numerical Gaussian process, because the covariance functions $k^{n,n}_{u,u}$ and $k^{n,n-1}_{u,u}$ result from a numerical scheme, in this case, the Euler method. Essentially, this introduces a structured prior that explicitly encodes the physical law modeled by the partial differential equation (1). In the following, we will generalize the framework outlined above to arbitrary linear multi-step methods, originally proposed by Bashforth and Adams , as well as Runge-Kutta methods, generally attributed to Runge . The biggest challenge here is the proper placement of the Gaussian process prior (see e.g., equation (3)) in order to avoid inversion of differential operators and to bypass the classical need for spatial discretization of such operators. For instance, in the above example (see equations (2) and (3)), it would have been an inappropriate choice to start by placing a Gaussian process prior on $u^{n}$, rather than on $u^{n-1}$, as obtaining the numerical Gaussian process (4) would then involve inverting operators of the form $I+\Delta t\mathcal{L}_{x}$ corresponding to the Euler method. Moreover, propagating the uncertainty associated with the noisy initial observations $\{\bm{x}^{0},\bm{u}^{0}\}$ through time is another major challenge addressed in the following.

Linear Multi-step Methods

Let us start with the most general form of the linear multi-step methods applied to equation (1); i.e.,

Different choices for the parameters $\alpha_{i}$ and $\beta_{i}$ result in specific schemes. For instance, in table 1, we present some specific members of the family of linear multi-step methods (5).

We encourage the reader to keep these special cases in mind while reading the rest of this section. Linear multi-step methods (5) can be equivalently written as

where $\mathcal{P}_{x}u:=u-\Delta t\beta_{0}\mathcal{L}_{x}u$ and $\mathcal{Q}^{i}_{x}u:=\alpha_{i}u+\Delta t\beta_{i}\mathcal{L}_{x}u$. Some special cases of equation (6) are given in table 2.

For every $j=0,1,\ldots,m$ and some $\tau\in$ which depends on the specific choices for the values of the parameters $\alpha_{i}$ and $\beta_{i}$, we define $u^{n-j+\tau}$ to be given by

Definition (7) takes the specific forms given in table 3 for some example schemes.

Shifting every term involved in the above definition (7) by $-\tau$ yields

To give an example, for the trapezoidal rule we obtain $\mathcal{P}_{x}u^{n+1/2}=u^{n}=\mathcal{Q}_{x}u^{n-1/2}$ and $\mathcal{P}_{x}u^{n-1/2}=u^{n-1}=\mathcal{Q}_{x}u^{n-3/2}$. Therefore, as a direct consequence of equation (8) we have

This, in the special case of the trapezoidal rule, translates to $u^{n}=\mathcal{Q}_{x}u^{n-1/2}$ and $u^{n-1}=\mathcal{P}_{x}u^{n-1/2}$. It is worth noting that by assuming $u^{n-1/2}(x)\sim\mathcal{GP}(0,k(x,x^{\prime};\theta))$, we can capture the entire structure of the trapezoidal rule in the resulting joint distribution of $u^{n}$ and $u^{n-1}$. This proper placement of the Gaussian process prior is key to the proposed methodology as it allows us to avoid any spatial discretization of differential operators since no inversion of such operators is necessary. We will capitalize on this idea in the following.

are $m$ independent processes, we obtain the following numerical Gaussian process

It is worth noting that the entire structure of linear multi-step methods (5) is captured by the kernels given in equations (11). Note that although we start from an independence assumption in equation (10), the resulting numerical Gaussian process exhibits a fully correlated structure as illustrated in equations (11). Moreover, the information on the boundary $\partial\Omega$ of the domain $\Omega$ can often be summarized by noisy observations $\{\bm{x}^{n}_{b},\bm{u}^{n}_{b}\}$ of a linear transformation $\mathcal{B}_{x}$ of $u^{n}$; i.e., noisy data on

Using this, we obtain the following covariance functions involving the boundary

The numerical examples accompanying this manuscript are designed to showcase different special treatments of boundary conditions, including Dirichlet, Neumann, mixed, and periodic boundary conditions.

2 Work flow and computational cost

The proposed work flow is summarized below:

Starting from the initial data $\{\bm{x}^{0},\bm{u}^{0}\}$ and the boundary data $\{\bm{x}^{1}_{b},\bm{u}^{1}_{b}\}$, we train the kernel hyper-parameters as outlined in section 2.3. This step carries the main computational burden as it scales cubically with the total number of training points since it involves Cholesky factorization of full symmetric positive-definite covariance matrices .

Having identified the optimal set of kernel hyper-parameters, we utilize the conditional posterior distribution to predict the solution at the next time-step and generate the artificial data $\{\bm{x}^{1},\bm{u}^{1}\}$. Note that $\bm{x}^{1}$ is randomly sampled in the spatial domain according to a uniform distribution, and $\bm{u}^{1}$ is a normally distributed random vector, as outlined in section 2.4.

Given the artificial data $\{\bm{x}^{1},\bm{u}^{1}\}$ and boundary data $\{\bm{x}^{2}_{b},\bm{u}^{2}_{b}\}$ we proceed with training the kernel hyper-parameters for the second time-stepTo be precise, we are using the mean of the random vector $\bm{u}^{1}$ for training purposes. (see section 2.3).

Having identified the optimal set of kernel hyper-parameters, we utilize the conditional posterior distribution to predict the solution at the next time-step and generate the artificial data $\{\bm{x}^{2},\bm{u}^{2}\}$, where $\bm{x}^{2}$ is randomly sampled in the spatial domain according to a uniform distribution. However, since $\bm{u}^{1}$ is a random vector, we have to marginalize it out in order to obtain consistent uncertainty estimates for $\bm{u}^{2}$. This procedure is outlined in section 2.5.

Steps 3 and 4 are repeated until the final integration time is reached.

In summary, the proposed methodology boils down to a sequence of Gaussian process regressions at every time-step. To accelerate training, one can use the optimal set of hyper-parameters from the previous time-step as an initial guess for the current one.

3 Training

In the following, for notational convenience and without loss of generalityThe reader should be able to figure out the details without much difficulty while generalizing to cases with $m>1$. Moreover, for the examples accompanying this manuscript, more details are also provided in the appendix., we will operate under the assumption that $m=1$ (see equation (5)). The hyper-parameters $\theta_{i},\ i=1,\ldots,m$, can be trained by employing the Negative Log Marginal Likelihood resulting from

where $\{\bm{x}^{n}_{b},\bm{u}^{n}_{b}\}$ are the (noisy) data on the boundary, $\{\bm{x}^{n-1},\bm{u}^{n-1}\}$ are artificially generated data to be explained later (see equation (16)), and

It is worth mentioning that the marginal likelihood provides a natural regularization mechanism that balances the trade-off between data fit and model complexity. This effect is known as Occam’s razor after William of Occam 1285–1349 who encouraged simplicity in explanations by the principle: “plurality should not be assumed without necessity”.

4 Posterior

In order to predict $u^{n}(x^{n}_{*})$ at a new test point $x^{n}_{*}$, we use the following conditional distribution

5 Propagating Uncertainty

However, to properly propagate the uncertainty associated with the initial data through time, one should not stop here. Since $\{\bm{x}^{n-1},\bm{u}^{n-1}\}$ are artificially generated data (see equation (16)) we have to marginalize them out by employing

Now, one can use the resulting posterior distribution (14) to obtain the artificially generated data $\{\bm{x}^{n},\bm{u}^{n}\}$ for the next time step with

Here, $\bm{\mu}^{n}=\mu^{n}(\bm{x}^{n})$ and $\bm{\Sigma}^{n,n}=\Sigma^{n,n}(\bm{x}^{n},\bm{x}^{n})$.

6 Example: Burgers’ equation (Backward Euler)

Burgers’ equation is a fundamental partial differential equation arising in various areas of applied mathematics, including fluid mechanics, nonlinear acoustics, gas dynamics, and traffic flow . In one space dimension the equation reads as

along with Dirichlet boundary conditions $u(t,-1)=u(t,1)=0$, where $u(t,x)$ denotes the unknown solution and $\nu$ is a viscosity parameter. Let us assume that all we observe are noisy measurements $\{\bm{x}^{0},\bm{u}^{0}\}$ of the black-box initial function $u(0,x)=-\sin(\pi x)$. Given such measurements, we would like to solve the Burgers’ equation (17) while propagating through time the uncertainty associate with the noisy initial data (see figure 1).

This example is important because it involves solving a non-linear partial differential equation. To illustrate how one can encode the structure of the physical laws expressed by Burgers’ equation in a numerical Gaussian process let us apply the backward Euler scheme to equation (17). This can be written as

We would like to place a Gaussian process prior on $u^{n}$. However, the nonlinear term $u^{n}\frac{d}{dx}u^{n}$ is causing problems simply because the product of two Gaussian processes is no longer Gaussian. Hence, we will approximate the nonlinear term with $\mu^{n-1}\frac{d}{dx}u^{n}$, where $\mu^{n-1}$ is the posterior mean of the previous time step. Therefore, the backward Euler scheme (18) can be approximated by

is a Gaussian process with a neural network covariance function

where $\theta=\left(\sigma^{2}_{0},\sigma^{2}\right)$ denotes the hyper-parameters. Here we have chosen a non-stationary prior motivated by the fact that the solution to the Burgers’ equation can develop discontinuities for small values of the viscosity parameter $\nu$. This enables us to obtain the following Numerical Gaussian Process

with covariance functions $k^{n,n}_{u,u}$, $k^{n,n-1}_{u,u}$, and $k^{n-1,n-1}_{u,u}$ given in section 5.1 of the appendix. Training, prediction, and propagating the uncertainty associated with the noisy initial observations can be performed as in sections 2.3, 2.4, and 2.5, respectively. Figure 1 depicts the noisy initial data along with the posterior distribution of the solution to the Burgers’ equation (17) at different time snapshots. It is remarkable that the proposed methodology can effectively propagate an infinite collection of correlated Gaussian random variables (i.e., a Gaussian process) through the complex nonlinear dynamics of the Burgers’ equation.

6.2 Numerical Study

It must be re-emphasized that numerical Gaussian processes, by construction, are designed to deal with cases where: (1) all we observe is noisy data on black-box initial conditions, and (2) we are interested in quantifying the uncertainty associated with such noisy data in our solutions to time-dependent partial differential equations. In fact, we recommend resorting to other alternative classical numerical methods such as Finite Differences, Finite Elements, and Spectral methods in cases where: (1) the initial function is not a black-box function and we have access to noiseless data, or (2) we are not interested in quantifying the uncertainty in our solutions. However, in order to be able to perform a systematic numerical study of the proposed methodology and despite the fact that this defeats the whole purpose of the current work, sometimes we will operate under the assumption that we have access to noiseless initial data. For instance, concerning the Burgers’ equation, if we had access to such noiseless data, we would obtain results similar to the ones reported in figure 2.

Moreover, in order to make sure that the numerical Gaussian process resulting from the backward Euler scheme (20) applied to the Burgers’ equation is indeed first-order accurate in time, we perform the numerical experiments reported in figures 5 and 5. Specifically, in figure 5 we report the time-evolution of the relative spatial $\mathcal{L}^{2}$ error until the final integration time $T=1.0$. We observe that the error indeed grows as $\mathcal{O}(\Delta{t})$, and its resulting behavior reveals both the shock development region as well as the energy dissipation due to diffusion at later times. Moreover, in figure 5 we fix the final integration time to $T=0.1$ and the number of initial and artificial data to 50, and vary the time-step size $\Delta{t}$ from $10^{-1}$ to $10^{-4}$. As expected, we recover the first-order convergence properties of the backward Euler scheme, except for a saturation region arising when we further reduce the time-step size below approximately $10^{-3}$. This behavior is not a result of the time stepping scheme but is attributed to the underlying Gaussian process regression and the finite number of spatial data points used for training and prediction. To investigate the accuracy of the posterior mean in predicting the solution as the number of training points is increased, we perform the numerical experiment reported in figure 5. Here we have considered two cases for which we fix the time step size to $\Delta{t}=10^{-2}$ and $\Delta{t}=10^{-3}$, respectively, and increase the number of initial as well as artificial data points. A similar accuracy saturation is also observed here as the number of training points is increased. In this case, this is attributed to the error accumulated due to time-stepping with the relatively large time step sizes for the first-order accurate Euler scheme. If we further keep decreasing the time-step, this saturation behavior will occur for higher numbers of total training points. The key point here is that although Gaussian processes can yield satisfactory accuracy, they, by construction, cannot force the approximation error down to machine precision. This is due to the fact that Gaussian processes are suitable for solving regression problems. This is exactly the reason why we recommend other alternative classical numerical methods for solving partial differential equations in cases where one has access to noiseless data. In such cases, it is desirable to use numerical schemes that are capable of performing exact interpolation on the data, rather than just a mere regression.

7 Example: Wave Equation (Trapezoidal Rule)

The wave equation is an important second-order linear partial differential equation for the description of wave propagation phenomena, including sound waves, light waves, and water waves. It arises in many scientific fields such as acoustics, electromagnetics, and fluid dynamics. In one space dimension the wave equation reads as

The function $u(t,x)=\frac{1}{2}\sin(\pi x)\cos(\pi t)+\frac{1}{3}\sin(3\pi x)\sin(3\pi t)$ solves this equation and satisfies the following initial and homogeneous Dirichlet boundary conditions

Now, let us assume that all we observe are noisy measurements $\{\bm{x}^{0}_{u},\bm{u}^{0}\}$ and $\{\bm{x}^{0}_{v},\bm{v}^{0}\}$ of the black-box initial functions $u^{0}$ and $v^{0}$, respectively. Given this data, we are interested in solving the wave equation (23) and quantifying the uncertainty in our solution associated with the noisy initial data (see figure 6).

To proceed, let us define $v:=u_{t}$ and rewrite the wave equation as a system of equations given by

This example is important because it involves solving a system of partial differential equations. One could rewrite the system of equations (25) in matrix-vector notations and obtain

This form is now amenable to the previous analysis provided for general linear multi-step methods. However, for pedagogical purposes, let us walk slowly through the trapezoidal rule and apply it to the system of equations (25). This can be written as

Now, let us define $u^{n-1/2}$ and $v^{n-1/2}$ to be given by

As outlined in section 2 this is a key step in the proposed methodology as it hints at the proper location to place the Gaussian process prior. Shifting the terms involved in the above equations by $-1/2$ and $+1/2$ we obtain

respectively. Now we can proceed with encoding the structure of the wave equation into a numerical Gaussian process prior for performing Bayesian machine learning of the solution $\{u(t,x),v(t,x)\}$ at any $t>0$.

are two independent Gaussian processes with squared exponential covariance functions

where $\theta_{u}=\left(\gamma_{u}^{2},w_{u}\right)$ and $\theta_{v}=\left(\gamma_{v}^{2},w_{v}\right)$. From a theoretical point of view, each covariance function gives rise to a Reproducing Kernel Hilbert space that defines a class of functions that can be represented by this kernel. In particular, the squared exponential covariance function chosen above implies smooth approximations. More complex function classes can be accommodated by appropriately choosing kernels (see e.g., equation (22)). This enables us to obtain the following numerical Gaussian process

which captures the entire structure of the trapezoidal rule (26), applied to the wave equation (23), in its covariance functions given in section 5.2 of the appendix. Training, prediction, and propagating the uncertainty associated with the noisy initial observations can be performed as in section 5.2 of the appendix. Figure 6 depicts the noisy initial data along with the posterior distribution (60) of the solution to the wave equation (23) at different time snapshots.

7.2 Numerical Study

In the case where we have access to noiseless initial data we obtain the results depicted in figure 7.

Moreover, we perform a numerical study similar to the one reported in section 2.6.2. This is to verify that the numerical Gaussian process resulting from the trapezoidal rule (26) applied to the wave equation is indeed second-order accurate in time. In particular, the numerical experiment shown in figure 10 illustrates the time evolution of the relative spatial $\mathcal{L}^{2}$ until the final integration time $T=1.5$. The second-order convergence of the algorithm is also demonstrated in figure 10 where we have fixed the number of noiseless initial and artificially generated data, while decreasing the time step size. We also investigate the convergence behavior of the algorithm for a fixed time-step $\Delta{t}=10^{-2}$ and as the number of training points is increased. The results are summarized in figure 10. The analysis of both temporal and spatial convergence properties yield qualitatively similar conclusions to the ones reported in section 2.6.2. One thing worth mentioning here is that the error in $u$ is not always less than the error in $v$ (as seen in figures 7 and 10). This just happens to be the case at time $T=0.2$.

Runge-Kutta Methods

Let us now focus on the general form of Runge-Kutta methods with $q$ stages applied to equation (1); i.e.,

Here, $u^{n+\tau_{i}}(x)=u(t^{n}+\tau_{i}\Delta t,x)$. This general form encapsulates both implicit and explicit time-stepping schemes, depending on the choice of the weights $\{a_{ij},b_{i}\}$. An important feature of the proposed methodology is that it is oblivious to the choice of these parameters, hence the implicit or explicit nature of the time-stepping scheme is ultimately irrelevant. This is in sharp contrast to classical numerical methods in which implicit time-integration is burdensome due to the need for repeatedly solving linear or nonlinear systems. Here, for a fixed number of stages $q$, the cost of performing implicit or explicit time-marching is identical. This is attributed to the fact that the structure of the time-stepping scheme is encoded in the numerical Gaussian process prior, and the algorithm only involves solving a sequence of regression problems as outlined in section 2.2. This allows us to enjoy the favorable stability properties of fully implicit schemes at no extra cost, and thus perform long-time integration using very large time-steps. Equations (33) can be equivalently written as

are $q+1$ mutually independent Gaussian processes. Therefore, we can write the joint distribution of $u^{n+1},u^{n+\tau_{q}},\ldots,u^{n+\tau_{1}},u^{n}_{q+1},\ldots,u^{n}_{1}$ which will capture the entire structure of the Runge-Kutta methods in the resulting numerical Gaussian process. However, rather than getting bogged down into heavy notation, and without sacrificing any generality, we will present the main ideas through the lens of an example.

We have chosen this classical pedagogical example as a prototype benchmark problem for testing the limits of long-time integration. This example also highlights the implementation of periodic constraints at the domain boundaries (3.1). The advection equation in one space dimension takes the form

The function $u(t,x)=\sin(2\pi(x-t))$ solves this equation and satisfies the following initial and periodic boundary conditions

However, let us assume that all we observe are noisy measurements $\{\bm{x}^{0},\bm{u}^{0}\}$ of the black-box initial function $u^{0}$. Given this data, we are interested in encoding the structure of the advection operator in a numerical Gaussian process prior and use it to infer the solution $u(t,x)$ with quantified uncertainty for any $t>0$ (see figure 11).

Let us apply the Gauss-Legendre time-stepping quadrature with two stages (thus fourth-order accurate) to the advection equation (36). Referring to equations (34), we obtain

Here, $\tau_{1}=\frac{1}{2}-\frac{1}{6}\sqrt{3}$, $\tau_{2}=\frac{1}{2}+\frac{1}{6}\sqrt{3}$, $b_{1}=b_{2}=\frac{1}{2}$, $a_{11}=a_{22}=\frac{1}{4}$, $a_{12}=\frac{1}{4}-\frac{1}{6}\sqrt{3}$, and $a_{21}=\frac{1}{4}+\frac{1}{6}\sqrt{3}$.

are three independent Gaussian processes with squared exponential covariance functions similar to the kernels used in equations (31). This assumption yields the following numerical Gaussian process

where the covariance functions are given in section 5.3 of the appendix.

1.2 Training

The hyper-parameters $\theta_{n+1}$, $\theta_{n+\tau_{2}}$, and $\theta_{n+\tau_{1}}$ can be trained by minimizing the Negative Log Marginal Likelihood resulting from

Here, $u^{n+1}(1)-u^{n+1}(0)=0$, $u^{n+\tau_{2}}(1)-u^{n+\tau_{2}}(0)=0$, and $u^{n+\tau_{1}}(1)-u^{n+\tau_{1}}(0)=0$ correspond to the periodic boundary condition (3.1). Moreover, $\bm{u}_{3}^{n}=\bm{u}_{2}^{n}=\bm{u}_{1}^{n}=\bm{u}^{n}$ and $\{\bm{x}^{n},\bm{u}^{n}\}$ are the artificially generated data. This last equality reveals a key feature of this Runge-Kutta numerical Gaussian process, namely the fact that it inspects the same data through the lens of different kernels. A detailed derivation of the covariance matrix $\bm{K}$ is given in section 5.3 of the appendix. Prediction and propagation of uncertainty associated with the noisy initial observations can be performed as in section 5.3 of the appendix. Figure 11 depicts the noisy initial data along with the posterior distribution (74) of the solution to the advection equation (36) at different time snapshots.

1.3 Numerical Study

In the case where we have access to noiseless initial data we obtain the results depicted in figure 12.

Moreover, in order to make sure that the numerical Gaussian process resulting from the Gauss-Legendre method (38) applied to the advection equation is indeed fourth-order accurate in time, we perform the numerical experiment reported in figures 15 and 15. The qualitative analysis of the temporal as well as the spatial convergence properties (as seen in figure 15) closely follows the conclusions drawn in section 2.6.2.

2 Example: Heat Equation (Trapezoidal Rule)

Revisiting the trapezoidal rule, equipped with the machinery introduced for the Runge-Kutta methods, we obtain an alternative numerical Gaussian process to the one proposed in section 2. We will apply the resulting scheme to the heat equation in two space dimensions; i.e.,

The function $u(t,x_{1},x_{2})=e^{-\frac{5\pi^{2}t}{4}}\sin(\pi x_{1})\sin\left(\frac{\pi x_{2}}{2}\right)$ solves this equation and satisfies the following initial and boundary conditions

Equations (42) involve Dirichlet boundary conditions while equation (43) corresponds to a Neumann-type boundary. Let us assume that all we observe are noisy measurements $\{(\bm{x}_{1}^{0},\bm{x}_{2}^{0}),\bm{u}^{0}\}$ of the black-box initial function $u(0,x_{1},x_{2})$. Given such measurements, we would like to infer the latent scalar field $u(t,x_{1},x_{2})$ (i.e., the solution to the heat equation (41)), while quantifying the uncertainty associated with the noisy initial data (see figure 16).

This example showcases the ability of the proposed methods to handle multi-dimensional spatial domains and mixed boundary conditions (see equations (42) and (43)). Let us apply the trapezoidal scheme to the heat equation (41). The trapezoidal rule for the heat equation is given by

Rearranging the terms, we can write $u^{n}_{1}:=u^{n}$ and

In other words, we are just rewriting equations (38) for the heat equation (41) with $\tau_{1}=0$, $\tau_{2}=1$, $b_{1}=b_{2}=\frac{1}{2}$, $a_{11}=a_{12}=0$, and $a_{21}=a_{22}=1/2$.

Similar to the strategy (35) adopted for the Runge-Kutta methods, and as an alternative to the scheme used in section 2, we make the following prior assumptions:

Here, we employ anisotropic squared exponential covariance functions of the form

The hyper-parameters are given by $\theta_{n+1}=(\gamma_{n+1}^{2},w_{n+1,1},w_{n+1,2})$ and $\theta_{n}=(\gamma_{n}^{2},w_{n,1},w_{n,2})$. To deal with the mixed boundary conditions (42) and (43), let us define $v^{n+1}:=\frac{d}{dx_{2}}u^{n+1}$ and $v^{n}:=\frac{d}{dx_{2}}u^{n}$. We obtain the following numerical Gaussian process

where the covariance functions are given in section 5.4 of the appendix.

2.2 Training

The hyper-parameters $\theta_{n+1}$ and $\theta_{n}$ can be trained by minimizing the Negative Log Marginal Likelihood resulting from

where $\{(\bm{x}_{1,D}^{n+1},\bm{x}_{2,D}^{n+1}),\bm{u}^{n+1}_{D}\}$ and $\{(\bm{x}_{1,D}^{n},\bm{x}_{2,D}^{n}),\bm{u}^{n}_{D}\}$ denote the data on the Dirichlet (42) portion of the boundary, while

correspond to the Neumann (43) boundary data. Moreover, $\bm{u}^{n}_{1}=\bm{u}^{n}_{3}=\bm{u}^{n}$ and $\{(\bm{x}^{n}_{1},\bm{x}^{n}_{2}),\bm{u}^{n}\}$ are the artificially generated data. The exact form of the covariance matrix $\bm{K}$ is given in section 5.4 of the appendix. Prediction and propagation of uncertainty associated with the noisy initial observations can be performed as in section 5.4 of the appendix. Figure 16 depicts the noisy initial data along with the posterior distribution (87) of the solution to the Heat equation (41) at different time snapshots.

2.3 Numerical Study

In order to be able to perform a systematic numerical study of the proposed methodology, we will operate under the assumption that we have access to noiseless initial data. The corresponding results are reported in figure 17.

Moreover, in order to make sure that the numerical Gaussian process resulting from the Runge-Kutta version of the trapezoidal rule (3.2) applied to the Heat equation is indeed second-order accurate in time, we perform the numerical experiments reported in figures 20 and 20. Again, the qualitative analysis of the temporal as well as the spatial convergence properties (as seen in figure 20) closely follows the conclusions drawn in section 2.6.2.

Concluding Remarks

We have presented a novel machine learning framework for encoding physical laws described by partial differential equations into Gaussian process priors for nonparametric Bayesian regression. The proposed algorithms can be used to infer solutions to time-dependent and nonlinear partial differential equations, and effectively quantify and propagate uncertainty due to noisy initial or boundary data. Moreover, to the best of our knowledge, this is the first attempt to construct structured learning machines which are 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.

In contrast to classical deterministic numerical methods for solving partial differential equations (e.g., finite difference and finite-element methods), the proposed approach is by construction capable of propagating entire probability distributions in time. Although this provides a natural platform for learning from noisy data and computing under uncertainty, it comes with a non-negligible computational cost. Specifically, a limitation of this work in its present form stems from the cubic scaling with respect to the total number of training data points. In future work we plan to design more computationally efficient algorithms by exploring ideas including recursive Kalman updates and variational inference .

From a classical numerical analysis standpoint, it also becomes natural to ask questions on convergence, derivation of dispersion relations, quantification of truncation errors, comparison against classical schemes, etc. We must underline that these questions become obsolete in presence of noisy data and cannot be straightforwardly tackled using standard techniques from numerical analysis due to the probabilistic nature of the proposed work flow. In the realm of numerical Gaussian processes such questions translate into investigating theoretical concepts like prior consistency , posterior robustness , and posterior contraction rates . These define a vast territory for analysis and future developments that currently remains unexplored.

In terms of future work, we plan to leverage the proposed framework to study more complex physical systems (e.g., fluid flows via the Navier-Stokes prior), propose extensions that can accommodate parameter inference, inverse and model discovery problems , as well as incorporate probabilistic time integration schemes that allow for a natural quantification of uncertainty due to time-stepping errors .

Acknowledgements

This works received support by the DARPA EQUiPS grant N66001-15-2-4055, and the AFOSR grant 5210009.

References

References

Appendix

The covariance functions for the Burgers’ equation example are given by $k^{n,n}_{u,u}=k$,

The only non-trivial operations in the aforementioned kernel computations are the ones involving derivatives of the kernels which can be performed using any mathematical symbolic computation program like Wolfram Mathematica.

2 Wave Equation

The covariance functions for the wave equation example are given by

It is worth highlighting that the only non-trivial but straightforward operations involved in the aforementioned kernel computations are

2.2 Training

The hyper-parameters $\theta_{u}$ and $\theta_{v}$ can be trained by minimizing the Negative Log Marginal Likelihood resulting from

where $\{\bm{x}^{n}_{b},\bm{u}^{n}_{b}\}$ are the data on the boundary, $\{\bm{x}^{n-1}_{u},\bm{u}^{n-1}\}$, $\{\bm{x}^{n-1}_{v},\bm{v}^{n-1}\}$ are artificially generated data, and

Here, the data on the boundary are given by

which correspond to the Dirichlet boundary conditions (2.7).

2.3 Posterior

In order to predict $u^{n}(x^{n}_{*u})$ and $v^{n}(x^{n}_{*v})$ at new test points $x^{n}_{*u}$ and $x^{n}_{*v}$, respectively, we use the following conditional distribution

where $\bm{q}=[\bm{q}_{u}\ \bm{q}_{v}]$ and

2.4 Propagating Uncertainty

Since $\{\bm{x}^{n-1}_{u},\bm{u}^{n-1}\}$ and $\{\bm{x}^{n-1}_{v},\bm{v}^{n-1}\}$ are artificially generated data, to properly propagate the uncertainty associated with the initial data, we have to marginalize them out by employing

Now, we can use the resulting posterior distribution to obtain the artificially generated data $\{\bm{x}^{n}_{u},\bm{u}^{n}\}$ and $\{\bm{x}^{n}_{v},\bm{v}^{n}\}$ for the next time step with

3 Advection Equation

The covariance functions for the advection equation example are given by

3.2 Training

The matrix $\bm{K}$ used in the distribution (40) is given by

3.3 Posterior

In order to predict $u^{n+1}(x^{n+1}_{*})$ at a new test point $x^{n+1}_{*}$, we use

3.4 Propagating Uncertainty

To propagate the uncertainty associate with the noisy initial data through time we have to marginalize out the artificially generated data $\{\bm{x}^{n},\bm{u}^{n}\}$ by employing

Now, we can use the resulting posterior distribution (74) to obtain the artificially generated data $\{\bm{x}^{n+1},\bm{u}^{n+1}\}$ with

4 Heat equation

The covariance functions for the Heat equation are given by

4.2 Training

The matrix $\bm{K}$ used in the distribution (48) is given by

4.3 Posterior

In order to predict $u^{n+1}(x^{n+1}_{1,*},x^{n+1}_{2,*})$ at a new test point $(x^{n+1}_{1,*},x^{n+1}_{2,*})$, we use

4.4 Propagating Uncertainty

To propagate the uncertainty associate with the noisy initial data through time we have to marginalize out the artificially generated data $\{(\bm{x}^{n}_{1},\bm{x}^{n}_{2}),\bm{u}^{n}\}$ by employing

Now, we can use the resulting posterior distribution (87) to obtain the artificially generated data $\{(\bm{x}^{n+1}_{1},\bm{x}^{n+1}_{2}),\bm{u}^{n+1}\}$ with