---
title: "FEM Matérn methods: SPDE, B-splines, and covariance comparison"
author: "Patrick Brown"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{FEM Matérn methods: SPDE, B-splines, and covariance comparison}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup}
library(adlaplaceGrf)
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:

```r
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.

# 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:

```{r params}
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)
```

Evaluation sites (same for both grids):

```{r sites}
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.

```{r compare}
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)))
)
```

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

```{r plot-corr, fig.width=7, fig.height=4}
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).