Preconditioning Kernel Matrices

Introduction

Kernel machines, in enabling flexible feature space representations of data, comprise a broad and important class of tools throughout machine learning and statistics; prominent examples include support vector machines (Schölkopf & Smola, 2001) and Gaussian processes (GPs) (Rasmussen & Williams, 2006). At the core of most kernel machines is the need to solve linear systems involving the Gram matrix $K=\left\{k(\mathbf{x}_{i},\mathbf{x}_{j}\mid\boldsymbol{\theta})\right\}_{i,j=1,...,n}$, where the kernel function $k$, parameterized by $\boldsymbol{\theta}$, implicitly specifies the feature space representation of data points $\mathbf{x}_{i}$. Because $K$ grows with the number of data points $n$, a fundamental computational bottleneck exists: storing $K$ is $\mathcal{O}(n^{2})$, and solving a linear system with $K$ is $\mathcal{O}(n^{3})$. As the need for large-scale kernel machines grows, much work has been directed towards this scaling issue.

Standard approaches to kernel machines involve a factorization (typically Cholesky) of $K$, which is efficient and exact but maintains the quadratic storage and cubic runtime costs. This cost is particularly acute when adapting (or learning) hyperparameters $\boldsymbol{\theta}$ of the kernel function, as $K$ must then be factorized afresh for each $\boldsymbol{\theta}$. To alleviate this burden, numerous works have turned to approximate methods (Candela & Rasmussen, 2005; Snelson & Ghahramani, 2007; Rahimi & Recht, 2008) or methods that exploit structure in the kernel (Gilboa et al., 2015). Approximate methods can achieve attractive scaling, often through the use of low-rank approximations to $K$, but they can incur a potentially severe loss of accuracy. An alternative to factorization is found in the conjugate gradient method (CG), which is used to directly solve linear systems via a sequence of matrix-vector products. Any kernel structure can then be exploited to enable fast multiplications, driving similarly attractive runtime improvements, and eliminating the storage burden (neither $K$ nor its factors need be represented in memory). Unfortunately, in the absence of special structure that accelerates multiplications, CG performs no better than $\mathcal{O}(n^{3})$ in the worst case, and in practice finite numerical precision often results in a degradation of runtime performance compared to factorization approaches.

Throughout optimization, the typical approach to the slow convergence of CG is to apply preconditioners to improve the geometry of the linear system being solved (Golub & Van Loan, 1996). While preconditioning has been explored in domains such as spatial statistics (Chen, 2005; Stein et al., 2012; Ashby & Falgout, 1996), the application of preconditioning to kernel matrices in machine learning has received little attention. Here we design and study preconditioned conjugate gradient methods (PCG) for use in kernel machines, and provide a full exploration of the use of approximations of $K$ as preconditioners.

Our contributions are as follows. (i) Extending the work in (Davies, 2014), we apply a broad range of kernel matrix approximations as preconditioners. Interestingly, this step allows us to exploit the important developments of approximate kernel machines to accelerate the exact computation that PCG offers. (ii) As a motivating example used throughout the paper, we analyze and provide a general framework to both learn kernel parameters and make predictions in GPs. (iii) We extend stochastic gradient learning for GPs (Filippone & Engler, 2015; Anitescu et al., 2012) to allow any likelihood that factorizes over the data points by developing an unbiased estimate of the gradient of the approximate log-marginal likelihood. We demonstrate this contribution in making the first use of PCG for GP classification. (iv) We evaluate datasets over a range of problem size and dimensionality. Because PCG is exact in the limit of iterations (unlike approximate techniques), we demonstrate a tradeoff between accuracy and computational effort that improves beyond state-of-the-art approximation and factorization approaches.

In all, we show that PCG, with a thoughtful choice of preconditioner, is a competitive strategy which is possibly even superior than existing approximation and CG-based techniques for solving general kernel machinesCode to replicate all results in this paper is available at http://github.com/mauriziofilippone/preconditioned_GPs.

Motivating example – Gaussian Processes

Gaussian processes (GPs) are the fundamental building block of many probabilistic kernel machines that can be applied in a large variety of modeling scenarios (Rasmussen & Williams, 2006). Throughout the paper, we will denote by $X=\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\}$ a set of $n$ input vectors and use $\mathbf{y}=(y_{1},\ldots,y_{n})^{\top}$ for the corresponding labels. GPs are formally defined as collections of random variables characterized by the property that any finite number of them is jointly Gaussian distributed. The specification of a kernel function determines the covariance structure of such random variables

In this work we focus in particular on the popular Radial Basis Function (RBF) kernel

where $\boldsymbol{\theta}$ represents the collection of the kernel parameters $\sigma^{2}$ and $l_{r}^{2}$. Defining $f_{i}=f(\mathbf{x}_{i})$ and $\mathbf{f}=(f_{1},\ldots,f_{n})^{\top}$, and assuming a zero mean GP, we have

where $K$ is the $n\times n$ Gram matrix with elements $K_{ij}=k(\mathbf{x}_{i},\mathbf{x}_{j}\mid\boldsymbol{\theta})$. Note that the kernel above and many popular kernels in machine learning give rise to dense kernel matrices. Observations are then modeled through a transformation $h$ of a set of GP-distributed latent variables, specifying the model

The success of nonparametric models based on kernels hinges on the adaptation of kernel parameters $\boldsymbol{\theta}$. The motivation for preconditioning begins with an inspection of the log-marginal likelihood of GP models with prior $\mathcal{N}(\mathbf{f}\mid\mathbf{0},K)$. In Gaussian processes with a Gaussian likelihood $y_{i}\sim\mathcal{N}(y_{i}\mid f_{i},\lambda)$, we have analytic forms for

and its derivatives with respect to kernel parameters $\theta_{i}$,

where $K_{\mathbf{y}}=K+\lambda I$. The traditional approach involves factorizing the kernel matrix $K_{\mathbf{y}}=LL^{\top}$ using the Cholesky algorithm (Golub & Van Loan, 1996) which costs $\mathcal{O}(n^{3})$ operations. After that, all other operations cost $\mathcal{O}(n^{2})$ except for the trace term in the calculation of $g_{i}$ which once again requires $\mathcal{O}(n^{3})$ operations. Similar computations are required for computing mean and variance predictions for test data (Rasmussen & Williams, 2006). Note that the solution of a linear system is required for computing the variance at every test point.

This approach is not viable for large $n$ and, consequently, many approaches have been proposed to approximate these computations, thus leading to approximate optimal values for $\boldsymbol{\theta}$ and approximate predictions. Here we investigate the possibility of avoiding approximations altogether, by arguing that for parameter optimization it is sufficient to obtain an unbiased estimate of the gradient $g_{i}$. In particular, when such an estimate is available, it is possible to employ stochastic gradient optimization that has strong theoretical guarantees (Robbins & Monro, 1951). In the case of GPs, the problematic terms in eq. 2 are the solution of the linear system $K_{\mathbf{y}}^{-1}\mathbf{y}$ and the trace term. In this work we make use of a stochastic linear algebra result that allows for an approximation of the trace term,

This result shows that all it takes to calculate stochastic gradients is the ability to efficiently solve linear systems. Linear systems can be iteratively solved using conjugate gradient (CG) (Golub & Van Loan, 1996). The advantage of this formulation is that we can attempt to optimize kernel parameters using stochastic gradient optimization without having to store $K_{\mathbf{y}}$ and, given that the most expensive operation is now multiplying the kernel matrix by vectors, only $\mathcal{O}(n^{2})$ computations are required. However, it is well known that the convergence of the CG algorithm depends on the condition number $\kappa(K_{\mathbf{y}})$ (ratio of largest to smallest eigenvalues), so the suitability of this approach may also be curtailed if $K_{\mathbf{y}}$ is badly conditioned. To this end, a well-known approach for improving the conditioning of a matrix, which in turn accelerates convergence, is preconditioning. This necessitates the introduction of a preconditioning matrix, $P$, which should be chosen in such a way that $P^{-1}K_{\mathbf{y}}$ approximates the identity matrix, $I$. Intuitively, this can be obtained by setting $P=K_{\mathbf{y}}$; however, given that in Preconditioned CG (PCG) we are required to solve linear systems involving $P$, this choice would be no easier than solving the original system. Thus we must choose $P$ which approximates $K_{\mathbf{y}}$ as closely as possible, but which can also be easily inverted. The PCG algorithm is shown in Algorithm 1.

2 Non-Gaussian Likelihoods

When the likelihood $p(y_{i}\mid f_{i})$ is not Gaussian, it is no longer possible to analytically integrate out latent variables. Instead, techniques such as Gaussian approximations (see, e.g., (Kuss & Rasmussen, 2005; Nickisch & Rasmussen, 2008)) and methods attempting to characterize the full posterior $p(\mathbf{f},\boldsymbol{\theta}\mid\mathbf{y})$ (Murray et al., 2010; Filippone et al., 2013) may be required. Among the various schemes to recover tractability in the case of models with a non-Gaussian likelihood, we choose the Laplace approximation, as we can formulate it in a way that only requires the solution of linear systems. The GP models we consider assume that the likelihood factorizes across all data points $p(\mathbf{y}\mid\mathbf{f})=\prod_{i=1}^{n}p(y_{i}\mid f_{i})$. The use of CG for computing the Laplace approximation has been proposed elsewhere (Flaxman et al., 2015), but we make the first use of preconditioning and stochastic gradient estimation within the Laplace approximation to compute stochastic gradients for non-conjugate models.

Defining $W=-\nabla_{\mathbf{f}}\nabla_{\mathbf{f}}\log[p(\mathbf{y}\mid\mathbf{f})]$ (a diagonal matrix), carrying out the Laplace approximation algorithm, computing its derivatives wrt $\boldsymbol{\theta}$, and making predictions, all possess the same computational bottleneck: the solution of linear systems involving the matrix $B=I+W^{\frac{1}{2}}KW^{\frac{1}{2}}$ (Rasmussen & Williams, 2006). For a given $\boldsymbol{\theta}$, each iteration of the Laplace approximation algorithm requires solving one linear system involving $B$ and two matrix-vector multiplications involving $K$; the linear system involving $B$ can be solved using CG or PCG. The Laplace approximation yields the mode $\hat{\mathbf{f}}$ of the posterior over latent variables and offers an approximate log-marginal likelihood in the form:

which poses the same computational challenges as the regression case. Once again, we therefore seek an alternative way to learn kernel parameters by stochastic gradient optimization based on computing unbiased estimates of the gradient of the approximate log-marginal likelihood. This is complicated further by the inclusion of an additional “implicit” term accounting for the change in the solution given by the Laplace approximation for a change in $\boldsymbol{\theta}$. The full derivation of the gradient is rather lengthy and is deferred to the supplementary material. Nonetheless, it is worth noting that the calculation of the exact gradient involves trace terms similar to the regression case that cannot be computed for large $n$, and we unbiasedly estimate these using the stochastic approximation of the trace.

Preconditioning Kernel Matrices

Here we consider choices for kernel preconditioners, and for the sake of clarity we focus on preconditioners for $K_{\mathbf{y}}$. Unless stated otherwise, we shall consider standard left preconditioning, whereby the original problem of solving $K_{\mathbf{y}}\mathbf{z}=\mathbf{v}$ is transformed by applying a preconditioner, $P$, to both sides of this equation. This formulation may thus be expressed as $P^{-1}K_{\mathbf{y}}\mathbf{z}=P^{-1}\mathbf{v}.$

The Nyström method was originally proposed to approximate the eigendecomposition of kernel matrices (Williams & Seeger, 2000); as a result, it offers a way to obtain a low rank approximation of $K$. This method selects a subset of $m\ll n$ data (inducing points) collected in the set $U$ which are intended for approximating the spectrum of $K$. The resulting approximation is $\hat{K}=K_{XU}K_{UU}^{-1}K_{UX}$ where $K_{UU}$ denotes the evaluation of the kernel function over the inducing points, and $K_{XU}$ denotes the evaluation of the kernel function between the input points and the inducing points. The resulting preconditioner $P=K_{XU}K_{UU}^{-1}K_{UX}+\lambda I$ can be inverted using the matrix inversion lemma

which has $\mathcal{O}(m^{3})$ complexity.

The use of a subset of data for approximating a GP kernel has also been utilized in the fully and partially independent training conditional approaches (FITC and PITC, respectively) for approximating GP regression (Candela & Rasmussen, 2005). In the former case, the prior covariance of the approximation can be written as follows:

As the name implies, this formulation enforces that the latent variables associated with $U$ are taken to be completely conditionally independent. On the other hand, the PITC method extends on this approach by enforcing that although inducing points assigned to a designated block are conditionally dependent on each other, there is no dependence between points placed in different blocks:

For the FITC preconditioner, the diagonal resulting from the training conditional can be added to the diagonal noise matrix, and the inversion lemma can be invoked as for the Nyström case. Meanwhile, for the PITC preconditioner, the noise diagonal can be added to the block diagonal matrix, which can then be inverted block-by-block. Once again, matrix inversion can then be carried out as before, where the inverted block diagonal matrix takes the place of $\lambda I$ in the original formulation.

2 Approximate factorization of kernel matrices

This group of preconditioners relies on approximations to $K$ that factorize as $\hat{K}=\Phi\Phi^{\top}$. We shall consider different ways of determining $\Phi$ such that $P$ can be inverted at a lower cost than the original kernel matrix $K$. Once again, this enables us to employ the matrix inversion lemma, and express the linear system:

We now review a few methods to approximate the kernel matrix $K$ in the form $\Phi\Phi^{\top}$.

The spectral approach uses random Fourier features for deriving a sparse approximation of a GP (Rahimi & Recht, 2008). This approach for GPs was introduced in (Lázaro-Gredilla et al., 2010), and relies on the assumption that stationary kernel functions can be represented as the Fourier transform of non-negative measures. As such, the elements of $K$ can be approximated as follows:

In the equation above, the vectors $\mathbf{s}_{r}$ denote the spectral points (or frequencies) which in the case of the RBF kernel can be sampled from $\mathcal{N}\left(\textbf{0},\frac{1}{4\pi^{2}}\Lambda\right)$, where $\Lambda=\left[1/l_{1}^{2},\dots,1/l_{n}^{2}\right]$. To the best of our knowledge, this is the first time such an approximation has been considered for the purpose of preconditioning kernel matrices.

2.2 Partial SVD

Another factorization approach that we consider in this work is the partial singular value decomposition (SVD) method (Golub & Van Loan, 1996). The SVD method factorizes the original kernel matrix $K$ into $A\Lambda A^{\top}$, where $A$ is a unitary matrix and $\Lambda$ is a diagonal matrix of singular values. Here, we shall consider a variation of this technique called randomized truncated SVD (Halko et al., 2011), which constructs an approximate low rank SVD factorization of $K$ using random sampling to accelerate computations.

2.3 Structured Kernel Interpolation (SKI)

Some recent work on approximating GPs has exploited the fast computation of Kronecker matrix-vector multiplications when inputs are located on a Cartesian grid (Gilboa et al., 2015). Unfortunately, not all datasets meet this requirement, thus limiting the widespread application of Kronecker inference. To this end, SKI (Wilson & Nickisch, 2015) is an approximation technique which exploits the benefits of the Kronecker product without imposing any requirements on the structure of the training data. In particular, a grid of inducing points, $U$, is constructed, and the covariance between the training data and $U$ is then represented as $K_{XU}=WK_{UU}.$ In this formulation, $W$ denotes a sparse interpolation matrix for assigning weights to the elements of $K_{UU}$. In this manner, a preconditioner exploiting Kronecker structure can be constructed as $P=WK_{UU}W^{\top}+\lambda I$. If we consider $V=W/\sqrt{\lambda}$, we can rewrite the (inverse) preconditioner as $P^{-1}=\lambda^{-1}(VK_{\text{UU}}V^{\top}+I)^{-1}.$ Since this can no longer be solved directly, we solve this (inner-loop) linear system using the CG algorithm (all within one iteration of the outer-loop PCG). For badly conditioned systems, although the complexity of the required matrix-vector multiplications is now much less than $\mathcal{O}({n^{2}})$, the number of iterations to solve linear systems involving the preconditioner is potentially very large, and could diminish the benefits of preconditioning.

3 Other approaches

An alternative to using a single subset of data involves constructing local GPs over segments of the original data (Snelson & Ghahramani, 2007). An example of such an approach is the Block Jacobi approximation, whereby the preconditioner is constructed by taking a block diagonal of $K$ and discarding all other elements in the kernel matrix. In this manner, covariance is only expressed for points within the same block, as $P=\textbf{bldiag}\left(K_{\mathbf{y}}+\lambda I\right).$ The inverse of this block diagonal matrix is computationally cheap (also block diagonal). However, given that a substantial amount of information contained in the original covariance matrix is ignored, this choice is intrinsically a rather crude approach.

3.2 Regularization

An appealing feature shared by the aforementioned preconditioners (aside from SKI) is that their structure enables us to directly solve $P^{-1}\mathbf{v}$. An alternative technique for constructing a preconditioner involves adding a positive regularization parameter, $\delta I$, to the original kernel matrix, such that $P=K_{\mathbf{y}}+\delta I$ (Srinivasan et al., 2014). This follows from the fact that adding noise to the diagonal of $K_{\mathbf{y}}$ makes it better-conditioned, and the condition number is expected to decrease further as $\delta$ increases. Nonetheless, for the purpose of preconditioning, this parameter should be tuned in such a way that $P$ remains a sensible approximation of $K_{\mathbf{y}}$. As opposed to the previous preconditioners, this is an instance of right preconditioning, which has the following general form $K_{\mathbf{y}}P^{-1}(P\textbf{x})=\textbf{v}.$

Given that it is no longer possible to evaluate $P^{-1}\mathbf{v}$ analytically, this linear system is solved yet again using CG, such that a linear system of equations is solved at every outer iteration of the PCG algorithm. Due to the potential loss of accuracy incurred while solving the inner linear systems, a variation of the standard PCG algorithm, referred to as flexible PCG (Notay, 2000), is used instead. Using this approach, a re-orthogonalization step is introduced such that the search directions remain orthogonal even when the inner system is not solved to high precision.

Comparison of Preconditioners

In this section, we provide an empirical exploration of these preconditioners in a practical setting. We begin by considering three datasets for regression from the UCI repository (Asuncion & Newman, 2007), namely the Concrete dataset ($n=1030,d=8$), the Power Plant dataset ($n=9568,d=4$), and the Protein dataset ($n=45730,d=9$). In particular, we evaluate the convergence in solving $K_{\mathbf{y}}\mathbf{z}=\mathbf{y}$ using iterative methods, where $\mathbf{y}$ denotes the labels of the designated dataset, and $K_{\mathbf{y}}$ is constructed using different configurations of kernel parameters.

With this experiment, we aim to assess the quality of different preconditioners based on how many matrix-vector products they require, which, for most approaches, corresponds to the number of iterations taken by PCG to converge. The convergence threshold is set to $\epsilon^{2}=n\cdot 10^{-10}$ so as to roughly accept an average error of $10^{-5}$ on each element of the solution.

For every variation, we set the parameters of the preconditioners so as to have a complexity lower than the $\mathcal{O}(n^{2})$ cost associated with matrix-vector products; by doing so, we can assume that the latter computations are the dominant cost for large $n$. In particular, for Nyström-type methods, we set $m=\sqrt{n}$ inducing points, so that when we invert the preconditioner using the matrix inversion lemma, the cost is in $\mathcal{O}(m^{3})=\mathcal{O}(n^{3/2})$. Similarly, for the Spectral preconditioner, we set $m=\sqrt{n}$ random features. For the SKI preconditioner, we take an equal number of elements on the grid for each dimension; under this assumption, Kronecker products have $\mathcal{O}(dn^{\frac{d+1}{d}})$ cost (Gilboa et al., 2015), and we set the size of the grid so that the complexity of applying the preconditioner matches $\mathcal{O}(n^{3/2})$, so as to be consistent with the other preconditioners. For the Regularized approach, each iteration needed to apply the preconditioner requires one matrix-vector product, and we add this to the overall count of such computations. For this preconditioner, we add a diagonal offset $\delta$ to the original matrix, equivalent to two orders of magnitude greater than the noise of the process. In general, although the complexity of PCG is indeed no different from that of CG, we emphasize that experiencing a 2-fold or 5-fold (in some cases even an order of magnitude) improvement can be very substantial when plain CG takes very long to converge or when the dataset is large.

We focus on an isotropic RBF variant of the kernel in eq. 1, fixing the marginal variance $\sigma^{2}$ to one. We vary the length-scale parameter $l$ and the noise variance $\lambda$ in $\log_{10}$ scale. The top part of fig. 1 shows the number of iterations that the standard CG algorithm takes, where we have capped the number of iterations to 100,000.

The bottom part of the figure reports the improvement offered by various preconditioners measured as

It is worth noting that when both CG and PCG fail to converge within the upper bound, the improvement will be marked as 0, i.e. neither a gain or a loss within the given bound. The results plotted in fig. 1 indicate that the low-rank preconditioners (PITC, FITC and Nyström) achieve significant reductions in the number of iterations for each dataset, and all approaches work best when the lengthscale is longer, characterising smoother processes. In contrast, preconditioning seems to be less effective when the lengthscale is shorter, corresponding to a kernel matrix that is more sparse. However, for cases yielding positive results, the improvement is often in the range of an order of magnitude, which can be substantial when a large number of iterations is required by the CG algorithm.

The results also confirm that, as alluded to in the previous section, Block Jacobi preconditioning is generally a poor preconditioner, particularly when the corresponding kernel matrix is dense. The only minor improvements were observed when CG itself converges quickly, in which case preconditioning serves very little purpose either way.

The regularization approach with flexible conjugate gradient does not appear to be effective in any case either, particularly due to the substantial amount of iterations required for solving an inner system at every iteration of the PCG algorithm. This implies that introducing additional small jitter to the diagonal does not necessarily make the system much easier to solve, whilst adding an overly large offset would negatively impact convergence of the outer algorithm. One could assume that tuning the value of this parameter could result in slightly better results; however, preliminary experiments in this regard yielded only minor improvements.

The results for SKI preconditioning are similarly discouraging at face value. When the matrix $K_{\mathbf{y}}$ is very badly conditioned, an excessive number of inner iterations are required for every iteration of outer PCG. This greatly increases the duration of solving such systems, and as a result, this method was not included in the comparison for the Protein dataset, where it was evident that preconditioning the matrix in this manner would not yield satisfactory improvements. Notwithstanding that these experiments depict a negative view of SKI preconditioning, it must be said that we assumed a fairly simplistic interpolation procedure in our experiments, where each data point was mapped to nearest grid location. The size of the constructed grid is also hindered considerably by the constraint imposed by our upper bound on complexity. Conversely, more sophisticated interpolation strategies or even grid formulation procedures could possibly speed up the convergence of CG for the inner systems. In line with this thought, however, one could argue that the preconditioner would no longer be straightforward to construct, which goes against our innate preference towards easily derived preconditioners.

Impact of preconditioning on GP learning

One of the primary objectives of this work is to reformulate GP regression and classification in such a way that preconditioning can be effectively exploited. In section 2, we demonstrated how preconditioning can indeed be applied to GP regression problems, and also proposed a novel way of rewriting GP classification in terms of solving linear systems (where preconditioning can thus be employed). We can now evaluate how the proposed preconditioned GP techniques compare to other state of the art methods.

To this end, in this section, we empirically report on the generalization ability of GPs as a function of the time taken to optimize parameters $\boldsymbol{\theta}$ and compute predictions. In particular, for each of the methods featured in our comparison, we iteratively run the optimization of kernel parameters for a few iterations and predict on unseen data, and assess how prediction accuracy varies over time for different methods.

We can make use of stochastic gradients for GP models to optimize kernel parameters using off-the-shelf stochastic gradient optimization algorithms. In order to reduce the number of parameters to tune, we employ ADAGRAD (Duchi et al., 2011) – an optimization algorithm having a single step-size parameter. For the purpose of this experiment, we do not attempt to optimize this parameter, since this would require additional computations. Nonetheless, our experience with training GP models indicates that the choice of this parameter is not critical: we set the step-size to one.

Fig. 2 shows the two error measures over time for a selection of approaches. In the figure, PCG and CG refer to stochastic gradient optimization of kernel parameters using ADAGRAD, where linear systems are solved with PCG and CG, respectively. In view of the results obtained in our comparison of preconditioners, we decide to proceed with the Nyström preconditioning method. Furthermore, we construct the preconditioner with $m=4\sqrt{n}$ points randomly selected from the input data at each iteration, such that the overall complexity of the PCG method matches plain CG. For these methods, stochastic estimates of trace terms are carried out using $N_{\mathbf{r}}=4$ random vectors. The baseline CHOL method refers to the optimization of kernel parameters using the L-BFGS algorithm, where the exact log-marginal likelihood and its gradient are calculated using the full Cholesky decomposition of $K_{\mathbf{y}}$ or $B$.

Alongside these approaches for optimizing kernel parameters without approximation, we also evaluate the performance of approximate GP methods. For this experiment, we chose to compare against approximations found in the software package GPstuff (Vanhatalo et al., 2013), namely the fully and partial independent training conditional approaches (FITC, PITC), and the sparse variational GP (VAR) (Titsias, 2009). In order to match the computational cost of CG/PCG, which is in $\mathcal{O}(n^{2})$, we set the number of inducing points for the approximate methods to be $n^{2/3}$.

All methods are initialized from the same set of kernel parameters, and the curves are averaged over $5$ folds ($3$ for the Protein and EEG datasets). For the sake of integrity, we ran each method in the comparison individually on a workstation with Intel Xeon E5-2630 CPU having 16 cores and 128GB RAM. We also ensured that all methods reported in the comparison used optimized linear algebra routines exploiting the multi-core architecture. This diligence for ensuring fairness gives credence to our assumption that the timings are not affected by external factors other than the actual implementation of the algorithms. The CG, PCG and CHOL approaches have been implemented in R; the fact that the approximate methods were implemented in a different environment (GPstuff is written in Matlab/Octave) and by a different developer may cast some doubt on the correctness of directly comparing results. However, we believe that the key point emerging from this comparison is that preconditioning feasibly enables the use of iterative approaches for optimization of kernel parameters in GPs, and the results are competitive with those achieved using popular GP software packages.

For the reported experiments, it was possible to store the kernel matrix $K$ for all datasets, making it possible to compare methods against the baseline GP where computations use Cholesky decompositions. We stress, however, that iterative approaches based on CG/PCG can be implemented without the need to store $K$, whereas this is not possible for approaches that attempt to factorize $K$ exactly. It is also worth noting that for the CG/PCG approach, calculating the log-likelihood on test data requires solving one linear system for each test point; this clearly penalizes the speed of these methods given the set-up of the experiment, where predictions are carried out every fixed number of iterations.

Discussion and Conclusions

Careful attention to numerical properties is essential in scaling machine learning to large and realistic datasets. Here we have introduced the use of preconditioning to the implementation of kernel machines, specifically, prediction and learning of kernel parameters for GPs. Our novel scheme permits the use of any likelihood that factorizes over the data points, allowing us to tackle both regression and classification. We have shown robust performance improvements, in both accuracy and computational cost, over a host of state-of-the-art approximation methods for kernel machines. Notably, our method is exact in the limit of iterations, unlike approximate alternatives. We have also shown that the use of PCG is competitive with exact Cholesky decomposition in modestly sized datasets, when the Cholesky factors can be feasibly computed. When data and thus the kernel matrix grow large enough, Cholesky factorization becomes unfeasible, leaving PCG as the optimal choice.

One of the key features of a PCG implementation is that it does not require storage of any $\mathcal{O}(n^{2})$ objects. We plan to extend our implementation to compute the elements of $K$ on the fly in one case, and in another case store $K$ in a distributed fashion (e.g. in TensorFlow/Spark). Furthermore, while we have focused on solving linear systems, we can also use preconditioning for other iterative algorithms involving the $K$ matrix, e.g., those to solve $\log(K)\mathbf{v}$ and $K^{1/2}\mathbf{v}$ (Chen et al., 2011), as is often useful in estimating marginal likelihoods for probabilistic kernel models like GPs.

Acknowledgements

KC and MF are grateful to Pietro Michiardi and Daniele Venzano for assisting the completion of this work by providing additional computational resources for running the experiments. JPC acknowledges support from the Sloan Foundation, The Simons Foundation (SCGB#325171 and SCGB#325233), and The Grossman Center at Columbia University.

References

Appendix A Other results not included in the paper

In fig. 3 we report some of the runs that we did not include in the main text for lack of space. The figure reports plots on the error vs. time for the same regression cases considered in the main text but with an isotropic kernel, and results on the concrete dataset with isotropic and ARD kernels.

Appendix B Gaussian Processes with non-Gaussian likelihood functions

In this section we report the derivations of the quantities needed to compute an unbiased estimate of the log-marginal likelihood given by the Laplace approximation for GP models with non-Gaussian likelihood functions. Throughout this section, we assume a factorizing likelihood

and we specialize the equations to the probit likelihood

where $\Phi$ denotes the cumulative function of the Gaussian density. The latent variables $\mathbf{f}$ are given a zero mean GP prior $\mathbf{f}\sim\mathcal{N}(\mathbf{f}|\mathbf{0},K)$.

For a given value of the hyperparameters $\boldsymbol{\theta}$, define

as the logarithm of the posterior density over $\mathbf{f}$. Performing a Laplace approximation amounts in defining a Gaussian $q(\mathbf{f}\mid\mathbf{y},\boldsymbol{\theta})=\mathcal{N}(\mathbf{f}\mid\hat{\mathbf{f}},\hat{\Sigma})$, such that

As it is not possible to directly solve the maximization problem in equation 5, an iterative procedure based on the following Newton-Raphson formula is usually employed,

starting from some initial $\mathbf{f}$ until convergence. The gradient and the Hessian of the log of the target density are

where we have defined $W=-\nabla_{\mathbf{f}}\nabla_{\mathbf{f}}\log[p(\mathbf{y}\mid\mathbf{f})]$, which is diagonal because the likelihood factorizes over observations. Note that if $\log[p(\mathbf{y}\mid\mathbf{f})]$ is concave, such as in probit classification, $\Psi(\mathbf{f})$ has a unique maximum.

We can rewrite the inverse of the negative Hessian using the matrix inversion lemma:

We can define $\mathbf{b}=(W\mathbf{f}+\nabla_{\mathbf{f}}\log[p(\mathbf{y}\mid\mathbf{f})])$ and rewrite this expression as:

As we will see later, the definition of $\mathbf{a}$ is useful for the calculation of the gradient and for predictions.

Proceeding with the calculations from right to left we see that in order to complete a Newton-Raphson iteration the expensive operations are: (i) carry out one matrix-vector multiplication $K\mathbf{b}$, (ii) solve a linear system involving the $B$ matrix, and (iii) carry out one matrix-vector multiplication involving $K$ and the vector in the parenthesis. Calculating $\mathbf{b}$ and performing any multiplications of $W^{\frac{1}{2}}$ with vectors cost $\mathcal{O}(n)$.

All these operations can be carried out without the need to store $K$ or any other $n\times n$ matrices. The linear system in (ii) can be solved using the CG algorithm that involves repeatedly multiplying $B$ (and therefore $K$) with vectors.

The Laplace approximation yields an approximate log-marginal likelihood in the following form:

Handy relationships that we will be using in the remainder of this section are:

The gradient of the log-marginal likelihood with respect to the kernel parameters $\boldsymbol{\theta}$ requires differentiating the terms that explicitly depend on $\boldsymbol{\theta}$ and those that implicitly depend on it because a change in the parameters reflects in a change in $\hat{\mathbf{f}}$. Denoting by $g_{i}$ the $i$th component of the gradient of $\frac{\partial\log[\hat{p}(\mathbf{y}|\boldsymbol{\theta})]}{\partial\theta_{i}}$, we obtain

The trace term cannot be computed exactly for large $n$ so we propose a stochastic estimate:

By noticing that the derivative of $B$ is $W^{\frac{1}{2}}\frac{\partial K}{\partial\theta_{i}}W^{\frac{1}{2}}$, this simplifies to

so we need to solve $N_{\mathbf{r}}$ linear systems involving $B$.

The second term contains the linear system $K^{-1}\hat{\mathbf{f}}$ that we already have from the Laplace approximation and is $\mathbf{a}$.

The third term is slightly more involved and will be dealt with in the next sub-section.

The last (implicit) term in the last equation can be simplified by noticing that:

and that the derivative of the first term wrt $\hat{\mathbf{f}}$ is zero because $\hat{\mathbf{f}}$ maximizes $\Psi(\hat{\mathbf{f}})$. Therefore:

The components of $\left[\nabla_{\hat{\mathbf{f}}}\log|B|\right]$ can be obtained by considering the identity $\log|B|=\log|I+KW|$, so differentiating $\log|B|$ wrt the components of $\hat{\mathbf{f}}$ becomes:

We can rewrite this by gathering $K$ inside the inverse and, due to the inversion of the matrix product, $K$ cancels out:

We notice here that the resulting trace contains the inverse of the same matrix needed in the iterations of the Laplace approximation and that the matrix $\frac{\partial W}{\partial(\hat{\mathbf{f}})_{j}}$ is zero everywhere except in the $j$th diagonal element where it attains the value:

For this reason, it would be possible to simplify the trace term as the product between the $j$th diagonal element of $(K^{-1}+W)^{-1}$ and $\frac{\partial^{3}\log[p(\mathbf{y}\mid\hat{\mathbf{f}})]}{\partial(\hat{\mathbf{f}})^{3}_{j}}$. Bearing in mind that we need $n$ of these quantities, we could define

which is the standard way to proceed when computing the gradient of the approximate log-marginal likelihood using the Laplace approximation (Rasmussen & Williams, 2006). However, this would be difficult to compute exactly for large $n$, as this would require inverting $K^{-1}+W$ first and then compute its diagonal. Using the matrix inversion lemma would not simplify things as there would still be an inverse of $B$ to compute explicitly. We therefore aim for a stochastic estimate of this term starting from:

which requires solving $N_{\mathbf{r}}$ linear systems involving the $B$ matrix:

The derivative of $\hat{\mathbf{f}}$ wrt $\theta_{i}$ can be obtained by differentiating the expression $\hat{\mathbf{f}}=K\nabla_{\hat{\mathbf{f}}}\log[p(\mathbf{y}\mid\hat{\mathbf{f}})]$:

Given that $\nabla_{\hat{\mathbf{f}}}\nabla_{\hat{\mathbf{f}}}\log[p(\mathbf{y}\mid\hat{\mathbf{f}})]=-W$ we can rewrite:

So an unbiased estimate of the implicit term in the gradient of the approximate log-marginal likelihood becomes:

Rewriting the inverse in terms of $B$ yields:

Putting everything together, the components of the stochastic gradient are:

B.2 Predictions

To obtain an approximate predictive distribution, conditioned on a value of the hyperparameters $\boldsymbol{\theta}$, we can compute:

Given the properties of multivariate normal variables, $f_{*}$ is distributed as $\mathcal{N}(f_{*}\mid\mu_{*},\beta^{2}_{*})$ with $\mu_{*}=\mathbf{k}_{*}^{\top}K^{-1}\mathbf{f}$ and $\beta^{2}_{*}=k_{**}-\mathbf{k}_{*}^{\top}K^{-1}\mathbf{k}_{*}$. Approximating $p(\mathbf{f}\mid\mathbf{y},\boldsymbol{\theta})$ with a Gaussian $q(\mathbf{f}\mid\mathbf{y},\boldsymbol{\theta})=\mathcal{N}(\mathbf{f}\mid\boldsymbol{\mu}_{q},\Sigma_{q})$ makes it possible to analytically perform integration with respect to $\mathbf{f}$ in eq. 14. In particular, the integration with respect to $\mathbf{f}$ yields $\mathcal{N}(f_{*}\mid m_{*},s^{2}_{*})$ with

This shows that the mean is cheap to compute, whereas the variance requires solving another linear system involving $B$ for each test point.

The univariate integration with respect to $f_{*}$ follows exactly in the case of a probit likelihood, as it is a convolution of a Gaussian and a cumulative Gaussian

B.3 Low rank preconditioning

When a low rank approximation of the matrix $K$ is available, say $\hat{K}=\Phi\Phi^{\top}$, the inverse of the preconditioner can be rewritten as:

By using the matrix inversion lemma we obtain:

Similarly to the GP regression case, the application of this preconditioner is in $\mathcal{O}(m^{3})$, where $m$ is the rank of $\Phi$.