Understanding SPDEs

The Matérn Covariance

The Matérn covariance is the most widely used covariance function in spatial statistics. It describes how the similarity between observations decreases with distance.

Correlation Function

For two locations \(\mathbf{s}_i\) and \(\mathbf{s}_j\), the Matérn correlation is:

Covariance Function

The covariance is the correlation scaled by the marginal variance \(\sigma^2\):

The Parameters

When distance is zero, the Bessel function term equals 1, so \(\text{Cor} = 1\) and \(\text{Cov} = \sigma^2\).

Understanding Smoothness (\(\nu\))

The smoothness parameter has a direct physical interpretation:

• \(\nu = 0.5\): Exponential covariance. Rough, non-differentiable surfaces (like Brownian motion).
• \(\nu = 1\): Once differentiable. Moderately smooth surfaces. Default in INLA.
• \(\nu = 1.5\): 1.5 times differentiable. Smoother surfaces.
• \(\nu = 2\): Twice differentiable. Very smooth surfaces.
• \(\nu \to \infty\): Gaussian (squared exponential) covariance. Infinitely smooth.

Practical Range

The scale parameter \(\kappa\) is hard to interpret directly. Instead, we use the practical range, which is the distance at which correlation drops to approximately 0.14 (i.e., \(e^{-1} \approx 0.13\)):

Interpretation: If the practical range is 50 km, then two locations 50 km apart have correlation ~0.14. Beyond this distance, observations are nearly uncorrelated.

Marginal Standard Deviation (\(\sigma\))

The marginal standard deviation \(\sigma\) controls how much the spatial field varies. For the SPDE parameterization:

where \(\tau\) is an internal precision parameter. In 2D with \(\nu=1\), this simplifies to \(\sigma = 1/(\tau \sqrt{4\pi} \kappa)\).

Interpretation: If \(\sigma = 2\), then at any location, the spatial effect has standard deviation 2 (in the units of your response variable).

The SPDE-Matérn Link

Why This Matters

Computing Matérn covariance directly requires building and inverting large dense matrices, which becomes computationally infeasible for thousands of locations. The SPDE approach provides a shortcut: instead of working with the covariance matrix, we work with a sparse precision matrix that is much faster to compute.

The Key Insight (No Differential Equations Required)

Think of a Matérn random field as a smooth, spatially correlated surface. Lindgren et al. (2011) discovered that this surface can be created through a simple recipe:

1. Start with random noise: Imagine sprinkling random values across your spatial domain, each location getting an independent random value.
2. Smooth the noise: Apply a smoothing operation that blends nearby values together. The amount of smoothing determines how correlated nearby locations become.
3. The result: After smoothing, you get exactly a Matérn covariance structure.

The Mathematical Statement

For those interested, the formal result is that a Gaussian field with Matérn covariance is the solution to this stochastic partial differential equation (SPDE):

The goal: We want to find \(u(\mathbf{s})\), the spatial random field. This is the unknown we're solving for. The SPDE tells us that \(u(\mathbf{s})\) (the result of smoothing white noise) has exactly Matérn covariance structure.

Understanding Each Part

Connection to Matérn Parameters

The smoothing parameters connect directly to the Matérn covariance parameters:

• \(\alpha = \nu + d/2\) where \(d\) is the spatial dimension (usually 2)
• Higher \(\alpha\) means more smoothing passes, giving smoother fields (higher \(\nu\))
• The scale parameter \(\kappa\) controls the correlation range

Summary: What You Need to Remember

You don't need to understand the differential equation. Just remember:

• The SPDE approach is a computational trick to get Matérn covariance efficiently
• It uses the mesh to represent spatial smoothing with sparse matrices
• This makes fitting models with thousands of locations fast and practical
• The results are statistically equivalent to fitting Matérn covariance directly

From SPDE to FEM: The Bridge

The Problem

The SPDE defines \(u(\mathbf{s})\), the spatial field we want to find. But it's defined continuously over all of space: we'd need to know \(u(\mathbf{s})\) at infinitely many points. That's impossible to compute.

The Solution: Discretize with a Mesh

This is where the triangular mesh comes in. Instead of finding \(u(\mathbf{s})\) everywhere, we approximate it using values only at mesh vertices:

1. Cover the domain with triangles: The mesh divides your spatial region into triangles.
2. Estimate \(u\) at vertices only: We find \(u\) at each of the \(m\) mesh vertices (call these values \(w_1, w_2, \ldots, w_m\)).
3. Interpolate between vertices: Within each triangle, \(u(\mathbf{s})\) is linearly interpolated from the three corner values (like stretching a rubber sheet over the vertices).

Why This Works

The SPDE says "smooth this noise to get \(u(\mathbf{s})\)" continuously. The FEM approximation says "smooth this noise, but only track values at mesh vertices, and interpolate in between." As the mesh gets finer (more vertices), the approximation gets closer to the true continuous solution.

FEM Matrices and Precision

Approximating the Spatial Field with Basis Functions

Remember: our goal is to find \(u(\mathbf{s})\), the spatial field from the SPDE. FEM approximates this continuous field as a weighted sum of basis functions defined on the triangular mesh:

Where:

• \(\psi_k(\mathbf{s})\): Basis functions (one per mesh vertex), defined by the triangles
• \(w_k\): Random weights at each vertex. These are what we estimate.
• \(m\): Number of mesh vertices

So instead of finding \(u(\mathbf{s})\) at every point in space, we only need to estimate the \(m\) weights \(\{w_1, w_2, \ldots, w_m\}\). The basis functions handle interpolation everywhere else.

What Are Tent Functions?

Think of the mesh as lying flat on the ground (the x-y plane). Each basis function \(\psi_k\) adds a third dimension (height) that rises above the mesh:

For each vertex \(k\), the tent function \(\psi_k\) has these properties:

• Height = 1 directly above vertex \(k\) (the peak of the tent)
• Height = 0 at all neighboring vertices (the tent touches the ground there)
• Slopes linearly from height 1 down to height 0 within the triangles touching vertex \(k\)
• Height = 0 everywhere outside those triangles

When you sum all these tent functions (each multiplied by its weight \(w_k\)), you get a piecewise linear surface that passes through the weight values at each vertex.

Plugging Back into the SPDE

Here's the key step: we substitute our approximation \(u(\mathbf{s}) \approx \sum_k \psi_k(\mathbf{s}) w_k\) back into the SPDE equation. This converts the continuous differential equation into a discrete system of equations for the weights \(\{w_1, \ldots, w_m\}\).

The process (called the Galerkin method) involves:

1. Substitute the basis function approximation into the SPDE
2. Multiply by each basis function \(\psi_i\) and integrate over the domain
3. This produces integrals of products of basis functions

The FEM Matrices

From plugging the basis functions into the SPDE, we get two matrices:

Key point: Both C and G are sparse because each basis function \(\psi_k\) is only non-zero near vertex \(k\). Two basis functions only interact if their vertices share a triangle.

The Precision Matrix

The precision matrix \(\mathbf{Q}\) is built from C and G. The exact formula depends on the smoothness parameter α you choose:

Why different formulas? Higher α means more smoothing, which requires a more complex formula. Think of it like applying a blur filter: α = 1 is one pass, α = 2 is two passes (smoother result). The formula for α = 2 is essentially "apply the α = 1 formula twice."

Worked Example: A Single Triangle

Let's work through a complete example with the simplest possible mesh: 3 vertices forming one triangle.

Step 1: The Mesh

Consider a triangle with vertices at:

• Vertex 1: \((0, 0)\)
• Vertex 2: \((1, 0)\)
• Vertex 3: \((0, 1)\)

This is a right triangle with area = 0.5.

Step 2: The Basis Functions

Each vertex gets a tent function. For this triangle:

• \(\psi_1(\mathbf{s})\): equals 1 at vertex 1, linearly decreases to 0 at vertices 2 and 3
• \(\psi_2(\mathbf{s})\): equals 1 at vertex 2, linearly decreases to 0 at vertices 1 and 3
• \(\psi_3(\mathbf{s})\): equals 1 at vertex 3, linearly decreases to 0 at vertices 1 and 2

At any point inside the triangle: \(\psi_1 + \psi_2 + \psi_3 = 1\) (they form a partition of unity).

Step 3: The Spatial Field Approximation

We approximate \(u(\mathbf{s})\) as:

Where \(w_1, w_2, w_3\) are the random weights we estimate. At vertex 1, \(u = w_1\). At vertex 2, \(u = w_2\). Inside the triangle, \(u\) is linearly interpolated.

Step 4: Computing the Mass Matrix C

Now let's compute the actual values of \(C_{ik}\).

Let's compute step by step:

The formula (standard result for linear basis functions on triangles):

Step 5: Computing the Stiffness Matrix G

The stiffness matrix measures how the slopes (gradients) of the tent functions interact.

For our right triangle:

Step 6: Building the Precision Matrix Q

For \(\alpha = 1\) (simplest case), with scale parameter \(\kappa\):

Let's compute Q step by step for κ = 2:

Step 7: What Q Tells Us

The precision matrix \(\mathbf{Q}\) is the inverse of the covariance matrix for the weights \(\mathbf{w} = (w_1, w_2, w_3)\):

Summary of the Full Chain

Deriving the Random Walk 1 (RW1) Precision Matrix

The same FEM/Galerkin approach works in 1D! Let's derive the Random Walk of order 1 precision matrix using the exact same steps we used for the 2D Matérn field.

The 1D Differential Equation

Step 1: The 1D Mesh

Consider n = 4 equally-spaced points on a line with spacing h = 1:

Step 2: The 1D Basis Functions (Hat Functions)

Each point gets a "hat function" - the 1D equivalent of tent functions:

Step 3: The Galerkin Derivation (Full Detail)

Step 4: Computing the Stiffness Matrix G

Deriving the Random Walk 2 (RW2) Precision Matrix

Following the same approach, let's derive the Random Walk of order 2 precision matrix. RW2 produces smoother curves than RW1 by penalizing changes in slope rather than changes in level.

The Galerkin Approach for RW2

Step 1: The Second-Difference Matrix D

Step 2: Computing Q = D^TD

Step 3: Computing Q Step by Step (n = 5)

From u(s) to Weights: What Are We Estimating?

The Continuous Field u(s)

The SPDE defines a continuous spatial field \(u(s)\) at every location \(s\) in the domain. But we can't estimate infinitely many values! We need a finite representation.

The Basis Function Approximation

We approximate \(u(s)\) using basis functions centered at mesh vertices:

Where:

• \(\psi_k(s)\): Basis function for vertex \(k\) (piecewise linear "tent" function)
• \(w_k\): Weight at vertex \(k\) (what INLA estimates)
• \(n\): Number of mesh vertices

Connecting the Notation

The weights \(\mathbf{w}\) are the same as \(\mathbf{u}_{\text{mesh}}\):

Once we have the weights, we can compute \(u(s)\) at any location by evaluating the basis functions.

The Projector Matrix A

To get \(u(s)\) at observation locations, we evaluate the basis function expansion. The projector matrix \(\mathbf{A}\) does this efficiently:

In matrix form:

Where \(A_{ik} = \psi_k(s_i)\) is the value of basis function \(k\) at observation location \(s_i\).

Structure of A

Each row of \(\mathbf{A}\) has at most 3 non-zero values because each point lies in one triangle with 3 vertices. These values are the barycentric coordinates of the point within its triangle, and each row sums to 1.

The Complete Model

Putting it all together, a spatial model with SPDE has three levels:

General Hierarchical Model

Where \(g(\cdot)\) is the link function connecting the mean \(\mu_i\) to the linear predictor \(\eta_i\):

Where Does Q Come From?

The precision matrix \(\mathbf{Q}\) comes from the SPDE discretization. For the Matérn field with α = 2:

Where \(\mathbf{C}\) (mass matrix) and \(\mathbf{G}\) (stiffness matrix) are computed from integrating the basis functions over the mesh triangles. The matrix \(\mathbf{Q}\) encodes the spatial correlation structure and defines the prior on the weights \(\mathbf{w}\).

Notation Summary

• \(\mathbf{y}\): Observed response vector (\(m \times 1\))
• \(\eta_i\): Linear predictor at observation \(i\)
• \(g(\cdot)\): Link function (identity, log, logit, etc.)
• \(\beta_0\): Intercept (and other fixed effects)
• \(\mathbf{A}\): Projector matrix (\(m \times n\)), evaluates basis functions at observation locations
• \(\mathbf{w}\): Weights at mesh vertices (\(n \times 1\)), what INLA estimates
• \(\mathbf{A}\mathbf{w}\): The spatial field \(u(s)\) evaluated at observation locations
• \(\mathbf{Q}\): Precision matrix from SPDE (\(n \times n\)), prior precision of weights \(\mathbf{w}\)

Next Steps

Now that you understand the theory, see how to fit SPDE models in practice:

• Modeling SPDE - Step-by-step workflow for fitting SPDE models
• Creating Meshes - Detailed guide on mesh construction
• Priors - Setting PC priors for range and sigma

References

• Lindgren, F., Rue, H., & Lindström, J. (2011). An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. JRSS-B, 73(4), 423-498.
• Fuglstad, G.-A., Simpson, D., Lindgren, F., & Rue, H. (2019). Constructing priors that penalize the complexity of Gaussian random fields. JASA, 114(525), 445-452.
• Krainski, E., et al. (2018). Advanced Spatial Modeling with Stochastic Partial Differential Equations Using R and INLA. Chapman & Hall/CRC.