Non-Stationary Gaussian Process Regression with Hamiltonian Monte Carlo

Introduction

Gaussian process regression has emerged as a powerful, yet practical class of non-parametric Bayesian models that quantify the uncertainties of the underlying process using Gaussian distributions (Rasmussen and Williams, 2006). Gaussian processes are commonly applied to time-series interpolation, regression and classification, where the GP can provide predictive distributions (Rasmussen and Williams, 2006).

Several authors have proposed extending GPs by learning a latent noise variance as another GP, and by inferring the unknown function and the noise model in a maximum likelihood (ML) (Kersting et al., 2007) or maximum a posteriori (MAP) fashion (Quadrianto et al., 2009). Fully Bayesian inference methods include MCMC sampling (Goldberg, 1998) and variational and expectation propagation approximations of the posterior (Lazaro-Gredilla and Titsias, 2011; Tolvanen et al., 2014). Non-stationarities can also be included in the signal variance or lengthscale with Gaussian process priors. Nonstationary lengthscale was introduced by Gibbs (1997) and further extended by Paciorek and Schervish (2004) with MCMC inference. Recently, Tolvanen et al. (2014) introduced a non-stationary signal variance using expectation propagation and approximate variational inference.

In this paper we introduce the first fully non-stationary and heteroscedastic GP regression framework, in which all three main components (noise variance, signal variance and the lengthscale) can be simultaneously input-dependent, with GP priorsMatlab implementation available from github.com/markusheinonen/adaptivegp. We propose an inference method for the exact joint posterior of the underlying signal and all three latent functions, avoiding the need for introducing variational or expectation propagation approximations (Lazaro-Gredilla and Titsias, 2011; Tolvanen et al., 2014). We use HMC-NUTS, which can effectively sample the posterior guided by the model gradients, which we derive analytically. Furthermore, an exact MAP solution arises as a simple gradient ascent of the posterior. We enhance both approaches by posterior whitening using Cholesky decompositions of the latent function priors. Our experiments demonstrate the necessity of non-stationary GPR to model realistic input-dependent dynamics, while the proposed method performs comparably to conventional stationary or previous non-stationary GPR models otherwise.

In Section 2 we introduce the fully nonstationary GP model. In its subsections we first introduce MAP and HMC inference, discuss model whitening and finally define the predictive distributions. Section 3 presents experimental results on several synthetic and one real biological datasets, and we conclude in Section 4.

Heteroscedatic nonstationary GP model

where both the underlying signal $f(x)$ and the zero-mean observation noise variance $\omega(x)^{2}$ are unknown functions to be learned. We proceed by first placing a zero mean GP prior on the unknown function $f(x)$,

which assumes that the output covariances $\operatorname{\textbf{cov}}(f(x),f(x^{\prime}))=K_{f}(x,x^{\prime})$ depend on the input covariance through a kernel function. We use a nonstationary generalisation of the squared exponential kernel (Gibbs, 1997)

We model the lengthscale, signal variance and noise variance with latent functions. We are interested in smoothly varying latent functions and thus we place separate GP priors on them as well:

where we set the priors on the logarithms to ensure their positivity. We select separate standard squared exponential covariances for each,

We note that by placing a GP prior on just the noise and restricting the other two to be constants, we arrive at the heteroscedastic model studied in several earlier works (Goldberg, 1998; Kersting et al., 2007; Quadrianto et al., 2009). Setting a prior on only the lengthscale retrieves the models of Gibbs (1997); Paciorek and Schervish (2004), and setting a prior on both the signal variance and the noise gives the model of Tolvanen et al. (2014). Out method is the first to combine heteroscedatic noise and nonstationary lengthscale, while also allowing the signal variance to vary over the inputs.

where $\mathbf{f}$ has been marginalised out. Using Bayes’ theorem this is equivalent to maximizing the marginal likelihood

whose logarithm we denote as the marginal log likelihood (MLL).

The partial derivatives of the log of marginal likelihood (3) with respect to the latent functions are analytical:

We perform gradient ascent over the MLL, $\log\mathcal{L}$. The solution is only guaranteed to converge to a local optimum, and hence we perform multiple restarts from random initial conditions. The MAP solution is adequate when the posterior is close to unimodal.

2 HMC inference

As a second approach we sample the full posterior using Hamiltonian Monte Carlo (HMC) (Hoffman and Gelman, 2014; Neal, 2011). In HMC, we introduce an additional momentum variable for each of the model variables and interpret the extended model as a Hamiltonian system. We simulate time evolution of the Hamiltonian dynamics to produce proposals for the Metropolis algorithm. This simulation step makes use of the gradient of the log joint density

3 Posterior whitening

The posterior of the latent vectors is by definition highly correlated due to Gaussian priors, leading to inefficient Monte Carlo sampling. To ease the sampling, we perform the sampling over the whitened latent vectors (Kuss and Rasmussen, 2005)

4 Making predictions

Experiments

We assess the performance of the proposed method on several synthetic and real datasets. We experiment with 8 synthetic datasets and a gene expression time series dataset (Heinonen et al., 2015). Of the synthetic datasets, three datasets are from the literature: the motorcycle dataset M (Silverman, 1985), the ‘jump’ dataset J (Paciorek and Schervish, 2004) and a nonstationary dataset T from GPstuff (demo_epinf in Vanhatalo et al. (2013)). We also generated five synthetic datasets with different types of nonstationarities (See Table 1). We expect datasets exhibiting specific types of input-dependent characteristics to require a model with the corresponding input-dependencies.

With MSE the stationary GP performs slightly better, albeit still worse than the optimal non-stationary method. This applies for every dataset.

Adding ‘unnecessary’ nonstationarities retains or only slightly worsens the performance, with the major exception being the ‘jump’ dataset J. Here, the lengthscale is clearly input-dependent (NLPD $-0.72$, optimal), while in contrast the nonstationary signal variance $\sigma$ is unable to model the data (NLPD $1.04$). Adding Heteroscedatic noise to any of the models of this dataset weakens the model.

2 HMC performance

3 Biological dataset

We showcase the method with a biological dataset of $205$ gene expression time series measurements of human endothelial cells after irradiation at time $t=0$. Due to the irradiation the dataset exhibits nonstationary dynamics as the cells try to repair themselves and revert back to steady states. The gene expressions are measured over 8 days $(0.5,1,2,3,4,7,14,21)$ in three replicates (Heinonen et al., 2015). The goal is to construct a realistic model of the underlying gene expression process and the underlying dynamics with no knowledge of the ‘true’ expression levels, given only the small number of sparse measurements.

4 Latent function reconstruction

Figure 5 bottom highlights the MAP and HMC solutions given $150$ datapoints, and compares them to the state-of-the-art $\sigma,\omega$-GP model of Tolvanen et al. (2014). We are able to accurately estimate the latent functions.

Discussion

In this paper, we have proposed a fully non-stationary Gaussian process regression framework, where all three key components can be input-dependent. Our approach uses analytical gradient-based techniques to perform inference with HMC sampling and MAP estimation. We are able to effectively sample from the exact posterior of the latent functions. We have shown that the method is able to infer the underlying latent functions and improve regression performance when the datasets truly are nonstationary, and achieve equivalent performance to a stationary model when they are not.

The interplay between the signal variance and the lengthscale is an interesting topic (Diggle et al., 1998; Zhang, 2004). When modeling the ‘jump’ dataset the signal variance was unable to model dynamics, while nonstationary lengthscale produced a good model. This is natural since the signal variance serves as a linear amplitude, while the lengthscale has a possibly non-linear effect on the function model.

The gradient-based HMC is a powerful inference tool for Gaussian processes, and could be further enhanced by utilizing natural gradients or position-dependent mass matrices with Riemannian Manifold HMC (Girolami and Calderhead, 2011). We note that the method could be extended by also inferring the hyperparameters $\boldsymbol{\theta}$ using HMC. However, proper care has to be taken to set their priors.

Supplemental information

Kernel SDP proof

is a one-dimensional Gaussian SDP kernel $C^{NS}(\mathbf{x}_{i},\mathbf{x}_{j})$ of Paciorek and Schervish (2004). The kernel $K_{f}$ can be stated as

which is positive definite for any function $\sigma(\cdot)$ (Shawe-Taylor and Christianini, 2004).

Conditional distributions

The conditional distributions of the latent functions at target timepoints $\mathbf{x}_{\star}$ given the latent functions at observed timepoints $\mathbf{x}$ are

Partial derivatives

The partial derivatives of the unconstrained latent functions against the marginal log likelihood

where $\boldsymbol{\alpha}=K_{y}^{-1}\mathbf{y}$.

References