pyrobolearn/tutorials/machine_learning/GP.ipynb

Gaussian Process

In this tutorial, we expose what gaussian processes are, and how to use the GPy library. We first provide a gentle reminder about Gaussian distributions and their properties.

In [ ]:

# Import all the important libraries
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.stats import multivariate_normal
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.mlab as mlab
from ipywidgets import widgets as wg
from matplotlib import cm

import GPy
#%matplotlib inline
%matplotlib notebook

1D Gaussian distribution

In [ ]:

# Plot
def plot_gaussian(mu=0, sigma=1):
    x = np.linspace(-3, 3, 100)
    plt.plot(x, norm.pdf(x, mu, sigma))
    plt.xlabel('x')
    plt.ylabel('p(x)')

wg.interact(plot_gaussian, mu=(-2,2,0.1), sigma=(-2,2,0.1))
plt.show()

Multivariate Gaussian distribution (2D)

The multivariable Gaussian distribution is a generalization of the Gaussian distribution to vectors. See wikipedia for more info.

In [ ]:

# moments
mu = np.array([0,0])
Sigma = np.array([[1,0], 
                  [0,1]])
Sigma1 = np.array([[1,0.5],
                   [0.5,1]])
Sigma2 = np.array([[1,-0.5],
                   [-0.5,1]])
Sigmas = [Sigma, Sigma1, Sigma2]

pts = []
for S in Sigmas:
    pts.append(np.random.multivariate_normal(mu, S, 1000).T)

# Plotting
width = 16
height = 4
plt.figure(figsize=(width, height))

# make plot
for i in range(len(Sigmas)):
    plt.subplot(1,3,i+1)
    plt.title('Plot '+str(i+1))
    plt.ylim(-4,4)
    plt.xlim(-4,4)
    plt.xlabel('x1')
    plt.ylabel('x2')
    plt.plot(pts[i][0], pts[i][1],'o')
    #plt.scatter(pts[i][0], pts[i][1])
    
plt.show()

The 1st plot above is described by: \begin{equation} \left[ \begin{array}{c} x_1\x_2 \end{array}\right] \sim \mathcal{N} \left(\left[ \begin{array}{c} 0\0 \end{array}\right], \left[ \begin{array}{cc} 1 & 0\ 0 & 1 \end{array}\right] \right) \end{equation}

The 2nd plot is given by: \begin{equation} \left[ \begin{array}{c} x_1\x_2 \end{array}\right] \sim \mathcal{N}\left(\left[ \begin{array}{c} 0\0 \end{array}\right], \left[ \begin{array}{cc} 1 & 0.5\ 0.5 & 1 \end{array}\right]\right) \end{equation}

Finally, the 3rd plot is given by: \begin{equation} \left[ \begin{array}{c} x_1\x_2 \end{array}\right] \sim \mathcal{N}\left(\left[ \begin{array}{c} 0\0 \end{array}\right], \left[ \begin{array}{cc} 1 & -0.5\ -0.5 & 1 \end{array}\right]\right) \end{equation}

The covariance (and the dot product) measures the similarity.

For the 2nd and 3rd plots, x_1 is correlated with x_2, i.e. knowing x_1 gives us information about x_2.

Joint distribution p(x_1,x_2)

The joint distribution p(x_1, x_2) is given by:

\begin{equation} \left[ \begin{array}{c} x_1\x_2 \end{array}\right] \sim \mathcal{N}\left(\left[ \begin{array}{c} \mu_1 \ \mu_2 \end{array}\right], \left[ \begin{array}{cc} \Sigma_{11} & \Sigma_{12} \ \Sigma_{21} & \Sigma_{22} \end{array}\right] \right) = \mathcal{N}(\pmb{\mu}, \pmb{\Sigma}) \end{equation}

In [ ]:

# Reference: http://stackoverflow.com/questions/38698277/plot-normal-distribution-in-3d

# moments
mu = np.array([0,0])
Sigma = np.array([[1,0], [0,1]])

# Create grid and multivariate normal
step = 500
bound = 10
x = np.linspace(-bound,bound,step)
y = np.linspace(-bound,bound,step)
X, Y = np.meshgrid(x,y)
pos = np.empty(X.shape + (2,))
pos[:, :, 0] = X; pos[:, :, 1] = Y
pdf = multivariate_normal(mu, Sigma).pdf(pos)

# Plot
fig = plt.figure(figsize=plt.figaspect(0.5)) # Twice as wide as it is tall.

# 1st subplot (3D)
ax = fig.add_subplot(1, 2, 1, projection='3d')
ax.plot_surface(X, Y, pdf, cmap='viridis', linewidth=0)
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_zlabel('p(x1, x2)')

# 2nd subplot (2D)
ax = fig.add_subplot(1, 2, 2)
ax.contourf(x, y, pdf)
#ax.colorbar()
ax.set_xlabel('x1')
ax.set_ylabel('x2')

fig.tight_layout()

plt.show()

Normalization

In order to be a valid probability distribution, the volume under the surface should equal to 1.

\begin{equation} \int \int p(x_1,x_2) dx_1 dx_2 = 1 \end{equation}

In [ ]:

# p(x1,x2) = pdf
# dx = 2.*bound/step
# dx1 dx2 = (2.*bound/step)**2
print("Summation: {}".format((2.*bound/step)**2 * pdf.sum()))

Conditional distribution p(x_2|x_1)

What is the mean \mu_{2|1} and the variance \Sigma_{2|1} of the conditional distribution p(x_2|x_1) = \mathcal{N}(\mu_{2|1}, \Sigma_{2|1})?

We know the mean \pmb{\mu} and the covariance \pmb{\Sigma} of the joint distribution p(x_1,x_2). Using the Schur complement, we obtain:

\begin{align} \mu_{2|1} &= \mu_{2} + \Sigma_{21}\Sigma_{22}^{-1}(x_2 - \mu_2) \ \Sigma_{2|1} &= \Sigma_{22} - \Sigma_{21}\Sigma_{22}^{-1}\Sigma_{12} \end{align}

For the demo, check Murphy's book "Machine Learning: A Probabilistic Perspective", section 4.3.4

In [ ]:

fig = plt.figure(figsize=plt.figaspect(0.5)) # Twice as wide as it is tall.

x1_value = 0
z_max = pdf.max()

# 1st subplot
ax = fig.add_subplot(1, 2, 1, projection='3d')
ax.plot_surface(X, Y, pdf, cmap='viridis', linewidth=0)
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_zlabel('p(x1, x2)')
y1 = np.linspace(-bound,bound,2)
z = np.linspace(0,z_max,2)
Y1, Z = np.meshgrid(y1,z)
ax.plot_surface(x1_value, Y1, Z, color='red', alpha=0.2)
#cset = ax.contourf(X, Y, pdf, zdir='x', offset=-bound, cmap=cm.coolwarm)

# 2nd subplot
ax = fig.add_subplot(1, 2, 2)
ax.plot(x, pdf[step//2 + x1_value*step//(2*bound)])
ax.set_xlabel('x2')
ax.set_ylabel('p(x2|x1)')

fig.tight_layout()

plt.show()

Marginal distribution p(x_1) and p(x_2)

\begin{align} p(x_1) &= \int p(x_1, x_2) dx_2 = \mathcal{N}(\mu_1, \Sigma_{11}) \ p(x_2) &= \int p(x_1, x_2) dx_1 = \mathcal{N}(\mu_2, \Sigma_{22}) \end{align}

In [ ]:

fig = plt.figure(figsize=plt.figaspect(0.5))
plt.subplot(1,2,1)
plt.title('By summing')
dx = 2. * bound / step
plt.plot(x, pdf.sum(0) * dx, color='blue')
plt.xlabel('x2')
plt.ylabel('p(x2)')

plt.subplot(1,2,2)
plt.title('by using the normal distribution')
plt.plot(x, norm.pdf(x, mu[1], Sigma[1,1]), color='red')
plt.xlabel('x2')
plt.ylabel('p(x2)')

fig.tight_layout()
plt.show()

Gaussian Processes (GPs)

A Gaussian process is a Gaussian distribution over functions. That is, it is a generalization of the multivariable Gaussian distribution to infinite vectors.

It will become clearer with an example.

In [ ]:

x = [0.5,0.8,1.4]
f = [1,2,6]

plt.plot(x,f,'o')
for i in range(len(x)):
    plt.annotate('f'+str(i+1), (x[i],f[i]))
plt.xlim(0,2)
plt.ylim(0,6.5)
plt.ylabel('f(x)')
plt.xlabel('x')
plt.xticks(x, ['x'+str(i+1) for i in range(len(x))])
plt.show()

\begin{align} \left[ \begin{array}{c} f_1 \ f_2 \ f_3 \end{array}\right] &\sim \mathcal{N}\left( \left[ \begin{array}{c} 0 \ 0 \ 0 \end{array}\right], \left[ \begin{array}{ccc} K_{11} & K_{12} & K_{13} \ K_{21} & K_{22} & K_{23} \ K_{31} & K_{32} & K_{33} \end{array}\right] \right) \ &\sim \mathcal{N}\left( \left[ \begin{array}{c} 0 \ 0 \ 0 \end{array}\right], \left[ \begin{array}{ccc} 1 & 0.7 & 0.2 \ 0.7 & 1 & 0.6 \ 0.2 & 0.6 & 1 \end{array}\right] \right) \end{align}

Similarity measure: K_{ij} = \exp(- ||x_i - x_j||^2) = \left\{ \begin{array}{ll} 0 & ||x_i - x_j|| \rightarrow \infty \\ 1 & x_i = x_j \end{array} \right.

Prediction (noiseless GP regression): given data \mathcal{D} = \{(x_1,f_1), (x_2,f_2), (x_3,f_3)\}, and new point x_* (e.g. $x_*$=1.4), what is the value of f_*?

\begin{equation} \pmb{f} \sim \mathcal{N}(\pmb{0}, \pmb{K}) \qquad \mbox{and} \qquad f_* \sim \mathcal{N}(0, K(x_,x_)) = \mathcal{N}(0, K_{**}) \end{equation}

In this case, K_{**} = K(x_*,x_*) = \exp(- ||x_* - x_*||^2) = 1.

Now, we can write:

\begin{equation} \left[ \begin{array}{c} \pmb{f} \ f_* \end{array} \right] \sim \mathcal{N}\left(\pmb{0}, \left[\begin{array}{cc} \left[ \begin{array}{ccc} K_{11} & K_{12} & K_{13} \ K_{21} & K_{22} & K_{23} \ K_{31} & K_{32} & K_{33} \end{array} \right] & \left[ \begin{array}{c} K_{1*} \ K_{2*} \ K_{3*} \end{array} \right] \ \left[ \begin{array}{ccc} K_{*1} & K_{2} & K_{3} \end{array} \right] & \left[\begin{array}{c} K_{**} \end{array} \right] \end{array} \right] \right) = \mathcal{N}\left(\pmb{0}, \left[\begin{array}{cc} \pmb{K} & \pmb{K}_ \ \pmb{K}_^T & \pmb{K}_{**} \end{array}\right]\right) \end{equation}

Using the formula for the conditional probability p(f_*|f), we have:

\begin{align} \mu_* &= \mathbb{E}[f_] = \pmb{K}_^T \pmb{K}^{-1}\pmb{f} \ c_* &= K_{**} - \pmb{K}*^T \pmb{K}^{-1}\pmb{K}* \end{align}

We can thus predict the mean \mu_* and the variance c_* for the test point x_*.

In [ ]:

x = [0.5,0.8,1.4]
f = [1,2,6]

x_new = 1.3
f_new = 5.2

plt.plot(x+[x_new],f+[f_new],'o')
for i in range(len(x)):
    plt.annotate('f'+str(i+1), (x[i],f[i]))
plt.errorbar(x_new, f_new, yerr=1)
plt.annotate('f*', (x_new+0.02, f_new))
plt.xlim(0,2)
plt.ylim(0,6.5)
plt.ylabel('f(x)')
plt.xlabel('x')
plt.xticks(x+[x_new], ['x'+str(i+1) for i in range(len(x))]+['x*'])
plt.show()

Generalization

A GP defines a distribution over functions p(f) (i.e. it is the joint distribution over all the infinite function values).

Definition: p(f) is a GP if for any finite subset \{x_1,...,x_n\} ⊂ X, the marginal distribution over that finite subset p(f) has a multivariate Gaussian distribution.

Prior on f: \begin{equation} \pmb{f}|\pmb{x} \sim \mathcal{GP}(\pmb{\mu}(\pmb{x}), \pmb{K}(\pmb{x}, \pmb{x})) \end{equation} with \begin{align*} \pmb{\mu}(\pmb{x}) &= \mathbb{E}_f \lbrack \pmb{x} \rbrack \ k(\pmb{x}, \pmb{x'}) &= \mathbb{E}_f \lbrack (\pmb{x} - \pmb{\mu}(\pmb{x})) (\pmb{x'} - \pmb{\mu}(\pmb{x'})) \rbrack \end{align*}

Often written as: \begin{equation} \pmb{f} \sim \mathcal{GP}(\pmb{0}, \pmb{K}) \end{equation}

Concretely, assume \pmb{x} \in \mathbb{R}^{50}, then \pmb{K}(\pmb{x}, \pmb{x}) \in \mathbb{R}^{50 \times 50}, then \pmb{f} \sim \mathcal{GP}(\pmb{0}, \pmb{K}) means:

\begin{equation} \left[ \begin{array}{c} f_1 \ \vdots \ f_{50} \end{array}\right] := \left[ \begin{array}{c} f(x_1) \ \vdots \ f(x_{50}) \end{array}\right] \sim \mathcal{N}\left( \left[ \begin{array}{c} 0 \ \vdots \ 0 \end{array}\right], \left[ \begin{array}{ccc} k(x_1,x_1) & \cdots & k(x_1, x_{50}) \ \vdots & \ddots & \vdots \ k(x_{50},x_1) & \cdots & k(x_{50}, x_{50}) \end{array} \right] \right) \end{equation}

RBF kernel

Let's choose a RBF (a.k.a Squared Exponential, Gaussian) kernel:

\begin{equation} \pmb{K} = \left[ \begin{array}{ccc} k(x_1,x_1) & \cdots & k(x_1, x_d) \ \vdots & \ddots & \vdots \ k(x_d,x_1) & \cdots & k(x_d, x_d) \end{array} \right] \end{equation} with \begin{equation} k(x_i, x_j) = \alpha^2 \exp \left( - \frac{(x_i - x_j)^2}{2l} \right) \qquad \mbox{ and hyperparameters } \pmb{\Phi} = \left{ \begin{array}{l} \alpha \mbox{: amplitude} \ l \mbox{: the lengthscale} \end{array} \right. \end{equation}

This function k is infinitely differentiable.

In [ ]:

# Reference: https://www.youtube.com/watch?v=4vGiHC35j9s&t=51s

# Hyperparameters
alpha = 1
l = 2

# Parameters
n = 50 # nb of points
n_func = 10 # nb of fct to draw
x_bound = 5 # bound on the x axis

def RBF_kernel(a,b):
    sqdist = np.sum(a**2,1).reshape(-1,1) + np.sum(b**2,1) - 2*np.dot(a,b.T)
    return alpha**2 * np.exp(-1/l * sqdist)

n = 50
X = np.linspace(-x_bound, x_bound, n).reshape(-1,1)
K = RBF_kernel(X, X) # dim(K) = n x n

L = np.linalg.cholesky(K + 1e-6 * np.eye(n))
f_prior = np.dot(L, np.random.normal(size=(n, n_func)))

# Plotting
width = 16
height = 4
plt.figure(figsize=(width, height))

# plot f_prior
plt.subplot(1,3,1)
plt.title('GP: prior on f')
plt.plot(X, f_prior)
plt.plot(X, f_prior.mean(1), linewidth=3, color='black')
plt.ylabel('f(x)')
plt.xlabel('x')

# plot Kernel
plt.subplot(1,3,2)
plt.title('Kernel matrix')
plt.pcolor(K[::-1])
plt.colorbar()

plt.subplot(1,3,3)
plt.title('Kernel function')
plt.plot(X, RBF_kernel(X, np.array([[1.0]])))
plt.show()

Kernel (prior knowledge)

By choosing a specific kernel, we can incorporate prior knowledge that we have about the function f, such as, if the function is:

• periodic
• smooth
• symmetric
• etc.

The hyperparameters for each kernel are also very intuitive/interpretable.

Note: kernels can be combined!

Indeed, if k(x,y), k_1(x,y) and k_2(x,y) are valid kernels then:

• \alpha k(x,y) with \alpha \geq 0
• k_1(x,y) + k_2(x,y)
• k_1(x,y) k_2(x,y)
• p(k(x,y)) with p being a polynomial function with non-negative coefficients
• exp(k(x,y))
• f(x) k(x,y) \overline{f(y)} with \overline{f} = complex conjugate
• k(\phi(x),\phi(y))

are all valid kernels!

Periodic Exponential kernel

In [ ]:

variance = 1.
lengthscale = 1.
period = 2.*np.pi

#K = periodic_kernel(X, X) # dim(K) = n x n
kern = GPy.kern.PeriodicExponential(variance=variance, lengthscale=lengthscale, period=period)
K1 = kern.K(X)

L = np.linalg.cholesky(K1 + 1e-6 * np.eye(n))
f_prior = np.dot(L, np.random.normal(size=(n, 1)))

# Plotting
width = 16
height = 4
plt.figure(figsize=(width, height))

# plot f_prior
plt.subplot(1,3,1)
plt.title('GP: prior on f')
plt.plot(X, f_prior)
plt.plot(X, f_prior.mean(1), linewidth=3, color='black')
plt.ylabel('f(x)')
plt.xlabel('x')

# plot Kernel
plt.subplot(1,3,2)
plt.title('Kernel matrix')
plt.pcolor(K1[::-1])
plt.colorbar()

plt.subplot(1,3,3)
plt.title('Kernel function')
plt.plot(X, kern.K(X, np.array([[1.0]])))
plt.show()

addition and multiplication of 2 kernels (SE and PE)

In [ ]:

K_add = K + K1

L = np.linalg.cholesky(K_add + 1e-6 * np.eye(n))
f_prior = np.dot(L, np.random.normal(size=(n, n_func)))

# Plotting
width = 16
height = 8
plt.figure(figsize=(width, height))

# plot f_prior
plt.subplot(2,2,1)
plt.title('GP: prior on f with K_add')
plt.plot(X, f_prior)
plt.plot(X, f_prior.mean(1), linewidth=3, color='black')
plt.ylabel('f(x)')
plt.xlabel('x')

# plot Kernel
plt.subplot(2,2,2)
plt.title('Kernel matrix: K_add')
plt.pcolor(K_add[::-1])
plt.colorbar()

K_prod = K * K1

L = np.linalg.cholesky(K_prod + 1e-6 * np.eye(n))
f_prior = np.dot(L, np.random.normal(size=(n, n_func)))

# plot f_prior
plt.subplot(2,2,3)
plt.title('GP: prior on f with K_prod')
plt.plot(X, f_prior)
plt.plot(X, f_prior.mean(1), linewidth=3, color='black')
plt.ylabel('f(x)')
plt.xlabel('x')

# plot Kernel
plt.subplot(2,2,4)
plt.title('Kernel matrix: K_prod')
plt.pcolor(K_prod[::-1])
plt.colorbar()
plt.show()

GP Posterior

Given \mathcal{D}=\{(x_i, y_i)\}_{i=1}^{i=N} = (\pmb{X}, \pmb{y}), we have:

\begin{equation} p(f|\mathcal{D}) = \frac{p(\mathcal{D}|f)p(f)}{p(\mathcal{D})} \end{equation}

GP Regression

\begin{equation} y_i = f(\pmb{x}_i) + \epsilon_i \qquad \left{ \begin{array}{l} f \sim \mathcal{GP}(\pmb{0}, \pmb{K}) \ \epsilon_i \sim \mathcal{N}(0, \sigma^2) \end{array} \right. \end{equation}

• Prior f is a GP \Leftrightarrow p(\pmb{f}|\pmb{X}) = \mathcal{N}(\pmb{0}, \pmb{K})
• Likelihood is Gaussian \Leftrightarrow p(\pmb{y}|\pmb{X},\pmb{f}) = \mathcal{N}(\pmb{f}, \sigma^2\pmb{I})
• \rightarrow p(f|\mathcal{D}) is also a GP.

Predictive distribution:

p(\pmb{y}_*|\pmb{x}_*, \pmb{X}, \pmb{y}) = \int p(\pmb{y}_{*}| \pmb{x}_{*}, \pmb{f}, \pmb{X}, \pmb{y}) p(f|\pmb{X}, \pmb{y}) d\pmb{f} = \mathcal{N}(\pmb{\mu}_*, \pmb{\Sigma}_*)

\begin{align} \pmb{\mu}* &= \pmb{K}{N} (\pmb{K}N + \sigma^2 \pmb{I})^{-1} \pmb{y} \ \pmb{\Sigma} &= \pmb{K}{**} - \pmb{K}{N} (\pmb{K}N + \sigma^2 \pmb{I})^{-1} \pmb{K}{N} \end{align}

Learning a GP

Marginal likelihood:

\begin{equation} p(\pmb{y}|\pmb{X}) = \int p(\pmb{y}|\pmb{f},\pmb{X}) p(\pmb{f}|\pmb{X}) d\pmb{f} = \mathcal{N}(\pmb{0}, \pmb{K} + \sigma^2\pmb{I}) \end{equation}

By taking the logarithm, and setting \pmb{K}_y = (\pmb{K} + \sigma^2\pmb{I}), we have:

\begin{equation} \mathcal{L} = \log p(\pmb{y}|\pmb{X}; \pmb{\Phi}) = \underbrace{-\frac{1}{2} \pmb{y}^T \pmb{K}y^{-1} \pmb{y}}{\mbox{data fit}} \underbrace{-\frac{1}{2} \log |\pmb{K}y^{-1}|}{\mbox{complexity penalty}} - \frac{n}{2} \log 2\pi \end{equation}

The marginal likelihood (i.e. ML-II) is used to optimize the hyperparameters \pmb{\Phi} that defines the covariance function and thus the GP.

\begin{equation} \pmb{\Phi}^* = argmax_{\pmb{\Phi}} \log p(\pmb{y}|\pmb{X}; \pmb{\Phi}) \end{equation}

Optimizing the marginal likelihood is more robust than the likelihood as it tries to optimize the complexity of the model, and the fitting of this last one to the observed data.

GPy

In [ ]:

# GP Regression
# Based on the tutorial: https://github.com/SheffieldML/notebook/blob/master/GPy/GPyCrashCourse.ipynb

# Create dataset
X = np.random.uniform(-3.0, 3.0, (20,1))
Y = np.sin(X) + np.random.randn(20,1) * 0.05 

# Create the kernel
# Reminder 1: The sum of valid kernels gives a valid kernel.
# Reminder 2: The product of valid kernels gives a valid kernel.
# Available kernels: RBF, Exponential, Matern32, Matern52, Brownian, Bias, Linear, PeriodicExponential, White.
kernel = GPy.kern.RBF(input_dim=1, variance=1.0, lengthscale=1.0)

# Create the model
gp_model = GPy.models.GPRegression(X, Y, kernel)

# Display and plot
print("Before optimization: ", gp_model)
gp_model.plot()
plt.show()

# Optimize the model (that is find the 'best' hyperparameters of the kernel matrix)
# By default, the optimizer is a 2nd order algo: lbfgsb. Others are available such as the scg, ...
gp_model.optimize(messages=False)

# Display and plot
print("After optimization: ", gp_model)
gp_model.plot()
plt.show()

Gaussian Process Latent Variable Model (GP-LVM)

In [ ]:

# GPLVM
# Based on the tutorials: 
# http://nbviewer.jupyter.org/github/SheffieldML/notebook/blob/master/GPy/MagnificationFactor.ipynb
# https://github.com/SheffieldML/notebook/blob/master/lab_classes/gprs/lab4-Copy0.ipynb

# Create dataset
N = 100
k1 = GPy.kern.RBF(5, variance=1, lengthscale=1./np.random.dirichlet(np.r_[10,10,10,0.1,0.1]), ARD=True)
k2 = GPy.kern.RBF(5, variance=1, lengthscale=1./np.random.dirichlet(np.r_[0.1,10,10,10,0.1]), ARD=True)
X = np.random.normal(0, 1, (N,5))
A = np.random.multivariate_normal(np.zeros(N), k1.K(X), 10).T
B = np.random.multivariate_normal(np.zeros(N), k2.K(X), 10).T

Y = np.vstack((A,B))

# latent space dimension
latent_dim = 2

# Create the kernel
kernel = GPy.kern.RBF(input_dim=latent_dim, variance=1.0, lengthscale=1.0)

# Create the GPLVM model
gplvm_model = GPy.models.GPLVM(Y, latent_dim, init='PCA', kernel=kernel)

# Display and plot
print("Before optimization: ", gplvm_model)
gplvm_model.plot_latent()
plt.show()

# Optimize the model (that is find the 'best' hyperparameters of the kernel matrix)
# By default, the optimizer is a 2nd order algo: lbfgsb. Others are available such as the scg, ...
gplvm_model.optimize(messages=False)

# Display and plot
print("After optimization: ", gplvm_model)
gplvm_model.plot_latent()
plt.show()