Gaussian Process Regression Networks

Introduction

Gaussian process models have become exceptionally popular for solving non-linear regression and classification problems. They are expressive, interpretable, avoid over-fitting, and have impressive predictive performance in many thorough empirical comparisons (Rasmussen, , 1996; Kuss and Rasmussen, , 2005; Rasmussen and Williams, , 2006).

In machine learning, Gaussian process regression developed out of neural networks research. Neal, (1996) showed that Bayesian neural networks became Gaussian processes as the number of hidden units approached infinity, and conjectured that “there may be simpler ways to do inference in this case.” These simple inference techniques became the cornerstone of subsequent Gaussian process models. However, neural networks had been motivated in part by their ability to capture correlations between multiple outputs (responses), by using adaptive hidden units that were shared between the outputs. In the infinite limit, this ability was lost.

Recently there has been an explosion of interest in extending the Gaussian process regression framework to account for fixed correlations between output variables (Alvarez and Lawrence, , 2011; Yu et al., , 2009; Alvarez and Lawrence, , 2008; Bonilla et al., , 2008; Osborne et al., , 2008; Teh et al., , 2005; Boyle and Frean, , 2004). These are often called ‘multi-task’ learning or ‘multiple output’ regression models. Capturing correlations between outputs (response variables) can be used to make better predictions. Imagine we wish to predict cadmium concentrations in a region of the Swiss Jura, where geologists are interested in heavy metal concentrations. A standard Gaussian process regression model would only be able to use cadmium training measurements. With a multi-task method, we can also make use of correlated heavy metal measurements to enhance cadmium predictions (Goovaerts, , 1997). We could further enhance predictions if we make use of how these (signal) correlations change with geographical location.

There has similarly been great interest in extending Gaussian process (GP) regression to account for input dependent noise variances (Goldberg et al., , 1998; Kersting et al., , 2007; Adams and Stegle, , 2008; Turner and Sahani, , 2008; Turner, , 2010; Wilson and Ghahramani, 2010a, ; Wilson and Ghahramani, 2010b, ; Lázaro-Gredilla and Titsias, , 2011). Wilson and Ghahramani, 2010b ; Wilson and Ghahramani, (2011) and Fox and Dunson, (2011) further extended the GP framework to accommodate input dependent noise correlations between multiple output (response) variables.

Other extensions include Gaussian process regression with non-stationary covariance function amplitudes (Turner and Sahani, , 2008; Adams and Stegle, , 2008) and length-scales (Gibbs, , 1997; Schmidt and O’Hagan, , 2003), and with heavy tailed predictive distributions (Neal, , 1997; Vanhatalo et al., , 2009) for outlier rejection (De Finetti, , 1956; Dawid, , 1973; O’Hagan, , 1979).

In this paper, we introduce a new regression framework, Gaussian Process Regression Networks (GPRN), which combines the structural properties of Bayesian neural networks (Neal, , 1996) with the nonparametric flexibility of Gaussian processes. This network is an adaptive mixture of Gaussian processes, which naturally accommodates input dependent signal and noise correlations between multiple output variables, input dependent length-scales and amplitudes, and heavy tailed predictive distributions, without expensive or numerically unstable computations.

We start by introducing the GPRN framework, and show how to perform efficient inference using both Markov chain Monte Carlo (MCMC) and variational Bayes (VB). Carefully following Alvarez and Lawrence, (2011), we compare to eight multiple output GP models on gene expression and geostatistics datasets. We then compare to multivariate volatility models on several benchmark financial datasets, following Wilson and Ghahramani, 2010b . In the Appendix, we review Gaussian process regression and the notation of Rasmussen and Williams, (2006).

Gaussian Process Regression Networks

where $\bm{\epsilon}$ and $\bm{z}$ are i.i.d. $\mathcal{N}(0,I)$ white noiseThe distribution of $\bm{z}$ could be Student-$t$, Laplace, or something different. Using diagonal noise would also be a straightforward extension., $W(x)$ is a $p\times q$ matrix of independent Gaussian processes such that $W(x)_{ij}\sim\mathcal{GP}(0,k_{w})$, and $\bm{f}(x)=(f_{1}(x),\dots,f_{q}(x))^{\top}$ is a $q\times 1$ vector of independent GPs with $f_{i}(x)\sim\mathcal{GP}(0,k_{f_{i}})$. Each of the latent Gaussian processes in $\bm{f}(x)$ have additive Gaussian noise. Changing variables to include the noise $\sigma_{f}\bm{\epsilon}$ we let $\hat{f}_{i}(x)\sim\mathcal{GP}(0,k_{\hat{f}_{i}})$, where

and $\delta_{xx^{\prime}}$ is the Kronecker delta.

We represent this Gaussian process regression network (GPRN) in Figure 1, labelling the length-scale hyperparameters for the kernels $k_{w}$ and $k_{f}$ as $\bm{\theta}_{w}$ and $\bm{\theta}_{f}$ respectively. We see the latent node functions $\hat{\bm{f}}(x)$ are connected together to form the outputs $\bm{y}(x)$. The strengths of the connections change as a function of $x$; the weights themselves – the entries of $W(x)$ – are functions. Old connections can break and new connections can form. This is an adaptive network, where the signal and noise correlations between the components of $\bm{y}(x)$ vary with $x$. Coincidentally, there is an unrelated paper called “Gaussian process networks” (Friedman and Nachman, , 2000), which is about learning the structure of Bayesian networks – e.g. the direction of dependence between random variables.

To explicitly separate the dynamic signal and noise correlations, we re-write (1) as

Conditioning on $W(x)$ in (3), we can better understand the signal correlations. In this case, each of the outputs $\bm{y}_{i}(x)$, $i=1,\dots,p$, is a Gaussian process with kernel

Even if $\sigma_{f}^{2}$ and $\sigma_{y}^{2}$ are zero, so that this is a noise free regression, there are still signal correlations; the components of $\bm{y}$ are coupled through the matrix $W(x)$. Once the network has been trained, $W(x)$ is conditioned on the data $\mathcal{D}$, and so the predictive covariances of $\bm{y}(x_{*})|\mathcal{D}$ are now influenced by the values of the observations themselves, and not just distances between the test point $x_{*}$ and the observed points $x_{1},\dots,x_{N}$ as is the case for independent GPs; we can view (4) as an adaptive kernel learned from the data. There are three other interesting features in equation (4): 1) the amplitude of the covariance function $\sum_{j=1}^{q}W_{ij}(x)W_{ij}(x^{\prime})$ is non-stationary (input dependent); 2) even if each of the kernels $k_{f_{j}}$ has different stationary length-scales, the mixture of the kernels $k_{f_{j}}$ is input dependent and so the effective overall length-scale is non-stationary; 3) the kernels $k_{f_{j}}$ may be entirely different: some may be periodic, others squared exponential, others Brownian motion, etc. . So the overall covariance function (kernel) may be continuously switching between regions of entirely different covariance structures.

In addition to modelling signal correlations, we can see from equation (3) that the GPRN is simultaneously a multivariate volatility model. The noise covariance is $\sigma_{f}^{2}W(x)W(x)^{\top}+\sigma_{y}^{2}I$. Since the entries of $W(x)$ are GPs, this noise model is an example of a generalised Wishart process (Wilson and Ghahramani, 2010b, ; Wilson and Ghahramani, , 2011).

The number of nodes $q$ influences how the model accounts for signal and noise correlations. If $q$ is smaller than $p$, the dimension of $\bm{y}(x)$, the model performs dimensionality reduction and matrix factorization as part of the regression on $\bm{y}(x)$ and $\text{cov}[\bm{y}(x)]$. However, we may want $q>p$, for instance if the output space were one dimensional ($p=1$). In this case we would need $q>1$ to realise features 2 and 3 listed above. For a given dataset, we can vary $q$ and select the value which gives the highest marginal likelihood on training data. We can also use ‘automatic relevance determination’ (MacKay and Neal, , 1994) as a proxy for model selection for $q$ for a given dataset. This is achieved by introducing $\{a_{j}\}$, signal variances for each node function $j$, so that $k_{\hat{f}_{j}}\to a_{j}k_{\hat{f}_{j}}$, and comparing magnitudes of the trained $a_{j}$.

When $q=p=1$, the GPRN essentially becomes the nonstationary GP regression model of Adams and Stegle, (2008) and Turner and Sahani, (2008). Likewise, when the weight functions are constants the GPRN becomes the Semiparametric Latent Factor Model (SLFM) of Teh et al., (2005), except that the resulting GP regression network is less prone to over-fitting through its use of full Bayesian inference.In Teh et al., (2005), the weight constants are a large matrix of hyperparameters, determined through maximising a marginal likelihood. Indeed, in an implementation of GPRN, one can switch features on or off; to switch off changing correlations and multivariate volatility, set $\sigma_{f}^{2}=0$ and the length-scales for weight function kernels ($k_{w}$) to large fixed values.

Inference

As a first step, we re-write the prior in terms of $\bm{u}=(\hat{\mathbf{f}},\mathbf{W})$, a vector composed of all the node and weight Gaussian process functions, evaluated at the training points $\{x_{1},\dots,x_{N}\}$. There will be $q$ node functions and $p\times q$ weight functions. Therefore

where $C_{B}$ is an $Nq(p+1)\times Nq(p+1)$ block diagonal matrix, since the weight and node functions are independent in the prior. The way we have ordered $\bm{u}$, the first $q$ blocks are $N\times N$ covariance matrices $K_{\hat{f}}$ from the node kernel $k_{\hat{f}}$, and the last blocks are $N\times N$ covariance matrices $K_{w}$ from the weight kernel $k_{w}$.

Next we specify our likelihood function, so we can use Bayes’ theorem to find the posterior $p(\bm{u}|\mathcal{D},\bm{\gamma})$. From (1), our likelihood is

By incorporating noise on $\bm{f}$, the GP network accounts for multivariate volatility (as in (3)), without the need for costly or numerically unstable matrix inversions. For other multivariate volatility models, like multivariate GARCH (Bollerslev et al., , 1988), or multivariate stochastic volatility (Harvey et al., , 1994), the likelihood takes the form $p(\mathcal{D}|\beta)=\prod_{i=1}^{N}\mathcal{N}(\bm{\mu}_{i},\Sigma_{i})$, and requires inversions of $p\times p$ covariance matrices. There are three other notable advantages to the inference with GPRN: 1) it is easy to simultaneously estimate $\bm{\mu}_{i}$ and $\Sigma_{i}$. Usually in the multivariate volatility setting, $\bm{\mu}_{i}$ is assumed to be a constant; 2) we can use a Student-$t$ observation model instead of a Gaussian observation model, by letting $\bm{z}$ in (1) be $t$ distributed, with minimal changes to the inference procedures; 3) we can transform the components of the product $W(x_{i})\hat{\bm{f}}(x_{i})$ so that the priors on the components of $\bm{y}(x_{i})$ become copula processes (Wilson and Ghahramani, 2010a, ) and have whatever marginals we desire. We can also do this without significantly changing inference procedures.

Now that we have specified our prior and likelihood, we can apply Bayes’ theorem:

In the next sections, we discuss how to either sample from or use variational Bayes to approximate this posterior in (7), so that we can estimate $p(\bm{y}(x_{*})|\mathcal{D})$. We also use variational Bayes to learn the hyperparameters $\bm{\gamma}$.

To sample from (7), we could use a Gibbs sampling scheme which would have conjugate posterior updates, alternately conditioning on weight and node functions. However, this Gibbs cycle would mix poorly because of the tight correlations between the weights and the nodes. In general, MCMC samples from (7) mix poorly because of the strong correlations in the prior imposed by $C_{B}$. The sampling process is also often slowed by costly matrix inversions in the likelihood.

We use Elliptical Slice Sampling (Murray et al., , 2010), a recent MCMC technique specifically designed to sample from posteriors with tightly correlated Gaussian priors. It does joint updates and has no free parameters. We find that it mixes well. And since there are no costly or numerically unstable matrix inversions in the likelihood of (6) we also find sampling to be highly efficient.

With a sample from (7), we can sample from the predictive $p(W(x_{*}),{\bm{f}}(x_{*})|\bm{u},\sigma_{f},\mathcal{D})$. Let $W_{*}^{i},{\bm{f}}_{*}^{i}$ be the $i^{\text{th}}$ such joint sample. Using (3) we can then construct samples of $p(\bm{y}(x_{*})|W_{*}^{i},\bm{f}_{*}^{i},\sigma_{f},\sigma_{y})$, from which we can construct the predictive distribution

We see that even with a Gaussian observation model, the predictive distribution in (8) is an infinite mixture of Gaussians, and will generally be heavy tailed and therefore robust to outliers.

Mixing was assessed by looking at trace plots of samples, and the likelihoods of these samples. Specific information about how long it takes to sample a solution for a given problem is in the experiments section.

2 Variational Bayes

We perform variational EM (Jordan et al., , 1999) to fit an approximate posterior $q$ to the true posterior $p$, by minimising the Kullback-Leibler divergence $KL(q||p)=-H[q(\mathbf{v})]-\int q(\mathbf{v})\log p(\mathbf{v})d\mathbf{v},$ where $H[q(\mathbf{v})]=-\int q(\mathbf{v})\log q(\mathbf{v})d\mathbf{v}$ is the entropy and $\mathbf{v}=\{\mathbf{f},\mathbf{W},\sigma^{2}_{f},\sigma^{2}_{y},a_{j}\}$.

We use Variational Message Passing (Winn and Bishop, , 2006) under the Infer.NET framework (Minka et al., , 2010) to estimate the posterior over $\mathbf{v}=\{\mathbf{f},\mathbf{W},\sigma^{2}_{f},\sigma^{2}_{y},a_{j}\}$. We specify inverse Gamma priors on $\{\sigma^{2}_{f},\sigma^{2}_{y},a_{j}\}$:

For mathematical and computational convenience we introduce the following variables which are deterministic functions of the existing variables in the model:

Note that the observations $y_{i}(x_{n})\sim\mathcal{N}(s_{in},\sigma_{y}^{2})$ and that $\hat{f}_{nj}\sim\mathcal{N}(f^{\prime}_{nj},\sigma^{2}_{f_{j}})$. Variational message passing uses these deterministic factors and the associated “pseudo-marginals” as conduits to pass appropriate moments, resulting in the same updates as standard VB (Winn and Bishop, , 2006). The full model can now be written as

We use a variational posterior of the following form:

where $q_{\sigma^{2}_{y}},q_{\sigma^{2}_{fj}}$ and $q_{a_{j}}$ are inverse Gamma distributions; $q_{w_{nij}},q_{f^{\prime}_{nj}},q_{\hat{f}_{nj}},q_{t_{nij}}$ and $q_{s_{in}}$ are univariate normal distributions; and $q_{\mathbf{f}_{j}}$ and $q_{\mathbf{W}_{ij}}(\mathbf{W}_{ij})$ are multivariate normal distributions.

The updates for $\mathbf{f},\mathbf{W},\sigma^{2}_{f},\sigma^{2}_{y}$ are standard VB updates and are available in Infer.NET. The update for the ARD parameters $a_{j}$ however required specific implementation. The factor itself is

where $\stackrel{{\scriptstyle c}}{{=}}$ denotes equality up to an additive constant. Taking expectations with respect to $\mathbf{f}$ under $q$ we obtain the VMP message to $a_{j}$ as being $\text{IG}\left(a_{j};\frac{N}{2}-1,\frac{1}{2}\langle\mathbf{f}_{j}^{T}K_{f}^{-1}\mathbf{f}_{j}\rangle\right)$. Since the variational posterior on $\mathbf{f}$ is multivariate normal the expectation $\langle\mathbf{f}_{j}^{T}K_{f}^{-1}\mathbf{f}_{j}\rangle$ is straightforward to calculate.

M-step.

In the M-step we optimise the variational lower bound with respect to the log length scale parameters $\{\theta_{f},\theta_{w}\}$, using gradient descent with line search. When optimising $\theta_{f}$ we only need to consider the contribution to the lower bound of the factor $\mathcal{N}(\mathbf{f}_{j};0,a_{j}K_{f_{j}})$ (see $\eqref{eqn:gpfactor}$), which is straightforward to evaluate and differentiate (see Appendix). For $\theta_{w}$ we consider the contribution of $\mathcal{N}(\mathbf{W}_{pq};0,K_{W})$.

3 Computational Considerations

GPRN is mainly limited by taking the Cholesky decomposition of the block diagonal $C_{B}$, an $Nq(p+1)\times Nq(p+1)$ matrix. But $pq$ of these blocks are the same $N\times N$ covariance matrix $K_{w}$ for the weight functions, and $q$ of these blocks are the covariance matrices $K_{\hat{f}_{i}}$ associated with the node functions, and $\text{chol}(\text{blkdiag}(A,B,\dots))=\text{blkdiag}(\text{chol}(A),\text{chol}(B),\dots)$. Therefore assuming the node functions share the same covariance function (which they do in our experiments), the complexity of this operation is only $\mathcal{O}(N^{3})$, the same as for regular Gaussian process regression. At worst it is $\mathcal{O}(qN^{3})$, assuming different covariance functions for each node.

Sampling also requires likelihood evaluations. Since there are input dependent noise correlations between the elements of the $p$ dimensional observations $\bm{y}(x_{i})$, multivariate volatility models would normally require inverting a $p\times p$ covariance matrix $N$ times, like MGARCH (Bollerslev et al., , 1988) or multivariate stochastic volatility models (Harvey et al., , 1994). This would lead to a total complexity of $\mathcal{O}(Nqp+Np^{3})$. However, by working directly with the noisy $\hat{\bm{f}}$ instead of the noise free $\bm{f}$, evaluating the likelihood requires no costly or numerically unstable inversions, and thus has a complexity of only $\mathcal{O}(Nqp)$. This allows GPRN to scale to high dimensions; indeed we have a 1000 dimensional gene expression experiment in Section 4.

The computational complexity of VB is dominated by the $\mathcal{O}(N^{3})$ inversions required to calculate the covariance of the node and weight functions in the E-step. Naively $q$ and $qp$ such inversions are required per iteration for the node and weight functions respectively, giving a total complexity of $\mathcal{O}(qpN^{3})$. However, under VB the covariances of the weight functions for the same $p$ are all equal, reducing the complexity to $\mathcal{O}(qN^{3})$. If $p$ is large the $\mathcal{O}(pqN^{2})$ cost of calculating the weight function means may become significant. Although the per iteration cost of VB is actually higher than for MCMC far fewer iterations are typically required to reach convergence.

Overall, even though the GPRN accounts for input dependent signal correlations (rather than fixing the correlations like other multi-task methods), the computational demands of GPRN compare favourably to most multi-task GP models, which commonly have a complexity of $\mathcal{O}(p^{3}N^{3})$ (Alvarez and Lawrence, , 2011).

Experiments

We compare the GPRN to multi-task learning and multivariate volatility models, and we also use the GPRN to gain new scientific insights into the data we model. Furthermore, we compare between variational Bayes (VB) and Markov chain Monte Carlo (MCMC) inference within the GPRN framework. To keep our comparisons up to date, we exactly reproduce many of the experiments in recent papers by Alvarez and Lawrence, (2011) and Wilson and Ghahramani, 2010b on benchmark datasets. In the multi-task setting, there are $p$ dimensional observations $\bm{y}(x)$, and the goal is to use the correlations between the elements of $\bm{y}(x)$ to make better predictions of $\bm{y}(x_{*})$, for a test input $x_{*}$, than if we were to treat the dimensions independently. A major difference between GPRN and alternative multi-task models is that the GPRN accounts for signal correlations that change with $x$, rather than fixed correlations. It also accounts for changing noise correlations (multivariate volatility).

We compare to the following multi-task GP methods: 1) the linear model of coregionalisation (LMC) (Journel and Huijbregts, , 1978; Goovaerts, , 1997), 2) the intrinsic coregionalisation model (ICM) (Goovaerts, , 1997), 3) ordinary co-kriging (Cressie, , 1993; Goovaerts, , 1997; Wackernagel, , 2003), 4) the semiparametric latent factor model (SLFM) (Teh et al., , 2005), 5) convolved multiple output Gaussian processes (CMOGP) (Barry and Jay, , 1996; Ver Hoef and Barry, , 1998; Boyle and Frean, , 2004), 6) standard independent Gaussian processes (GP), 7) and the DTC (Csató and Opper, , 2001; Seeger et al., , 2003; Quiñonero-Candela and Rasmussen, , 2005; Rasmussen and Williams, , 2006) , 8) FITC (Snelson and Ghahramani, , 2006), and 9) PITC (Quiñonero-Candela and Rasmussen, , 2005) sparse approximations for CMOGP (Alvarez and Lawrence, , 2011), which we respectively label as MDTC, MFITC and MPITC. Detail about each of these methods is in Alvarez and Lawrence, (2011). We compare on a 3 dimensional geostatistics heavy metal dataset from the Swiss Jura, where 28% of the observations for one of the outputs (response variables) is missing, and on gene expression datasets with 50 and 1000 dimensional time dependent outputs $\bm{y}(t)$.

In the multi-task experiments, the GPRN accounts for input dependent noise covariance matrices $\text{cov}[\bm{y}(x)]=\Sigma(x)$. To specifically test GPRN’s ability to model input dependent noise covariances (multivariate volatility), we also compare predictions of $\Sigma(x)$ to those made by popular multivariate volatility models – full BEKK MGARCH (Engle and Kroner, , 1995), generalised Wishart processes (Wilson and Ghahramani, 2010b, ), the original Wishart process (Bru, , 1991; Gouriéroux et al., , 2009), and empirical estimates – on benchmark return series datasets which are especially suited to MGARCH (Poon and Granger, , 2005; Hansen and Lunde, , 2005; Brownlees et al., , 2009; McCullough and Renfro, , 1998; Brooks et al., , 2001).

In all experiments, GPRN uses a squared exponential covariance function for its node functions, and another squared exponential covariance function for its weight functions.

Tomancak et al., (2002) measured gene expression levels every hour for 12 hours during Drosophila embryogenesis; they then repeated this experiment for an independent replica (a second independent time series). Gene expression is activated and deactivated by transcription factor proteins. We focus on genes which are thought to at least be regulated by the transcription factor twi, which influences mesoderm and muscle development in Drosophila (Zinzen et al., , 2009). The assumption is that these gene expression levels are all correlated. We would like to use how these correlations change over time to make better predictions of time varying gene expression in the presence of transcription factors. In total there are 1621 genes (outputs) at $N=12$ time points (inputs), on two independent replicas. For training, $p=50$ random genes were selected from the first replica, and the corresponding 50 genes in the second replica were used for testing. We then repeated this experiment 10 times with a different set of genes each time, and averaged the results. We then repeated the whole experiment, but with $p=1000$ genes. We used exactly the same training and testing sets as Alvarez and Lawrence, (2011).

We use a relatively small $p=50$ dataset so that we are able to compare with popular alternative multi-task methods (LMC, CMOGP, SLFM) which have a complexity of $\mathcal{O}(N^{3}p^{3})$ and would not scale to $p=1000$ (Alvarez and Lawrence, , 2011). For $p=1000$, we compare to the sparse convolved multiple output GP methods (MFITC, MDTC, and MPITC) of Alvarez and Lawrence, (2011). In both of these regressions, the GPRN is accounting for multivariate volatility; this is the first time a multivariate stochastic volatility model has been estimated for $p>50$ (Chib et al., , 2006). We assess performance using standardised mean square error (SMSE) and mean standardized log loss (MSLL), as defined in Rasmussen and Williams, (2006) on page 23, and discussed in Alvarez and Lawrence, (2011) on page 1469. Using the empirical mean and variance to fit the data would give an SMSE and MSLL of $1$ and respectively. The smaller the SMSE and more negative the MSLL the better.

The results are in Table 4.3, under the headings GENE (50D) and GENE (1000D). For SET 1 of the 50D dataset, we used Replica 1 and Replica 2 in Alvarez and Lawrence, (2011) respectively as training and testing replicas. We follow Alvarez and Lawrence, (2011) and reverse training and testing replicas to create SET 2. The results for LMC, CMOGP, MFITC, MPITC, and MDTC are reproduced from Alvarez and Lawrence, (2011). GPRN significantly outperforms all of the other models, with between 46% and 68% of the SMSE, and similarly strong results on the MSLL error metric. On the 50D dataset the MCMC and VB results are comparable. However, on the 1000D dataset GPRN with VB noticeably outperforms GPRN with MCMC, likely because MCMC is not mixing as well in high dimensions. Indeed VB may have an advantage over MCMC in high dimensions. On the other hand, GPRN with MCMC is still robust, outperforming all the other methods on the 1000 dimensional dataset.

On both the 50 and 1000 dimensional datasets, the marginal likelihood for the network structure is sharply peaked at $q=1$. This is evidence for the hypothesis that there is only the one transcription factor twi controlling the expression levels of the genes in question. These datasets can be found in Neil Lawrence’s GPSIM toolbox: http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/gpsim/

Typical GPRN (VB) runtimes for the 50D and 1000D datasets were respectively 12 seconds and 330 seconds.

2 Jura Geostatistics

Here we are interested in predicting concentrations of cadmium at 100 locations within a 14.5 km2 region of the Swiss Jura. For training, we have access to measurements of cadmium at 259 neighbouring locations. We also have access to nickel and zinc concentrations at these 259 locations, as well as at the 100 locations we wish to predict cadmium. While a standard Gaussian process regression model would only be able to make use of the cadmium training measurements, a multi-task method can use the correlated nickel and zinc measurements to enhance predictionsThis can be seen as a multivariate missing data problem, with $p=3$ outputs.. With GPRN we can also make use of how the correlations between nickel, zinc, and cadmium change with location to further enhance predictions.

Here the network structure with by far the highest marginal likelihood has $q=2$ latent node functions. The node and weight functions learnt using VB for this setting are shown in Figure 2. Since there are $p=3$ output dimensions, the result $q<p$ suggests that heavy metal concentrations in the Swiss Jura are correlated. Indeed, using our model we can observe the spatially varying correlations between heavy metal concentrations, as shown for cadmium and zinc in Figure 3. Although the correlation between cadmium and zinc is generally positive (with values around 0.6), there is a region where the correlations drop of noticeably, perhaps corresponding to a geological structure. The quantitative results in Table 1 confirm that the ability of GPRN to learn these spatially varying correlations is beneficial in terms of being able to predict cadmium concentrations.

We assess performance quantitatively using mean absolute error (MAE) between the predicted and true cadmium concentrations. We restart the experiment 10 times with different initialisations of the parameters, and average the MAE. The results are marked by JURA in Table 4.3. This experiment follows Goovaerts, (1997) and Alvarez and Lawrence, (2011). The results for SLFM, ICM and CMOGP are from Alvarez and Lawrence, (2011), and the results for co-kriging are from Goovaerts, (1997). It is unclear what preprocessing was performed for these methods, but we found log transforming and normalising each dimension to have zero mean and unit variance to be beneficial due to the skewed distribution of the $y$-values (but we also include results on untransformed data, marked with *). All of the multiple output methods give lower MAE than using an independent GP, and GPRN outperforms SLFM and the other multiple output methods.

For the JURA dataset, the improved performance of GPRN is at the cost of a slightly greater runtime. However, GPRN is accounting for input dependent signal and noise correlations, unlike the other methods. Moreover, the complexity of GPRN scales with $p$ as $\mathcal{O}(Nqp)$, unlike the other methods which scale as $\mathcal{O}(N^{3}p^{3})$ (Alvarez and Lawrence, , 2011). This is why GPRN runs relatively quickly on the 1000 dimensional gene expression dataset, for which the other methods are intractable. This data is available from http://www.ai-geostats.org/.

3 Multivariate Volatility

In the previous experiments the GPRN implicitly accounted for multivariate volatility (input dependent noise covariance) in making predictions of $\bm{y}(x_{*})$. The GPRN incorporates a generalised Wishart process (Wilson and Ghahramani, 2010b, ; Wilson and Ghahramani, , 2011) noise model into a more general model which can also account for signal correlations and other nonstationarities. Although our focus is the ability of GPRN to model input dependent correlations in a multiple output regression setting, we here test the GPRN explicitly as a model of multivariate volatility, and assess predictions of $\Sigma(t)=\text{cov}[\bm{y}(t)]$, where the observations $\bm{y}$ are time dependent. We make historical predictions at observed time points, and one day ahead forecasts. Historical predictions can be used, for example, to understand a past financial crisis. We follow Wilson and Ghahramani, 2010b exactly, and predict $\Sigma(t)$ for returns on three currency exchanges (EXCHANGE) and five equity indices (EQUITY) processed exactly as in Wilson and Ghahramani, 2010b . These datasets are especially suited to MGARCH, the most popular multivariate volatility model, and have become a benchmark for assessing GARCH models (Poon and Granger, , 2005; Hansen and Lunde, , 2005; Brownlees et al., , 2009; McCullough and Renfro, , 1998; Brooks et al., , 2001). We make 200 historical predictions of $\Sigma(x)$ and 200 one step ahead forecasts. The forecasts are assessed using the log likelihood of the new observations under the predicted covariance, denoted $\mathcal{L}$ Forecast. We compare to full BEKK MGARCH (Engle and Kroner, , 1995), the generalised Wishart process (Wilson and Ghahramani, 2010b, ), the original Wishart process (Bru, , 1991; Gouriéroux et al., , 2009), and using the empirical covariance of the training set. We see in Table 4.3 that GPRN (VB) is competitive with MGARCH, even though these datasets are particularly suited to MGARCH. The Historical MSE for EXCHANGE is between the learnt covariance $\Sigma(x)$ and $\mathbf{y}(x)\mathbf{y}(x)^{\top}$, so the high MSE values for GPRN on EXCHANGE are essentially training error, and are less meaningful than the encouraging step ahead forecast likelihoods. The historical predictions are more relevant in EQUITY, where we can compare to the true covariances. See Wilson and Ghahramani, 2010b for details.

GPRN and GWP are both highly flexible but fully Bayesian models for multivariate volatility, so it understandable that their performance is comparable. While GPRN (MCMC) sometimes outperforms MGARCH, the GWP, and the WP, it is often outperformed by GPRN (VB) on the multivariate volatility, perhaps suggesting convergence problems. These data were obtained using Bloomberg (http://www.bloomberg.com/).

A Gaussian process regression network (GPRN) has a simple and interpretable structure, and generalises many of the recent extensions to the Gaussian process regression framework. The model naturally accommodates input dependent signal and noise correlations between multiple output variables, heavy tailed predictive distributions, input dependent length-scales and amplitudes, and adaptive covariance functions. Furthermore, GPRN has scalable inference procedures, and strong empirical performance on several benchmark datasets.

In the future, it would be enlightening to use GPRN with different types of adaptive covariance structures, particularly in the case where $p=1$ and $q>1$; in one dimensional output space it would be easy, for instance, to visualise a process gradually switching between brownian motion, periodic, and smooth covariance functions. It would also be interesting to apply this adaptive network to classification. We hope the GPRN will inspire further research into adaptive networks, and further connections between different areas of machine learning and statistics.

Acknowledgements

Thanks to Mauricio Álvarez and Neil Lawrence for their valuable feedback about the gene expression datasets.

Appendix

Since we extend the Gaussian process framework, we briefly review Gaussian process regression, some notation, and expand on some of the points in the introduction. For more detail see Rasmussen and Williams, (2006).

A Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution. Using a Gaussian process, we can define a distribution over functions $w(x)$:

where $w$ is the output variable, $x$ is an arbitrary (potentially vector valued) input variable, and the mean $m(x)$ and covariance function (or kernel) $k(x,x^{\prime})$ are respectively defined as

This means that any collection of function values has a joint Gaussian distribution:

where the $N\times N$ covariance matrix $K$ has entries $K_{ij}=k(x_{i},x_{j})$, and the mean $\bm{\mu}$ has entries $\bm{\mu}_{i}=m(x_{i})$. The properties of these functions (smoothness, periodicity, etc.) are determined by the kernel function. The squared exponential kernel is popular:

The Ornstein-Uhlenbeck kernel is also widely applied:

In one dimension it is the covariance function of an Ornstein-Uhlenbeck process (Uhlenback and Ornstein, , 1930), which was introduced to model the velocity of a particle undergoing Brownian motion. With this kernel, the corresponding GP is a continuous time AR(1) process with Markovian dynamics: $w(x+a)$ is independent of $w(x-a)$ given $w(x)$ for any constant $a$. Indeed the OU kernel belongs to a more general class of Matérn kernels,

where $K_{\alpha}$ is a modified Bessel function (Abramowitz and Stegun, , 1964). In one dimension the corresponding GP is a continuous time AR($p$) process, where $p=\alpha+1/2$.Discrete time autoregressive processes such as $w(t+1)=w(t)+\epsilon(t)$, where $\epsilon(t)\sim\mathcal{N}(0,1)$, are widely used in time series modelling and are a particularly simple special case of Gaussian processes. The OU kernel is recovered by setting $\alpha=1/2$.

There are many other useful kernels, like the periodic kernel (with a period that can be learned from data), or the Gibbs kernel (Gibbs, , 1997) which allows for input dependent length-scales. Kernels can be combined together, e.g. $k=a_{1}k_{1}+a_{2}k_{2}+a_{3}k_{3}$, and the relative importance of each kernel can be determined from data (e.g. from estimating $a_{1},a_{2},a_{3}$). Rasmussen and Williams, (2006) and Bishop, (2006) have a discussion about how to create and combine kernels.

Suppose we are doing a regression using points $\{y(x_{1}),\dots,y(x_{N})\}$ from a noisy function $y=w(x)+\epsilon$, where $\epsilon$ is additive i.i.d Gaussian noise, such that $\epsilon\sim\mathcal{N}(0,\sigma_{n}^{2})$. Letting $\bm{y}=(y(x_{1}),\dots,y(x_{N}))^{\top}$, and $\bm{w}=(w(x_{1}),\dots,w(x_{N})^{\top}$, we have $p(\bm{y}|\bm{w})=\mathcal{N}(\bm{w},\sigma_{n}^{2}I)$ and $p(\bm{w})=\mathcal{N}(\bm{\mu},K)$ as above. For notational simplicity, we assume $\bm{\mu}=0$. For a test point $w(x_{*})$, the joint distribution $p(w(x_{*}),\bm{y})$ is Gaussian:

where $K$ is defined as above, and $(\bm{k}_{*})_{i}=k(x_{*},x_{i})$ with $i=1,\dots,N$. We can therefore condition on $\bm{y}$ to find $p(w(x_{*})|\bm{y})=\mathcal{N}(\mu_{*},v_{*})$ where

We can find this more laboriously by noting that $p(\bm{w}|\bm{y})$ and $p(w(x_{*})|\bm{w})$ are Gaussian and integrating, since $p(w(x_{*})|\bm{y})=\int p(w(x_{*})|\bm{w})p(\bm{w}|\bm{y})d\bm{w}$.

We see that (23) doesn’t depend on the data $\bm{y}$, just on how far away the test point $x_{*}$ is from the training inputs $\{x_{1},\dots,x_{N}\}$.

In regards to the introduction, we also see that for this standard Gaussian process regression, the observation model $p(y|w)$ is Gaussian, the predictive distribution in (22) and (23) is Gaussian, the marginals in the prior (from marginalising equation (15)) are Gaussian, the noise is constant, and in the popular covariance functions given, the amplitude and length-scale are constant. A brief discussion of multiple outputs, noise models with dependencies, and non-Gaussian observation models can be found in sections 9.1, 9.2 and 9.3 on pages 190-191 of Rasmussen and Williams, (2006), available free online at the book website www.gaussianprocess.org/gpml. An example of an input dependent length-scale is in section 4.2 on page 43.

2 Constraining W𝑊W

It is possible to reduce the number of modes in the posterior by somehow constraining the weights $W$ to be positive. For MCMC it is straightforward to do this by exponentiating the weights, as in Adams and Stegle, (2008) and Adams et al., (2010). For VB it is more straightforward to explicitly constrain the weights to be positive using a truncated Gaussian representation. We found that these extensions did not significantly improve empirical performance, although exponentiating the weights sometimes improved numerical stability for MCMC on the multivariate volatility experiments. For Adams and Stegle, (2008) exponentiating the weights will have been more valuable because they use Expectation Propagation which is known to perform badly in the presence of multimodality. MCMC and VB approaches are more robust to this problem.

3 VB M-step

We will need the gradient with respect to $\theta_{f}$:

The expectations here are straightforward to compute analytically.

4 VB predictive distributions

The predictive distribution is calculated as

VB fits the approximation $p(W(x),f(x)|\mathcal{D})=q(W)q(f)$, so the approximate predictive is

We can calculate the mean and covariance of this distribution analytically:

It is also of interest to calculate the noise covariance. Recall our model can be written as

Let $\bm{n}=\sigma_{f}W(x)\bm{\epsilon}+\sigma_{y}\bm{z}$ be the noise. The covariance of $\mathbf{n}$ is then

References