Wrong choice of covariance function? I am applying Gaussian process regression (GPR) model on some data assuming covSEiso (a.k.a. RBF) covariance function.  I made sure there are no identical data points (but there are similar data points) and gpml still runs into numerical problems when applying Cholesky decomposition on $K+\sigma_{noise}^2$ I matrix.  With other choice of covariance function (e.g. Matern with d parameter 1) things run smoothly.  Does this mean the GRP with covSEiso was a wrong model to begin with?  I want to decide weather to learn how to deal with the numerical problem (e.g. by removing similar data points which is a hassle IMO) or search for a suitable covariance function.
 A: The rbf kernel is extremely smooth and despite it being possibly the most popular choice of kernel it is rarely the best choice (not just IMO but this also seems to be the opinion of most well know researchers in the GP community - watch talks from GPSS for example). 
However, no matter what kernel you use there is always a chance that your covariance matrix will be poorly conditioned due to the distribution of the input values - not anything to do with how well the true function is modelled! This is kind of annoying in practice but there are loads of ways to get around it. 
The first I don't recommend but is by far the most common hack which is to increase the iid noise (like jitter chol which is build into GPML, GPy, PyGPs... etc). They refer to this extra noise as jitter.
What makes a lot more sense is to use a low rank representation of the covariance matrix. Numerical problems occur when eigen values are close to zero and hence the conditioning is poor. So the first approach adds a constant to every eigen value and the second truncates the spectrum and hence removes the smallest values.
