Scikit-learn's Gaussian Processes: How to include multiple hyperparameters in kernel/cov function?

I'm using the scikit-learn's implementation of Gaussian processes. A simple thing to do is to combine multiple kernels as a linear combination to describe your time series properly. So I'd like to include both the squared exponential kernel and the periodic kernel. Linear combinations of valid kernels produce valid kernels, and same goes for multiplying valid kernels (given by Rasmussen and Williams).

Unfortunately I haven't figured out how to give the theta parameters properly to the model. For example, if we have:

$$k_{Gauss}(x,x') = \exp{(\theta (x-x')^2)}$$

then it is alright (this is how the squared-exponential kernel is defined in scikit-learn). But if I wanted:

$$k_{Gauss}(x,x') = \theta_0 \exp{(\theta_1 (x-x')^2)}$$

then it is impossible, it seems. The $\mathbf{\theta}$ thing is supposed to be an array, in case you have multiple dimensions/features (even though scikit-learn doesn't support multidimensional GPs, someone developed it, and it will be merged soon). So there is one row with the columns being the parameter in such-and-such dimension. But you cannot have more rows, otherwise it screams at you.

So question: has anyone actually been able to use kernels that use more than one hyperparameter? If so, what am I doing wrong? And if it is indeed not possible with the current code in scikit, does anyone have some tips on how to extend it so that it can? This is a really important feature that I need. Thanks.

• Thank you Dougal, I will keep that in mind for the future. Oct 15, 2013 at 3:05

On scikit-learn==0.14.1.

$\theta_0$ can be a vector. The following code works for me.

import numpy as np
from sklearn.gaussian_process import GaussianProcess
from sklearn.datasets import make_regression
X, y = make_regression()
model.fit(X,y)


You can pass any kernel you want as the parameter corr. The following is the radial basis function that sklearn uses for Gaussian processes.

def squared_exponential(theta, d):
"""
Squared exponential correlation model (Radial Basis Function).
(Infinitely differentiable stochastic process, very smooth)::

n
theta, dx --> r(theta, dx) = exp(  sum  - theta_i * (dx_i)^2 )
i = 1

Parameters
----------
theta : array_like
An array with shape 1 (isotropic) or n (anisotropic) giving the
autocorrelation parameter(s).

dx : array_like
An array with shape (n_eval, n_features) giving the componentwise
distances between locations x and x' at which the correlation model
should be evaluated.

Returns
-------
r : array_like
An array with shape (n_eval, ) containing the values of the
autocorrelation model.
"""

theta = np.asarray(theta, dtype=np.float)
d = np.asarray(d, dtype=np.float)

if d.ndim > 1:
n_features = d.shape[1]
else:
n_features = 1

if theta.size == 1:
return np.exp(-theta[0] * np.sum(d ** 2, axis=1))
elif theta.size != n_features:
raise ValueError("Length of theta must be 1 or %s" % n_features)
else:
return np.exp(-np.sum(theta.reshape(1, n_features) * d ** 2, axis=1))


It looks like you're doing something pretty interesting, btw.

• Thanks Jacob. Theta can indeed be a vector, however the elements of the vector represent that theta parameter in D-dimensional space. For example, if I wanted a hyperparam $\sigma$ for that SE kernel to be its variance, then that works. However, say for some reason I also needed another hyperparameter $A$ that scaled the entire thing, i.e. I had Aexp(\sigmad**2), then it doesn't work. What if I wanted this kernel with those two hyperparameters but with two features? So we have $\mathbf{\theta} = \{A,\sigma\}$, and in the 1-D case $A$ and $\sigma$ are both scalars, but in a 2-D case they Oct 15, 2013 at 3:10
• must be 2-dimensional vectors themselves. So in terms of scikits code, see the final return value for the SE kernel, it reshapes theta just to make sure it has the proper formatting for a single row vector, where each element is a feature (but that kernel restricts $|\mathbf{\theta}|$ to 1). In a 2-D case I would expect $\mathbf{\theta}$ to be 2x2 matrix, where the first row corresponds to $A$ and the second row corresponds to $\sigma$, where the columns are the dimensions/features. Anyway, it's not actually like that, and it screams at me. I don't know the correct way to do this. Any insight? Oct 15, 2013 at 3:13
• Write the function for corr with the parameters hardcoded. From a programmatic view you probably want to do it as a partial application. From a statistics view you lose the ability to specify upper and lower bounds of theta, which would require a rewrite. Oct 15, 2013 at 3:37
• OK, that will probably be really tedious but maybe do-able for proof of concept at least. I really just wanted to make sure there wasn't another way and I was just not seeing it. Thanks. Oct 15, 2013 at 4:03
• The most obvious way to me, would be write the two separate corr functions that are partially applied on theta. E.g. def square_exp(d), then write a function to combine the results via multiplication or addition with the argument def function(theta, d), where theta doesn't do anything in the function. Oct 15, 2013 at 17:21

For future users, in sklearn 0.19.1 (and probably earlier) you can use combined kernels.

I think you can create the kernel you want by

kernel = ConstantKernel(1.0) * RBF(np.ones(nrOfFeatures))


The ConstantKernel would be your theta0. In case you provide the RBF kernel with a single float, this will be your theta1