I'm trying to figure out how to calculate the cross-terms of the covariance matrix between the function values and derivatives of a Gaussian Process. For context, this is needed for Gaussian Process regression with positivity and monotonicity constraints, where the derivative and/or function values are constrained (see Da Veiga and Marrel 2020).
Specifically, this means calculating the posterior (predictive) covariance between the function values and derivative at a set of 'constraint points' $X_z$: $$\mathrm{Cov}(f(X_z), \frac{\partial f(X_z)}{\partial x}) $$
So far I have found this link (which references a pdf written by McHutchon for which the URL is now broken) for calculating the variance of the derivative of the Gaussian Process $\mathrm{Cov}(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial x})$, but I've struggled to derive any expression for the combined covariance of $f(x)$ and $\frac{\partial f}{\partial x}$