# Differential Equations and Gaussian Processes¶

## Gaussian Process Winter School, Genova, Italy¶

### Cristian Guarnizo Lemus, Mu Niu and Neil D. Lawrence¶

Gaussian process models can deal with linear operations and retain their nature as a Gaussian process. One example of this type of linear operation is a convolution. Such convolutions are often the result of driving a differential equaition with a function or setting that function as the initial condition. In this notebook we demonstrate the combination of Gaussian processes with differential equations, an approach that we refer to as latent force models. We consider a second order dynamical system and a spatial differential equation

## Simulation Study¶

First, we generate some data points from a second order differential equations. We start by importing some libraries

In [1]:
import numpy as np
from scipy import integrate
from IPython.display import display
import matplotlib.pyplot as plt


Now we define a function that allow us to simulate a dynamical system

In [2]:
def calc_deri(yvec, time, damper, spring, s):
return (yvec[1], 1.*s - damper*yvec[1] - spring*yvec[0])


Using the previous function, We simulate the dynamical system regarding the following equation: $$\frac{\text{d}^2 y_d(t)}{\text{d} t^2} + C_d \frac{\text{d} y_d(t)}{\text{d} t} + B_d y_d(t) = \sum_{q=1}^{Q} S_{d,q}u_q(t)$$ In the function coded before, u(t) is the unit step.

In [3]:
# System parameters:
# Spring constants
B=[1.,3.]
# Damper Constants
C=[4.,1.]
# Sensitivities
S=[[1.],[1.]]
# Numer of inputs
q=1
# Number of outputs
d=2
# Initialization of output values
Nd = 100
ti = 0.
tf = 10.
t = np.linspace(ti, tf, Nd)
start = 0
Y = np.empty([d*Nd, 1])
T = np.empty([d*Nd, 1])
index = np.empty(Y.shape, dtype = np.int16)
for i in range(d):
T[start:start+Nd, 0] = t
yarr = integrate.odeint(calc_deri, (0, 0), t, args=(C[i], B[i], S[i][0]))
Y[start:start+Nd, 0] = yarr[:, 0] + 0.01*np.random.randn(yarr.shape[0])
index[start:start+Nd, 0] = i*np.ones((Nd,), dtype = np.int16)
start = start + Nd

X = np.hstack((T, index))


Using the minimum and maximum value of time inputs we generate the inducing inputs (here we make use of a trick to differentiate between ouput and latent time inputs)

In [4]:
Mq = 50
indexz = np.empty([Mq*q, 1], dtype = np.int16)
z = np.empty([Mq*q, 1])
start = 0
zi = np.linspace(ti, tf, Mq)
for i in range(q):
z[start:start+Mq] = zi.reshape((zi.size, 1))
indexz[start:start+Mq] = i*np.ones((zi.size, 1), dtype = np.int16)
start = start + Mq
#Index number greater than the number of outputs are considered as latent time inputs
indexz += d
Z = np.hstack((z, indexz))


## Model identification using GPy's sparse inference¶

Kernel is created using the second order differential equation covariance function. Then, the model is trained using GPy's SparseGPRegression model.

In [5]:
import GPy
kern = GPy.kern.EQ_ODE2(2, output_dim=d, rank=q)
model = GPy.models.SparseGPRegression(X=X, Y=Y, Z=Z, kernel=kern)
model.inducing_inputs.fix()
model.optimize('scg', max_iters = 100)


Now we review the outcome from the sparse inference.

In [6]:
display(model)

Name                 : sparse gp mpi
Log-likelihood       : 586.791497527
Number of Parameters : 108
Parameters:
sparse_gp_mpi.           |        Value        |  Constraint  |  Prior  |  Tied to
inducing inputs          |            (50, 2)  |    fixed     |         |
eq_ode2.lengthscale      |      1.61296828862  |     +ve      |         |
eq_ode2.C                |               (2,)  |     +ve      |         |
eq_ode2.B                |               (2,)  |     +ve      |         |
eq_ode2.W                |             (2, 1)  |              |         |
Gaussian_noise.variance  |  0.000112833006706  |     +ve      |         |


In [7]:
display(model.kern.B)
display(model.kern.C)
display(model.kern.W)

  Index  |  sparse_gp_mpi.eq_ode2.B  |  Constraint  |   Prior   |  Tied to
[0]   |               0.53967786  |     +ve      |           |    N/A
[1]   |                2.3962038  |     +ve      |           |    N/A       Index  |  sparse_gp_mpi.eq_ode2.C  |  Constraint  |   Prior   |  Tied to
[0]   |                2.0954422  |     +ve      |           |    N/A
[1]   |                1.0402882  |     +ve      |           |    N/A       Index  |  sparse_gp_mpi.eq_ode2.W  |  Constraint  |   Prior   |  Tied to
[0 0]  |               0.35545259  |              |           |    N/A
[1 0]  |               0.53738722  |              |           |    N/A


Finally, we perform predictions from the training data, and compare them with the real output values.

In [8]:
Yp_mean = m.predict(X)
Yp = Yp_mean[0][:, 0]

In [9]:
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
plt.plot(X[:100, 0], Yp[:100], X[:100, 0], Y[:100, 0])

Out[9]:
[<matplotlib.lines.Line2D at 0x7f914416b210>,
<matplotlib.lines.Line2D at 0x7f914416b490>]
In [10]:
plt.plot(X[100:, 0], Yp[100:], X[100:, 0], Y[100:, 0])

Out[10]:
[<matplotlib.lines.Line2D at 0x7f914409d6d0>,
<matplotlib.lines.Line2D at 0x7f914409d950>]