Background and problem
I am using Gaussian Processes (GP) for regression and subsequent Bayesian optimization (BO). For regression I use the gpml package for MATLAB with several custom-made modifications, but the problem is general.
It is a well-known fact that when two training inputs are too close in input space, the covariance matrix may become not-positive definite (there are several questions about it on this site). As a result, the Cholesky decomposition of the covariance matrix, necessary for various GP computations, may fail due to numerical error. This happened to me in several cases when performing BO with the objective functions I am using, and I'd like to fix it.
AFAIK, the standard solution to alleviate ill-conditioning is to add a ridge or nugget to the diagonal of the covariance matrix. For GP regression, this amounts to adding (or increasing, if already present) observation noise.
So far so good. I modified the code for exact inference of gpml so that whenever the Cholesky decomposition fails, I try to fix the covariance matrix to the closest symmetric positive definite (SPD) matrix in Frobenius norm, inspired by this MATLAB code by John d'Errico. The rationale is to minimize intervention on the original matrix.
This workaround does the job, but I noticed that the performance of BO reduced substantially for some functions -- possibly whenever the algorithm would need to zoom-in in some area (e.g., because it's getting nearer to the minimum, or because the length scales of the problem become non-uniformly small). This behaviour makes sense since I am effectively increasing the noise whenever two input points get too close, but of course it's not ideal. Alternatively, I could just remove problematic points, but again, sometimes I need the input points to be close.
I don't think that numerical issues with Cholesky factorization of GP's covariance matrices is a novel problem, but to my surprise I couldn't find many solutions so far, aside of increasing the noise or removing points that are too close to each other. On the other hand, it is true that some of my functions are pretty badly behaved, so perhaps my situation is not so typical.
Any suggestion/reference that could be useful here?