Thoughts on Massively Scalable Gaussian Processes

Andrew Gordon Wilson, Christoph Dann, Hannes Nickisch

Introduction

Every minute of the day, users share hundreds of thousands of pictures, videos, tweets, reviews, and blog posts. More than ever before, we have access to massive datasets in almost every area of science and engineering, including genomics, robotics, and climate science. This wealth of information provides an unprecedented opportunity to automatically learn rich representations of data, which allows us to greatly improve performance in predictive tasks, but also provides a mechanism for scientific discovery.

Expressive non-parametric methods, such as Gaussian processes (GPs) (Rasmussen and Williams,, 2006), have great potential for large-scale structure discovery; indeed, these methods can be highly flexible, and have an information capacity that grows with the amount of available data. However, large data problems are mostly uncharted territory for GPs, which can only be applied to at most a few thousand training points $n$ , due to the $\mathcal{O}(n^{3})$ computations and $\mathcal{O}(n^{2})$ storage required for inference and learning.

Even more scalable approximate GP approaches, such as inducing point methods (Quiñonero-Candela and Rasmussen, 2005a, ), typically require $\mathcal{O}(m^{2}n+m^{3})$ computations and $\mathcal{O}(m^{2}+mn)$ storage, for $m$ inducing points, and are hard to apply massive datasets, containing $n>10^{5}$ examples. Moreover, for computational tractability, these approaches require $m\ll n$ , which can severely affect predictive performance, limit representational power, and the ability for kernel learning, which is most needed on large datasets (Wilson,, 2014). New directions for scalable Gaussian processes have involved mini-batches of data through stochastic variational inference (Hensman et al.,, 2013) and distributed learning (Deisenroth and Ng,, 2015). While these approaches are promising, inference can undergo severe approximations, and a small number of inducing points are still required. Indeed, stochastic variational approaches scale as $\mathcal{O}(m^{3})$ .

In this paper, we introduce a new framework for massively scalable Gaussian processes (MSGP), which provides near-exact $\mathcal{O}(n)$ inference and learning and $\mathcal{O}(1)$ test time predictions, and does not require distributed learning or severe assumptions. Our approach builds on the recently introduced KISS-GP framework (Wilson and Nickisch,, 2015), with several significant advances which enable its use on massive datasets. In particular, we provide:

Near-exact $\mathcal{O}(1)$ mean and variance predictions. By contrast, standard GPs and KISS-GP cost $\mathcal{O}(n)$ for the predictive mean and $\mathcal{O}(n^{2})$ for the predictive variance per test point. Moreover, inducing point and finite basis expansions (e.g., Quiñonero-Candela and Rasmussen, 2005a, ; Lázaro-Gredilla et al.,, 2010; Yang et al.,, 2015) cost $\mathcal{O}(m)$ and $\mathcal{O}(m^{2})$ per test point.

Circulant approximations which (i) integrate Kronecker and Toeplitz structure, (ii) enable extremely and accurate fast log determinant evaluations for kernel learning, and (iii) increase the speed of Toeplitz methods on problems with 1D predictors.

The ability to exploit more general block-Toeplitz-Toeplitz-block (BTTB) structure, which enables fast and exact inference and learning in cases where multidimensional Kronecker structure is not present.

Projections which help alleviate the limitation of Kronecker methods to low-dimensional input spaces.

Code will be available as part of the GPML package (Rasmussen and Nickisch,, 2010), with demonstrations at http://www.cs.cmu.edu/~andrewgw/pattern.

We begin by briefly reviewing Gaussian processes, structure exploiting inference, and KISS-GP, in sections 2 - 4. We then introduce our MSGP approach in section 5. We demonstrate the scalability and accuracy of MSGP in the experiments of section 6. We conclude in section 7.

Gaussian Processes

We briefly review Gaussian processes (GPs), and the computational requirements for predictions and kernel learning. Rasmussen and Williams, (2006) contains a full treatment of GPs.

We assume a dataset $\mathcal{D}$ of $n$ input (predictor) vectors $X=[{\bm{\mathbf{x}}}_{1},\dots,{\bm{\mathbf{x}}}_{n}]$ , each of dimension $D$ , corresponding to a $n\times 1$ vector of targets ${\bm{\mathbf{y}}}=[y({\bm{\mathbf{x}}}_{1}),\dots,y({\bm{\mathbf{x}}}_{n})]^{\top}$ . If $f({\bm{\mathbf{x}}})\sim\mathcal{GP}(\mu,k_{{\bm{\mathbf{\theta}}}})$ , then any collection of function values ${\bm{\mathbf{f}}}$ has a joint Gaussian distribution,

with mean vector and covariance matrix defined by the mean vector and covariance function of the Gaussian process: $({\bm{\mathbf{\mu}}}_{X})_{i}=\mu(x_{i})$ , and $(K_{X,X})_{ij}=k_{\bm{\mathbf{\theta}}}({\bm{\mathbf{x}}}_{i},{\bm{\mathbf{x}}}_{j})$ . The covariance function $k_{\bm{\mathbf{\theta}}}$ is parametrized by ${\bm{\mathbf{\theta}}}$ . Assuming additive Gaussian noise, $y({\bm{\mathbf{x}}})|f({\bm{\mathbf{x}}})\sim\mathcal{N}(y({\bm{\mathbf{x}}});f({\bm{\mathbf{x}}}),\sigma^{2})$ , then the predictive distribution of the GP evaluated at the $n_{*}$ test points indexed by $X_{*}$ , is given by

$K_{X_{*},X}$ represents the $n_{*}\times n$ matrix of covariances between the GP evaluated at $X_{*}$ and $X$ , and all other covariance matrices follow the same notational conventions. ${\bm{\mathbf{\mu}}}_{X_{*}}$ is the $n_{*}\times 1$ mean vector, and $K_{X,X}$ is the $n\times n$ covariance matrix evaluated at training inputs $X$ . All covariance matrices implicitly depend on the kernel hyperparameters ${\bm{\mathbf{\theta}}}$ .

The marginal likelihood of the targets ${\bm{\mathbf{y}}}$ is given by

where we have used $K_{{\bm{\mathbf{\theta}}}}$ as shorthand for $K_{X,X}$ given ${\bm{\mathbf{\theta}}}$ . Kernel learning is performed by optimizing Eq. (3) with respect to ${\bm{\mathbf{\theta}}}$ .

The computational bottleneck for inference is solving the linear system $(K_{X,X}+\sigma^{2}I)^{-1}{\bm{\mathbf{y}}}$ , and for kernel learning is computing the log determinant $\log|K_{X,X}+\sigma^{2}I|$ . Standard procedure is to compute the Cholesky decomposition of the $n\times n$ matrix $K_{X,X}$ , which requires $\mathcal{O}(n^{3})$ operations and $\mathcal{O}(n^{2})$ storage. Afterwards, the predictive mean and variance of the GP cost respectively $\mathcal{O}(n)$ and $\mathcal{O}(n^{2})$ per test point ${\bm{\mathbf{x}}}_{*}$ .

Structure Exploiting Inference

Structure exploiting approaches make use of existing structure in $K_{X,X}$ to accelerate inference and learning. These approaches benefit from often exact predictive accuracy, and impressive scalability, but are inapplicable to most problems due to severe grid restrictions on the data inputs $X$ . We briefly review Kronecker and Toeplitz structure.

Kronecker (tensor product) structure arises when we have multidimensional inputs i.e. $P>1$ on a rectilinear grid, ${\bm{\mathbf{x}}}\in\mathcal{X}_{1}\times\dots\times\mathcal{X}_{P}$ , and a product kernel across dimensions $k({\bm{\mathbf{x}}}_{i},{\bm{\mathbf{x}}}_{j})=\prod_{p=1}^{P}k({\bm{\mathbf{x}}}_{i}^{(p)},{\bm{\mathbf{x}}}_{j}^{(p)})$ . In this case, $K=K_{1}\otimes\dots\otimes K_{P}$ . One can then compute the eigendecomposition of $K=QVQ^{\top}$ by separately taking the eigendecompositions of the much smaller $K_{1},\dots,K_{P}$ . Inference and learning then proceed via $(K+\sigma^{2}I)^{-1}{\bm{\mathbf{y}}}=(QVQ^{\top}+\sigma^{2}I)^{-1}{\bm{\mathbf{y}}}=Q(V+\sigma^{2}I)^{-1}Q^{\top}{\bm{\mathbf{y}}}$ , and $\log|K+\sigma^{2}I|=\sum_{i}\log(V_{ii}+\sigma^{2})$ , where $Q$ is an orthogonal matrix of eigenvectors, which also decomposes a Kronecker product (allowing for fast MVMs), and $V$ is a diagonal matrix of eigenvalues and thus simple to invert. Overall, for $m$ grid data points, and $P$ grid dimensions, inference and learning cost $\mathcal{O}(Pm^{1+\frac{1}{P}})$ operations (for $P>1$ ) and $\mathcal{O}(Pm^{\frac{2}{P}})$ storage (Saatchi,, 2011; Wilson et al.,, 2014). Unfortunately, there is no efficiency gain for 1D inputs (e.g., time series).

2 Toeplitz Structure

A covariance matrix constructed from a stationary kernel $k(x,z)=k(x-z)$ on a 1D regularly spaced grid has Toeplitz structure. Toeplitz matrices $T$ have constant diagonals, $T_{i,j}=T_{i+1,j+1}$ . Toeplitz structure has been exploited for GP inference (e.g., Zhang et al.,, 2005; Cunningham et al.,, 2008) in $\mathcal{O}(n\log n)$ computations. Computing $\log|T|$ , and the predictive variance for a single test point, requires $\mathcal{O}(n^{2})$ operations (although finite support can be exploited (Storkey,, 1999)) thus limiting Toeplitz methods to about $n<10,000$ points when kernel learning is required. Since Toeplitz methods are limited to problems with 1D inputs (e.g., time series), they complement Kronecker methods, which exploit multidimensional grid structure.

KISS-GP

Recently, Wilson and Nickisch, (2015) introduced a fast Gaussian process method called KISS-GP, which performs local kernel interpolation, in combination with inducing point approximations (Quiñonero-Candela and Rasmussen, 2005b, ) and structure exploiting algebra (e.g., Saatchi,, 2011; Wilson,, 2014).

for any single inputs ${\bm{\mathbf{x}}}$ and ${\bm{\mathbf{z}}}$ . The $n\times n$ training covariance matrix $K_{X,X}$ thus has the approximation

Wilson and Nickisch, (2015) propose to perform local kernel interpolation, in a method called KISS-GP, for extremely sparse interpolation matrices. For example, if we are performing local cubic (Keys,, 1981) interpolation for $d$ -dimensional input data, $W_{X}$ and $W_{Z}$ contain only $4^{d}$ non-zero entries per row.

Furthermore, Wilson and Nickisch, (2015) show that classical inducing point methods can be re-derived within their SKI framework as global interpolation with a noise free GP and non-sparse interpolation weights. For example, the subset of regression (SoR) inducing point method effectively uses the kernel ${k}_{\text{SoR}}({\bm{\mathbf{x}}},{\bm{\mathbf{z}}})=K_{{\bm{\mathbf{x}}},U}K_{U,U}^{-1}K_{U,{\bm{\mathbf{z}}}}$ (Quiñonero-Candela and Rasmussen, 2005b, ), and thus has interpolation weights ${\bm{\mathbf{w}}}_{\text{SoR}({\bm{\mathbf{x}}})}=K_{{\bm{\mathbf{x}}},U}K_{U,U}^{-1}$ within the SKI framework.

GP inference and learning can be performed in $\mathcal{O}(n)$ using KISS-GP, a significant advance over the more standard $\mathcal{O}(m^{2}n)$ scaling of fast GP methods (Quiñonero-Candela and Rasmussen, 2005a, ; Lázaro-Gredilla et al.,, 2010). Moreover, Wilson and Nickisch, (2015) show how – when performing local kernel interpolation – one can achieve close to linear scaling with the number of inducing points $m$ by placing these points $U$ on a rectilinear grid, and then exploiting Toeplitz or Kronecker structure in $K_{U,U}$ (see, e.g., Wilson,, 2014), without requiring that the data inputs $X$ are on a grid. Such scaling with $m$ compares favourably to the $\mathcal{O}(m^{3})$ operations for stochastic variational approaches (Hensman et al.,, 2013). Allowing for large $m$ enables near-exact performance, and large scale kernel learning.

In particular, for inference we can solve $(K_{\text{SKI}}+\sigma^{2}I)^{-1}{\bm{\mathbf{y}}}$ , by performing linear conjugate gradients, an iterative procedure which depends only on matrix vector multiplications (MVMs) with $(K_{\text{SKI}}+\sigma^{2}I)$ . Only $j\ll n$ iterations are required for convergence up to machine precision, and the value of $j$ in practice depends on the conditioning of $K_{\text{SKI}}$ rather than $n$ . MVMs with sparse $W$ (corresponding to local interpolation) cost $\mathcal{O}(n)$ , and MVMs exploiting structure in $K_{U,U}$ are roughly linear in $m$ . Moreover, we can efficiently approximate the eigenvalues of $K_{\text{SKI}}$ to evaluate $\log|K_{\text{SKI}}+\sigma^{2}I|$ , for kernel learning, by using fast structure exploiting eigendecompositions of $K_{U,U}$ . Further details are in Wilson and Nickisch, (2015).

Massively Scalable Gaussian Processes

We introduce massively scalable Gaussian processes (MSGP), which significantly extend KISS-GP, inducing point, and structure exploiting approaches, for: (1) $\mathcal{O}(1)$ test predictions (section 5.1); (2) circulant log determinant approximations which (i) unify Toeplitz and Kronecker structure; (ii) enable extremely fast marginal likelihood evaluations (section 5.2); and (iii) extend KISS-GP and Toeplitz methods for scalable kernel learning in D=1 input dimensions, where one cannot exploit multidimensional Kronecker structure for scalability. (3) more general BTTB structure, which enables fast exact multidimensional inference without requiring Kronecker (tensor) decompositions (section 5.3); and, (4) projections which enable KISS-GP to be used with structure exploiting approaches for $D\gg 5$ input dimensions, and increase the expressive power of covariance functions (section 5.4).

We note that the methodology here for fast test predictions does not rely on having performed inference and learning in any particular way: it can be applied to any trained Gaussian process model, including a full GP, inducing points methods such as FITC (Snelson and Ghahramani,, 2006), or the Big Data GP (Hensman et al.,, 2013), or finite basis models (e.g., Yang et al.,, 2015; Lázaro-Gredilla et al.,, 2010; Rahimi and Recht,, 2007; Le et al.,, 2013; Williams and Seeger,, 2001).

1.2 Predictive Variance

Practical applications typically do not require the full predictive covariance matrix $\text{cov}({\bm{\mathbf{f}}}_{*})$ of Eq. (2) for a set of $n_{*}$ test inputs $X_{*}$ , but rather focus on the predictive variance,

where ${\bm{\mathbf{\nu}}}_{*}=\text{diag}(K_{X_{*},X}[K_{X,X}+\sigma^{2}I]^{-1}K_{X,X_{*}})$ , the explained variance, is approximated by local interpolation from the explained variance on the grid $U$

The overall (unbiased) estimate (Papandreou and Yuille,, 2011) is obtained by clipping

2 Circulant Approximation

Kronecker and Toeplitz methods (section 3) are greatly restricted by requiring that the data inputs $X$ are located on a grid. We lift this restriction by creating structure in $K_{U,U}$ , with the unobserved inducing variables $U$ , as part of the structured kernel interpolation framework described in section 4.

Toeplitz methods apply only to 1D problems, and Kronecker methods require multidimensional structure for efficiency gains. Here we present a circulant approximation to unify the complementary benefits of Kronecker and Toeplitz methods, and to greatly scale marginal likelihood evaluations of Toeplitz based methods, while not requiring any grid structure in the data inputs $X$ .

If $U$ is a regularly spaced multidimensional grid, and we use a stationary product kernel (e.g., the RBF kernel), then $K_{U,U}$ decomposes as a Kronecker product of Toeplitz (section 3.2) matrices:

Because the algorithms which leverage Kronecker structure in Gaussian processes require eigendecompositions of the constituent matrices, computations involving the structure in Eq. (11) are no faster than if the Kronecker product were over arbitrary positive definite matrices: this nested Toeplitz structure, which is often present in Kronecker decompositions, is wasted. Indeed, while it is possible to efficiently solve linear systems with Toeplitz matrices, there is no particularly efficient way to obtain a full eigendecomposition.

Fast operations with an $m\times m$ Toeplitz matrix $T$ can be obtained through its relationship with an $a\times a$ circulant matrix. Symmetric circulant matrices $C$ are Toeplitz matrices where the first column ${\bm{\mathbf{c}}}$ is given by a circulant vector: ${\bm{\mathbf{c}}}=[c_{1},c_{2},c_{3},..,c_{3},c_{2}]^{\top}$ . Each subsequent column is shifted one position from the next. In other words, $C_{i,j}={\bm{\mathbf{c}}}_{|j-i|\text{ mod }a}$ . Circulant matrices are computationally attractive because their eigendecomposition is given by

where $F$ is the discrete Fourier transform (DFT): $F_{jk}=\exp(-2jk\pi i/a)$ . The eigenvalues of $C$ are thus given by the DFT of its first column, and the eigenvectors are proportional to the DFT itself (the $a$ roots of unity). The log determinant of $C$ – the sum of log eigenvalues – can therefore be computed from a single fast Fourier transform (FFT) which costs $\mathcal{O}(a\log a)$ operations and $\mathcal{O}(a)$ memory. Fast matrix vector products can be computed at the same asymptotic cost through Eq. (12), requiring two FFTs (one FFT if we pre-compute $F{\bm{\mathbf{c}}}$ ), one inverse FFT, and one inner product.

In a Gaussian process context, fast MVMs with Toeplitz matrices are typically achieved through embedding a $m\times m$ Toeplitz matrix $K$ into a $(2m-1)\times(2m-1)$ circulant matrix $C$ (e.g., Zhang et al.,, 2005; Cunningham et al.,, 2008), with first column ${\bm{\mathbf{c}}}=[k_{1},k_{2},..,k_{m-1},k_{m},k_{m-1},..,k_{2}]$ . Therefore $K=C_{i=1\dots m,j=1\dots m}$ , and using zero padding and truncation, $K{\bm{\mathbf{y}}}=[C{\bm{\mathbf{y}}},{\bm{\mathbf{0}}}]^{\top}_{i=1\dots m}$ , where the circulant MVM $C{\bm{\mathbf{y}}}$ can be computed efficiently through FFTs. GP inference can then be achieved through LCG, solving $K^{-1}{\bm{\mathbf{y}}}$ in an iterative procedure which only involves MVMs, and has an asymptotic cost of $\mathcal{O}(m\log m)$ computations. The log determinant and a single predictive variance, however, require $\mathcal{O}(m^{2})$ computations.

To speed up LCG for solving linear Toeplitz systems, one can use circulant pre-conditioners which act as approximate inverses. One wishes to minimise the distance between the preconditioner $C$ and the Toeplitz matrix $K$ , $\text{arg min}_{C}d(C,K)$ . Three classical pre-conditioners include $C_{\text{Strang}}=\arg\min_{C}\left\|C-K\right\|_{1}$ (Strang,, 1986), $C_{\text{T. Chan}}=\arg\min_{C}\left\|C-K\right\|_{F}$ (Chan,, 1988), and $C_{\text{Tyrtyshnikov}}=\arg\min_{C}\left\|I-C^{-1}K\right\|_{F}$ (Tyrtyshnikov,, 1992).

A distinct line of research was explored more than 60 years ago in the context of statistical inference over spatial processes (Whittle,, 1954). The circulant Whittle approximation $\text{circ}_{\text{Whittle}}({\bm{\mathbf{k}}})$ is given by truncating the sum

i.e. we retain $\sum_{j=-w}^{w}k_{i+jm}$ only.

Positive definiteness of $C=\text{toep}({\bm{\mathbf{c}}})$ for the complete sum is guaranteed by construction (Guinness and Fuentes,, 2014, Section 2). For large lattices, the approach is often used due to its accuracy and favourable asymptotic properties such as consistency, efficiency and asymptotic normality (Lieberman et al.,, 2009). In fact, the circulant approximation $c_{i}$ is asymptotically equivalent to the initial covariance $k_{i}$ , see Gray, (2005, Lemma 4.5), hence the logdet approximation inevitably converges to the exact value.

where we threshold ${\bm{\mathbf{c}}}=F^{\textsf{H}}\max(F\text{circ}({\bm{\mathbf{k}}}),{\bm{\mathbf{0}}})$ . ${\bm{\mathbf{1}}}$ denotes a vector ${\bm{\mathbf{1}}}=[1,1,\dots,1]^{\top}$ , $F^{H}$ is the conjugate transpose of the Fourier transform matrix. $\text{circ}({\bm{\mathbf{k}}})$ denotes one of the T. Chan, Tyrtyshnikov, Strang, Helgason or Whittle approximations.

3 BCCB Approximation for BTTB

We note that blocks in the BTTB matrix need not be symmetric. Moreover, symmetric BCCB matrices – in contrast to symmetric BTTB matrices – are fully characterised by their first column.

4 Projections

Even if we do not exploit structure in $K_{U,U}$ , our framework in section 5 provides efficiency gains over conventional inducing point approaches, particularly for test time predictions. However, if we are to place $U$ onto a multidimensional (Cartesian product) grid so that $K_{U,U}$ has Kronecker structure, then the total number of inducing points $m$ (the cardinality of $U$ ) grows exponentially with the number of grid dimensions, limiting one to about five or fewer grid dimensions for practical applications. However, we need not limit the applicability of our approach to data inputs $X$ with $D\leq 5$ input dimensions, even if we wish to exploit Kronecker structure in $K_{U,U}$ . Indeed, many inducing point approaches suffer from the curse of dimensionality, and input projections have provided an effective countermeasure (e.g., Snelson,, 2007).

The resulting kernel generalises the RBF and ARD kernels, which respectively have spherical and diagonal covariance matrices, with a full covariance matrix $\Sigma=PP^{\top}$ , which allows for richer models (Vivarelli and Williams,, 1998). But in our context, this added flexibility has special meaning. Kronecker methods typically require a kernel which separates as a product across input dimensions (section 3.1). Here, we can capture sophisticated correlations between the different dimensions of the data inputs ${\bm{\mathbf{x}}}$ through the projection matrix $P$ , while preserving Kronecker structure in $K_{U,U}$ . Moreover, learning $P$ in a supervised manner, e.g., through the Gaussian porcess marginal likelihood, has immediate advantages over unsupervised dimensionality reduction of ${\bm{\mathbf{x}}}$ ; for example, if only a subset of the data inputs were used in producing the target values, this structure would not be detected by an unsupervised method such as PCA, but can be learned through $P$ . Thus in addition to the critical benefit of allowing for applications with $D>5$ dimensional data inputs, $P$ can also enrich the expressive power of our MSGP model.

The entries of $P$ become hyperparameters of the Gaussian process marginal likelihood of Eq. (3), and can be treated in exactly the same way as standard kernel hyperparameters such as length-scale. One can learn these parameters through marginal likelihood optimisation.

Computing the derivatives of the log marginal likelihood with respect to the projection matrix requires some care under the structured kernel interpolation approximation to the covariance matrix $K_{X,X}$ . We provide the mathematical details in the appendix A.

For practical reasons, one may wish to restrict $P$ to be orthonormal or have unit scaling. We discuss this further in section 6.2.

Experiments

In these preliminary experiments, we stress test MSGP in terms of training and prediction runtime, as well as accuracy, empirically verifying its scalability and predictive performance. We also demonstrate the consistency of the model in being able to learn supervised projections, for higher dimensional input spaces.

We compare to exact Gaussian processes, FITC (Snelson and Ghahramani,, 2006), Sparse Spectrum Gaussian Processes (SSGP) (Lázaro-Gredilla et al.,, 2010), the Big Data GP (BDGP) (Hensman et al.,, 2013), MSGP with Toeplitz (rather than circulant) methods, and MSGP not using the new scalable approach for test predictions (scalable test predictions are described in section 5).

The experiments were executed on a workstation with Intel i7-4770 CPU and 32 GB RAM. We used a step length of $0.01$ for BDGP based on the values reported by Hensman et al., (2013) and a batchsize of $300$ . We also evaluated BDGP with a larger batchsize of $5000$ but found that the results are qualititively similar. We stopped the stochastic optimization in BDGP when the log-likelihood did not improve at least by $0.1$ within the last $50$ steps or after $5000$ iterations.

One cannot exploit Kronecker structure in one dimensional inputs for scalability, and Toeplitz methods, which apply to 1D inputs, are traditionally limited to about $10,000$ points if one requires many marginal likelihood evaluations for kernel learning. Thus to stress test the value of the circulant approximation most transparently, and to give the greatest advantage to alternative approaches in terms of scalability with number of inducing points $m$ , we initially stress test in a 1D input space, before moving onto higher dimensional problems.

In particular, we sample 1D inputs $x$ uniform randomly in $ $, so that the data inputs have no grid structure. We then generate out of class ground truth data$ f(x)=\sin(x)\exp\left(\frac{-x^{2}}{2\times 5^{2}}\right) $with additive Gaussian noise to form$ {\bm{\mathbf{y}}} $. We distribute inducing points on a regularly spaced grid in$ $.

In Figure 2 we show the runtime for a marginal likelihood evaluation as a function of training points $n$ , and number inducing points $m$ . MSGP quickly overtakes the alternatives as $n$ increases past this point. Moreover, the runtimes for MSGP with different numbers of inducing points converge quickly with increases in $m$ . By $n=10^{7}$ , MSGP requires the same training time for $m=10^{3}$ inducing points as it does for $m=10^{6}$ inducing points! Indeed, MSGP is able to accommodate an unprecedented number of inducing points, overtaking the alternatives, which are using $m=10^{3}$ inducing points, when using $m=10^{6}$ inducing points. Such an exceptionally large number of inducing points allows for near-exact accuracy, and the ability to retain model structure necessary for large scale kernel learning. Note that $m=10^{3}$ inducing points (resp. basis functions) is a practical upper limit in alternative approaches (Hensman et al., (2013), for example, gives $m\in$ as a practical upper bound for conventional inducing approaches for large $n$ ).

We emphasize that although the stochastic optimization of the BDGP can take some time to converge, the ability to use mini-batches for jointly optimizing the variational parameters and the GP hyper parameters, as part of a principled framework, makes BDGP highly scalable. Indeed, the methodology in MSGP and BDGP are complementary and could be combined for a particularly scalable Gaussian process framework.

In Figure 3, we show that the prediction runtime for MSGP is practically independent of both $m$ and $n$ , and for any fixed $m,n$ is much faster than the alternatives, which depend at least quadratically on $m$ and $n$ . We also see that the local interpolation strategy for test predictions, introduced in section 5.1, greatly improves upon MSGP using the exact predictive distributions while exploiting Kronecker and Toeplitz algebra.

As we increase the number of inducing points $m$ the quality of predictions are improved, as we expect. The fast predictive variances, based on local kernel interpolation, are not as sensitive to the number of inducing points as the alternatives, but nonetheless have reasonable performance. We also see that the the fast predictive variances are improved by having an increasing number of samples $n_{s}$ in the stochastic estimator of section 5.1.2, which is most noticeable for larger values numbers of inducing points $m$ . Notably, the fast predictive mean, based on local interpolation, is essentially indistinguishable from the predictive mean (‘slow’) using the standard GP predictive equations without interpolation. Overall, the error of the fast predictions is much less than the average variability in the data.

2 Projections

Here we test the consistency of our approach in section 5.4 for recovering ground truth projections, and providing accurate predictions, on $D\gg 5$ dimensional input spaces.

To generate data, we begin by sampling the entries of a $d\times D$ projection matrix $P$ from a standard Gaussian distribution, and then project $n=3000$ inputs ${\bm{\mathbf{x}}}$ (of dimensionality $D\times 1$ ), with locations randomly sampled from a Gaussian distribution so that there is no input grid structure, into a $d=2$ dimensional space: ${\bm{\mathbf{x}}}^{\prime}=P{\bm{\mathbf{x}}}$ . We then sample data ${\bm{\mathbf{y}}}$ from a Gaussian process with an RBF kernel operating on the low dimensional inputs ${\bm{\mathbf{x}}}^{\prime}$ . We repeat this process $30$ times for each of $D=1,2,3,4,5,6,7,8,9,10,15,20,25,\dots,100$ .

We now index the data ${\bm{\mathbf{y}}}$ by the high dimensional inputs ${\bm{\mathbf{x}}}$ and attempt to reconstruct the true low dimensional subspace described by ${\bm{\mathbf{x}}}^{\prime}=P{\bm{\mathbf{x}}}$ . We learn the entries of $P$ jointly with covariance parameters through marginal likelihood optimisation (Eq. (3)). Using a $(d=2)$ $50\times 50$ Cartesian product grid for $2500$ total inducing points $U$ , we reconstruct the projection matrix $P$ , with the subspace error,

shown in Figure 5. Here $G_{i}$ is the orthogonal projection onto the $d$ -dimensional subspace spanned by the rows of $P_{i}$ . This metric is motivated by the one-to-one correspondence between subspaces and orthogonal projections. It is bounded in $ $, where the maximum distance$ 1$ indicates that the subspaces are orthogonal to each other and the minimum is achieved if the subspaces are identical. More information on this metric including how to compute dist of Eq. (13) is available in chapter 2.5.3 of Golub and Van Loan, (2012).

We implement variants of our approach where $P$ is constrained to be 1) orthonormal, 2) to have unit scaling, e.g., $P_{\text{unit}}=\text{diag}(\sqrt{(PP^{\top})}^{-1})P$ . Such constraints prevent degeneracies between $P$ and kernel hyperparameters from causing practical issues, such as length-scales growing to large values to shrink the marginal likelihood log determinant complexity penalty, and then re-scaling $P$ to leave the marginal likelihood model fit term unaffected. In practice we found that unit scaling was sufficient to avoid such issues, and thus preferable to orthonormal $P$ , which are more constrained.

Discussion

We demonstrate these advantages, comparing to state of the art alternatives. In particular, we show MSGP is exceptionally scalable in terms of training points $n$ , inducing points $m$ , and number of testing points. The ability to handle large $m$ will prove important for retaining accuracy in scalable Gaussian process methods, and in enabling large scale kernel learning.

This document serves to report substantial developments regarding the SKI and KISS-GP frameworks introduced in Wilson and Nickisch, (2015).

References

Appendix A APPENDIX

apply chain rule to obtain from $\frac{\partial\psi}{\partial Q}$ from $\frac{\partial\psi}{\partial P}$

A.2 Derivatives For Orthonormal Projections