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. 
 A: 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
y_pred_grad = 0.0*y_pred;

# 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
    k_val=gp.kernel_(X, np.atleast_2d(x_star), eval_gradient=False).ravel()

    # 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
    y_pred_grad[key] = final_result

A: 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

A: 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))

A: 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).
Expression for the expected gradient
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).
Taking the gradient, we have:
$$\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)$.

