Directional Cross Derivatives

[daseyb]

I'm working with Gaussian processes where the kernel $k(x, y)$ is a function of two position vectors. I would like to take "directional cross derivatives" with respect to each position. That is, I want to compute something like $D_{v_y} D_{v_x} k(x,y) = D^2_{v_x,v_y} k(x,y)$ where $D_{v_*} k(x,y)$ is the directional derivative of $k(x,y)$ with respect to the change of $*$ in direction $v_*$.

I at first naively did something like
autodiff::derivatives( k(x, y), autodiff::along(dirX, dirY), at(x, y))[2]
but that's subtly different as far as I can tell (i.e. the second derivative along the combined direction $[v_x, v_y]$, instead of the cross derivative).

Is there something built in that would allow me to do this reasonably efficiently?
Thanks!

[daseyb]

Update: I did some working out by hand and found that I can do this relatively cleanly (though a bit wastefully) using the hessian.

$$ D^2_{v_x,v_y} k(x,y) = \sum_{j=1}^d ( \sum_{i=1}^d ( \frac{\partial^2 k(x,y)}{\partial x_i \partial y_j} \cdot v_{x_i} ) \cdot v_{y_j} ) \\ = v_x^T H_{k(x,y)}^{d:2d, 0:d} v_y $$

Where $H_{k(x,y)}^{d:2d, 0:d}$ is the off-diagonal dxd "block" of the full 2dx2d hessian matrix, which is where the "actual" (across the x and y vectors, instead of within) cross derivatives are.

In code, this looks pretty neat, but of course, requires dual2nd instead of real2nd:

Eigen::Matrix3d hess = autodiff::hessian(k(x, y), wrt(x, y), at(x, y)).block(3, 0, 3, 3);
double res = dirX.transpose() * hess * dirY;

The only thing that makes this not 100% satisfying is that I'm throwing away a quarter of the hessian, so there might be a more efficient way to compute the required cross derivatives.