[FIXED] Vectorized evaluation of sklearn.gaussian_process.kernels.Matern

Issue

It is unclear for me from the documentation how the __call__ function of the Matern kernel in sklearn.gaussian_process.kernels works. In particular, I would like to build the matrix K = k(x_i, x_j) for all couples of elements x_i, x_j in a grid x. The intuitive way, for me, would be to do a meshgrid operation and then feed it to the kernel. Apparently this does not work (see pictures), but with two nested loops I obtain what I expect (again see pictures). How do I avoid the double loop in the code below?

from sklearn.gaussian_process.kernels import Matern
import numpy as np
import matplotlib.pyplot as plt

n = 51
x = np.linspace(0, 1, n)
kernel = Matern(length_scale=1, nu=1.5)
xx, yy = np.meshgrid(x, x)

# evaluate the kernel vectorially
k1 = kernel(xx, yy)

# evaluate the kernel with for loops
k2 = np.zeros([n, n])
for i in range(n):
    for j in range(n):
        k2[i, j] = kernel([[x[i]]], [[x[j]]])

plt.matshow(k1)
plt.matshow(k2)

Edit:

I (almost) made this work by using

x_eval = np.array([x, x]).T
k1 = kernel(x_eval)

Still there is some non-negligible difference between the vectorized and loop-based versions. Some ideas why?

Edit 2:

The code to get k1 in the Edit above is wrong, I checked with a self-made implementation of the Gaussian kernel K(x_i, x_j) = exp(-(x_i - x_j)**2 / (2 * scale**2)) which is the same as the Matern covariance with parameter nu=np.inf. The value I obtain coincides with the for loop version, and is different than the other.

Solution

Correct implementation:

x_eval = np.reshape(x, [n, 1])
k1 = kernel(x_eval)

Answered By - G. Gare