Scalable Log Determinants for Gaussian Process Kernel Learning

11/09/2017 ∙ by Kun Dong, et al. ∙ 0

For applications as varied as Bayesian neural networks, determinantal point processes, elliptical graphical models, and kernel learning for Gaussian processes (GPs), one must compute a log determinant of an n × n positive definite matrix, and its derivatives - leading to prohibitive O(n^3) computations. We propose novel O(n) approaches to estimating these quantities from only fast matrix vector multiplications (MVMs). These stochastic approximations are based on Chebyshev, Lanczos, and surrogate models, and converge quickly even for kernel matrices that have challenging spectra. We leverage these approximations to develop a scalable Gaussian process approach to kernel learning. We find that Lanczos is generally superior to Chebyshev for kernel learning, and that a surrogate approach can be highly efficient and accurate with popular kernels.



page 1

page 2

page 3

page 4

Code Repositories


Scalable Log Determinants for Gaussian Process Kernel Learning ( (NIPS 2017)

view repo
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

There is a pressing need for scalable machine learning approaches to extract rich statistical structure from large datasets. A common bottleneck — arising in determinantal point processes

(Kulesza et al., 2012), Bayesian neural networks (MacKay, 1992), model comparison (MacKay, 2003), graphical models (Rue and Held, 2005), and Gaussian process kernel learning (Rasmussen and Williams, 2006) — is computing a log determinant over a large positive definite matrix. While we can approximate log determinants by existing stochastic expansions relying on matrix vector multiplications (MVMs), these approaches make assumptions, such as near-uniform eigenspectra (Boutsidis et al., 2015)

, which are unsuitable in machine learning contexts. For example, the popular RBF kernel gives rise to rapidly decaying eigenvalues. Moreover, while standard approaches, such as stochastic power series, have reasonable asymptotic complexity in the rank of the matrix, they require too many terms (MVMs) for the precision necessary in machine learning applications.

Gaussian processes (GPs) provide a principled probabilistic kernel learning framework, for which a log determinant is of foundational importance. Specifically, the marginal likelihood

of a Gaussian process is the probability of data given only kernel hyper-parameters. This utility function for kernel learning compartmentalizes into automatically calibrated model fit and complexity terms — called

automatic Occam’s razor — such that the simplest models which explain the data are automatically favoured (Rasmussen and Ghahramani, 2001; Rasmussen and Williams, 2006)

, without the need for approaches such as cross-validation, or regularization, which can be costly, heuristic, and involve substantial hand-tuning and human intervention. The automatic complexity penalty, called the

Occam’s factor (MacKay, 2003), is a log determinant of a kernel (covariance) matrix, related to the volume of solutions that can be expressed by the Gaussian process.

Many current approaches to scalable Gaussian processes (e.g., Quiñonero-Candela and Rasmussen, 2005; Le et al., 2013; Hensman et al., 2013) focus on inference assuming a fixed kernel, or use approximations that do not allow for very flexible kernel learning (Wilson, 2014), due to poor scaling with number of basis functions or inducing points. Alternatively, approaches which exploit algebraic structure in kernel matrices can provide highly expressive kernel learning (Wilson et al., 2014), but are essentially limited to grid structured data.

Recently, Wilson and Nickisch (2015) proposed the

structured kernel interpolation

(SKI) framework, which generalizes structuring exploiting methods to arbitrarily located data. SKI works by providing accurate and fast matrix vector multiplies (MVMs) with kernel matrices, which can then be used in iterative solvers such as linear conjugate gradients for scalable GP inference. However, evaluating the marginal likelihood and its derivatives, for kernel learning, has followed a scaled eigenvalue approach (Wilson et al., 2014; Wilson and Nickisch, 2015) instead of iterative MVM approaches. This approach can be inaccurate, and relies on a fast eigendecomposition of a structured matrix, which is not available in many consequential situations where fast MVMs are available, including: (i) additive covariance functions, (ii) multi-task learning, (iii) change-points (Herlands et al., 2016), and (iv) diagonal corrections to kernel approximations (Snelson and Ghahramani, 2006). Fiedler (Fiedler, 1984) and Weyl (Weyl, 1912) bounds have been used to extend the scaled eigenvalue approach (Flaxman et al., 2015; Herlands et al., 2016), but are similarly limited. These extensions are often very approximate, and do not apply beyond sums of two and three matrices, where each matrix in the sum must have a fast eigendecomposition.

In machine learning there has recently been renewed interest in MVM based approaches to approximating log determinants, such as the Chebyshev (Han et al., 2015) and Lanczos (Ubaru et al., ) based methods, although these approaches go back at least two decades in quantum chemistry computations (Bai et al., 1998). Independently, several authors have proposed various methods to compute derivatives of log determinants (MacKay and Gibbs, 1997; Stein et al., 2013). But both the log determinant and the derivatives are needed for efficient GP marginal likelihood learning: the derivatives are required for gradient-based optimization, while the log determinant itself is needed for model comparison, comparisons between the likelihoods at local maximizers, and fast and effective choices of starting points and step sizes in a gradient-based optimization algorithm.

In this paper, we develop novel scalable and general purpose Chebyshev, Lanczos, and surrogate approaches for efficiently and accurately computing both the log determinant and its derivatives simultaneously. Our methods use only fast MVMs, and re-use the same MVMs for both computations. In particular:

  • We derive fast methods for simultaneously computing the log determinant and its derivatives by stochastic Chebyshev, stochastic Lanczos, and surrogate models, from MVMs alone. We also perform an error analysis and extend these approaches to higher order derivatives.

  • These methods enable fast GP kernel learning whenever fast MVMs are possible, including applications where alternatives such as scaled eigenvalue methods (which rely on fast eigendecompositions) are not, such as for (i) diagonal corrections for better kernel approximations, (ii) additive covariances, (iii) multi-task approaches, and (iv) non-Gaussian likelihoods.

  • We illustrate the performance of our approach on several large, multi-dimensional datasets, including a consequential crime prediction problem, and a precipitation problem with training points. We consider a variety of kernels, including deep kernels (Wilson et al., 2016a), diagonal corrections, and both Gaussian and non-Gaussian likelihoods.

  • We have released code and tutorials as an extension to the GPML library (Rasmussen and Nickisch, 2010) at A Python implementation of our approach is also available through the GPyTorch library:

When using our approach in conjunction with SKI (Wilson and Nickisch, 2015) for fast MVMs, GP kernel learning is , for inducing points and training points, where . With algebraic approaches such as SKI we also do not need to worry about quadratic storage in inducing points, since symmetric Toeplitz and Kronecker matrices can be stored with at most linear cost, without needing to explicitly construct a matrix.

Although we here use SKI for fast MVMs, we emphasize that the proposed iterative approaches are generally applicable, and can easily be used in conjunction with any method that admits fast MVMs, including classical inducing point methods (Quiñonero-Candela and Rasmussen, 2005), finite basis expansions (Le et al., 2013), and the popular stochastic variational approaches (Hensman et al., 2013). Moreover, stochastic variational approaches can naturally be combined with SKI to further accelerate MVMs (Wilson et al., 2016b).

We start in §2 with an introduction to GPs and kernel approximations. In §3 we introduce stochastic trace estimation and Chebyshev (§3.1) and Lanczos (§3.2) approximations. In §4, we describe the different sources of error in our approximations. In §5 we consider experiments on several large real-world data sets. We conclude in §6. The supplementary materials also contain several additional experiments and details.

2 Background

A Gaussian process (GP) is a collection of random variables, any finite number of which have a joint Gaussian distribution

(e.g., Rasmussen and Williams, 2006). A GP can be used to define a distribution over functions , where each function value is a random variable indexed by , and and are the mean and covariance functions of the process.

The covariance function is often chosen to be an RBF or Matérn kernel (see the supplementary material for more details). We denote any kernel hyperparameters by the vector

. To be concise we will generally not explicitly denote the dependence of and associated matrices on .

For any locations , where and represent the vectors of function values for and evaluated at each of the , and is the matrix whose entry is . Suppose we have a vector of corresponding function values

, where each entry is contaminated by independent Gaussian noise with variance

. Under a Gaussian process prior depending on the covariance hyperparameters , the log marginal likelihood is given by


where and . Optimization of (1) is expensive, since the cheapest way of evaluating and its derivatives without taking advantage of the structure of involves computing the Cholesky factorization of . computations is too expensive for inference and learning beyond even just a few thousand points.

A popular approach to GP scalability is to replace the exact kernel by an approximate kernel that admits fast computations Quiñonero-Candela and Rasmussen (2005). Several methods approximate via inducing points . An example is the subset of regressor (SoR) kernel:

which is a low-rank approximation Silverman (1985). The SoR matrix has rank at most , allowing us to solve linear systems involving and to compute in time. Another popular kernel approximation is the fully independent training conditional (FITC), which is a diagonal correction of SoR so that the diagonal is the same as for the original kernel Snelson and Ghahramani (2006). Thus kernel matrices from FITC have low-rank plus diagonal structure. This modification has had exceptional practical significance, leading to improved point predictions and much more realistic predictive uncertainty Quiñonero-Candela and Rasmussen (2005); Quinonero-Candela et al. (2007), making FITC arguably the most popular approach for scalable Gaussian processes.

Wilson and Nickisch (2015) provides a mechanism for fast MVMs through proposing the structured kernel interpolation (SKI) approximation,


where is an -by- matrix of interpolation weights; the authors of Wilson and Nickisch (2015) use local cubic interpolation so that is sparse. The sparsity in makes it possible to naturally exploit algebraic structure (such as Kronecker or Toeplitz structure) in when the inducing points are on a grid, for extremely fast matrix vector multiplications with the approximate even if the data inputs are arbitrarily located. For instance, if is Toeplitz, then each MVM with the approximate costs only . By contrast, placing the inducing points on a grid for classical inducing point methods, such as SoR or FITC, does not result in substantial performance gains, due to the costly cross-covariance matrices and .

3 Methods

Our goal is to estimate, for a symmetric positive definite matrix ,

where is the matrix logarithm Higham (2008). We compute the traces involved in both the log determinant and its derivative via stochastic trace estimators Hutchinson (1990), which approximate the trace of a matrix using only matrix vector products.

The key idea is that for a given matrix and a random probe vector with independent entries with mean zero and variance one, then ; a common choice is to let the entries of the probe vectors be Rademacher random variables. In practice, we estimate the trace by the sample mean over independent probe vectors. Often surprisingly few probe vectors suffice.

To estimate , we need to multiply by probe vectors. We consider two ways to estimate : by a polynomial approximation of or by using the connection between the Gaussian quadrature rule and the Lanczos method Han et al. (2015); Ubaru et al.

. In both cases, we show how to re-use the same probe vectors for an inexpensive coupled estimator of the derivatives. In addition, we may use standard radial basis function interpolation of the log determinant evaluated at a few systematically chosen points in the hyperparameter space as an inexpensive surrogate for the log determinant.

3.1 Chebyshev

Chebyshev polynomials are defined by the recursion

For the Chebyshev interpolant of degree is

where is the Kronecker delta and for ; see Gil et al. (2007). Using the Chebyshev interpolant of , we approximate by

when has eigenvalues .

For stochastic estimation of , we only need to compute for each given probe vector . We compute vectors and via the coupled recurrences

This gives the estimators

Thus, each derivative of the approximation costs two extra MVMs per term.

3.2 Lanczos

We can also approximate via a Lanczos decomposition; see Golub and Meurant (2010) for discussion of a Lanczos-based computation of and Ubaru et al. ; Bai et al. (1998) for stochastic Lanczos estimation of log determinants. We run steps of the Lanczos algorithm, which computes the decomposition

where is a matrix with orthonormal columns such that , is tridiagonal, is the residual, and is the th Cartesian unit vector. We estimate


where is the first column of the identity. The Lanczos algorithm is numerically unstable. Several practical implementations resolve this issue Cullum and Willoughby (2002); Saad (1992). The approximation (3) corresponds to a Gauss quadrature rule for the Riemann-Stieltjes integral of the measure associated with the eigenvalue distribution of . It is exact when is a polynomial of degree up to . This approximation is also exact when has at most distinct eigenvalues, which is particularly relevant to Gaussian process regression, since frequently the kernel matrices only have a small number of eigenvalues that are not close to zero.

The Lanczos decomposition also allows us to estimate derivatives of the log determinant at minimal cost. Via the Lanczos decomposition, we have

This approximation requires no additional matrix vector multiplications beyond those used to compute the Lanczos decomposition, which we already used to estimate ; in exact arithmetic, this is equivalent to steps of CG. Computing in this way takes additional time; subsequently, we only need one matrix-vector multiply by for each probe vector to estimate .

3.3 Diagonal correction to SKI

The SKI approximation may provide a poor estimate of the diagonal entries of the original kernel matrix for kernels with limited smoothness, such as the Matérn kernel. In general, diagonal corrections to scalable kernel approximations can lead to great performance gains. Indeed, the popular FITC method (Snelson and Ghahramani, 2006) is exactly a diagonal correction of subset of regressors (SoR).

We thus modify the SKI approximation to add a diagonal matrix ,


such that the diagonal of the approximated is exact. In other words, substracts the diagonal of and adds the true diagonal of . This modification is not possible for the scaled eigenvalue method for approximating log determinants in (Wilson and Nickisch, 2015), since adding a diagonal matrix makes it impossible to approximate the eigenvalues of from the eigenvalues of .

However, Eq. (4) still admits fast MVMs and thus works with our approach for estimating the log determinant and its derivatives. Computing with SKI costs only flops since is sparse for local cubic interpolation. We can therefore compute in flops.

3.4 Estimating higher derivatives

We have already described how to use stochastic estimators to compute the log marginal likelihood and its first derivatives. The same approach applies to computing higher-order derivatives for a Newton-like iteration, to understand the sensitivity of the maximum likelihood parameters, or for similar tasks. The first derivatives of the full log marginal likelihood are

and the second derivatives of the two terms are

Superficially, evaluating the second derivatives would appear to require several additional solves above and beyond those used to estimate the first derivatives of the log determinant. In fact, we can get an unbiased estimator for the second derivatives with no additional solves, but only fast products with the derivatives of the kernel matrices. Let

and be independent probe vectors, and define and . Then

Hence, if we use the stochastic Lanczos method to compute the log determinant and its derivatives, the additional work required to obtain a second derivative estimate is one MVM by each second partial of the kernel for each probe vector and for , one MVM of each first partial of the kernel with , and a few dot products.

3.5 Radial basis functions

Another way to deal with the log determinant and its derivatives is to evaluate the log determinant term at a few systematically chosen points in the space of hyperparameters and fit an interpolation approximation to these values. This is particularly useful when the kernel depends on a modest number of hyperparameters (e.g., half a dozen), and thus the number of points we need to precompute is relatively small. We refer to this method as a surrogate, since it provides an inexpensive substitute for the log determinant and its derivatives. For our surrogate approach, we use radial basis function (RBF) interpolation with a cubic kernel and a linear tail. See e.g. Buhmann (2000); Fasshauer (2007); Schaback and Wendland (2006); Wendland (2004) and the supplementary material for more details on RBF interpolation.

4 Error properties

In addition to the usual errors from sources such as solver termination criteria and floating point arithmetic, our approach to kernel learning involves several additional sources of error: we approximate the true kernel with one that enables fast MVMs, we approximate traces using stochastic estimation, and we approximate the actions of

and on probe vectors.

We can compute first-order estimates of the sensitivity of the log likelihood to perturbations in the kernel using the same stochastic estimators we use for the derivatives with respect to hyperparameters. For example, if is the likelihood for a reference kernel , then

and we can bound the change in likelihood at first order by . Given bounds on the norms of , we can similarly estimate changes in the gradient of the likelihood, allowing us to bound how the marginal likelihood hyperparameter estimates depend on kernel approximations.

If , the Hutchinson trace estimator has known variance Avron and Toledo (2011)

If the eigenvalues of the kernel matrix without noise decay rapidly enough compared to , the variance will be small compared to the magnitude of . Hence, we need fewer probe vectors to obtain reasonable accuracy than one would expect from bounds that are blind to the matrix structure. In our experiments, we typically only use 5–10 probes — and we use the sample variance across these probes to estimate a posteriori the stochastic component of the error in the log likelihood computation. If we are willing to estimate the Hessian of the log likelihood, we can increase rates of convergence for finding kernel hyperparameters.

The Chebyshev approximation scheme requires steps to obtain an approximation error in computing , where is the condition number of  Han et al. (2015). This behavior is independent of the distribution of eigenvalues within the interval , and is close to optimal when eigenvalues are spread quasi-uniformly across the interval. Nonetheless, when the condition number is large, convergence may be quite slow. The Lanczos approach converges at least twice as fast as Chebyshev in general (Ubaru et al., , Remark 1), and converges much more rapidly when the eigenvalues are not uniform within the interval, as is the case with log determinants of many kernel matrices. Hence, we recommend the Lanczos approach over the Chebyshev approach in general. In all of our experiments, the error associated with approximating by Lanczos was dominated by other sources of error.

5 Experiments

We test our stochastic trace estimator with both Chebyshev and Lanczos approximation schemes on: (1) a sound time series with missing data, using a GP with an RBF kernel; (2) a three-dimensional space-time precipitation data set with over half a million training points, using a GP with an RBF kernel; (3) a two-dimensional tree growth data set using a log-Gaussian Cox process model with an RBF kernel; (4) a three-dimensional space-time crime datasets with a log-Gaussian Cox model with Matérn 3/2 and spectral mixture kernels; and (5) a high-dimensional feature space using the deep kernel learning framework (Wilson et al., 2016a). In the supplementary material we also include several additional experiments to illustrate particular aspects of our approach, including kernel hyperparameter recovery, diagonal corrections (Section 3.3), and surrogate methods (Section 3.5). Throughout we use the SKI method (Wilson and Nickisch, 2015) of Eq. (2) for fast MVMs. We find that the Lanczos and surrogate methods are able to do kernel recovery and inference significantly faster and more accurately than competing methods.

5.1 Natural sound modeling

Here we consider the natural sound benchmark in Wilson and Nickisch (2015), shown in Figure 1(a). Our goal is to recover contiguous missing regions in a waveform with training points. We exploit Toeplitz structure in the matrix of our SKI approximate kernel for accelerated MVMs.

The experiment in Wilson and Nickisch (2015) only considered scalable inference and prediction, but not hyperparameter learning, since the scaled eigenvalue approach requires all the eigenvalues for an Toeplitz matrix, which can be computationally prohibitive with cost . However, evaluating the marginal likelihood on this training set is not an obstacle for Lanczos and Chebyshev since we can use fast MVMs with the SKI approximation at a cost of .

In Figure 1(b), we show how Lanczos, Chebyshev and surrogate approaches scale with the number of inducing points

compared to the scaled eigenvalue method and FITC. We use 5 probe vectors and 25 iterations for Lanczos, both when building the surrogate and for hyperparameter learning with Lanczos. We also use 5 probe vectors for Chebyshev and 100 moments. Figure

1(b) shows the runtime of the hyperparameter learning phase for different numbers of inducing points , where Lanczos and the surrogate are clearly more efficient than scaled eigenvalues and Chebyshev. For hyperparameter learning, FITC took several hours to run, compared to minutes for the alternatives; we therefore exclude FITC from Figure 1(b). Figure 1(c) shows the time to do inference on the 691 test points, while 1(d) shows the standardized mean absolute error (SMAE) on the same test points. As expected, Lanczos and surrogate make accurate predictions much faster than Chebyshev, scaled eigenvalues, and FITC. In short, Lanczos and the surrogate approach are much faster than alternatives for hyperparameter learning with a large number of inducing points and training points.

(a) Sound data
(b) Recovery time
(c) Inference time
(d) SMAE
Figure 1: Sound modeling using 59,306 training points and 691 test points. The intensity of the time series can be seen in (a). Train time for RBF kernel hyperparameters is in (b) and the time for inference is in (c). The standardized mean absolute error (SMAE) as a function of time for an evaluation of the marginal likelihood and all derivatives is shown in (d). Surrogate is ( ——), Lanczos is ( - - -), Chebyshev is (), scaled eigenvalues is ( — + —), and FITC is ( — o —).

5.2 Daily precipitation prediction

This experiment involves precipitation data from the year of 2010 collected from around weather stations in the US111 The hourly precipitation data is preprocessed into daily data if full information of the day is available. The dataset has entries in terms of precipitation per day given the date, longitude and latitude. We randomly select data points as test points and use the remaining points for training. We then perform hyperparameter learning and prediction with the RBF kernel, using Lanczos, scaled eigenvalues, and exact methods.

For Lanczos and scaled eigenvalues, we optimize the hyperparameters on the subset of data for January 2010, with an induced grid of points per spatial dimension and in the temporal dimension. Due to memory constraints we only use a subset of entries for training with the exact method. While scaled eigenvalues can perform well when fast eigendecompositions are possible, as in this experiment, Lanczos nonetheless still runs faster and with slightly lower MSE.

Method MSE Time [min]
Lanczos 528k 3M 0.613 14.3
Scaled eigenvalues 528k 3M 0.621 15.9
Exact 12k - 0.903 11.8
Table 1: Prediction comparison for the daily precipitation data showing the number of training points , number of induced grid points , the mean squared error, and the inference time.

Incidentally, we are able to use 3 million inducing points in Lanczos and scaled eigenvalues, which is enabled by the SKI representation (Wilson and Nickisch, 2015) of covariance matrices, for a a very accurate approximation. This number of inducing points is unprecedented for typical alternatives which scale as .

5.3 Hickory data

In this experiment, we apply Lanczos to the log-Gaussian Cox process model with a Laplace approximation for the posterior distribution. We use the RBF kernel and the Poisson likelihood in our model. The scaled eigenvalue method does not apply directly to non-Gaussian likelihoods; we thus applied the scaled eigenvalue method in (Wilson and Nickisch, 2015) in conjunction with the Fiedler bound in (Flaxman et al., 2015) for the scaled eigenvalue comparison. Indeed, a key advantage of the Lanczos approach is that it can be applied whenever fast MVMs are available, which means no additional approximations such as the Fiedler bound are required for non-Gaussian likelihoods.

This dataset, which comes from the R package spatstat, is a point pattern of hickory trees in a forest in Michigan. We discretize the area into a grid and fit our model with exact, scaled eigenvalues, and Lanczos. We see in Table 2 that Lanczos recovers hyperparameters that are much closer to the exact values than the scaled eigenvalue approach. Figure 2 shows that the predictions by Lanczos are also indistinguishable from the exact computation.

Method Time [s]
Scaled eigenvalues
Table 2: Hyperparameters recovered on the Hickory dataset.
(a) Point pattern data
(b) Prediction by exact
(c) Scaled eigenvalues
(d) Lanczos
Figure 2: Predictions by exact, scaled eigenvalues, and Lanczos on the Hickory dataset.

5.4 Crime prediction

In this experiment, we apply Lanczos with the spectral mixture kernel to the crime forecasting problem considered in Flaxman et al. (2015). This dataset consists of incidents of assault in Chicago from January 1, 2004 to December 31, 2013. We use the first years for training and attempt to predict the crime rate for the last years. For the spatial dimensions, we use the log-Gaussian Cox process model, with the Matérn-5/2 kernel, the negative binomial likelihood, and the Laplace approximation for the posterior. We use a spectral mixture kernel with components and an extra constant component for the temporal dimension. We discretize the data into a spatial grid corresponding to -by- mile grid cells. In the temporal dimension we sum our data by weeks for a total of weeks. After removing the cells that are outside Chicago, we have a total of observations.

The results for Lanczos and scaled eigenvalues (in conjunction with the Fiedler bound due to the non-Gaussian likelihood) can be seen in Table 3. The Lanczos method used Hutchinson probe vectors and Lanczos steps. For both methods we allow iterations of LBFGS to recover hyperparameters and we often observe early convergence. While the RMSE for Lanczos and scaled eigenvalues happen to be close on this example, the recovered hyperparameters using scaled eigenvalues are very different than for Lanczos. For example, the scaled eigenvalue method learns a much larger than Lanczos, indicating model misspecification. In general, as the data become increasingly non-Gaussian the Fiedler bound (used for fast scaled eigenvalues on non-Gaussian likelihoods) will become increasingly misspecified, while Lanczos will be unaffected.

Method [s] [s]
Scaled eigenvalues
Table 3: Hyperparameters recovered, recovery time and RMSE for Lanczos and scaled eigenvalues on the Chicago assault data. Here and are the length scales in spatial dimensions and is the noise level. is the time for recovering hyperparameters. is the time for prediction at all observations (including training and testing).

5.5 Deep kernel learning

To handle high-dimensional datasets, we bring our methods into the deep kernel learning framework Wilson et al. (2016a) by replacing the final layer of a pre-trained deep neural network (DNN) with a GP. This experiment uses the gas sensor dataset from the UCI machine learning repository. It has instances with dimensions. We pre-train a DNN, then attach a Gaussian process with RBF kernels to the two-dimensional output of the second-to-last layer. We then further train all parameters of the resulting kernel, including the weights of the DNN, through the GP marginal likelihood. In this example, Lanczos and the scaled eigenvalue approach perform similarly well. Nonetheless, we see that Lanczos can effectively be used with SKI on a high dimensional problem to train hundreds of thousands of kernel parameters.

Method DNN Lanczos Scaled eigenvalues
Time [s]
Table 4: Prediction RMSE and per training iteration runtime.

6 Discussion

There are many cases in which fast MVMs can be achieved, but it is difficult or impossible to efficiently compute a log determinant. We have developed a framework for scalable and accurate estimates of a log determinant and its derivatives relying only on MVMs. We particularly consider scalable kernel learning, showing the promise of stochastic Lanczos estimation combined with a pre-computed surrogate model. We have shown the scalability and flexibility of our approach through experiments with kernel learning for several real-world data sets using both Gaussian and non-Gaussian likelihoods, and highly parametrized deep kernels.

Iterative MVM approaches have great promise for future exploration. We have only begun to explore their significant generality. In addition to log determinants, the methods presented here could be adapted to fast posterior sampling, diagonal estimation, matrix square roots, and many other standard operations. The proposed methods only depend on fast MVMs—and the structure necessary for fast MVMs often exists, or can be readily created. We have here made use of SKI (Wilson and Nickisch, 2015) to create such structure. But other approaches, such as stochastic variational methods (Hensman et al., 2013), could be used or combined with SKI for fast MVMs, as in (Wilson et al., 2016b). Moreover, iterative MVM methods naturally harmonize with GPU acceleration, and are therefore likely to increase in their future applicability and popularity. Finally, one could explore the ideas presented here for scalable higher order derivatives, making use of Hessian methods for greater convergence rates.


  • Kulesza et al. (2012) Alex Kulesza, Ben Taskar, et al. Determinantal point processes for machine learning. Foundations and Trends® in Machine Learning, 5(2–3):123–286, 2012.
  • MacKay (1992) David JC MacKay. Bayesian methods for adaptive models. PhD thesis, California Institute of Technology, 1992.
  • MacKay (2003) David JC MacKay. Information theory, inference and learning algorithms. Cambridge university press, 2003.
  • Rue and Held (2005) Havard Rue and Leonhard Held. Gaussian Markov random fields: theory and applications. CRC Press, 2005.
  • Rasmussen and Williams (2006) C. E. Rasmussen and C. K. I. Williams. Gaussian processes for Machine Learning. The MIT Press, 2006.
  • Boutsidis et al. (2015) Christos Boutsidis, Petros Drineas, Prabhanjan Kambadur, Eugenia-Maria Kontopoulou, and Anastasios Zouzias. A randomized algorithm for approximating the log determinant of a symmetric positive definite matrix. arXiv preprint arXiv:1503.00374, 2015.
  • Rasmussen and Ghahramani (2001) Carl Edward Rasmussen and Zoubin Ghahramani. Occam’s razor. In Neural Information Processing Systems (NIPS), 2001.
  • Quiñonero-Candela and Rasmussen (2005) Joaquin Quiñonero-Candela and Carl Edward Rasmussen. A unifying view of sparse approximate gaussian process regression. Journal of Machine Learning Research, 6(Dec):1939–1959, 2005.
  • Le et al. (2013) Q. Le, T. Sarlos, and A. Smola. Fastfood-computing Hilbert space expansions in loglinear time. In Proceedings of the 30th International Conference on Machine Learning, pages 244–252, 2013.
  • Hensman et al. (2013) J Hensman, N Fusi, and N.D. Lawrence. Gaussian processes for big data. In

    Uncertainty in Artificial Intelligence (UAI)

    . AUAI Press, 2013.
  • Wilson (2014) Andrew Gordon Wilson. Covariance kernels for fast automatic pattern discovery and extrapolation with Gaussian processes. PhD thesis, University of Cambridge, 2014.
  • Wilson et al. (2014) Andrew Gordon Wilson, Elad Gilboa, Nehorai Arye, and John P Cunningham. Fast kernel learning for multidimensional pattern extrapolation. In Advances in Neural Information Processing Systems, pages 3626–3634, 2014.
  • Wilson and Nickisch (2015) Andrew Gordon Wilson and Hannes Nickisch. Kernel interpolation for scalable structured Gaussian processes (KISS-GP). International Conference on Machine Learning (ICML), 2015.
  • Herlands et al. (2016) William Herlands, Andrew Wilson, Hannes Nickisch, Seth Flaxman, Daniel Neill, Wilbert Van Panhuis, and Eric Xing. Scalable Gaussian processes for characterizing multidimensional change surfaces. Artificial Intelligence and Statistics, 2016.
  • Snelson and Ghahramani (2006) Edward Snelson and Zoubin Ghahramani. Sparse Gaussian processes using pseudo-inputs. In Advances in neural information processing systems (NIPS), volume 18, page 1257. MIT Press, 2006.
  • Fiedler (1984) M. Fiedler. Hankel and Loewner matrices. Linear Algebra and Its Applications, 58:75–95, 1984.
  • Weyl (1912) Hermann Weyl. Das asymptotische verteilungsgesetz der eigenwerte linearer partieller differentialgleichungen (mit einer anwendung auf die theorie der hohlraumstrahlung). Mathematische Annalen, 71(4):441–479, 1912.
  • Flaxman et al. (2015) Seth Flaxman, Andrew Wilson, Daniel Neill, Hannes Nickisch, and Alex Smola. Fast kronecker inference in gaussian processes with non-gaussian likelihoods. In International Conference on Machine Learning, pages 607–616, 2015.
  • Han et al. (2015) Insu Han, Dmitry Malioutov, and Jinwoo Shin. Large-scale log-determinant computation through stochastic Chebyshev expansions. In ICML, pages 908–917, 2015.
  • (20) Shashanka Ubaru, Jie Chen, and Yousef Saad. Fast estimation of via stochastic Lanczos quadrature.
  • Bai et al. (1998) Zhaojun Bai, Mark Fahey, Gene H Golub, M Menon, and E Richter. Computing partial eigenvalue sums in electronic structure calculations. Technical report, Tech. Report SCCM-98-03, Stanford University, 1998.
  • MacKay and Gibbs (1997) D MacKay and MN Gibbs. Efficient implementation of gaussian processes. Neural Computation, 1997.
  • Stein et al. (2013) Michael L Stein, Jie Chen, Mihai Anitescu, et al. Stochastic approximation of score functions for gaussian processes. The Annals of Applied Statistics, 7(2):1162–1191, 2013.
  • Wilson et al. (2016a) Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P Xing. Deep kernel learning. In Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, pages 370–378, 2016a.
  • Rasmussen and Nickisch (2010) Carl Edward Rasmussen and Hannes Nickisch. Gaussian processes for machine learning (GPML) toolbox. Journal of Machine Learning Research (JMLR), 11:3011–3015, Nov 2010.
  • Wilson et al. (2016b) Andrew G Wilson, Zhiting Hu, Ruslan R Salakhutdinov, and Eric P Xing. Stochastic variational deep kernel learning. In Advances in Neural Information Processing Systems, pages 2586–2594, 2016b.
  • Silverman (1985) Bernhard W Silverman. Some aspects of the spline smoothing approach to non-parametric regression curve fitting. Journal of the Royal Statistical Society. Series B (Methodological), pages 1–52, 1985.
  • Quinonero-Candela et al. (2007) Joaquin Quinonero-Candela, Carl Edward Rasmussen, and Christopher KI Williams. Approximation methods for Gaussian process regression. Large-scale kernel machines, pages 203–223, 2007.
  • Higham (2008) Nicholas J Higham. Functions of matrices: theory and computation. SIAM, 2008.
  • Hutchinson (1990) Michael F Hutchinson. A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Communications in Statistics-Simulation and Computation, 19(2):433–450, 1990.
  • Gil et al. (2007) Amparo Gil, Javier Segura, and Nico Temme. Numerical Methods for Special Functions. SIAM, 2007.
  • Golub and Meurant (2010) Gene Golub and Gérard Meurant. Matrices, Moments and Quadrature with Applications. Princeton University Press, 2010.
  • Cullum and Willoughby (2002) Jane K Cullum and Ralph A Willoughby. Lanczos algorithms for large symmetric eigenvalue computations: Vol. I: Theory. SIAM, 2002.
  • Saad (1992) Youcef Saad. Numerical methods for large eigenvalue problems. Manchester University Press, 1992.
  • Buhmann (2000) Martin Dietrich Buhmann. Radial basis functions. Acta Numerica 2000, 9:1–38, 2000.
  • Fasshauer (2007) Gregory E Fasshauer. Meshfree approximation methods with MATLAB, volume 6. World Scientific, 2007.
  • Schaback and Wendland (2006) Robert Schaback and Holger Wendland. Kernel techniques: from machine learning to meshless methods. Acta Numerica, 15:543–639, 2006.
  • Wendland (2004) Holger Wendland. Scattered data approximation, volume 17. Cambridge university press, 2004.
  • Avron and Toledo (2011) Haim Avron and Sivan Toledo. Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix. J. ACM, 58(2):8:1–8:34, 2011. doi: 10.1145/1944345.1944349. URL
  • Wathen and Zhu (2015) Andrew J. Wathen and Shengxin Zhu. On spectral distribution of kernel matrices related to radial basis functions. Numer. Algor., 70:709–726, 2015.

Appendix A Background

Two popular covariance kernels are the RBF kernel

and the Matérn kernel

where , , and are popular choices for to model heavy-tailed correlations between function values. The spectral behavior of these and other kernels has been well-studied for years, and we recommend [Wathen and Zhu, 2015] for recent results. Particularly relevant to our discussion is a theorem due to Weyl, which says that if a symmetric kernel has continuous derivatives, then the eigenvalues of the associated integral operator decay like . Hence, the eigenvalues of kernel matrices for the smooth RBF kernel (and of any given covariance matrix based on that kernel) tend to decay much more rapidly than those of the less smooth Matérn kernel, which has two derivatives at zero for , one derivative at zero for , and no derivatives at zero for . This matters to the relative performance of Chebyshev and Lanczos approximations of the log determinant for large values of and small values of on the exact and approximate RBF kernel.

Appendix B Methods

b.1 Scaled eigenvalue method

The scaled eigenvalue method was introduced in  [Wilson et al., 2014] to estimate , where consists of points. The eigenvalues of can be approximated using the largest eigenvalues of a covariance matrix on a full grid with points such that . Specifically,

The induced kernel plays the role of when the scaled eigenvalue method is applied to SKI and the eigenvalues of can be efficiently computed. Assuming that the eigenvalues can be computed efficiently is a much stronger assumption than our fast MVM based approach.

b.2 Radial basis function surrogates

Radial basis function (RBF) interpolation is one of the most popular approaches to approximating scattered data in a general number of dimensions [Buhmann, 2000, Fasshauer, 2007, Schaback and Wendland, 2006, Wendland, 2004]. Given distinct interpolation points , the RBF model takes the form


where the kernel is a one-dimensional function and , the space of polynomials with variables of degree no more than . There are many possible choices for such as the cubic kernel and the thin-plate spline kernel . The coefficients are determined by imposing the interpolation conditions for and the discrete orthogonality condition


For appropriate RBF kernels, this linear system is nonsingular provided that polynomials in are uniquely determined by their values on the interpolation set.

b.3 Comparison to a reference kernel

Suppose more generally that is an approximation to a reference kernel matrix , and let . Let and be the log likelihood functions for the two kernels; then

If we are willing to pay the price of a few MVMs with , we can use these expressions to improve our maximum likelihood estimate. Let and be independent probe vectors with and . To estimate the trace in the derivative computation, we use the standard stochastic trace estimation approach together with the observation that :

This linearization may be used directly (with a stochastic estimator); alternately, if we have an estimates for and , we can substitute these in order to get estimated bounds on the magnitude of the derivatives. Coupled with a similar estimator for the Hessian of the likelihood function (described in the supplementary materials), we can use this method to compute the maximum likelihood parameters for the fast kernel, then compute a correction to estimate the maximum likelihood parameters of the reference kernel.

Appendix C Additional experiments

This section contains several experiments with synthetic data sets to illustrate particular aspects of the method.

c.1 1D cross-section plots

In this experiment we compare the accuracy of Lanczos and Chebyshev for 1-dimensional perturbations of a set of true hyper-parameters, and demonstrate how critical it is to use diagonal replacement for some approximate kernels. We choose the true hyper-parameters to be and consider two different types of datasets. The first dataset consists of equally spaced points in the interval in which case the kernel matrix of a stationary kernel is Toeplitz and we can make use of fast matrix-vector multiplication. The second dataset consists of data points drawn independently from a distribution. We use SKI with cubic interpolation to construct an approximate kernel based on equally spaced points. The function values are drawn from a GP with the true hyper-parameters, for both the true and approximate kernel. We use iterations for Lanczos and Chebyshev moments in order to assure convergence of both methods. The results for the first dataset with the RBF and Matérn kernels can be seen in Figure 3(a)-3(d). The results for the second dataset with the SKI kernel can be seen in Figure 4(a)-4(d).

(a) log marginal likelihood for the RBF kernel
(b) log marginal likelihood for the Matérn kernel
(c) log determinant for the RBF kernel
(d) log determinant for the Matérn kernel
Figure 3: 1-dimensional perturbations for the exact RBF and Matérn 1/2 kernel where the data is equally spaced points in the interval . The exact values are (), Lanczos is (—–), Chebyshev is (—–). The error bars of Lanczos and Chebyshev are standard deviation and were computed from runs with different probe vectors
(a) log marginal likelihood for the RBF kernel
(b) log marginal likelihood for the Matérn kernel
(c) log determinant for the RBF kernel
(d) log determinant for the Matérn kernel
Figure 4: 1-dimensional perturbations with the SKI (cubic) approximations of the RBF and Matérn 1/2 kernel where the data is 1000 points drawn from . The exact values are (), Lanczos with diagonal replacement is (—–), Chebyshev with diagonal replacement is (—–), Lanczos without diagonal replacement is (—–), Chebyshev without diagonal replacement is (—–), and scaled eigenvalues is ( ). Diagonal replacement makes no perceptual difference for the RBF kernel so the lines are overlapping in this case. The error bars of Lanczos and Chebyshev are standard deviation and were computed from runs with different probe vectors

Lanczos yields an excellent approximation to the log determinant and its derivatives for both the exact and the approximate kernels, while Chebyshev struggles with large values of and small values of on the exact and approximate RBF kernel. This is expected since Chebyshev has issues with the singularity at zero while Lanczos has large quadrature weights close to zero to compensate for this singularity. The scaled eigenvalue method has issues with the approximate Matérn 1/2 kernel.

c.2 Why Lanczos is better than Chebyshev

In this experiment, we study the performance advantage of Lanczos over Chebyshev. Figure 5 shows that the Ritz values of Lanczos quickly converge to the spectrum of the RBF kernel thanks to the absence of interior eigenvalues. The Chebyshev approximation shows the expected equioscillation behavior. More importantly, the Chebyshev approximation for logarithms has its greatest error near zero where the majority of the eigenvalues are, and those also have the heaviest weight in the log determinant.

Another advantage of Lanczos is that it requires minimal knowledge of the spectrum, while Chebyshev needs the extremal eigenvalues for rescaling. In addition, with Lanczos we can get the derivatives with only one MVM per hyper-parameter, while Chebyshev requires an MVM at each iteration, leading to extra computation and memory usage.

(a) True spectrum
(b) Lanczos weights
(c) Chebyshev weights
(d) Chebyshev absolute error
Figure 5: A comparison between the true spectrum, the Lanczos weights (, and the Chebyshev weights () for the RBF kernel with , , and . All weights and counts are on a log-scale so that they are easier to compare. Blue bars correspond to positive weights while red bars correspond to negative weights.

c.3 The importance of diagonal correction

This experiment shows that diagonal correction of the approximate kernel can be very important. Diagonal correction cannot be used efficiently for some methods, such as the scaled eigenvalue method, and this may hurt its predictive performance. Our experiment is similar to [Quiñonero-Candela and Rasmussen, 2005]. We generate uniformly distributed points in the interval , and we choose a small number of inducing points in such a way that there is a large chunk of the interval where there is no inducing point. We are interested in the behavior of the predictive uncertainties on this subinterval. The function values are given by

and normally distributed noise with standard deviation

is added to the function values. We find the optimal hyper-parameters of the Matérn using the exact method and use these hyper-parameters to make predictions with Lanczos, Chebyshev, FITC, and the scaled eigenvalue method. We consider Lanczos both with and without diagonal correction in order to see how this affects the predictions. The results can be seen in Figure 6.

(a) Lanczos with diagonal correction
(b) Lanczos without diagonal correction
(c) Chebyshev with diagonal correction
(d) Chebyshev without diagonal correction
(e) FITC
(f) Scaled eigenvalue method
Figure 6: Example that shows how important diagonal correction can be for some kernels. The Matérn kernel was used to fit the data given by the black dots. This data was generated from the function to which we added normally distributed noise with standard deviation

. We used the exact method to find the optimal hyper-parameters and used these hyper-parameters to study the different behavior of the predictive uncertainties when the inducing points are given by the green crosses. The solid blue line is the predictive mean and the dotted red lines shows a confidence interval of two standard deviations.

It is clear that Lanczos and Chebyshev are too confident in the predictive mean when diagonal correction is not used, while the predictive uncertainties agree well with FITC when diagonal correction is used. The scaled eigenvalue method cannot be used efficiently with diagonal correction and we see that this leads to predictions similar to Lanczos and Chebyshev without diagonal correction. The flexibility of being able to use diagonal correction with Lanczos and Chebyshev makes these approaches very appealing.

c.4 Surrogate log determinant approximation

The point of this experiment is to illustrate how accurate the level-curves of the surrogate model are compared to the level-curves of the true log determinant. We consider the RBF and the Matérn kernels and the same datasets that we considered in C.1. We fix and study how the level curves compare when we vary and . Building the surrogate with all three hyper-parameters produces similar results, but requires more design points. We use design points to construct a cubic RBF with a linear tail. The values of the log determinant and its derivatives are computed with Lanczos. It is clear from Figure 7 that the surrogate model does a good job approximating the log determinant for both kernels.

(a) RBF exact
(b) Matérn exact
(c) RBF surrogate
(d) Matérn surrogate
Figure 7: Level curves of the exact and surrogate approximation of the log determinant as a function of and for the RBF and Matérn kernels. We used and the dataset consisted of 1000 equally spaced points in the interval . The surrogate model was constructed from the points shown with () and the log determinant values were computed using stochastic Lanczos.

c.5 Kernel hyper-parameter recovery

This experiments tests how well we can recover hyper-parameters from data generated from a GP. We compare Chebyshev, Lanczos, the surrogate, the scaled eigenvalue method, and FITC. We consider a dataset of points generated from a distribution. We use SKI with cubic interpolation and a total of inducing points for Lanczos, Chebyshev, and then scaled eigenvalue method. FITC was used with equally spaced points because it has a longer runtime as a function of the number of inducing points. We consider the RBF kernel and the Matérn kernel and sample from a GP with ground truth parameters . The GPs for which we try to recover the hyper-parameters were generated from the original kernel. It is important to emphasize that there are two sources of errors present: the error from the kernel approximation errors and the stochastic error from Lanczos and Chebyshev. We saw in Figure 3 and 4 that the stochastic error for Lanczos is relatively small, so this follow-up experiment helps us understand how Lanczos is influenced by the error incurred from an approximate kernel. We show the true log marginal likelihood, the recovered hyper-parameters, and the run-time in Table 5.

max width= RBF Matérn True Hypers Exact Hypers Time (s) Lanczos Hypers Time (s) Chebyshev Hypers Time (s) Surrogate Hypers Time (s) Scaled eigenvalues Hypers Time (s) FITC Hypers Time (s)

Table 5: Hyper-parameter recovery for the RBF and Matérn kernels. The data was generated from normally distributed points. Lanczos, surrogate, and scaled eigenvalues all used 2000 inducing points while FITC used 750. These numbers where chosen to make their run times close to equal. Diagonal correction was applied to the Matérn approximate kernel. The value of the log marginal likelihood was was computed from the exact kernel and shows the value of the hyper-parameters recovered by each method. We ran Lanczos 5 times and averaged the values.

It is clear from Table 5 that most methods are able to recover parameters close to the ground truth for the RBF kernel. The results are more interesting for the Matérn kernel where FITC struggles and the parameters recovered by FITC have a value of the log marginal likelihood that is much worse than the other methods.