FEM Matérn methods: SPDE, B-splines, and covariance comparison

Patrick Brown

2026-10-11

library(adlaplaceFem)
library(Matrix)

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\):

α ν SPDE Precision operator
2 1 \(\tau(\kappa^2-\Delta)u=\mathcal{W}\) \(\tau^2(\kappa^2-\Delta)^2\)
3 2 \(\tau(\kappa^2-\Delta)^{3/2}u=\mathcal{W}\) \(\tau^2(\kappa^2-\Delta)^3\)

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\).

Quantity Meaning Units (loaloa)
Coordinates / knots Spatial location m
range \(=\rho=\sqrt{8\nu}/\kappa\) Practical range (≈0.1 correlation) m
sd \(=\sigma\) Marginal SD of \(u(\mathbf{s})\) same as \(u\) (e.g. logit)
\(\tau\) Internal SPDE precision scale derived from \(\kappa,\sigma\)

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

Matrix Entries
Mass \(\mathbf{C}\) \(C_{ij}=\langle\psi_i,\psi_j\rangle\)
Stiffness \(\mathbf{G}\) \(G_{ij}=\langle\nabla\psi_i,\nabla\psi_j\rangle\)
Biharmonic \(\mathbf{G}_2\) \((G_2)_{ij}=\langle\Delta\psi_i,\Delta\psi_j\rangle\)
Third-order \(\mathbf{G}_3\) \((G_3)_{ij}=\langle\nabla\Delta\psi_i,\nabla\Delta\psi_j\rangle\)

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.

Derivatives

The Gaussian prior on the weights \(\mathbf{w}\) is

\[ \log\pi(\mathbf{w}\mid\theta) = \frac12\log\det\mathbf{Q}(\theta) -\frac12\mathbf{w}^\top\mathbf{Q}(\theta)\mathbf{w} -\frac{n}{2}\log(2\pi), \]

with \(\theta=(\log\mathrm{range},\log\mathrm{sd})\) and \(n\) the number of weights. The two terms are separate AD groups. Both see the same sparse \(\mathbf{Q}(\theta)=\sum_j c_j(\theta)\,\mathbf{M}_j\), where \(\mathbf{M}_j\) runs over \((\mathbf{C},\mathbf{G},\mathbf{G}_2)\) for \(\alpha=2\) and also \(\mathbf{G}_3\) for \(\alpha=3\). The scalars \(c_j(\theta)\) are the coefficients in the formulas for \(\mathbf{Q}_2\) and \(\mathbf{Q}_3\) above, after substituting \(\kappa\) and \(\tau\) from \((\mathrm{range},\mathrm{sd})\).

The Laplace approximation differentiates this density in three places. The inner optimization needs the value, the gradient, and the Hessian in \(\mathbf{w}\). The outer optimization needs first and second derivatives in \(\theta\). The profile trace needs the contraction of \(H^{-1}\) with third derivatives of the joint density; for this prior those third derivatives are nonzero only in \(\theta\).

Quadratic form

random_fem_ssq_2 and random_fem_ssq_3 contribute \(f=-\frac12\mathbf{w}^\top\mathbf{Q}\mathbf{w}\). \(\mathbf{Q}\) is linear in the \(c_j\) and \(f\) is exactly quadratic in \(\mathbf{w}\), so the derivatives are

\[ \begin{aligned} \frac{\partial f}{\partial\mathbf{w}} &= -\mathbf{Q}\mathbf{w},\\ \frac{\partial f}{\partial\theta_k} &= -\frac12\sum_j\frac{\partial c_j}{\partial\theta_k}\, (\mathbf{w}^\top\mathbf{M}_j\mathbf{w}),\\ \frac{\partial^2 f}{\partial\mathbf{w}\,\partial\mathbf{w}^\top} &= -\mathbf{Q},\\ \frac{\partial^2 f}{\partial\mathbf{w}\,\partial\theta_k} &= -\frac{\partial\mathbf{Q}}{\partial\theta_k}\mathbf{w},\\ \frac{\partial^2 f}{\partial\theta_a\,\partial\theta_b} &= -\frac12\sum_j \frac{\partial^2 c_j}{\partial\theta_a\,\partial\theta_b}\, (\mathbf{w}^\top\mathbf{M}_j\mathbf{w}). \end{aligned} \]

Every one of these is a pass over the stored upper triangle of the Grams: combine the \(c_j\) (or their \(\theta\)-derivatives) into \(\mathbf{Q}\) or \(\partial\mathbf{Q}/\partial\theta_k\), then one quadratic form or one sparse matrix-vector product. A two-input tape maps \(\theta\) to the \(c_j\) and supplies \(\partial c_j/\partial\theta\) and \(\partial^2 c_j/\partial\theta^2\). The products along the stored pattern of \(\mathbf{Q}\) are not operations on that tape.

Third derivatives in \(\mathbf{w}\) are identically zero. In a direction \(\mathbf{v}\) in the weight space the only surviving third-order term is

\[ \frac{\partial}{\partial\theta_k} \Bigl(\mathbf{v}^\top\frac{\partial^2 f}{\partial\mathbf{w}\,\partial\mathbf{w}^\top}\mathbf{v}\Bigr) = -\frac12\mathbf{v}^\top\frac{\partial\mathbf{Q}}{\partial\theta_k}\mathbf{v}. \]

Log-determinant atomic

random_fem_det_2 and random_fem_det_3 contribute \(\frac12\log\det\mathbf{Q}-\frac{n}{2}\log(2\pi)\). This term does not depend on \(\mathbf{w}\), so its inner gradient and Hessian are zero. All of its derivatives are in \(\theta\), obtained by differentiating in the coefficients \(x_j=c_j(\theta)\) and chaining back through the same two-input map.

The atomic fem_logdet accepts the coefficient vector \(\mathbf{x}\) (length 3 or 4) and returns \(\log\det\mathbf{Q}(\mathbf{x})\) for \(\mathbf{Q}=\sum_j x_j\mathbf{M}_j\). The sparse LDL and the selected inverse stay inside the atomic; the tape records only the call and the chain rule from \(\mathbf{x}\) to \(\theta\).

One evaluation assembles the numeric values of \(\mathbf{Q}\), factors it with the symbolic LDL pattern stored on the payload, and, when a derivative is required, forms the Takahashi selected inverse \(\boldsymbol{\Sigma}\) (\(\mathbf{Q}^{-1}\) on the sparsity pattern of the factor). The identity \(d\log\det\mathbf{Q}=\operatorname{Tr}(\mathbf{Q}^{-1}\,d\mathbf{Q})\) then gives

\[ \frac{\partial\log\det\mathbf{Q}}{\partial x_j} = \operatorname{Tr}(\mathbf{Q}^{-1}\mathbf{M}_j) = \sum_k \omega_k (M_j)_k \Sigma_k, \]

where the sum is over stored upper-triangle entries and \(\omega_k\) is \(1\) on the diagonal and \(2\) off it.

That selected inverse is takahashi_davis. For each column of \(\mathbf{L}\), from last to first, the column of \(\boldsymbol{\Sigma}\) is scattered into a work vector of length \(n\), and the Takahashi update reads the vector by row index. The stored result is \(\mathbf{Q}^{-1}\) on the sparsity pattern of \(\mathbf{L}\). fem_logdet calls this function. A dual scalar in the work vector is what produces the directional derivative below.

The same subset can be computed without a FEM model. From R, adlaplace::takahashi_davis(Q) returns a dgCMatrix of those entries in the original variable order. Pass chol = fem_chol_pattern(Q) when that symbolic factor is already known. From C++, a package with LinkingTo: adlaplace includes adlaplace/takahashi_impl.hpp. adlaplace::chol::takahashi_davis takes \(\mathbf{L}\) and \(\mathbf{D}\). A caller that has only \(\mathbf{Q}\) factors it with chol_update_csc first.

Because \(\mathbf{Q}\) is linear in \(\mathbf{x}\), the second derivative has no \(\operatorname{Tr}(\mathbf{Q}^{-1}\partial^2\mathbf{Q})\) term:

\[ \frac{\partial^2\log\det\mathbf{Q}}{\partial x_i\,\partial x_j} = -\operatorname{Tr}(\mathbf{Q}^{-1}\mathbf{M}_i\mathbf{Q}^{-1}\mathbf{M}_j). \]

CppAD asks the atomic for this as a Hessian-vector product. The same LDL and Takahashi code runs with a dual scalar \((a,\dot a)\), so one factorization returns both \(\partial\log\det\mathbf{Q}/\partial x_j\) and its directional derivative. The prior density uses half of these derivatives.

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. \]

Hierarchical B-splines

For spatial domains where data are concentrated in subregions, build hierarchical knot lines with [hb_knots()] and pass the result to matern():

kn <- hb_knots(
  outer = list(                          # level 0: padded domain
    xmin = ..., xmax = ..., ymin = ..., ymax = ...,
    resolution = c(dx, dy)
  ),
  inner = list(
    c(xmin, xmax, ymin, ymax),           # level 1: study region
    list(west_ext, east_ext)             # level 2: sibling rectangles
  ),
  fact = c(2, 4)
)
matern(geometry, knots = kn, shape = 2L)

fact[i] is the cell-size ratio from the previous level. Each inner extent is snapped to that level’s nested cell grid. This gives nested spaces \(V^0 \subset V^1 \subset \cdots \subset V^L\). A truncated hierarchical (THB) basis keeps level-\(\ell\) tensor B-splines active when their support lies in the level-\(\ell\) region but not in finer regions. Every active function embeds in the finest space via a sparse map \(\mathbf{S}\), so the Gram matrices become

\[ \mathbf{C}_H = \mathbf{S}^\top \mathbf{C}_L \mathbf{S},\quad \mathbf{G}_H = \mathbf{S}^\top \mathbf{G}_L \mathbf{S},\ \ldots \]

with \(\mathbf{C}_L,\mathbf{G}_L,\ldots\) the usual Kronecker Grams on the finest level. The SPDE precision formulas and C++ random_fem_* kernels are unchanged. Inspect the layout with hb_summary(hb_basis(kn, degree)).

Cost scales with the finest global knot grid; local refinement adds resolution where needed without refining the entire padded domain.

Numerical comparison: Matérn vs FEM

Fix \(\nu=1\) (SPDE order \(\alpha=2\)), \(\kappa=3\), and marginal \(\mathrm{sd}=1\). The practical range is \(\rho=\sqrt{8\nu}/\kappa\). The analytic covariance below uses the Matérn kernel at that \((\kappa,\mathrm{sd})\). The FEM covariances use the same \((\rho,\mathrm{sd})\) through [matern_fem_covariance()], which returns \(\mathbf{A}\mathbf{Q}^{-1}\mathbf{A}^\top\) at the geometry in data.

nu <- 1
kappa <- 3
sd <- 1
rho <- sqrt(8 * nu) / kappa
c(kappa = kappa, range = rho, sd = sd)
##    kappa    range       sd 
## 3.000000 0.942809 1.000000

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 <- data.frame(
  x = runif(100),
  y = runif(100)
)
sites$geometry <- as.matrix(sites[c("x", "y")])

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[c("x", "y")]))
Sigma <- matern_cov(D, nu, kappa, sigma2 = sd^2)
diag(Sigma) <- sd^2

FEM covariance at the same sites, for two knot resolutions:

term_coarse <- matern("geometry", knots = knots_coarse, shape = 1L)
term_fine <- matern("geometry", knots = knots_fine, shape = 1L)
Sig_coarse <- matern_fem_covariance(term_coarse, sites, range = rho, sd = sd)
Sig_fine <- matern_fem_covariance(term_fine, sites, range = rho, sd = sd)

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).