A Sparse Covariance Function for Exact Gaussian Process Inference in Large Datasets

1 Introduction

Gaussian processes (GPs) are a useful and powerful tool for regression in supervised machine learning [Rasmussen and Williams, 2006]. The range of applications includes geophysics, mining, hydrology, reservoir engineering and robotics. Despite its increasing popularity, modelling largescale spatial stochastic processes is still challenging. The difficulty comes from the fact that inference in GPs is usually computationally expensive due to the need to invert a potentially large covariance matrix during inference time, which has O (N 3 ) cost. For problems with thousands of observations, exact inference in normal GPs is intractable and approximation algorithms are required.

Most of the approximation algorithms employ a subset of points to approximate the posterior distribution at a new point given the training data and hyper-parameters. These approximations rely on heuristics to select the subset of points [Lawrence et al., 2003; Seeger et al., 2003], or use pseudo targets obtained during the optimisation of the log-marginal likelihood of the model [Snelson and Ghahramani, 2006].

In this work, we address the complexity problem differently. Instead of relying on sparse GP approximations, we

Fabio Ramos

Australian Centre for Field Robotics The University of Sydney NSW 2006, Australia f.ramos@acfr.usyd.edu.au

propose a new covariance function which provides intrinsically sparse covariance matrices. This allows exact inference in GPs using conventional methods. As the new sparse covariance function can be multiplied by any other valid covariance function and the result is a sparse covariance matrix, a lot of flexibility is given to practitioners to accurately model their problems while still preserving sparseness properties. We call the GPs constructed using our sparse covariance function Exact Sparse Gaussian Processes (ESGPs). The main idea behind is the formulation of a valid and smooth covariance function whose output equals to zero whenever the distance between input observations is larger than a hyperparameter. As with other hyper-parameters, this can be estimated by maximising the marginal likelihood to better model the properties of the data such as smoothness, characteristic length-scale, and noise. Additionally, the proposed covariance function in much resembles the popular squared exponential in terms of smoothness, being four times continuously differentiable. We empirically compare ESGP with local approximation techniques and demonstrate how other covariance functions can be integrated in the same framework. Our method results in very sparse covariance matrices (up to 90% of the elements are zeros in in-ground grade estimation problems) which requires significantly less memory while providing similar performance.

This paper is organised as follows. In Section 2 we review the basics of GP regression and introduce notation. Section 3 summarises previous work on approximate inference with GPs. Section 4 presents our new intrinsically sparse covariance function and its main properties. We evaluate the framework providing experimental results in both artificial and real data in Section 6. Finally, Section 7 concludes the paper and discusses further developments.

2 Gaussian Processes

In this section we briefly review Gaussian Processes for regression and introduce notation. We consider the supervised learning problem where given a training set D = {xi, yi}Ni=1 consisting of N input points xi ∈ RD and the corresponding outputs yi ∈ R the objective is to compute the predictive distribution f (x∗) at a new test point x∗. A Gaussian process model places a multivariate Gaussian distribution over the space of function variables f(x) mapping input to output spaces. The model is specified by defining a mean function m(x) and the covariance function k(x, x′ ) resulting in the Gaussian process written as

Denoting groups of these points as (X, f, y) = ({xi}, {fi}, {yi})Ni=1 for the training set and (X∗, f∗, y∗) = ({x∗,i}, {f∗,i}, {y∗,i})Ni=1 for the testing points, the joint Gaussian distribution with m(x) = 0 is

where N (μ,Σ) is a multivariate Gaussian distribution with mean μ and covariance Σ, and K is used to denote the covariance matrix computed between all points in the set. If we assume observations with Gaussian noise ϵ and variance σ2 such that y = f(x) + ϵ, the joint distribution becomes

A popular choice for the covariance function is the squared exponential used in this paper for comparisons in the experiment section:

with M = diag(l)-2 where l is a vector of positive numbers representing the length-scales in each dimension.

2.1 Inference for New Points

By conditioning on the observed training points, the predictive distribution can be obtained as

where

From Equation 5, it can be observed that the predictive mean is a linear combination of N kernel functions each centred on a training point, μ∗ = ∑Ni=1 αik (xi, x∗), where α = (K (X, X) + σ2 I )-1 y. A GP is also a best unbiased linear estimator [Cressie, 1993; Kitanidis, 1997] in the mean squared error sense. During inference, most of the computational cost takes place while computing the inversion in Equation 5, which is O (N 3 ) if implemented naïvely.

2.2 Learning Hyper-Parameters

Commonly, the covariance function k(x, x′ ) is parametrised by set of hyper-parameters θ, and we can write k (x, x′ ; θ). These parameters allow for more flexibility in modelling the properties of the data. Thus, learning a GP model is equivalent to determining the hyper-parameters of the covariance function from some training dataset. In a Bayesian framework this can be performed by maximising the log of the marginal likelihood w.r.t. θ:

where Ky = K(X,X) + σ2 I is the covariance matrix for the targets y. The marginal likelihood has three terms (from left to right), the first accounts for the data fit; the second is a complexity penalty term (encoding the Occam's Razor principle) and the last is a normalisation constant.

Eq. (6) is a non-convex function on the hyper-parameters θ and therefore only local maxima can be obtained. In practice, this is not a major issue since good local maxima can be obtained with gradient descent techniques by using multiple starting points. However, this requires the computation of partial derivatives resulting in:

Note that this expression requires the computation of partial derivatives of the covariance function w.r.t θ.

3 RelatedWork

Recently, there have been several methods proposed to tackle the problem of GP inference in large datasets. However, most of these approaches rely on approximation techniques. A common and simple procedure is to select a subset of data points and perform inference using only these points. This is equivalent to ignoring part of the data which makes the selection of the subset very important. In [Lawrence et al., 2003], the selection of points is based on the differential entropy. Similarly, [Seeger et al., 2003] suggest the use of another information theory quantity, the information gain.

Another interesting procedure is to select a subset of data points to act as an inducing set and project this set up to all the data points available. This is known as sparse GP approximation [Williams and Seeger, 2001; Smola and Bartlett, 2001; Candela and Rasmussen, 2005] which usually performs better than simply selecting a subset of data points. However, the definition of the inducing set is difficult and can involve non-convex optimisations [Snelson and Ghahramani, 2006].

Local methods have been applied in geostatistics for a long time [Wackernagel, 2003]. The idea is to perform inference by evaluating the covariance function only at points in the neighbourhood of a query point. This method can be effective but the definition of the neighbourhood is crucial. Our method is inspired by this approach but rather than defining the neighbourhood manually, we obtain it automatically during learning. An interesting idea on combining local and global methods (such as the sparse Gaussian process) was proposed in [Snelson and Ghahramani, 2007]. We compare our method to theirs in Section 6.

This work differs from other studies by not addressing the GP inference problem through an approximation technique. Rather, it proposes a new covariance function that naturally generates sparse covariance matrices. This idea was used in [Wendland, 2005] with piecewise polynomials but extensions to multiple dimensions is difficult due to the need to guarantee positive definiteness. A similar formulation to ours was proposed in [Storkey, 1999]. However, there is no hyperparameter learning and the main properties are not analysed.

To the best of our knowledge, this work is the first to demonstrate with real examples how the complexity problem can be addressed through the construction of a new sparse covariance function allowing for exact GP inference in large datasets.

4 Exactly Sparse Gaussian Processes

For very large datasets, the inversion or even storage of a full matrix K(X, X) + σ2 I can be prohibitive. In geology problems for example, it is not uncommon to have datasets with 100K points or more. To deal with such large problems while still being able to perform exact inference in the GP model, we develop the covariance function below. First, note that the mean prediction in Eq. 5 can be rewritten as a linear combination of N evaluations of the covariance function, each one centred on a training point, μ∗ = ∑Ni=1 αik (x∗, xi), where α = (K (X, X) + σ2 I )-1 y. To avoid the inversion of the full matrix, we can instead develop a covariance function whose output vanishes outside some region R, so that k(x∗, xi) = 0 when xi is outside a region R. In this way, only a subset of α would need to be computed which effectively means that only few columns of (K (X, X) + σ2 I )-1 need to be computed, significantly reducing the computational and storage costs as R diminishes. As we shall see, the region can be specified automatically during learning.

4.1 Intrinsically Sparse Covariance Function

The covariance function we are looking for must vanish out of some finite region R for exact sparse GP regression. It must produce smooth curves but it should not be infinitely differentiable so that it can be applicable to problems with some discontinuities. For our derivation, the function g (x) = cos2 (πx)H (0.5- |x|) was chosen, which due to cos2 (πx) = (cos (2πx) + 1) /2 is actually the cosine function shifted up, normalised and set to zero out of the interval x ∈ (-0.5, 0.5). The cosine function was selected as the basis function due to the following reasons: 1) it is analytically well tractable; 2) integrals with finite limits containing combinations of its basic form can be calculated in closed form; and 3) the cosine function usually provides good approximations for different functions, being the core element for Fourier analysis. Here and afterwards H (·) represents the Heaviside unit step function. As it stands, the chosen basis function g (x) is smooth on the whole real axis, vanishes out of the interval x ∈ (-0.5, 0.5) and has discontinuities in the second derivative. To derive a valid covariance function we conduct calculations analogous to presented in [Rasmussen and Williams, 2006]. Using the transfer function h (x; u) = g (x - u) the following 1D covariance function is obtained:

x-x'

Figure 1: Plot showing the output of the covariance function for different values of Δx.

Due to the chosen form of the basis functions, the integral in Eq. (8) can be analytically evaluated (see Appendix A for details) to result in:

where σ0 > 0 is a constant coefficient, l > 0 is a given scale and d is the distance between the points:

From Eq. (8) it follows that for any points xi and any real numbers ai where i = 1, 2, ... , n the inequality

holds, so that the constructed covariance function is positive semi-definite.

Based on Eqs. (9)-(10) we calculate that

which shows that the covariance function k1 is continuous and has continuous 4th derivative at d = l.

The function k1 (Δx, 0; l, 1) is compared with squared exponential in Figure 1. Note that it follows the squared exponential covariance function closely but vanishes when |Δx| ≥ l.

4.2 Extending toMultiple Dimensions

This covariance function can be extended to multiple dimensions in the following ways:

1. Using direct products for all axes:

where D is the dimensionality of the points and l is the vector of the characteristic lengths, l = (l1, l2, ..., lD)T . 2. Using Mahalanobis distance:

where σ0 > 0, Ω is positive semi-definite and

After this point we will frequently use the short notation kM for the function k(2) M (r; σ0,Ω).

4.3 Important properties of the new covariance function

The developed multi-dimensional covariance function in both forms k(1) M and k(2) M has the following remarkable properties:

1. It vanishes out of some finite region R:

Sizes of the regions R1 and R2 can be controlled via the characteristic lengths li. Moreover, these sizes can be learnt from data by maximising the marginal likelihood as common in the GP framework.

1. All the derivatives up to (and including) the fourth order derivative are continuous, which guarantees mean square differentiability up to the corresponding order of the sample curves in GPs. There are discontinuities for the fifth order gradient.
2. The region R1 is a D dimensional rectangle and the region R2 is a D dimensional ellipsoid.
3. In the case of 1 dimension k(1) M and k(2) M become identical.
4. The covariance function is anisotropic, i.e. has different inner properties for different dimensions.
5. This covariance function leads to sparse covariance matrices and allows GP inference in large datasets without the need for approximations.

5 Partial Derivatives for Learning

Learning the GP requires the computation of the covariance function partial derivatives w.r.t. the hyper-parameters (Eq. (7)). Based on Eqs.(9), (12), the following expressions for the partial derivatives of k(1) M (x, x′ ; l, σ0) can be calculated:

where i = 1, 2, ..., D.

For the second case, if Ω is diagonal and positive definite, it can be expressed via the characteristic lengths as follows:

From Eqs. (14), (19) it follows that

Based on Eqs. (13), (19)-(20) the following gradient components of this multi-dimensional covariance function kM can be obtained:

In Eq. (22), r is in the denominator, so that direct calculations cannot be carried out using Eq. (22) when r = 0. However, using the equality

one can directly show that

Based on Eq. (25) it must be taken directly

Eqs. (21)-(23), (26) fully define the gradient of the new covariance function kM at every point and can be directly used in the learning procedure.

6 Experiments

This section provides empirical comparisons between exact sparse GP, conventional GP with squared exponential and approximation procedures.

6.1 Artificial Dataset

In this experiment we compare the exact GP with the proposed covariance function against the approach proposed in [Snelson and Ghahramani, 2007]. The data is essentially the same as presented in the experiment section in [Snelson and Ghahramani, 2007]. As can be observed in Figure 2, the sparse covariance function provides a much smoother prediction for the underlying function than the combination of local and global approximations. This example shows qualitatively that in some situations approximation methods can lead to discontinuities. The same does not occur in the exact sparse GP inference.

6.2 Rainfall Dataset

In this experiment we compare the exact sparse GP with the exact GP with the squared exponential covariance function and the covariance function obtained by the multiplication of both of them. The dataset used is a popular dataset in geostatistics for comparing inference procedures and is known as the Spatial Interpolation Comparison dataset [Dubois et al., 2003](SIC) 1 . The dataset consists of 467 points measuring rainfall in 2D space. We divide these points into two sets, inference and testing. The inference set contains the points used to perform inference on the testing points. For each case the experiment is repeated 1500 times with randomly selected inference and testing sets. Figure 3 shows the normalised squared error for the different covariance functions and the standard deviation (one sigma for each part of the bar) as a function of the number of inference points. As the number of inference points increases, so does the size of the covariance matrix. The results demonstrate that very similar errors are obtained for the different covariances. However, the sparse covariance function produces sparse matrices thus requiring much less floating point operations. Figure 4 shows the percentage of zeros in the covariance matrix as a function of the number of inference points. As can be observed, the percentage of zeros grows quickly as more inference points are added and it reaches its limit around 100 inference points. Although the percentage of zeros reaches its limit, the number of zeros in the covariance matrix continues to increase because the size of the covariance matrix increases with the number of inference points. Also worth noticing is the performance of the multiplication between the two covariance functions. The error is essentially the same as for the sparse covariance function alone but the percentage of zeros is significantly smaller. This example demonstrates the benefits of the proposed approach in reducing storage and number of operations for similar accuracy.

1 The SIC dataset can be downloaded at: http://www.aigeostats.org

6.3 Iron Ore Dataset

In this dataset the goal is to estimate iron ore grade in 3D space over a region of 2.5 cubic kilometres. The dataset is from an open pit iron mine in Western Australia. About 17K samples were collected and the iron concentration measured with X-Ray systems. We divide the 17K dataset points into inference and testing sets. The inference set is taken arbitrarily from the dataset and from the remaining points the testing points are arbitrarily chosen. The experiments are repeated 500 times. Figure 5 shows the normalised mean squared error and the standard deviation (one sigma for each part of the bar) in the cases of squared exponential, sparse covariance functions and their product. The results demonstrate that all the three lead to similar errors with the sparse covariance function performing slightly better with the increase of the number of inference points. Figure 6 shows that although they result in similar errors, the sparse and the product lead to about 48% of zeros in the covariance matrix, which is 120K to 12M cells exactly equal to zero when the number of inference points varies from 500 to 5000. This example demonstrates that the proposed method provides greater savings for bigger datasets.

6.4 Speed Comparison

This experiment demonstrates the computational gains in using the proposed sparse covariance function. Synthetic datasets containing 1000, 2000, 3000 and 4000 points were generated by sampling a polynomial function with white noise. We compare the speed of a GP with the sparse covariance function to a GP with the squared exponential covariance function for different length scales and corresponding number of non-zero elements in the covariance matrix. The results are presented in Figure 7. The code is implemented in Matlab and uses the sparse matrix implementation package provided. Further gains could be obtained in more efficient sparse matrix packages. As the number of points in the datasets increases, the speed up becomes more evident. With 4000 points, the sparse covariance function is faster than the squared exponential for up to 70% of non-zeros elements in the covariance matrix. After this point, the computational cost of the sparse matrix implementation becomes dominant. As in general the sparse covariance function provides covariance matrices much sparser, speed gains can be quite substantial (in addition to storage gains).

7 Conclusions

This paper proposed a new covariance function constructed upon the cosine function for analytical tractability that naturally provides sparse covariance matrices. The sparseness of the data is controlled by a hyper-parameter that can be learnt from data. The sparse covariance function enables exact inference in GPs even for large datasets, providing both storage and computational benefits. Although the main focus of this paper was on GPs, it is important to emphasise that the covariance function proposed is also a Mercer kernel and therefore can be applied to kernel machines such as support vector machines, kernel principal component analysis and others [Schölkopf and Smola, 2002]. The use of the sparse covariance function in other kernel methods is objective of our future work.

Acknowledgements

We thank Edward Snelson for providing the pictures for the comparison with local and global sparse Gaussian process approximation. This work has been supported by the Rio Tinto Centre for Mine Automation and the ARC Centre of Excellence programme, funded by the Australian Research Council (ARC) and the New South Wales State Government.