Continuous model: Matérn via Whittle SPDE

A zero-mean Gaussian field \(u(\mathbf{s})\) on \(\mathbb{R}^d\) with Matérn covariance

\[ \operatorname{Cov}\bigl(u(\mathbf{s}),u(\mathbf{s}')\bigr) = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)} (\kappa\|\mathbf{h}\|)^\nu K_\nu(\kappa\|\mathbf{h}\|), \quad \mathbf{h}=\mathbf{s}-\mathbf{s}', \]

is a stationary solution of

\[ \tau\,(\kappa^2-\Delta)^{\alpha/2}\,u(\mathbf{s})=\mathcal{W}(\mathbf{s}), \qquad \alpha=\nu+d/2, \]

with Gaussian white noise \(\mathcal{W}\) (Whittle; Lindgren, Rue & Lindström 2011). The continuous precision operator is \(\tau^2(\kappa^2-\Delta)^\alpha\).

In \(d=2\):

On a bounded domain, homogeneous Neumann boundary conditions are used; pad the domain by about half a correlation range in applications.

Parameterization: range and sd

The FEM term exposes hyperparameters \((\mathrm{range},\mathrm{sd})\) (optimized as \(\log\mathrm{range}\), \(\log\mathrm{sd}\)), not the raw SPDE pair \((\kappa,\tau)\). The public range is the same practical range as geostatsp/INLA, \(\rho=\sqrt{8\nu}/\kappa\). Internally the AD kernels convert

\[ \kappa = \sqrt{8\nu}/\mathrm{range}, \qquad \tau = \begin{cases} 1\big/\bigl(\kappa\,\sigma\sqrt{4\pi}\bigr) & \alpha=2\ (\nu=1),\\ 1\big/\bigl(\kappa^2\,\sigma\sqrt{8\pi}\bigr) & \alpha=3\ (\nu=2), \end{cases} \]

with \(\sigma=\mathrm{sd}\), so that \(\mathrm{sd}\) is the marginal standard deviation of the Matérn field \(u(\mathbf{s})\) on the linear-predictor scale. In 2D, \(\nu=\alpha-1\).

So with coordinates in metres and \(\nu=2\), range = 1e5 means a practical range of \(100\,\mathrm{km}\) and \(\kappa=4\times10^{-5}\,\mathrm{m}^{-1}\).

geostatsp / INLA glgm uses the same (range, sd) meaning, so coefficient tables can be compared directly:

cf <- coef(fit)
c(
  range = unname(cf[["geometry_matern_log_range"]]),
  sd = unname(cf[["geometry_matern_log_sd"]])
)

(glgm may label the range row range/1000 when the summary is in kilometres; convert by \(\times 1000\) first.)

The Gram matrices \(\mathbf{C}\), \(\mathbf{G}\), … use the same coordinate units as the knots, so range is already on that scale.

Finite-element weights and Grams

Write \(\widetilde{u}(\mathbf{s})=\boldsymbol{\psi}(\mathbf{s})^\top\mathbf{w}\). The weight precision is assembled from

Working formulas (higher-order B-splines; no mass lumping):

\[ \mathbf{Q}_2=\tau^2\bigl(\kappa^4\mathbf{C}+2\kappa^2\mathbf{G}+\mathbf{G}_2\bigr), \]

\[ \mathbf{Q}_3=\tau^2\bigl(\kappa^6\mathbf{C}+3\kappa^4\mathbf{G}+3\kappa^2\mathbf{G}_2+\mathbf{G}_3\bigr). \]

These replace Galerkin expressions involving \(\mathbf{C}^{-1}\), which destroy sparsity. Recommended B-spline degrees: \(p=2\) for α = 2; \(p=3\) for α = 3.

Tensor-product B-splines

On a rectangle, with 1-D B-splines \(B_i^x\) and \(B_j^y\), \(\psi_{ij}(x,y)=B_i^x(x)\,B_j^y(y)\) and the Grams factor as Kronecker sums

\[ \mathbf{C}=\mathbf{C}_x\otimes\mathbf{C}_y, \]

\[ \mathbf{G}=\mathbf{G}_x\otimes\mathbf{C}_y+\mathbf{C}_x\otimes\mathbf{G}_y, \]

\[ \mathbf{G}_2=\mathbf{G}_{2,x}\otimes\mathbf{C}_y+\mathbf{C}_x\otimes\mathbf{G}_{2,y}+2\,\mathbf{G}_x\otimes\mathbf{G}_y. \]

Numerical comparison: Matérn vs FEM

Fix α = 2, \(\nu=1\), \(\kappa=2\), \(\tau\) chosen so the Matérn marginal variance is 1:

nu <- 1
kappa <- 3
# sigma^2 = Gamma(nu) / ( Gamma(alpha) (4*pi)^{d/2} kappa^{2*nu} tau^2 ) with alpha=2, d=2
# => 1 = 1 / (1 * 4*pi * kappa^2 * tau^2)  => tau^2 = 1/(4*pi*kappa^2)
tau <- 1 / sqrt(4 * pi * kappa^2)
rho <- sqrt(8 * nu) / kappa
c(kappa = kappa, tau = tau, range = rho)

##     kappa       tau     range 
## 3.0000000 0.0940316 0.9428090

Evaluation sites (same for both grids):

knots_coarse <- list(
  x = seq(-0.2, 1.2, by = 0.2),
  y = seq(-0.2, 1.2, by = 0.2)
)
knots_fine <- list(
  x = seq(-1, 2, by = 0.1),
  y = seq(-1, 2, by = 0.1)
)

set.seed(0)
sites_eval <- data.frame(
  x = runif(100),
  y = runif(100)
)


matern_cov <- function(h, nu, kappa, sigma2 = 1) {
  h <- pmax(h, 1e-12)
  r <- kappa * h
  sigma2 * (2^(1 - nu) / gamma(nu)) * (r^nu) * besselK(r, nu)
}

D <- as.matrix(dist(sites_eval))
Sigma <- matern_cov(D, nu, kappa, sigma2 = 1)
diag(Sigma) <- 1

FEM covariance at the same sites: \(\mathbf{A}\mathbf{Q}^{-1}\mathbf{A}^\top\) for two knot resolutions.

fem_cov_at_sites <- function(sites_eval, knots_list) {
  fem <- grf_bspline(
    sites_eval,
    knots_list,
    degree = 2L
  )
  Q <- fem_precision(kappa, tau, fem$C, fem$G, fem$G2, alpha = 2L)
  A <- fem$A
  X <- Matrix::solve(Q, Matrix::t(A))
  as.matrix(A %*% X)
}

Sig_coarse <- fem_cov_at_sites(sites_eval, knots_coarse)
Sig_fine <- fem_cov_at_sites(sites_eval, knots_fine)

rel_err <- function(S_hat, S) {
  sqrt(sum((S_hat - S)^2)) / sqrt(sum(S^2))
}

data.frame(
  grid = c("coarse (by = 0.2)", "fine (by = 0.1)"),
  frobenius_rel = c(rel_err(Sig_coarse, Sigma), rel_err(Sig_fine, Sigma)),
  max_abs = c(max(abs(Sig_coarse - Sigma)), max(abs(Sig_fine - Sigma))),
  mean_diag = c(mean(diag(Sig_coarse)), mean(diag(Sig_fine)))
)

##                grid frobenius_rel    max_abs mean_diag
## 1 coarse (by = 0.2)    0.62707120 1.16299161 1.4525017
## 2   fine (by = 0.1)    0.01084377 0.02165488 0.9876545

Correlation versus distance (exact Matérn vs FEM):

upper <- upper.tri(D)
df <- data.frame(
  dist = D[upper],
  matern = Sigma[upper] / sqrt(diag(Sigma)[row(Sigma)[upper]] * diag(Sigma)[col(Sigma)[upper]]),
  coarse = Sig_coarse[upper] / sqrt(pmax(diag(Sig_coarse)[row(Sig_coarse)[upper]], 1e-8) *
    pmax(diag(Sig_coarse)[col(Sig_coarse)[upper]], 1e-8)),
  fine = Sig_fine[upper] / sqrt(pmax(diag(Sig_fine)[row(Sig_fine)[upper]], 1e-8) *
    pmax(diag(Sig_fine)[col(Sig_fine)[upper]], 1e-8))
)
df <- df[order(df$dist), ]
plot(df$dist, df$matern,
  type = "l", lwd = 2,
  ylim = c(0, 1),
  xlim = c(0, 1),
  xlab = "distance", ylab = "correlation",
  main = "Matern vs FEM correlations"
)
lines(df$dist, df$coarse, col = "tomato", lty = 2, lwd = 4)
lines(df$dist, df$fine, col = "steelblue", lty = 3, lwd = 2)
legend("topright", c("Matern", "FEM coarse", "FEM fine"),
  col = c("black", "tomato", "steelblue"), lty = c(1, 2, 3), lwd = c(2, 1, 1), bty = "n"
)

The finer knot grid tracks the Matérn correlation more closely (especially away from the boundary).

References

Whittle (1954, 1963); Lindgren, Rue & Lindström (2011); Bolin & Lindgren (CSDA 2013); Liu, Guillas & Lai (2016); Lindgren, Bolin & Rue (2022); Bakka (arXiv:1803.03765).