# Computing gradients via Gaussian Process Regression

I have a set of noisy data that I am fitting using Gaussian Process Regression via Python's sklearn package. The posterior mean of the GP is essentially my output with an associated error. Based on either the posterior mean or the original data itself, is there a systematic or recommended routine to calculate gradients (i.e. derivative of y respect to x) of the original data via GP?

I was planning to simply apply a basic finite difference approximation of the fit, but am wondering if this is a good idea or if there are better techniques (preferably compatible with sklearn) that permit accurate computation of the gradient and its associated propagating error. I am only seeking a solution in 1 dimension (i.e. gradient in x). But suggestions for calculating gradients in multidimensional space via GP are also welcome.

• Is a 'Gaussian Process Regression' with a posterior mean any different from Bayesian linear regression? You might attract more useful replies if you use simpler language. – Frans Rodenburg Oct 24 '18 at 7:40
• You should not use finite differences because there are better ways to do this. If your kernel is twice differentiable then your GP is differentiable. Have a look at this question and its answer: stats.stackexchange.com/questions/180823/… – g g Oct 24 '18 at 20:24
• Ah and I forgot here you find more theory: link.springer.com/chapter/10.1007/978-3-642-21738-8_20 – g g Oct 24 '18 at 20:36
• @gg Awesome, that's almost exactly what I was searching for. If that can be supplemented with a simple example code implementation (e.g. computing derivatives of a GP with SE kernel on some random data set), that would answer my question. – Mathews24 Oct 25 '18 at 12:20

Gaussian process regression (GPR) gives a posterior distribution over functions mapping input to output. We can differentiate to obtain a distribution over the gradient. Below, I'll derive an expression for the expected gradient. There's no need to use finite differencing, as it can be computed in closed form (as long as the covariance function is differentiable; otherwise it doesn't exist).

Assume the model:

$$y = f(\mathbf{x}) + \epsilon, \quad \epsilon \underset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma_n^2)$$

where the observed output $$y \in \mathbb{R}$$ is a function of input $$\mathbf{x} \in \mathbb{R}^d$$, plus i.i.d. Gaussian noise with variance $$\sigma_n^2$$. Say we fit a GPR model with differentiable covariance function $$k$$. Let $$X = \{\mathbf{x_1}, \dots, \mathbf{x_n}\}$$ denote the training inputs, and let $$\mathbf{y} = [y_1, \dots, y_n]^T$$ denote the corresponding training outputs. Let $$\mathbf{x_*}$$ denote a new input, and let $$f_*$$ be a random variable representing the function value at $$\mathbf{x_*}$$.

We want to compute $$E[\nabla f_* \mid X, \mathbf{y}, \mathbf{x^*}]$$, the expected gradient of the function evaluated at $$\mathbf{x_*}$$ (where the gradient is taken w.r.t. the input and the expectation is over the GPR posterior distribution). Because differentiation is a linear operation, this is equivalent to $$\nabla E[ f_* \mid X, \mathbf{y}, \mathbf{x_*}]$$, the gradient of the expected function value (i.e. posterior mean) at $$\mathbf{x_*}$$.

The expected function value at $$\mathbf{x_*}$$ is:

$$E[f_* \mid X, \mathbf{y}, \mathbf{x_*}] = \sum_{i=1}^n \alpha_i k(\mathbf{x_i}, \mathbf{x_*})$$

where $$\mathbf{\alpha} = (K + \sigma_n^2 I)^{-1} \mathbf{y}$$, $$I$$ is the identity matrix, and matrix $$K$$ contains the covariance for all pairs of training points ($$K_{ij} = k(\mathbf{x_i}, \mathbf{x_j})$$). For details, see chapter 2 of Rasmussen and Williams (2006).

$$\nabla E[f_* \mid X, \mathbf{y}, \mathbf{x_*}] = \nabla \sum_{i=1}^n \alpha_i k(\mathbf{x_*}, \mathbf{x_i})$$

$$= \sum_{i=1}^n \alpha_i \nabla k(\mathbf{x_*}, \mathbf{x_i})$$

Note that the weights $$\mathbf{\alpha}$$ are the same as used to compute the expected function value at $$\mathbf{x^*}$$. So, to compute the expected gradient, the only extra thing we need is the gradient of the covariance function.

For the squared exponential covariance function

As an example, the squared exponential (a.k.a. RBF) covariance function with signal variance $$\sigma_f^2$$ and length-scale $$\ell$$ is:

$$k(\mathbf{x}, \mathbf{x'}) = \sigma_f^2 \exp \left[ -\frac{\|\mathbf{x}-\mathbf{x'}\|^2}{2\ell^2} \right]$$

Taking $$k(\mathbf{x_*}, \mathbf{x_i})$$ and differentiating w.r.t. $$\mathbf{x_*}$$ gives:

$$\nabla k(\mathbf{x_*}, \mathbf{x_i}) = k(\mathbf{x_*}, \mathbf{x_i}) \frac{\mathbf{x_i} - \mathbf{x_*}}{\ell^2}$$

This can be plugged into the expression above for the expected gradient.

Example

Here's an example for the 1d function $$f(x) = \sin(2 \pi x)$$. I fit a GPR model with squared exponential covariance function to 200 noisy observations. The noise variance and kernel parameters (signal variance and length-scale) were estimated by maximizing the marginal likelihood. The expected gradient (computed as above) is similar to the true gradient $$\nabla f(x) = 2 \pi \cos (2 \pi x)$$.

• Just to clarify a couple points: this is assuming a stationary kernel, correct? For example, if the covariance was dependent upon the specific $\bf x$, then the weights $\alpha$ would be dependent upon $\bf x$, correct? Also, I believe you may have placed an extra $\nabla$ when you first take the gradient of the expected value function. – Mathews24 Feb 10 '19 at 19:25
• Is it possible to provide the code you used to generate the above example? I am having trouble implementing the derivative of the covariance function. The issue I have is xi and x* have a different number of rows so can't be subtracted. It is probably just a case of me no understanding the notation. – Robin Dec 9 '19 at 5:21
• @Robin I doubt I even saved the code; it was just a quick example from more than a year ago. $x_i$ is a single data point from the training set and $x_*$ is a single new data point (where we'd like to evaluate the gradient). So, they should have the same dimensions. In the example, where input space is one-dimensional, they should both be scalars. – user20160 Dec 9 '19 at 6:33
• @Mathews24 I'm seeing your comment pretty late, but if you're still around...The expression for the gradient should be valid for any kernel; no need to assume it's stationary. The gradient is taken w.r.t. the new input $x_*$, and the weights $\alpha$ don't depend on this; they only depend on the training data. – user20160 Dec 9 '19 at 6:39
• Thanks @user20160 out of curiosity could you calculate the variance of the gradient following the same method you use to calculate the variance of a regular GP? – Robin Dec 19 '19 at 7:32

I don't have enough karma to comment on the above solution by @user20160 , so I'm posting this here. This provides the source code to implement the definition given by @user20160 for the gradient using GPR in sklearn.

Here is a basic working example using an RBF kernel:

gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)
gp.fit(X, y)

# gets 'l' used in denominator of expected value of gradient for RBF kernel
k2_l = gp.kernel_.get_params()['k2__length_scale']

# not necessary to do predict, but now y_pred has correct shape
y_pred, sigma = gp.predict(x, return_std=True)

# allocate array to store gradient

# set of points where gradient is to be queried
x = np.atleast_2d(np.linspace(-5, 0.8, 1000)).T

# loop over each point that a gradient is needed
for key, x_star in enumerate(x):
# eval_gradient can't be true when eval site doesn't match X
# this gives standard RBF kernel evaluations

# x_i - x_star / l^2
x_diff_over_l_sq = ((X-x_star)/np.power(k2_l,2)).ravel()

# pair-wise multiply
intermediate_result = np.multiply(k_val, x_diff_over_l_sq)

# dot product intermediate_result with the alphas
final_result = np.dot(intermediate_result, gp.alpha_)

# store gradient at this point


Re kuberry’s implementation of user20160’s answer: isn’t de-normalization missing here? See this line in sklearn’s GaussianProcess.predict():

 y_mean = self._y_train_std * y_mean + self._y_train_mean


So I guess the final_result should also be multiplied by the normalization constant:

 final_result *= gp._y_train_std


I would like to add in my code as well. It computes the the first and second derivates as well as antiderivates of the Process.

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF,ConstantKernel
from scipy.special import erf

#Gives a fitted Gaussian Process object that can then be used for predictions.
#The Input is of the Form x.shape = (n), y.shape = (n,t) where both x and y
#are np.ndarrays.
#The normalisation has to be set to False for now since it didn't work with
#my current version of sklearn. Could be added in customary by normalizing the
#input data and denormalizing the output directly.
#The Kernel types (not their parameters though) have to stay this way since the derivates
#and antiderivates are computed for this setup. Should no constant kernel be
#desired its parameters can be set to constant_value = 1.0 and
#constant_value_bounds = 'fixed'.
#All other values, as n_restarts, the RBF kernel and Constant kernel parameters
#have to be selected according to the input data.

class GPR:
def __init__(self,x,y):
normalize = False #hardcoded, don't change.
n_restarts = 2

k1 = ConstantKernel(constant_value=1.0,constant_value_bounds=(1e-5,1e5))
k2 = RBF(length_scale=0.1,length_scale_bounds=(1e-5,1e5))

self.gp = GaussianProcessRegressor(k1*k2,
n_restarts_optimizer=n_restarts,
normalize_y=normalize).fit(x.reshape(-1,1),y)

def predict(self,x,k=0):
#x of shape (m)

#returns the gp predictions where f is the true function and
#df, ddf, If, IIf are its first and second derivate respectively antiderivates
#the outputs are the predictions f_p,df_p,ddf_p,If_p,IIf_p where
#f(x) = f_p(x), df(x) = df_p(x), ddf(x) = ddf_p(x), If(x) = If_p(x) + C1,
#IIf(x) = IIf_p(x) + C1*x + C2 with some constants C1,C2
#set k = 0 for the normal prediction, K = 1,2 for the first or second derivates
#and k = -1,-2 for the first or second antiderivates

x = x.reshape(-1,1)

X = x - self.gp.X_train_.reshape(1,-1)
c = self.gp.kernel_.k1.constant_value
l = self.gp.kernel_.k2.length_scale
A = self.gp.alpha_

f = np.exp(-(X)**2 / (2*l**2))
df = (f * (-X / l ** 2))
ddf = f * ((-X / l ** 2)**2 + -1/l**2)
If = np.sqrt(np.pi/2) * l * erf(X/(np.sqrt(2)*l))
IIf = X * If + l**2 * f

if k == 0:
return c * f @ A
elif k == 1:
return c * df @ A
elif k == 2:
return c * ddf @ A
elif k == -1:
return c * If @ A
elif k == -2:
return c * IIf @ A
else:
raise Exception('Unknown parameter k: {}'.format(k))