Inferring the unknown parameters in Differential Equation by Gaussian Process Regression with Constraint

Specifically, the solution $u({\bm{x}})$ is assumed as a Gaussian process with zero mean and the covariance function $k_{uu}({\bm{x}},{\bm{x}}^{\prime},\bm{\gamma})$, which is denoted as

$$u({\bm{x}})\sim\mathcal{GP}(0,k_{uu}({\bm{x}},{\bm{x}}^{\prime};\bm{\gamma})),$$

(3.1)

where, $\bm{\gamma}$ denotes the hyper-parameters of the kernel function. One property of the Gaussian process in our favor is that any linear transformation, such as differentiation and integration, of a Gaussian process is still a Gaussian process. With the assumption (3.1), we consider the case when $\mathcal{F}(\cdot;{\bm{\theta}})$ is a linear operator or a combination of different linear operators, denoted as $\mathcal{L}^{\bm{\theta}}$ with given ${\bm{\theta}}$. We introduce a random function $r({\bm{x}})$ as the residual of the linear differential equation $\mathcal{L}^{\bm{\theta}}u=0$ for convenience:

$$r({\bm{x}}):=\mathcal{L}^{\bm{\theta}}u(x).$$

Then the function $r({\bm{x}})$ is also a mean-zero Gaussian process

$$r({\bm{x}})=\mathcal{GP}(0,k_{rr}({\bm{x}},{\bm{x}}^{\prime};\bm{\gamma})),$$

where $k_{rr}({\bm{x}},{\bm{x}}^{\prime})=\text{Cov}(\mathcal{L}^{\bm{\theta}}u({\bm{x}}),\mathcal{L}^{\bm{\theta}}u({\bm{x}}^{\prime}))$ denotes the covariance function of $\mathcal{L}^{\bm{\theta}}u$ between ${\bm{x}}$ and ${\bm{x}}^{\prime}$. The following fundamental relationship between the kernels $k_{uu}$ and $k_{rr}$, $k_{ur}$, $k_{ru}$ is well-known (see e.g. [21, 9]),

$$\begin{split}k_{rr}({\bm{x}},{\bm{x}}^{\prime};\bm{\gamma})&=\mathcal{L}_{{\bm{x}}}^{\bm{\theta}}\mathcal{L}_{{\bm{x}}^{\prime}}^{\bm{\theta}}k_{uu}({\bm{x}},{\bm{x}}^{\prime};\bm{\gamma}),\\ k_{ur}({\bm{x}},{\bm{x}}^{\prime};\bm{\gamma})&=\mathcal{L}_{{\bm{x}}^{\prime}}^{\bm{\theta}}k_{uu}({\bm{x}},{\bm{x}}^{\prime};\bm{\gamma}),\\ k_{ru}({\bm{x}},{\bm{x}}^{\prime};\bm{\gamma})&=\mathcal{L}^{\bm{\theta}}_{{\bm{x}}}k_{uu}({\bm{x}},{\bm{x}}^{\prime};\bm{\gamma}),\end{split}$$

(3.2)

Based on the Gaussian assumption and the covariance expressions between $u$ and $r$ in (3.2), a joint inference framework of Gaussian process regression for the available observation data of $u$ and $r$ can be naturally constructed. By interpreting the equation $r:=\mathcal{L}^{\bm{\theta}}u=0$ as the constraint of the function $r$, GPRC will significantly improve the accuracy of estimation of solution and its derivatives due to the additional equation information.

Training a Gaussian Process

The training process is to find the optimal parameters $\bm{\beta}=\{\bm{\gamma},\sigma_{obs}^{2}\}$ by maximum likelihood estimation (MLE). Given $n$ noisy observations of the state $u$ at $n$ points $\{{\bm{x}}_{i}:1\leq i\leq n\}$, we denote the state vector ${\bm{y}}\equiv[y_{1},y_{2},\dots,y_{n}]^{T}\in{\mathbb{R}}^{n\times 1}$, the residual vector $\bm{r}\equiv[r_{1},r_{2},\dots,r_{n}]^{T}\in{\mathbb{R}}^{n\times 1}$ and the training point matrix $X\equiv[{\bm{x}}_{1},{\bm{x}}_{2},\dots,{\bm{x}}_{n}]^{T}\in{\mathbb{R}}^{n\times D}$. Let $Y=\begin{bmatrix}{\bm{y}}\\ \bm{r}\end{bmatrix}\in{\mathbb{R}}^{2n}$. Then the negative log marginal likelihood of $p(Y|\bm{\beta})$ has the following expression

$$\begin{split}-\log p(Y|\bm{\beta})=\frac{1}{2}\log(\det{K})+\frac{1}{2}Y^{T}K^{-1}Y+\frac{n}{2}\log 2\pi,\end{split}$$

(3.3)

where the $2n\times 2n$ matrix $K$ is defined by

$$K=\begin{bmatrix}&K_{uu}+\sigma^{2}_{u}I&K_{ur}\\ &K_{ru}&K_{rr}+\sigma^{2}_{r}I\end{bmatrix},$$

The matrices $K_{uu},K_{ur}$, $K_{ru}$ and $K_{rr}$ correspond to the kernel functions $k_{uu},k_{ur},k_{ru}$ and $k_{rr}$ respectively in (3.1) and (3.2), evaluated at the $n$ points $X$.
The minimisation of the negative log marginal likelihood is a non-convex optimisation task. However gradients are easily obtained

$$\begin{split}\frac{\partial L}{\partial\beta_{k}}=\frac{1}{2}tr(K\frac{\partial K^{-1}}{\partial\beta_{k}})+\frac{1}{2}Y^{T}\frac{\partial K}{\partial\beta_{k}}K^{-1}\frac{\partial K}{\partial\beta_{k}}Y,\end{split}$$

(3.4)

where $\beta_{k}$ indicates the $i$th hyper-parameter of $\bm{\beta}$ and therefore standard gradient optimisers can be used to minimize the negative log marginal likelihood (3.3) within the gradient (3.4).

Prediction After the hyper-parameters in the GPRC are computed, the prediction of the function $u({\bm{x}})$ or its derivatives of interest, denoted as a function $l({\bm{x}})$ (e.g., $l({\bm{x}})=\partial_{{\bm{x}}}u({\bm{x}})$), at a new test point ${\bm{x}}_{*}$ is described below. The covariance function of the GP $l({\bm{x}})$ is denoted by $k_{ll}({\bm{x}},{\bm{x}}^{\prime})$ by the convenience.

With a given point ${\bm{x}}_{*}$, the differential equation provides the fact $r({\bm{x}}_{*})=0$, which is a useful observation to incorporate the Bayesian inference. To enhance this condition locally, $m$ equally-spaced points around ${\bm{x}}_{*}$, which are referred as the extended set $\chi:=\{{\bm{x}}_{*}^{j}:1\leq j\leq m\}$ are actually considered in the prediction. The set $\chi$ is similar to the window/scale concept in local polynomial regression and thus adaptive strategy is possible. Therefore, the value of the equation constraints in the extended set should be zero, i.e., $\bm{r}_{\chi}=\bm{0}$. These points ${\bm{x}}_{*}^{j}$ are supposed to resolve a characteristic length for the residual process $r({\bm{x}})$.

This is the posterior distribution for a specific set of test case,

		$\displaystyle\bar{u}({\bm{x}}_{*})=K_{u_{*}\bullet}\widehat{K}_{u\bm{r}_{\chi}}^{-1}\begin{bmatrix}{\bm{y}}\\ \bm{r}_{\chi}\end{bmatrix},S_{u}({\bm{x}}_{*})=K_{uu}({\bm{x}}_{*},{\bm{x}}_{*})-K_{u_{*}\bullet}\widehat{K}_{u\bm{r}_{\chi}}^{-1}K_{u_{*}\bullet}^{T}.$		(3.5)
		$\displaystyle\bar{l}({\bm{x}}_{*})=K_{l_{*}\bullet}\widehat{K}_{u\bm{r}_{\chi}}^{-1}\begin{bmatrix}{\bm{y}}\\ \bm{r}_{\chi}\end{bmatrix},S_{l}({\bm{x}}_{*})=K_{ll}({\bm{x}}_{*},{\bm{x}}_{*})-K_{l_{*}\bullet}\widehat{K}_{u\bm{r}_{\chi}}^{-1}K_{l_{*}\bullet}^{T},$		(3.6)

where

$$\widehat{K}_{u\bm{r}_{\chi}}=\begin{bmatrix}K_{u}&K_{u\bm{r}_{\chi}}\\ K_{u\bm{r}_{\chi}}^{T}&K_{\bm{r}_{\chi}}\end{bmatrix}\in{\mathbb{R}}^{(n+m)\times(n+m)},$$

$K_{u}=K_{uu}+\sigma^{2}_{u}\bm{I}_{u}\in{\mathbb{R}}^{n\times n}$, $K_{\bm{r}_{\chi}}=K_{\bm{r}_{\chi}\bm{r}_{\chi}}+\sigma^{2}_{r}\bm{I}_{r}\in{\mathbb{R}}^{m\times m}$, $K_{u_{*}}=k_{uu}({\bm{x}}_{*},{\bm{x}}_{*})$, $K_{l_{*}}=k_{ll}({\bm{x}}_{*},{\bm{x}}_{*})$ and $K_{u_{*}\bullet}=[K_{u_{*}u},K_{u_{*}\bm{r}_{\chi}}]\in{\mathbb{R}}^{1\times(n+m)}$ and $K_{l_{*}\bullet}=[K_{l_{*}u},K_{l_{*}\bm{r}_{\chi}}]\in{\mathbb{R}}^{1\times(n+m)}$. The posterior variances $S_{u}({\bm{x}}_{*})$ and $S_{l}({\bm{x}}_{*})$ can be used as good indicators of how confident the predictions are.

For nonlinear differential equations, the equation constraint can’t be formulated as a Gaussian process since the product of Gaussian processes is not a Gaussian process anymore. Here an iterative linearization method is introduced for nonlinear equations, motivated by the Picard iteration method[7] , which is an extension of method used in Van der Pol equation example in [24].

Assume there have an initial guess of the solution $\hat{u}_{0}$ and its derivative $\hat{l}^{i}_{0}$ which can be fitted by any classical method, like Gaussian process regression and basis function expansion. Here the hat notation $\hat{\cdot}$ represents an estimation and we are only interest in certain derivative terms which is useful for lineariztion. So the $i$th derivative term $\hat{l}^{i}$ is used to represent all derivatives of interest for simplicity. Then any nonlinear term of the nonlinear equation, which is composed of $K$ solution or its derivative terms, can be simply approximated by replacing $K-1$ terms with the corresponding given initial guess of the solution $\hat{u}_{0}$ and its derivatives $\hat{l}^{i}_{0}$, and only keep one term. This is a quite straightforward linearlization strategy and can be easily applied. There is a trick that replacing the solution or low order derivatives, would be preferred than the one who replacing the high order derivatives of the nonlinear term in the approximation. That is because the instability of numerical differentiation of the fitted function solely from the data which is the measurement of solution with adding noise, will be more serious in estimating higher order derivative. The linearization of $\mathcal{F}(u(\bm{x});{\bm{\theta}})$ is expressed as

$$\mathcal{F}(u,l^{i},\dots;\hat{u}_{0},\hat{l}^{i}_{0},\dots;{\bm{\theta}})=0.$$

(3.7)

Based on linearization in the Eq. (3.7), GPRC introduced in Section 3.1 can be used to estimate a new solution $\hat{u}$ and derivative $\hat{l}^{i}$. Then an iterative scheme can be constructed with the new pair $\{\hat{u},\hat{l}^{i}\}$. Roughly speaking in the $n$th cycle, given the current guess of solution $\hat{u}_{n}$ and derivatives $\hat{l}^{i}_{n}$, we can compute a new (and possibly better) approximation $\{\hat{u}_{n+1},\hat{l}^{i}_{n+1}\}$ using

$$\{\hat{u}_{n+1},\hat{l}^{i}_{n+1}\}=\textbf{GPRC}(\mathcal{F}(u,l^{i},\dots;\hat{u}_{n},\hat{l}^{i}_{n},\dots;{\bm{\theta}}),\mathbf{y}),$$

(3.8)

where $\textbf{GPRC}(\cdot)$ indicates the estimator of GPRC whose output is the posterior mean functions of the solution and its derivatives. Our stopping criterion is the sum of RMSE values of solution and the RMSE values of all derivatives of interest, which is expressed by

$$e=\int(||\hat{u}_{n+1}-\hat{u}_{n}||^{2}+\sum_{i}||\hat{l}^{i}_{n+1}-\hat{l}^{i}_{n}||^{2})p({\bm{x}})d{\bm{x}},$$

where $||\cdot||^{2}$ denotes the L2 normalization. If $e$ is small than a give small positive threshold $\epsilon_{rmse}$, then we think the iteration reaches stability.

First we consider the following $1$-dimension linear ordinary differential equation,

$$u^{\prime\prime}(t)+\theta_{1}u^{\prime}(t)+\theta_{2}u(t)=0.$$

(5.1)

where $u$ is a function of time $t$, ^′ and ^′′ denotes first and second order derivative operators, respectively. ${\bm{\theta}}=[\theta_{1},\theta_{2}]$ is the unknown parameters which need to refer. Our true parameter values are taken as $\theta^{*}_{1}=1$,$\theta^{*}_{2}=3$. Then $21$ observations are evenly measured in the time interval $[0,10]$ with Gaussian noise $\mathcal{N}(0,\sigma^{2}_{obs}=0.1)$. The observations and the true solution $u$ are shown in Fig. 1.

We assume the prior of the sate function $u$ is a zero-mean Gaussian process expressed in (3.1). As discussed in Section 3, the equation constraint $r=u^{\prime\prime}+\theta_{1}u^{\prime}+\theta_{2}u$ is also a Gaussian process, $r\sim\mathcal{GP}(0,k_{rr}({\bm{x}},{\bm{x}}^{\prime};\bm{\gamma})).$ With the property of covariance, the covariance between $r$ and $u$ can be expanded as

	$\displaystyle\mathrm{Cov}(r,u)$	$\displaystyle=\mathrm{Cov}(u^{\prime\prime}+\theta_{1}u^{\prime}+\theta_{2}u,u)$
		$\displaystyle=\mathrm{Cov}(u^{\prime\prime},u)+\theta_{1}\mathrm{Cov}(u^{\prime},u)+\theta_{2}\mathrm{Cov}(u,u).$

So, $k_{ru}=k_{u^{\prime\prime}u}+\theta_{1}k_{u^{\prime}u}+\theta_{2}k_{u,u}$. Similarly, the kernel functions corresponding to covariances $\mathrm{Cov}(r,u^{\prime}),\mathrm{Cov}(r,u^{\prime\prime})$ and $\mathrm{Cov}(r,r)$ are expressed as

	$\displaystyle k_{ru^{\prime}}$	$\displaystyle=k_{u^{\prime\prime}u^{\prime}}+\theta_{1}k_{u^{\prime}u^{\prime}}+\theta_{2}k_{u,u^{\prime}},$
	$\displaystyle k_{ru^{\prime\prime}}$	$\displaystyle=k_{u^{\prime\prime}u^{\prime\prime}}+\theta_{1}k_{u^{\prime}u^{\prime\prime}}+\theta_{2}k_{u,u^{\prime\prime}},$
	$\displaystyle k_{rr}$	$\displaystyle=k_{u^{\prime\prime}u^{\prime\prime}}+{\theta_{1}}^{2}k_{u^{\prime}u^{\prime}}+{\theta_{2}}^{2}k_{uu}+2{\theta_{1}}k_{u^{\prime\prime}u^{\prime}}+2{\theta_{2}}k_{u^{\prime\prime}u}+2{\theta_{1}}{\theta_{2}}k_{u^{\prime}u}.$

By equation (3.5) and (3.6), the posterior distributions $p(u_{*}|{\bm{y}},\bm{r}_{\chi})={\mathcal{N}}(m_{u},\Sigma_{u})$, $p(u_{*}^{\prime}|{\bm{y}},\bm{r}_{\chi})={\mathcal{N}}(m_{u^{\prime}},\Sigma_{u^{\prime}})$ and $p(u_{*}^{\prime\prime}|{\bm{y}},\bm{r}_{\chi})={\mathcal{N}}(m_{u^{\prime\prime}},\Sigma_{u^{\prime\prime}})$ are given below

	$\displaystyle m_{u_{*}}$	$\displaystyle=K_{\bullet u_{*}}^{T}\widehat{K}_{u\bm{r}_{\chi}}^{-1}Y,\quad\Sigma_{u_{*}}=K_{u_{*}}-K_{\bullet u_{*}}^{T}\widehat{K}_{u\bm{r}_{\chi}}^{-1}K_{\bullet u_{*}}$
	$\displaystyle m_{u^{\prime}_{*}}$	$\displaystyle=K_{\bullet u_{*}^{\prime}}^{T}\widehat{K}_{u\bm{r}_{\chi}}^{-1}Y,\quad\Sigma_{u^{\prime}_{*}}=K_{u^{\prime}_{*}}-K_{\bullet u_{*}^{\prime}}^{T}\widehat{K}_{u\bm{r}_{\chi}}^{-1}K_{\bullet u_{*}^{\prime}}$
	$\displaystyle m_{u^{\prime\prime}_{*}}$	$\displaystyle=K_{\bullet u_{*}^{\prime\prime}}^{T}\widehat{K}_{u\bm{r}_{\chi}}^{-1}Y,\quad\Sigma_{u^{\prime\prime}_{*}}=K_{u^{\prime\prime}_{*}}-K_{\bullet u_{*}^{\prime\prime}}^{T}\widehat{K}_{u\bm{r}_{\chi}}^{-1}K_{\bullet u_{*}^{\prime\prime}},$

where $Y=[{\bm{y}},\bm{r}_{\chi}]^{T}$. This is the posterior estimation of state and derivative functions. In GPRC algorithm, we set $\sigma_{r}^{2}=10$ and $\chi=[t_{*}]$, which means only thee predicted point is included in the extend set $\chi$. Although we would sacrifice some precision of prediction due to the small extend set, this setting can efficiently speed up the calculation. Then we can obtain the posterior estimation of solution and its derivatives within any given ${\bm{\theta}}$. A Gaussian proposal distribution $q({\bm{\theta}}^{*}|{\bm{\theta}}^{(i)})=\mathcal{N}({\bm{\theta}}^{(i)},0.6^{2})$ is used as the proposal distribution in our MH sampler with $\alpha=100$ in our likelihood function (4.1). $5000$ samples are drawn from the posterior. As a comparison same number of samples are obtained by two-stage method and the both posterior results are shown in Figure 2. Obviously from the figure on the left, we can find that the samples (red scattered points) computed by our method have a wider variance, which means a bigger uncertainty of the posterior is estimated. In the marginal pdf figures, the first dimension marginal pdf value at $\theta_{1}^{*}$ computed by our method (red line) are larger than one computed by two-stage method (blue dash) and the second dimension marginal pdf values computed by two method at $\theta_{1}^{*}$ are equal. These results indicate the posterior estimation computed by our method is more credible. This point would be better clear in a small-noise case because the bias caused by the numerical error of derivatives estimation would be highlights when the measure noise is small.

We then consider the small-noise case where we use the same implementation configurations as in the large-noise case. We change the noise level from $0.1$ to $0.05$ and collect the same number of observations at same locations. The samples from posterior are shown in Figure 2 (right). Similar to the large-noise case, the figure shows that the results of the two-stage are of low accuracy, while our method yield rather good results.

The Van der Pol model is described by 1-dimensional nonlinear ODE equation,

$$u^{\prime\prime}-\mu(1-u^{2})u^{\prime}+u=0,$$

(5.2)

where $\mu$ is the unknown parameter. For this nonlinear differential equation, we first use the GPR method to obtain initial guess of the solution $\hat{u}_{0}$ only from data. Then we can obtain a linearizaiton form of (5.2), by using the approximation $\hat{u}_{0}$ to replace the original $u$ of nonlinear term, expressed as

$$u^{\prime\prime}+\mu(1-\hat{u}_{0}^{2})u^{\prime}+u=0.$$

In this example 41 observations are measured at every $0.5$ time units on the interval [0,20] with given true parameter $\mu^{*}=0.5$. The clean data are corrupted with additive Gaussian noise with variance $\sigma_{obs}^{2}=0.01$. The figure on the left in Fig 3 shows the true solution $u$ and the observations. In GPRC, $\sigma_{r}^{2}=0.1$ is set and $10$ points which are evenly selected in the domain $[t_{*}-1,t_{*}+1]$, form the extend set $\chi$ . In the MH sampling procedure, the design of the proposal distribution $q(\theta^{*}|\theta^{(i)})$ is $\mathcal{N}(\theta^{(i)},0.6^{2})$. The histogram of the posterior samples computed by both our method and two-stage method are shown in the Fig. 3 (right). In the figure, the true parameter (black line) are closer to the mean of red histogram of samples computed by our method than the one of blue histogram, computed by two-stage method.

The KdV equation is a nonlinear, dispersive partial differential equation for function $u$ of two real variables, space $x\in\mathbb{R}^{1\times 1}$ and time $t$:

$$\frac{\partial u}{\partial t}+\theta_{1}u\frac{\partial u}{\partial x}+\theta_{2}\frac{\partial^{3}u}{\partial x^{3}}=0,$$

It is an asymptotic simplification of Euler equations used to model waves in shallow water. It can also be viewed as Burger’s equation with added dispersive term. In this equation two parameters need to be estimated from some noisy observations of function $u$. The true solution is created using a spectral method with $512$ spatial points and $200$ timesteps. More details about this model refer [20]. This is a highly non-linear PDE. So $200$ solution values are randomly selected from the given discrete solution in the time domain $[15,20]$ and spatial domain $[-21,26]$, within our true parameters are $\theta^{*}_{1}=6$ and $\theta^{*}_{2}=1$. A Gaussian noise is added on the selected solution values with zero mean and $0.1$ variance. The true solution and the observations are shown in the Fig. 4 (left).

In MH sampler procedure, $5000$ samples are obtained by the given proposal distribution $\mathcal{N}(0,0.6I)$. The posterior samples and marginal distributions are shown in right figure of Fig. 4. As shown, the true parameter (white star) is closer to the center of the posterior samples computed by our method than the one of two-stage method. In addition, at the true parameters $\theta^{*}_{1}$ and $\theta^{*}_{2}$, the values of marginal PDF computed by our method are larger than the one computed by two-stage method, express that our method yields a better posterior estimation.