# 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