Discrete Exterior Calculus-based Whittle-Matérn Fields on Meshes

Introduction

Gaussian Processes (GPs) are a commonly-used tool for spatial modelling under uncertainty. A GP is a multivariate distribution where the joint distribution of every finite set of variables is Gaussian [16]. A GP is wholly defined by its mean and covariance functions, and the mean and variance of an unobserved spatial index is defined by the GP and any observed data. An isotropic GP is one where the covariance function is a function of the Euclidean distance between two spatial locations, and common GPs, such as the Matérn or squared exponential GPs, are isotropic.

However, on curved manifolds like spheres, Euclidean distance functions are not sufficient for defining spatial covariance, and using geodesic distance instead generally does not produce positive definite covariance. Borovitskiy et al. [6] introduce an alternative formulation for scalar covariance kernels defined on an arbitrary Riemannian manifold $(M, g)$ through the spectral decomposition of the Laplace-Beltrami operator $-\Delta_g$. Both Hutchinson et al. [12] and Robert-Nicoud et al. [17] extend this idea to modelling vector fields, the former through projected kernels, and the latter through applying derivative operators to the eigenfunctions of the Laplace-Beltrami to get vector fields. In particular, Robert-Nicoud et al. [17] show that Hodge-compositional processes can be constructed in this way, allowing separate parametrisations of the curling and diverging velocities in a vector field. Borovitskiy et al. [5] show that the spectral formulation of the Matérn GP can be applied to graphs and triangular meshes through the eigenpairs of the graph Laplacian, and show empirically that it converges to its smooth counterpart in the Riemannian setting. Yang et al. [21] propose a methodology for modelling functions on the edges of meshes, extending the formulation of intrinsic vector GPs with the graph Laplacian. Gillan et al. [9] show that the Hodge-compositional vector field model of Robert-Nicoud et al. [17] has a discrete analogue by substituting the discrete exterior calculus (DEC) cotangent Laplacian for $-\Delta_g$ and showing that it can be used to model vector fields defined on the vertices of a manifold mesh. However, this precludes a sparse formulation. Alain et al. [2] develop a model for Matérn GPs over the components of a cellular complex, analogous to the simplicial complex of DEC, and propose a kernel for propagating information between different orders of cells. Alain and Takao [1] develop a similar model to Gillan et al. [9] using finite element exterior calculus, eliciting a similar covariance-based Hodge-compositional model without the restriction to mesh vertices.

Standard dense Gaussian-process calculations require $O(N^2)$ factorisation and $O(N^2)$ storage. The SPDE approach Lindgren et al. [14] explicitly constructs Matérn precision matrices for Gaussian Markov Random Fields (GMRFs) with respect to a mesh by discretising the SPDE using piece-wise linear basis functions, enabling sparse linear-algebra methods. As highlighted in a later review [13], this result can be generalised to smooth manifolds by replacing the Euclidean $\Delta$ with the equivalent manifold operator $\Delta_g$. Nambiarajan (2026) show that a Gaussian Markov Random Field (GMRF) approximation for Whittle-Matérn fields can be built from DEC operators for fields on vertices and faces.

Background

Discrete Exterior Calculus

Exterior calculus is the calculus defined by generalising the derivative to operate on higher-degree differential forms [8]. A differential $k$-form is a quantity that can be integrated over a $k$-dimensional manifold; for example, the expression $h(x)dx$ is a 1-form and can be integrated over a 1-dimensional line segment by

The discretisation of exterior calculus is called discrete exterior calculus (DEC) [11]. An implementation of most of the core theory is provided in the PyDEC library [3].

In the following, $\omega$, $\eta$ are smooth $k$-forms, $\alpha$, $\beta$ are discretised $k$-forms, cochains, $c$ is a $k$-chain, $M$ is the Riemannian manifold that is triangulated by the simplicial complex $K$, and $\sigma^k$ is a simplex of dimension $k$. For an in-depth treatment see, e.g. Desbrun et al. [7].

Simplicial Complexes, Chains and Cochains

A smooth Riemannian manifold $(M,g)$ of dimension $n$ can be discretised as a finite collection of 0-, 1-, $\ldots$, $(n-1)$-, and $n$-dimensional piecewise linear components. For a 2-dimensional manifold, these would be vertices (0-dimensional), edges (1-dimensional), and faces (2-dimensional). This collection of components is called a simplicial complex, $K$, and the individual $k$-dimensional components, $\sigma^k$, are called $k$-simplices. We assume $K$ triangulates $M$ and that $K$ is embedded in Euclidean space with vertex positions chosen so that its lengths, areas and angles approximate those of $(M,g)$. We write $K^{(k)}$ for the set of $k$-simplices (e.g. $K^{(1)}$ is the collection of edges on a triangulated surface) and $N_k = |K^{(k)}|$.

The orientation of a $k$-simplex is the ordering of its vertices $[v_0,\ldots,v_k]$, up to even permutations. The chain group $C_k(K)$ is the abelian group generated by the oriented $k$-simplices, where opposite orientations are negative. Fixing one orientation per simplex, $C_k(K)$ is the free abelian group on $K^{(k)}$,

and we set $C_k(K) = 0$ for $k < 0$ or $k > n$. The orientations of individual simplices are an arbitrary choice.

A discrete $k$-form, or $k$-cochain, is a homomorphism $\alpha: C_k(K) \to\mathbb{R}$, and these form the cochain group

A cochain is determined by its values $\alpha(\sigma^k)$ on the simplices, and $C^k(K) \cong\mathbb{R}^{N_k}$; a 0-cochain is a function on vertices, a 1-cochain a value on each oriented edge, and so on. We write

for the evaluation of a cochain on a chain.

The circumcentric dual of $K$ associates each $k$-simplex $\sigma^k$ with an $(n-k)$-dimensional dual cell $\star\sigma^k$ built from the circumcentres of the simplices containing $\sigma^k$. On a triangulated surface, the dual of a face is its circumcentre, the dual of an edge is the segment joining the circumcentres of its adjacent faces, and the dual of a vertex is the polygon joining the circumcentres of its incident faces. We write $|\sigma^k|$ and $|\star\sigma^k|$ for the $k$- and $(n-k)$-dimensional volumes of a simplex and its dual cell, taking the volume of a point to be 1. On well-centred meshes (every simplex contains its circumcentre) all dual volumes are positive.

de Rham and Whitney Maps

Smooth and discrete forms are related through the de Rham and Whitney maps [19]. The de Rham map $\psi^k$ discretises a smooth $k$-form $\omega$ by integration over $k$-simplices,

$\psi^k\omega$ is a well-defined cochain in $C^k(K)$, and by linearity $\langle\psi^k\omega,c\rangle=\int_c\omega$ for every chain $c\in C_k(K)$.

In the opposite direction, the Whitney map $\phi^k$ interpolates a cochain $\alpha\in C^k(K)$ as a piecewise-smooth $k$-form,

Writing $\lambda_i$ for the barycentric coordinate of vertex $v_i$, equal to one at $v_i$ and zero at all other vertices, the Whitney forms of degree zero, one and two are

where $f=[v_i,v_j,v_l]$ is a face on a triangulated surface and $\mathrm{d}A_f$ is its area form, oriented by the vertex ordering of $f$. These satisfy $\int_{\sigma^k}W_{\sigma^{\prime k}}=\delta_{\sigma^k\sigma^{\prime k}}$ for all $\sigma,\sigma'\in K^{(k)}$, so that

i.e. interpolating a cochain and then discretising recovers it exactly.

Exterior Derivative

The discrete boundary operator $\partial_k:C_k(K)\to C_{k-1}(K)$ is defined by

where $[v_0,\ldots,\hat{v}_i,\ldots,v_k]$ is the $(k-1)$-dimensional simplex obtained by omitting the vertex $v_i$ [7]. $\partial_k$ is represented by a sparse signed incidence matrix with entries in $\{0,\pm1\}$, and $\partial_k\partial_{k+1}=0$.

The discrete exterior derivative, $\mathrm{d}_k:C^k(K)\to C^{k+1}(K)$, is the discrete difference operator acting on $k$-cochains. For any $(k+1)$-chain $c$, Stokes’ theorem can be written as an adjoint relationship

We define $\mathrm{d}_k$ as the adjoint of the discrete boundary operator,

It follows from $\partial_k\partial_{k+1}=0$ that $\mathrm{d}_{k+1}\mathrm{d}_k=0$. Additionally, we can see that $\mathrm{d}_k$ commutes with the de Rham map $\psi^k$,

for any smooth $k$-form $\omega$ and $c\in C_{k+1}(K)$. Since a cochain is determined by its values on chains, $\mathrm{d}_k\psi^k=\psi^{k+1}\mathrm{d}$, and $\mathrm{d}_k$ introduces no discretisation error. In matrix form, $\mathrm{d}_k$ is the transpose of the signed incidence matrix of $\partial_{k+1}$, and is simply a sparse difference matrix as a result of the sparsity of $\partial_{k+1}$.

Hodge Star

In the DEC framing, all metric information enters through an inner product on cochains. For smooth $k$-forms on $(M,g)$, the $L^2$ inner product is

where $\{e_i\}$ is any local orthonormal frame. We seek a symmetric positive definite matrix $\star_k$ such that

and define the discrete inner product

The natural choice is the Gram matrix of the Whitney forms, $(\star_k^W)_{\sigma\sigma'} = (W_\sigma,W_{\sigma'})_{L^2}$, for which $(\alpha,\beta)_k = (\phi^k\alpha,\phi^k\beta)_{L^2}$ exactly; $\star_k^W$ is then the finite-element mass matrix of $k$-cochains. In practice, DEC uses a diagonal (lumped) approximation of this mass matrix. For $\alpha= \psi^k\omega$, $\alpha(\sigma^k) = \int_{\sigma^k}\omega\approx\bar{\omega}_{\sigma^k}|\sigma^k|$, and the ratio $\alpha(\sigma^k)/|\sigma^k|$ approximates the average density of $\omega$ along $\sigma^k$. Giving $\sigma^k$ a region of measure proportional to $|\sigma^k||\star\sigma^k|$ and summing density squared times volume gives

where $\star_k$ is the diagonal Hodge star,

$\star_k$ is symmetric positive definite on well-centred meshes, where all $|\star\sigma^k|$ are positive.

Co-derivative

The discrete co-derivative $\delta_{k+1}: C^{k+1}(K) \to C^k(K)$ is defined as the adjoint of $d_k$ with respect to the inner product in Eq. (14),

Since $\star_k$ is symmetric positive definite, and hence invertible,

so

Just as $d_k$ is the adjoint of $\partial_{k+1}$ under the pairing $\langle\cdot,\cdot\rangle$, $\delta_{k+1}$ is the adjoint of $d_k$ under the inner product $(\cdot,\cdot)_k$.

Hodge-Laplacian

The discrete Hodge-Laplacian² on $k$-cochains is

with $\delta_0 = 0$ and $d_n = 0$, since $C^{-1}(K) = C^{n+1}(K) = 0$. $\Delta_k$ is self-adjoint and positive semidefinite with respect to $(\cdot,\cdot)_k$, since $(\Delta_k\alpha,\alpha)_k = (d_k\alpha,d_k\alpha)_{k+1} + (\delta_k\alpha,\delta_k\alpha)_{k-1} \ge0$. $\Delta_k$ as a matrix is not symmetric, but $\star_k\Delta_k$ is.

Whittle-Matérn Fields

The Matérn covariance function for locations $\mathbf{u},\mathbf{v}\in\mathbb{R}^n$ and Euclidean distance $\lVert\cdot\rVert$ is defined as

where $K_\nu$ is the modified Bessel function of the second kind with order $\nu> 0$, $\kappa> 0$ is a scaling parameter and $\sigma^2$ is the marginal variance. The Matérn random field has mean-square derivatives of integer orders $m < \nu$. Lindgren et al. [14] suggests a more natural interpretation of $\kappa$ as a range parameter $\rho$, with an empirical reparametrisation $\rho= \frac{\sqrt{8\nu}}{\kappa}$ giving correlations near 0.1 when $\lVert v-u\rVert= \rho$. Whittle [20] formalises the Matérn GP as the solution to the SPDE

with the same $\kappa> 0$ and $\alpha= \nu+ \frac{n}{2}$, and a white noise process $\mathcal{W}$. $\Delta$ is the Euclidean Laplacian and the marginal variance of the field is given by

Methodology

Whittle-Matérn Fields on $k$-Cochains

When a Gaussian vector is defined implicitly by a linear system driven by Gaussian noise, its precision is defined by the system matrix and the noise covariance. Suppose

with $S$ invertible. Then $x=S^{-1}b$ is Gaussian with

We can translate the Whittle-Matérn SPDE formulation in Eq. (22) to DEC by replacing $-\Delta$ with the $k$-form Hodge-Laplacian $\Delta_k$ as defined in Eq. (20).

Note the change of sign between Eqs. (22) and (26) since $\Delta$ is negative semi-definite and $\Delta_k$ is positive semi-definite.

To discretise Eq. (26) for $\alpha=2$, we use its weak form: we require $((\kappa^2+\Delta_k)x,v)_k=(\mathcal{W},\phi_v^k)_{L^2}$ for all $v\in C^k(K)$. For each $\sigma\in K^{(k)}$, the left-hand side is the $\sigma$-th entry of $\star_k(\kappa^2+\Delta_k)x=Sx$, and the right-hand side is $b_\sigma=(\mathcal{W},W_\sigma)_{L^2}$, giving $Sx=b$. Since $\kappa^2+\Delta_k$ is self-adjoint with respect to $(\cdot,\cdot)_k$, $S$ is symmetric, and its explicit form follows from the quadratic form $u^{\mathsf{T}}Su=((\kappa^2+\Delta_k)u,u)_k$,

where

$S\in\mathbb{R}^{N_k\times N_k}$ is sparse since $\star_k^{-1}$ is diagonal. Using $\star_k^W$ would make the term $\star_kd_{k-1}(\star_{k-1}^W)^{-1}d_{k-1}^{\mathsf{T}}\star_k$ dense, which justifies our use of the lumped mass matrix. The white noise process $\mathcal{W}$ is a random linear functional such that

Using Whitney forms $\{W_\sigma\}$ as a basis, $b_\sigma:=(\mathcal{W},W_\sigma)_{L^2}$, and

Approximating $\star_k^W$ by the lumped mass matrix $\star_k$ gives $b \sim\mathcal{N}(0,\star_k)$. Thus, when $\alpha=2$, the weak form of Eq. (26) is

and the precision matrix of $x$ is given by

For general $\alpha\ne2$, the fractional power in Eq. (26) is most naturally defined after whitening. The discrete inner product satisfies

so the change of variables $\tilde{x}=\star_k^{1/2}x$ transforms the $\star_k$ inner product to the Euclidean one. The self-adjoint operator $\kappa^2+\Delta_k$ is then represented by the symmetric positive definite matrix

using $S=\star_k(\kappa^2+\Delta_k)$. The noise process becomes

For $\alpha=2$, $Sx=b$ is equivalent to $\tilde{S}\tilde{x}=\tilde{b}$. For general $\alpha$ the whitened form of Eq. (26) is

Hence $\operatorname{Cov}(\tilde{x})=\tilde{S}^{-\alpha}$. Transforming back via $x=\star_k^{-1/2}\tilde{x}$ gives $\operatorname{Cov}(x)=\star_k^{-1/2}\tilde{S}^{-\alpha}\star_k^{-1/2}$, and so

For integer $\alpha$, we can use telescoping to define $Q_\alpha$,

For $\alpha=1$ this is $Q_1=S$, and for $\alpha=2$ it recovers $Q_2=S\star_k^{-1}S$. For non-integer $\alpha$, $\tilde{S}^{\alpha}$ remains well defined through spectral powers but is dense, so $Q_\alpha$ loses both its sparsity and the Markov property. We therefore restrict attention to integer $\alpha$; fractional values can be handled by rational approximation [4].

The smoothness parameter is related to $\alpha$ by $\nu=\alpha-n/2$, and $\nu>0$ requires $\alpha>n/2$. With integer $\alpha$, $\nu$ takes integer values when $n$ is even and half-integer when $n$ is odd. On a surface ($n=2$), $\alpha=1$ gives $\nu=0$, and so does not define a Matérn covariance function, but is still a well-defined random measure [14].

Implementation

Prior Sampling

To sample from $\mathcal{N}(0,Q_\alpha^{-1})$ we require $L$ such that $LL^{\mathsf{T}}=Q_\alpha$. We consider the cases where $\alpha=2m$ and $\alpha=2m+1$ separately. When $\alpha= 2m$,

When $\alpha= 2m+1$,

where

A sample $x$ is then given by $x = L^{-\mathsf{T}}z$, $z \sim\mathcal{N}(0,I)$, i.e.

Figure 1 shows samples of $k$-cochains on a unit sphere.

Posterior Conditioning

The posterior covariance $\operatorname{Cov}(x \mid y) = Q_{\mathrm{post}}^{-1}$ where

$A$ is an observation operator that maps the latent cochain to the observations, $y = Ax+\epsilon$, and $R$ is the covariance of the observation noise, $\epsilon\sim\mathcal{N}(0,R)$ e.g. for i.i.d noise with variance $\tau^2$, $R = \tau^2I$. For $R = \tau^2I$, the posterior mean of $x$ is given by

$Q_{\mathrm{post}}^{-1}$ is a dense matrix, but we can find the mean and take posterior samples without ever forming it explicitly. For the mean we employ two triangular solves

For posterior samples, we employ one additional triangular solve,

We can find the posterior marginal variance approximately through the Monte Carlo sample variance, or exactly by solving Takahashi equations to find only $\operatorname{diag}(Q_{\mathrm{post}}^{-1})$ [18]. Figure 2 shows posterior inference on a 2-cochain on a sphere.

Observation Operators

The latent field is a cochain $x \in C^k(K)$ with prior $\mathcal{N}(0,Q_\alpha^{-1})$. Any observation that is a linear functional of the underlying form (e.g. an integral, a point value, or a vector) enters through the observation operator $A$ in $y = Ax+\epsilon$, and the posterior of § 4.2 applies unchanged.

Mesh-aligned observations. If $y_i$ is an integral over a chain $c \in C_k(K)$, then by the de Rham map $\int_c \omega= \langle\psi^k\omega,c\rangle$, so the corresponding row of $A$ is the coefficient vector of $c$. For a single simplex this is a signed indicator, with sign given by the relative orientation. Note that this includes scalar values observed at vertices of the mesh.

Point observations. For non-mesh aligned observations, we evaluate the Whitney map $\phi_x^k$ at the observation location $p$ as in Eq. (5), using the barycentric coordinates $\lambda_v(p)$ of $p$. Only the Whitney forms of the face $f$ containing $p$ contribute, so each row of $A$ is sparse. For $k = 1$, when observations are vectors rather than 1-forms, we use a vector proxy for $\phi_x^1$,

where $\nabla\lambda_i$ is constant on $f$. An observed vector is compared with $E_p^{\mathsf{T}}\mathbf{v}(p)$, where $E_p \in\mathbb{R}^{3 \times2}$ is the observation’s local frame.

Vector-field priors. For a collection of points $p_1,\ldots,p_m$, we can define the operator $B : C^1(K) \to\mathbb{R}^{2m}$ by evaluating $E_{p_i}^{\mathsf{T}}\mathbf{v}(p_i)$ at each $p_i$. The prior on vector fields is the push-forward of the cochain prior through $B$, $\mathbf{v} = Bx \sim\mathcal{N}(0, BQ_\alpha^{-1}B^{\mathsf{T}})$. We can avoid constructing this dense covariance explicitly; samples are obtained by sampling $x$ as in Section 4.1 and applying $B$, and vector observations enter inference through the corresponding rows of $B$ in $A$.

Experiments

As a demonstration, we apply this model to interpolation of a vector field on a sphere. The data are monthly averaged $u$ and $v$ components of wind at the 500hPa pressure level for the month of July 2020, taken from the ERA5 atmospheric reanalysis dataset [10]³, scaled to have an average norm of 1. For observations we subsample the data to a $10^\circ\times10^\circ$ grid, and infer back to the $0.25^\circ\times0.25^\circ$ grid, giving 1,035,360 output locations. The ground truth data is shown in Figure 3 and the observations are shown in Figure 4.

We begin by generating an icosphere mesh with 5 subdivisions, producing a well-centred mesh with 10242 vertices, 30720 edges, and 20480 faces, and using PyDEC we assemble a simplicial complex. For hyperparameter values, we fix $\alpha= 2$ ($\nu= 1$) and estimate $\rho$ by empirically measuring the spatial correlation of the observations. For our data this gives $\rho= 0.615$ rad and so $\kappa= 4.597$.

Using the DEC operators defined by the simplicial complex, we construct $Q_2$ for $k = 1$ as in Eq. (31), giving a prior precision over the edges of the mesh. The sparsity pattern of this matrix is shown in Figure 5. We then construct our observation operator $A$ as the Whitney map for vector field priors described in § 4.3. Note that since we use an icosphere, neither the observations nor the ground truth are (necessarily) aligned with the vertices of the mesh. We set $\tau= 0.001$ and use $R = \tau^2 I$ for the observation noise.

Using the procedure detailed in § 4.2, we obtain the posterior mean and standard deviation, shown in Figure 6 and Figure 7, respectively. The posterior estimate has a mean square error of 0.046 with respect to the ground truth, and the spatial distribution of square error is shown in Figure 8.

Conclusion

This work presents a sparse precision formulation for Whittle–Matérn fields for $k$-cochains using DEC operators. By applying Whitney maps and introducing the Whitney observation operators, this model can be applied to any point, edge, areal, or vector field data, with observation and prediction locations anywhere on the underlying manifold. Future development will focus on boundary conditions and non-stationarity, as well as sparse approximations of Hodge-compositional covariance kernels.