Gaussian process latent variable models for visualisation of high dimensional data

Gaussian Process Latent Variable Models for Visualisation of High Dimensional Data
Neil D. Lawrence Department of Computer Science, University of Shef?eld, Regent Court, 211 Portobello Street, Shef?eld, S1 4DP, U.K. neil@https://www.wendangku.net/doc/516542550.html,
Abstract
In this paper we introduce a new underlying probabilistic model for principal component analysis (PCA). Our formulation interprets PCA as a particular Gaussian process prior on a mapping from a latent space to the observed data-space. We show that if the prior’s covariance function constrains the mappings to be linear the model is equivalent to PCA, we then extend the model by considering less restrictive covariance functions which allow non-linear mappings. This more general Gaussian process latent variable model (GPLVM) is then evaluated as an approach to the visualisation of high dimensional data for three different data-sets. Additionally our non-linear algorithm can be further kernelised leading to ‘twin kernel PCA’ in which a mapping between feature spaces occurs.
1 Introduction
Visualisation of high dimensional data can be achieved through projecting a data-set onto a lower dimensional manifold. Linear projections have traditionally been preferred due to the ease with which they can be computed. One approach to visualising a data-set in two dimensions is to project the data along two of its principal components. If we were forced to choose a priori which components to project along, we might sensibly choose those associated with the largest eigenvalues. The probabilistic reformulation of principal component analysis (PCA) also informs us that choosing the ?rst two components is also the choice that maximises the likelihood of the data [11]. 1.1 Integrating Latent Variables, Optimising Parameters Probabilistic PCA (PPCA) is formulated as a latent variable model: given a set centred of ? -dimensional data Yò ò ? and denoting the latent variable associated with each datapoint üò we may write the likelihood for an individual data-point under the PPCA model as
? ′Yò ? ? μ
? ′Yò üò
?
? μ ? ′üò μ
üò

where ? ′üò μ is Gaussian distributed with unit covariance, ? ′üò μ ? ′üò ? áμ, and ? ? ′Yò üò ? ? μ ? Yò ?üò ? ? á . The solution for ? can then be found1 by assuming that Yò is i.i.d. and maximising the likelihood of the data-set, ? ? ′ ? ?μ ? ′Yò ? ? μ ò ? T where Y? Y? ? is the ? ￠ design matrix. Probabilistic principal component analysis and other latent variable models, such as factor analysis (FA) or independent component analysis (ICA), require a marginalisation of the latent variables and optimisation of the parameters. In this paper we consider the dual approach of marginalising ? and optimising each üò . This probabilistic model also turns out to be equivalent to PCA. 1.2 Integrating Parameters, Optimising Latent Variables
?? By ?rst specifying a prior distribution, ? ′?μ ? ? ? ? ? á where ? is the th row of the matrix ?, and integrating over ? we obtain a marginalised likelihood for ,
é
?
?
?′
where ? then
?μ
T
? ′? μ ?
?
?
·
?
? ? á and ?
? ü?
?
ü?
T
? ? tr
?
???
T
?
(1)
ü? ? . The corresponding log-likelihood is
Dò ? ? ? tr
?
?
Dò ′? μ
?
?
?
???
T
?
(2)
Now that the parameters are marginalised we may focus on optimisation of the likelihood with respect to the . The gradients of (2) with respect to may be found as,
which implies that at our solution
?
T
? 
????
T
???
??
???
??? í? ??
T
some algebraic manipulation of this formula [11] leads to where í? is an ? ￠ are eigenvectors of
? ?
?
T ? , and ? is an arbitrary ? ￠ ? , where is the th eigenvalue of ? ?? 2 orthogonal matrix . Note that the eigenvalue problem we have developed can easily be shown to be equivalent to that solved in PCA (see e.g. [10]), indeed the formulation of T is PCA in this manner is a key step in the development of kernel PCA [9] where replaced with a kernel. Our probabilistic PCA model shares an underlying structure with [11] but differs in that where they optimise we marginalise and where they marginalise we optimise. The marginalised likelihood we are optimising in (1) is recognised as the product of independent Gaussian processes where the (linear) covariance function is given by ? T · ? ? á. Therefore a natural extension is the non-linearisation of the mapping from latent space to the data space through the introduction of a non-linear covariance function.
1
matrix (? is the dimension of the latent space) whose columns T , ? is a ? ￠ ? diagonal matrix whose th element is D
it.
As can the solution for ? but since the solution for ? is not dependent on ? we will disregard
2 For independent component analysis the correct rotation matrix ? must also be found, here we have placed no constraints on the orientation of the axes so this matrix cannot be recovered.

2 Gaussian Process Latent Variable Models
We saw in the previous section how PCA can be interpreted as a Gaussian process ‘mapping3 ’ from a latent space to a data space where the locale of the points in latent space is determined by maximising the Gaussian process likelihood with respect to . We will refer to models of this class as Gaussian process latent variable models (GPLVM). Principal component analysis is a GPLVM where the process prior is based on the ? ￠ ? inner product matrix of , in this section we develop an alternative GPLVM by considering a prior which allows for non-linear processes, speci?cally we focus on the popular ‘RBF kernel’ which takes the form
ò?
?
ü?
- ′üò ü? μT ′üò ü? μ ?
·
?ò? ? ?
where ò ? is the element in the òth row and ?th column of ?, - is a scale parameter and ?ò? denotes the Kronecker delta. Gradients of (2) with respect to the latent points can be found through combining
 with ü? via the chain rule. These gradients may be used in combination with (2) in a ò non-linear optimiser such as scaled conjugate gradients (SCG) [7] to obtain a latent variable representation of the data. Furthermore gradients with respect to the parameters of the kernel matrix may be computed and used to jointly optimise , ?, - and ? . The solution for will naturally not be unique; even for the linear case described above the solution is subject to an arbitrary rotation, here we may expect multiple local minima.
2.1 Illustration of GPLVM via SCG To illustrate a simple Gaussian process latent variable model we turn to the ‘multi-phase oil ?ow’ data [2]. This is a twelve dimensional data-set containing data of three known classes corresponding to the phase of ?ow in an oil pipeline: strati?ed, annular and homogeneous. In this illustration, for computational reasons, the data is sub-sampled to 100 data-points. Figure 1 shows visualisations of the data using both PCA and our GPLVM algorithm which required 766 iterations of SCG. The positions for the GPLVM model were initialised using PCA (see https://www.wendangku.net/doc/516542550.html,/~neil/gplvm/ for the MATLAB code used). The gradient based optimisation of the RBF based GPLVM’s latent space shows results which are clearly superior (in terms of greater separation between the different ?ow domains) to those achieved by the linear PCA model. Additionally the use of a Gaussian process to perform our ‘mapping’ means that there is uncertainty in the positions of the points in the data space. For our formulation the level of uncertainty is shared across all4 dimensions and thus may be visualised in the latent space. In Figure 1 (and subsequently) this is done through varying the intensity of the background pixels. Unfortunately, a quick analysis of the complexity of the algorithm shows that each gradient ? step requires an inverse of the kernel matrix, an ? ? ? operation, rendering the algorithm impractical for many data-sets of interest.
3 Strictly speaking the model does not represent a mapping as a Gaussian process ‘maps’ to a distribution in data space rather than a point. 4 This apparent weakness in the model may be easily recti?ed to allow different levels of uncertainty for each output dimension, our more constrained model allows us to visualise this uncertainty in the latent space and is therefore preferred for this work.
? ?
???
T
???
?
???

0.3
2.4 0.2
0.2
0.15
2.2 2 0.1 1.8
0.1
0
0.05 1.6 0
?0.1
?0.05
1.4 1.2 ?0.1 1
?0.2
?0.3
?0.15 0.8 ?0.2
?0.4
0.6 ?0.25 0.4
?0.5 ?0.8
?0.6
?0.4
?0.2
0
0.2
0.4
0.6
?0.2
?0.1
0
0.1
0.2
0.3
Figure 1: Visualisation of the Oil data with (a) PCA (a linear GPLVM) and (b) A GPLVM which uses an RBF kernel. Crosses, circles and plus signs represent strati?ed, annular and homogeneous ?ows respectively. The greyscales in plot (b) indicate the precision with which the manifold was , expressed in data-space for that latent point. The optimised parameters of the kernel were - ? and ? .
? ??
??
??
2.2 A Practical Algorithm for GPLVMs There are three main components to our revised, computationally ef?cient, optimisation process: Sparsi?cation. Kernel methods may be sped up through sparsi?cation, i.e. representing the data-set by a subset, á , of points known as the active set. The remainder, the inactive set, is denoted by ? . We make use of the informative vector machine [6] which selects points sequentially according to the reduction in the posterior process’s entropy that they induce. Latent Variable Optimisation. A point from the inactive set, , can be shown to project into the data space as a Gaussian distribution
? ? ′Y ü ? ? - μ ? Y á (3) T ? ? á where ?á á denotes the kernel matrix developed from
?
?
whose mean is áá the active set and á is a column vector consisting of the elements from the th column of ? ? that correspond to the active set. The variance is ? ′ü ü μ T ?á á á Note á that since ü does not appear in the inverse, gradients with respect to ü do not depend on other data in ? . We can therefore independently optimise the likelihood of each Y with respect to each ü . Thus the full set ? can be optimised with one pass through the data. Kernel Optimisation. The likelihood of the active set is given by
? ′ áμ ?
?
? ′? μ ?
?á á
? ?
ü?
?? ?
-
á ?á á á
T
??
(4)
which?can be optimised5 with respect to
?, ?
and
with gradient evaluations costing
Algorithm 1 summarises the order in which we implemented these steps. Note that whilst we never optimise points in the active set, we repeatedly reselect the active set so it is
5 In practice we looked for MAP solutions for all our optimisations, specifying a unit covariance Gaussian prior for the matrix and using ?, ? and - for ?, ? and - respectively.
?
?
?

Algorithm 1 An algorithm for modelling with a GPLVM. Require: A size for the active set, . A number of iterations, ì . Initialise through PCA. for ì iterations do Select a new active set using the IVM algorithm. Optimise (4) with respect to the parameters of ? using scaled conjugate gradients. Select a new active set. for Each point not in active set, . do Optimise (3) with respect to ü using scaled conjugate gradients. end for end for
unlikely that many points remain in their original location. For all the experiments that follow we used ì ? iterations and an active set of size ???. The experiments were run on a ‘one-shot’ basis6 so we cannot make statements as to the effects that signi?cant modi?cation of these parameters would have. We present results on three data-sets: for the oil ?ow data (Figure 2) from the previous section we now make use of all 1000 available points and we include a comparison with the generative topographic mapping (GTM) [4].
1
2
0.8
1.5
1.8
0.6
1 1.6
0.4
0.5
0.2
0
1.4
0
?0.5 1.2
?0.2
?1 1 ?1.5 0.8
?0.4
?0.6
?0.8
?2
?1 ?1
?0.8
?0.6
?0.4
?0.2
0
0.2
0.4
0.6
0.8
1
?2.5 ?3
?2
?1
0
1
2
3
Figure 2: The full oil ?ow data-set visualised with (a) GTM with 225 latent points laid out on a
grid and with 16 RBF nodes and (b) an RBF based GPLVM. The parameters of the latent ,? and - . Notice how the GTM arti?cially variable model were found to be ? ‘discretises’ the latent space around the locations of the 225 latent points.
? ￠?
? ??
??
We follow [5] in our 2-D visualisation of a sub-set of 3000 of the digits 0-4 (600 of each digit) from a ? ￠ ? greyscale version of the USPS digit data-set (Figure 3). Finally we modelled a face data-set [8] consisting of 1965 images from a video sequence digitised at ? ￠ ??. Since the images are originally from a video sequence we might expect the underlying dimensionality of the data to be one — the images are produced in a smooth way over time which can be thought of as a piece of string embedded in a high (560) dimensional pixel space. We therefore present ordered results from a 1-D visualisation in Figure 4 . All the code used for performing the experiments is available from http://www.dcs.
By one-shot we mean that, given the algorithm above, each experiment was only run once with one setting of the random seed and the values of ì and given. If we were producing a visualisation for only one dataset this would leave us open to the criticism that our one-shot result was ‘lucky’. However we present three data-sets in what follows and using a one-shot approach in problems with multiple local minima removes the temptation of preferentially selecting ‘prettier’ results.
6

0.8
4
0.75 0.7
2
0.65 0.6
0
0.55
?2
0.5
0.45
?4 ?4 ?2 0 2 4
0.4
Figure 3: The digit images visualised in the 2-D latent space. We followed [5] in plotting images in
a random order but not plotting any image which would overlap an existing image. 538 of the 3000 digits are plotted. Note how little space is taken by the ‘ones’ (the thin line running from (-4, -1.5) to (-1, 0)) in our visualisation, this may be contrasted with the visualisation of a similar data-set in [5]. We suggest this is because ‘ones’ are easier to model and therefore do not require a large region in latent space.
https://www.wendangku.net/doc/516542550.html,/~neil/gplvm/ along with avi video ?les of the 1-D visualisation and results from two further experiments on the same data (a 1-D GPLVM model of the digits and a 2-D GPLVM model of the faces).
3 Discussion
Empirically the RBF based GPLVM model gives useful visualisations of a range of datasets. Strengths of the method include the ability to optimise the kernel parameters and to generate fantasy data from any point in latent space. Through the use of a probabilistic process we can obtain error bars on the position of the manifolds which can be visualised by imposing a greyscale image upon the latent space. When Kernels Collide: Twin Kernel PCA The eigenvalue problem which provides the maxima of (2) with respect to for the linear kernel is exploited in kernel PCA. One could T T · ? ? á and consider a ‘twin kernel’ PCA where both ? are replaced by kernel functions. Twin kernel PCA could no longer be undertaken with an eigenvalue decomposition but Algorithm 1 would still be a suitable mechanism with which to determine the values of and the parameters of ’s kernel.

Figure 4: Top: Fantasy faces from the 1-D model for the face data. These faces were created by taking 64 uniformly spaced and ordered points from the latent space and visualising the mean of their distribution in data space. The plots above show this sequence unfolding (starting at the top left and moving right). Ideally the transition between the images should be smooth. Bottom: Examples from the data-set which are closest to the corresponding fantasy images in latent space. Full sequences of 2000 fantasies and the entire dataset are available on the web as avi ?les. Stochastic neighbor embedding. Consider that (2) could be written as ? ê ?Y ′T ? ?Y μ Dò ?ü ′T ? ?ü μ T where we have introduced a vector, T, of length ? , T ? ?Y and we have rede?ned ? as ?ü . The entropy of ?Y ′T ? ?Y μ is constant7 in , we therefore may add it to ? to obtain KL ′?Y
?ü μ
?Y ′T ?
?Y μ Dò
?ü ′T ? ?Y ′T ?
?ü μ ?Y μ
T
(5)
which is recognised Kullback-Leibler (KL) divergence between the two distributions. Stochastic neighbor embedding (SNE) [5] also minimises this KL divergence to visualise data. However in SNE the vector T is discrete. Generative topographic mapping. The Generative topographic mapping [3] makes use of a radial basis function network to perform the mapping from latent space to observed space. Marginalisation of the latent space is achieved with an expectation-maximisation
7 Computing the entropy requires ?Y to be of full rank, this is not true in general but can be forced by adding ‘jitter’ to ?Y , e.g. ?Y ?Y - ? á.
·

(EM) algorithm. A radial basis function network is a special case of a generalised linear model and can be interpreted as a Gaussian process. Under this interpretation the GTM becomes GPLVM with a particular covariance function. The special feature of the GTM is the manner in which the latent space is represented, as a set of uniformly spaced delta functions. One could view the GPLVM as having a delta function associated with each data-point: in the GPLVM the positions of the delta functions are optimised, in the GTM each data point is associated with several different ?xed delta functions.
4 Conclusions
We have presented a new class of models for probabilistic modelling and visualisation of high dimensional data. We provided strong theoretical grounding for the approach by proving that principal component analysis is a special case. On three real world data-sets we showed that visualisations provided by the model cluster the data in a reasonable way. Our model has an advantage over the various spectral clustering algorithms that have been presented in recent years in that, in common with the GTM, it is truly generative with an underlying probabilistic interpretation. However it does not suffer from the arti?cial ‘discretetisation’ suffered by the GTM. Our theoretical analysis also suggested a novel non-linearisation of PCA involving two kernel functions. Acknowledgements We thank Aaron Hertzmann for comments on the manuscript.
References
[1] S. Becker, S. Thrun, and K. Obermayer, editors. Advances in Neural Information Processing Systems, volume 15, Cambridge, MA, 2003. MIT Press. [2] C. M. Bishop and G. D. James. Analysis of multiphase ?ows using dual-energy gamma densitometry and neural networks. Nuclear Instruments and Methods in Physics Research, A327:580– 593, 1993. [3] C. M. Bishop, M. Svensén, and C. K. I. Williams. GTM: a principled alternative to the SelfOrganizing Map. In Advances in Neural Information Processing Systems, volume 9, pages 354–360. MIT Press, 1997. [4] C. M. Bishop, M. Svensén, and C. K. I. Williams. GTM: the Generative Topographic Mapping. Neural Computation, 10(1):215–234, 1998. [5] G. Hinton and S. Roweis. Stochastic neighbor embedding. In Becker et al. [1], pages 857–864. [6] N. D. Lawrence, M. Seeger, and R. Herbrich. Fast sparse Gaussian process methods: The informative vector machine. In Becker et al. [1], pages 625–632. [7] I. T. Nabney. Netlab: Algorithms for Pattern in Pattern Recognition. Springer, Berlin, 2001. https://www.wendangku.net/doc/516542550.html,/netlab/. Recognition. Advances Code available from
[8] S. Roweis, L. K. Saul, and G. Hinton. Global coordination of local linear models. In T. G. Dietterich, S. Becker, and Z. Ghahramani, editors, Advances in Neural Information Processing Systems, volume 14, pages 889–896, Cambridge, MA, 2002. MIT Press. [9] B. Sch?lkopf, A. J. Smola, and K.-R. Müller. Kernel principal component analysis. In Proceedings 1997 International Conference on Arti?cial Neural Networks, ICANN’97, page 583, Lausanne, Switzerland, 1997. [10] M. E. Tipping. Sparse kernel principal component analysis. In T. K. Leen, T. G. Dietterich, and V. Tresp, editors, Advances in Neural Information Processing Systems, volume 13, pages 633–639, Cambridge, MA, 2001. MIT Press. [11] M. E. Tipping and C. M. Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Society, B, 6(3):611–622, 1999.