Consider a covariance function of the form
$$K_{i,j}=\alpha\times exp(-0.5 (x_i-x_j)^2/l^2)$$
This is a very common function used in Gaussian processes. How to show that this covariance is non-negative definite ?
Solved – Prove that the squared exponential covariance is non-negative definite
covariancegaussian process
Related Solutions
You could try following first the standard way of deriving the solution for $k(s,t) = e^{-|s-t|}$ given for instance in [1] and then try to reapply this to your case. If you don't have the access to [1] I can give an outline here. Generally you have to carefully differentiate
$\int k_{SE}(s,t)f_i(s)dt = \lambda_i f_i(s)$
with respect to $t$ and see if it simplifies to an ODE for $f_i$. Your ODE is going to be more complicated but should be solvable (haven't done the calculation myself!).
[1] Ghanem, R. and Spanos, P. "Stochastic Finite Elements: A Spectral Approach", 1991 Springer, pp 29-33
I think the comment given in the linked answer is incorrect or might be a misunderstanding. In general covariance matrices just need to be positive semi-definite. But the covariance matrix $\Sigma$ constructed in the specific way you did from the RBF kernel function will always be strictly positive definite. This means $x^T \Sigma x > 0$ unless $x=0$. This fact is crucial. Because if your $\Sigma$ were just semi-definite, it would not be invertible. But you need the inverse to solve the interpolation or prediction equation in Gaussian Process regression or interpolation.
So, if $\Sigma$ is positive definite, why did you run into trouble? Most likely you have a numerical problem. The RBF kernel is notorious for being numerically unstable. The reason is its fast (exponential!) decay of eigenvalues. Very small eigenvalues lead to bad condition numbers and numerical problems. Explicit formulas for eigenvalues and eigenfunctions of the RBF kernel can be found in the book Rasmussen-Williams Chapter 4.3.1. In the chapter before they also explain the connection between rapidly decaying eigenvalues and smoothness.
I do not know the algorithm used by Mathematica to determine positive definiteness but any algorithm is doomed to fail if the condition number of $\Sigma$ becomes too large. You can check this yourself, by decreasing the number of points in your grid while checking the condition number or plotting the eigenvalues.
This also explains the widely used "trick" to add a small $\epsilon$ to the diagonal. By increasing the smallest eigenvalue, the condition number of $\Sigma + \epsilon I$ will be (much) smaller than the condition number of $\Sigma$ alone. Again this is something you can easily check while varying grid granularity.
Other kernels are better behaved in this respect. If you really want to do interpolation you should probably choose a different kernel. One example would be a kernel from the Matern family (see Rasmussen-Williams, Formula 4.14)
Best Answer
I am not an expert but I'll sketch a standard argument which is explained in more detail in Rasmussen and Williams, Chapter 4 Section 2.1 (that book has answered a ton of my question about GPs). So we are working with the squared exponential function right? We have: $$K_{i,j}= \alpha \cdot \mathrm{exp}(\frac{-(x_i-x_j)^2}{2\ell^2}) = \alpha \cdot \mathrm{exp}(\frac{-(|x_i-x_j|)^2}{2\ell^2})$$
Since kernel can be written as a function of $|x_i-x_j|$, it is stationary (isotropic, even). Since it is stationary, the trick is we can apply Bochner's theorem to $K_{i,j}$. In this case, showing positive semidefiniteness of the square exponential reduces to finding a suitable function $S(s)$ which we can take the Fourier transform $\mathcal{F}_s$ such that $\mathcal{F}_s(S(s))=K_{i,j}$. Now the Fourier transform of a Gaussian is another Gaussian, so the $S(s)$ function that we are looking for turns out to be $$ S(s) = \alpha (2\pi \ell^2)^{D/2} \mathrm{exp}(−2\pi^2 \ell^2s^2). $$ If you calculate the Fourier transform of this function you will get your kernel, thus showing it's positive semidefinite. I'm sorry if that's too terse, but I can try to derive more details if that helps.