The SPDE Approach: From Gaussian Fields to GMRFs

Introduction

The SPDE approach proposed by Lindgren, Rue, and Lindström (2011) provides a powerful link between continuous Gaussian Fields (GFs) and discrete Gaussian Markov Random Fields (GMRFs). The key to this approach is the parameter mapping that connects the SPDE parameters to the interpretable parameters of the Matèrn covariance function.

The SPDE is given by:

\[ (\kappa^2 - \Delta)^{\alpha/2}(\tau \xi(s)) = W(s), \tag{1} \]

And its stationary solution has the Matèrn covariance function:

\[ \text{Cov}(\xi(s_i), \xi(s_j)) = \frac{\sigma^2}{\Gamma(\lambda) 2^{\lambda-1}} (\kappa \|s_i - s_j\|)^\lambda K_\lambda(\kappa \|s_i - s_j\|), \tag{2} \]

The link between the SPDE in Eq. (1) and the Matèrn parameters is given by the following equations:

\[ \begin{cases} \lambda = \alpha - \frac{d}{2} \\[1.2em] \sigma^2 = \frac{\Gamma(\lambda)}{\Gamma(\alpha) (4\pi)^{d/2} \kappa^{2\lambda} \tau^2} \end{cases} \]

For the case of \(s \in \mathbb{R}^2\) (\(d = 2\)), this simplifies to:

\[ \begin{cases} \lambda = \alpha - 1 \\[1.2em] \sigma^2 = \frac{\Gamma(\lambda)}{\Gamma(\alpha) (4\pi) \kappa^{2\lambda} \tau^2} \end{cases} \]

In this document, we explain the meaning of each parameter, derive the mapping, and show how to use it in practice.

Understanding the Parameters

The SPDE Parameters

The SPDE in Eq. (1) uses three parameters:

Parameter	Symbol	Role	Interpretation
Fractional power	\(\alpha\)	Controls smoothness	Higher \(\alpha\) = smoother field
Scale parameter	\(\kappa\)	Controls correlation range	Larger \(\kappa\) = faster decay
Precision scaling	\(\tau\)	Controls variance	Larger \(\tau\) = smaller variance

These parameters are convenient for mathematical operations such as:

• Fourier transforms
• Differential equations
• Finite element discretization

The Matèrn Parameters

The Matèrn covariance in Eq. (2) uses three parameters:

Parameter	Symbol	Role	Interpretation
Smoothness	\(\lambda\) (or \(\nu\))	Controls differentiability	\(\lambda = 0.5\): Exponential (rough) \(\lambda = 1\): Once differentiable \(\lambda \to \infty\): Infinitely smooth
Scale parameter	\(\kappa\)	Controls correlation range	Same as in SPDE
Marginal variance	\(\sigma^2\)	Controls amplitude	\(\text{Var}(\xi(s)) = \sigma^2\)

These parameters are convenient for statistical interpretation:

• \(\lambda\) tells us how wiggly the field is
• \(\sigma^2\) tells us how much the field fluctuates
• \(\kappa\) tells us how far the correlation extends

The General Parameter Mapping (\(d\) dimensions)

1. The Smoothness Mapping: \(\lambda = \alpha - d/2\)

Why is there a \(-d/2\) term?

This arises from the Fourier transform relationship. The Matèrn spectral density has the form:

\[ f(\omega) \propto \frac{1}{(\kappa^2 + \|\omega\|^2)^{\lambda + d/2}} \]

While the SPDE spectral density has:

\[ f(\omega) \propto \frac{1}{(\kappa^2 + \|\omega\|^2)^{\alpha}} \]

For these to match exactly, we must have:

\[ \alpha = \lambda + \frac{d}{2} \]

Rearranging gives:

\[ \lambda = \alpha - \frac{d}{2} \]

Why is this important?

• If we want a field with smoothness \(\lambda = 1\) (once differentiable) in \(d = 2\) dimensions, we must set \(\alpha = 1 + 2/2 = 2\) in the SPDE.
• If we want a field with smoothness \(\lambda = 0.5\) (Exponential covariance) in \(d = 2\), we must set \(\alpha = 0.5 + 1 = 1.5\).
• For \(d = 3\) dimensions, the same \(\lambda\) would require \(\alpha = \lambda + 1.5\).

Smoothness Interpretation

\(\lambda\)	Differentiability	Shape of Covariance	Example Process
0.5	Not differentiable	Exponential decay	Rough, fractal-like surfaces
1.0	Once differentiable	Smooth at origin	Typical environmental data
1.5	Twice differentiable	Very smooth at origin	Smooth physical fields
\(\to \infty\)	Infinitely differentiable	Gaussian covariance	Extremely smooth processes

2. The Variance Mapping: \(\sigma^2 = \frac{\Gamma(\lambda)}{\Gamma(\alpha) (4\pi)^{d/2} \kappa^{2\lambda} \tau^2}\)

Derivation

The marginal variance \(\sigma^2\) is the integral of the spectral density over all frequencies:

\[ \sigma^2 = \int_{\mathbb{R}^d} f(\omega) d\omega \]

From the SPDE solution, the spectral density is:

\[ f(\omega) = \frac{1}{\tau^2} \frac{1}{(\kappa^2 + \|\omega\|^2)^\alpha} \]

Therefore:

\[ \sigma^2 = \frac{1}{\tau^2} \int_{\mathbb{R}^d} \frac{1}{(\kappa^2 + \|\omega\|^2)^\alpha} d\omega \]

Solving the Integral

The integral \(\int_{\mathbb{R}^d} \frac{1}{(\kappa^2 + \|\omega\|^2)^\alpha} d\omega\) can be evaluated analytically. In \(d\)-dimensions, it equals:

\[ \int_{\mathbb{R}^d} \frac{1}{(\kappa^2 + \|\omega\|^2)^\alpha} d\omega = \frac{\Gamma(\alpha - d/2)}{\Gamma(\alpha)} \frac{\pi^{d/2}}{\kappa^{2\alpha - d}} \]

Using \(\lambda = \alpha - d/2\), we get:

\[ \sigma^2 = \frac{1}{\tau^2} \cdot \frac{\Gamma(\lambda)}{\Gamma(\alpha)} \frac{\pi^{d/2}}{\kappa^{2\lambda}} \]

Rearranging:

\[ \sigma^2 = \frac{\Gamma(\lambda)}{\Gamma(\alpha) \pi^{d/2} \kappa^{2\lambda} \tau^2} \]

The \((4\pi)^{d/2}\) Convention

In the Lindgren et al. (2011) paper, they use a specific Fourier transform convention that introduces a factor of \(4\):

\[ \sigma^2 = \frac{\Gamma(\lambda)}{\Gamma(\alpha) (4\pi)^{d/2} \kappa^{2\lambda} \tau^2} \]

The factor \(4\) comes from different normalization conventions for the Fourier transform. Different authors use different constants (\(2\pi\), \(4\pi\), etc.), but the key relationships remain the same:

Relationship	Interpretation
\(\sigma^2 \propto 1/\tau^2\)	Larger precision \(\tau\) = smaller variance
\(\sigma^2 \propto 1/\kappa^{2\lambda}\)	Larger scale \(\kappa\) = smaller variance (for fixed \(\tau\))
\(\Gamma\) ratios	Ensure correct normalization as \(\lambda\) changes

The 2D Case (\(d = 2\))

In most spatial statistics applications, we work in \(s \in \mathbb{R}^2\) (2D space). For this special case, the mapping simplifies:

\[ \begin{cases} \lambda = \alpha - 1 \\[1.2em] \sigma^2 = \frac{\Gamma(\lambda)}{\Gamma(\alpha) (4\pi) \kappa^{2\lambda} \tau^2} \end{cases} \]

Validity Constraint

For the field to be valid (positive variance and finite smoothness), we require:

\[ \lambda > 0 \Rightarrow \alpha > 1 \]

This means:

• \(\alpha\) must be strictly greater than 1 in 2D
• \(\alpha = 1\) would give \(\lambda = 0\), which is the limiting case of infinite roughness (not a proper field)

Examples

Example 1: Once Differentiable Field (\(\lambda = 1\))

Goal: Create a field that is once differentiable (smooth enough for most environmental applications).

Parameters: - \(d = 2\), \(\lambda = 1\) - \(\alpha = \lambda + d/2 = 1 + 1 = 2\) - \(\kappa = 0.5\) (moderate correlation range) - \(\tau = 1\) (moderate precision)

Variance Calculation:

\[ \sigma^2 = \frac{\Gamma(1)}{\Gamma(2) (4\pi) (0.5)^{2(1)} (1)^2} \]

Since \(\Gamma(1) = 1\) and \(\Gamma(2) = 1\):

\[ \sigma^2 = \frac{1}{4\pi \cdot 0.25} = \frac{1}{\pi} \approx 0.318 \]

Example 2: Rough Field (\(\lambda = 0.5\))

Goal: Create a rough field with Exponential covariance (not differentiable).

Parameters: - \(d = 2\), \(\lambda = 0.5\) - \(\alpha = \lambda + 1 = 1.5\) - \(\kappa = 0.5\) - \(\tau = 1\)

Variance Calculation:

\[ \sigma^2 = \frac{\Gamma(0.5)}{\Gamma(1.5) (4\pi) (0.5)^{2(0.5)} (1)^2} \]

Since \(\Gamma(0.5) = \sqrt{\pi}\) and \(\Gamma(1.5) = \sqrt{\pi}/2\):

\[ \sigma^2 = \frac{\sqrt{\pi}}{(\sqrt{\pi}/2) \cdot (4\pi) \cdot 0.5} = \frac{2}{4\pi \cdot 0.5} = \frac{2}{2\pi} = \frac{1}{\pi} \approx 0.318 \]

Example 3: Smooth Field (\(\lambda = 2\))

Goal: Create a very smooth field (twice differentiable).

Parameters: - \(d = 2\), \(\lambda = 2\) - \(\alpha = \lambda + 1 = 3\) - \(\kappa = 0.5\) - \(\tau = 1\)

Variance Calculation:

\[ \sigma^2 = \frac{\Gamma(2)}{\Gamma(3) (4\pi) (0.5)^{2(2)} (1)^2} \]

Since \(\Gamma(2) = 1\) and \(\Gamma(3) = 2\):

\[ \sigma^2 = \frac{1}{2 \cdot 4\pi \cdot (0.5)^4} = \frac{1}{8\pi \cdot 0.0625} = \frac{1}{0.5\pi} = \frac{2}{\pi} \approx 0.637 \]

Notice that increasing \(\lambda\) (smoothness) while keeping \(\kappa\) and \(\tau\) fixed increases the variance.

The Helmholtz Operator: A Comprehensive Explanation

What is the Helmholtz Operator?

The Helmholtz operator (also called the Helmholtz equation operator) is a fundamental differential operator that appears in the SPDE approach. It is defined as:

\[ (\kappa^2 - \Delta) \]

where: - \(\kappa^2\) is a positive constant (the square of the scale parameter) - \(\Delta\) is the Laplacian operator

Breaking Down the Components

1. The Laplacian Operator (\(\Delta\))

The Laplacian is the sum of all unmixed second partial derivatives:

In 2D space (\(s = (x, y)\)): \[ \Delta f(s) = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} \]

In 3D space (\(s = (x, y, z)\)): \[ \Delta f(s) = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} + \frac{\partial^2 f}{\partial z^2} \]

What Does the Laplacian Measure?

The Laplacian measures the local curvature or “bumpiness” of a function:

• At a local peak (maximum): The function curves downward in all directions. The second derivatives are negative, so \(\Delta f < 0\).

• At a local valley (minimum): The function curves upward in all directions. The second derivatives are positive, so \(\Delta f > 0\).

• At a flat region: The function is constant, so \(\Delta f = 0\).

• At a saddle point: The function curves up in one direction and down in another, so \(\Delta f\) can be positive, negative, or zero depending on the balance.

Visual Intuition

Imagine a rubber sheet: - High \(\Delta\): The sheet has sharp bumps or dips (high curvature) - Low \(\Delta\): The sheet is smooth and flat (low curvature) - Zero \(\Delta\): The sheet is perfectly flat or has no curvature (harmonic)

2. The Zero-Order Term (\(\kappa^2\))

The term \(\kappa^2 f(s)\) multiplies the function by a positive constant. This term:

• Ensures invertibility: Without it, the operator would have a zero eigenvalue (for constant functions), making it non-invertible.
• Controls the scale: Larger \(\kappa\) means faster decay of correlations (shorter range).
• Prevents division by zero: In the frequency domain, the operator becomes \((\kappa^2 + \|\omega\|^2)\), which is always positive.

The Helmholtz Operator in Action

When applied to a function \(f(s)\):

\[ (\kappa^2 - \Delta)f(s) = \kappa^2 f(s) - \Delta f(s) \]

Physical Interpretation

The Helmholtz operator combines two effects:

1. The \(\kappa^2\) term: A restoring force that pulls values toward zero (damping)
2. The \(-\Delta\) term: A smoothing force that penalizes curvature

Together, they create a “stiffness” operator:

• At a sharp peak (\(\Delta f < 0\)): \(-\Delta f > 0\), so the operator amplifies the peak
• At a sharp valley (\(\Delta f > 0\)): \(-\Delta f < 0\), so the operator suppresses the valley
• At a smooth region (\(\Delta f \approx 0\)): The operator simply scales by \(\kappa^2\)

The Helmholtz Operator in the Frequency Domain

The true power of the Helmholtz operator becomes clear in the frequency domain. Using Fourier transforms:

\[ \mathcal{F}\{(\kappa^2 - \Delta)f(s)\} = (\kappa^2 + \|\omega\|^2)\hat{f}(\omega) \]

where: - \(\omega\) is the frequency vector - \(\|\omega\|\) is the magnitude of the frequency - \(\hat{f}(\omega)\) is the Fourier transform of \(f\)

What Does This Mean?

The operator multiplies each frequency component by the scalar \((\kappa^2 + \|\omega\|^2)\):

• Low frequencies (small \(\|\omega\|\)): The multiplication factor is approximately \(\kappa^2\). These smooth, long-range variations are preserved.

• High frequencies (large \(\|\omega\|\)): The multiplication factor is approximately \(\|\omega\|^2\). These rapid, short-range variations are amplified.

Why Is the Helmholtz Operator Important in SPDE?

1. Local Behavior

The Helmholtz operator only involves local derivatives. This means: - It only depends on the function value at a point and its immediate neighbors - It can be discretized using sparse matrices - It enables efficient computation with the Finite Element Method (FEM)

2. Connection to the Matèrn Covariance

In the SPDE approach, the Helmholtz operator is raised to a fractional power:

\[ (\kappa^2 - \Delta)^{\alpha/2} \]

The spectral representation becomes:

\[ (\kappa^2 + \|\omega\|^2)^{\alpha/2} \]

This directly leads to the Matèrn spectral density:

\[ f(\omega) \propto \frac{1}{(\kappa^2 + \|\omega\|^2)^{\nu + d/2}} \]

where \(\alpha = \nu + d/2\).

3. Interpretation as a Filter

The Helmholtz operator acts as a spatial filter:

• Low-pass filter: When \(\kappa\) is large, the operator suppresses high frequencies less, allowing more rapid variation.
• High-pass filter: When \(\kappa\) is small, the operator suppresses high frequencies more, producing smoother fields.
• Scale control: \(\kappa\) determines the correlation length. Larger \(\kappa\) = shorter range = faster decay.

Examples

Example 1: 1D Analog

In 1D, the Helmholtz operator becomes:

\[ (\kappa^2 - \frac{d^2}{dx^2})f(x) \]

Applied to \(f(x) = \sin(kx)\):

\[ (\kappa^2 - \frac{d^2}{dx^2})\sin(kx) = (\kappa^2 + k^2)\sin(kx) \]

Example 2: 2D Gaussian Peak

Consider \(f(x,y) = e^{-(x^2 + y^2)}\) at the origin \((0,0)\):

• \(\Delta f(0,0) = -4\) (negative, since it’s a peak)
• \(\kappa^2 f(0,0) = \kappa^2\)
• \((\kappa^2 - \Delta)f(0,0) = \kappa^2 + 4\)

Example 3: 2D Valley

Consider \(f(x,y) = -e^{-(x^2 + y^2)}\) at the origin \((0,0)\):

• \(\Delta f(0,0) = 4\) (positive, since it’s a valley)
• \(\kappa^2 f(0,0) = -\kappa^2\)
• \((\kappa^2 - \Delta)f(0,0) = -\kappa^2 - 4\)

Connection to Physical Systems

The Helmholtz operator appears in many physical contexts:

1. Wave Equation

The Helmholtz equation arises from the wave equation: \[ \nabla^2 u + k^2 u = 0 \] This is the same as \((\Delta + k^2)u = 0\) or \((\kappa^2 - \Delta)u = 0\) with \(\kappa = ik\).

2. Heat Equation

The heat equation \(\frac{\partial u}{\partial t} = \alpha \Delta u\) has solutions that decay as \(e^{-\alpha \kappa^2 t}\), where \(\kappa\) determines the spatial scale.

3. Elasticity

In solid mechanics, the operator \(\kappa^2 - \Delta\) describes the stiffness of a material under deformation.

Summary Table

Aspect	Description
Definition	\((\kappa^2 - \Delta)\)
Components	Zero-order term \(\kappa^2\) + second-order term \(-\Delta\)
What it measures	Local curvature plus damping
Frequency domain	Multiplication by \((\kappa^2 + \|\omega\|^2)\)
Effect on smooth fields	Minimal amplification
Effect on rough fields	Strong amplification
Role in SPDE	Determines correlation range and smoothness
Computational advantage	Local operator → sparse matrices

Visualizing the Helmholtz Operator

# Create a simple 1D function
x <- seq(-5, 5, length.out = 100)
f <- exp(-x^2/2)  # Gaussian bump

# Compute derivatives numerically
dx <- x[2] - x[1]
df2 <- c(0, diff(diff(f)/dx)/dx, 0)  # Second derivative

# Apply Helmholtz operator for different kappa values
kappa_vals <- c(0.5, 1, 2)
colors <- c("blue", "red", "green")

# Plot
plot(x, f, type = "l", lwd = 2, col = "black", 
     ylim = c(-1, 2), 
     main = "Effect of Helmholtz Operator on a Gaussian Bump",
     xlab = "x", ylab = "Value")
legend("topright", legend = c("Original", "kappa=0.5", "kappa=1", "kappa=2"),
       col = c("black", colors), lwd = 2)

for (i in 1:length(kappa_vals)) {
  kappa <- kappa_vals[i]
  helmholtz <- kappa^2 * f - df2
  lines(x, helmholtz, col = colors[i], lwd = 2, lty = 2)
}

grid()

This plot shows how the Helmholtz operator transforms a smooth Gaussian bump for different values of \(\kappa\). Notice that: - Larger \(\kappa\) values produce a larger overall magnitude - The operator amplifies the peak and creates negative regions around it - The shape reflects the balance between the damping term (\(\kappa^2 f\)) and the curvature term (\(-\Delta f\))

Breaking Down the Operator: \((\kappa^2 - \Delta)^{\alpha/2}(\tau \xi(s))\)

Layer 1: The Input - \(\tau \xi(s)\)

• \(\xi(s)\) is the spatial Gaussian field we want to model.
• \(\tau > 0\) is a scalar that scales the field.

Purpose: Controls the overall precision (inverse variance) of the process. A larger \(\tau\) means smaller variance.

Layer 2: The Base Operator - \((\kappa^2 - \Delta)\)

As explained in detail above, this is the Helmholtz operator. Applied to a function \(f(s)\):

\[ (\kappa^2 - \Delta)f(s) = \kappa^2 f(s) - \Delta f(s) \]

• \(\kappa^2 f(s)\): The zero-order term. It multiplies the field by a positive constant, ensuring the operator is invertible.
• \(-\Delta f(s)\): The second-order term. The Laplacian \(\Delta\) measures local curvature or convexity:
  • At a local peak, \(\Delta f(s)\) is negative, so \(-\Delta f(s)\) is positive.
  • At a local valley, \(\Delta f(s)\) is positive, so \(-\Delta f(s)\) is negative.

Physical Intuition: This operator measures the “stiffness” of a rubber sheet, penalizing sharp local wiggles while rewarding smooth slopes.

Layer 3: The Fractional Power - \((\cdot)^{\alpha/2}\)

This is the most abstract part. The exponent \(\alpha/2\) applies a fractional power of the Helmholtz operator, which is defined using Fourier transforms:

In the frequency domain, the operator \((\kappa^2 - \Delta)\) becomes multiplication by the scalar \((\kappa^2 + \|\omega\|^2)\).

Raising to the power \(\alpha/2\) means:

\[ (\kappa^2 - \Delta)^{\alpha/2} f(s) = \text{InverseFT}\left\{(\kappa^2 + \|\omega\|^2)^{\alpha/2} \cdot \hat{f}(\omega)\right\} \]

What does the fractional power do?

• If \(\alpha = 2\): The operator becomes exactly \((\kappa^2 - \Delta)\). This produces a field with standard smoothness.
• If \(\alpha = 4\): The operator becomes \((\kappa^2 - \Delta)^2\), involving fourth-order derivatives, producing a much smoother field.
• If \(\alpha = 1\) (or any non-integer): This involves fractional derivatives, which are non-local operators that look at the entire spatial domain, not just immediate neighbors.

The Crucial Trick: By choosing a fractional \(\alpha\), we can achieve any desired level of smoothness for the field, not just integer derivatives.

Layer 4: The Full Expression

When the compound operator acts on \(\tau \xi(s)\), it outputs a new spatial function representing the residual forcing or local innovation:

• If \(\xi(s)\) is smooth (low curvature), the output is small.
• If \(\xi(s)\) is rough (high curvature), the output is large.

The SPDE states that this output must equal spatial white noise \(W(s)\):

\[ (\kappa^2 - \Delta)^{\alpha/2}(\tau \xi(s)) = W(s) \]

Since white noise has no spatial correlation, the equation means: “Apply this differential filter to the field to strip away all spatial structure; the residuals must be completely independent in space.”

Time-Series Analogy

If you are familiar with time-series, this equation is the spatial equivalent of an Autoregressive (AR) model:

• Time-series: \((1 - \phi B)^p y_t = \epsilon_t\) (where \(B\) is the backshift operator, and \(\epsilon_t\) is white noise).
• Spatial: \((\kappa^2 - \Delta)^{\alpha/2}(\tau \xi(s)) = W(s)\).

The operator \((\kappa^2 - \Delta)^{\alpha/2}\) acts like the polynomial \((1 - \phi B)^p\). It is a differential filter that transforms a correlated spatial field \(\xi(s)\) into uncorrelated white noise \(W(s)\). The parameters \(\kappa\), \(\tau\), and \(\alpha\) control the shape and strength of this filter, ultimately determining the correlation range, variance, and smoothness of the resulting field.

Understanding Equation (2): The Matèrn Covariance

What Does the Equation Say?

Equation (2) describes the covariance between two points in the spatial field:

\[ \text{Cov}(\xi(s_i), \xi(s_j)) = \frac{\sigma^2}{\Gamma(\lambda) 2^{\lambda-1}} (\kappa \|s_i - s_j\|)^\lambda K_\lambda(\kappa \|s_i - s_j\|) \]

• \(\xi(s_i)\) and \(\xi(s_j)\): Values of the spatial process at two distinct locations.
• \(\sigma^2\): The marginal variance (variance at a single point).
• \(\Gamma(\lambda)\): The Gamma function, a normalizing constant.
• \(\|s_i - s_j\|\): The Euclidean distance between the two points.
• \(K_\lambda(\cdot)\): The modified Bessel function of the second kind, which gives the Matèrn its characteristic shape:
  • For large distances, it decays exponentially.
  • For small distances, it behaves like a power-law.

Key Property: The covariance depends only on the distance between the points (stationarity and isotropy). The correlation starts at 1 at distance 0, decays to 0 as distance goes to infinity, and the rate of decay is governed by \(\kappa\) and \(\lambda\).

The “Famous” Spectral Density of the Matèrn

The Matèrn covariance has a famous spectral density (its Fourier transform):

\[ f(\omega) \propto \frac{1}{(\kappa^2 + \|\omega\|^2)^{\nu + d/2}} \]

This is famous because it has three extraordinary properties:

1. Power-Law Tails: For high frequencies, it decays as \(1/\|\omega\|^{2\nu+d}\), which matches the fractal behavior of many real-world spatial phenomena.

2. Rational Form: It is a rational function of \(\|\omega\|^2\), making it the spatial equivalent of Autoregressive (AR) processes in time-series analysis. This is the key that allows connection to the SPDE.

3. Continuity in Smoothness: It spans the gap between the Exponential covariance (\(\nu = 0.5\)) and the Gaussian covariance (\(\nu \to \infty\)), covering all levels of smoothness.

Derivation: From the SPDE to the Matèrn Covariance

Step 1: Fourier Transform the SPDE

Apply the Fourier transform to both sides of Equation (1):

\[ \mathcal{F}\{(\kappa^2 - \Delta)^{\alpha/2}(\tau \xi(s))\} = \mathcal{F}\{W(s)\} \]

The Fourier transform of the Laplacian is: \(\mathcal{F}\{\Delta \xi(s)\} = -\|\omega\|^2 \hat{\xi}(\omega)\).

Therefore, the Fourier transform of the operator is: \((\kappa^2 + \|\omega\|^2)^{\alpha/2} \hat{\xi}(\omega)\).

The Fourier transform of white noise \(W(s)\) is a constant (usually 1).

Thus:

\[ (\kappa^2 + \|\omega\|^2)^{\alpha/2} \tau \hat{\xi}(\omega) = 1 \]

Solving for \(\hat{\xi}(\omega)\):

\[ \hat{\xi}(\omega) = \frac{1}{\tau} \frac{1}{(\kappa^2 + \|\omega\|^2)^{\alpha/2}} \]

The spectral density \(f(\omega)\) (the squared magnitude of the Fourier coefficients) is:

\[ f(\omega) = \frac{1}{\tau^2} \frac{1}{(\kappa^2 + \|\omega\|^2)^\alpha} \]

Step 2: Parameter Substitution

The classic Matèrn spectral density is:

\[ f(\omega) = \frac{\sigma^2 2^d \pi^{d/2} \Gamma(\nu + d/2) \kappa^{2\nu}}{\Gamma(\nu)} \cdot \frac{1}{(\kappa^2 + \|\omega\|^2)^{\nu + d/2}} \]

For the SPDE solution to match this, we must set:

\[ \alpha = \nu + \frac{d}{2} \]

Substituting this into the SPDE spectrum:

\[ \hat{\xi}(\omega) = \frac{1}{(\kappa^2 + \|\omega\|^2)^{\nu/2 + d/4}} \]

And the spectral density becomes:

\[ f(\omega) = \frac{1}{(\kappa^2 + \|\omega\|^2)^{\nu + d/2}} \]

This perfectly matches the denominator of the classic Matèrn spectrum.

Step 3: Inverse Fourier Transform

By the Wiener-Khinchin theorem, the covariance is the inverse Fourier transform of the spectral density:

\[ C(h) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} \frac{1}{(\kappa^2 + \|\omega\|^2)^{\nu + d/2}} e^{i\omega \cdot h} d\omega \]

where \(h = \|s_i - s_j\|\) and \(\omega\) is the frequency vector.

Reducing to a 1D Integral

Because the spectrum is radially symmetric, the \(d\)-dimensional Fourier transform reduces to a 1D Hankel transform. In hyperspherical coordinates:

\[ C(h) = \frac{1}{(2\pi)^{d/2}} \int_0^\infty \frac{r^{d-1}}{(\kappa^2 + r^2)^{\nu + d/2}} (rh)^{1-d/2} J_{d/2-1}(rh) dr \]

where \(r = \|\omega\|\) and \(J\) is the Bessel function of the first kind.

The Standard Integral Solution

This integral is a known form of the Hankel transform of a Bessel potential. The solution is:

\[ \int_0^\infty \frac{r^{d-1}}{(\kappa^2 + r^2)^{\nu + d/2}} (rh)^{1-d/2} J_{d/2-1}(rh) dr = \frac{\pi^{d/2} 2^{1-\nu}}{\Gamma(\nu)} h^\nu \kappa^{d/2-\nu} K_\nu(\kappa h) \]

where \(K_\nu\) is the modified Bessel function of the second kind.

Substituting Back and Simplifying

Plugging this back into the covariance integral:

\[ C(h) = \frac{1}{(2\pi)^{d/2}} \cdot \frac{\pi^{d/2} 2^{1-\nu}}{\Gamma(\nu)} h^\nu \kappa^{d/2-\nu} K_\nu(\kappa h) \]

Simplify the constants:

\[ C(h) = \frac{2^{1-\nu-d/2}}{\Gamma(\nu)} \kappa^{d/2-\nu} h^\nu K_\nu(\kappa h) \]

To get the form with \((\kappa h)^\nu\), multiply and divide by \(\kappa^\nu\):

\[ C(h) = \frac{2^{1-\nu-d/2}}{\Gamma(\nu)} \kappa^{d/2-2\nu} (\kappa h)^\nu K_\nu(\kappa h) \]

Normalizing to Get the Marginal Variance

Define the marginal variance \(\sigma^2\) such that:

\[ \frac{2^{1-\nu-d/2}}{\Gamma(\nu)} \kappa^{d/2-2\nu} = \frac{\sigma^2}{\Gamma(\nu) 2^{\nu-1}} \]

Then:

\[ C(h) = \frac{\sigma^2}{\Gamma(\nu) 2^{\nu-1}} (\kappa h)^\nu K_\nu(\kappa h) \]

Replacing \(\nu\) with \(\lambda\) and writing \(h = \|s_i - s_j\|\), we obtain exactly Equation (2):

\[ \text{Cov}(\xi(s_i), \xi(s_j)) = \frac{\sigma^2}{\Gamma(\lambda) 2^{\lambda-1}} (\kappa \|s_i - s_j\|)^\lambda K_\lambda(\kappa \|s_i - s_j\|) \]

Visualizing the Effects of Parameters

Effect of \(\alpha\) on Smoothness

# Function to compute Matèrn covariance
matern_cov <- function(h, kappa, nu, sigma2 = 1) {
  if (h == 0) {
    return(sigma2)
  } else {
    Kv <- besselK(kappa * h, nu)
    return(sigma2 / (gamma(nu) * 2^(nu - 1)) * (kappa * h)^nu * Kv)
  }
}

# Create distance grid
h_vals <- seq(0, 10, length.out = 200)

# Define different alpha values for d=2
# alpha = nu + 1
alpha_vals <- c(1.5, 2, 2.5, 3)
nu_vals <- alpha_vals - 1  # lambda

# Parameters
kappa <- 1
tau <- 1

# Compute sigma^2 for each alpha
sigma2_vals <- sapply(nu_vals, function(nu) {
  alpha <- nu + 1
  gamma(nu) / (gamma(alpha) * (4*pi) * kappa^(2*nu) * tau^2)
})

# Compute covariance for each combination
results <- data.frame()
for (i in 1:length(alpha_vals)) {
  temp <- data.frame(
    h = h_vals,
    covariance = sapply(h_vals, matern_cov, 
                        kappa = kappa, 
                        nu = nu_vals[i],
                        sigma2 = sigma2_vals[i]),
    alpha = paste("alpha =", alpha_vals[i]),
    nu = paste("nu =", round(nu_vals[i], 2))
  )
  results <- rbind(results, temp)
}

# Plot
ggplot(results, aes(x = h, y = covariance, color = alpha, group = alpha)) +
  geom_line(size = 1.2) +
  labs(
    title = "Effect of Alpha (Smoothness) on the Covariance Function",
    subtitle = "d = 2, kappa = 1, tau = 1",
    x = "Distance (h)",
    y = "Covariance C(h)",
    color = "Alpha"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 12),
    legend.position = "bottom"
  ) +
  scale_color_viridis(discrete = TRUE)