Gaussian Process predictive gradient with GPML

Given a Gaussian Process $$f(x) = \mathcal{GP}\big( m(x), k(x,x') \big)$$ with $$x\in\mathbb{R}^d$$, I would like to calculate $$\mathbb{E}\left[ \frac{\partial f}{\partial x_i} \right]$$; that is, the $$i^\mathrm{th}$$ coordinate of $$\mathrm{grad} f$$.

For short, I use the notation $$\partial_i \bullet$$ for the $$i^\mathrm{th}$$ gradient coordinate below. Following the GPML book by Rasmussen and Williams p.16, I think this can be done as follows: $$\mathbb{E}\left[ \partial_i f_*\ \big| X, f, X_* \right] = \partial_i \mathbb{E}\left[ f_*\ \big| X, f, X_* \right] = \partial_i K(X_*, X) K(X,X)^{-1} f$$ where $$X$$ are coordinates with known function values $$f$$, $$X_*$$ is a matrix of coordinates with unknown function values $$f_*$$, and we want the $$i^\mathrm{th}$$ coordinate of the gradient at these points $$\partial_i f_*$$.

In short, we simply need the derivative of the covariance function with respect to that coordinate. For example with the (isotropic) squared-exponential covariance: $$\partial_i k_\mathrm{SE}(x,x') = \partial_i k_\mathrm{SE}(r = \|x -x'\|_2) = -\frac{2r (x_i - x_i')}{L^2}\exp\left(-\frac{r^2}{2L^2}\right)$$

I was able to program this, but I don't want to reinvent the wheel and it seems likely that there would be a way to compute this with the GPML toolbox already. The help for covFunctions indicates that the second output dK is the "linear directional derivative function", which can apparently be called with dK(Q); but I am not sure what this means, or how to use it. Hence my questions:

Q1. Is it correct that the gradient of a GP with respect to its variables (not its parameters, as often found elsewhere in the context of GP training) is simply obtained as above, given the gradient of the covariance function?

Q2. How can one obtain these predictions using the GPML toolbox (if at all possible)?

• have you found your answer how to do this? Or did you find that you just don't need it? :-) – Curious Dec 6 '19 at 9:43

First, let's start with a note: you call X matrix of "coordinates"; these are rather covariates, in general.

Next, you came to the $$\partial_i K(X_*, X) K(X,X)^{-1} f$$, which seems correct. But in the further text you are just considering the derivation of $$K$$, while you still need the derivation of $$K^{-1}$$ and $$f$$ as well, which I guess will be much more complicated.

So the answer to Q1 is NO, you are still missing much more pieces of the equation.

Regarding the Q2 and the $$\partial K \over \partial x_i$$ (i.e. your $$\partial_i k(x,x')$$) ... first, I cannot wrap my head around this, it seems like there is some principial problem there... since when you derive $$k(x,x')$$ by $$x$$ then by which $$x$$: $$x$$ or $$x'$$? Second, I am not user of gpml and matlab, but every GP software package needs to have implemented $$\partial K \over \partial l_i$$ (you can e.g. see the R implementation here) but I doubt it has implemented the $$\partial K \over \partial x_i$$, because of the difficulties I mentioned in the previous points.

EDIT: Q2: I actually found the computation of $$\partial K \over \partial x_i$$ ! Not in GPML, but in another Matlab software, GPstuff: see function gpcf_exp_ginput():

%    DKff = GPCF_EXP_GINPUT(GPCF, X) takes a covariance function
%    structure GPCF, a matrix X of input vectors and returns
%    DKff, the gradients of covariance matrix Kff = k(X,X) with
%    respect to X (cell array with matrix elements). This subfunction
%    is needed when computing gradients with respect to inducing
%    inputs in sparse approximations.