Cholesky, in my opinion, is not robust as you have discovered. When it encounters a non-psd matrix it gives up. Maybe that is the correct course of action, but I still think that when you are using the getCovariance to estimate the curvature in the neighborhood of a point it would be okay to take the generalized inverse and not worry too much about why your matrix is bordering on non-PSDness... If your optimization stops on that point, that is another story and should be flagged.
> A possibly more robust option here is to use Cholesky decomposition, > > which is known to be stable for symmetric positive definite > > matrices, which the covariance matrix being inverted here should > > be. The exceptions thrown will be different; but they will give > > more specific information about what is wrong with the covariance > > matrix. > > I've tried it with my problem, and it also throws an exception. > However, I would like to obtain the covariance matrix anyway, because I've > no other clue as to what might be wrong. > So I think that, at least, users should be able to set the positive > definiteness threshold in order to avoid raising an exception. > > >