SPDE Matérn Model

Overview

The SPDE approach represents a continuous spatial random field with Matérn covariance as the solution of a stochastic partial differential equation (SPDE). Solving the SPDE on a triangular mesh by finite elements produces a sparse Gaussian Markov random field approximation whose precision matrix is small-banded and cheap to factorize. This makes it possible to fit smooth Gaussian processes on tens or hundreds of thousands of locations with INLA in seconds rather than hours.

• Continuous space: for point-referenced (geostatistical) data: temperature, pollution, species counts at irregular locations.
• Matérn covariance: controlled by two parameters, the practical range ρ (distance at which the correlation drops to roughly 0.1; about 0.139 at the default smoothness ν=1) and the marginal standard deviation σ.
• PC priors: Penalised-Complexity priors on (ρ, σ) parameterised in interpretable units, "P(ρ < r₀) = p" and "P(σ > s₀) = p".
• Mesh-based: you build a triangular mesh once and project both observations and predictions through a sparse projector matrix A.

Mathematical Formulation

Continuous Matérn Field

The latent field is a stationary, isotropic Gaussian process \(u(\mathbf{s})\) with Matérn covariance

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

where \(K_\nu\) is the modified Bessel function of the second kind, \(\nu > 0\) is the smoothness, and \(\kappa > 0\) controls the correlation range.

SPDE Representation

Lindgren, Rue & Lindström (2011) showed that this Matérn field is the stationary solution of the SPDE

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

where \(\mathcal{W}\) is Gaussian white noise on \(\mathbb{R}^d\) and the SPDE-order parameter is

\[\alpha = \nu + d/2.\]

In two dimensions (\(d=2\)) the default alpha=2 gives smoothness \(\nu = 1\). The marginal variance in this case is

\[\sigma^2 = \frac{1}{4\pi\,\kappa^2\,\tau^2}.\]

Practical Range

The SPDE inverse-length scale \(\kappa\) is unintuitive. The doc parameterises everything in terms of the practical range, the distance at which the correlation has dropped to roughly 0.1:

\[\rho = \frac{\sqrt{8\nu}}{\kappa}.\]

The exact correlation at distance \(\rho\) drifts with the smoothness \(\nu\); at the default \(\nu = 1\) it is approximately 0.139, with \(\rho = 2\sqrt{2}/\kappa \approx 2.83/\kappa\). The PC prior is placed directly on \(\rho\) so the user reasons in distance units.

Finite-Element Approximation

The continuous SPDE is solved on a triangular mesh \(\mathcal{T}\) with \(n\) vertices using piecewise-linear basis functions \(\{\psi_i\}_{i=1}^n\). The latent field is approximated by

\[u(\mathbf{s}) \approx \sum_{i=1}^{n} \psi_i(\mathbf{s})\, u_i, \qquad \mathbf{u} = (u_1,\ldots,u_n)^\top.\]

The mesh-vertex weights \(\mathbf{u}\) form a sparse Gaussian Markov random field. With \(\alpha = 2\) the precision matrix has the closed form

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

where \(\mathbf{C}\) is the (lumped) mass matrix and \(\mathbf{G}\) is the stiffness matrix produced by the finite-element discretization. For this \(\alpha = 2\) case, the \(\mathbf{G}\,\mathbf{C}^{-1}\,\mathbf{G}\) term gives \(\mathbf{Q}\) the sparsity pattern of the mesh adjacency squared (two-hop neighbourhood); for \(\alpha = 1\) it would be one-hop only.

Projector Matrix

Observation locations \(\mathbf{s}_1, \ldots, \mathbf{s}_m\) need not coincide with mesh vertices. The sparse projector \(\mathbf{A} \in \mathbb{R}^{m\times n}\) maps mesh weights to observation values via barycentric interpolation: each row of \(\mathbf{A}\) has at most three non-zeros (the barycentric weights of the triangle containing the observation). The model is then

\[y_j \mid \boldsymbol{\eta} \sim \pi(y_j;\, \eta_j), \qquad \boldsymbol{\eta} = \mathbf{X}\boldsymbol{\beta} + \mathbf{A}\,\mathbf{u}.\]

Points outside the mesh have an all-zero row in \(\mathbf{A}\) (they contribute the zero spatial effect).

When to Use SPDE

The SPDE Matérn model is the right choice for geostatistical data: observations attached to continuous spatial coordinates (latitude/longitude, x/y in a planar map) where you expect smooth spatial variation that decays with distance.

Specification in pyINLA

An SPDE random effect is specified by passing an SPDE object in the spec's 'model' slot, plus a sparse projector matrix in 'A.local'. This shape differs from the string-keyed specs of iid/rw1/besag: the safety layer recognises an SPDE object via its class and only allows the three keys below.

Allowed Parameters

No other keys are accepted (no hyper, constr, diagonal, etc.). Hyperparameters and priors are configured at the SPDE-object construction step instead, via spde2_pcmatern's prior_range, prior_sigma, and alpha arguments.

Per-Construction-Function Runnable Examples

Most of the SPDE API surface is in the construction functions you call before pyinla(...). Click any header below to expand a minimal runnable snippet. All snippets assume the shared setup block.

# shared setup for every snippet below
import numpy as np
import pandas as pd
from pyinla import pyinla
from pyinla.fmesher import fm_mesh_2d, spde2_pcmatern, spde_make_A

rng = np.random.default_rng(42)
nobs = 200
coords = rng.uniform(0, 10, size=(nobs, 2))   # observation locations on a 10x10 square
df = pd.DataFrame({
    'y': rng.standard_normal(nobs),
    'x': rng.standard_normal(nobs),
})

# Two key knobs:
#   max_edge = (inner, outer): max triangle side length, inside the data
#              region and in the outer "buffer" zone.
#   offset   = (inner, outer): how far the buffer extends beyond the data.
# Smaller max_edge = more vertices = finer field but slower fit.
# cutoff merges any input points closer than this distance.
mesh = fm_mesh_2d(
    loc      = coords,
    max_edge = [0.5, 2.0],
    offset   = [0.5, 2.0],
    cutoff   = 0.1,
)
print(f"Mesh vertices: {mesh.n}, triangles: {mesh.n_triangle}")

# prior_range = (r0, p): P(range < r0) = p
# prior_sigma = (s0, p): P(sigma > s0) = p
# alpha=2 (default) gives Matern smoothness nu = 1 in 2D.
spde = spde2_pcmatern(
    mesh        = mesh,
    prior_range = (1.0, 0.5),     # 50% probability that range < 1
    prior_sigma = (1.0, 0.5),     # 50% probability that sigma > 1
    alpha       = 2,
)
print(spde.summary())

# Sparse matrix of shape (n_obs, n_vertices). Each row has at most 3
# non-zeros (the barycentric weights of the containing triangle).
# Points outside the mesh get an all-zero row (zero spatial effect).
A = spde_make_A(mesh, loc=coords)
print(f"A shape: {A.shape}, nnz: {A.nnz}")

# Pass the SPDE object (not a string) as 'model', and the sparse
# projector matrix as 'A.local'. The id is the latent-component name.
formula = {
    'response': 'y',
    'fixed':    ['1', 'x'],
    'random':   [{
        'id':      'spatial',
        'model':   spde,
        'A.local': A,
    }],
}
result = pyinla(formula=formula, family='gaussian', data=df)

# In 2D, smoothness nu = alpha - d/2 = alpha - 1.
# Currently only alpha = 2 (nu = 1) is supported by spde2_pcmatern in
# pyinla; other values raise NotImplementedError ("B matrices not
# implemented for alpha=...") because the FEM B-matrix derivation is
# only wired up for that case. alpha=2 produces a smooth Matern field
# and is the standard default in the SPDE literature.
spde_smooth = spde2_pcmatern(mesh, prior_range=(1.0, 0.5),
                                prior_sigma=(1.0, 0.5), alpha=2)

# For comparison, the (currently unsupported) alpha=1 call would be:
# spde_rough = spde2_pcmatern(mesh, prior_range=(1.0, 0.5),
#                             prior_sigma=(1.0, 0.5), alpha=1)
# -> NotImplementedError: B matrices not implemented for alpha=1.0, d=2.

# Tuple (r0, p) interpreted as P(range < r0) = p.
# Tighter prior: smaller p (e.g. 0.05) puts more weight on "range >= r0".
# Choose r0 in the units of your coordinates, around what you expect.
spde = spde2_pcmatern(
    mesh,
    prior_range = (2.0, 0.05),     # P(range < 2) = 0.05  -> "range likely >= 2"
    prior_sigma = (1.0, 0.5),
)

# Tuple (s0, p) interpreted as P(sigma > s0) = p.
# Tighter prior: smaller p shrinks sigma toward 0 more strongly.
# Choose s0 around the magnitude of variation you expect in the response.
spde = spde2_pcmatern(
    mesh,
    prior_range = (1.0, 0.5),
    prior_sigma = (0.5, 0.05),     # P(sigma > 0.5) = 0.05  -> "sigma likely small"
)

from pyinla.fmesher import fm_segm

# Define a closed boundary as a sequence of (x, y) points.
# fm_segm(is_bnd=True) marks the segment as a domain boundary.
bnd_pts = np.array([
    [0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0],
])
boundary = fm_segm(loc=bnd_pts, is_bnd=True)
mesh = fm_mesh_2d(
    loc      = coords,
    boundary = boundary,
    max_edge = [0.5, 2.0],
)

from pyinla.fmesher import fm_hexagon_lattice

# Generate a regular hexagonal lattice of points inside a polygon.
# edge_len is the hex side length (in coordinate units). The returned
# points are ideal as seed `loc` for fm_mesh_2d on irregular regions.
poly_coords = np.array([
    [0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0],
])
hex_pts = fm_hexagon_lattice(poly_coords, edge_len=1.0)
print(f"Hex lattice: {len(hex_pts)} points")

# Build the mesh from the hex seeds (more uniform than from raw observations).
mesh = fm_mesh_2d(loc=hex_pts, max_edge=[0.5, 2.0], offset=[0.5, 2.0])

from pyinla.fmesher import fm_nonconvex_hull

# Concave hull around scattered points. Use as a boundary when the data
# lie on an irregular region (e.g. coastlines) and a convex hull would
# over-cover empty space. Requires the R `sf` package.
#   convex  ~ tightness of the hull (more negative = tighter)
#   concave ~ minimum gap to consider a concavity
boundary = fm_nonconvex_hull(
    coords,
    convex  = -0.15,
    concave = -0.15,
)
mesh = fm_mesh_2d(
    loc      = coords,
    boundary = boundary,
    max_edge = [0.5, 2.0],
)

from pyinla.fmesher import spde_grid_projector

# Build a sparse projector from mesh vertices onto a regular dims[0] x dims[1]
# grid covering xlim x ylim. Use it to draw posterior-mean / sd maps of the
# spatial effect, or to set up a prediction grid (see Typical workflow below).
grid = spde_grid_projector(
    mesh = mesh,
    xlim = (0.0, 10.0),
    ylim = (0.0, 10.0),
    dims = (200, 200),
)
print(f"Grid: {grid.dims} = {grid.n_grid} points; A shape {grid.proj.shape}")

# grid.proj is a sparse (n_grid, n_vertices) matrix. To visualise the
# posterior-mean spatial effect after fitting:
#   spatial_mean = result.summary_random['spatial']['mean'].values
#   z = grid.project(spatial_mean).reshape(grid.dims[1], grid.dims[0])
#   plt.imshow(z, origin='lower', extent=[*grid.xlim, *grid.ylim])

Hyperparameters

The SPDE2-PC-Matérn model has two free hyperparameters:

1. Practical range ρ, internally parameterised as \(\theta_1 = \log(\rho)\).
2. Marginal standard deviation σ, internally parameterised as \(\theta_2 = \log(\sigma)\).

The internal SPDE coefficients \(\kappa\) and \(\tau\) are derived from \((\rho, \sigma)\) and the smoothness \(\nu\) by

\[\kappa = \frac{\sqrt{8\nu}}{\rho}, \qquad \tau = \sqrt{\frac{\Gamma(\nu)}{\Gamma(\alpha)\,(4\pi)^{d/2}\,\kappa^{2\nu}\,\sigma^2}}.\]

PC Prior Defaults

Both priors are Penalised-Complexity priors (Simpson et al., 2017), which shrink the field toward the simpler "no spatial effect" base model unless the data argue otherwise.

Notes

• Mesh quality drives accuracy and cost. Smaller max_edge → finer mesh → better SPDE approximation but larger \(\mathbf{Q}\) and slower fit. A common rule of thumb is max_edge <= prior_range[0] / 5 in the inner region.
• Outer buffer matters. The Matérn covariance produces edge artefacts near the domain boundary. offset controls how far the mesh extends beyond your data; the buffer should be roughly one practical range.
• Coordinate units. prior_range uses the same units as your coordinates. If you pass lat/lon, ranges are in degrees; for kilometres or metres, project first.
• Only alpha = 2 (smoothness \(\nu = 1\)) is supported by spde2_pcmatern in pyinla today. Any other alpha raises NotImplementedError from SPDE2PcMatern._compute_b_matrices at construction time ("B matrices not implemented for alpha=…"). The internal precision builder has stub branches for alpha = 1, but they're currently unreachable because the B-matrix derivation is wired up only for alpha = 2.
• Points outside the mesh get an all-zero row in A from spde_make_A. This means they contribute the zero spatial effect (not NaN); make sure your mesh covers all observation and prediction locations or the predictions will be flat.
• Constraints. SPDE is a proper (non-intrinsic) prior, so no sum-to-zero constraint is added. The intercept in fixed remains identifiable.