# Implicit Gaussian process representation of vector fields over arbitrary
latent manifolds

Robert L. Peach*
University Hospital Würzburg
<EMAIL_ADDRESS>&Matteo Vinao-Carl*, Nir Grossman, Michael David
Imperial College London
(* Indicates equal contribution) Emma Mallas, David Sharp, Paresh A. Malhotra
Imperial College London
&Pierre Vandergheynst, Adam Gosztolai
EPFL
<EMAIL_ADDRESS>
###### Abstract

Gaussian processes (GPs) are popular nonparametric statistical models for
learning unknown functions and quantifying the spatiotemporal uncertainty in
data. Recent works have extended GPs to model scalar and vector quantities
distributed over non-Euclidean domains, including smooth manifolds appearing
in numerous fields such as computer vision, dynamical systems, and
neuroscience. However, these approaches assume that the manifold underlying
the data is known, limiting their practical utility. We introduce RVGP, a
generalisation of GPs for learning vector signals over latent Riemannian
manifolds. Our method uses positional encoding with eigenfunctions of the
connection Laplacian, associated with the tangent bundle, readily derived from
common graph-based approximation of data. We demonstrate that RVGP possesses
global regularity over the manifold, which allows it to super-resolve and
inpaint vector fields while preserving singularities. Furthermore, we use RVGP
to reconstruct high-density neural dynamics derived from low-density EEG
recordings in healthy individuals and Alzheimer’s patients. We show that
vector field singularities are important disease markers and that their
reconstruction leads to a comparable classification accuracy of disease states
to high-density recordings. Thus, our method overcomes a significant practical
limitation in experimental and clinical applications.

##  1 Introduction

A cornerstone of statistical learning theory is the manifold assumption, which
posits that high-dimensional datasets are often distributed over low-
dimensional smooth manifolds – topological spaces characterised by locally
Euclidean structure. For instance, images of an object from varying camera
angles or diverse renditions of a written letter can all be viewed as samples
from a smooth manifold (Tenenbaum, 2000). Further, the common approximation of
data by a proximity graph, based on a notion of affinity or similarity between
data points, induces a Riemannian structure that is instrumental in geometric
learning theories. For example, the analogy between the graph Laplacian matrix
and the Laplace-Beltrami operator associated with a Riemannian manifold
(Chung, 1997) has been widely exploited in manifold learning (Belkin & Niyogi,
2003; Coifman et al., 2005), shape analysis (Taubin, 1995), graph signal
processing (Ortega et al., 2018), discrete geometry (Gosztolai & Arnaudon,
2021), graph neural networks (Defferrard et al., 2016; Kipf & Welling, 2017;
Peach et al., 2020) and Gaussian processes (Borovitskiy et al., 2020; 2021).

However, many datasets contain a richer structure comprising a smoothly
varying vector field over the manifold. Prime examples are dissipative
dynamical systems where post an initial transient phase, trajectories converge
to a manifold in state space (Fefferman et al., 2016). Likewise, in
neuroscience, smooth vector fields arise from the firing rate trajectories of
neural populations evolving over neural manifolds, which is instrumental in
neural information coding (Sussillo & Barak, 2013; Khona & Fiete, 2022;
Gardner et al., 2022). Smooth vector fields are also pertinent in areas like
gene expression profiling during development (La Manno et al., 2018) and
multireference rotational alignment in cryoelectron microscopy (Singer & Wu,
2012). The breadth of these applications emphasises the need to generalise
current learning paradigms to capture both the manifold structure and its
associated vector field.

To address this need, a promising avenue is to consider the Laplace-Beltrami
operator as a hierarchy of Laplacians that act on tensor bundles of a manifold
with increasing order. The first member of this hierarchy is the Laplace-
Beltrami operator, which acts on rank-$0$ tensors, i.e., scalar signals.
Similarly, higher-order signals, including vector fields, have associated
Laplacian operators, which can encode their spatial regularity. Among these,
the connection Laplacian (Barbero et al., 2022), defined on vector bundles,
and the related sheaf Laplacian (Knöppel et al., 2013; Bodnar et al., 2022),
which allows the vector spaces on nodes to have different dimensions, are
emerging as leading operators in machine learning (Bronstein et al., 2017;
Gosztolai et al., 2023). Both of these operators are related to heat diffusion
of higher-order signals over manifolds (Singer & Wu, 2012; Sharp et al., 2019)
and thus intrinsically encode the signals’ smoothness. The connection
Laplacian is particularly appealing because it can be constructed, even when
the manifold is unknown, from graph-based data descriptions (Singer & Wu,
2012; Budninskiy et al., 2019). We, therefore, asked how one could use this
discrete approximation to derive continuous functions that implicitly
represent the vector field over the manifold. Such representation could take
account of the global regularity of the vector field to reconstruct intricate
vector field structures lost in data sampling.

Gaussian processes – a renowned family of nonparametric stochastic processes –
offer an excellent framework for learning implicit functional descriptions of
data. While GPs are traditionally defined on Euclidean spaces (Rasmussen &
Williams, 2006), several studies have extended them to Riemannian manifolds.
However, these studies have either considered scalar signals (Wilson et al.,
2021; Mallasto & Feragen, 2018; Mallasto et al., 2020; Borovitskiy et al.,
2020; 2021; Jensen et al., 2020) or vector signals Hutchinson et al. (2021)
but only in cases where the underlying manifold is known and is analytically
tractable, such as spheres and tori. In this work, we generalise GPs to vector
fields on arbitrary latent manifolds, which can only be approximated based on
local similarities between data points, making them applicable in practice to
real-world datasets.

Our contributions are as follows. (i) We generalise GPs to vector-valued data
using the connection Laplacian operator, assuming that the data originates
from a stationary stochastic process. (ii) We show that the resulting
Riemannian manifold vector field GP (RVGP) method encodes the manifold and
vector field’s smoothness as inductive biases, enabling out-of-sample
predictions from sparse or obscured data. (iii) To underscore the practical
implications of RVGP, we apply it to electroencephalography (EEG) recordings
from both healthy individuals and Alzheimer’s disease patients. The global
spatial regularity learnt by our method significantly outperforms the state-
of-the-art approaches for reconstructing high-density electrical fields from
low-density EEG arrays. This enables our method to better resolve vector field
singularities and dramatically increase the classification power of disease
states.

Taken together, our technique opens the door to a differential geometric
formulation of kernel-based operators and demonstrates a direct relevance for
fundamental and clinical neuroscience.

##  2 Comparison with related works

Let us begin by carefully comparing our method to related works in the
literature.

##### Implicit neural representations (INRs)

There has been an increasing interest in defining signals implicitly as
parametrised functions from an input domain to the space of the signal. In
Euclidean spaces, INRs have been a breakthrough in replacing pixel-wise
description images or voxel-wise descriptions of 3D shapes by neural networks
(Sitzmann et al., 2020; Lipman, 2021; Gosztolai et al., 2021; Koestler et al.,
2022; Mildenhall et al., 2020). INRs have also been extended to signals over
manifolds by graph Laplacian positional encoding (Grattarola & Vandergheynst,
2022). However, INRs are data-intensive due to the lack of explicit spatial
regularisation. Further, they have not been extended to handle vector-valued
data.

##### Gaussian processes over specific Riemannian manifolds

Several closely related works in the GP literature have provided various
definitions of GPs on Riemannian manifolds. One line of works defined GPs as
manifold-valued processes $f:{\mathbb{X}}\to\mathcal{M}$ (Mallasto & Feragen,
2018; Mallasto et al., 2020) by using the exponential map of the manifold to
perform regression in the tangent space. However, these works require that the
manifold $\mathcal{M}$ be known in order to define the exponential map. More
notable are studies which define GPs as scalar-valued functions
$f:{\mathbb{X}}\to\mathbb{R}$. For example, considering data points as samples
from a stationary stochastic process, the domain ${\mathbb{X}}$ can be defined
based on positional encoding using eigenfunctions of the Laplace-Beltrami
operator (Solin & Särkkä, 2020; Borovitskiy et al., 2020) or the graph
Laplacian (Borovitskiy et al., 2021). However, these works cannot be directly
applied to vector-valued signals by treating vector entries as scalar
channels. This is because these channels are generally not independent but
related through the curvature of the manifold. To address this gap, Hutchinson
et al. (2021) defined GPs as functions $f:{\mathbb{X}}\to\mathcal{TM}$ over
the tangent bundle $\mathcal{TM}$ by first mapping the manifold isometrically
to Euclidean space and then using a multi-input, multi-output GP to learn the
projected signal. However, Hutchinson et al. (2021) has focused on cases when
the manifold is explicitly known, and its mapping to Euclidean space can be
explicitly defined. Here, we are specifically interested in the case where the
manifold is unknown – a common scenario in many scientific domains. To achieve
this, we generalise the works of Borovitskiy et al. (2020; 2021) to vector
fields on Riemannian manifolds using the connection Laplacian operator and its
eigenvectors as positional encoding.

##### Gaussian processes in neuroscience

GPs have also been widely used in the neuroscience literature, particularly
combined with linear dimensionality reduction to discover latent factors
underlying neural dynamics. One popular method is Gaussian Process Factor
Analysis (GPFA) (Yu et al., 2009), which defines GPs in the temporal domain
and does not encode spatial regularity over ensembles of trajectories as
inductive bias. GPFA has been used to define time-warping functions to align
neural responses across trials, i.e., individual presentation of a stimulus or
task (Duncker & Sahani, 2018). Likewise, GPFA has been extended to non-
Euclidean spaces by simultaneously identifying the latent dynamics and the
manifold over which it evolves (Jensen et al., 2020). However, this model is
limited to manifolds built from $SO(3)$ symmetry groups and requires them to
be enumerated to perform Bayesian model selection. We instead seek a
constructive framework that requires no assumption on the manifold topology.

##  3 Background

Here, we introduce the function-space view of GPs in Euclidean domains and
their stochastic partial differential equation formulation developed for
scalar-valued signals, which will then allow us to extend them to vector-
valued signals.

###  3.1 Gaussian processes in Euclidean spaces

A GP is a stochastic process $f:{\mathbb{X}}\to\mathbb{R}^{d}$ defined over a
set ${\mathbb{X}}$ such that for any finite set of samples
$X=({\bm{x}}_{1},\dots,{\bm{x}}_{n})\in{\mathbb{X}}^{n}$, the random variables
$(f({\bm{x}}_{1}),\dots,f({\bm{x}}_{n}))\in\mathbb{R}^{n\times d}$ are jointly
multivariate Gaussian, $\mathcal{N}({\bm{x}};{\bm{\mu}},{\bm{K}})$, with mean
vector ${\bm{\mu}}=m(X)$ and covariance matrix ${\bm{K}}=k(X,X)$.
Consequentially, a GP is fully characterised by its mean function
$m({\bm{x}}):=\mathbb{E}(f({\bm{x}}))$ and covariance function
$k({\bm{x}},{\bm{x}}^{\prime}):=\mathrm{Cov}(f({\bm{x}}),f({\bm{x}}^{\prime}))$,
also known as kernel, and denoted as
$f\sim\mathcal{GP}(m({\bm{x}}),k({\bm{x}},{\bm{x}}^{\prime}))$. It is typical
to assume that ${\bm{m}}({\bm{x}})=0$, which does not reduce the expressive
power of GPs (Rasmussen & Williams, 2006).

One may obtain the best-fit GP to a set of training data points
$(X,{\bm{y}})=\\{({\bm{x}}_{i},{\bm{y}}_{i})|i=1,...,n\\}$ by Bayesian linear
regression, assuming that the observations ${\bm{y}}_{i}$ differ from the
predictions of $f$ by some Gaussian measurement noise, i.e.,
${\bm{y}}_{i}=f({\bm{x}}_{i})+{\mathbf{\epsilon}}_{i}$, where
$f\sim\mathcal{GP}(0,k)$ and
${\mathbf{\epsilon}}_{i}\sim\mathcal{N}(0,\sigma_{n}^{2})$ for some standard
deviation $\sigma_{n}$. Then, the distribution of training outputs ${\bm{y}}$
and model outputs ${\bm{f}}_{*}:=f({\bm{x}}_{*})$ at test points
${\bm{x}}_{*}$, is jointly Gaussian, namely

$\begin{bmatrix}{\bm{y}}\\\
{\bm{f}}_{*}\end{bmatrix}\sim\mathcal{N}\left(0,\begin{bmatrix}k(X,X)+\sigma_{n}^{2}\mathbf{I}&k(X,X_{*})\\\
k(X_{*},X)&k(X_{*},X_{*})\end{bmatrix}\right).$ (1)

To generate predictions for a test set
$X_{*}=({\bm{x}}_{1},\dots,{\bm{x}}_{n^{*}})$, one can derive the posterior
predictive distribution conditioned on the training set (Rasmussen & Williams,
2006), namely
$f_{*}|X_{*},X,\mathbf{y}\sim\mathcal{N}({\bm{\mu}}_{|\mathbf{y}},{\bm{K}}_{|\mathbf{y}})$
whose mean vector and covariance matrix are given by the expressions

$\displaystyle{\bm{\mu}}_{|\mathbf{y}}({\bm{x}}_{*})$
$\displaystyle=k(X_{*},X)(k(X,X)+\sigma_{n}^{2}\mathbf{I})^{-1}{\bm{y}},$ (2)
$\displaystyle{\bm{K}}_{|\mathbf{y}}({\bm{x}}_{*},{\bm{x}}_{*})$
$\displaystyle=k({\bm{x}}_{*},{\bm{x}}_{*})-k(X_{*},X)(K+\sigma^{2}\mathbf{I})^{-1}k(X,X_{*}).$
(3)

The advantage of GPs is that the smoothness of the training set regularises
their behaviour, which is controlled by the kernel function. We focus on
kernels from the Matérn family, stationary kernels of the form:

$k_{\nu}({\bm{x}},{\bm{x}}^{\prime})\equiv
k_{\nu}({\bm{x}}-{\bm{x}}^{\prime})=\sigma^{2}\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\frac{||{\bm{x}}-{\bm{x}}^{\prime}||}{\kappa}\right)K_{\nu}\left(\sqrt{2\nu}\frac{||{\bm{x}}-{\bm{x}}^{\prime}||}{\kappa}\right)$
(4)

for $\nu<\infty$, where $\Gamma(\nu)$ is the Gamma function and $K_{\nu}$ is
the modified Bessel function of the second kind. Matérn family kernels are
favoured due to their interpretable behaviour with respect to their
hyperparameters. Specifically, $\sigma,\kappa,\nu$ control the GP’s
variability, smoothness and mean-squared differentiability. Moreover, the
well-known squared exponential kernel, also known as radial basis function
$k_{\infty}({\bm{x}}-{\bm{x}}^{\prime})=\sigma^{2}\exp\left(-||{\bm{x}}-{\bm{x}}^{\prime}||^{2}/2\kappa^{2}\right)$
is obtained in the limit as $\nu\to\infty$.

###  3.2 Scalar-valued GPs on Riemannian manifolds

In addition to their interpretable hyperparameters, Matérn GPs lend themselves
to generalisation over non-Euclidean domains. A formulation that will allow
extension to the vector-valued case is the one by Whittle (1963) who has shown
that in Euclidean domains ${\mathbb{X}}=\mathbb{R}^{d}$ Matérn GPs can be
viewed as a stationary stochastic process satisfying the stochastic partial
differential equation

$\left(\frac{2\nu}{\kappa^{2}}-\Delta\right)^{\frac{\nu}{2}+\frac{d}{4}}f=\mathcal{W},$
(5)

where $\Delta$ is the Laplacian and $\mathcal{W}$ is the Gaussian white noise.
Likewise, for $\nu\to\infty$, the limiting GP satisfies
$\exp(-\kappa^{2}\Delta/4)f=\mathcal{W}$, where the left-hand side has the
form of the heat kernel.

As shown by Borovitskiy et al. (2020), Eq. 5 readily allows generalising
scalar-valued ($d=1$) Matérn GPs to compact Riemannian manifolds
${\mathbb{X}}=\mathcal{M}$ by replacing $\Delta$ by the Laplace-Beltrami
operator $\Delta_{\mathcal{M}}$. The corresponding GPs are defined by
$f_{\mathcal{M}}\sim\mathcal{GP}(0,k_{\mathcal{M}})$, with kernel

$k_{\mathcal{M}}({\bm{x}},{\bm{x}}^{\prime})=\frac{\sigma^{2}}{C_{\nu}}\sum_{i=0}^{\infty}\left(\frac{2\nu}{\kappa^{2}}+\lambda_{i}\right)^{-\nu-\frac{d}{2}}f_{i}({\bm{x}})f_{i}({\bm{x}}^{\prime}),$
(6)

where $(\lambda_{i},f_{i})$ are the eigenvalue-eigenfunction pairs of
$\Delta_{\mathcal{M}}$ and $C_{\nu}$ is a normalisation factor.

###  3.3 Scalar-valued GPs on graphs

Analogously, scalar-valued Matérn GPs can be defined on graph domains,
${\mathbb{X}}=G$ (Borovitskiy et al., 2021) by using the graph Laplacian in
place of $\Delta_{\mathcal{M}}$ as its discrete approximation. Specifically,
the graph Laplacian is ${\bm{L}}:={\bm{D}}-{\bm{W}}$ with weighted adjacency
matrix ${\bm{W}}\in\mathbb{R}^{n\times n}$ of and diagonal node degree matrix
${\bm{D}}=\text{diag}({\bm{W}}\bm{1}^{T})$. The graph Laplacian admits a
spectral decomposition,

${\bm{L}}={\bm{U}}{\bm{\Lambda}}{\bm{U}}^{T},$ (7)

where ${\bm{U}}$ is the matrix of eigenvectors and ${\bm{\Lambda}}$ is the
diagonal matrix of eigenvalues and, by the spectral theorem,
$\Phi({\bm{L}})={\bm{U}}\Phi({\bm{\Lambda}}){\bm{U}}^{T}$ for some function
$\Phi:\mathbb{R}\to\mathbb{R}$. Therefore, choosing
$\Phi(\lambda)=\left(2\nu/\kappa^{2}+\lambda\right)^{\nu/2}$ obtains the
operator on the left-hand side of Eq. 5 up to a scaling factor111Note the
different signs due to the opposite sign convention of the Laplacians.. Using
this, one may analogously write

$\left(\frac{2\nu}{\kappa^{2}}-{\bm{L}}\right)^{\frac{\nu}{2}}{\bm{f}}=\mathcal{W},$
(8)

for a vector ${\bm{f}}\in\mathbb{R}^{n}$. Thus, analogously to Eq. 6, the
scalar-valued GP on a graph becomes $f_{G}\sim\mathcal{GP}(0,k_{G})$, with
kernel

$k_{G}(p,q)=\sigma^{2}{\bm{u}}(p)\Phi({\bm{\Lambda}})^{-2}{\bm{u}}(q)^{T}$ (9)

where ${\bm{u}}(i),{\bm{u}}(j)$ are the columns of ${\bm{U}}$ corresponding to
nodes $i,j\in V$, respectively.

##  4 Intrinsic representation of vector fields over arbitrary latent
manifolds

We may now construct GPs on unknown manifolds and associated tangent bundles.

Figure 1: Construction of vector-valued Gaussian processes on unknown
manifolds. A Input data consists of samples from a vector field over a latent
manifold. B The manifold is approximated by a proximity graph. Black circles
mark two sample points, $i$ and $j$ and their graph neighbourhood. C The
tangent bundle is a collection of locally Euclidean vector spaces over the
manifold. It is approximated by parallel transport maps between local tangent
space approximations. D The eigenvectors of the connection Laplacian are used
as positional encoding to define the GP that leans the vector fields. E The GP
is evaluated as unseen points to predict the smoothest possible vector field
that is consistent with the training data. We use this GP to accurately
predict singularities, near which sampling is typically sparse.

###  4.1 Vector-valued GPs on unknown manifolds

We consider training data consisting of pairs
$\\{({\bm{x}}_{i},{\bm{v}}_{i})|i=1,\dots,n\\}=(X,V)$, where ${\bm{x}}_{i}$
are samples from a manifold $\mathcal{M}\subset\mathbb{R}^{d}$ and
${\bm{v}}_{i}$ are sampled from the tangent bundle
$\mathcal{TM}=\cup_{i}\mathcal{T}_{i}\mathcal{M}$. If the dimension of the
manifold is $m\leq d$ the tangent spaces
$\mathcal{T}_{i}\mathcal{M}:=\mathcal{T}_{{\bm{x}}_{i}}\mathcal{M}$ anchored
to ${\bm{x}}_{i}$ are isomorphic as a vector space to $\mathbb{R}^{m}$.
Importantly, we assume that both $\mathcal{M}$ and
$\mathcal{T}_{i}\mathcal{M}$ are unknown and seek a GP to provide an implicit
description of the vector field over $\mathcal{M}$ that agrees with the
training set and provides a continuous interpolation at out-of-sample test
points with controllable smoothness properties.

##### Approximating the manifold and the tangent bundle

We first fit a proximity graph $G=(V,E)$ to $X$, defined based on some notion
of similarity (spatial or otherwise) in the data. While $G$ approximates
$\mathcal{M}$ it will not restrict the domain to $V$ as in Borovitskiy et al.
(2021). Then, to approximate $\mathcal{TM}$, note that the tangent spaces do
not come with preferred coordinates. However, being isomorphic to
$\mathbb{R}^{m}$, $\mathcal{T}_{i}\mathcal{M}$ can be parametrised by $m$
orthogonal vectors in the ambient space $\mathbb{R}^{d}$, to form a local
frame, or gauge, $\mathbb{T}_{i}$. To obtain this frame, we take vectors
$\mathbf{e}_{ij}\in\mathbb{R}^{d}$ from $i$ to $N$ nearest nodes $j$, assuming
that they span $\mathcal{T}_{i}\mathcal{M}$222In practice, we pick $N$ closest
nodes to $i$ on the proximity graph where $N$ is a hyperparameter. Larger $N$
increases the overlaps between the nearby tangent spaces. We find that
$N=2D_{ii}$ is often a good compromise between locality and robustness to
noise of the tangent space approximation. and form a matrix by stacking them
column-wise. The left singular vectors corresponding to the $m$ largest
singular values yield the desired frame

$\mathbb{T}_{i}=(\mathbf{t}_{1}^{(1)},\dots\mathbf{t}_{i}^{(m)})\in\mathbb{R}^{d\times
m}.$ (10)

Then, $\hat{{\bm{v}}}_{i}=\mathbb{T}^{T}_{i}{\bm{v}}_{i}$ acts as a projection
of the signal to the tangent space in the $\ell_{2}$ sense.

##### Constraining the vector field over the manifold

Armed with the approximation of $\mathcal{M}$ and $\mathcal{TM}$, by $G$ and
$\\{{\mathbb{T}}_{i}\\}$, we may define the connection Laplacian operator
${\bm{L}}_{c}$ that will regularise the GP’s behaviour by using the global
smoothness of the vector field.

The notion of vector field smoothness is formalised by the parallel transport
map $\mathcal{P}_{j\to i}$ that aligns $\mathcal{T}_{j}\mathcal{M}$ with
$\mathcal{T}_{i}\mathcal{M}$ to allow the comparison of vectors in a common
space. While parallel transport is generally path dependent, we assume that
$i,j$ are close enough such that $\mathcal{P}_{j\to i}$ is the unique smallest
rotation. Indeed, the construction of the nearest neighbour proximity graph
limits pairs $i,j$ to be close in space. This is known as the Lévy-Civita
connection and can be computed as a matrix $\mathbf{O}_{ji}\in O(m)$ in the
orthogonal group (rotation and reflection)

$\mathbf{O}_{ji}=\arg\min_{\mathbf{O}\in
O(m)}||\mathbb{T}_{i}-\mathbb{T}_{j}\mathbf{O}||_{F},$ (11)

where $||\cdot||_{F}$ is the Frobenius norm and can be uniquely computed by
the Kabsch algorithm in $\mathcal{O}(m)$-time (Kabsch, 1976).

Using the parallel transport maps, we may define the connection Laplacian
(Singer & Wu, 2012), a block matrix ${\bm{L}}_{c}\in\mathbb{R}^{nd\times nd}$,
whose $(i,j)$ block entry is given by

${\bm{L}}_{c}(i,j)=\begin{cases}D_{ii}{\bm{I}}_{m\times
m}\;&\text{for}\;i=j\\\ W_{ij}{\bm{O}}_{ij}&\text{for}\;i,j\text{
adjacent}.\end{cases}$ (12)

Let us remark that ${\bm{L}}_{c}$ prescribes the smoothness of the vector
field over an underlying continuous manifold that agrees with the available
training data. In fact, as $n\to\infty$, the eigenvectors of ${\bm{L}}_{c}$
converge to the eigenfunctions of the connection Laplacian over the tangent
bundle (Singer & Wu, 2017), such that the corresponding continuous signal
satisfies the associated vector diffusion process (Berline N., 1996). The
solution of this diffusion process, minimises the vector Dirichlet energy
$\sum_{ij\in E}w_{ij}|{\bm{v}}_{i}-{\bm{O}}_{ij}{\bm{v}}_{j}|^{2}$, which
quantifies the smoothness of the vector field (Knöppel et al., 2015).

##### Vector-field GP on arbitrary latent manifolds

We are now ready to define a GP to regress the vector field over
$\mathcal{M}$. To this end, we consider a positional encoding of points on the
tangent bundle $\mathcal{TM}$ based on the spectrum of the connection
Laplacian, ${\bm{L}}_{c}={\bm{U}}_{c}{\bm{\Lambda}}_{c}{\bm{U}}_{c}^{T}$,
where ${\bm{\Lambda}}_{c},{\bm{U}}_{c}\in\mathbb{R}^{nm\times nm}$. Comparing
with Eq. 7 one notices that each node represents a vector space of dimension
$m$. Using the eigenvectors corresponding to the $k$ smallest eigenvalues, and
for some point ${\bm{x}}$ corresponding to node $i$, we consider the matrix

$({\bm{U}}_{c})_{i}=\sqrt{nm}\begin{pmatrix}u_{im,1}&\dots&u_{im,k}\\\
\vdots&&\vdots\\\
u_{(i+1)m,1}&\dots&u_{(i+1)m,k}\end{pmatrix}\in\mathbb{R}^{m\times k},$ (13)

obtained by keeping $k$ columns of ${\bm{U}}_{c}$. This allows us to define
the positional encoding of ${\bm{x}}$ by mapping the eigencoordinates defined
in the respective tangent spaces using ${\mathbb{T}}_{i}$ back into the
ambient space

${\bm{P}}_{{\bm{x}}}={\mathbb{T}}_{i}({\bm{U}}_{c})_{i}\in\mathbb{R}^{d\times
k}.$ (14)

Then, by analogy to Eqs. 5-6, we define the scalar-valued Matérn GP that
represents the vector field over the manifold as
$f_{\mathcal{TM}}:\mathcal{TM}\to\mathbb{R}^{d}$ with a $\mathbb{R}^{d\times
d}$-valued kernel

$k_{\mathcal{TM}}({\bm{v}},{\bm{v}}^{\prime})=\sigma^{2}{\bm{P}}_{{\bm{v}}}\Phi({\bm{\Lambda}}_{c})^{-2}{\bm{P}}_{{\bm{v}}^{\prime}}^{T},$
(15)

where ${\bm{v}},{\bm{v}}^{\prime}\in\mathcal{TM}$. Note that we recover
$f_{\mathcal{M}}$ for scalar signals ($m=1$), where Laplace-Beltrami operator
equals the connection Laplacian for a trivial line bundle on $\mathcal{M}$.
Therefore, the tangent spaces become trivial (scalar), recovering the well-
known Laplacian eigenmaps
${\bm{P}}_{{\bm{x}}}=\sqrt{n}(u_{i,1},\dots,u_{i,k})$. Likewise,
$k_{\mathcal{TM}}({\bm{v}},{\bm{v}}^{\prime})$ reduces to the scalar-valued
kernel $k_{\mathcal{M}}({\bm{x}},{\bm{x}}^{\prime})$ for
${\bm{x}},{\bm{x}}^{\prime}\in\mathcal{M}$ in Eq. 6. Note that the kernel in
Eq. 15 is a member of a general class defined in Hutchinson et al. (2021,
Proposition 7) arising as an isometric projection of the tangent bundle into
Euclidean space. Our construction makes this projection explicit using the
connection Laplacian positional encoding (Eq. 13) and tangent spaces (Eq. 10).
In Fig. 2, we show that this construction leads to specific performance
advantages in recovering vector field singularities.

###  4.2 Scalable training via inducing point methods

A notable drawback of GPs is their computational inefficiency, largely driven
by the need to compute an $n\times n$ covariance matrix. Therefore, inducing
point methods, which reduce the effective number of data points to a set of
$n^{\prime}<n$ inducing points, have become a mainstay of GPs in practice
(Titsias, 2009). Hutchinson et al. (2021) showed that the inducing point
methods such as that of Titsias (2009) are readily extendible to vector fields
over Riemannian manifolds, provided the covariance matrix of inducing points
is represented in tangent space coordinates. Since the kernel of RVGP is
constructed from a positional encoding using connection Laplacian
eigenvectors, which are expressed in local coordinates, our method is readily
amenable to computation by inducing point methods which we provide an
implementation of.

##  5 Experiments

###  5.1 Manifold-consistent interpolation of vector fields

Figure 2: Superresolution and inpainting. A Sparse samples randomly
distributed over the Stanford bunny and torus surface mesh are interpolated to
a higher resolution ($k=50$). B Ablation studies for the Stanford Bunny,
showing the dependence of mean alignment of the superresolved vector fields in
the test set against the average distance between manifold points $\alpha$
(for $k=50$ fixed) and the number of eigenvectors $k$. C Prediction in a
masked area containing vector field singularity. Predictions using an RBF
kernel with Laplacian eigenvectors as positional encoding produce
discontinuities along the masked boundary (triangle), vectors that protrude
the mesh surface (star) and do not converge to zero magnitude at the
singularity. RVGP using connection Laplacian eigenvectors produces smoothly
varying inpainting ($k=50$).

We expected that RVGP fitted to sparse vector samples over an unknown manifold
would leverage the global regularity of the vector field to provide accurate
out-of-sample predictions. Thus, we conducted two experiments on the Stanford
Bunny and toroidal surface mesh to test the RVGP’s ability to super-resolve
sparse samples and inpaint missing regions containing singularities on diverse
manifold topologies. In both experiments, we first generated a ground truth
vector signal. Next, we drew $n$ evenly sampled anchor points
$\\{({\bm{x}}_{i}\\})$ using the furthest point sampling algorithm (Qi et al.,
2017), which has a parameter $\alpha$ quantifying the average pairwise
distance of points relative to the manifold diameter. We then sampled
corresponding vectors $\\{{\bm{v}}_{i}\\}$ from a uniform distribution on the
sphere $S^{3}$, projected them onto their respective tangent spaces
$\hat{{\bm{v}}}_{i}={\mathbb{T}}_{i}^{T}{\bm{v}}_{i}$ and used the vector heat
method (Sharp et al., 2019) to find the smoothest vector field. Specifically,
concatenating signals as
$\hat{{\bm{v}}}=\mathbin{\|}_{i=0}^{n}\hat{{\bm{v}}}_{i}\in\mathbb{R}^{nd\times
1}$ the vector heat method obtains
$\hat{{\bm{v}}}(\tau)=\hat{{\bm{v}}}(\tau)/u(\tau)$, where
$\hat{{\bm{v}}}(\tau)$ is the solution of the vector heat equation
$\hat{{\bm{v}}}(\tau)=\hat{{\bm{v}}}\exp{\left(-{\bm{L}}_{c}\tau\right)}$ and
$u(\tau)$ is the solution of the scalar heat equation
$u(\tau)=|\hat{{\bm{v}}}|\exp{\left(-{\bm{L}}\tau\right)}$. To obtain a smooth
field, we ran the process until diffusion time $\tau=100$. On surfaces of
genus $0$ this process will lead to at least one singularity.

##### Superresolution

In the first experiment, we asked if RVGP can smoothly interpolate (also known
as super-resolve) the vector field from sparse samples. To this end, we fitted
a graph over the bunny mesh and trained RVGP using vectors over $50\%$ of the
nodes, reserving the remaining $50\%$ for testing. We found that RVGP
predictions were in excellent visual alignment with the training vector field
(Fig. 2A). To investigate the sensitivity of predictions to sampling density
and manifold description, we conducted ablation studies by varying both
$\alpha$ and the number of eigenvectors $k$. As quantified by the inner
product between predicted and test vectors, we found that RVGP produces
accurate alignment with only $10$ eigenvectors over a broad range of sampling
densities (Fig. 2B).

##### Inpainting

In the second experiment, we tested RVGP’s ability to inpaint whole vector
field regions containing singularities. This experiment is more challenging
than classical inpainting because it requires our method to infer the
smoothest topologically consistent vector field. To this end, we masked off
the vortex singularity and used the remaining points to train RVGP consisting
of a positional encoding of the manifold based on graph Laplacian
eigenvectors. For comparison, we also trained a radial basis function kernel
that treats vector entries independently. We found that vectors predicted by
RVGP closely followed the mesh surface, aligned with the training set on the
mask boundary and smoothly resolved the singularity by gradually reducing the
vector amplitudes to zero. By contrast, the RBF kernel did not produce vectors
with smooth transitions on the mask boundary and often protruded from the mesh
surface. This example thus shows that the connection Laplacian positional
encoding with vector field kernel provides sufficient regularity to learn
vector fields over complex shapes.

###  5.2 Superresolution of EEG data

Finally, as a biologically and clinically relevant use case, we applied RVGP
to superresolve electroencephalography (EEG) recordings from human subjects.
EEG recordings measure spatiotemporal wave patterns associated with neural
dynamics, which play a fundamental role in human behaviour (Sato et al., 2012;
Xu et al., 2023). Accurately resolving these dynamics requires high-density
EEG setups in excess of 200 channels Robinson et al. (2017); Seeber et al.
(2019); Siclari et al. (2018). However, due to long setup times during which
signal quality can rapidly degrade, experimentalists and clinicians commonly
resort to low-density recordings with 32 or 64 channels (Chu, 2015).

Thus, we asked whether superresolving low-density 64-channel recordings (Fig.
3A) using RVGP into high-density 256-channel recordings (Fig. 3C) can
facilitate biological discovery and clinical diagnostics. As a ground truth,
we collected 256-channel EEG recordings of resting-state brain activity from
33 Alzheimer’s patients (AD) and 28 age-matched healthy controls (see Appendix
A.1). We focused on low-frequency alpha waves spanning the 8-15 Hz range,
which represent the dominant rhythm at rest (Berger, 1934) and exhibits
impaired dynamics in AD (Moretti et al., 2004; Besthorn et al., 1994; Dauwels
et al., 2010). After preprocessing of EEG time series (see Appendix A.1), we
used a triangulated mesh of the head surface, with vertices given by 256
electrode locations, and finite differencing to compute the vector field at
mesh vertices (Fig. 3B,C (ground truth), Appendix A.2). Then, we constructed
RVGP kernels using the connection Laplacian eigenvectors derived from a
$k$-nearest neighbour graph ($k=5$) fit to mesh vertices. Finally, we trained
RVGP using vectors at 64 training vertices (Fig. 3B) and then used it to infer
the vector field at the remaining 192 test vertices (Fig. 3C). As a benchmark,
we also used a channel-wise linear and spline interpolation of the vector
field.

Figure 3: Reconstruction of spatiotemporal wave patterns in the human EEG. A
Snapshot of an alpha wave pattern (8-15 Hz) recorded on low density (64
channel) EEG from a healthy subject projected in two dimensions. B Phase field
of an alpha wave. Vector field denotes the the spatial gradient of the voltage
signal. C Ground-truth and reconstructed high-density phase field (256
channel) using RVGP, linear and spline interpolation. Streamlines, computed
based on the vector field, highlight features of the phase field. RVGP
significantly better preserves singularities, i.e., sources, sinks and
vortices. D Reconstruction accuracy, measured by the preservation of
singularities. E Receiver operating characteristic (ROC) for binary
classification of patients with Alzheimer’s disease against healthy controls
using a linear support vector machine trained on the divergence and vorticity
fields.

To assess the quality of reconstructions, we computed the divergence and curl
of the predicted EEG vector field and computed the mean absolute error (MAE)
relative to the ground truth 256-channel EEG. The divergence and curl have
previously been used to identify singularities in neural wave patterns from
human neuroimaging recordings and have been linked to cognitive function and
behaviour (Roberts et al., 2019; Xu et al., 2023). We found that the RVGP
reconstruction closely approximated the divergence and curl of the high-
density EEG (Fig. 3C). For visual comparison, we show a single time snapshot
with characteristic vector field singularities such as sources, sinks and
vortices. In particular, RVGP reconstruction substantially better preserved
these singularities than linear interpolation, which had the tendency to
underfit the vector field, and spline interpolation, which introduced spurious
local structures. This observation is corroborated by significantly lower
angular error and curl-divergence MAE for all subjects for RVGP compared with
benchmark ($n=61$, Fig. 3D).

Given the superior reconstruction accuracy of our approach, we sought to
determine if it could enhance the classification accuracy between patients
with Alzheimer’s Disease (AD) and healthy controls. Contemporary diagnostic
methods for AD are costly, invasive, and lengthy in procedure (Zetterberg &
Bendlin, 2021). To underscore the potential clinical significance and
scalability of our methodology, we employed a linear support vector machine
for the classification of AD patients versus age-matched controls. This
classification was based on the reconstructed divergence and curl fields
derived from a brief, one-minute resting-state low-density EEG recording – a
procedure that is feasibly integrated into clinical settings due to its non-
invasive nature and cost-efficiency. Our results show significantly higher
than state-of-the-art classification accuracy, approaching that derived from
the ground truth high-density EEG signal (Figure 3E).

##  6 Conclusion

In this study, we introduced RVGP, a novel extension of Gaussian processes
designed to model vector fields on latent Riemannian manifolds. Utilising the
spectrum of the connection Laplacian operator, RVGP intrinsically captures the
manifold’s geometry and topology and the vector field’s smoothness. This
enables the method to learn global patterns while preserving singularities,
thereby filling a significant gap in existing approaches that are limited to
known, analytically tractable manifolds. A key strength of RVGP is its data-
driven, intrinsic construction via a proximity graph, which enhances its
practical utility by making it highly applicable to real-world datasets where
the manifold is often not explicitly known. Demonstrated across diverse
scientific domains such as neuroscience and geometric data analysis, RVGP not
only advances the field of geometrically-informed probabilistic modelling but
also offers a statistical tool for various high-impact applications, including
clinical neuroscience, that is robust to sampling density and noisy data.

#### Author Contributions

AG and PV devised the RVGP method and developed the theoretical foundations
with contributions from RLP. AG and RLP wrote the RVGP code. AG, RLP and MVC
designed the experiments. MVC and NG conceptualised the EEG preprocessing
pipeline. MD, EM, DS, and PAM collected the EEG data. MVC, RLP and NG
performed the analysis of the EEG data with contributions from AG. AG, RLP,
MVC and PV wrote the paper.

#### Acknowledgments

We thank Nicolas Aspert for the critical reading and feedback on the
manuscript. AG acknowledges support from an HFSP Cross-disciplinary
Postdoctoral Fellowship (LT000669/2020-C). PAM acknowledges funding from the
NIHR Biomedical Research Centre at Imperial College London. RLP acknowledges
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project-
ID 424778381-TRR 295.

### Code availability

The Python code to reproduce our results is openly available at
https://github.com/agosztolai/RVGP.

## References

  * Barbero et al. (2022) Federico Barbero, Cristian Bodnar, Haitz Sáez de Ocáriz Borde, Michael Bronstein, Petar Veličković, and Pietro Liò.  Sheaf neural networks with connection laplacians, 2022.
  * Belkin & Niyogi (2003) Mikhail Belkin and Partha Niyogi.  Laplacian eigenmaps for dimensionality reduction and data representation.  _Neural Comput._ , 15(6):1373–1396, 2003.  ISSN 08997667.  doi: 10.1162/089976603321780317.
  * Berger (1934) Hans Berger.  Über das elektrenkephalogramm des menschen.  _DMW-Deutsche Medizinische Wochenschrift_ , 60(51):1947–1949, 1934.
  * Berline N. (1996) Vergne M. Berline N., Getzler E.  _Heat kernels and Dirac operators_.  Springer, 2nd. edition, 1996.
  * Besthorn et al. (1994) Christoph Besthorn, Hans Förstl, Claudia Geiger-Kabisch, Heribert Sattel, Theo Gasser, and Ursula Schreiter-Gasser.  Eeg coherence in alzheimer disease.  _Electroencephalography and clinical neurophysiology_ , 90(3):242–245, 1994.
  * Bodnar et al. (2022) Cristian Bodnar, Francesco Di Giovanni, Benjamin Paul Chamberlain, Pietro Liò, and Michael M. Bronstein.  Neural sheaf diffusion: A topological perspective on heterophily and oversmoothing in GNNs.  In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho (eds.), _Advances in Neural Information Processing Systems_ , 2022.
  * Borovitskiy et al. (2020) Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Deisenroth.  Matérn gaussian processes on riemannian manifolds.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (eds.), _Advances in Neural Information Processing Systems_ , volume 33, pp. 12426–12437. Curran Associates, Inc., 2020.
  * Borovitskiy et al. (2021) Viacheslav Borovitskiy, Iskander Azangulov, Alexander Terenin, Peter Mostowsky, Marc Deisenroth, and Nicolas Durrande.  Matérn gaussian processes on graphs.  In Arindam Banerjee and Kenji Fukumizu (eds.), _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , volume 130 of _Proceedings of Machine Learning Research_ , pp. 2593–2601. PMLR, 13–15 Apr 2021.
  * Bronstein et al. (2017) Michael M. Bronstein, Joan Bruna, Yann Lecun, Arthur Szlam, and Pierre Vandergheynst.  Geometric Deep Learning: Going beyond Euclidean data.  _IEEE Signal Process. Mag._ , 34(4):18–42, 2017\.  ISSN 10535888.  doi: 10.1109/MSP.2017.2693418.
  * Budninskiy et al. (2019) Max Budninskiy, Gloria Yin, Leman Feng, Yiying Tong, and Mathieu Desbrun.  Parallel transport unfolding: A connection-based manifold learning approach.  _SIAM Journal on Applied Algebra and Geometry_ , 3(2):266–291, 2019.  doi: 10.1137/18M1196133.
  * Chu (2015) Catherine J Chu.  High density eeg—what do we have to lose?  _Clinical neurophysiology: official journal of the International Federation of Clinical Neurophysiology_ , 126(3):433, 2015\.
  * Chung (1997) Fan Chung.  _Spectral Graph Theocy_ , volume 92 of _American Mathematical Soc._ American Mathematical Soc., 1997.  ISBN 0821803158.
  * Coifman et al. (2005) Ronald R Coifman, Stéphane Lafon, Anne B Lee, Mauro Maggioni, Boaz Nadler, Frederick Warner, and Steven W Zucker.  Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps.  _Proc. Natl. Acad. Sci. U. S. A._ , 102(21):7426–7431, 2005.  doi: 10.1073/pnas.0500334102.
  * Dauwels et al. (2010) Justin Dauwels, Francois Vialatte, and Andrzej Cichocki.  Diagnosis of alzheimer’s disease from eeg signals: where are we standing?  _Current Alzheimer Research_ , 7(6):487–505, 2010\.
  * Defferrard et al. (2016) Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst.  Convolutional Neural Networks on Graphs with.  In _NIPS’ 2016_ , pp. 3844–3852, 2016.  ISBN 9781510838819.
  * Delorme & Makeig (2004) Arnaud Delorme and Scott Makeig.  Eeglab: an open source toolbox for analysis of single-trial eeg dynamics including independent component analysis.  _Journal of neuroscience methods_ , 134(1):9–21, 2004.
  * Duncker & Sahani (2018) Lea Duncker and Maneesh Sahani.  Temporal alignment and latent Gaussian process factor inference in population spike trains.  _Adv. Neural Inf. Process. Syst._ , 2018-December(NeurIPS):10445–10455, 2018.  ISSN 10495258.
  * Fefferman et al. (2016) Charles Fefferman, Sanjoy Mitter, and Hariharan Narayanan.  Testing the manifold hypothesis.  _Journal of the American Mathematical Society_ , 29(4):983–1049, 2016.
  * Gardner et al. (2022) Richard J. Gardner, Erik Hermansen, Marius Pachitariu, Yoram Burak, Nils A. Baas, Benjamin A. Dunn, May Britt Moser, and Edvard I. Moser.  Toroidal topology of population activity in grid cells.  _Nature_ , 602(7895):123–128, 2022.  doi: 10.1038/s41586-021-04268-7.
  * Gosztolai & Arnaudon (2021) Adam Gosztolai and Alexis Arnaudon.  Unfolding the multiscale structure of networks with dynamical Ollivier-Ricci curvature.  _Nat. Commun._ , 12(1):1–11, 2021.  ISSN 20411723.  doi: 10.1038/s41467-021-24884-1.
  * Gosztolai et al. (2021) Adam Gosztolai, Semih Günel, Victor Lobato-Ríos, Marco Pietro Abrate, Daniel Morales, Helge Rhodin, Pascal Fua, and Pavan Ramdya.  LiftPose3D, a deep learning-based approach for transforming two-dimensional to three-dimensional poses in laboratory animals.  _Nat. Methods_ , 18(8):975–981, 2021.  ISSN 1548-7091.  doi: 10.1038/s41592-021-01226-z.
  * Gosztolai et al. (2023) Adam Gosztolai, Robert L. Peach, Alexis Arnaudon, Mauricio Barahona, and Pierre Vandergheynst.  Interpretable statistical representations of neural population dynamics and geometry.  _arXiv:2304.03376_ , 2023.
  * Grattarola & Vandergheynst (2022) Daniele Grattarola and Pierre Vandergheynst.  Generalised implicit neural representations.  _Advances in Neural Information Processing Systems_ , 35:30446–30458, 2022.
  * Hutchinson et al. (2021) Michael Hutchinson, Alexander Terenin, Viacheslav Borovitskiy, So Takao, Yee Teh, and Marc Deisenroth.  Vector-valued gaussian processes on riemannian manifolds via gauge independent projected kernels.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan (eds.), _Advances in Neural Information Processing Systems_ , volume 34, pp. 17160–17169. Curran Associates, Inc., 2021.
  * Illoul & Lorong (2011) Lounes Illoul and Philippe Lorong.  On some aspects of the cnem implementation in 3d in order to simulate high speed machining or shearing.  _Computers & Structures_, 89(11-12):940–958, 2011.
  * Jensen et al. (2020) Kristopher T. Jensen, Ta Chu Kao, Marco Tripodi, and Guillaume Hennequin.  Manifold GPLVMs for discovering non-Euclidean latent structure in neural data.  _Adv. Neural Inf. Process. Syst._ , 2020-Decem:1–22, 2020\.  ISSN 10495258.
  * Kabsch (1976) Wolfgang Kabsch.  A solution for the best rotation to relate two sets of vectors.  _Acta Crystallographica Section A_ , 32(5):922–923, September 1976.  doi: 10.1107/S0567739476001873.
  * Khona & Fiete (2022) Mikail Khona and Ila R. Fiete.  Attractor and integrator networks in the brain.  _Nat. Rev. Neurosci._ , 23(12):744–766, 2022\.  doi: 10.1038/s41583-022-00642-0.
  * Kipf & Welling (2017) Thomas N. Kipf and Max Welling.  Semi-Supervised Classification with Graph Convolutional Networks.  _ICLR_ , 2017.  ISSN 10963626.  URL arXiv.org.
  * Knöppel et al. (2013) Felix Knöppel, Keenan Crane, Ulrich Pinkall, and Peter Schröder.  Globally optimal direction fields.  _ACM Trans. Graph._ , 32(4), 2013.  doi: 10.1145/2461912.2462005.
  * Knöppel et al. (2015) Felix Knöppel, Keenan Crane, Ulrich Pinkall, and Peter Schröder.  Stripe patterns on surfaces.  _ACM Trans. Graph._ , 34(4), 2015.  ISSN 15577368.  doi: 10.1145/2767000.
  * Koestler et al. (2022) Lukas Koestler, Daniel Grittner, Michael Moeller, Daniel Cremers, and Zorah Lähner.  Intrinsic Neural Fields: Learning Functions on Manifolds.  _Lect. Notes Comput. Sci. (including Subser. Lect. Notes Artif. Intell. Lect. Notes Bioinformatics)_ , 13662 LNCS:622–639, 2022.  ISSN 16113349.  doi: 10.1007/978-3-031-20086-1˙36.
  * La Manno et al. (2018) Gioele La Manno, Ruslan Soldatov, Amit Zeisel, Emelie Braun, Hannah Hochgerner, Viktor Petukhov, Katja Lidschreiber, Maria E. Kastriti, Peter Lönnerberg, Alessandro Furlan, Jean Fan, Lars E. Borm, Zehua Liu, David van Bruggen, Jimin Guo, Xiaoling He, Roger Barker, Erik Sundström, Gonçalo Castelo-Branco, Patrick Cramer, Igor Adameyko, Sten Linnarsson, and Peter V. Kharchenko.  RNA velocity of single cells.  _Nature_ , 560(7719):494–498, 2018.  ISSN 14764687.  doi: 10.1038/s41586-018-0414-6.
  * Lipman (2021) Yaron Lipman.  Phase Transitions, Distance Functions, and Implicit Neural Representations.  _Proc. Mach. Learn. Res._ , 139:6702–6712, 2021.  ISSN 26403498.
  * Mallasto & Feragen (2018) Anton Mallasto and Aasa Feragen.  Wrapped gaussian process regression on riemannian manifolds.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2018.
  * Mallasto et al. (2020) Anton Mallasto, Søren Hauberg, and Aasa Feragen.  Probabilistic Riemannian submanifold learning with wrapped Gaussian process latent variable models.  _AISTATS 2019 - 22nd Int. Conf. Artif. Intell. Stat._ , 89, 2020.
  * Mildenhall et al. (2020) Ben Mildenhall, Pratul P. Srinivasan, Matthew Tancik, Jonathan T. Barron, Ravi Ramamoorthi, and Ren Ng.  Nerf: Representing scenes as neural radiance fields for view synthesis.  In _ECCV_ , 2020.
  * Moretti et al. (2004) Davide V Moretti, Claudio Babiloni, Giuliano Binetti, Emanuele Cassetta, Gloria Dal Forno, Florinda Ferreric, Raffaele Ferri, Bartolo Lanuzza, Carlo Miniussi, Flavio Nobili, et al.  Individual analysis of eeg frequency and band power in mild alzheimer’s disease.  _Clinical Neurophysiology_ , 115(2):299–308, 2004\.
  * Ortega et al. (2018) Antonio Ortega, Pascal Frossard, Jelena Kovacevic, Jose M.F. Moura, and Pierre Vandergheynst.  Graph Signal Processing: Overview, Challenges, and Applications.  _Proc. IEEE_ , 106(5):808–828, 2018.  ISSN 15582256.  doi: 10.1109/JPROC.2018.2820126.
  * Peach et al. (2020) Robert L Peach, Alexis Arnaudon, and Mauricio Barahona.  Semi-supervised classification on graphs using explicit diffusion dynamics.  _Foundations of Data Science_ , 2(1):19–33, 2020\.
  * Perrin et al. (1989) François Perrin, Jacques Pernier, Olivier Bertrand, and Jean Francois Echallier.  Spherical splines for scalp potential and current density mapping.  _Electroencephalography and clinical neurophysiology_ , 72(2):184–187, 1989.
  * Qi et al. (2017) Charles R. Qi, Hao Su, Kaichun Mo, and Leonidas J. Guibas.  PointNet: Deep learning on point sets for 3D classification and segmentation.  _Proc. - 30th IEEE Conf. Comput. Vis. Pattern Recognition, CVPR 2017_ , 2017-January:77–85, 2017.  doi: 10.1109/CVPR.2017.16.
  * Rasmussen & Williams (2006) Carl Edward Rasmussen and Christopher K. I. Williams.  _Gaussian Processes for Machine Learning_.  2006\.  ISBN 026218253X.
  * Roberts et al. (2019) James A Roberts, Leonardo L Gollo, Romesh G Abeysuriya, Gloria Roberts, Philip B Mitchell, Mark W Woolrich, and Michael Breakspear.  Metastable brain waves.  _Nature communications_ , 10(1):1056, 2019.
  * Robinson et al. (2017) Amanda K Robinson, Praveen Venkatesh, Matthew J Boring, Michael J Tarr, Pulkit Grover, and Marlene Behrmann.  Very high density eeg elucidates spatiotemporal aspects of early visual processing.  _Scientific reports_ , 7(1):16248, 2017.
  * Sato et al. (2012) Tatsuo K Sato, Ian Nauhaus, and Matteo Carandini.  Traveling waves in visual cortex.  _Neuron_ , 75(2):218–229, 2012.
  * Seeber et al. (2019) Martin Seeber, Lucia-Manuela Cantonas, Mauritius Hoevels, Thibaut Sesia, Veerle Visser-Vandewalle, and Christoph M Michel.  Subcortical electrophysiological activity is detectable with high-density eeg source imaging.  _Nature communications_ , 10(1):753, 2019.
  * Sharp et al. (2019) Nicholas Sharp, Yousuf Soliman, and Keenan Crane.  The vector heat method.  _ACM Trans. Graph._ , 38(3), 2019.  doi: 10.1145/3243651.
  * Siclari et al. (2018) Francesca Siclari, Giulio Bernardi, Jacinthe Cataldi, and Giulio Tononi.  Dreaming in nrem sleep: a high-density eeg study of slow waves and spindles.  _Journal of Neuroscience_ , 38(43):9175–9185, 2018.
  * Singer & Wu (2012) Amit Singer and Hau Tieng Wu.  Vector diffusion maps and the connection Laplacian.  _Commun. Pure Appl. Math._ , 65(8):1067–1144, 2012.  doi: 10.1002/cpa.21395.
  * Singer & Wu (2017) Amit Singer and Hau Tieng Wu.  Spectral convergence of the connection Laplacian from random samples.  _Inf. Inference_ , 6(1):58–123, 2017.  doi: 10.1093/imaiai/iaw016.
  * Sitzmann et al. (2020) Vincent Sitzmann, Julien N.P. Martel, Alexander W. Bergman, David B. Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  _Adv. Neural Inf. Process. Syst._ , 2020-December(NeurIPS):1–12, 2020.  ISSN 10495258.
  * Solin & Särkkä (2020) Arno Solin and Simo Särkkä.  Hilbert space methods for reduced-rank Gaussian process regression.  _Stat. Comput._ , 30(2):419–446, 2020.  ISSN 15731375.  doi: 10.1007/s11222-019-09886-w.
  * Sussillo & Barak (2013) David Sussillo and Omri Barak.  Opening the black box: low-dimensional dynamics in high-dimensional recurrent neural networks.  _Neural Comput._ , 25(3):626–649, 2013.  doi: 10.1162/NECO˙a˙00409.
  * Taubin (1995) Gabriel Taubin.  A signal processing approach to fair surface design.  In _Proceedings of the 22nd Annual Conference on Computer Graphics and Interactive Techniques_ , SIGGRAPH ’95, pp. 351–358, New York, NY, USA, 1995. Association for Computing Machinery.  ISBN 0897917014.  doi: 10.1145/218380.218473.
  * Tenenbaum (2000) Joshua B Tenenbaum.  A Global Geometric Framework for Nonlinear Dimensionality Reduction.  _Science (80-. )._ , 290(5500):2319–2323, 2000\.  doi: 10.1126/science.290.5500.2319.
  * Titsias (2009) Michalis Titsias.  Variational learning of inducing variables in sparse gaussian processes.  In David van Dyk and Max Welling (eds.), _Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics_ , volume 5 of _Proceedings of Machine Learning Research_ , pp. 567–574, Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA, 16–18 Apr 2009\. PMLR.
  * Whittle (1963) Peter Whittle.  Stochastic-processes in several dimensions.  _Bulletin of the International Statistical Institute_ , 40(2):974–994, 1963.
  * Wilson et al. (2021) James T. Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Peter Deisenroth.  Pathwise conditioning of gaussian processes.  _J. Mach. Learn. Res._ , 22:1–47, 2021.
  * Xu et al. (2023) Yiben Xu, Xian Long, Jianfeng Feng, and Pulin Gong.  Interacting spiral wave patterns underlie complex brain dynamics and are related to cognitive processing.  _Nature Human Behaviour_ , pp. 1–20, 2023.
  * Yu et al. (2009) Byron M. Yu, John P. Cunningham, Gopal Santhanam, Stephen I. Ryu, Krishna V. Shenoy, and Maneesh Sahani.  Gaussian-process factor analysis for low-dimensional single-trial analysis of neural population activity.  _J. Neurophysiol._ , 102(1):614–635, 2009.  ISSN 15221598.  doi: 10.1152/jn.90941.2008.
  * Zetterberg & Bendlin (2021) Henrik Zetterberg and Barbara B Bendlin.  Biomarkers for alzheimer’s disease—preparing for a new era of disease-modifying therapies.  _Molecular psychiatry_ , 26(1):296–308, 2021\.

##  Appendix A EEG analysis

The following analyses were scripted in MATLAB R2023a:

###  A.1 Experimental data, and pre-processing

Experimental data: We analysed 10-minute resting state EEG recordings from a
cohort of 61 subjects containing 33 patients with a clinical diagnosis of AD
and 28 age-matched healthy control subjects. Data was acquired from the
clinical imaging facility at Hammersmith Hospital in London. The AD group
consisted of 15 Females and 18 Males with a mean age of 76 years (std = 7.8
years), healthy control group contained 15 Females and 13 Males with a mean
age of 77 (std = 4.9 years). 10 participants were excluded from the original
full cohort of 71 subjects due to poor data quality. During pre-processing,
the EEG was bandpass filtered from 2-40 Hz with a second-order Butterworth
filter and downsampled to 250Hz. Channels exceeding a kurtosis threshold of 3
were rejected before average re-referencing. We performed an independent
component analysis with the Picard algorithm and applied an automatic artefact
rejection algorithm in eeglab with a maximum of 5 % of components removed from
each dataset. The pre-processed EEG was then filtered in the alpha band
(8-15Hz). The low-density EEG consisted of 64 channels sampled evenly across
the scalp. These channels were then used to reconstruct the full 256-channel
high-density recording.

###  A.2 Wave velocities

We calculated the velocity vector field at each time point using a method
similar to Roberts et al. (2019). We compute the instantaneous phase at each
channel using the Hilbert transform and estimate the wave velocity $v$ from
the spatial and temporal derivatives of the unwrapped phase $\phi(x,y,z,t)$,
as $v=-||\partial t\partial\phi||/||\nabla\phi||_{2}$ implemented using the
constrained natural element method (CNEM) (Illoul & Lorong, 2011). CNEM is a
mesh-free calculus method for solving partial differential equations that
avoids artefacts that can arise from mesh-based interpolation of nodes
positioned on the outer boundary of the brain’s convex hull.

###  A.3  Linear and spherical spline interpolation

Spherical spline interpolation was carried out on the raw EEG for each subject
using 64 electrodes sampled evenly across the scalp and implemented in eeglab
(Delorme & Makeig, 2004). Each channel was mapped onto a unit sphere, and the
electrical potential at the missing channel locations was interpolated using
the spline function, which weights the contributions of the neighbouring
electrode based on its distance to the interpolation point (Perrin et al.,
1989).

Linear interpolation was performed on the 64-node 3-dimensional phase flow
field. The gradient vector for each node was projected onto the tangent plane
prior to linear interpolation.

###  A.4 Feature Extraction

To determine the behaviourally relevant regions of the cortex, we computed the
time-resolved divergence and vorticity fields for each subject. Next, we
estimated the probability density of sources, sinks and spirals by computing
the mean positive/negative vorticity and mean positive/negative divergence
fields thresholded at $1$ and $<-1$ averaged over time for each node. This
produced four cortical probability maps per subject; mean source probability
(positive divergence, $D>1$), mean sink probability (negative divergence
$D<-1$), clockwise spiral probability (positive vorticity, $C>1$),
anticlockwise spiral probability (negative vorticity, $C<-1$). We extracted
two predictors from these four probability maps. The first predictor was the
probability ratio of sources to sinks at a given node. The second predictor
was the probability ratio of clockwise to anti-clockwise spirals at a given
node:

$R_{n}=\log_{10}\left(\frac{P_{n}(A)}{P_{n}(B)}\right)$ (16)

where $P_{n}(A)$ represent either the probability of a source or a clockwise
spiral at node $n$ and $P_{n}(B)$ represents the probability of a sink or an
anti-clockwise at node $n$. We ran a group-level cluster permutation analysis
to identify nodes where the source-to-sink ratio or clockwise-to-anti-
clockwise spiral ratio differed between the patients with AD and healthy
controls in the reference EEG (alpha = 0.05, 50000 permutation). The three
clusters with the lowest p-value were then used as seed regions where we
extracted the source-to-sink probability ratio and clockwise to anti-clockwise
spiral probability ratio in the interpolated EEG for each reconstruction
method (RVGP, spline and linear interpolation).

###  A.5 Binary Classification

To evaluate whether the RVGP reconstruction contained relevant information
about cognitive function, we tested whether healthy controls and Alzheimer’s
disease patients could be accurately classified from the ratio of source to
sinks and the ratio of clockwise to anti-clockwise spirals. We tested the
classification accuracy using 3 seed brain regions and compared the accuracy
from each reconstructed approach (RVGP, spline, linear interpolation) against
the high-density reference EEG. For binary classification, we used a linear
support vector machine and computed the classification accuracy and ROC curves
for each approach separately using 10-fold cross-validation. We trained two
one-vs-all classifiers to binarily categorize samples pooled from all 3 brain
regions and across subjects into each class (AD or healthy controls). For each
fold of training, we fit the posterior distribution of the scores to 90% of
the data and used this to determine the probability of each sample belonging
to each class in the test set (10%). To handle class imbalance, we repeated
the classification 500 times, with each iteration trained on a random subset
of 100 AD and 100 healthy control samples. The AUC and ROC reported in the
results were calculated from the distribution estimated over 500 iterations.

###  A.6 Visualisation

For visualisation, the $x,y,z$ coordinates for each channel location were
flattened onto a 2D grid using multi-dimensional scaling to most accurately
preserve the local distances between nodes across the scalp. For plotting, we
use the 2D tangent vectors of the phase field and streamlines computed using
the streamslice function in MATLAB.

###  A.7 Data Collection

Data was collected in the UKDRI Care Research & Technology Centre at the
Micheal Uren Hub in London from 2022-2023. AD exclusion criteria were a
previous history of head injury or any other neurological condition. All
patients had a clinical diagnosis of AD, with one subject retrospectively
excluded after a change in diagnosis to frontal-temporal lobe dementia (FTLD).
Subjects on AD medication were asked to abstain on the day of testing and to
avoid caffeine. 11 subjects were excluded from EEG analysis due to low signal
quality leaving 61 subjects total (33 AD, 28 controls).