Slice sampling covariance hyperparameters of latent Gaussian models

Introduction

Many probabilistic models incorporate multivariate Gaussian distributions to explain dependencies between variables. Gaussian process (GP) models and generalized linear mixed models are common examples. For non-Gaussian observation models, inferring the parameters that specify the covariance structure can be difficult. Existing computational methods can be split into two complementary classes: deterministic approximations and Monte Carlo simulation. This work presents a method to make the sampling approach easier to apply.

In recent work Murray et al. developed a slice sampling variant, elliptical slice sampling, for updating strongly coupled a-priori Gaussian variates given non-Gaussian observations. Previously, Agarwal and Gelfand demonstrated the utility of slice sampling for updating covariance parameters, conventionally called hyperparameters, with a Gaussian observation model, and questioned the possibility of slice sampling in more general settings. In this work we develop a new slice sampler for updating covariance hyperparameters. Our method uses a robust representation that should work well on a wide variety of problems, has very few technical requirements, little need for tuning and so should be easy to apply.

We consider generative models of data that depend on a vector of latent variables $\mathbf{f}$ that are Gaussian distributed with covariance $\Sigma_{\theta}$ set by unknown hyperparameters $\theta$. These models are common in the machine learning Gaussian process literature [e.g. 4] and throughout the statistical sciences. We use standard notation for a Gaussian distribution with mean $\mathbf{m}$ and covariance $\Sigma$,

and use ${\mathbf{f}\sim\mathcal{N}(\mathbf{m},\Sigma)}$ to indicate that $\mathbf{f}$ is drawn from a distribution with the density in (1).

The generic form of the generative models we consider is summarized by

The methods discussed in this paper apply to covariances $\Sigma_{\theta}$ that are arbitrary positive definite functions parameterized by $\theta$. However, our experiments focus on the popular case where the covariance is associated with $N$ input vectors $\{\mathbf{x}_{n}\}_{n=1}^{N}$ through the squared-exponential kernel,

For non-Gaussian likelihoods we wish to sample from the joint posterior over unknowns,

We would like to avoid implementing new code or tuning algorithms for different covariances $\Sigma_{\theta}$ and conditional likelihood functions $\mathcal{L}(\mathbf{f})$.

Markov chain inference

Recent work has focused on transition operators for updating the latent variables $\mathbf{f}$ given data and a fixed covariance $\Sigma_{\theta}$ . Updates to the hyperparameters for fixed latent variables $\mathbf{f}$ need to leave the conditional posterior,

invariant. The simplest algorithm for this is the Metropolis–Hastings operator, see Algorithm 1. Other possibilities include slice sampling and Hamiltonian Monte Carlo .

Alternately fixing the unknowns $\mathbf{f}$ and $\theta$ is appealing from an implementation standpoint. However, the resulting Markov chain can be very slow in exploring the joint posterior distribution. Figure 1a shows latent vector samples using squared-exponential covariances with different lengthscales. These samples are highly informative about the lengthscale hyperparameter that was used, especially for short lengthscales. The sharpness of $P(\theta\kern 1.0pt|\kern 1.0pt\mathbf{f})$, Figure 1b, dramatically limits the amount that any Markov chain can update the hyperparameters $\theta$ for fixed latent values $\mathbf{f}$.

Often the conditional likelihood is quite weak; this is why strong prior smoothing assumptions are often introduced in latent Gaussian models. In the extreme limit in which there is no data, i.e. $\mathcal{L}$ is constant, the target distribution is the prior model, $P(\mathbf{f},\theta)\!=\!\mathcal{N}(\mathbf{f};0,\Sigma_{\theta})\,p_{h}\kern-1.0pt(\theta)$. Sampling from the prior should be easy, but alternately fixing $\mathbf{f}$ and $\theta$ does not work well because they are strongly coupled. One strategy is to reparameterize the model so that the unknown variables are independent under the prior.

Independent random variables can be identified from a commonly-used generative procedure for the multivariate Gaussian distribution. A vector of independent normals, $\boldsymbol{\nu}$, is drawn independently of the hyperparameters and then deterministically transformed:

We can choose to update the hyperparameters $\theta$ for fixed $\boldsymbol{\nu}$ instead of fixed $\mathbf{f}$. As the original latent variables $\mathbf{f}$ are deterministically linked to the hyperparameters $\theta$ in (5), these updates will actually change both $\theta$ and $\mathbf{f}$. The samples in Figure 1a resulted from using the same whitened variable $\boldsymbol{\nu}$ with different hyperparameters. They follow the same general trend, but vary over the lengthscales used to construct them.

The posterior over hyperparameters for fixed $\boldsymbol{\nu}$ is apparent by applying Bayes rule to the generative procedure in (5), or one can laboriously obtain it by changing variables in (3):

Algorithm 2 is the Metropolis–Hastings operator for this distribution. The acceptance rule now depends on the latent variables through the conditional likelihood $\mathcal{L}(\mathbf{f})$ instead of the prior $\mathcal{N}(\mathbf{f};0,\Sigma_{\theta})$ and these variables are automatically updated to respect the prior. In the no-data limit, new hyperparameters proposed from the prior are always accepted.

Surrogate data model

Neither of the previous two algorithms are ideal for statistical applications, which is illustrated in Figure 2. Algorithm 2 is ideal in the “weak data” limit where the latent variables $\mathbf{f}$ are distributed according to the prior. In the example, the likelihoods are too restrictive for Algorithm 2’s proposal to be acceptable. In the “strong data” limit, where the latent variables $\mathbf{f}$ are fixed by the likelihood $\mathcal{L}$, Algorithm 1 would be ideal. However, the likelihood terms in the example are not so strong that the prior can be ignored.

We augment the latent Gaussian model with auxiliary variables, $\mathbf{g}$, a noisy version of the true latent variables:

For now $S_{\theta}$ is an arbitrary free parameter that could be set by hand to either a fixed value or a value that depends on the current hyperparameters $\theta$. We will discuss how to automatically set the auxiliary noise covariance $S_{\theta}$ in Section 3.2.

The original model, $\mathbf{f}\sim\mathcal{N}(0,\Sigma_{\theta})$ and (7) define a joint auxiliary distribution $P(\mathbf{f},\mathbf{g}\kern 1.0pt|\kern 1.0pt\theta)$ given the hyperparameters. It is possible to sample from this distribution in the opposite order, by first drawing the auxiliary values from their marginal distribution

and then sampling the model’s latent values conditioned on the auxiliary values from

That is, under the auxiliary model the latent variables of interest are drawn from their posterior given the surrogate data $\mathbf{g}$. Again we can describe the sampling process via a draw from a spherical Gaussian:

We then condition on the “whitened” variables $\boldsymbol{\eta}$ and the surrogate data $\mathbf{g}$ while updating the hyperparameters $\theta$. The implied latent variables $\mathbf{f}(\theta,\boldsymbol{\eta},\mathbf{g})$ will remain a plausible draw from the surrogate posterior for the current hyperparameters. This is illustrated in Figure 2.

We can leave the joint distribution (3) invariant by updating the following conditional distribution derived from the above generative model:

The Metropolis–Hastings Algorithm 3 contains a ratio of these terms in the acceptance rule.

The Metropolis–Hastings algorithms discussed so far have a proposal distribution $q(\theta^{\prime};\theta)$ that must be set and tuned. The efficiency of the algorithms depend crucially on careful choice of the scale $\sigma$ of the proposal distribution. Slice sampling is a family of adaptive search procedures that are much more robust to the choice of scale parameter.

Algorithm 4 applies one possible slice sampling algorithm to a scalar hyperparameter $\theta$ in the surrogate data model of this section. It has a free parameter $\sigma$, the scale of the initial proposal distribution. However, careful tuning of this parameter is not required. If the initial scale is set to a large value, such as the width of the prior, then the width of the proposals will shrink to an acceptable range exponentially quickly. Stepping-out procedures could be used to adapt initial scales that are too small. We assume that axis-aligned hyperparameter moves will be effective, although reparameterizations could improve performance [e.g. 9].

For likelihood terms that factorize, $\mathcal{L}(\mathbf{f})\!=\!\prod_{i}\mathcal{L}_{i}(f_{i})$, we can measure how much the likelihood restricts each variable individually:

A Gaussian can be fitted by moment matching or a Laplace approximation (matching second derivatives at the mode). Such fits, or close approximations, are often possible analytically and can always be performed numerically as the distribution is only one-dimensional. Given a Gaussian fit to the site-posterior (12) with variance $v_{i}$, we can set the auxiliary noise to a level that would result in the same posterior variance at that site alone: ${(S_{\theta})_{ii}\!=\!(v_{i}^{-1}\!-\!{(\Sigma_{\theta})_{ii}}^{-1})^{-1}}$. (Any negative $(S_{\theta})_{ii}$ must be thresholded.) The moment matching procedure is a grossly simplified first step of “assumed density filtering” or “expectation propagation” , which are too expensive for our use in the inner-loop of a Markov chain.

Related work

We have discussed samplers that jointly update strongly-coupled latent variables and hyperparameters. The hyperparameters can move further in joint moves than their narrow conditional posteriors (e.g., Figure 1b) would allow. A generic way of jointly sampling real-valued variables is Hamiltonian/Hybrid Monte Carlo (HMC) . However, this method is cumbersome to implement and tune, and using HMC to jointly update latent variables and hyperparameters in hierarchical models does not itself seem to improve sampling .

Christensen et al. have also proposed a robust representation for sampling in latent Gaussian models. They use an approximation to the target posterior distribution to construct a reparameterization where the unknown variables are close to independent. The approximation replaces the likelihood with a Gaussian form proportional to $\mathcal{N}(\mathbf{f};\hat{\mathbf{f}},\Lambda(\hat{\mathbf{f}}))$:

where $\Lambda$ is often diagonal, or it was suggested one would only take the diagonal part. This Taylor approximation looks like a Laplace approximation, except that the likelihood function is not a probability density in $\mathbf{f}$. This likelihood fit results in an approximate Gaussian posterior $\mathcal{N}(\mathbf{f};\mathbf{m}_{\theta\kern-1.0pt,\mathbf{g}=\hat{\mathbf{f}}},R_{\theta})$ as found in (9), with noise $S_{\theta}\!=\!\Lambda(\hat{\mathbf{f}})^{-1}$ and data $\mathbf{g}\!=\!\hat{\mathbf{f}}$.

When the likelihood is Gaussian, the reparameterized variables $\boldsymbol{\omega}$ are independent of each other and the hyperparameters. The hope is that approximating non-Gaussian likelihoods will result in nearly-independent parameterizations on which Markov chains will mix rapidly.

Taylor expanding some common log-likelihoods around the maximum is not well defined, for example approximating probit or logistic likelihoods for binary classification, or Poisson observations with zero counts. These Taylor expansions could be seen as giving flat or undefined Gaussian approximations that do not reweight the prior. When all of the likelihood terms are flat the reparameterization approach reduces to that of Section 2.1. The alternative $S_{\theta}$ auxiliary covariances that we have proposed could be used instead.

The surrogate data samplers of Section 3 can also be viewed as using reparameterizations, by treating $\boldsymbol{\eta}\!=\!L_{\kern-1.0ptR_{\theta}}^{-1}(\mathbf{f}-\mathbf{m}_{\theta\kern-1.0pt,\mathbf{g}})$ as an arbitrary random reparameterization for making proposals. A proposal density $q(\boldsymbol{\eta}^{\prime},\theta^{\prime};\boldsymbol{\eta},\theta)$ in the reparameterized space must be multiplied by the Jacobian $|L_{\kern-1.0ptR_{\theta^{\prime}}}^{-1}|$ to give a proposal density in the original parameterization. The probability of proposing the reparameterization must also be included in the Metropolis–Hastings acceptance probability:

A few lines of linear algebra confirms that, as it must do, the same acceptance ratio results as before. Alternatively, substituting (3) into (15) shows that the acceptance probability is very similar to that obtained by applying Metropolis–Hastings to (14) as proposed by Christensen et al. . The differences are that the new latent variables $\mathbf{f}^{\prime}$ are computed using different pseudo-posterior means and the surrogate data method has an extra term for the random, rather than fixed, choice of reparameterization.

The surrogate data sampler is easier to implement than the previous reparameterization work because the surrogate posterior is centred around the current latent variables. This means that 1) no point estimate, such as the maximum likelihood $\hat{\mathbf{f}}$, is required. 2) picking the noise covariance $S_{\theta}$ poorly may still produce a workable method, whereas a fixed reparameterized can work badly if the true posterior distribution is in the tails of the Gaussian approximation. Christensen et al. pointed out that centering the approximate Gaussian likelihood in their reparameterization around the current state is tempting, but that computing the Jacobian of the transformation is then intractable. By construction, the surrogate data model centers the reparameterization near to the current state.

Experiments

We empirically compare the performance of the various approaches to GP hyperparameter sampling on four data sets: one regression, one classification, and two Cox process inference problems. Further details are in the rest of this section, with full code as supplementary material. The results are summarized in Figure 3 followed by a discussion section.

In each of the experimental configurations, we ran ten independent chains with different random seeds, burning in for 1000 iterations and sampling for 5000 iterations. We quantify the mixing of the chain by estimating the effective number of samples of the complete data likelihood trace using R-CODA , and compare that with three cost metrics: the number of hyperparameter settings considered (each requiring a small number of covariance decompositions with $O(n^{3})$ time complexity), the number of likelihood evaluations, and the total elapsed time on a single core of an Intel Xeon 3GHz CPU.

The experiments are designed to test the mixing of hyperparameters $\theta$ while sampling from the joint posterior (3). All of the discussed approaches except Algorithm 1 update the latent variables $\mathbf{f}$ as a side-effect. However, further transition operators for the latent variables for fixed hyperparameters are required. In Algorithm 2 the “whitened” variables $\boldsymbol{\nu}$ remain fixed; the latent variables and hyperparameters are constrained to satisfy $\mathbf{f}\!=\!L_{\kern-1.0pt\Sigma_{\theta}}\boldsymbol{\nu}$. The surrogate data samplers are ergodic: the full joint posterior distribution will eventually be explored. However, each update changes the hyperparameters and requires expensive computations involving covariances. After computing the covariances for one set of hyperparameters, it makes sense to apply several cheap updates to the latent variables. For every method we applied ten updates of elliptical slice sampling to the latent variables $\mathbf{f}$ between each hyperparameter update. One could also consider applying elliptical slice sampling to a reparameterized representation, for simplicity of comparison we do not. Independently of our work Titsias has used surrogate data like reparameterizations to update latent variables for fixed hyperparameters.

We implemented six methods for updating Gaussian covariance hyperparameters. Each method used the same slice sampler, as in Algorithm 4, applied to the following model representations. fixed: fixing the latent function $\mathbf{f}$ . prior-white: whitening with the prior. surr-site: using surrogate data with the noise level set to match the site posterior (12). We used Laplace approximations for the Poisson likelihood. For classification problems we used moment matching, because Laplace approximations do not work well . surr-taylor: using surrogate data with noise variance set via Taylor expansion of the log-likelihood (13). Infinite variances were truncated to a large value. post-taylor and post-site: as for the surr- methods but a fixed reparameterization based on a posterior approximation (14).

We evaluated four different methods for performing binary GP classification: fixed, prior-white, surr-site and post-site. We applied these methods to the Ionosphere dataset , using 200 training data and 34 dimensions. We used a logistic likelihood with zero-mean prior, inferring lengthscales as well as signal variance. The -taylor methods reduce to other methods or don’t apply because the maximum of the log-likelihood is at plus or minus infinity.

We tested all six methods on an inhomogeneous Poisson process with a Gaussian process prior for the log-rate. We sampled the hyperparameters in (2) and a mean offset to the log-rate. The model was applied to two point process datasets: 1) a record of mining disasters with 191 events in 112 bins of 365 days. 2) 195 redwood tree locations in a region scaled to the unit square split into $25\!\times\!25\!=\!625$ bins. The results for the mining problem were initially highly variable. As the mining experiments were also the quickest we re-ran each chain for 20,000 iterations.

Discussion

On the Ionosphere classification problem both of the -site methods worked much better than the two baselines. We slightly prefer surr-site as it involves less problem-specific derivations than post-site.

On the synthetic test the post- and surr- methods perform very similarly. We had expected the existing post- method to have an advantage of perhaps up to 2–3$\times$, but that was not realized on this particular dataset. The post- methods had a slight time advantage, but this is down to implementation details and is not notable.

On the mining problem the Poisson likelihoods are often close to Gaussian, so the existing post-taylor approximation works well, as do all of our new proposed methods. The Gaussian approximations to the Poisson likelihood fit most poorly to sites with zero counts. The redwood dataset discretizes two-dimensional space, leading to a large number of bins. The majority of these bins have zero counts, many more than the mining dataset. Taylor expanding the likelihood gives no likelihood contribution for bins with zero counts, so it is unsurprising that post-taylor performs similarly to prior-white. While surr-taylor works better, the best results here come from using approximations to the site-posterior (12). For unreasonably fine discretizations the results can be different again: the site- reparameterizations do not always work well.

Our empirical investigation used slice sampling because it is easy to implement and use. However, all of the representations we discuss could be combined with any other MCMC method, such as recently used for Cox processes. The new surrogate data and post-site representations offer state-of-the-art performance and are the first such advanced methods to be applicable to Gaussian process classification.

An important message from our results is that fixing the latent variables and updating hyperparameters according to the conditional posterior — as commonly used by GP practitioners — can work exceedingly poorly. Even the simple reparameterization of “whitening the prior” discussed in Section 2.1 works much better on problems where smoothness is important in the posterior. Even if site approximations are difficult and the more advanced methods presented are inapplicable, the simple whitening reparameterization should be given serious consideration when performing MCMC inference of hyperparameters.

We thank an anonymous reviewer for useful comments. This work was supported in part by the IST Programme of the European Community, under the PASCAL2 Network of Excellence, IST-2007-216886. This publication only reflects the authors’ views. RPA is a junior fellow of the Canadian Institute for Advanced Research.

References