Gaussian Process Road School, 16th February 2014 written by Max Zwiessele, Neil Lawrence
This lab session will focus on three aspects of GPs: sampling, the design of Experiments and uncertainty propagation.
As for the last two lab sessions, the first thing to do is set the plots to appear inline and import relevant modules for the lab session.
In [1]:
%pylab inline
import numpy as np
import pylab as pb
import GPy
import string
For this lab, we've created a dataset digits.npy containing all handwritten digits from $0 \cdots 9$ handwritten, provided by deCampos et al. [2009]. All digits were cropped and scaled down to an appropriate format.
You can retrieve the dataset as follows:
In [2]:
import urllib
urllib.urlretrieve('http://staffwww.dcs.sheffield.ac.uk/people/J.Hensman/gpsummer/Lab3.zip', 'Lab3.zip')
import zipfile
zip = zipfile.ZipFile('Lab3.zip', 'r')
for name in zip.namelist():
zip.extract(name, '.')
from load_plotting import * # for plotting
We will only use some of the digits for the demonstrations in this lab class, but you can edit the code below to select different subsets of the digit data as you wish.
In [3]:
digits = np.load('digits.npy')
which = [0,1,2,6,7,9] # which digits to work on
digits = digits[which,:,:,:]
num_classes, num_samples, height, width = digits.shape
labels = np.array([[str(l)]*num_samples for l in which])
You can try to plot some sample using pb.matshow
Principal component analysis (PCA) finds a rotation of the observed outputs, such that the rotated principal component (PC) space maximizes the variance of the data observed, sorted from most to least important (most to least variable in the corresponding PC).
In order to apply PCA in an easy way, we have included a PCA module in pca.py. You can import the module by import <path.to.pca> (without the ending .py!). To run PCA on the digits we have to reshape (Hint: np.reshape ) digits .
We will call the reshaped observed outputs $\mathbf{Y}$ in the following.
In [4]:
Y = digits.reshape((digits.shape[0]*digits.shape[1],digits.shape[2]*digits.shape[3]))
Yn = Y-Y.mean()
Now let’s run PCA on the reshaped dataset $\mathbf{Y}$:
In [5]:
import pca
p = pca.PCA(Y) # create PCA class with digits dataset
The resulting plot will show the lower dimensional representation of the digits in 2 dimensions.
In [6]:
p.plot_fracs(20) # plot first 20 eigenvalue fractions
p.plot_2d(Y,labels=labels.flatten(), colors=colors)
pb.legend()
Out[6]:
The Gaussian Process Latent Variable Model (GP-LVM) embeds of PCA into a Gaussian process framework, where the latent inputs $\mathbf{X}$ are learnt as hyperparameters and the mapping variables $\mathbf{W}$ are integrated out. The advantage of this interpretation is it allows PCA to be generalized in a non linear way by replacing the resulting linear covariance witha non linear covariance. But first, let's see how GPLVM is equivalent to PCA using an automatic relevance determination (ARD, see e.g. Bishop et al. [2006]) linear kernel:
In [36]:
input_dim = 4 # How many latent dimensions to use
kernel = GPy.kern.Linear(input_dim, ARD=True) # ARD kernel
#kernel += GPy.kern.White(input_dim) + GPy.kern.Bias(input_dim)
m = GPy.models.GPLVM(Yn, input_dim=input_dim, kernel=kernel)
#m.ensure_default_constraints()
m['.*noise'] = m.Y.var()/20. # start noise is 5% of datanoise
m.optimize(messages=1, max_iters=1000) # optimize for 1000 iterations
print 'done'
m.kern.plot_ARD()
plot_model(m, m['.*linear.variances'].argsort()[-2:], labels.flatten())
pb.legend()
Out[36]:
In [38]:
print m
As you can see the solution with a linear kernel is the same as the PCA solution with the exception of rotational changes and axis flips.
For the sake of time, the solution you see was only running for 1000 iterations, thus it might not be converged fully yet. The GP-LVM proceeds by iterative optimization of the inputs to the covariance. As we saw in the lecture earlier, for the linear covariance, these latent points can be optimized with an eigenvalue problem, but generally, for non-linear covariance functions, we are obliged to use gradient based optimization.
a) How do your linear solutions differ between PCA and GPLVM with a linear kernel? Look at the plots and also try and consider how the linear ARD parameters compare to the eigenvalues of the principal components.
b) The next step is to use a non-linear mapping between inputs $\mathbf{X}$ and ouputs $\mathbf{Y}$ by selecting the exponentiated quadratic (GPy.kern.rbf) covariance function.
In [39]:
kern = GPy.kern.RBF(input_dim)
c) How does the nonlinear model differe from the linear model? Are there digits that the GPLVM with an exponentiated quadratic covariance can separate, which PCA is not able to?
d) Try modifying the covariance function and running the model again. For example you could try a combination of the linear and exponentiated quadratic covariance function or the Matern 5/2. If you run into stability problems try initializing the covariance function parameters differently.
In GP-LVM we use a point estimate of the distribution of the input $\mathbf{X}$. This estimate is derived through maximum likelihood or through a maximum a posteriori (MAP) approach. Ideally, we would like to also estimate a distribution over the input $\mathbf{X}$. In the Bayesian GPLVM we approximate the true distribution $p(\mathbf{X}|\mathbf{Y})$ by a variational approximation $q(\mathbf{X})$ and integrate $\mathbf{X}$ out.
Approximating the posterior in this way allows us to optimize a lower bound on the marginal likelihood. Handling the uncertainty in a principled way allows the model to make an assessment of whether a particular latent dimension is required, or the variation is better explained by noise. This allows the algorithm to switch off latent dimensions. The switching off can take some time though, so below in Section 6 we provide a pre-learnt module, but to complete section 6 you'll need to be working in the IPython console instead of the notebook.
For the moment we'll run a short experiment applying the Bayesian GP-LVM with an exponentiated quadratic covariance function.
In [41]:
# Model optimization
input_dim = 5 # How many latent dimensions to use
kern = GPy.kern.RBF(input_dim,ARD=True) # ARD kernel
m = GPy.models.BayesianGPLVM(Yn, input_dim=input_dim, kernel=kern, num_inducing=25)
# initialize noise as 1% of variance in data
m['*.noise.variance'] = m.likelihood.Y.var()/100.
m.optimize('scg', messages=1, max_iters=1000)
In [ ]:
# Plotting the model
plot_model(m, m['rbf_len'].argsort()[:2], labels.flatten())
pb.legend()
m.kern.plot_ARD()
# Saving the model:
m.pickle('bgplvm_rbf.pickle')
Because we are now also considering the uncertainty in the model, this optimization can take some time. However, you are free to interrupt the optimization at any point selecting Kernel->Interupt from the notepad menu. This will leave you with the model, m in the current state and you can plot and look into the model parameters.
a) How does the Bayesian GP-LVM compare with the standard model?
A good way of working with latent variable models is to interact with the latent dimensions, generating data. This is a little bit tricky in the notebook, so below in section 6 we provide code for setting up an interactive demo in the standard IPython shell. If you are working on your own machine you can try this now. Otherwise continue with section 5.
In Manifold Relevance Determination we try to find one latent space, common for $K$ observed output sets (modalities) $\{\mathbf{Y}_{k}\}_{k=1}^{K}$. Each modality is associated with a separate set of ARD parameters so that it switches off different parts of the whole latent space and, therefore, $\mathbf{X}$ is softly segmented into parts that are private to some, or shared for all modalities. Can you explain what happens in the following example?
Again, you can stop the optimizer at any point and explore the result obtained with the so far training:
In [ ]:
m = GPy.examples.dimensionality_reduction.mrd_simulation(optimize = False)
m.optimize(messages = True, max_iters=1000, optimizer = 'bfgs')
In [ ]:
m.plot_scales()
m.plot_X_1d()
The simulated data set is a sinusoid and a double frequency sinusoid function as input signals.
a) Which signal is shared across the three datasets?
b) Which are private?
c) Are there signals shared only between two of the three datasets?
The module below loads a pre-optimized Bayesian GPLVM model (like the one you just trained) and allows you to interact with the latent space. Three interactive figures pop up: the latent space, the ARD scales and a sample in the output space (corresponding to the current selected latent point of the other figure). You can sample with the mouse from the latent space and obtain samples in the output space. You can select different latent dimensions to vary by clicking on the corresponding scales with the left and right mouse buttons. This will also cause the latent space to be projected on the selected latent dimensions in the other figure.
In [42]:
run load_bgplvm_dimension_select.py
In [43]:
import cPickle as pickle
with open('./digit_bgplvm_demo.pickle', 'rb') as f:
m = pickle.load(f)
Prepare for plotting of this model. If you run on a webserver the interactive plotting will not work. Thus, you can skip to the next codeblock and run it on your own machine, later.
In [ ]:
fig = pb.figure('Latent Space & Scales', figsize=(16,6))
ax_latent = fig.add_subplot(121)
ax_scales = fig.add_subplot(122)
fig_out = pb.figure('Output', figsize=(1,1))
ax_image = fig_out.add_subplot(111)
fig_out.tight_layout(pad=0)
data_show = GPy.util.visualize.image_show(m.likelihood.Y[0:1, :], dimensions=(16, 16), transpose=0, invert=0, scale=False, axes=ax_image)
lvm_visualizer = GPy.util.visualize.lvm_dimselect(m.X.copy(), m, data_show, ax_latent, ax_scales, labels=labels.flatten())
Observations
Confirm the following observations by interacting with the demo:
You can also run the dimensionality reduction example
In [ ]:
GPy.examples.dimensionality_reduction.bgplvm_simulation()
Questions
C. M. Bishop. Pattern recognition and machine learning, volume 1. springer New York, 2006.
T. de Campos, B. R. Babu, and M. Varma. Character recognition in natural images. VISAPP 2009.
N. D. Lawrence. Probabilistic non-linear principal component analysis with Gaussian process latent variable models. In Journal of Machine Learning Research 6, pp 1783--1816, 2005