Multi-Kernel Gaussian Processes

1 Introduction

Over the past years Gaussian processes (GPs) have become an important tool for machine learning. Initially proposed under the name kriging in the geostatistical literature [Cressie, 1993], its formulation as a non-parametric Bayesian regression technique boosted the application of these models to problems beyond spatial stochastic process modeling [MacKay, 1997; Rasmussen and Williams, 2006].

Gaussian process inference is usually formulated for a single output, i.e. given a training set consisting of inputs x and outputs y compute the mean and variance of the predictive distribution for a new point x∗ . However, in many machine learning problems the objective is to infer multiple tasks jointly, possibly using the dependencies between them to improve results. Real world examples of this problem include ore mining where the objective is to infer the concentration of several chemical components to assess the ore quality. Similarly, in robotics and control problems there are several actuators and the understanding and accurate modeling of the dependencies between the control outputs can significantly improve the controller.

Given its importance for mining applications, the geostatistical community has studied the multiple output case for several years under the name co-kriging [Wackernagel, 2003].

Fabio Ramos

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

The main difficulty in cokriging or multiple output GPs is the definition of a valid covariance function, able to model dependencies between different outputs in the cross terms. Recently, new approaches were proposed where the cross terms are obtained from a convolution process [Boyle and Frean, 2005; Alvarez and Lawrence, 2009]. These approaches use the same covariance function to model different tasks. In many problems however, tasks can have very different behaviors and still be dependent. For example, in mining a chemical component might have a very variable concentration but still be correlated to another component with a smooth concentration in the same area.

In this paper we generalize the multi-task Gaussian process through convolutional processes to allow the use of multiple covariance functions, possibly having a different covariance function per task. We develop a general mathematical framework to build valid cross covariance terms and demonstrate the applicability to real world problems. As examples, we provide closed form solutions to cross covariance terms between Mat´ern and squared exponential, Mat´ern and sparse, and sparse and squared exponential covariance functions. The sparse covariance function was recently proposed for exact GP inference in large datasets [Melkumyan and Ramos, 2009]. This property can be naturally incorporated in the multiple output case with the definition of valid cross sparse terms as described in this paper.

The paper is organized as follows: in Section 2 we discuss the previous work. Section 3 reviews multiple output Gaussian process inference and learning. In Section 4 we propose a new methodology for constructing multi-task covariance functions. Section 5 presents analytical solutions for six cross covariance terms and Section 6 demonstrates experimental results. Finally, Section 7 concludes the paper.

2 RelatedWork

Multi-task learning has received a lot of attention recently. In regression problems, this can be achieved in the linear case with multidimensional linear regression [Hastie et al., 2001]. For non-linear cases, neural networks were employed in [Caruana, 1997] where hidden units represent the sharing knowledge between multiple tasks. In [Evgeniou et al., 2005] correlation between multiple tasks is obtained by specifying a correlated prior over linear regression parameters for support vector machines. The kernels obtained can model linear dependencies aptly but are not suitable for more complex (non-linear) dependencies between tasks.

An interesting multi-task formulation combining Gaussian processes and parametric models was proposed in [Teh et al., 2005]. Each task is modeled with a different covariance function and correlated with a positive definite matrix learned by optimizing a variational lower bound on the marginal likelihood.

In geostatistics the multi-task inference problem has been investigated under the name of co-kriging [Wackernagel, 2003]. Different forms of co-kriging have been proposed. Ordinary co-kriging assumes zero mean while simple co-kriging computes an arbitrary mean in the same inference process. In both cases the covariance function is usually assumed to have the same form and parametrization for the different tasks. The linear model of coregionalization proposes the introduction of positive definite matrix multiplying the covariance function similar to the prior model proposed in [Bonilla et al., 2008].

A convolution process between a smoothing kernel and a latent function was used by Boyle and Frean to specify cross covariance terms in [Boyle and Frean, 2005]. Alvarez and Lawrence provide a sparse approximation for that approach in [Alvarez and Lawrence, 2009]. Our model can also be seen as a convolution process of two smoothing kernels (basis functions) assuming the influence of one latent function. Extensions to multiple latent functions are also possible by following the procedure in [Alvarez and Lawrence, 2009].

Our method differs from previous approaches by directly modeling the dependencies between multiple tasks through new cross covariance terms. A mathematical procedure is presented to obtain basis functions for general stationary covariance functions which can be used to construct new cross covariance terms. This is expected to provide much more flexibility in representing complex dependencies between the different outputs. It also generalizes the use of kernels with different forms for the joint multiple output prediction.

3 Multiple Output Gaussian Processes

Consider the supervised learning problem of estimating M tasks y∗ for a query point x∗ given a set X of inputs x11, . . . , xN11, x12, . . . , xN22, . . . , x1M, . . . , xNMM and corresponding noisy outputs y = (y11, . . . , yN11, y12, . . . , yN22, . . . , y1M, . . . , yNMM)T , where xil and yil correspond to the ith input and output for task l respectively, and Nl is the number of training examples for task l.

The Gaussian processes approach to this problem is to place a Gaussian prior over the latent functions fl mapping inputs to outputs. Assuming zero mean for the outputs we define a covariance matrix over all latent functions in order to explore the dependencies between different tasks

where Klk with l, k = 1 : M define the positive semi-definite (PSD) block matrix K. In this work we allow Klk to be computed with multiple covariance functions (or kernels) resulting in a final PSD matrix. To fully define the model we need to specify the auto covariance terms klk with l = k and the cross covariance terms klk with l ̸=k. The main difficulty in this problem is to define cross covariance terms that provide PSD matrices. A general framework for this is proposed in the next section.

With these definitions, inference can be computed using the conventional GP equations for mean and variance:

¯ fl(x∗ ) = kTl K-1 y y, V[fl(x∗ )] = kl∗ -kTl K-1 y kl, (2) where Ky = K+σ2 I is the covariance matrix for the targets y and

Similarly, learning can be performed by maximizing the log marginal likelihood

where Θ is a set of hyper-parameters.

4 ConstructingMulti-Kernel Covariance Functions

4.1 General Form

To construct valid cross covariance terms between M covariance functions k11, k22, . . . , kMM we need to go back to their basis functions and construct cross covariance terms between all the kernel pairs. The proposition below states that if the M covariance functions kii, i = 1 : M can be written as convolution of their basis functions (or smoothing kernels) gi, defining the cross covariance terms as the convolutions of gi with gj where i, j = 1 : M results in a PSD multi-task covariance function.

Proposition 1. Suppose kii (x, x′ ), i = 1 : M are singletask stationary covariance functions and can be written in the following form:

Then the M task covariance function defined as

where xi and x′j belong to the tasks i and j, respectively, is a PSD multi-task covariance function.

This preposition can be proved by considering the quadratic form that it generates. After some algebraic manipulations this quadratic form can be rearranged into a sum of squares which proves that K is a M task PSD covariance function. The details of the proof can be found in the technical report ACFR-TR-2011-002 1 . The covariance functions kii, i = 1 : M can have the same form with different hyperparameters or can have completely different forms. When the covariance functions can be written as in Eq. (4) the cross covariance terms can be calculated as in Eq. (5). The main difficulty is finding gi (smoothing kernel) for popular covariance functions and computing the integrals in Eq. (5). The following section demonstrates how to obtain smoothing kernels for stationary covariance functions through the Fourier analysis.

1 The technical report ACFR-TR-2011-002 is available at: http://www.acfr.usyd.edu.au/techreports/

4.2 Constructing Cross Covariance Terms with Fourier Analysis

Consider the covariance function k (x, x′ ) which can be represented in the form

where g (u) ≡ g (-u):

Changing the variable of integration in Eq. (6) we obtain

where ∗ stands for convolution and τ = x - x′ .

Applying the Fourier transformation

to Eq. (7) and using the equality

one has that

Using Eqs. (9), (8) one can calculate the basis function g (τ) of arbitrary stationary covariance function k (τ) via the formula:

5 Examples

In this section we provide analytical solutions for three cross covariance functions using the framework described above. Proofs can be found in the technical report ACFR-TR-2011002 referenced above. We have included the sparse covariance function proposed in [Melkumyan and Ramos, 2009] as this provides intrinsically sparse matrices which can be inverted efficiently. The definitions for the squared exponential, Mat´ern (ν = 3/2 see [Rasmussen and Williams, 2006], p.85) and sparse covariance functions are given below

where lSE, lM and lS are the corresponding length scales, H (x) is the Heaviside unit step function and r = |x - x′ |.

5.1 Squared Exponential ×Matern 3/2

Given the squared exponential and Mat´ern 3/2 covariance functions defined in Eq. (11) and Eq. (12) respectively, the cross covariance term can be analytically calculated via the proposed methodology resulting in:

5.2 Matern 3/2 ×Matern 3/2

A combination of two Mat´ern 3/2 kernels results in the following cross covariance term:

where σ12 = 2√l1l2/ (l2 1 - l2 2 ), and l1 and l2 are the length scales of the first and second Mat´ern 3/2 kernels, respectively.

5.3 Sparse × Sparse

One of the main challenges in multi-task learning with Gaussian processes is the computational burden of inverting a large matrix of size ∑Mi=1 Ni where M is the number of tasks, and Ni is the number of points in task i. To tackle this problem in the case of single task GPs, an intrinsically sparse covariance function was proposed in [Melkumyan and Ramos, 2009]. This covariance function has the property of vanishing to zero after a certain distance - represented as a hyperparameter estimated during the learning phase. As the produced covariance matrix is sparse, significant speed ups can be obtained during inversion. Here, we derive an extension of that covariance function for the multi-task case, in combination with itself and other covariance functions. Note that the resulting multitask sparse×sparse covariance function has compact support.

Combination of two sparse kernels with characteristic length-scales l1 and l2 results in

and

where H (x) is the Heaviside unit step function, l = (l1 + l2)/2, lmin = min (l1, l2) and lmax = max (l1, l2).

5.4 Squared Exponential × Sparse

The cross covariance term between squared exponential and sparse kernels can be calculated analytically resulting in

where Re [γ] stands for the real part of γ and

σSE×S = (2π)1/4 √lSE/(6lS).

5.5 Matern 3/2 × Sparse

The cross covariance term between Mat´ern 3/2 and sparse kernels can be calculated analytically resulting in

where

and lM and lS are the length scales of the corresponding covariance functions.

5.6 Squared Exponential × Squared Exponential

The cross covariance term between two squared exponential kernels has an analytical form given by

For the general anisotropic case:

Multidimensional and anisotropic extensions of the other models are possible by taking the product of the cross covariance terms defined for each input dimension.

Note that the examples above do not consider parameters for the amplitude (signal variance) of the covariance functions. This, however, can be added by multiplying blocks of the multi-task covariance matrix by coefficients from a PSD matrix as in [Bonilla et al., 2008].

6 Experiments

6.1 1D Simulation

The first experiment demonstrates the benefits of using the multi-kernel methodology in an artificial 1D problem for two dependent tasks. The observations for the first task are generated from a minus sine function corrupted with Gaussian noise. Only the observations for the second part of the function are used and the objective is to infer the first part from observations of the second task. Observations for the second task were generated from a sine function with some additional complexity to make the function less smooth and corrupted by Gaussian noise. A comparison between independent GP predictions, multi-task GP with squared exponential kernel for both tasks, and the multi-kernel GP (squared exponential kernel for the first task and Mat´ern 3/2 for the second) is presented in Figure 1. It can be observed in Figure 1(c) that the multi-kernel GP models the second function more accurately. This helps in providing a better prediction for the first task.

Despite the simplicity of this experiment it simulates a very common phenomenon in grade estimation for mining. Some elements have a much higher concentration variability but follow the same trend as others. Being able to aptly model these dependencies from noisy x-ray lab samples is essential for an accurate final product. This is empirically demonstrated in the second experiment.

6.2 Iron Ore

1363 samples from an iron ore mine were collected and analyzed in laboratory with x-ray instruments to determine the concentration of three components: iron, silica and alumina. Iron is the main product but equally important is to asses the concentration of the contaminants silica and alumina. The samples were collected from exploration holes of about 200m deep, distributed in an area of 6 km2 . Each hole is divided into 2 meter sections for laboratory assessment, the lab result for each section is then an observation in our dataset. The final dataset consists of 4089 data points representing 31 exploration holes. We separate two holes to use as testing data. For these holes we predict the concentration of silica given iron and alumina. The experiment is repeated employing different multi-task covariance functions with either squared exponential or Mat´ern kernel for each task combined with the cross-covariance terms presented in Section 5. The results are summarized in Table 1 which demonstrates that the dependencies between iron, silica and alumina are better captured by the Mat´ern 3/2 × Mat´ern 3/2 × SqExp multi-kernel covariance function.

6.3 Jura

In the third experiment GPs with different multi-kernel covariance functions were applied to the Jura dataset. The Jura dataset is a benchmark dataset in geostatistics 2 . It consists of a training set with 259 samples in an area of 14.5km2 and a testing set with 100 samples. The task is to predict the concentration of cadmium (Cd), lead (Pb) and zinc (Zn) at new locations. The multi-kernel covariance functions proposed in this paper enable considering different kernels for each of the materials thus maximizing the predictive qualities of the GP. The 259 training samples were used at the learning stage and the 100 testing samples were used to evaluate the predictive qualities of the models. The square root mean square error (SMSE) for all possible triplet combinations of SqExp and Mat´ern 3/2 kernels are presented in Table 2. The results demonstrate that the dependencies between cadmium, lead and zinc are better captured by the Mat´ern 3/2 × SqExp × SqExp triplet-kernel.

2 The Jura dataset is available at: http://goovaerts.pierre.googlepages.com/

6.4 Concrete Slump

In the last experiment the concrete slump dataset 3 is considered. This dataset contains 103 data points with seven input dimensions and 3 outputs describing the influence of the constituent parts of concrete on the overall properties of the concrete. The seven input dimensions are cement, slag, fly ash, water, SP, coarse aggr. and fine aggr., and the outputs are slump, flow and 28-day compressive strength of concrete. 83 data points are used for learning and 20 data points are used for testing. The square root mean square error (SMSE) for all possible triplet combinations of SqExp and Mat´ern 3/2 kernels for this dataset are presented in Table 3. The results demonstrate that the dependencies between slump, flow and 28-day compressive strength of concrete are better captured by the SqExp × Mat´ern 3/2 × Mat´ern 3/2 triplet-kernel.

3 The concrete slump dataset is available at: http://archive.ics.uci.edu/ml/datasets/Concrete+Slump+Test

7 Conclusions

This paper presented a novel methodology to construct cross covariance terms for multi-task Gaussian process. This methodology allows the use of multiple covariance functions for the same multi-task prediction problem. We prove that if a stationary covariance function can be written as a convolution of two identical basis functions, a cross covariance term can always be defined resulting in a positive semi-definite multitask covariance matrix. A general methodology to find the basis function is then developed based on Fourier analysis.

We provide analytical solutions for six combinations of covariance functions, three of them combining different covariance functions. The analytical forms for the cross covariance terms can be directly applied to GP prediction problems but are useful for other kernel machines.

We presented a multi-task sparse covariance function which provides computationally efficient (and exact) way of performing inference in large datasets [Melkumyan and Ramos, 2009]. Note however that approximate techniques such as [Williams and Seeger, 2001; Lawrence et al., 2003; Snelson and Ghahramani, 2006] can also be used.

As a future work we plan to extend the approach to nonstationary covariance functions, possibly combining nonstationary and stationary kernels. This can be useful in applications involving space and time domains such as pollution estimation or weather forecast.

Acknowledgements

This work was supported by the Rio Tinto Centre for Mine Automation and the Australian Centre for Field Robotics.