GPclust
Collapsed Variational Bayes
view repo
In this publication, we combine two Bayesian nonparametric models: the Gaussian Process (GP) and the Dirichlet Process (DP). Our innovation in the GP model is to introduce a variation on the GP prior which enables us to model structured timeseries data, i.e. data containing groups where we wish to model inter and intragroup variability. Our innovation in the DP model is an implementation of a new fast collapsed variational inference procedure which enables us to optimize our variationala pproximation significantly faster than standard VB approaches. In a biological time series application we show how our model better captures salient features of the data, leading to better consistency with existing biological classifications, while the associated inference algorithm provides a twofold speedup over EMbased variational inference.
READ FULL TEXT VIEW PDFCollapsed Variational Bayes
Clustering time series using Gaussian processes and Variational Bayes.
Collapsed Variational Bayes
We consider the problem of modelling and clustering structured time series. We turn to two tools from the field of Bayesian nonparametrics, using Gaussian Processes (GPs) to model time series and Dirichlet Processes to model clusterings.
Our model is constructed as follows. Given some data which is partitioned into disjoint groups (or batches), we construct a hierarchical GP model, using a GP to model each group, and a single additional GP to model the prior mean for the whole set. The general behaviour of the groups is governed by the last GP, and the deviations of each group from this mean behaviour is also modelled as a GP. We envisage that for many applications, these groupings may have subdivisions, which we model using additional layers of hierarchy. Further, we construct a model where the top level partition is unknown apriori, i.e. in the case of clustering, using a Dirichlet Process (DP) prior with a GP base distribution, each atom of which becomes the prior mean for a hierarchical GP. This allows us to perform inference over clusters of hierarchically grouped data.
Our model is inspired by the analysis of geneexpression timeseries data. Previous models for clustering timeseries using GP models (Dunson, 2010; Cooke et al., 2011) failed to account for structure in the data and previously proposed inference procedures (Gibbs sampling, agglomerative clustering) do not scale well. In this biological application, values of gene expression are measured at regular or irregular time intervals spanning some phenomenon such as development or disease progression. The high cost of measurement or limited temporal resolution of the underlying system usually dictates a small number of time points, and the measurement process is subject to both technical and biological variation. Groups in the data occur naturally: time series may be taken for different patients in a clinical trial, thus grouping them by patient number; measurements can be taken during development of related species or subject to replicated experiments. We can envisage more than one level of grouping: we might have different patients with measurements taken at different hospitals, say, or developmental time series taken by different laboratories using different technologies.
Whilst inference in a GP is tractable^{1}^{1}1subject to known or optimised hyperparameters and a Gaussian likelihood, inference in a DP requires some numerical procedure such as Gibbs sampling or variational approaches. Whilst variational methods are widely acknowledged to be faster than sampling, there is a need for even faster inference in these methods. In particular, faster inference allows for the exploration of larger datasets given the same computational resources and avoids the common practice of applying a crude filtering to reduce the size of data set prior to modelling. We propose a fast inference scheme based on recent work by Hensman et al. (2012b).
This novel derivation of variational Bayes (VB) amalgamates several key ideas from the literature. First, we construct a KLcorrected (King and Lawrence, 2006) bound on the marginal likelihood, where our objective function depends only on the approximate distribution of the clustering (latent) variables. The other model parameters are marginalised after constructing a lower bound on the conditional likelihood, and are not explicitly parameterised in optimization. This makes the parameterspace of the optimization significantly smaller. On this reduced parameter space, we make use of the Riemann structure of the approximation to derive the natural gradient, which is closely linked to the VBEM procedure. Using approximate geometrical conjugacy on the manifold (Honkela et al., 2010b), we implement a conjugate natural gradient approach that outperforms VBEM and freeform optimization.
Gaussian Process methods have been applied to gene expression timeseries with several aims, such as to infer transcription regulation (Honkela et al., 2010a), and to find dynamic differential expression (Kalaitzis and Lawrence, 2011; Stegle et al., 2010). Recently, we have proposed hierarchical Gaussian processes (Hensman et al., 2012a) for modelling gene expression, and showed that this simple proposal led to much improved modelling of the timeseries with various applications. Park and Choi (2010) also proposed a method based on hierarchical Gaussian processes. Their presentation is conceptually similar, but with the objective of saving some computation in performing Gaussian process regression, and Behseta et al. (2005)
proposed a hierarchical Gaussian process model with application to neuronal Poissonprocess intensities.
Clustering gene expression time series is an application which has attracted a lot of interest. Analysis of time series clusters is an important tool in exploring and understanding gene networks, whilst incorporating knowledge of the timeseries into the model has the potential to improve the ability of the method to discern clusters. Dunson (2010) proposed a DPGP model much like ours, but we make some important additions. In Dunson’s model, a series of GP functions are drawn from a DPGP prior, and each observation is then assigned to one of the functions. However Dunson makes no use of structure
in the model: observations differ from the latent function draws only by white noise. In our model, we use further GPs to describe how genes differ from the latent function, and
further GPs to describe more structure in the experiment, such as replications. Rasmussen and Ghahramani (2002) Also presented a method which combined GPs using DPs, but this method used a gating approach to produce a mixture of experts model, with different aims to that presented here.Outside the nonparametric framework, Medvedovic et al. (2004)
has described a hierarchical clustering model, but our approach is novel in applying structure to a GP model of the time series within the clustering.
Cooke et al. (2011) proposed a GP based clustering approach for gene expression, where replicates were used beforeclustering to estimate the level of noise in the experiment. Our richer model explicitly accounts for replicate structure in the experiment.
Our method for inferring structure and clustering in the model also offers an improvement over the aforementioned approaches. We show that our inference method is considerably faster than the usual variational method (VBEM), which is widely acknowledged to be faster than the Gibbs sampling approach adopted by Dunson (2010). The agglomerative clustering method proposed by Cooke et al. (2011) will suffer significant scalability issues: in the first round of agglomeration one needs to compute marginal likelihoods for GP models of every pair of genes. We shall derive a collapsed variational procedure for inference of clustering.
Collapsed variational Bayes (CVB) (Kurihara et al., 2007; Teh et al., 2007), and the latent variable method of Sung et al. (2008) share many properties of our approach, see (Hensman et al., 2012b) for details. Yet our method improves on CVB by considering the Riemannian structure of the collapsed problem, and relating it to the VBEM method. We can derive extremely efficient gradient methods, and apply conjugate gradient algorithms which account for the Riemman structure of the collapsed problem. Previously, Honkela et al. (2010b) proposed natural gradient based variational methods, but without connection to the collapsed approach they were unable to achieve a speedup over standard VBEM.
A related model which uses the KLcorrected approach is Overlapping mixtures of Gaussian Processes (LázaroGredilla et al., 2011). In this model, a series of latent GP functions are assumed, to which each observation is then assigned, with the objective of tracking. Our model can be reduced to this by removing the structured GP element, as well as the DP prior. Further (LázaroGredilla et al., 2011) used freeform variational optimization, which we shall show to be much slower than our approach.
In this section, we briefly review Gaussian Process regression and introduce our notation. We then extend the GPs to model structured time series before introducing our notation for mixture models.
Gaussian process (GP) regression is perhaps the most widely applied of Bayesian nonparametric methods, particularly since the publication of Rasmussen and Williams (2006)
. The idea is to place a prior over the space of functions, and use Bayesian inference to update one’s belief in the function by observing noisy realisations of the function at a finite number of points.
The GP is specified by a mean function and a covariance function . The mean function is often assumed to be zero everywhere and the covariance function takes a parametric form. In this publication we make use of the squareexponential covariance function:
. The parameters of the covariance function control the type of function permitted: the variance of the function is controlled by the parameter
, and the lengthscale of the function by. In this publication where we may have several separate covariance functions, we use subscripts to identify the function to which the parameters belong, and the vector
to collect appropriate covariance function parameters. The prior over functions is written(1) 
The critical property of a GP is that the distribution of any finite set of function values has a Gaussian distribution. If the vector
contains the values of the function at times , the prior over is(2) 
where is a covariance matrix constructed from the covariance function: . In the regression setting, we are usually presented with a noise corrupted version of , . Assuming that the noise is Gaussian, writing
(3) 
it is trivial to marginalise the values due to the conjugate nature of the Gaussian noise and the Gaussian process prior:
(4) 
One interpretation of this is that there is a function with Gaussian process prior covariance , which we observe directly. It is this conjugate relationship that we will use to construct structured models, using not i.i.d white noise, but further Gaussian process functions.
Consider a set of time series which we wish to model. We have groups of data taken at times . For example, each group could represent an experimental replication under different conditions. Under our model, there is a latent GP function which governs all the time series, which we denote . Given a draw for , each group of data is then drawn from a GP: . The additional covariance can account for both correlated structure in and noise.
It is possible to marginalise the latent function and thus introduce covariance between the groups, using the same conjugate properties as for the noise discussed above. This covariance amongst the groups then depends on the group index , and we can write a compound covariance function
(5) 
Considering the group index as an input to the function, we can write
(6) 
To obtain a likelihood for the grouped data, we concatenate as and similarly for (constructing ) and the group indexes (). We construct a kernel matrix on the concatenated vectors such that and write
(7) 
where is a vector of all covariance function parameters which we can then infer by typeII maximum likelihood or MCMC sampling.
Mixture models involve the assignment of data into clusters, as well as the inference of the properties of each cluster. The popular Gaussian mixture model involves the assignment of each data vector to one of
Gaussian components, whilst simultaneously inferring the mean and covariance of each of the components. The EM algorithm for Gaussian mixture models is widely known: it treats the assignment of data to clusters as a latent variable problem, and estimates the cluster means and variances by maximising the likelihood. The algorithm alternates between inference of the latent variables and maximisation of the likelihood with respect to the parameters.In Bayesian mixture models, the cluster parameters are treated as random variables, and inference is performed by computing (or approximating) the joint posterior distribution of the latent variables and cluster parameters. A variational approach to inference in a mixture model can be achieved by approximating the posterior with a factorising distribution, which is updated using the VBEM algorithm. Using a set of conjugate priors, the VBEM algorithm alternates between computing the optimal approximating distribution for the assignment variables (the VBE step), and finding that for the cluster parameters and mixing proportions (the VBM step). The variational procedure is better than the EM method since it avoids the pitfalls of maximum likelihood estimation, however both methods require the specification of the number of components,
.Using a Dirichlet process prior for the mixing proportions avoids the problem of selecting the number of components in the model. It also offers a convenient inference procedure via Gibbs sampling, though in this paper we focus on a variational approach. The Dirichlet process can be seen as an infinite mixture model (Rasmussen, 2000): a normal mixture model where the atoms correspond to (parameters of) Gaussian densities, and the number of clusters has been allowed to tend to infinity. In the most general case, DPs can be used for infinite mixture models, not simply Gaussian densities. The mixing proportions of the clusters (and hence the expected number of clusters) are controlled by a concentration parameter .
We propose a Dirichlet Process Gaussian Process (DPGP) mixture model, using the stick breaking construction as follows. Let be a space of functions mapping , and let be a DP on that space, with a GP base distribution and concentration . We draw a series of atoms and associated stickbreaking lengths independently such that
(8) 
From the stick breaking weights we define a series of mixing proportions , and the distribution can be written
(9) 
thus each atom of our DP is a function drawn from a GP. This construction is similar to that provided by Dunson (2010), though we have innovated in using additional structure in the model (which we will show empirically to be very effective), and we also propose a novel inference procedure based on variational Bayes.
To use this construction in clustering functional data, we use the hierarchical GP developed in the previous section. Each group of data is then associated with a single atom of the DP by the variables , and varies from the atom by an independent draw from another GP. Each atom becomes the mean function in a hierarchical GP described above.
In our applications, we have further levels of this hierarchy, with unknown groups (clusters) at the highest level and known groups (replicates) at lower levels. It is simple to extend the model with a series of levels with known and unknown groupings at each, depending on the application.
The generative procedure for our model is then:
Select a concentration parameter
, and GP hyperparameters for the DPGP construct.
Draw an infinite series of stickbreaking lengths and associated atomic functions , compute the infinite mixing proportions .
For each group of data, draw the random variable , thus selecting the atomic function associated with the group.
For each subgroup in that group, draw a function for a GP which defines the deviation from the selected atomic function (this may itself be a hierarchical GP over subsubgroups).
Evaluate the functions at a finite set of points , reporting the values .
Presented with a set of data , perhaps with some known structure, we are tasked with inferring the unknown groupings (via the variables ), the latent functions , the GP parameters and the all the functions which occur in the data.
Variational Bayes is a method for probabilistic inference where the posterior distribution is approximated using some simpler distribution. The usual assumption is that the posterior factorises in some way which yields a tractable lower bound on the marginal likelihood, which then serves as an objective in optimization. The factorising assumption leads naturally to the VBEM algorithm, where each of the factors is updated in turn.
Recently, various forms of collapsed variational Bayes have been proposed for specific models (King and Lawrence, 2006; LázaroGredilla and Titsias, 2011; LázaroGredilla et al., 2011; Teh et al., 2007; Kurihara et al., 2007; Sung et al., 2008), where some of the variables are marginalised analytically. In Hensman et al. (2012b), we showed that many of these schemes are equivalent. We also proposed a Riemann optimization scheme similar to Honkela et al. (2010b), and showed that VBEM is in fact a steepest ascent algorithm upon a Riemannian manifold. Thus introducing geometrically conjugate gradient directions serves to increase the speed of convergence. Here we consider the application of these ideas to the DPGP model.
For Dirichlet process mixture models as we consider here, Kurihara et al. (2007) considered forms of a collapsed stickbreaking prior. Although their approach differs from our derivation in that they marginalised stick breaking lengths before making a variational approximation, we will show that we end up with similar expressions. Further, we collapse all of the parameters from our model aside from the cluster allocation variables. This leads to greater simplification in computing gradients and proposing nongradient moves in optimization such as mergesplit and reordering of the clusters.
The steepest ascent direction on a Riemann manifold is given by the natural gradient (Amari, 1998). For our model, and for related mixture models, we show that the necessary informationgeometric quantities can be computed in closed form without the expensive matrix inverse which hindered previous approaches (such as Honkela et al., 2010b). We demonstrate empirically that our algorithm converges faster than VBEM and the freeform approach of LázaroGredilla and Titsias (2011).
We note that the variables are each of infinite dimension with a single unitary element such that , in our approximate posterior, we shall truncate the number of components, selecting a truncation parameter as in (Blei and Jordan, 2006; Kurihara et al., 2007). To select this we adopt a mergesplit approach which has been applied before in maximum likelihood clustering (Ueda et al., 2000) and also as a MetropolisHastings step in a collapsed Gibbs sampler (Jain and Neal, 2004).
First, we follow the procedure outlined by Hensman et al. (2012b) to derive a collapsed lower bound on the marginal likelihood, which serves as an objective function in optimization.
We briefly formalise the notations for our model, summarising them in Table I. We have a DPGP construction as described above, whose atoms we denote . Suppose we have groups of data which we wish to cluster. In gene expression data, the expression of all genes are necessarily gathered at the same time, so each of the vectors are the same, simplifying our exposition somewhat. Let the values of the function at times be gathered into the vector , and denote the collection of these . Collect the stick breaking lengths similarly . From here, the vectors represent the data that we wish to cluster, and any subgroupings have been concatenated. This simplifies the model as illustrated by Figures 2 and 3. Accordingly, we can define the likelihood in the usual mixture form:
(10) 
where is a covariance matrix constructed according to any subgroups in as per equation (5), the prior for is as defined in equation (8), and the prior for occurs through the usual GP methodology as
(11) 
In the above, the GP hyperparameters have been omitted for brevity: in practise we make point estimates for the parameters interleaving gradientbased optimization with the VB procedure.
symbol  type  description 

Allocates the data group to a latent function  
collection of allocation variables  
the latent function  
realisations of the function at the points  
collection of all realised function values  
the group of observed data  
collection of all observed data  
the stickbreaking length in the DP construction  
collection of all stickbreaking lengths  
mixing proportions defined by stick breaking construction  
concentration parameter of the stick breaking process  
(approximate) posterior probability of assigning the datum to the component 

effective number of assignment to the component in the approximate posterior 
The first step in deriving a collapsed bound is to select which variables should be used for parameterisation, and which should be collapsed from the problem: this can be done using a dseparation test. Given the observed variables (data) and treating the variables we wish to parameterise as observed, the collapsed variables must factorise as in the prior (Hensman et al., 2012b). Examining the graphical representation of our model in Figure 3, we can see that if the latent variables were observed, then the model would dseparate appropriately. We shall use a variational distribution , and introduce the variational parameters and the truncation level such that . To ensure that is a valid distribution, we shall reparameterize it through the softmax function. This will also assist in the computation of natural gradients, as we shall show.
An important difference between our variational method and the VBEM procedure is that we do not introduce parameterisation for the distributions of collapsed variables: they will be analytically marginalised from the problem. The first step is to use Jensen’s inequality to derive a lower bound on the data likelihood, conditioned on the variables that we shall be collapsing, and :
(12) 
Since provides a lower bound on , we have trivially that provides a lower bound on . Thus we also have
(13) 
The second integral is tractable because of the conjugacy between and the prior. We now have a lower bound on the marginal likelihood without specifying any form for the approximate distribution of or .
The integral in (13) separates in and . The integral for follows easily by completing the square and using the Gaussian identity. The integral for is also straightforward, but reveals some relations to previous studies of collapsed stick breaking priors (Kurihara et al., 2007).
The integrals separate as follows:
(14) 
the middle term can be solved as follows. First, the expectation of is trivially
(15) 
where we have defined . The stick breaking lengths all have beta priors with parameters . Since and we have:
(16) 
Substituting these two results back into the main expression, we are left with
(17) 
Note that the part of the prior beyond is trivially marginalised. Substituting the definition for and refactoring gives
(18) 
where we have defined . Finally, we recognise a series of simple betaintegrals and write
(19) 
Note the similarity to the proposed collapsed stick breaking prior of Kurihara et al:
(20) 
where , . In Kurihara et al’s approach, the stick breaking lengths are marginalised from the expression before a variational approximation is made, leading to (20). To make this variational approximation tractable, a first order Taylor expansion of is used around the point . This approximate ’marginalisation’ of (20) leads to a similar expression to (19).
The approximating distribution factorises into a multinomial distributions , with parameters . Since , we use the softmax reparameterisation to avoid constrained optimization: . Denoting the gradient in as , the natural gradient is given by
(21) 
where is the Fisher information matrix of in the parameterisation , which is given by . This matrix is singular due to the overparameterised nature of the softmax function, which makes computing the natural gradient problematic. Kuusela et al. (2009) suggests omitting the first element of , writing , though this can be avoided as follows.
First note that inverse of can be calculated through the ShermanMorrison inversion as
(22) 
where is an appropriately sized matrix of ones. Since , the natural gradient is
(23) 
We note that the symmetry of the softmax parameterisation constrains , thus the gradients are . Taking a step of this length in the variable is equivalent to taking a step of length in the variable , thus the natural gradient can be computed simply by dividing by , with no matrix inversions required.
Since
often contains many elements which are close to zero, this division may cause numerical problems. This can be avoided by considering the chainrule which is used to compute the gradients with respect to
from those for :(24) 
Dividing through by obtains the following expression for the natural gradient, which be find to be stable in computation:
(25) 
This expression for the natural gradient is applicable for the multinomial distribution wherever the softmax parameterisation is used.
We have presented the KLcorrected bound and its natural gradient for a Dirichlet process mixture of Hierarchical Gaussian processes. In the following we discuss some relations between our optimization approach and the VBEM method.
First, it can be shown that the optimal distribution for the collapsed variables () factorises. Our approximate posterior thus takes the same form as the standard mean field approximation.
Next, the meanfield bound can be shown to be exactly the same as the KLcorrected bound. If we set to the same distribution in each, and then update the meanfield bound with a single step of VBEM, the bounds will be equivalent (Hensman et al., 2012b; Sung et al., 2008).
Furthermore, the VBEM procedure is effectively a gradient method (Sato, 2001; Hoffman et al., 2012), taking unit length steps in the natural gradient direction in each of the coordinates in turn (where coordinates correspond to the approximation to each node of the graph). An important result is that the natural gradient on the KLcorrected bound is the same as that on the mean field bound (so long as the other variables are all updated). This means that we can recover exactly the VBEM algorithm by taking unit steps in the natural gradient direction on the KLcorrected bound.
The KLcorrected bound has then brought about the following advantages: it is simpler to derive (there are fewer variables to deal with), and it has a lowerdimensional space for optimization. Surprisingly, this more compact representation does not complicate the optimization: natural gradient steps on the KLcorrected bound have the same effect as a full set of steps on the meanfield bound (or a full round of updates).
Perhaps the most important advantage of the KLcorrected bound is that it enables conjugate gradient descent to be preformed easily. One only needs to consider the conjugate computations (as presented by Honkela et al. (2010b)) in a small number of variables. If the conjugate gradient step fails to improve the bound, then reverting to a unit step in the natural gradient direction will recover the VBEM update again, and the conjugate part of the algorithm can be ’restarted’.
A mergesplit approach has been suggested for mixture models using EM (Ueda et al., 2000) and in MCMC (Jain and Neal, 2004). In this approach, the current solution is reinitialised by either redefining two clusters as one, or one cluster as two new. The collapsed nature of the KLcorrected bound is particularly helpful in performing mergesplit, since we only have the parameters (or equivalently ) to deal with. Since we have a lower bound on the marginal likelihood, we also have a natural method for accepting proposed moves, depending on whether they increase the bound.
To perform a split, we select a cluster component and find the data which are currently associated by examining . We increase the truncation parameter , adding a new cluster, and move half of the probabilistic mass for to this new cluster . After optimising to convergence, we accept the move if the bound increases. We found empirically that merge procedures were not necessary: that optimization naturally managed to merge clusters appropriately. We simply removed empty clusters as appropriate. We also make use of the reordering move (Kurihara et al., 2007) which reorder the solution so that the largest cluster is first: this increases the bound under the DP prior.
We note that the collapsed parameterisation of the model makes these procedures simple to implement: with only a matrix containing to deal with, columns (corresponding to clusters) can be deleted, added, moved or adjusted arbitrarily, so long as the bound on the loglikelihood increases.
We present the application of our model and variational inference procedure to three data sets. In all the experiments, we initialised the allocation of clusters randomly. The effect of the covariance function hyperparameters can play quite a strong role in the results. We initialised using the following rules of thumb: lengthscales were initialised to half of the span of the input data; the toplevel variance (cluster variance) was set to account for 60% of the signal variance; the hierarchical variance was set to 30% and noise was set to account for 10% of the data variance. Optimization of the hyperparameters (against the lower bound on the marginal likelihood) was interleaved with the variational optimization. Unless otherwise stated, the Dirichlet process concentration parameter was fixed to 1.
To demonstrate our model, we generated a synthetic data set as follows. We selected 12 time points randomly in the region , and defined clusters by evaluating the sine function with uniformly randomised phase and randomised frequency around . We randomly selected data per cluster in the interval , and for each datum in each cluster selected a correlated offset from the mean for each cluster using further randomised sine functions, and added a small amount of i.i.d. noise. The data are illustrated in Figure 5.
We present results of clustering data from Kalinka et al. (2010). In this paper, the gene expression of six species of Drosophila was presented, measured at two hour intervals during embryonic development. The data contains a natural structure: aside from structure across species, the experiments was performed in replicate. Pools of embryos were used, with measurements taken every two hours from each pool. Not all of the pools were measured at all time points. We use a hierarchical GP to model this replicate structure, accounting for correlated differences between replicates. We use a further level of the hierarchy for clustering the genes. Using a method similar to Kalaitzis and Lawrence (2011) to eliminate silent genes, we selected 1000 genes for clustering.
An advantage of using GP models for clustering is that we can incorporate prior information into the model. In Gossan et al. (2013) gene expression in mouse cartilage was measured at four hour intervals in duplicate, following a 12h12h lightdark cycle. 304 genes corresponding to circadian rhythms were identified by fitting sine functions to the data. Here we propose a clustering model for these circadian genes. We use a periodic covariance function based on a projection of the Matern covariance (Durrande et al., 2013), drawing functions from the DPGP construct which are periodic (not necessarily sinusoidal) in nature. The next layer of our hierarchical structure uses the standard RBF covariance with i.i.d. noise. In this model, the genes are able to share only a periodic component: deviation from this periodicity is accounted for on a genebygene basis. Clusters inferred by our model are shown in Figure 8. Further analysis of the discovered cluster structure showed that it reflected known groupings of established clock genes as well as providing insights into the regulation of cartilage specific genes by the circadian clock. For more details, see Gossan et al. (2013).
We tested our model on the synthetic data, and compared to GPDP construct without a hierarchical structure, using only i.i.d. noise to model the difference in clustered signals, and a Dirichlet process Gaussian mixture model, which attempts to infer all the covariance structure in the data. For comparison, we set the concentration parameter to give the correct number of clusters a priori, and inferred using the same scheme for each.
Figure 4 shows the ground truth cluster allocation and the inferred cluster allocation for each. From the Figure we see that our model has correctly inferred most of the correct structure: only two true clusters are confused. The GPDP construct without a structured model is poor in finding the clusters: since it cannot account for correlations amongst signals using a noise model, it attempts to introduce extra components into the model to account for variance. The Gaussian mixture model also fails to infer the correct structure: without prior knowledge of signal correlations, it is unable to separate the groups. The clusters inferred by our model are shown in Figure 5.
To optimize the parameters of the covariance functions, we interleave standard our Riemannian method with standard optimization of the covariance functions parameters, keeping the variational distributions fixed. Using the synthetic data, we set up a simple trial to examine the sensitivity of the method to the initial conditions of the covariance function parameters. We created 20 initial conditions, drawing parameters values from the standard lognormal distribution, and optimized the models using the Riemannian procedure interleaved with standard conjugategradient optimization of the parameters and a mergesplit routine. In 16 of the 20 cases, the optimal structure (as shown in Figure
4) was discovered; in the remaining 4 cases, two of the clusters were conflated, which was reflected by a lower bound on the marginal likelihood, whilst the remaining structure was inferred correctly. In all cases, the lengthscales and variances of the covariance functions were estimated correctly, to with 2 decimal places of the best solution.The results of clustering the Drosophila data are shown in the supplementary material. An example of the structure inferred in a single cluster is illustrated in Figure 6. The use of our hierarchical model allows us to find biologically meaningful clusterings that may not be possible without the structure. For example, the first three clusters appear to be very similar: signals in the second cluster rise only slightly faster than in the first. Using the online DAVID tool, we found that the cellular component gene ontology (CCGO) terms were enriched differently in each. In the first cluster, the top CCGO terms were extracelular matrix, extracelluar region part and proteinaceous extracellular matrix, with pvalue belows 0.005. In the second cluster, the top CCGO terms were plasma membrane part and cell junction, with similar pvalues. The third cluster, whose signals arrive slightly later still, was also enriched for genes intrinsic to membrane, but also showed enrichment for cell adhesion and biological adhesion.
For comparative purposes, we also applied our code to clustering the data without the hierarchical structure. This is similar in spirit to that proposed in Dunson (2010), though we maintain our variational framework. In this model, all variation from the cluster mean is modelled as independent Gaussian noise. The result is that many more clusters have to be used to model the data. Noise in the measurement process is not simply i.i.d. as this model enforces: varying sensitivities of the microarray system as well as true biological variation in the genes mean that those genes activated by similar pathways – thus having similar temporal patterns – will vary in a correlated manner. We summarise the results of the two models in Table II.
Hierarchical  Nonhierarchical
(similar to Dunson, 2010) 


N. clusters  52  245 
4256.2  58.9  
Signal variance  1.26  1.41 
Noise variance  0.03  0.04 
hierarchical variance  0.06  n/a 
We first note that the lower bound on the marginal likelihood is dramatically higher for the hierarchical model. Even accounting for the few extra parameters required, the difference is extremely significant. From the table we see that the hierarchical model uses a smaller noise level to model the data, and uses twice the variance of the noise to model hierarchical structure (including both replicate and genewise variance). Plots of the discovered clusters can be found in the supplementary material – it is clear that the nonhierarchical version discovers many small clusters with very similar profiles.
To perform inference, we used the VB procedure described with a mergesplit, and used a gradientbased method to find pointestimates of the kernel parameters, based on maximising .
Subsequently, to compare our Riemann procedure with VBEM and to test the effectiveness of the mergesplit approach, we set the kernel parameters to sensible values found using several optimization runs.
Previous uses of the collapsed VB method have reported that it finds superior solutions (Sung et al., 2008). We empirically found that this occurred in our algorithm: the VBEM procedure, which is a steepest ascent method on the Riemann manifold, becomes stuck on plateaus where there is little gradient. This is illustrated by Figure 7: both the conjugate Riemann method and VBEM pause at the same levels of likelihood, but the nature of our algorithm allows it to pass through this solution, whilst the VBEM algorithm is stuck.
To monitor the effects of this on the ability of the algorithm to find a good solution, we ran 200 restarts for each of our data sets, using the same initial conditions for each without using the mergesplit approach. We then monitored, over all the restarts, how many times the algorithm came to good solution (which we defined as being within 10 nats of the bestfound solution). We divided the total time taken (or iterations used) for all 200 restarts by the number of runs which found such a solution. This statistic then asses how well the algorithms perform in not only speeding up convergence but also escaping the plateaus as discussed. The results are shown in Table III.
Dataset (N. genes)  VBEM iters.  Riemann iters.  VBEM (s)  Riemann (s) 

Synthetic (241)  304  234  1.16  0.90 
Drosophila (600)  680  381  153  88 
mouse cartilage (896)  232  102  34  14 
The effect of expedited convergence using the Riemann approach is amplified when optimising the hyperparameters and using the mergesplit method. In practise, it is necessary to run the variational optimization many times, interleaved between mergesplit trials and optimization of the hyperparameters. Since the Riemann procedure often finds better local solutions for the variational parameters, it does not need to use as many splitprocedures to find a good global optimum. When optimising the hyperparameters interleaved with the VB parameters, the Riemann procedure is also particularly effective, finding local solutions more quickly.
We have presented a method for clustering of structured time series. The method is based on hierarchical Gaussian process. This simple idea allows us to combine related time series groups in a natural way. We introduced a clustering model which allows us to infer the groupings, modelling further structure within subgroups.
Inspired by gene expression time series, we applied our model to three datasets. We showed how the Gaussian process methodology allows us to incorporate knowledge of the problem into the model such as periodicity of the shared timeseries function. The model has many applications in clustering time series, and we are currently exploring the application to motion capture data.
We performed inference in the model using a recent modification of variational Bayes. This not only provided a speed improvement, but also allowed for extremely simple implementations of a mergesplit approach. We related the collapsed expression to that used in collapsed variational Bayes, and showed how to compute the natural gradient for a set of softmaxparameterised variables, a derivation which has wide application in clustering models.
A python implementation of the algorithm and code for running all the experiments can be found on our website at http://staffwww.dcs.shef.ac.uk/people/J.Hensman/ .
Bayesian hierarchical clustering for microarray time series data with replicates and outlier measurements.
BMC Bioinformatics, 12(1):399, 2011.The Journal of Machine Learning Research
, 9999:3235–3268, 2010b.A splitmerge markov chain monte carlo procedure for the dirichlet process mixture model.
Journal of Computational and Graphical Statistics, 13(1):158–182, 2004.Proceedings of the International Joint Conference on Artificial Intelligence
, volume 20, page 19, 2007.Variational heteroscedastic Gaussian process regression.
In Proceedings of the International Conference on Machine Learning (ICML), 2011, 2011.
Comments
There are no comments yet.