Hierarchical GPs with varying lengthscales

That was my question :slight_smile:

I’ll post my faster method when I have it more rigorously quantified, but it involves computing the cholesky decomposition more efficiently through knowledge of the symmetric-Toeplitz structure of the correlation matrix, then computing the matrix-multiply more efficiently through knowledge that it’s a lower-tri and using the foward-substitution algorithm.

Thanks for the heads-up on vectorize(); I hadn’t come across it yet and should clearly add it to my toolbox.