What is the rationale of the Matérn covariance function?

Question

The Matérn covariance function is commonly used as kernel function in Gaussian Process. It is defined like this

$$ {\displaystyle C_{\nu }(d)=\sigma ^{2}{\frac {2^{1-\nu }}{\Gamma (\nu )}}{\Bigg (}{\sqrt {2\nu }}{\frac {d}{\rho }}{\Bigg )}^{\nu }K_{\nu }{\Bigg (}{\sqrt {2\nu }}{\frac {d}{\rho }}{\Bigg )}} $$

where $d$ is a distance function (such as Euclidean distance), $\Gamma$ is the gamma function, $K_\nu$ is the modified Bessel function of the second kind, $\rho$ and $\nu$ are positive parameters. $\nu$ is a lot of time chosen to be $\frac{3}{2}$ or $\frac{5}{2}$ in practice.

A lot of time this kernel works better than the standard Gaussian kernel as it is 'less smooth', but except that, are there any other reason why one would prefer this kernel? Some geometric intuition about how it behaves, or some explanation of the seemingly cryptic formula would be highly appreciated.

Answer 1

In addition to @Dahn's nice answer, I thought I would try to say a little bit more about where the Bessel and Gamma functions come from. One starting point for arriving at the covariance function is Bochner's theorem.

Theorem (Bochner) A continuous stationary function $k(x, y) = \widetilde{k}(|x − y|)$ is positive definite if and only if $\widetilde{k}$ is the Fourier transform of a finite positive measure: $$\widetilde{k}(t) = \int_{\mathbb{R}} e^{−iωt}\mathrm{d}µ(ω) .$$

From this you can deduce that the Matérn covariance matrix is derived as the Fourier transform of $\frac{1}{(1+\omega^2)^p}$ (Source: Durrande). That's all good but it doesn't really tell us how you arrive at this finite positive measure given by $\frac{1}{(1+\omega^2)^p}$. Well, it's the (power) spectral density of a stochastic process $f(x)$.

Which stochastic process? It's known that a random process on $\mathbb{R}^d$ with a Matérn covariance function is a solution to the stochastic partial differential equation (SPDE) $$ (κ^2 − ∆)^{α/2} X(s) = φW(s), $$ where $W(s)$ is Gaussian white noise with unit variance, $$\Delta = \sum_{i=1}^d \frac{\partial^2}{\partial x^2_i}$$ is the Laplace operator, and $α =ν + d/2$ (I think this is in Cressie and Wikle).

Why pick this particular SPDE/stochastic process? The origin is in spatial statistics where it's argued that this is the simplest and natural covariance that works well in $\mathbb{R}^2$:

The exponential correlation function is a natural correlation in one dimension, since it corresponds to a Markov process. In two dimensions this is no longer so, although the exponential is a common correlation function in geostatistical work. Whittle (1954) determined the correlation corresponding to a stochastic differential equation of Laplace type:

$$ \left[ \left(\frac{\partial}{\partial t_1}\right)^2 + \left(\frac{\partial}{\partial t_2}\right)^2 - \kappa^2 \right] X(t_1, t_2) = \epsilon(t_1 , t_2) $$ where $\epsilon$ is white noise. The corresponding discrete lattice process is a second order autoregression. (Source: Guttorp&Gneiting)

The family of processes included in the SDE associated with the Matérn equation includes the $AR(1)$ Ornstein–Uhlenbeck model of the velocity of a particle undergoing Brownian motion. More generally, you can define a power spectrum for a family of $AR(p)$ processes for every integer $p$ which also have a Matérn family covariance. This is in the appendix of Rasmussen and Williams.

This covariance function is not related to Matérn cluster process.

References

Cressie, Noel, and Christopher K. Wikle. Statistics for spatio-temporal data. John Wiley & Sons, 2015.

Guttorp, Peter, and Tilmann Gneiting. "Studies in the history of probability and statistics XLIX On the Matern correlation family." Biometrika 93.4 (2006): 989-995.

Rasmussen, C. E. and Williams, C. K. I. Gaussian Processes for Machine Learning. the MIT Press, 2006.

Answer 2

I do not know, but I found this question very interesting and here's what I got after a bit of reading on it.

For certain values of $\nu$, the Matérn covariance function can be expressed as a product of an exponential and a polynomial. E.g. for $\nu = 5/2$: $$C_{5/2}(d) = \sigma^2\left(1 + \frac{\sqrt 5 d}{\rho} + \frac{5d^2}{3\rho^2} \right) \exp \left(- \frac{\sqrt 5 d}{\rho}\right)$$ It is then not too surprising that, as $\nu \to \infty$, $C_\nu$ actually converges to the Gaussian RBF: $$\lim_{\nu \to \infty} C_\nu(d) = \sigma^2 \exp \left( -\frac{d^2}{2\rho^2}\right)$$ For $\nu = 1/2$, the Matérn covariance function gives the absolute exponential kernel $$C_{1/2}(d) = \sigma^2 \exp\left( -\frac{d}{\rho} \right)$$

Furthermore, a Gaussian process with the Matérn covariance function with parameter $\nu$ is $\lceil \nu \rceil -1$-time differentiable .

This is quite nicely demonstrated on a picture taken from Rasmussen & Williams (2006)

In Interpolation of Spatial Data, Stein (who actually proposed the name of the Matérn covariance function), argues (pg. 30) that the infinite differentiability of the Gaussian covariance function yields unrealistic results for physical processes, since observing only a small continuous fraction of space/time should, in theory, yield the whole function. He thus proposed the Matérn version as a generalization that is able to match physical processes more realistically.

Summary

The Matérn covariance function can be seen as a generalization of the Gaussian radial basis function. It contains even the absolute exponential kernel, which gives radically different results, and is better able to capture physical processes due to its finite differentiability (for finite $\nu$).

As for the mysteriousness of the appearance of the Bessel function, I'd love to see further intuition behind that, but I would guess that it is precisely its (asymptotic) behaviour in $\nu$ that made it useful in this context and lead Stein to define the Matérn covariance function. That of course does not rule out the possibility that there's a beautiful argument as to why all of that is true.

Answer 3

There is one aspect of Matérn covariance functions that makes them very useful for physical systems:

It describes an electrical signal with white Gaussian noise passing through an RC low-pass filter. The output signal is time-correlated according to the Matérn covariance function $\nu= 1/2$. When this output signal passes a second low-pass filter, the new output is the Matérn covariance function $\nu=3/2$.

In general, a series of $n$ low-pass filters on white Gaussian noise has the effect of correlating it according to the Matérn function $\nu=(2n-1)/2$.

In physical systems, one often finds influences according to an exponential decay due to one or more independent physical mechanisms, leading to the Matérn covariance functions.