If L is the L from a Cholesky decomposition, yeah, that’s lower triangular, but I dunno if you can easily get to them from non-square matrices from the expansion.
Oooh, I see. Do you need all the eigenvalues/eigenvectors? Or just the first few?
There’s a section on using SVDs in that Fasshauer book, but I dunno if they’re coming at it from the same angle as you. SVDs certainly aren’t cheap.
How many hyperparameters you have on these GPs? If by some accidental luck you only have like one or maybe two, then maybe you could try building lookup tables to interpolate on the quantities you want?
I did this once for some 1D RBF GPs. The input to the lookup table was the length scale and the output were the Cholesky factors (so you don’t have to run Choleskys[ies?] at runtime – just interpolations). I just stole some eqs. for cubic interpolation off Wikipedia and added some custom autodiff.
Just going by this: Bad LKJ Initialization in model from "Fast hierarchical Gaussian processes" - #3 by andre.pfeuffer, it looks like you have one hyperparameter? bw1? var1 could come outside the eigenvalue decomposition? I dunno. You can do a better job than me of knowing if this makes sense or not with what you’re doing.
I talked to Dan and Aki about this and they basically said this sorta trick isn’t too uncommon, it’s just really limited because interpolation is hard beyond a couple dimensions.
Anyway here are some plots of the lookup table GP vs. brute forcing the full GP:
And here is the stan model: https://github.com/bbbales2/gp/blob/master/models/cubic_interpolated_gp.stan
I did the interpolation with some custom C++ stuff. Matrix stuff can really clog up the autodiff stack: https://github.com/bbbales2/gp/blob/master/models/cubic_interpolated_gp.hpp
And here is the R test code: https://github.com/bbbales2/gp/blob/master/test_interpolate.R