Unstable factorization by default used to sample Gaussian process

[Toby Isaac]

Hi MUQ team,

I'm going through a notebook demonstrating Gaussian processes in MUQ. If I take the example demonstrating the squared exponential kernel and simply change the marginal variance, the samples should look the same, only rescaled. Changing the variance from 1 to 3, however, gives samples that looks like this:

This is pretty clearly due to an unstable factorization. My best from looking at the source is that llt from eigen is used.

Is there an interface to switch to a more stable approach already? If not, there should be.

[Matthew Parno]

Thanks for pointing this out. You’re correct that we’re using llt. It will likely be more stable to use ldlt instead or even an eigenvalue decomposition of the covariance instead. I’ll do some speed tests to see what can be done to make this more stable without harming performance. I agree that maybe we should include an option for more stable factorizations if there is a significant performance hit with the more stable factorizations.

[Matthew Parno]

This has recently been fixed in the GaussianProcess class, but still needs to be fixed in the Gaussian class.