Description

The RBF kernel corresponds to the form:
k(x_i, x_j) = exp(-1/2 ||x_i - x_j||^2 / length_scale^2)

Therefore, the gradient with respect to the parameter length_scale should be:
gradient = k(x_i, x_j) * ( ||x_i - x_j||^2 / length_scale^3)

However, the current implementation seems to use the form:
gradient = k(x_i, x_j) * ( ||x_i - x_j||^2 / length_scale^2)

Steps/Code to Reproduce

Example:

import numpy as np
from sklearn.gaussian_process.kernels import RBF
np.random.seed(1)
X = np.array([[1,2], [3,4], [5,6]])
sk_kernel = RBF(2.0)
K_grad = sk_kernel(X, eval_gradient=True)[1][:,:,0]

Expected Results

K_grad = 
array([[ 0.        ,  0.36787944,  0.07326256],
       [ 0.36787944,  0.        ,  0.36787944],
       [ 0.07326256,  0.36787944,  0.        ]])

Actual Results

K_grad = 
array([[ 0.        ,  0.73575888,  0.14652511],
       [ 0.73575888,  0.        ,  0.73575888],
       [ 0.14652511,  0.73575888,  0.        ]])

Versions

Darwin-14.5.0-x86_64-i386-64bit
Python 3.6.3 (default, Oct 8 2017, 15:07:13)
[GCC 4.2.1 Compatible Apple LLVM 7.0.0 (clang-700.0.72)]
NumPy 1.13.3
SciPy 0.19.1
Scikit-Learn 0.19.0