Efficient Computation of Hamiltonian Derivatives in RHMC

[THargreaves]

I thought I'd move this conversation out into its own issue since there are quite a few ideas to talk about.

Originally posted by @nsiccha in #439 (comment)

The trace of the product of two (square) matrices can be computed in O(n^2) without computing the O(n^3) matrix product, see e.g. https://en.wikipedia.org/wiki/Trace_(linear_algebra)#Trace_of_a_product.

Yeah, I guess if we're going to be using G^{-1} for this trace product for each θ it's probably worth to invert it directly and use the Kronecker trick.

Furthermore, do we think that Julia is clever enough to reorder this expression: r' * invG * ∂G∂θᵢ * invG * r? It's in any case unncessary to compute it as it's implemented there. Do the tmp = invG * r once, then compute the tmp' * ∂G∂θᵢ * tmp in one go - I'm sure there's a function for that.

I'd be quite confident that it wouldn't. A tmp variable is what I used form the SSM-RHMC implementation.

[nsiccha]

I do wonder about numerical stability though.

Obviously, "you should never actually compute the inverse of a matrix", but we do. If we do that anyways, we may as well make use of that inverse.

If we were not computing the inverse (and we have a general dense G), is there a better way than actually computing the matrix product d times? That would make the overall complexity of ∂H∂θ O(d^4), no?