Covariance of derivative of Gaussian Process Regression

Question

There are a quite a few questions and answers which discuss how to calculate the gradients/derivatives of the posterior of Gaussian Process Regression (see here, here). These include the equations for calculating the mean but the calculation of the covariance is less clear. How can we calculate the covariance of the gradient?

Answer 1

You may find McHutchon (2013) useful; everything you need is there, but in case the link goes dead, I'll put a streamlined version here.

As you noted, multiple answers here cover

$$ \mathbb{E} \left[ \frac{\partial \mathbf{f}_\ast}{\partial \mathbf{x}_\ast} \right] = \frac{\partial k\left(\mathbf{x}_\ast, \mathbf{X}\right)}{\partial \mathbf{x}_\ast} K^{-1} \mathbf{y}, $$

but how do we get

$$ \mathbb{V} \left[ \frac{\partial \mathbf{f}_\ast}{\partial \mathbf{x}_\ast} \right]? $$

We consider an additional test point $\mathbf{x}_\ast + \boldsymbol\delta$. Then

\begin{align} f \left( \mathbf{x}_\ast \right) & = \bar{f} \left( \mathbf{x}_\ast \right) + \mathbf{z}_\ast \\ f \left( \mathbf{x}_\ast + \boldsymbol\delta \right) & = \bar{f} \left( \mathbf{x}_\ast + \boldsymbol\delta \right) + \mathbf{z}_\delta \\ \end{align}

and

$$ \begin{bmatrix} \mathbf{z}_\ast \\ \mathbf{z}_\delta \end{bmatrix} \sim \mathcal{N} \left( \mathbf{0}, \begin{bmatrix} k_{\ast\ast} - \mathbf{k}_\ast^T K^{-1} \mathbf{k}_\ast & k_{\ast\delta} - \mathbf{k}_\ast^T K^{-1} \mathbf{k}_\delta \\ k_{\delta\ast} - \mathbf{k}_\delta^T K^{-1} \mathbf{k}_\ast & k_{\delta\delta} - \mathbf{k}_\delta^T K^{-1} \mathbf{k}_\delta \\ \end{bmatrix} \right). $$

Taking the limit as $\boldsymbol\delta \to \mathbf{0}$,

\begin{align} \frac{\partial \mathbf{f}_\ast}{\partial \mathbf{x}_\ast} & = \lim_{\boldsymbol\delta \to \mathbf{0}} \frac{f \left( \mathbf{x}_\ast + \boldsymbol\delta \right) - f \left( \mathbf{x}_\ast \right)}{\mathbf{x}_\ast + \boldsymbol\delta - \mathbf{x}_\ast} \\ & = \frac{\partial \bar{\mathbf{f}}_\ast}{\partial \mathbf{x}_\ast} + \lim_{\boldsymbol\delta \to \mathbf{0}} \frac{\mathbf{z}_\delta - \mathbf{z}_\ast}{\boldsymbol\delta}, \end{align}

we find

$$ \mathbb{V} \left[ \lim_{\boldsymbol\delta \to \mathbf{0}} \frac{\mathbf{z}_\delta - \mathbf{z}_\ast}{\boldsymbol\delta} \right] = \mathbb{V} \left[ \frac{\partial \mathbf{f}_\ast}{\partial \mathbf{x}_\ast} \right] = \frac{\partial^2 k \left(\mathbf{x}_1^\ast, \mathbf{x}_2^\ast \right)}{\partial \mathbf{x}_1^\ast \partial \mathbf{x}_2^\ast} - \frac{\partial k \left(\mathbf{x}_\ast, \mathbf{X} \right)}{\partial \mathbf{x}_\ast} K^{-1} \frac{\partial k \left(\mathbf{x}_\ast, \mathbf{X} \right)}{\partial \mathbf{x}_\ast}^T . $$

Please note that in

$$ \frac{\partial^2 k \left(\mathbf{x}_1^\ast, \mathbf{x}_2^\ast \right)}{\partial \mathbf{x}_1^\ast \partial \mathbf{x}_2^\ast}, $$

$\mathbf{x}_1^\ast = \mathbf{x}_2^\ast = \mathbf{x}_\ast$, but we must do it in this cross-partial way to avoid negatives on the diagonal. All credit for this presentation of the derivation goes to McHutchon (2013), like I say, I simply reproduce relevant parts here for completeness of the answer.