Inverting a Gaussian (RBF) Kernel Matrix Without Numerical Trouble

The Gaussian (RBF) kernel matrix for points $x_1, \dots, x_n$ is

You need its inverse, or more precisely $K^{-1}y$, in kernel ridge regression, Gaussian process regression and RBF interpolation. There are two separate questions: is $K$ invertible in theory, and can you invert it reliably in practice?

In theory: always invertible for distinct points

The Gaussian kernel is strictly positive definite. For any distinct points $x_1, \dots, x_n$ and any non-zero vector $c$,

This follows from Bochner’s theorem: the Gaussian’s Fourier transform is again a Gaussian, which is positive everywhere. So $K$ is symmetric positive definite, all its eigenvalues are positive, and $K^{-1}$ exists. Duplicate points produce identical rows and make $K$ singular.

In practice: often numerically singular

The smallest eigenvalue of $K$ can be astronomically small. That happens when:

• points are close together compared with $\sigma$, so nearby rows are nearly identical, or
• $\sigma$ is large compared with the spread of the data, so every entry is close to 1.

A condition number of $10^{16}$ or more is common, at which point double-precision arithmetic can’t distinguish $K$ from a singular matrix.

Two-point example

With two points at distance $d$, write $k = e^{-d^2/(2\sigma^2)}$:

For $d = \sigma = 1$: $k = e^{-1/2} \approx 0.6065$ and $1 - k^2 = 1 - e^{-1} \approx 0.6321$, so

As $d \to 0$ (or $\sigma \to \infty$), $k \to 1$ and the entries of $K^{-1}$ blow up like $1/(1-k^2)$. The eigenvalues are $1 \pm k$, so the condition number is $(1+k)/(1-k)$. With $n$ points the effect compounds quickly.

What to do instead

1. Add a ridge (nugget, jitter) term. Replace $K$ by $K + \lambda I$ with $\lambda > 0$. Every eigenvalue moves up by $\lambda$, so the condition number is at most $(\lambda_{\max} + \lambda)/\lambda$.

• In kernel ridge regression and GP regression, $\lambda$ is part of the model: the regularisation strength or the noise variance $\sigma_n^2$.
• For pure interpolation, use a tiny jitter such as $10^{-8}$ to $10^{-6}$ times the mean of the diagonal, accepting a tiny loss of exactness.

2. Solve, don’t invert. You almost never need $K^{-1}$ itself. Factor once with Cholesky, $K + \lambda I = LL^\top$, then solve two triangular systems:

import numpy as np
from scipy.linalg import cho_factor, cho_solve
from scipy.spatial.distance import cdist

def rbf_kernel(X, Y, sigma):
    return np.exp(-cdist(X, Y, "sqeuclidean") / (2 * sigma**2))

K = rbf_kernel(X, X, sigma)
c, low = cho_factor(K + lam * np.eye(len(X)))  # fails loudly if not positive definite
alpha = cho_solve((c, low), y)                 # alpha = (K + lam*I)^{-1} y
y_pred = rbf_kernel(X_new, X, sigma) @ alpha

Cholesky costs about $n^3/3$ operations, roughly half the cost of an LU factorisation. It is backward stable for positive definite matrices, and its diagonal gives you $\log\det K = 2\sum_i \log L_{ii}$ for free (needed for GP marginal likelihoods).

3. Choose $\sigma$ sensibly. A common starting point is the median pairwise distance (the “median heuristic”). Then tune $\sigma$ and $\lambda$ by cross-validation or by maximising the marginal likelihood.

4. For large $n$, use low-rank approximations (Nyström, random Fourier features) or iterative solvers (conjugate gradients), which avoid the $O(n^3)$ factorisation.

Common mistakes

• Calling np.linalg.inv(K) and multiplying. It is slower and less accurate than a solve, and it hides ill-conditioning.
• Treating a Cholesky failure as a bug. It is a signal that $K$ is numerically singular. Increase $\lambda$ or check for duplicate points.
• Using the same $\lambda$ after rescaling the data. The appropriate $\lambda$ depends on the scale of $K$ and $y$.