When doing a Cholesky decomposition of a covariance matrix with very low eigenvalues, numpy.linalg.cholesky and scipy.linalg.cholesky say the matrix is not positive definite given that it is. I'm guessing this is a numerical precision issue.
In the following code, I create $25$ correlated variables as follows:
Each random variable corresponds to a cell in a $5\times 5$ grid, and the covariance between $X_i$ and $X_j$ is $\sigma^2 e^{-d^2/2\epsilon^2}$ where $\sigma$ and $\epsilon$ are given, and $d$ is the distance on the grid between cell $i$ and cell $j$. For example, the distance between cells $1$ and $25$ (assuming left-right, top-bottom numbering, i.e. cell $1$ is $(1,1)$, and cell $25$ is $(5,5)$) is $d=\sqrt{(5-1)^2+(5-1)^2}\approx 5.66$. In the code, I start by constructing my $25\times 25$ matrix $\textbf D$ where element $D_{ij}$ is the distance between cells $i$ and $j$. Then I just apply the gaussian element-wise.
import numpy as np
from scipy.spatial.distance import cdist
from scipy.linalg import cholesky
def gauss_rbf(d, epsilon, sigma):
return sigma**2*np.exp(-0.5*(float(d)/epsilon)**2)
x,y = np.meshgrid(range(5), range(5))
coords = zip(x.ravel(),y.ravel())
dists = cdist(coords,coords)
epsilon = 10
sigma = 5
cov = gauss_rbf(dists, epsilon, sigma)
cholesky(cov)
I get this: LinAlgError: 25-th leading minor not positive definite
This is because np.linalg.eigvalsh(cov) shows there's a very small eigenvalue (-8e-15). But this is a numerical error. The covariance matrix is really positive definite. How can I fix it?
