# gaussian_processes_on_cellular_complexes__8a92009a.pdf

Gaussian Processes on Cellular Complexes

Mathieu Alain * 1 So Takao * 1 2 Brooks Paige 1 Marc Peter Deisenroth 1

In recent years, there has been considerable interest in developing machine learning models on graphs to account for topological inductive biases. In particular, recent attention has been given to Gaussian processes on such structures since they can additionally account for uncertainty. However, graphs are limited to modelling relations between two vertices. In this paper, we go beyond this dyadic setting and consider polyadic relations that include interactions between vertices, edges and one of their generalisations, known as cells. Specifically, we propose Gaussian processes on cellular complexes, a generalisation of graphs that captures interactions between these higher-order cells. One of our key contributions is the derivation of two novel kernels, one that generalises the graph Mat ern kernel and one that additionally mixes information of different cell types.

1. Introduction

The abundance of graph-structured problems in science and engineering stimulates the development of machine learning (ML) models, such as graph neural networks (GNNs) (Scarselli et al., 2008) and graph kernel machines (Smola & Kondor, 2003). The former has achieved great success in a broad range of tasks, from molecular docking (Corso et al., 2023) to text summarisation (Fernandes et al., 2019). However, GNNs do not provide predictive uncertainty, which is an essential feature in decision-making applications. Recent work on Gaussian processes (GPs) defined on graphs (Borovitskiy et al., 2020; Nikitin et al., 2022; Opolka et al., 2022; Zhi et al., 2023) takes the latter approach, which naturally accounts for uncertainty quantification, but may lack the expressibility of GNNs.

*Equal contribution 1Centre for Artificial Intelligence, University College London, London, UK. 2Department of Computing and Mathematical Sciences, California Institute of Technology, Pasadena, CA, USA. Correspondence to: Mathieu Alain <mathieu.alain.21@ucl.ac.uk>.

Proceedings of the 41 st International Conference on Machine Learning, Vienna, Austria. PMLR 235, 2024. Copyright 2024 by the author(s).

(b) Simplicial complex

(c) Cellular complex

Figure 1. Graph, simplicial complex, and cellular complex (specifically: polyhedral complex). A simplicial complex cannot represent arbitrary polygons like the pentagon in (c).

Although graphs are an invaluable data structure, one of their main limitations is that they cannot represent interactions beyond the dyadic setting (i.e., between two vertices). However, these interactions do exist and have important applications, such as group interactions in social networks (Alvarez-Rodriguez et al., 2021), neuronal dynamics in cortexes (Yu et al., 2011), and trigenic interactions in gene networks (Kuzmin et al., 2018). Cellular complexes are a generalisation of graphs that have the ability to model such polyadic interactions (Hatcher, 2001) (see Figure 1). For this reason, they are gradually being used in ML (Hajij et al., 2020; Bodnar et al., 2021) and signal processing (Barbarossa & Sardellitti, 2020; Roddenberry et al., 2022). Although these models have helped to expand the variety of problems one can tackle with ML, as with GNNs, they typically do not quantify uncertainty.

In this paper, we fill this gap by proposing GPs defined on cellular complexes, which enables predictions on different types of cells, such as vertices and edges, while equipping them with uncertainty that is consistent with prior assumptions about the model and data. In particular, we propose the cellular Mat ern kernel, a generalisation of the graph Mat ern kernel (Borovitskiy et al., 2020), which enables modelling of signals on arbitrary cell types. By fixing an orientation on each cell (i.e., a reference direction ), this produces predictions on cells that are directed, allowing us to consider signals that have an associated direction.

Furthermore, prompted by settings with strong correlations

Gaussian Processes on Cellular Complexes

between data supported on different cell types, we also develop a new kernel, the reaction-diffusion kernel, which leverages the Dirac operator (Bianconi, 2021; Calmon et al., 2023) to mix information between different types of cells. This enables us to model various types of signal jointly so that inference on one type of cell can help the other, and vice versa. Interestingly, the Dirac operator has also been used in higher-order neural networks (Battiloro et al., 2023).

Our main contribution is to define GPs on cellular complexes, which allows us to extend the modelling capability of graph GPs in the following novel ways:

Model quantities supported on the edges, and on higherorder cells such as volumes.

Handles orientation on cells naturally. In particular, we can make directed predictions.

Allows for joint modelling of data supported on vertices, edges and higher-order cells.

2. Gaussian Processes

A Gaussian process (GP) on a set X is a random function f : X R, such that for any finite collection of points x1, . . . , x N X, the random vector (f(x1), . . . , f(x N)) RN is jointly Gaussian. GPs are defined by a mean function µ : X R and a kernel κ : X X R, which satisfy µ(x) = E[f(x)] and κ(x, x ) = Cov[f(x), f(x )], respectively, for any x, x X.

Let (x, y) := {(xi, yi)}N i=1 be training data, where yi = f(xi) + ϵi, ϵi N(0, σ2). We can make inference about f := f(x ) at arbitrary test points x by computing the posterior predictive mean and covariance

µf |y = µf + Kf f(Kff + σ2I) 1(y µf) (1)

Σf |y = Kf f Kf f(Kff + σ2I) 1Kff , (2)

respectively, where µf = µ(x), µf = µ(x ) denotes the mean and Kff = κ(x, x), Kf f = κ(x , x ), Kff = κ(x, x ) denotes the covariance and cross-covariance of f at x and x .

2.1. Gaussian processes on graphs

While most existing GPs are defined on continuous domains, in this work, we are interested in graphs G = (V, E), where we take the input set X to be the set of vertices V , and the edges E to model the proximity between two vertices. As a mathematical object, a GP in this sense is identical to a multivariate Gaussian f R|V |, whose indices are the vertices of the graph. However, the extra information provided by the edges allows one to impose a more rigid correlation structure, whereby two vertices are more strongly correlated

Figure 2. A cellular complex is constructed by attaching boundaries of k-cells ek α to the (k 1)-skeleton Xk 1 via a continuous map ϕk α.

if they are closer in the graph. Thus, this enables f to take on the characteristics of a continuous GP while being discrete.

A common way of encoding graph structures into multivariate Gaussians is by imposing specific sparsity patterns in the precision matrix, as seen in Gaussian Markov random fields (Rue & Held, 2005). For example, the Mat ern GP on graphs (Borovitskiy et al., 2021) is formally defined as f N(0, ( 2ν

ℓ2 + ) ν), where is the graph Laplacian matrix, a discrete analogue of the standard Laplacian operator on Rn. If we observe the sparsity pattern of the corresponding precision matrix ( 2ν

ℓ2 + )ν for ν N, we see that each vertex x is connected to vertices that are within a radius ν of x in the graph, controlling the smoothness of the process. More generally, we can impose a graph structure by penalising frequencies in the spectral domain of the kernel (Smola & Kondor, 2003). However, graphs are limited in the data they can support. If the data contains polyadic interactions, then a more general structure than graphs is needed.

3. Modelling with Cellular Complexes

Cellular complexes generalise graphs by incorporating higher-order interactions via cells , extending the dyadic relation modelled by edges on a graph. Concretely, a k-cell is a topological space that is homeomorphic to the unit disk Dk := {x Rk : x < 1}. We will employ the standard notation to refers to the boundary of a topological space. A finite cellular complex X of dimension n < is constructed iteratively, as follows (Hatcher, 2001):

(Step 0) Start with a collection X0 = {e0 α}N0 α=1 of 0-cells (i.e., points), called the 0-skeleton.

(Steps k = 1, . . . , n) Take a collection {ek α}Nk α=1 of k-cells and glue their boundaries to points in Xk 1 via a continuous attaching map ϕk α : ek α Xk 1. The resulting space Xk

is the k-skeleton (see Figure 2).

(Step n+1) Define the cellular complex to be the topological space X := n k=0Xk.

Gaussian Processes on Cellular Complexes

Figure 3. The path going from vertex v0 to v9 (in red) can be expressed as a 1-chain c = e1 + e2 + e5 + e7 + e9, while c represents the reverse path going from v9 to v0. Boundaries of cells can be expressed as chains, whose direction is consistent with the cell orientation. For example, h0 = e1 +e2 +e5 e4 e3 e0.

In the special case where the attaching maps ϕk α : ek α Xk 1 are embeddings, we say that the cellular complex is regular. This omits pathological cases, such as when boundaries of a k-cell collapse to a point or fold into itself. Hereafter, when we refer to cellular complexes, we will always assume that they are regular and finite. Further, when we fix an orientation on each cell (viewed as a topological manifold), we say that X is oriented, which we will also assume hereafter.

Example 1. A directed/undirected graph (see Figure 1a) is a one-dimensional oriented/non-oriented cellular complex, where the vertices are the 0-cells, the edges are the 1-cells, and the attaching maps associate the two endpoints of an edge to a pair of vertices.

Example 2. Another important class are the simplicial complexes (see Figure 1b), where the k-cells are taken to be the k-simplices and the attaching maps glue the boundaries of two k-simplices homeomorphically. Simplices are cells that contain all their sub-cells.

3.1. Chains and Cochains

Chains and cochains are key concepts on cellular complexes that formalise the notion of paths and functions over k-cells, respectively. Given an n-dimensional cellular complex X, we define a k-chain ck for 0 k n as a formal sum of k-cells

α=1 nαek α, nα Z. (3)

Intuitively, this generalises the notion of directed paths on a graph, as we show in Figure 3. We denote the set of all k-chains on X by Ck(X), which has the algebraic structure of a free Abelian group1 with basis {ek α}Nk α=1.

1One may interpret this as a vector space with coefficients restricted to the integers.

For any cellular complex X, there is a canonical operation defined on chains called the boundary operator k : Ck(X) Ck 1(X), associating the boundary of a k-chain to a (k 1)-chain. This is defined as a linear map2

α=1 na ek α, (4)

where ek α is the boundary of the k-cell ek α, expressed as a (k 1)-chain. The cycle direction of ek α must agree with the orientation of ek α, as illustrated in Figure 3 (see Appendix A for more details).

Parallel to chains, there exist dual objects known as cochains on a cellular complex X, defined as follows. Definition 3. A k-cochain on X is a linear map2 f : Ck(X) R assigning real numbers to k-chains, i.e.,

α=1 nαek α =

α=1 nαf(ek α), (5)

where f(ek α) R is the value of f at cell ek α.

The space of k-cochains, denoted Ck(X) (with a superscript), forms a real vector space with dual basis {(ek α) }Nk α=1, satisfying (ek α) ek β = δαβ. The boundary operators on chains naturally induce analogous operators on the space of cochains, referred to as the coboundary operators, defined as follows. Definition 4. For 0 k < n, the coboundary operator dk : Ck(X) Ck+1(X) is defined via the relation

dkf(c) = f( k+1c) (6)

for all f Ck(X) and c Ck+1(X). To simplify our presentation, we also define dkf 0 for k { 1, n}.

3.2. Hodge Laplacian and Dirac Operator

We introduce a generalisation of the graph Laplacian to cellular complexes, which will be used in our construction of kernels. For each k, let {wk α}Nk α=1 be a set of real, positive weights. Then, for any f, g Ck(X), we define the weighted L2-inner product as

f, g L2(wk) :=

α=1 wk α f(ek α) g(ek α). (7)

The inner product induces an adjoint of the coboundary operator, denoted d k : Ck+1(X) Ck(X), i.e.,

d kf, g L2(wk) = f, dkg L2(wk+1) (8)

for any f Ck+1(X) and g Ck(X). This leads to the following definition.

2 A group homomorphism to be more precise.

Gaussian Processes on Cellular Complexes

Object Representation

Chain c = PNk α=1 nαek α c = (n1, . . . , n Nk)

Cochain f = PNk α=1 fα(ek α) f = (f1, . . . , f Nk)

Boundary operator k Bk ZNk 1 Nk

Coboundary operator dk B k+1 ZNk+1 Nk

Inner-product f, g L2(wk) f Wkg

Hodge Laplacian k k := B k W 1 k 1Bk Wk +W 1 k Bk+1Wk+1B k+1

Table 1. Numerical representation of key objects and operations defined on cellular complexes. Here, Wk = diag(wk 1, . . . , wk Nk) is the matrix of cell-weights and Bk is the order-k incidence matrix, whose j-th column corresponds to the vector representation of the cell boundary ek j , viewed as a (k 1)-chain.

Definition 5. The Hodge Laplacian k : Ck(X) Ck(X) on the space of k-cochains is defined as

k := dk 1 d k 1 + d k dk. (9)

We defer to Section 4 to see that this generalises the graph Laplacian by considering its numerical representation.

The Dirac operator

δk := d k 1 dk : Ck(X) Ck 1 Ck+1(X) (10)

maps a k-cochain onto a direct sum of (k 1) and (k + 1)- cochains. One can verify that the relation k = δ kδk holds. Hence, the Dirac operator is understood as a formal square root of the Hodge Laplacian.

3.3. Numerical Representation

In Table 1, we display how the objects considered above can be represented as matrices and vectors in order to make explicit computations with them. We refer to Appendix B for the full derivation. In the case k = 0 and taking W0 = I (i.e., no vertex weights), we see that the expression for the Hodge Laplacian reduces to the expression for the weighted graph Laplacian 0 = B1W1B 1 . Furthermore, taking W1 = I (no edge weights), we obtain the expression for the standard graph Laplacian 0 = B1B 1 .

4. Gaussian Processes on Cellular Complexes

In this section, we establish our notion of Gaussian processes on cellular complexes, which we take to be a direct sum of Gaussian random cochains. Given a sample space Ω, an event space F and a probability measure P, we denote by (Ω, F, P) the underlying probability triple for which our random variables will be defined over.

4.1. Gaussian Random Cochains

Definition 6. A random variable f : Ω Ck(X) is called a Gaussian random cochain if for any k-chain c Ck(X), the random variable f(c) : Ω R is Gaussian.

To facilitate computations, we characterise them using a mean function and a kernel. Defining a mean is straightforward this is just a fixed cochain. For the kernel, we consider the following definition.

Definition 7. A kernel on Ck(X) is a symmetric bilinear form3 κ : Ck(X) Ck(X) R such that for any c1, . . . , cm Ck(X), we have

i,j=1 κ(ci, cj) 0. (11)

This is an appropriate notion of the kernel:

Theorem 8. A Gaussian random cochain f : Ω Ck(X) is fully characterised by a mean µ Ck(X) and a kernel κ : Ck(X) Ck(X) R.

Proof: Appendix C.1.

The vector representation f of f is simply a multivariate Gaussian f N(µ, K), whose covariance

κ(ek 1, ek 1) κ(ek 1, ek Nk) ... ... ... κ(ek Nk, ek 1) κ(ek Nk, ek Nk)

is the matrix representation of the kernel. Due to the bilinearity of κ, for any c, d Ck(X), and their vector representations c, d ZNk, we can write

κ(c, d) = c Kd. (13)

By fixing orientations on each k-cell, we also have a notion of direction for the signal f for a k-cell ek α, depending on whether the sign of f(ek α) is positive or negative, the direction of f at ek α is aligned to, or runs counter to the orientation of ek α, respectively.

4.2. GPs on Cellular Complexes

Next, we extend our notion of Gaussian random cochains to direct sums of cochains of different orders. We take this as our definition of Gaussian processes on cellular complexes.

Definition 9. Let X be an n-dimensional cellular complex. We define a Gaussian process on X as a random variable f : Ω Ln k=0 Ck(X) such that for any c = (c0, . . . , cn) Ln k=0 Ck(X), the random variable f(c) : Ω R is univariate Gaussian.

3A group bi-homomorphism to be more precise.

Gaussian Processes on Cellular Complexes

As before, we have an appropriate notion of a kernel on this space as a symmetric bilinear form3

k=0 Ck(X) R (14)

satisfying P

i,j κ(ci, cj) 0 for ci, cj Ln k=0 Ck(X). We also have the following result characterising GPs on cellular complexes via a mean and a kernel.

Theorem 10. A GP on a cellular complex X is fully characterised by a mean µ Ln k=0 Ck(X) and a kernel κ : Ln k=0 Ck(X) Ln k=0 Ck(X) R.

Proof: Appendix C.2.

Again, we can represent (14) by a matrix

K11 K1n ... ... ... Kn1 Knn

with [Knm]ij = κ(en i , em j ), which defines the covariance of f RN1+ +Nn, the vector representation of f.

Remark 11. We emphasise that our notion of a GP on a cellular complex X is not defined as a random function X R (i.e. a GP on the topological space X), but rather as a direct sum of Gaussian random cochains, i.e., a random function Ln k=0 Ck(X) R. This will allow us to define covariance structures between cells themselves, instead of between points on cells.

4.3. Kernels on Cellular Complexes

We provide some concrete examples of kernels defining GPs on cellular complexes, which encompass existing kernels in the literature. For simplicity, we take unit cell-weights wk α = 1 here and defer the treatment of the general case to Appendix D.3. For any chain c Ck(X), we also denote by c Ck(X) the cochain defined by f(c) = f, c

L2(wk) for any f Ck(X).

4.3.1. MAT ERN KERNEL

We first consider a generalisation of the Mat ern kernel on cellular complexes. Following (Borovitskiy et al., 2020; 2021), consider the stochastic system 2ν

ν/2 f = W, (16)

where f Ck(X), k is the Hodge Laplacian (Definition 5), and W : Ω Ck(X) is a Gaussian random cochain satisfying E[W(c0)] = 0 and E[W(c1)W(c2)] = c 1, c 2

L2 for any c0, c1, c2 Ck(X). The operator 2ν

ℓ2 + k ν/2 is defined rigorously in Appendix D.1.

The Mat ern kernel κ : Ck(X) Ck(X) R is then defined as a solution to the system

ν κ(c, ) = c , c Ck(X), (17)

which is a kernel in the sense that it satisfies the following property.

Proposition 12. The solution to (17) is related to the solution f of the system (16) by

κ(c, c ) = E[f(c)f(c )], c, c Ck(X). (18)

Thus, κ solving (17) is the kernel of the Gaussian random cochain f.

Proof: Appendix D.1.1.

We can express the kernel in (17) by a matrix

ℓ2 I + Λ2 ν U , (19)

where k = UΛ2U is the eigendecomposition of the Hodge Laplacian. In the case k = 0, (19) recovers exactly the graph Mat ern kernel by (Borovitskiy et al., 2021). We may also extend this construction to direct sums of cochains by replacing the Hodge Laplacian in (16) and (17) by the super Laplacian L := Ln k=0 k, resulting in kernel (19), where UΛ2U is now the eigendecomposition of the super Laplacian matrix L = blockdiag( 0, . . . , n). While this defines a GP on direct sums of cochains, due to the block-diagonal structure of the super-Laplacian, no mixing occurs between different cochains.

In the special case of GPs defined over the edges of a graph (i.e., Gaussian 1-cochains), a similar construction of the Mat ern kernel has been explored in the concurrent work Yang et al. (2024), where they also employ the Hodge decomposition to add further flexibility of the model. Another related work (Pinder et al., 2021) explores the construction of Mat ern GPs on hypergraphs to model higher-order interactions. However, the work only considers inference on the vertices, whereas our method focuses on making inferences on signals supported on the interactions themselves (i.e., the cells).

4.3.2. REACTION-DIFFUSION KERNEL

Next, we introduce a new type of kernel, which we term the reaction-diffusion kernel, that enables mixing of information between cochains of different orders. We will operate here entirely with the vector representation of cochains for ease of presentation.

Gaussian Processes on Cellular Complexes

Consider the Dirac matrix

B 1 ... ... ... ... ... ... Bn 0 . . . B n 0

whose k-th column is a numerical representation of the k-th Dirac operator (10). We can check that D2 = L holds. Thus, the Dirac matrix D and the super-Laplacian matrix L share a common eigenbasis U with eigenvalues Λ and Λ2, respectively. Now consider the stochastic system

(r I c D + d L)

ν 2 f = w, w N(0, I) (21)

for some constants r, c, d, ν R 0. Then, the corresponding kernel matrix is given by

K = U(r I cΛ + dΛ2) νU , (22)

which we term the reaction-diffusion kernel.4 This is a kernel in the sense that it satisfies the following result.

Proposition 13. The kernel defined by (22) is related to the solution f of the system (21) by

[K]ij = E[fifj], i, j. (23)

Proof: Appendix D.2.1.

Remark 14. Our naming of the reaction-diffusion kernel derives from the similarity of system (21) with the multicomponent reaction-diffusion equation

t = (r c D + d L)f, (24)

where the first and third term model the reaction and diffusion of a quantity respectively, and the second term models the cross-diffusion of multiple quantities.

To interpret this kernel, we look at the corresponding precision matrix P = K 1, which encodes the probabilistic graphical model (PGM) representation of the model f, wherein variables fi and fj are linked if and only if [P]ij = 0 (Rue & Held, 2005). For simplicity, taking ν = 1 and restricting to a 1-dimensional cellular complex, the precision takes the expression

P = r I + d 0 c B1 c B 1 r I + d 1

4In general, K is indefinite. To fix this, we set ν to be an even integer in (22), which will make it positive definite for all r, c, d R 0, excluding the set X := S

i{(r, c, d) : r cλi + dλ2 i = 0}. This set has Lebesgue measure zero. If (r, c, d) X, K becomes positive semi-definite. Therefore it defines a degenerate Gaussian measure. Inference using degenerate Gaussian measures is, however, still valid, provided σ > 0 is strictly positive in (1) (2).

(a) 1D cellular complex (b) PGM representation

Figure 4. Probabilistic graphical structure of the reaction-diffusion GP. Interactions between vertices (green) and between edges (blue) are shown as well as the mixing between cochains of different orders (red). The cellular Mat ern kernel does not have this mixing property.

For f = (f 0, f 1) in (21), the graphical structure of the k-th component f k is summarised by the matrix r I + d k for k {0, 1}, and the dependence between f 0 and f 1 is represented by the incidence matrix B1 (dotted red in 4b). Since no communication between f 0 and f 1 occurs when c = 0, the Dirac term is essential for allowing information to propagate between cochains of different orders.

Let us now consider two special cases of the kernel (22). In the first case, taking r = 2ν/ℓ2, c = 0, d = 1 and ν = ν, we see that (22) recovers the Mat ern kernel (19). Since c = 0, there is no flow of information between cochains of different orders, resulting in independence between the random cochains f 0, . . . , f n.

For r = m2, c = 1, d = 0 and ν = 2, we obtain

K = U(m I Λ) 2U = (m I D) 2. (26)

This kernel is considered by (Calmon et al., 2023) (in the form of a regulariser) for retrieving mixed topological signals supported on the k-cells for k 2.

In this section, we demonstrate the results of our GP model defined over cellular complexes (hereafter referred to as CCGP) on two examples. First, we demonstrate that CC-GPs can make directed predictions on the edges of a graph by considering the problem of ocean current interpolation. In the second example, we investigate the effect of inter-signal mixing in the reaction-diffusion kernel. We provide details of the experimental setups in Appendix E.

5.1. Directed Edge Prediction

In the numerical simulation of fluids and especially in finite element methods (FEMs), it is common to treat vector fields as signals supported on the edges of a mesh (Arnold et al., 2006) to give them the flexibility for dealing with com-

Gaussian Processes on Cellular Complexes

Figure 5. Prediction of geostrophic current around the Southern tip of Africa using the CC-Mat ern GP on edges. (Top left) Ground truth. (Top right) Predicted mean. (Bottom left) Absolute error. (Bottom right) Standard deviation. Orange dots are observed edges.

plex geometries, e.g., arising from the coastlines in ocean modelling. As there is an increasing adoption of FEM for weather and climate modelling (for example, the UK Met Office s Gung Ho model uses FEM with cubical elements (Staniforth et al., 2013)), it is of interest to consider methods that propagate information from observations directly onto the finite element vertices and edges.

Our CC-GP model can naturally be applied in this setting: Consider the geostrophic current data from the NOAA Coast Watch (2023) database. First, we project the geostrophic current around the southern tip of Africa onto the oriented edges of a two-dimensional cubical complex by averaging the flow along each edge (Desbrun et al., 2006). This yields directed signals on the edges representing the vector field. We then train our edge-based Mat ern GP (see Section 4.3.1 defined over a 1-cochain) on 30% of the data, selected randomly. The main objective of this experiment is to demonstrate that our approach can capture the directional information of the vector field, which would otherwise be difficult with existing approaches.

Figure 5 shows resulting predictions. For ease of visualisation, we display the magnitude of the predicted signals by colours on the edges; arrows indicate the predicted direction on each cell, computed by averaging the signals on its boundaries and then taking the resulting direction. We see that our CC-Mat ern GP on edges captures the general characteristics of the ground-truth vector field with similar magnitudes and directions, indicating that it can correctly diffuse information onto neighbouring edges.

We compare the results with a graph Mat ern GP baseline

MSE ( ) NLL ( )

Graph Mat ern 0.030 0.000 684.54 4.20 Edge Mat ern (ours) 0.029 0.001 703.42 5.10

Table 2. Mean square error (MSE) and negative log-likelihood (NLL) of ocean current magnitude predictions using (a) a graph Mat ern baseline, and (b) the edge Mat ern GP. Mean and standard error are shown.

MSE ( ) CC-Mat ern Reaction-diffusion

Vertices 0.165 0.005 0.076 0.004 Edges 0.335 0.014 0.200 0.010 Triangles 0.272 0.005 0.166 0.005

Vertices 28.81 1.75 -9.31 1.61 Edges 136.77 4.32 71.07 6.29 Triangles 82.84 1.56 39.78 2.47

Table 3. Mean square error (MSE) and negative log-likelihood (NLL) of predictions on the synthetic data (mean and standard error across 20 random seeds). Overall, the performance of the reaction-diffusion GP is superior to that of the Mat ern GP on the cellular complex, highlighting the benefits of mixing information across different cell types.

defined over the corresponding line graph5. Since this baseline cannot infer directions on the edge signals, we use it to only predict its magnitude. The results are shown in Table 2, where we report the mean square error (MSE) and negative log-likelihood (NLL) scores on the magnitude of the ocean current. MSE results for both models are comparable; however, our edge-based Mat ern GP performs better than the graph Mat ern GP on the NLL. This suggests that in addition to being able to infer the directions on edges, predictive uncertainties are better and the topological inductive bias contained in the edge-based Mat ern GP also helps to improve predictions for the magnitudes.

5.2. Signal Mixing

We illustrate the benefits of mixing signals using the reaction-diffusion kernel on synthetic data, where we constructed a 2D-simplicial mesh consisting of 10 10 vertices. We generated artificial signals on the edges by considering a random 1-cochain f = PK i=k ξiui, where ξi N(0, λ 1 i ), {(λi, ui)}i are the eigenpairs of the Hodge Laplacian 1 and 0 < k < K are the minimal and maximal wavenumbers controlling the smoothness of the edge field. We then generated data supported on the triangles and vertices of the mesh by applying the coboundary operator d1 and its adjoint d 1 respectively to f (numerically, this corresponds to applying

5This is the graph constructed by treating the edges as vertices and connecting them if they share a vertex.

Gaussian Processes on Cellular Complexes

(a) Vertex prediction, RD

(b) Edge prediction, RD

(c) Triangle prediction, RD

(d) Vertex prediction, CC-Mat ern

(e) Edge prediction, CC-Mat ern

(f) Triangle prediction, CC-Mat ern

Figure 6. Predictive distributions on vertices, edges, and triangles using the reaction-diffusion kernel (RD, top row) and the CC-Mat ern kernel (bottom row). Left panels: differences between mean prediction and ground truth (values close to 0 are good); right panels: corresponding predictive standard deviations. By taking correlation properties into account, the RD kernel produces better predictions on vertices, edges, and triangles.

the matrices B 2 and B1). We randomly selected a third of the data supported on each type of cell (vertices, edges and triangles) and perturbed them by noise for training. The aim is to recover the signals on the remaining two thirds of the cells.

We compare the results of the Mat ern GP on the cellular complex (CC-Mat ern) and the reaction-diffusion (RD) GP, used to make predictions across various cell-types. Figure 6 shows an example result with k = 20 and K = 100, where we display the difference (left panels) from the ground truth and the predictive standard deviation (right panels) for vertex (left column), edge (center column), and triangle predictions (right column) for the RD kernel (top row) and the CC-Mat ern kernel (bottom row). Overall, we see how the mixing of information across different cell types in the RD kernel helps to improve predictions, as evidenced by the overall smaller errors in each cell type (with small predictive uncertainty), compared to the predictions made by the CC-Mat ern kernel. In particular, we see how data supported on one type of cell can be used to enhance the predictions on another cell type as we can infer from the standard deviation plots. The standard deviation of the CC-Mat ern GP is more localised around the data points than that of the reactiondiffusion GP, suggesting that information on different cell types is being mixed in the latter.

This behaviour is also quantitatively reflected in Table 3, which reports the mean squared error (MSE) and the negative log-likelihood (NLL) scores of the predictions for both models, computed across 20 random seeds (mean and standard errors). We see how the results for the RD-GP are on average significantly better than that for the CC-Mat ern kernel on both metrics, highlighting the benefits of mixing for both prediction and uncertainty quantification.

5.3. Modelling Electromagnetic Fields

Another potential application area of GPs defined over CCs is to model electromagnetic fields, which have natural representations as cochains. In particular, the scalar potential, electric field and magnetic field generated by point charges on a 2D plane can be modelled geometrically as scalar, vector and pseudovector fields, respectively. Upon discretisation, these can then be represented by 0, 1 and 2-cochains, respectively, using the continuous-discrete correspondence between tensor fields (more precisely, differential forms) and cochains (Desbrun et al., 2006).

In this experiment, we use simulations of the scalar potential (V ), electric field (E) and magnetic field (B) generated by 10 randomly sampled point charges. These are then projected onto 0, 1 and 2-cochains, respectively, from which we extract noisy observations at a sixth of randomly selected cells in each skeleton of the generated CC. Then, similar to the experiment in Section 5.2, the objective is to recover the signals on the remaining cells. In contrast to the previous experiment, the correlation structure between signals supported on different cell types are more complex, provided indirectly through Maxwell s equations. Our goal is thus to see whether the RD kernel is still useful in this setting where the correlation between fields exist, but are not artificially imposed, as with the previous experiment. The electromagnetic simulation was performed using the Python package Py Charge (Filipovich & Hughes, 2022).

In Table 4, we compare the results of the CC-Mat ern GP and the reaction-diffusion (RD) GP on this experiment. Generally, we observe that the RD-GP yields slightly better results than the CC-Mat ern GP on both the MSE and the NLL, with the exception of the MSE on the scalar potential. For predictions of the electric and magnetic fields, the RD kernel

Gaussian Processes on Cellular Complexes

Mat ern MSE RD MSE Mat ern NLL RD NLL

V 0.102 0.008 0.11 0.011 45.75 7.20 45.41 7.35 E 0.130 0.008 0.128 0.009 104.3 14.0 100.7 15.8 B 0.190 0.026 0.178 0.023 124.7 18.8 117.3 17.2

Table 4. Comparison of independent Mat ern GPs and the reactiondiffusion (RD) GP for joint modelling of the 2D potential (V), electric (E) and magnetic (B) fields. We compare the mean square error (MSE) and negative log-likelihood (NLL). For the E and B fields, RD ourperforms independent Mat ern on 8/10 seeds.

outperformed the Mat ern kernel on eight out of ten random seeds that were used to generate the results in Table 4. Our results indicate that interestingly, mixing of information between the different cell types using the RD kernel can still be useful in this setting. However, we also note that the improvements that we see here are much less pronounced than what we observed in the previous experiment, where correlations between signals on different cell types were imposed more directly.

6. Discussion

Our experiments demonstrate the benefits of incorporating the structure of cellular complexes to model data that are naturally supported on higher-order networks, enabling us to treat (a) directed signals and (b) mixed signals with ease, which a standard graph GP cannot handle. This opens up new possibilities for modelling data on non-Euclidean domains. For instance, as a promising direction, we believe that our approach could be useful for modelling vectorial or higher-order tensorial quantities supported on arbitrary manifolds, by relying on the structural parallels between CCs and differential forms, typically employed in numerical simulations of vectorial/tensorial quantities supported on manifolds (Arnold et al., 2006). This idea is explored in our experiments in Sections 5.1 and 5.3 on simple 2D domains. We aim to extend this to more complex domains in future work, possibly incorporating boundary conditions.

Indeed, vectorial GPs on manifolds have been considered before for example in (Hutchinson et al., 2021) and more recently in (Robert-Nicoud et al., 2024). However, the former construction relies on embedding the manifold in a larger ambient space in order to make use of scalar kernels, and the latter relies on the Helmholtz decomposition to make use of scalar kernels. On the other hand, by directly encoding the vector information on the edges of a CC, one can easily model vector fields over arbitrary manifolds, possibly with boundaries (see Figure 7). Adopting a CC perspective thus makes modelling with vectorial GPs easier, and can be further generalised to the tensorial setting naturally. Further, the reaction-diffusion GP provides a topologically consistent extension to multitask GPs to the CC setting, where the individual GPs can now live on different skeletons of the

Figure 7. Left: Vector field on a manifold-with-boundaries; Right: Vector field s encoding on the edges of an oriented CC. Modelling such fields is straightforward with CC-GPs.

same CC. This will be useful in situations for modelling multiple fields (cochains) that are correlated.

We note that at present, our reaction-diffusion kernel is controlled only by three hyperparameters r, c, d, which may be too restrictive for some modelling purposes. To make our model more flexible, one possible direction that we can consider is to make use of the Hodge decomposition to split the cochains into exact, co-exact and harmonic components. This will allow us to model each of these components separately, for added flexibility. In a recent work (Yang et al., 2024), the authors consider GPs to model edge signals on a graph and observed that modelling the exact, co-exact and harmonic signals separately can lead to improved results. For signals supported on a cellular complex, a natural analogue of the Hodge decomposition can be provided through the Dirac matrix (Calmon et al., 2023). Hence, a promising direction would be to combine this decomposition with our reaction-diffusion kernel to develop a more flexible extension that could fit more complex mixing of information between different cell-types.

Another current limitation of our model is the computational cost associated with computing the eigendecomposition of the Laplacian/Dirac matrix, which can grow very quickly with the size of the network. While this limits the size of networks we can work with, this expensive computation needs to be performed only once and not during training, so we expect one to be able to work with reasonably large networks, consisting of hundreds of thousands of cells.

7. Conclusion

We introduced Gaussian processes on cellular complexes as tools for probabilistic modelling on higher-order networks. We identify these as GPs defined over random chains or direct sums thereof, enabling inference on vertices, edges, and higher-order cells. We constructed kernels appropriate for practical modelling. In particular, we generalise the Mat ern kernel to cellular complexes and propose the reaction-diffusion kernel, which allows for propagation of information between cells of different orders.

Gaussian Processes on Cellular Complexes

Acknowledgements

MA is supported by a Mathematical Sciences Doctoral Training Partnership held by Prof. Helen Wilson, funded by the Engineering and Physical Sciences Research Council (EPSRC), under Project Reference EP/W523835/1. ST is supported by a Department of Defense Vannevar Bush Faculty Fellowship held by Prof. Andrew Stuart, and by the Sci AI Center, funded by the Office of Naval Research (ONR), under Grant Number N00014-23-1-2729.

Impact Statement

This paper presents work whose goal is to advance the field of machine learning. However, due to its theoretical nature, we do not identify any potential societal consequences of our work.

Alvarez-Rodriguez, U., Battiston, F., de Arruda, G. F., Moreno, Y., Perc, M., and Latora, V. Evolutionary Dynamics of Higher-Order Interactions in Social Networks. Nature Human Behaviour, 2021. (Cited on page 1.)

Arnold, D. N., Falk, R. S., and Winther, R. Finite Element Exterior Calculus, Homological Techniques, and Applications. Acta Numerica, 2006. (Cited on pages 6 and 9.)

Barbarossa, S. and Sardellitti, S. Topological Signal Processing over Simplicial Complexes. IEEE Transactions on Signal Processing, 2020. (Cited on page 1.)

Battiloro, C., Testa, L., Giusti, L., Sardellitti, S., Lorenzo, P. D., and Barbarossa, S. Generalized Simplicial Attention Neural Networks. ar Xiv preprint ar Xiv:2309.02138, 2023. (Cited on page 2.)

Bianconi, G. The Topological Dirac Equation of Networks and Simplicial Complexes. Journal of Physics: Complexity, 2021. (Cited on page 2.)

Bodnar, C., Frasca, F., Otter, N., Wang, Y. G., Li o, P., Mont ufar, G., and Bronstein, M. Weisfeiler and Lehman Go Cellular: CW Networks. In Advances in Neural Information Processing Systems, 2021. (Cited on page 1.)

Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. Mat ern Gaussian Processes on Riemannian Manifolds. In Advances in Neural Information Processing Systems, 2020. (Cited on pages 1, 5, and 23.)

Borovitskiy, V., Azangulov, I., Terenin, A., Mostowsky, P., Deisenroth, M. P., and Durrande, N. Mat ern Gaussian Processes on Graphs. In International Conference on Artificial Intelligence and Statistics, 2021. (Cited on pages 2, 5, and 17.)

Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., Vander Plas, J., Wanderman-Milne, S., and Zhang, Q. JAX: Composable Transformations of Python+Num Py Programs, 2018. URL http://github.com/google/jax. (Cited on page 23.)

Calmon, L., Schaub, M. T., and Bianconi, G. Dirac Signal Processing of Higher-Order Topological Signals. New Journal of Physics, 2023. (Cited on pages 2, 6, and 9.)

Corso, G., St ark, H., Jing, B., Barzilay, R., and Jaakkola, T. Diff Dock: Diffusion Steps, Twists, and Turns for Molecular Docking. In International Conference on Learning Representations, 2023. (Cited on page 1.)

Desbrun, M., Kanso, E., and Tong, Y. Discrete Differential Forms for Computational Modeling. In ACM SIGGRAPH 2006 Courses. 2006. (Cited on pages 7, 8, and 23.)

Fernandes, P., Allamanis, M., and Brockschmidt, M. Structured Neural Summarization. In International Conference on Learning Representations, 2019. (Cited on page 1.)

Filipovich, M. J. and Hughes, S. Py Charge: an Opensource Python Package for Self-Consistent Electrodynamics Simulations of Lorentz Oscillators and Moving Point Charges. Computer Physics Communications, 274: 108291, 2022. (Cited on page 8.)

Hajij, M., Istvan, K., and Zamzmi, G. Cell Complex Neural Networks. In Workshop on Topological Data Analysis and Beyond at Neur IPS, 2020. (Cited on page 1.)

Hatcher, A. Algebraic Topology. Cambridge University Press, 2001. (Cited on pages 1, 2, and 13.)

Hutchinson, M., Terenin, A., Borovitskiy, V., Takao, S., Teh, Y., and Deisenroth, M. Vector-Valued Gaussian Processes on Riemannian Manifolds via Gauge Independent Projected Kernels. In Advances in Neural Information Processing Systems, 2021. (Cited on pages 9 and 15.)

Kuzmin, E., Vander Sluis, B., Wang, W., Tan, G., Deshpande, R., Chen, Y., Usaj, M., Balint, A., Usaj, M. M., van Leeuwen, J., Koch, E., Pons, C., Dagilis, A., Pryszlak, M., Wang, J. Z. Y., Hanchard, J., Riggi, M., Xu, K., Heydari, H., Luis, B.-J. S., Shuteriqi, E., Zhu, H., Dyk, N. V., Sharifpoor, S., Costanzo, M., Loewith, R., Caudy, A., Bolnick, D., Brown, G., Andrews, B., Boone, C., and Myers, C. Systematic Analysis of Complex Genetic Interactions. Science, 2018. (Cited on page 1.)

Nikitin, A., John, S., Solin, A., and Kaski, S. Non-separable Spatio-temporal Graph Kernels via SPDEs. In International Conference on Artificial Intelligence and Statistics, 2022. (Cited on page 1.)

Gaussian Processes on Cellular Complexes

NOAA Coast Watch. Sea Level Anomaly and Geostrophic Currents, Multi-mission, Global, Optimal Interpolation, Gridded. https://coastwatch.noaa.gov/cwn/products/sealevel-anomaly-and-geostrophic-currents-multi-missionglobal-optimal-interpolation.html, 2023. Altimetry data are provided by the NOAA Laboratory for Satellite Altimetry. (Cited on pages 7 and 23.)

Opolka, F., Zhi, Y.-C., Li o, P., and Dong, X. Adaptive Gaussian Processes on Graphs via Spectral Graph Wavelets. In International Conference on Artificial Intelligence and Statistics, 2022. (Cited on page 1.)

Pinder, T. and Dodd, D. GPJax: a Gaussian Process Framework in JAX. Journal of Open Source Software, 2022. (Cited on page 23.)

Pinder, T., Turnbull, K., Nemeth, C., and Leslie, D. Gaussian Processes on Hypergraphs. ar Xiv preprint ar Xiv:2106.01982, 2021. (Cited on page 5.)

Robert-Nicoud, D., Krause, A., and Borovitskiy, V. Intrinsic Gaussian Vector Fields on Manifolds. In International Conference on Artificial Intelligence and Statistics, 2024. (Cited on page 9.)

Roddenberry, T. M., Schaub, M. T., and Hajij, M. Signal Processing on Cell Complexes. In IEEE International Conference on Acoustics, Speech and Signal Processing, 2022. (Cited on page 1.)

Rue, H. and Held, L. Gaussian Markov Random Fields: Theory and Applications. CRC Press, 2005. (Cited on pages 2 and 6.)

Scarselli, F., Gori, M., Tsoi, A. C., Hagenbuchner, M., and Monfardini, G. The Graph Neural Network Model. IEEE Transactions on Neural Networks, 2008. (Cited on page 1.)

Smola, A. J. and Kondor, R. Kernels and Regularization on Graphs. In Conference on Learning Theory, 2003. (Cited on pages 1 and 2.)

Staniforth, A., Melvin, T., and Wood, N. Gungho! a New Dynamical Core for the Unified Model. In ECMWF Workshop on Recent Developments in Numerical Methods for Atmosphere and Ocean Modelling, 2013. (Cited on page 7.)

Yang, M., Borovitskiy, V., and Isufi, E. Hodge Compositional Edge Gaussian Processes. In International Conference on Artificial Intelligence and Statistics, 2024. (Cited on pages 5 and 9.)

Yu, S., Yang, H., Nakahara, H., Santos, G., Nikoli c, D., and Plenz, D. Higher-Order Interactions Characterized in Cortical Activity. Journal of Neuroscience, 2011. (Cited on page 1.)

Zhi, Y.-C., Ng, Y. C., and Dong, X. Gaussian Processes on Graphs via Spectral Kernel Learning. IEEE Transactions on Signal and Information Processing over Networks, 2023. (Cited on page 1.)

Gaussian Processes on Cellular Complexes

A. Cell Orientation and Boundaries

In this appendix, we provide more details on the computation of cell boundaries. To this end, we first require a notion of orientation on cells. For low dimensional k-cells, we have the following intuitive definitions of what we mean by an orientation:

k = 0 (i.e., a point): A choice of either + or .

k = 1 (i.e., a line segment): A choice of direction from one endpoint to the other.

k = 2 (i.e., a 2D-disk): A choice of rotation direction (clockwise or anticlockwise).

Generally, the orientation of a k-dimensional unit disk Dk is the choice of a continuously varying unit normal vector field ˆn : Dk Rk+1 on Dk, viewed as a surface embedded in Rk+1. Here, a unit normal vector ˆn(s) for s Dk is a vector in Rk+1 such that ˆn(s) = 1 and ˆn(s) v = 0 for any v Ts Dk Rk+1, where Ts Dk is the tangent space of Dk at point s (i.e., the space of all vectors in Rk+1 that are tangential to Dk). Since orientation is a topological invariant, we can define the orientation on a k-cell generally to be the orientation of the k-dimensional disk that it is homeomorphic to.

A useful way of thinking about the orientation of Dk is in terms of how it is embedded in Rk+1. In particular, we may

identify a point x Dk as a vector (x1, . . . , xk, 0) Rk+1 such that q Pk i=1 x2 i < 1. We can thus choose the unit normal field to be given by ˆn = (x1, . . . , xk, 1) Rk+1 for any x1, . . . , xk parameterising Dk. From this perspective, we can define an induced orientation ˆn on its boundary Dk = Sk 1, by taking the unit normal field pointing outwards from the disk, with respect to this embedding.

(a) Orientation and induced orientation on a 2-cell (b) Orientation and induced orientation on a 1-cell

Figure 8. Illustration of orientations (blue arrows) on (a) a 2-cell, and (b) a 1-cell, along with the induced orientation (red arrows) on the corresponding boundaries. The orientation can be understood as a choice of a continuous vector field in the ambient Euclidean space that points in the direction normal to the k-cell (in this case, pointing upwards). The induced orientation on the boundaries can be understood as the outward pointing normal field with respect to the embedding.

Example 15. Consider a 2D-disk D2 embedded in R3. The orientation on D2 is then given by the unit normal field ˆn(x, y, 0) = (x, y, 1) for p

x2 + y2 < 1 and the induced orientation on its boundary S1 = D2 is given by the outward unit normal field ˆn(x, y, 0) = (x, y, 0) for p

x2 + y2 = 1. Intuitively, the former can be thought of as anticlockwise rotation of the disk, deduced by aligning one s thumb with the unit normal direction and applying the right-hand rule. Now, aligning the thumb with ˆn and the index finger with ˆn, the middle finger, according to the right-hand rule, points in the anticlockwise direction around the boundary, giving us a more intuitive interpretation of the induced orientation (see Figure 8a).

Example 16. Consider a line-segment (i.e., a 1D disk) embedded in R2 with respect to the parameterisation {(x, 0) R2 : x ( 1, 1)}. The orientation of the line segment is determined by the unit normal field ˆn(x, 0) = (x, 1) for x ( 1, 1) and the induced orientation is given by the outward normal field ˆn(x, 0) = (x, 0) for x { 1, 1}. Again, we can use the right-hand rule to determine the intuitive interpretation of orientations here: pointing the thumb towards the page and the index finger aligned with ˆn, the middle finger points in the direction going from left endpoint x = 1 to the right endpoint x = 1. Since the middle finger as a result aligns with ˆn(x, 0) at x = 1, we think of the induced orientation at (1, 0) as having the sign + . However, it does not align with ˆn(x, 0) at x = 1, hence the induced orientation at ( 1, 0) has the sign (see Figure 8b).

Gaussian Processes on Cellular Complexes

Next, the boundary of a k-cell ek α, is defined generally as a (k 1)-chain

β=1 deg(χαβ)ek 1 β , (27)

where deg(χαβ) is the Brouwer degree of the surjection χαβ : Sk 1 = ek α ϕk α Xk 1 / Xk 1/(Xk 1 ek 1 β ) = Sk 1, mapping the boundary of ek α (homeomorphic to Sk 1) to the (k 1)-cell ek 1 β , identified with Sk 1 by collapsing its boundary to a single point (Hatcher, 2001). For a regular cellular complex, deg(χαβ) is either 0 or 1 depending on whether or not ek 1 β is a part of ek α under the attaching map. That is, if ek 1 β is not part of the boundary ek α, then we take deg(χαβ) = 0. Otherwise, if the induced orientation of ek α aligns with the orientation of ek 1 β (i.e., take an embedding Xk 1 , Rk and see how the unit normal fields of ek α and ek 1 β align), then we take deg(χαβ) = 1 and if it has opposite orientations, we take deg(χαβ) = 1. We will also adopt the notation deg(ek α ek 1 β ) to denote deg(χαβ) in order to make the cells involved more explicit.

The computation of the Brouwer degree can in general be challenging. However, in some special cases, the computation can be simplified, as we show in the following examples. Example 17. Consider a k-simplex, which can be identified as a collection of k + 1 vertices, say v0, . . . , vk without loss of generality. Its orientation is determined by the parity of the permutation of these k + 1 vertices. Hence, we represent it by an equivalence class [v0, . . . , vk], where equivalence is defined by the parity. We set ek α = [v0, . . . , vk]. Now consider an oriented face of this simplex, represented by an equivalence class of k vertices ek 1 β = [w1, . . . , wk], such that {w1, . . . , wk} = {v0, . . . , vℓ, . . . , vk} for some ℓ {0, . . . , k}. Naturally, this can be represented by a permutation σ Sk+1 on {0, . . . , k} with σ(0) = ℓand wi = vσ(i) for i = 1, . . . , k. Then, we take deg(ek α ek 1 β ) = sgn(σ), the signature of the permutation σ. Thus, in general, we can write

[v1, . . . , vk+1] =

ℓ=1 ( 1)ℓ[v1, . . . , vℓ, . . . , vk+1]. (28)

Example 18. Consider a 2D polygonal cell ek α with vertices v1, . . . , vk, and re-order them such that it revolves around the polygon in a clockwise or anticlockwise manner: i.e., vσ(1) vσ(2) vσ(k 1) vσ(k) vσ(1) for some permutation σ Sk. Whether the ordering here is clockwise or anticlockwise determines the orientation of the cell. Since the choice of the first vertex in the ordering is not important, we represent this polygon as the equivalence class ek α = [vσ(1), . . . , vσ(k)], where equivalence is defined by permutation with respect to the cyclic group Ck Sk. Now, an oriented edge of this polygon can be represented by an ordered tuple ek 1 β = (vi, vj) for i < j, such that either (i, j) = (σ(ℓ), σ(ℓ+1)) or (i, j) = (σ(ℓ+1), σ(ℓ)) for some ℓ {1, . . . , k} (we use the convention k +1 1, i.e., assume the indices are elements of Z/k Z). Then, we set deg(ek α ek 1 β ) = 1 in the former case and deg(ek α ek 1 β ) = 1 in the latter. Hence, using the convention (vi, vj) = (vj, vi), we can check that

[vσ(1), . . . , vσ(k)] =

ℓ=1 (vσ(ℓ), vσ(ℓ+1)). (29)

In the general case, we can compute the boundary of an arbitrary cell by considering its simplicial decomposition. That is, we discretise the k-cell ek α by a collection of N k-simplices. This can be expressed as a k-chain

i=1 [vi1, . . . , vik]. (30)

We can compute its boundary by first taking

i=1 [vi1, . . . , vik], (31)

then using (28) to give a simplicial decomposition of ek α, and finally collectivising with respect to the simplicial decomposition of ek 1 β .

Gaussian Processes on Cellular Complexes

B. Numerical Representation of Cellular Complexes

In this appendix, we provide a full derivation of the numerical representation of key concepts on cellular complexes, as displayed in Table 1. To make this appendix self-contained, we first re-introduce some definitions.

Definition 19 (k-chains). A k-chain on X is a free Abelian group whose generator is the set of all k-cells comprising X. The space of all k-chains on X is denoted Ck(X).

Definition 20 (Boundary operators). For k {1, . . . , n}, the boundary operator is a group homomorphism k : Ck(X) Ck 1(X) mapping k-chains to (k 1)-chains, i.e.,

α=1 na ek α, (32)

where ek α is the boundary of the cell eα k, viewed as a (k 1)-chain (see Appendix A). For convention, we also take kc 0 for k {0, n + 1}.

Definition 21 (k-cochains). A k-cochain on X is a group homomorphism f : Ck(X) R assigning real numbers to k-chains, i.e.,

α=1 nαek α =

α=1 nαf(ek α), (33)

where f(ek α) R is the value of f at cell ek α. The space of all k-cochains on X is denoted Ck(X) and forms a real vector space.

Definition 22 (Coboundary operators). For 0 k < n, the coboundary operator dk : Ck(X) Ck+1(X) is a linear map defined via the relation

dkf(c) = f( k+1c) (34)

for all f Ck(X) and c Ck+1(X). For convention, we also take dkf 0 for k { 1, n}.

Definition 23 (L2-inner product on cochains). For each k, let {wk α}Nk α=1 be a set of real, positive weights. Then, for any f, g Ck(X), we define the weighted L2-inner product as

f, g L2(wk) :=

α=1 wk α f(ek α) g(ek α). (35)

Definition 24 (L2-adjoint of the coboundary operator). For each k, let {wk α}Nk α=1 be a set of real, positive weights. The L2-adjoint of the coboundary operator, denoted d k : Ck+1(X) Ck(X) is defined by

d kf, g L2(wk) = f, dkg L2(wk+1) , (36)

for any f Ck+1(X) and g Ck(X).

Definition 25 (Hodge Laplacian). The Hodge Laplacian k : Ck(X) Ck(X) on the space of k-cochains is defined as

k := dk 1 d k 1 + d k dk. (37)

To make explicit computations with them, we wish to represent these using matrices and vectors. Fortunately, this is not difficult as the space of chains / cochains forms a free Abelian group / vector space, which is naturally isomorphic to Zn / Rn.

To this end, we fix a labelling α 7 ek α of the k-cells comprising a cellular complex X, which forms an ordered basis (ek 1, . . . , ek Nk). Then, an arbitrary k-chain c = PNk α=1 nαek α Ck(X) may be represented by a vector c = (n1, . . . , n Nk)

in ZNk. Similarly, a k-cochain f = PNk α=1 fα(ek α) Ck(X) can be represented by a vector f = (f1, . . . , f Nk) in RNk. Under this representation, cochain evaluation (33) can simply be expressed as a dot product f(c) = f c R.

Gaussian Processes on Cellular Complexes

Next, we consider the boundary and coboundary operators. The boundary operator can be expressed as a signed incidence matrix Bk ZNk 1 Nk, whose j-th column corresponds to the vector representation of the cell boundary ek j , viewed as a (k 1)-chain (see Appendix A). That is,

[Bk]ij = deg(ek j ek 1 i ) (38)

for k {1, . . . , n}, and by convention, we take Bk 0 for k {0, n + 1}. Similarly, the coboundary operator can be represented by a matrix Dk RNk+1 Nk. Using (34), we have

f D k c = f Bk+1c, (39)

Dk = B k+1. (40)

Thus, the coboundary operator is identified with the transpose of the signed incidence matrix. Finally, let Wk = diag(wk 1, . . . , wk Nk) be the weight matrix defining the L2-inner product (35). i.e.,

f, g L2(wk) = f Wkg. (41)

Then, letting D k RNk Nk+1 be the matrix representation of the adjoint of the coboundary, (36) implies

f (D k) Wkg = f Wk+1B k+1g (42)

D k = W 1 k Bk+1Wk+1. (43)

Putting this together, we find the matrix expression k RNk Nk for the Hodge Laplacian operator using (37):

k = Dk D k + D k+1Dk+1 (44)

= B k (W 1 k 1Bk Wk) + (W 1 k Bk+1Wk+1)B k+1. (45)

C. Characterisation of GPs on Cellular Complexes

Here, we prove the results in Section 4 characterising Gaussian random cochains and GPs on cellular complexes by a mean and a kernel.

C.1. Proof of Theorem 8

We restate Theorem 8 below for completeness.

Theorem 26. A Gaussian random cochain f : Ω Ck(X) is fully characterised by a mean µ Ck(X) and a kernel κ : Ck(X) Ck(X) R.

Proof. The proof is almost identical to that of Lemmas 9 10 in (Hutchinson et al., 2021).

( ) First, we show that given a Gaussian random cochain f, we can define a mean and a kernel object. For this, we simply set µ(c) := E[f(c)] and κ(c, c ) := Cov[f(c), f(c )]. The former is clearly a cochain since f is a (random) cochain. The latter can be easily checked to be a kernel in the sense of Definition 7:

(Symmetry) This follows from the symmetry of the covariance operator Cov[f(c), f(c )] = Cov[f(c ), f(c)].

(Group bi-homomorphism) This follows from the fact that by definition, f is a group homomorphism and using the bilinearity of the covariance operator.

(Positive semi-definiteness) Fixing c1, . . . , cm Ck(X), we have

α,β=1 κ(cα, cβ) =

α,β=1 Cov[f(cα), f(cβ)] = E

α=1 (f(cα) E[f(cα)]) 2 #

Gaussian Processes on Cellular Complexes

( ) Next, we show that given µ Ck(X) and a kernel κ : Ck(X) Ck(X) R, we can construct a Gaussian random cochain. Take

µ(ek 1) ... µ(ek Nk)

k(ek 1, ek 1) k(ek 1, ek Nk) ... ... ... k(ek Nk, ek 1) k(ek Nk, ek Nk)

RNk Nk. (47)

This uniquely defines a multivariate Gaussian random variable f N(µ, K). Now let φ : Ck(X) RNk be the group isomorphism identifying k-cochains by a vector in RNk via the labelling α 7 ek α. We also take the group isomorphism ψ : Ck(X) ZNk, defined by f(c) = φ(f) ψ(c) for any f Ck(X) and c Ck(X). Then, we set f := φ 1f, which defines a Gaussian random cochain, since for any c Ck(X), we have f(c) = [φ 1f](c) = f ψ(c), which is univariate Gaussian.

C.2. Proof of Theorem 10

We restate Theorem 10 below for completeness. Theorem 27. A Gaussian process on a cellular complex X (abbreviated as CC-GP) is fully characterised by a mean µ Ln k=0 Ck(X) and a kernel κ : Ln k=0 Ck(X) Ln k=0 Ck(X) R.

Proof. The proof is almost identical to the proof of Theorem 8 (see Appendix C.1), hence we will omit some details to avoid repetition.

( ) For a CC-GP f : Ω Ln k=0 Ck(X), we can define a mean and a kernel by setting µ(c) := E[f(c)] and κ(c, c ) := Cov[f(c), f(c )], for any c, c Ln k=0 Ck(X).

( ) Given µ Ln k=0 Ck(X) and a kernel κ : Ln k=0 Ck(X) Ln k=0 Ck(X) R, we take

K11 K1n ... ... ... Kn1 Knn

with [µi]j = µ(ei j) and [Knm]ij = κ(en i , em j ), which uniquely defines a multivariate Gaussian random variable f N(µ, K) in RN1+ +Nn. Now, consider the group isomorphism φ : Ln k=0 Ck(X) RN1+ +Nn by fixing a labelling α 7 ek α for each k = 0, . . . , n. Then, f := φ 1f defines a CC-GP.

D. Kernels on Cellular Complexes

In this appendix, we provide further details on the kernels that we consider in this paper, namely the Mat ern kernel on cellular complexes (CC-Mat ern) and the reaction-diffusion kernel (RD). We will need the following definition to formalise our kernels. Definition 28. Consider the weighted L2-inner product , L2(wk) : Ck(X) Ck(X) R on the space of k-cochains (Definition 23). We define a group homomorphism : Ck(X) Ck(X) by

f(c) = f, c

L2(wk), (49)

for any c Ck(X) and f Ck(X). The existence of c follows from the Riesz representation theorem.

Next, we introduce the concept of a Gaussian white-noise cochain, defined as follows. Definition 29 (Gaussian white-noise cochain). We define a zero-mean Gaussian white noise cochain as a Gaussian random cochain W : Ω Ck(X) satisfying

E[ W, f L2(wk)] = 0 for any f Ck(X).

E[ W, f L2(wk) W, g L2(wk)] = f, g L2(wk) for any f, g Ck(X)

Now, we are ready to define our notion of the Mat ern kernel on cellular complexes.

Gaussian Processes on Cellular Complexes

D.1. Mat ern Kernel

For simplicity, let us assume for now that the cell weights are wk α = 1 for all k, α. In this special case, we will use the shorthand notation , L2 , L2(wk). We will deal with the more general case in Appendix D.3.1.

Following the construction in (Borovitskiy et al., 2021), we define a Mat ern Gaussian random k-cochain as a solution to the stochastic system

ν/2 f = W, (50)

where f Ck(X) and W : Ω Ck(X) is the Gaussian white-noise cochain (Definition 29).

The operator 2ν

ℓ2 + k ν/2 is to be understood as an operation in frequency space, by the following construction. Let {(λ2 i , ui)}Nk i=1, be solutions to the eigenproblem kui = λ2 i ui such that the eigencochains {ui}Nk i=1 are orthonormal in L2. Representing f as f = P i f, ui L2(wk) ui, we define

ν/2 f, ui L2 ui, (51)

which is a linear operator on the space of k-cochains. We define the Mat ern kernel κ : Ck(X) Ck(X) R as as a solution to the linear system

ν κ(c, ) = c , c Ck(X). (52)

D.1.1. PROOF OF PROPOSITION 12

We restate Proposition 12 below for completeness.

Proposition 30. The solution to (52) is related to the solution f of the system (50) as:

κ(c, c ) = E[f(c)f(c )], c, c Ck(X) (53)

Thus, κ is the kernel corresponding to the Gaussian random cochain f.

Proof. We first claim that the unique solution f to (50) can be represented as

ν/2 W, ui L2 ui. (54)

This can be checked by simply substituting this expression inside (51) and using the L2-orthonormality of the eigencochains {ui}Nk i=1. The uniqueness can be checked by the linearity of the operator 2ν

ℓ2 + k ν/2 and the fact that the solution to the

ℓ2 + k ν/2 f = 0 is satisfied only by f 0.

Similarly, the solution to (52) is given by

L2ui( ) (55)

ν ui(c)ui( ) (56)

Gaussian Processes on Cellular Complexes

Next, we show that for arbitrary c, c Ck(X), we have

E[f(c)f(c )] = E

ν/2 W, ui L2 W, uj L2 ui(c)uj(c )

ν/2 E W, ui L2 W, uj L2 ui(c)uj(c ) (58)

ν/2 ui, uj L2 | {z } =δij

ui(c)uj(c ) (59)

ν ui(c)ui(c ) (60)

(56) = κ(c, c ), (61)

verifying property (53). Finally, we show that κ is indeed a kernel. The symmetry of κ can be easily verified from the explicit expression (56), that is, κ(c, c ) = κ(c , c) for any c, c Ck(X). Checking that κ is a group bi-homomorphism also follows easily from expression (56) using the fact that ui is a group homomorphism by the definition of cochains. Fixing c1, . . . , cm Ck(X) such that cα = 0 for some α, we also have

α,β=1 κ(cα, cβ) =

ν ui(cα)ui(cβ) (62)

β=1 ui(cβ) (63)

verifying the positive-definiteness of κ. Hence, κ is a kernel on Ck(X), as expected.

D.1.2. MATRIX REPRESENTATION

From (56), we can deduce the matrix representation of the Mat ern kernel as

ℓ2 I + Λ2 ν U . (66)

Another way to derive this representation is by directly considering the numerical representation of system (50):

Lf = w, (67)

ℓ2 I + Λ2 ν/2 U (68)

is the numerical representation of the linear operator (51), and w N(0, I). Then, we have

f = L 1w N(0, K), (69)

K = L 1L . (70)

Gaussian Processes on Cellular Complexes

We claim that

ℓ2 I + Λ2 ν/2 U . (71)

This can be verified by computing

L 1L = U 2ν

ℓ2 I + Λ2 ν/2 U U | {z } =I

ℓ2 I + Λ2 ν/2 U (72)

ℓ2 I + Λ2 ν/2 2ν

ℓ2 I + Λ2 ν/2

and similarly, LL 1 = I. Then, we can check that indeed we have

K = L 1L (76)

ℓ2 I + Λ2 ν/2 U U | {z } =I

ℓ2 I + Λ2 ν/2 U (77)

ℓ2 I + Λ2 ν U . (78)

D.2. Reaction-diffusion Kernel

Here, we provide further details on the reaction-diffusion kernel, presented in Section 4.3.2. Since many of the ideas are similar to the Mat ern kernel (Appendix D.1), we omit some details. We first lift the Dirac operator (10) to the direct sum space Ln k=1 Ck(X), where it is more natural as it defines a group homomorphism to itself (i.e., an endomorphism) D : Ln k=1 Ck(X) Ln k=1 Ck(X). This is given explicitly as

d 0f1 d0f0 + d 1f2 ... dn 2fn 2 + d n 1fn dn 1fn 1

Using the property dk+1 dk = 0 (equivalently, d k d k+1 = 0) of the coboundary operator, one can check that D2 = L (the super-Laplacian operator) holds. We also extend the L2 inner-product to the direct sum space Ln k=1 Ck(X), which we define by

f, g L2(w) :=

α=1 wk α f(ek α) g(ek α). (80)

This trivially lifts Definition 23 and Definition 29 to the direct sum setting, defining a group homomorphism : Ln k=1 Ck(X) Ln k=1 Ck(X) and a Gaussian white noise process W : Ω Ln k=1 Ck(X), respectively.

Now let W : Ω Ln k=1 Ck(X) be the white noise process on the direct sum space. We define the reaction-diffusion GP to be the solution to the stochastic system

(r + c D + d L)ν/2 f = W. (81)

As before, the operator (r + c D + d L)ν/2 is to be understood as an operation in frequency space, as follows. Let {(λi, ui)}N1+ +Nn i=1 , be solutions to the eigenproblem Dui = λiui such that {ui}N1+ +Nn i=1 are orthonormal in L2. Since

Gaussian Processes on Cellular Complexes

D2 = L, we have that D and L share the same eigenfunction ui with eigenvalues λi and λ2 i , respectively. Thus, we define

(r + c D + d L)ν/2 f :=

r + cλi + dλ2 i ν/2 f, ui L2 ui, (82)

which is a linear operator on the direct sum space Ln k=1 Ck(X). We then define the reaction-diffusion kernel to be the solution to the system

(r + c D + d L)ν k(c, ) = c , (83)

for any c Ln k=1 Ck(X). The solutions to (81) and (83) are given explicitly by

r + cλi + dλ2 i ν/2 W, ui L2 ui( ), (84)

r + cλi + dλ2 i ν ui(c)ui( ). (85)

Then, following the proof in Appendix D.1.1 line-by-line, one can verify that κ is indeed a kernel for the GP f, that is, κ(c, c ) = E[f(c)f(c )]. Below, we present a more explicit proof under the numerical representation of (84) and (85), which can be written in the form

f = U r I + cΛ + dΛ2 ν/2 U w (86)

K = U r I + cΛ + dΛ2 ν U . (87)

Here, we denoted by Λ = diag(λ1, . . . , λN1+ +Nn) the diagonal matrix of eigenvalues, w N(0, I) is the numerical representation of W, and U = (u1, . . . , u N1+ +Nn) is the matrix of eigenvectors.

D.2.1. PROOF OF PROPOSITION 13

Proposition 31. The kernel defined by (87) is related to f given by (86), as

[K]ij = E[fifj], i, j = 1, . . . , N1 + + Nn. (88)

Proof. We have

E[ff ] = U r I cΛ + dΛ2 ν/2 U E[ww ] | {z } =I

U r I cΛ + dΛ2 ν/2 U (89)

= U r I cΛ + dΛ2 ν/2 U U | {z } =I

r I cΛ + dΛ2 ν/2 U (90)

= U r I cΛ + dΛ2 ν U (91)

which proves the claim.

D.3. Generalisation to Arbitrary Cell Weights

Here, we consider the case of general cell weights, extending the results in Appendix D.1 and D.2. In particular, we demonstrate how we arrive at identical expressions for the kernels, only the eigenbasis must be orthonormal with respect to the weighted L2 inner product instead of the standard one.

Gaussian Processes on Cellular Complexes

D.3.1. MAT ERN KERNEL

Let {(λ2 i , ui)}Nk i=1 be the eigenpairs of the Hodge Laplacian operator k, defined with respect to the weighted L2-inner product , L2(wk). We claim that in this case, the eigencochains {ui}Nk i=1 can be set to be orthonormal under the weighted L2-inner product. To see this, we first show that k is self-adjoint with respect to the weighted L2-inner product:

f, kg L2(wk) = f, dk 1d k 1g L2(wk) + f, d kdkg L2(wk) (93)

= d k 1f, d k 1g L2(wk 1) + dkf, dkg L2(wk+1) (94)

= dk 1d k 1f, g L2(wk) + d kdkf, g L2(wk) (95)

= kf, g L2(wk). (96)

Now consider

ui, kuj L2(wk) = λj ui, uj L2(wk), kui, uj L2(wk) = λi ui, uj L2(wk). (97)

Since ui, kuj L2(wk) = kui, uj L2(wk) owing to the self-adjointness of k, we have λj ui, uj L2(wk) = λi ui, uj L2(wk) by (97). For λi = λj, this relation is true if and only if ui, uj L2(wk) = 0. In the case λi = λj = λ, letting E(λ) denote the corresponding eigenspace, we can simply take the orthonormal basis of E(λ) to be the elements of {ui}Nk i=1 corresponding to the eigenvalue λ. Thus, we have a choice of {ui}Nk i=1 such that ui, uj L2(wk) = δij for all i, j. Now, with this choice of the eigenbasis, we can follow the arguments in Appendix (D.1) almost identically to show that Proposition 12 still holds in the weighted setting, under the same definition for the Mat ern GP (50) and the Mat ern kernel (52). The only difference is that the operator and the Gaussian white noise W must take into consideration the cell weights, according to Definitions 23 and 29.

To illustrate this better, we consider its explicit representation in terms of a matrix-vector system. We first represent the weighted orthonormality condition of the eigenbasis by

U WU = I, (98)

where U is the matrix of eigenvectors of k and W = diag(wk 1, . . . , wk Nk). Due to the orthonormality, we also have that for any f Ck(X), we have the expression f = P

i f, ui L2(wk)ui. This has the vector expression

U(U Wf) = f, (99)

which implies

UU W = I, (100)

since f is arbitrary. We can check that the operator (51) takes the form (contrast this with (68) in the non-weighted case):

ℓ2 I + Λ2 ν/2 U W. (101)

Again, we can check that its inverse reads

L 1 := U 2ν

ℓ2 I + Λ2 ν/2 U W, (102)

by verifying

L 1L = U 2ν

ℓ2 I + Λ2 ν/2 U WU | {z } =I

ℓ2 I + Λ2 ν/2 U W (103)

ℓ2 I + Λ2 ν/2 2ν

ℓ2 I + Λ2 ν/2

= UU W (105)

Gaussian Processes on Cellular Complexes

and similarly, LL 1 = I. Now, we consider the numerical representation w of the white noise process W in the weighted setting. We claim that

w N(0, W 1). (107)

This can be checked by using its definition (Definition 29), we have

E[ W, f L2(wk) W, g L2(wk)] = f, g L2(wk) (108)

(f W)E[ww ](Wg) = f Wg (109)

E[ww ] = W 1. (110)

Hence by (67), we have

f = L 1w N(0, K), (111)

K = L 1W 1L (112)

ℓ2 I + Λ2 ν/2 U WW 1W | {z } W U 2ν

ℓ2 I + Λ2 ν/2 U (113)

ℓ2 I + Λ2 ν/2 U WU | {z } =I

ℓ2 I + Λ2 ν/2 U (114)

ℓ2 I + Λ2 ν U . (115)

The expression for the kernel is the same as in the unweighted case, except that now Λ, U are the eigenpairs of the weighted Hodge Laplacian (45), the latter being orthonormal with respect to the weighted L2-inner product instead of the standard one.

D.3.2. REACTION-DIFFUSION KERNEL

The same approach also applies to extending the reaction-diffusion kernel to the weighted setting. In this case, the matrix expression for the Dirac operator (79) reads

D0 ... ... ... ... ... ... D n 1 0 . . . Dn 1 0

0 W 1 0 B1W1 0

B 1 ... ... ... ... ... ... W 1 n 1Bn Wn 0 . . . B n 0

We can check that D is self-adjoint under the weighted L2-inner product, that is

f WDg = g WDf. (117)

Hence, we can choose an eigenbasis U of D that is orthonormal under the weighted L2-inner product, with eigenvalues given by Λ = diag(λ1, . . . , λN1+ +Nn). The operator (82) in this case can be expressed as

(r I c D + d L)ν/2 f = U r I cΛ + dΛ2 ν/2 (U 1U)U Wf (118)

= U r I cΛ + dΛ2 ν/2 U Wf. (119)

Then as before, taking

E[ww ] = W 1 (120)

and following the remaining steps in Appendix D.3.1, we arrive at the expression for the weighted reaction-diffusion kernel

K = U r I cΛ + dΛ2 ν U . (121)

This is essentially the same as in the non-weighted case, except Λ, U are now the eigenpairs of the weighted Dirac operator (116), the latter being orthonormal under the weighted L2-inner product.

Gaussian Processes on Cellular Complexes

D.3.3. AMPLITUDE OF THE PROCESS

As a special case, given a positive constant σ > 0, let us consider the weights wk α = σ 2 for all k, α. That is,

W = σ 2I. (122)

In this case, notice that conditions (98) and (100) become

U U = UU = σ2I. (123)

Thus, the normalised eigenbasis b U := σ 1U is orthonormal under the standard L2-inner product, i.e.,

b U b U = b U b U = I. (124)

Under this basis, the expressions for the Mat ern and reaction-diffusion kernels read

KMat ern = σ2 b U 2ν

ℓ2 I + Λ2 ν b U , (125)

Kr.d. = σ2 b U r I cΛ + dΛ2 ν b U . (126)

The extra parameter σ controls the amplitude of the process

Var(fi) = [K]ii = σ2ci, (127)

for ci = [ b UΦ(Λ) b U ]ii, where Φ(Λ) = 2ν

ℓ2 I + Λ2 ν in the case of the Mat ern kernel and Φ(Λ) = r I cΛ + dΛ2 ν

in the case of the reaction-diffusion kernel. This can be introduced as an extra hyperparameter in the model to fit the data more appropriately, which is recommended to obtain better results.

E. Experimental Details

In this paper, all Gaussian processes (graph Mat ern GP, CC-GP and RD-GP) are implemented using the GPJax library (Pinder & Dodd, 2022). The objective function is the conjugate marginal log-likelihood and the optimiser is an implementation of Adam from Optax (Bradbury et al., 2018) with a learning rate set at 0.1.

E.1. Directed Edge Prediction

This experiment compares our CC-GP on edges (Mat ern-CC kernel) and the graph Mat ern kernel (Borovitskiy et al., 2020). The task is to predict the edge flow constructed from the geostrophic current around the southern tip of Africa. Here, the geostrophic current refers to the dominant component of the ocean current derived by balancing the pressure gradient with the Coriolis effect. The unprocessed data is retrieved from the (NOAA Coast Watch, 2023) database, which comes in the form of two scalar fields: one representing the x-component and the other the y-component of the geostrophic current vector field. The current around the southern tip of Africa is then extracted (lat = [ 45.0, 15.0], lon = [20.0, 53.1]) and its components are rescaled to a 2D grid of dimension 20 20 (see Figure 9).

The next step in the pre-processing is to transform this data into edge signals of a cubical 1-complex. We adopt the method in (Desbrun et al., 2006) to generate these signals. To do so, a cubical mesh of the same resolution as the data (20 20) is first generated, where each edge e in the mesh is assigned an orientation. Here, the orientation is represented by a unit vector ˆte pointing from one endpoint to the other. Then for each edge e, we compute how much of the geostrophic current flows along e in the direction specified by its orientation.

More precisely, the value f0 on the edge e (illustrated in Figure 10) is computed according to

e v(s) ˆte ds v0 ˆte + v1 ˆte

where s : [0, 1] e is a parameterisation of the edge e and v is the geostrophic current. This yields directed edge signals on a cubical mesh, where we use the usual rule of setting the direction to be aligned with the orientation of e if f0 is positive

Gaussian Processes on Cellular Complexes

(a) Geostrophic current

(b) x-component (coarsened)

(c) y-component (coarsened)

Figure 9. The geostrophic current around the southern tip of Africa. (Left) Quiver plot of the geostrophic current, (Middle) The xcomponent of the geostrophic current coarsened to a 20 20 grid, (Right) The y-component of the geostrophic current coarsened to a 20 20 grid.

(a) Two vectors of a vector field. (b) Directed edge signals

Figure 10. The construction of the edge signals.

and opposite to it if f0 is negative. The training data is obtained by randomly selecting 30% of the generated edge signals and adding i.i.d. noise from a Gaussian N(0, 10 4).

For the training of the graph Mat ern GP, the smoothness hyperparameter ν is fixed at 2. The amplitude and lengthscale hyperparameters σ2, ℓare both initialised at 1.0 and optimised for 1000 iterations using Adam. The training took less than 30 seconds on a Mac Book Pro with M1 chip. In a similar way, when training CC-Mat ern GP on edges, the smoothness hyperparameter ν is set to 2, and the amplitude and lengthscale hyperparameters σ2, ℓare initialised at 1.0, before optimising them for 1000 iterations using Adam. The training for this model also takes less than 30 seconds.

E.2. Signal Mixing

This experiment compares the performance of the RD-GP and the Mat ern CC-GP in the task of predicting signals on the vertices, edges and triangles of a 2D simplicial mesh. The mesh is constructed by first defining a 10 10 grid, then subdividing this grid into triangles to transform it into a 2D simplicial mesh. The resulting complex is composed of 523 simplices: 100 vertices, 261 edges and 162 triangles.

The signals on the edges are created by taking inspiration from the Karhunen-Lo eve theorem, which states that a stochastic process can be expressed as a linear combination of L2-orthogonal basis functions with random coefficients (one may view this as a stochastic analogue of the Fourier expansion). Here, the orthogonal basis functions are the set of eigenfunctions {ui}i of the Hodge Laplacian 1. The orthogonality of the eigenfunctions is ensured by the symmetry of the operator 1. This forms a basis for edge signals (i.e. 1-cochains) that encodes the topology of the mesh through the information contained in 1. For the coefficients in the basis expansion, we use i.i.d. Gaussians ξi N(0, λ 1 i ), where λi is the eigenvalue of 1 corresponding to ui. This expansion is truncated to lie between 0 < k < K, which represent the minimal and maximal wavenumbers controlling the smoothness of the edge field. Putting this together yields the random 1-cochain

i=k ξiui, ξi N(0, λ 1 i ). (129)

Gaussian Processes on Cellular Complexes

Once the signals on the edges are obtained, the signals on the vertices and the triangles are computed by applying the coboundary operator d1 and its adjoint d 1 to f, respectively. Using the numerical representation of cochains, {ui}i becomes the set of eigenvectors of the Hodge Laplacian matrix 1, the coboundary operator becomes the matrix B1 and its adjoint becomes B 2 (see Appendix B). The vertex signals and the triangle signals are obtained by computing B1f and B 2 f, respectively. An example signal for k = 20 and K = 100 is displayed in Figure 11, which we use as the ground truth in our experiment.

(a) Vertex signals

(b) Edge signals

(c) Triangle signals

Figure 11. An example synthetic signal on the vertices, edges and triangles.

The training data is generated by randomly selecting a third of the vertices, a third of the edges and a third of the triangles from this ground truth field, and corrupting them by i.i.d. noise from a Gaussian N(0, 10 2).

For training the CC-Mat ern GP, the smoothness hyperparameter ν is fixed at 2, and the amplitude and lengthscale hyperparameters σ2, ℓare both initialised at 1.5, before optimising them for 1000 iterations using Adam. The training takes less than a minute on a Mac Book Pro equipped with a M1 Pro chip. The training of RD-GP is similar: The smoothness hyperparameter ν is fixed at 2, and the amplitude hyperparameter σ2, the reaction coefficient r, the diffusion coefficient d, and the cross-diffusion coefficient c are all initialised at 1.5. They are then optimised for 1000 iterations using Adam, again taking less than a minute to run.

E.3. Modelling Electromagnetism

In this experiment, we compare the performance of the RD-GP and the Mat ern GP on imputing signals on the vertices, edges and faces of a 20 20 square lattice. The signals come from simulations of electromagnetic fields. In particular, we used the Python package Py Charge to generate 2D electromagnetic fields on a square domain, generated by 10 oscillating point charges at randomly generated locations. The fields that were computed were the scalar potential (V ), electric field (E) and the magnetic field (B). Physically, these are a scalar field, a vector field and a two-form / pseudovector field (i.e., a field of vectors whose sign depends on the orientation of the manifold), respectively. An example of such fields is displayed in Figure 12.

(a) Scalar potential

(b) Electric field

(c) Magnetic flux into page

Figure 12. We plot an example scalar potential, electric field and magnetic field generated from ten randomly sampled oscillating point charges. For the electric field, we display the amplitudes of the vectors in colour in the background. For the magnetic field, we plot the magnetic flux going into the page, which becomes a scalar field.

The next step involves projecting these fields onto a cellular complex of dimension two, given by a 20 20 square lattice. Projecting the scalar potential on the vertices of a square lattice involves just extracting the point values of the field at

Gaussian Processes on Cellular Complexes

the vertex locations. To project the electric field onto the edges of the lattice, the procedure is similar to that described in Appendix E.1. Finally, projecting the magnetic field onto the square faces of the lattice involves averaging the magnetic flux (i.e., B ˆn) over the square cells, where the unit normal ˆn is given by the normal vector determining the orientation of the cell (see Appendix A). We also normalise the projected values, due to the large discrepancies of magnitudes between the different fields. The final projections of the fields in Figure 12 onto the cells of a square lattice are displayed in Figure 13.

(a) Scalar potential

(b) Electric field

(c) Magnetic field

Figure 13. Discrete representations of the scalar potential, electric field and magnetic field as 0, 1 and 2-cochains of a square lattice respectively. Each square cell is assigned clockwise orientation.

The training data is generated by randomly selecting a sixth of the vertices, edges and square faces from the projected fields and adding i.i.d Gaussian noise with standard deviation of 10 2.

For training the CC-Mat ern GP, the smoothness hyperparameter ν is fixed at 2, and the amplitude and lengthscale hyperparameters σ2, ℓare both initialised at 1.5, before optimising them for 1000 iterations using Adam. The training takes less than a minute on a Mac Book Pro equipped with a M1 Pro chip. The training of RD-GP is similar: The smoothness hyperparameter ν is fixed at 2, the amplitude hyperparameter σ2, the reaction coefficient r and the diffusion coefficient d are initialised at 1.5. The cross-diffusion coefficient c is initialised at 2.5. They are then optimised for 1000 iterations using Adam, again taking less than a minute to run.

The predictions made by the RD-GP is displayed in Figure 14 and those made by the CC-Mat ern GP is displayed in Figure 15. We see that both GPs recover the ground truth field (Figure 13) fairly accurately from the observations. While the metrics indicate that the RD-GP output is slightly better than those of CC-Mat ern (Table 5), perceptually, the differences are too small to see.

(a) Scalar potential

(b) Electric field

(c) Magnetic field

Figure 14. Predictions of the scalar potential, electric field and magnetic field made from the reaction-diffusion GP. The top row displays the predictive mean and the bottom row displays the standard deviations.

Gaussian Processes on Cellular Complexes

(a) Scalar potential

(b) Electric field

(c) Magnetic field

Figure 15. Predictions of the scalar potential, electric field and magnetic field made from the CC-Mat ern GP. The top row displays the predictive mean and the bottom row displays the standard deviations.

MSE ( ) CC-Mat ern Reaction-diffusion

Scalar potential 0.113 0.110 Electric field 0.125 0.108 Magnetic field 0.151 0.137

Scalar potential 70.8 68.5 Electric field 110.3 92.3 Magnetic field 118.4 110.1

Table 5. Mean square error (MSE) and negative log-likelihood (NLL) of predictions of the electromagnetic fields in Figure 14 (RD-GP) and Figure 15 (CC-Mat ern GP). The performance of the reaction-diffusion GP is slightly better than the Mat ern GP on the cellular complex, suggesting that mixing on this example has some positive impact on the predictions.