# Manifold GCN: Diffusion-based Convolutional Neural Network for Manifold-
valued Graphs

Martin Hanik, Gabriele Steidl, and Christoph von Tycowicz M. Hanik is with the
Freie Universität Berlin, Institute of Mathematics, Arnimallee 14, 14195
Berlin, Germany. E-mail<EMAIL_ADDRESS>Steidl is with the Technische
Universität Berlin, Institute of Mathematics, Straße des 17. Juni 136, 10623
Berlin, Germany. E-mail<EMAIL_ADDRESS>v. Tycowicz is with the
Zuse Institute Berlin, Takustr. 7, 14195 Berlin, Germany. E-mail:
<EMAIL_ADDRESS>author: M. Hanik.

###### Abstract

We propose two graph neural network layers for graphs with features in a
Riemannian manifold. First, based on a manifold-valued graph diffusion
equation, we construct a diffusion layer that can be applied to an arbitrary
number of nodes and graph connectivity patterns. Second, we model a tangent
multilayer perceptron by transferring ideas from the vector neuron framework
to our general setting. Both layers are equivariant with respect to node
permutations and isometries of the feature manifold. These properties have
been shown to lead to a beneficial inductive bias in many deep learning tasks.
Numerical examples on synthetic data as well as on triangle meshes of the
right hippocampus to classify Alzheimer’s disease demonstrate the very good
performance of our layers.

###### Index Terms:

graph neural networks, manifold-valued features, diffusion, hyperbolic
embeddings, shape classification

##  I Introduction

Graph neural networks (GNNs) have gained widespread popularity for the
analysis of graph-structured data. They have found applications in various
areas, such as bioinformatics [1, 2, 3], physics [4, 5, 6], and social
sciences [7, 8, 9, 10]. Important information that GNNs utilize are the
features attributed to the graph’s nodes. Here, the vast majority of the GNNs
relies on a Euclidean feature space. However, there are many applications
where the data lies in a non-flat manifold; for example, on the sphere [11,
12, 13], the special orthogonal group SO(3) [14, 15], the space of symmetric
positive definite (SPD) matrices [16, 17, 18, 19, 20, 21, 22, 23], the Stiefel
manifold [24, 25, 26, 27], and the Grassmannian manifold [28, 29, 30]. Also,
other Lie Group-valued features are of interest [31, 32, 33, 34, 35, 36]. So
far, there exist only networks that are designed for very specific types of
manifolds, such as (products of) spaces of constant curvature [37, 38, 39, 40]
or spaces of SPD matrices [41, 22, 42]. These architectures cannot deal with
general manifold-valued features.

In this paper, we present a GNN architecture that can work with data from any
Riemannian manifold for which geodesics can be efficiently approximated. The
core element is a convolution-type layer that discretizes graph diffusion in
manifolds. We thereby utilize a link between convolution and diffusion
processes that has already successfully been used by GNNs for Euclidean data
[43, 44, 45]. Combined with an explicit solution scheme for the diffusion
equation, the _new diffusion layer_ exhibits the same local aggregation of
features that makes convolution layers so successful [46]. Furthermore, it is
equivariant under node permutations as well as isometric transformations of
the feature space and, thus, preserves the symmetries inherent to both the
base and feature domains. Both properties hold at least approximately for many
relationships that deep learning techniques intend to preserve [46]. Indeed,
with our layer, we can build networks that are not disturbed, for example, by
permutations of a graph’s nodes (symmetry of the base domain) or different
feature representations (symmetry of the feature/signal domain). To the best
of our knowledge, our layers are the first GNN layers for manifolds that
exhibit equivariance under isometries.

Additionally, we propose a novel tangent multilayer perceptron that can be
seen as a generalization of a fully-connected multilayer perceptron to
manifolds, including nonlinearities between layers. It shares the symmetries
of the diffusion layer.

We test a network incorporating our novel layers on two different
classification tasks: a benchmark on synthetic graphs and the classification
of Alzheimer’s disease from triangle meshes of the right hippocampus. Our
method outperforms its competition in both cases. We find that our layers lead
to fast learning, even with smaller data sets. The code to reproduce the
experiments will be made public upon acceptance.

_Outline of the paper._ We start by recalling related work in Section II.
Then, in Section III, we provide the necessary background on Riemannian
manifolds and graph Laplacians with manifold-valued features. In Section IV,
we introduce our novel diffusion layer and prove its equivariance properties.
Our second layer, the tangent multilayer perceptron, is constructed in Section
V and shows the same desirable equivariance behavior as the diffusion layer.
Both layers are then combined within a generic graph convolutional neural
network (GCN) block in Section VI. Section VII contains our numerical results.
Finally, a summary and conclusions for future work are given in Section VIII.
The appendices contain further theoretical results on the diffusion layer.

##  II Related Work

In this section, we review related works and highlight the differences to our
contribution. We distinguish between two types of architecture: We start with
an overview of GNNs for graphs with Euclidean features that discretize a
continuous flow in the feature space. We refer to [47] for an overview of the
full landscape of Euclidean GNNs. Then, we collect deep learning architectures
that can handle manifold-valued features.

###  II-A Flow-based GNNs in the Euclidean Space

Here, two approaches were mainly applied, namely diffusion networks and neural
ordinary differential equations.

####  II-A1 Diffusion networks

Diffusion was utilized to construct powerful (graph) convolutional neural
networks in [43, 45]. Several papers improved and modified the basic idea.
While adaptive support for diffusion-convolutions was proposed in [48],
implicit nonlinear diffusion was introduced to capture long-range interactions
in [49]. A polynomial filter on the graph Laplacian was used in [50] that
corresponds to a multi-scale diffusion equation.

Based on a graph diffusion equation, a broad class of well-performing GNNs was
identified in [51], where it was shown in [52] that adding a source term to
the equation is useful in certain scenarios. A nonlinear extension of the
graph diffusion equation was considered in [53], and many known GNN
architectures are instances of the resulting class of GNNs. The framework of
[53] also allows the use of non-Euclidean distances for data encoding.
However, in contrast to our method, this architecture does not treat the
features as being intrinsic to a manifold. A different line of work uses
diffusion processes in the context of cellular sheaf theory for graph learning
[54, 55].

Finally, a diffusion-based architecture, especially for graphs that discretize
2-manifolds (for example, triangular surface meshes), was proposed in [56].

####  II-A2 Neural ordinary differential equations

An approach similar to diffusion-based GNNs is the neural ordinary
differential equations (neural ODE) framework [57]. Instead of diffusion-based
methods, which treat GNNs as particular instances of the discretization of a
partial differential equation (in both space and time) [53], it discretizes an
underlying ODE. The idea was transferred to graph learning in [58, 59].

###  II-B Deep Neural Networks for Manifold-valued Signals

Here we distinguish between three approaches.

####  II-B1 Networks for manifold-valued grids

ManifoldNet [60] is a convolutional neural network that can take manifold-
valued images as input. Convolutions in manifolds are generalized through
weighted averaging. However, in contrast to our approach, their layers can
only work on regular grids as underlying structure and not on general graphs.
A similar approach for convolutions using diffusion means was presented in
[61].

####  II-B2 GNNs for special manifolds

There have been two main motivations for building GNNs that can deal with
manifold-valued features: first, learning from interrelated measurements that
take values in some manifold, and second, embedding abstract graphs in curved
manifolds to utilize their geometric structure for downstream tasks.
Surprisingly, so far, most of the work concentrates on the second aim. The
resulting GNNs can only handle data from very restricted classes of manifolds.
We discuss the networks and their associated spaces in the following.

There is mounting empirical evidence that embedding graphs in non-Euclidean
spaces can help with various tasks. So the authors of [62] observed that a
hyperbolic representation is beneficial when dealing with graphs with a
hierarchical structure. Later, products of constant curvature spaces appeared
as appropriate embedding spaces for more classes of graphs [63]. In graph
representation learning, these observations were further backed up by certain
constructions that successfully embed (abstract) graphs into hyperbolic space
for downstream tasks [64, 65, 66, 67]. Other studies showed that other spaces
can be better suited than hyperbolic ones when learning, for example
classification tasks on wider classes of graphs: While spaces of constant
curvature (or products of them) were successfully used in [37, 38, 39, 40],
spaces of non-constant curvature were utilized in the form of Grassmannians
[68] and SPD [68, 42] manifolds. In [69], also pseudo-Riemmannian manifolds
were employed as embedding spaces in the form of pseudo-hyperboloids.

All the above approaches have in common that the initial embedding can be
chosen or learned freely. This allows for the choice of reference points and
the placement of nodes that capture the semantics of input features. Clearly,
we do not have this freedom if the data already comes with a fixed embedding.
Think, for example, of a network of cities that naturally live on the
2-sphere. The above works usually do not investigate the performance of their
networks on data whose embeddings are “far away” from the preferred ones. It
must thus be checked additionally if such applications are tackled.

An example where the graphs come with an initial embedding was investigated in
[22]. The authors considered SPD-valued graphs that naturally arise from EEG
imaging data, and developed network blocks that could navigate the space of
SPD matrices with the affine invariant metric. This led to very good
classification results.

Unfortunately, all the above GNN architectures have in common that they are
not equivariant or invariant under the isometries of the manifold. Whenever
the function to be learned (approximately) misses one of these properties,
this can lead to suboptimal performance; in such scenarios, non-invariant or
non-equivariant networks tend to have many “unnecessary” parameters that make
successful training more difficult. We refer to the discussion in [46]) on the
curse of dimensionality in machine learning and how shift-
invariant/equivariant CNNs helped to overcome it.

####  II-B3 Flow-based networks on manifolds

Neural ODEs were transferred to manifolds in [70, 71]. Besides that, there are
also deep learning approaches via discrete [72] and continuous [73, 74]
normalizing flows in Riemannian manifolds. For an overview of normalizing
flows and its generalizations for Euclidean data, we refer, for example, to
[75, 76]. In contrast to our approach, those authors aim to solve a given ODE
to learn intricate probability distributions instead of using the ODE to build
convolutional layers. The same is true for generative diffusion models [12,
77].

##  III Background

In this section, we recall the necessary background from differential geometry
(see, for example, [78, 79] for more), the notation of graph Laplacians for
manifold-valued graphs from [80], and the definitions of invariance and
equivariance.

###  III-A Riemannian Geometry

A Riemannian manifold is a $d$-dimensional manifold111Also when not mentioned
explicitly, we always assume that manifolds and all maps between manifolds are
smooth (that is, infinitely often differentiable). $M$ together with a
Riemannian metric $\langle\;\cdot\;,\;\cdot\;\rangle$. The latter assigns to
each tangent space $T_{p}M$ at $p\in M$ an inner product
$\langle\;\cdot\;,\;\cdot\;\rangle_{p}$ that varies smoothly in $p$. It
induces a norm $\|\;\cdot\;\|_{p}$ on each tangent space $T_{p}M$ and a
distance function $\textnormal{dist}:M\times M\to\mathbb{R}_{\geq 0}$.
Furthermore, it determines the so-called Levi-Civita connection $\nabla$,
which allows for the differentiation of vector fields. For vector fields $X$
and $Y$ on $M$, we denote the derivative of $Y$ along $X$ by $\nabla_{X}Y$.

A geodesic is a curve $\gamma:I\to M$ on an interval $I\subseteq\mathbb{R}$
without acceleration, that is,

$\nabla_{\gamma^{\prime}}\gamma^{\prime}=0,$

where $\gamma^{\prime}:=\frac{\textnormal{d}}{\textnormal{d}t}\gamma$. The
manifold $M$ is called complete if every geodesic can be defined on all of
$\mathbb{R}$. In the rest of the paper, we consider only complete, connected
Riemannian manifolds $M$. Locally, geodesics are shortest paths, and the
length of a geodesic connecting two points $p,q\in M$ is equal to the distance
$\textnormal{dist}(p,q)$. Importantly, every point in $M$ has a so-called
normal convex neighborhood $U\subseteq M$ in which any pair $p,q\in U$ can be
joined by a unique length-minimizing geodesic $\gamma:[0,1]\to M$ with
$\gamma(0)=p$ and $\gamma(1)=q$ that lies completely in $U$. We will need the
Riemannian exponential map: For $X\in T_{p}M$, let $\gamma_{X}$ be the
geodesic with $\gamma(0)=p$ and $\gamma^{\prime}(0)=X$. Since geodesics are
solutions to ordinary differential equations, there is a neighborhood
$W\subset T_{p}M$ of $0\in T_{p}M$ on which the so-called exponential map

$\displaystyle\exp_{p}$ $\displaystyle:W\to M,\quad X\mapsto\gamma_{X}(1)$

at $p$ is well defined. It can be shown that $\exp_{p}$ is a local
diffeomorphism. Let $D_{p}\subset T_{p}M$ be the maximal neighborhood of $0\in
D_{p}$, where this is the case, and set $D_{p}M:=\exp_{p}(D_{p})$. Then the
inverse $\log_{p}:D_{p}M\to D_{p}$ of $\exp_{p}$ is defined and is called the
Riemannian logarithm at $p$. The points $q$, for which $\log_{p}$ is not
defined, constitute the so-called cut locus of $p$. Products of Riemannian
manifolds are again Riemannian manifolds, where everything works component-
wise; this holds in particular for products $M^{n}$.

An important class of maps are the so-called isometries (or symmetries) of
$M$. These are diffeomorphisms $\Phi:M\to M$ that conserve the Riemannian
metric, that is, for every $X,Y\in T_{p}M$, $p\in M$, we have

$\langle\textnormal{d}_{p}\Phi(X),\textnormal{d}_{p}\Phi(Y)\rangle_{\Phi(p)}=\langle
X,Y\rangle_{p},$

where $\textnormal{d}_{p}\Phi:T_{p}M\to T_{\Phi(p)}M$ is the differential of
$\Phi$ at $p$. Isometries preserve the distance between any two points, map
geodesics into geodesics, and map normal convex neighborhoods to normal convex
neighborhoods. For every $p\in M$, the exponential and logarithm commute with
an isometry $\Phi$ as follows:

$\displaystyle\exp_{\Phi(p)}\big{(}\textnormal{d}_{p}\Phi(X)\big{)}$
$\displaystyle=\Phi\big{(}\exp_{p}(X)\big{)},\quad X\in T_{p}M,$ (1)
$\displaystyle\textnormal{d}_{p}\Phi\big{(}\log_{p}(q)\big{)}$
$\displaystyle=\log_{\Phi(p)}\big{(}\Phi(q)\big{)},\;q\in D_{p}M.$ (2)

The set of all isometries of $M$ constitutes the isometry group of $M$.

Another fundamental property of Riemannian geometry is curvature. The
geodesics define the real-valued sectional curvatures $K_{p}$ that intuitively
measure how much surfaces, which are images of planes under $\exp_{p}$, bend.
Important manifolds of constant sectional curvature $K:=K_{p}\equiv c$ are
spheres ($c>0$), Euclidean spaces ($c=0$), and hyperbolic spaces ($c<0$). The
curvature has a profound influence on the properties of the manifold. In
particular, if $M$ is a so-called Hadamard manifold [78, Ch. 6], that is, a
simply connected manifold with non-positive sectional curvature $K_{p}\leq 0$
everywhere, the whole manifold is a normal convex neighborhood. Apart from
Euclidean and hyperbolic spaces, the spaces of symmetric positive definite
matrices with the affine invariant metric belong to this class.

###  III-B Graph Laplacian in Manifolds

Let $G=(V,E,w)$ be a directed graph with vertex (node) set
$V=\\{v_{1},\dots,v_{n}\\}$, directed edge set $E\subset V\times V$, and
positive edge weights $w:E\to\mathbb{R}_{+}$. The ordering of the vertices
encodes the direction of the edge, where $(v,u)$ as an edge going from $v$ to
$u$, and we write $u\sim v$ to denote that there exists such an edge.

Starting with the set of functions

$\mathcal{F}(V,M):=\\{f:V\to M\\},$

we define the set of admissible $M$-valued vertex functions by

$\mathcal{H}(V,M):=\\{f\in\mathcal{F}(V,M)\ |f(u)\in
D_{f(v)}M\;\forall(v,u)\in E\\}.$

The definition ensures that $\log_{f(v)}f(u)$ exists whenever $(v,u)\in E$.
Consequently, if $\log_{p}(q)$ is defined for any $p,q\in M$, then
$\mathcal{H}(V,M)=\mathcal{F}(V,M)$. This holds, for instance, for Hadamard
manifolds. When a graph $G$ comes with a map $f\in\mathcal{H}(V,M)$, we call
$f$ (vertex) features of $G$. In this case, we also write $G=(V,E,w,f)$.

For $f\in\mathcal{H}(V,M)$, we denote the disjoint union of tangent spaces at
the values of $f$ by

$T_{f}M:=\bigcup_{v\in V}T_{f(v)}M,$

and the space of tangent space functions corresponding to $f$ by

$\mathcal{H}(V,T_{f}M):=\\{F:V\to T_{f}M\ |\ F(v)\in T_{f(v)}M\\}.$

Now the graph Laplacian $\Delta:\mathcal{H}(V,M)\to\mathcal{H}(V,T_{f}M)$ is
defined for $f\in\mathcal{H}(V,M)$ by

$\Delta f(v)=-\\!\\!\\!\sum_{u\sim v}w(v,u)\log_{f(v)}f(u)\in T_{f(v)}M,$ (3)

see [80]. If $M=\mathbb{R}^{d}$, this reduces to a well-known graph Laplacian
for Euclidean functions; see, for example, [81].

###  III-C Equivariance and Invariance

Equivariance and invariance of maps are essential properties of many neural
networks. Let $S$ be a set and $\mathcal{G}$ a group with neutral element by
$e$. Assume that $\mathcal{G}$ acts on $S$ (from the left), that is, there is
a (left) group action $\rho:\mathcal{G}\times S\to S$ such that $\rho(e,s)=s$
and $\rho(g,\rho(h,s))=\rho(gh,s)$. A function $\phi:S\to\mathbb{R}$ is called
invariant with respect to the group action of $\mathcal{G}$ if

$\phi(\rho(g,s))=\phi(s)$

for all $g\in\mathcal{G}$ and $s\in S$. A function $\Phi:S\to S$ is called
equivariant with respect to the group action of $\mathcal{G}$ if

$\Phi(\rho(g,s))=\rho(g,\Phi(s))$

for all $g\in\mathcal{G}$ and $s\in S$. Note that the concatenation of an
equivariant and invariant function is again an invariant one.

Equivariance of neural networks under some group operation is usually achieved
by stacking, that is, concatenating of layers that are all equivariant
themselves. If an invariant network is the goal, several equivariant layers
are followed by at least one invariant layer before the final output is
generated.

##  IV Diffusion Layer

Figure 1: Diffusion of a rectangle graph on the 2-sphere. Only the edges of
the initial graph are shown. Color indicates diffusion time.

In this section, we introduce a novel diffusion layer. Let $G=(V,E,w,f)$ be a
graph with feature map $f$. We augment $f$ with a time parameter $t$ and
consider for $a>0$, the $M$-valued graph diffusion equation

$\begin{cases}\frac{\partial}{\partial
t}\widetilde{f}(v,t)=-\Delta\widetilde{f}(v,t),&v\in V,\quad t\in(0,a),\\\
\widetilde{f}(v,0)=f(v),&v\in V.\end{cases}$ (4)

Figure 1 depicts the diffusion of a rectangle graph with constant edge
weights. In Appendix A, we show that there always exists a solution to
Equation (4) that lives for some time. We also have the following theorem.

######  Theorem IV.1.

Let $G=(V,E,w,f)$ be a graph with positive weights and features
$f\in\mathcal{H}(V,M)$ such that the smallest closed geodesic ball that
contains the features of the graph is convex. Assume $\sum_{u\sim v}w(v,u)\leq
1$ for all $v\in V$. Then, Equation (4) has a solution
$\widetilde{f}:V\times[0,\infty)\to M$ that is defined for all $t\geq 0$.

The proof is given in Appendix A.

As discussed in Section II, discrete approximation schemes were used to build
neural network layers from diffusion processes in Euclidean space. We will
transfer this to manifolds. To this end, we introduce the maps that are needed
for the explicit Euler discretization of Equation (4). Motivated by the
behavior of the Rectified Linear Unit (ReLU), we define for any $p\in M$,
nonlinear activation function $\sigma_{p}^{\alpha}:T_{p}M\to T_{p}M$ by

$\sigma_{p}^{\alpha}(X):=\begin{cases}X,&\|X\|_{p}\geq\alpha,\\\
0,&\text{otherwise.}\end{cases}$

Then we define for $f\in\mathcal{H}(V,M)$ and $t\geq 0$, the $1$-step map
$\textnormal{step}_{t,\alpha}:\mathcal{H}(V,M)\to\mathcal{F}(V,M)$ by

$(\textnormal{step}_{t,\alpha}f)(v):=\exp_{f(v)}\Big{(}\\!\\!-t\sigma_{f(v)}^{\alpha}\left(\Delta
f(v)\right)\Big{)},\quad v\in V.$

The map $\textnormal{step}_{t,\alpha}$ realizes a “nonlinearly activated”
explicit diffusion step of length $t$ of Equation (4). The activation thereby
allows only steps that are not shorter than $t\alpha$. Note that since the
graph Laplacian is only influenced by a node’s direct neighbors—its so-called
1-hop neighborhood—, the 1-step map only aggregates information over these
neighborhoods.

######  Remark IV.2.

As mentioned before, $\mathcal{H}(V,M)=\mathcal{F}(V,M)$ holds for example in
Hadamard manifolds. In this case, $\textnormal{step}_{t}^{\alpha}(f)$ cannot
lie outside the space of admissible functions. For spaces with positive
sectional curvatures,
$\textnormal{step}_{t}^{\alpha}(f)\notin\mathcal{H}(V,M)$ is in general
possible, when two vertices that share an edge lie in their respective cut
loci. However, any maximal neighborhood $D_{p}$, on which the exponential is
invertible, is dense in $M$ (see, for example, [82, Cor. 28.2]). Thus, every
cut locus has null measure under the volume measure of the manifold [83, Lem.
III.4.4]. Consequently, also in positively curved manifolds, the image
$\textnormal{step}_{t,\alpha}(f)$ lies in $\mathcal{H}(V,M)$ except for some
exceptional cases.

In Appendix B, we additionally show that the 1-step map always yields
admissible features even in positively curved spaces, whenever the data is
“local enough” and the step size is not too large. $\Box$

Since numerical approximation schemes usually employ several approximation
steps, we extend the 1-step map for $\ell>1$ and $f\in\mathcal{H}(V,M)$ with
$\textnormal{step}_{t,\alpha}^{l-1}(f)\in\mathcal{H}(V,M)$ as

$\displaystyle\textnormal{step}^{\ell}_{t,\alpha}(f):=\underbrace{\textnormal{step}_{t,\alpha}\circ\dots\circ\textnormal{step}_{t,\alpha}}_{\ell\text{
times}}(f).$ (5)

Again, $\textnormal{step}^{\ell}_{t,\alpha}$ realizes $\ell$ “nonlinearly
activated” explicit Euler steps of length $t>0$ for Equation (4). It directly
follows that the $\ell$-step map aggregates information over $\ell$-hop
neighborhoods, that is, a nodes new feature is influenced by (the features of)
all nodes that are connected to it through a sequence of not more than $\ell$
edges.

Now we can define our diffusion layer. Let $\Omega$ be the set of graphs with
vertex features.

######  Definition IV.3 (Diffusion Layer).

For $c\in\mathbb{N}$, let
$\boldsymbol{t}=(t_{1},\dots,t_{c})\in\mathbb{R}_{\geq 0}^{c}$ and
$\boldsymbol{\alpha}=(\alpha_{1},\dots,\alpha_{c})\in\mathbb{R}_{\geq 0}^{c}$.
A diffusion layer with $c\geq 1$ channels is a map
$\textnormal{Diff}_{\boldsymbol{t},\boldsymbol{\alpha}}:\Omega^{c}\to\Omega^{c}$
with

$\displaystyle\textnormal{Diff}_{\boldsymbol{t},\boldsymbol{\alpha}}\big{(}G_{1},\dots,G_{c}\big{)}=\Big{(}\big{(}V_{i},E_{i},w_{i},\textnormal{step}^{\ell}_{t_{i},\alpha_{i}}(f_{i})\big{)}\Big{)}_{i=1}^{c}.$

The diffusion layer takes $c$ graphs as input and diffuses them for certain
amounts of time. Thus, it consists of $c$ diffusion “channels.” The diffusion
times $\boldsymbol{t}$ and the step activations $\boldsymbol{\alpha}$ are the
learnable parameters. An important property of the new layer is that it can be
applied to graphs with varying number of nodes and connectivity patterns. It
is also clear that the information aggregation takes place in $\ell$-hop
neighborhoods. This local aggregation of information is thought to be one of
the keys for the success of convolutional (graph) neural networks in many
situations [46, 84].

Apart from the already discussed qualities, the diffusion layer possesses the
following equivariance properties.

######  Theorem IV.4.

The diffusion layer is equivariant under node permutations and isometric
transformations of the feature manifold.

###### Proof.

We can restrict ourselves to one channel. The equivariance of the diffusion
layer under permutations of the vertex set is already inherent in our
formulation as $f$ and Diff are defined on the unordered set $V$. Therefore,
if an order is imposed on the vertex set, the node features are sorted
accordingly, as well as the output features of the diffusion layer.

Next, we deal with the equivariance under isometries. Let $\Phi:M\to M$ be an
isometry. We must show that transforming features with $\Phi$ commutes with
applying the diffusion layer:

$\textnormal{step}^{\ell}_{t,\alpha}(\Phi\circ
f)=\Phi\circ\textnormal{step}^{\ell}_{t,\alpha}(f).$

Obviously, it suffices to show the relation for $\ell=1$. By definition and
Equation (2), we have for all $v\in V$ that

$\displaystyle\textnormal{step}_{t,\alpha}(\Phi\circ f)(v)$
$\displaystyle=\exp_{\Phi(f(v))}\bigg{(}\\!\\!\\!-t\sigma^{\alpha}_{\Phi(f(v))}\Big{(}\sum_{u\sim
v}w(v,u)\log_{\Phi(f(v))}\Phi(f(u))\Big{)}\bigg{)}$
$\displaystyle=\exp_{\Phi(f(v))}\bigg{(}\\!\\!\\!-t\sigma^{\alpha}_{\Phi(f(v))}\Big{(}\sum_{u\sim
v}w(v,u)(\textnormal{d}_{f(v)}\Phi)(\log_{f(v)}f(u))\Big{)}\bigg{)},$

and by linearity and norm preservation of $\textnormal{d}_{f(v)}\Phi$ further

$\displaystyle\textnormal{step}_{t,\alpha}(\Phi\circ f)(v)$
$\displaystyle=\exp_{\Phi(f(v))}\bigg{(}\textnormal{d}_{f(v)}\Phi\Big{(}-t\sigma_{f(v)}^{\alpha}\Big{(}\sum_{u\sim
v}w(v,u)\log_{f(v)}f(u)\Big{)}\bigg{)}.$

Finally, we use (1) to get

$\displaystyle\textnormal{step}_{t,\alpha}(\Phi\circ f)(v)$
$\displaystyle=\Phi\bigg{(}\exp_{f(v)}\Big{(}\\!\\!-t\sigma_{f(v)}^{\alpha}\Big{(}\sum_{u\sim
v}w(v,u)\log_{f(v)}f(u)\Big{)}\Big{)}\bigg{)}$
$\displaystyle=\Phi\circ\textnormal{step}_{t,\alpha}(f)(v).$

This proves the claim. ∎

Our diffusion layer is equivariant with respect to the symmetries of the graph
and the feature space. As already discussed earlier, restricting to layers and
networks that exhibit this behavior by default is a highly successful ansatz
that helps to counter the curse of dimensionality [46]. Thus, the diffusion
layer is a highly versatile building block for GNNs.

##  V Tangent Multilayer Perceptron

In this section, we introduce a “1x1 convolution” layer that can transform
node features and allows to alter the width, that is, the number of channels,
of the network. To this end, we transfer ideas from the vector neuron
framework [85] to the manifold setting to construct a linear layer on tangent
vectors with a subsequent nonlinearity. Due to the latter, the layer can be
stacked to increase the depth. Indeed, it works like the well-known multilayer
perceptron (MLP).

The layer transforms $c_{\textnormal{in}}$ features
$f_{1},\dots,f_{c_{\textnormal{in}}}\in\mathcal{H}(V,M)$ from
$c_{\textnormal{in}}$ channels into $c_{\textnormal{out}}$ output feature
channels $g_{1},\dots,g_{c_{\textnormal{out}}}\in\mathcal{H}(V,M)$. To this
end, the node features are mapped to a reference tangent space and two sets of
linear combinations thereof are learned. One set of vectors is needed to
define positive and negative half spaces in the tangent space. The
other—output—set is then (nonlinearly) transformed according to the half space
they lie in. Finally, the exponential map is applied to transfer the result
back to the manifold.

######  Definition V.1.

(Tangent MLP) For $f_{1},\dots,f_{c_{\textnormal{in}}}\in\mathcal{H}(V,M)$ and
each $v\in V$, let $\bar{f}(v)\in M$ be a point depending on
$\boldsymbol{f}_{v}:=\left(f_{i}(v)\right)_{i=1}^{c_{\textnormal{in}}}$ that
is equivariant under node permutations and isometries in $M$. Further assume
that $f_{1}(v),\dots,f_{c_{\textnormal{in}}}(v)\in D_{\bar{f}(v)}M$. The
tangent layer perceptron with weights $\omega^{i}_{j},\xi^{i}_{j}$,
$i=1,\dots,c_{\textnormal{in}}$, $j=1,\dots,c_{\textnormal{out}}$ and
nonlinear scalar function $\sigma:\mathbb{R}\to\mathbb{R}$ transforms the
features of each node into new features
$g_{1},\dots,g_{c_{\textnormal{in}}}\in\mathcal{H}(V,M)$ in three steps:
First, the directions

$\displaystyle X_{j}(\boldsymbol{f}_{v})$
$\displaystyle:=\sum_{i=1}^{c_{\textnormal{in}}}\omega_{j}^{i}\log_{\bar{f}(v)}\big{(}f_{i}(v)\big{)},$
$\displaystyle\widetilde{Y}_{j}(\boldsymbol{f}_{v})$
$\displaystyle:=\sum_{i=1}^{c_{\textnormal{in}}}\xi_{j}^{i}\log_{\bar{f}(v)}\big{(}f_{i}(v)\big{)},$
$\displaystyle Y_{j}(\boldsymbol{f}_{v})$
$\displaystyle:=\frac{\widetilde{Y}_{j}(\boldsymbol{f}_{v})}{\|\widetilde{Y}_{j}(\boldsymbol{f}_{v})\|_{\bar{f}(v)}}$

are computed, and the $X_{j}$ are orthogonally decomposed with respect to the
$Y_{j}$ as

$\displaystyle
X_{j}^{\ThisStyle{\abovebaseline[0pt]{\scalebox{-1.0}{$\SavedStyle\perp$}}}}(\boldsymbol{f}_{v}$
$\displaystyle:=\big{\langle}X_{j}(\boldsymbol{f}_{v}),Y_{j}(\boldsymbol{f}_{v})\big{\rangle}_{\bar{f}(v)}Y_{j}(\boldsymbol{f}_{v}),$
$\displaystyle X_{j}^{\perp}(\boldsymbol{f}_{v})$
$\displaystyle:=X_{j}(\boldsymbol{f}_{v})-X_{j}^{\ThisStyle{\abovebaseline[0pt]{\scalebox{-1.0}{$\SavedStyle\perp$}}}}(\boldsymbol{f}_{v}).$

Second, the nonlinear scalar function $\sigma$ is applied

$\displaystyle Z_{j}(\boldsymbol{f}_{v})$
$\displaystyle:=\frac{\sigma\big{(}\|X_{j}^{\ThisStyle{\abovebaseline[0pt]{\scalebox{-1.0}{$\SavedStyle\perp$}}}}(\boldsymbol{f}_{v})\|_{\bar{f}(v)}\big{)}}{\|X_{j}^{\ThisStyle{\abovebaseline[0pt]{\scalebox{-1.0}{$\SavedStyle\perp$}}}}(\boldsymbol{f}_{v})\|_{\bar{f}(v)}}X_{j}^{\ThisStyle{\abovebaseline[0pt]{\scalebox{-1.0}{$\SavedStyle\perp$}}}}(\boldsymbol{f}_{v})+X_{j}^{\perp}(\boldsymbol{f}_{v}),$

and third, the result is transformed back to $M$ via

$\displaystyle g_{j}(v;\boldsymbol{f}_{v})$
$\displaystyle:=\exp_{\bar{f}(v)}\big{(}Z_{j}(\boldsymbol{f}_{v})\big{)}.$

The tangent multilayer perceptron (tMLP) with $m$ layers is the concatenation
of $m$ tangent linear layers with nonlinearity. We give some further remarks
on the layer.

######  Remark V.2.

  * i)

For the reference point $\bar{f}$, different choices are possible. For
example, we can select one of the channel features, that is,
$\bar{f}(v):=f_{i}(v)$ for some $i\in\\{1,\dots,c_{\textnormal{in}}\\}$, or
the Fréchet mean [86] of some of the
$\left(f_{i}(v)\right)_{i=1}^{c_{\textnormal{in}}}$. The latter is more
expensive to compute and harder to differentiate, which is why we use the
former choice. Note that the tangent linear layer is not equivariant with
respect to permutations of the channels. Different from the ordering of the
nodes, the ordering of the channels is learned by the model, and we want to
give the network the freedom to extract information from it.

  * ii)

Classical options for the nonlinear function $\sigma$ are the rectified linear
unit (ReLU) or leaky ReLU.

  * iii)

For a tMLP with $m>1$, it is convenient to use the same reference points in
all layers. With this choice, the computation of
$g_{j}(v;\boldsymbol{f}_{v}))$ and the applications of the logarithm in the
next layer (for $X_{j}(\boldsymbol{f}_{v}))$ and $Y_{j}(\boldsymbol{f}_{v}))$)
cancel each other, so these operations can be left away for a higher
computational efficiency.

  * iv)

The arguments from Remark IV.2 concerning the cut locus of a non-Hadamard
manifold also apply to the tMLP.

  * v)

Instead of the tMLP, we could also use sequences of 1D-convolutions from [60].
However, when stacked, the latter require nonlinear activations to guarantee
non-collapsibility, and there are no isometry-equivariant manifold-valued
activation functions that we are aware of.

The following theorem shows that the tMLP shares the equivariance properties
of the diffusion layer.

######  Theorem V.3.

The tMLP is equivariant under node permutations and isometric transformations
of the feature manifold.

###### Proof.

We show the claim for a single tangent layer perceptron. This directly yields
the result for the tMLP since the latter is a concatenation of several linear
layers.

The equivariance under permutations of the nodes holds because the weights
$\omega^{i}_{j},\xi^{i}_{j}$ have the same value for every node, and
$\overline{f}$ is permuted just like the underlying vertices.

Let $\Phi:M\to M$ be an isometry. By (2), linearity of
$\textnormal{d}_{\bar{f}(v)}\Phi$, and since the reference points $\bar{f}$
transforms with the isometry, we obtain with
$\Phi\circ\boldsymbol{f}_{v}:=\big{(}\Phi\left(f_{i}(v)\right)\big{)}_{i=1}^{c_{\textnormal{in}}}$
that

$\displaystyle X_{j}(\Phi\circ\boldsymbol{f}_{v})$
$\displaystyle=\sum_{i=1}^{c_{\textnormal{in}}}\omega_{j}^{i}\log_{\Phi(\bar{f}(v))}\Big{(}\Phi\big{(}f_{i}(v)\big{)}\Big{)}$
$\displaystyle=\textnormal{d}_{\bar{f}(v)}\Phi\left(\sum_{i=1}^{c_{\textnormal{in}}}\omega_{j}^{i}\log_{\bar{f}(v)}\big{(}f_{i}(v)\big{)}\right)$
$\displaystyle=\textnormal{d}_{\bar{f}(v)}\Phi\big{(}X_{j}(\boldsymbol{f}_{v})\big{)}.$

Analogously, since $\textnormal{d}_{\bar{f}(v)}\Phi$ preserves norms, we
obtain

$\widetilde{Y}_{j}(\Phi\circ\boldsymbol{f}_{v})=\textnormal{d}_{\bar{f}(v)}\Phi\big{(}\widetilde{Y}_{j}(\boldsymbol{f}_{v})\big{)},$

and further

$Y_{j}(\Phi\circ\boldsymbol{f}_{v})=\frac{\textnormal{d}_{\bar{f}(v)}\Phi\big{(}\widetilde{Y}_{j}(\boldsymbol{f}_{v})\big{)}}{\|\textnormal{d}_{\bar{f}(v)}\Phi\big{(}\widetilde{Y}_{j}(\boldsymbol{f}_{v})\big{)}\|_{\Phi(\bar{f}(v))}}=\textnormal{d}_{\bar{f}(v)}\Phi\big{(}Y_{j}(\boldsymbol{f}_{v})\big{)}.$

Then, since isometries preserve angles as well, we get

$\displaystyle
X_{j}^{\ThisStyle{\abovebaseline[0pt]{\scalebox{-1.0}{$\SavedStyle\perp$}}}}(\Phi\circ\boldsymbol{f}_{v})$
$\displaystyle=\big{\langle}X_{j}(\Phi\circ\boldsymbol{f}_{v}),Y_{j}(\Phi\circ\boldsymbol{f}_{v})\big{\rangle}_{\Phi\circ\bar{f}(v)}Y_{j}(\Phi\circ\boldsymbol{f}_{v}),$
$\displaystyle=\textnormal{d}_{\bar{f}(v)}\Phi\big{(}X_{j}^{\ThisStyle{\abovebaseline[0pt]{\scalebox{-1.0}{$\SavedStyle\perp$}}}}(\boldsymbol{f}_{v})\big{)},$
$\displaystyle X_{j}^{\perp}(\Phi\circ\boldsymbol{f}_{v})$
$\displaystyle=\textnormal{d}_{\bar{f}(v)}\Phi\big{(}X_{j}(\boldsymbol{f}_{v})-X_{j}^{\ThisStyle{\abovebaseline[0pt]{\scalebox{-1.0}{$\SavedStyle\perp$}}}}(\boldsymbol{f}_{v})\big{)}$
$\displaystyle=\textnormal{d}_{\bar{f}(v)}\Phi\big{(}X_{j}^{\perp}(\boldsymbol{f}_{v})\big{)}.$

Thus, we conclude since norms are preserved,

$\displaystyle Z_{j}\big{(}\Phi\circ\boldsymbol{f}_{v}\big{)}$
$\displaystyle=\frac{\sigma\big{(}\|X_{j}^{\ThisStyle{\abovebaseline[0pt]{\scalebox{-1.0}{$\SavedStyle\perp$}}}}(\Phi\circ\boldsymbol{f}_{v})\|_{\Phi(\bar{f}(v))}\big{)}}{\|X_{j}^{\ThisStyle{\abovebaseline[0pt]{\scalebox{-1.0}{$\SavedStyle\perp$}}}}(\Phi\circ\boldsymbol{f}_{v})\|_{\Phi(\bar{f}(v))}}X_{j}^{\ThisStyle{\abovebaseline[0pt]{\scalebox{-1.0}{$\SavedStyle\perp$}}}}(\Phi\circ\boldsymbol{f}_{v})$
$\displaystyle\quad+X_{j}^{\perp}(\Phi\circ\boldsymbol{f}_{v})$
$\displaystyle=\textnormal{d}_{\bar{f}(v)}\Phi\big{(}Z_{j}(\boldsymbol{f}_{v})\big{)}.$

Finally, (1) yields

$\displaystyle g_{j}\big{(}v;\Phi\circ\boldsymbol{f}_{v}\big{)}$
$\displaystyle=\exp_{\Phi(\bar{f}(v))}\Big{(}Z_{j}(\Phi\circ\boldsymbol{f}_{v})\Big{)}$
$\displaystyle=\exp_{\Phi(\bar{f}(v))}\Big{(}\textnormal{d}_{\bar{f}(v)}\Phi\big{(}Z_{j}(\boldsymbol{f}_{v}\big{)}\Big{)}$
$\displaystyle=\Phi\Big{(}\exp_{\bar{f}(v)}\big{(}Z_{j}(\boldsymbol{f}_{v})\big{)}\Big{)}$
$\displaystyle=\Phi\big{(}g_{j}(v;\boldsymbol{f}_{v})\big{)},$

which finishes the proof. ∎

##  VI Manifold GCN

We now describe a generic graph convolutional neural network (GCN) block to
which only the last task-specific layers need to be added for a full
architecture: It consists of a sequence

$\textnormal{Diff}\rightarrow\textnormal{tMLP}\rightarrow\textnormal{Diff}\rightarrow\dots\rightarrow\textnormal{tMLP}\rightarrow\textnormal{Diff}\rightarrow\textnormal{tMLP},$

in which the numbers of channels of each layer can be chosen freely, as long
as the number of output channels of one layer equals the number of input
channels of the next. When a graph is fed into the network, $c$ copies of it
are given to the first diffusion layer, where $c$ is the number of channels it
has. While the diffusion layers aggregate information over the graph, the
tMLPs ensure that the architecture is not equivalent to a single diffusion
layer (that is, it does not collapse). The tMLPs increase the nonlinearity and
allow the network to transfer information between channels.

Later, when we speak of the “depth” of a manifold GCN block, the number of
diffusion layers is meant. When all layers have the same number of channels
$c$, we say that the block has “width” $c$.

When learning scalar outputs (as, for example, in classification), it is
usually beneficial to add layers after the manifold GCN block that make the
network invariant under node permutations and isometries. We found that the
invariant layer from [60, Sec. 2.3] works well; we apply it node-wise and add
graph-wise (max- and/or mean-)pooling afterward. While the invariant layer
makes the network invariant under isometries [60, Prop. 3], the global pooling
operation ensures invariance under node permutations. An MLP with a trailing
log-softmax function can finally be added to map, for example, to class
probabilities for classification.

We found that normalizing the graph weights can boost the performance of GNNs
that use a manifold GCN block. The reason, in all likelihood, is that if the
step size is too large, the $\ell$-step map will be far from the continuous
flow, so that the evolution of the graphs becomes “chaotic.” Before feeding a
graph $G=(V,E,w,f)$ with $b:=\max_{v\in V}\sum_{u\sim v}w(v,u)>1$ into the
network, we thus recommend a global scaling of the weights:

$\widetilde{w}:=\frac{w}{b}.$

Together with “short” initial diffusion times (usually distributed tightly
around some $t_{0}\in[0,1]$), we have observed the best performance. Appendix
B, which discusses how Diff behaves in certain scenarios, gives further
reasons why normalization is helpful.

##  VII Experiments

The results of our experiments are in this section. We used our manifold GCN
in two graph classification tasks.

###  VII-A Synthetic Graphs

A benchmark [65] for graph classification algorithms is to let it learn
whether a graph was created using the Erdős-Renyi [87], Barabasi-Albert [88],
or Watts-Strogatz [89] algorithm. Examples of graphs that were created with
each of them are depicted in Figure 2. The authors of [65, 66] showed that
classifying from a learned embedding in hyperbolic space is superior to
classifying from a Euclidean one. We used our manifold GCN for this task.
Since our network can work on arbitrary manifolds, we also compare the
hyperbolic space to the SPD space as an embedding space. The latter is thought
to be advantageous as it has a more complex geometric structure [42].

Figure 2: Graphs created with the Erdős-Renyi (left), Barabasi-Albert
(middle), and Watts-Strogatz (right) algorithms.

####  VII-A1 Data

We used differently-sized data sets of graphs built with the three creation
algorithms. Each set was balanced; that is, it contained the same number of
graphs for each algorithm. Every graph had 100 nodes. The number of edges was
chosen (as in [65]) in the following way: For Barabási-Albert graphs, the
number of edges to attach from a new node to existing nodes was discrete-
uniform distributed between 1 and 200. For Erdős-Rényi, the probability of
edge creation was chosen from a uniform distribution on $[0.1,1]$. For Watts-
Strogatz, each node was connected to a discrete-uniform distributed random
number between 1 and 200 neighbors in the ring topology, and the probability
of rewiring each edge was taken from a uniform distribution on $[0.1,1]$. All
edge-weights were set to 1/100.

####  VII-A2 Architecture and Loss

We used a manifold GCN block of depth 2. The initial flow layer had 5
channels; the following tMLP, which consisted of only one layer and used leaky
ReLU as nonlinearity, increased this number to 8. To map to scalar features,
we added the invariant layer from [60]. For each node, it calculated the
distances to two weighted means, outputting scalar features in 8 channels.
Graph-wise max- and mean-pooling were then used for aggregation, and a
standard 2-layered MLP with leaky ReLU mapped to the class probabilities using
a final log-softmax. The network was trained using the standard cross-entropy
loss.

####  VII-A3 Embedding Manifolds

In our experiments, we embedded the graphs in hyperbolic and SPD spaces. As
the model of hyperbolic space, we used the Lorentz model (see, e.g., [65] for
formulas), as it performed best on the given task in [65, Table 1]. When we
used the SPD space, we endowed the latter with the affine-invariant metric
[86].

####  VII-A4 Initial Encoding

For the initial embedding of a graph into hyperbolic space, we used the
following options:

  * (i)

The “degree” embedding from [65]: We assign to each node a one-hot vector $X$
of length 101, the index of the one-entry being the degree of the node. Since
the last entry of $X$ is 0, this vector is a tangent vector at the origin
$o=[0\ \cdots\ 0\ 1]^{T}$ of the 100-dimensional Lorentz-hyperboloid
$\mathcal{H}^{100}$. We can then use a (learned) linear function that maps to
the tangent space at the origin $o_{d}$ of a $d$-dimensional Lorentz space.
The final step is then to apply the exponential $\exp_{o_{d}}$ of the Lorentz
model to obtain a feature $\exp_{o_{d}}(X)\in\mathcal{H}^{d}$ for the node.

  * (ii)

The “one-hot” embedding: We choose an order for the graph’s nodes and assign
to each node the one-hot vector $X\in T_{o}\mathcal{H}^{100}$ of length 101
that indicates its place. The final feature is $\exp_{o}(X)$.

While the first embedding is equivariant under node permutations because of
the use of the degree, the one-hot embedding is not. Therefore, networks such
as HGNN [65] cannot use it without sacrificing invariance under permutations.
For our manifold GCN, this is different since two one-hot embeddings can be
transformed into each other by a series of rotations (about its symmetry axis)
and reflections (about hyperplanes defined by the coordinate axes). Since the
latter are isometries of the Lorentz model, the equivariance under isometries
of the proposed layers secures the permutation-invariance of the network, even
with the one-hot embedding.

For the SPD space, we used a “one-hot” embedding similar to the hyperbolic
space. There, the tangent space at the identity matrix consists of all
symmetric matrices of the same size. We embedded the graphs in the space of
$15$-by-$15$ SPD matrices, as this is the smallest such space whose number of
free off-diagonal entries is not smaller than the number of nodes. We assigned
a symmetric “one-hot” matrix to each node, that is, a symmetric matrix with
two one-entries off the diagonal and zero everywhere else. This leads again to
permutation-invariant networks for the following reason: We can always
transform a one-hot matrix into another via a congruence transformation with
the product of two elementary permutation matrices. This represents an
isometry since the affine-invariant metric is invariant under congruence with
orthogonal matrices [90].

TABLE I: Classification of synthetic graphs in Hyperbolic space. Method | Mean F1 Score | # parameters
---|---|---
# Graphs | 90 | 180 | 360 | 1080 | 2880 |
HGNN ($d=6$) | $0.518\pm 0.133$ | $0.567\pm 0.085$ | $0.605\pm 0.060$ | $0.654\pm 0.040$ | $0.681\pm 0.042$ | 7053
HGNN ($d=100$) | $0.418\pm 0.103$ | $0.413\pm 0.105$ | $0.404\pm 0.097$ | $0.547\pm 0.154$ | $0.766\pm 0.079$ | 161903
Ours (degree) | $0.594\pm 0.114$ | $0.638\pm 0.084$ | $0.645\pm 0.062$ | $0.679\pm 0.037$ | $0.756\pm 0.046$ | 10429
Ours (one-hot) | $\textbf{0.633}\pm 0.129$ | $\textbf{0.687}\pm 0.093$ | $\textbf{0.744}\pm 0.059$ | $\textbf{0.787}\pm 0.034$ | $\textbf{0.820}\pm 0.024$ | 429

####  VII-A5 Evaluation and Comparison Methods

As comparison method, we used the hyperbolic graph neural network (HGNN) from
[65].222It performed similarly well as the network from [66] on the task (in
particular, reported differences are small and have not been tested for
statistical significance). Furthermore, there is a tested code available
online. We tried two HGNN models, embedding the graphs in the Lorentz models
of 100-dimensional and 6-dimensional hyperbolic space, respectively. We also
used two of our networks: one using the degree and the other the one-hot
embedding, both in 100-dimensional space.

To compare how well the networks learn with varying amounts of training data,
we tested them on data sets with 90, 180, 360, 1800, and 2880 graphs. As
measure of classification accuracy, we used the (macro) F1 score. We always
split the data set into training, validation, and test sets using a 4:1:1
ratio. The training was done with ADAM. After each epoch, the accuracy on the
validation set was computed. Our final model was chosen as the last one with
the highest validation score, and its score on the test set was reported. For
all but the set with 2880 graphs, we repeated this process 100 times, each
time with a new (random) data set; because the standard deviation of the
results went down (and the computation times up), we only performed 25
repetitions for the former set.

To also check the performance for large data sets, we created a single data
set of 6000 graphs and trained our method with one-hot encoding and the HGNN
methods on it. To compare embedding manifolds, we also trained our network
using the SPD(15) space. As the evaluation method for this experiment, we
chose 3-fold cross-validation.

####  VII-A6 Software

All our experiments were performed using Python 3.11. For computations in the
hyperbolic and SPD spaces, we used Morphomatics 3.0.2 [91]. The graphs were
created with NetworkX 3.2.1 [92] and jraph 0.0.6. The network was implemented
in Haiku 0.0.11 using JAX 4.20. The training was done on a GPU using the ADAM
implementation of Optax 0.1.7. The parameters were updated using Optax’s
incremental update function with step size $0.1$. All experiments with HGNN
were performed using the code that is offered by the authors
online333https://github.com/facebookresearch/hgnn, and the software they
require.

####  VII-A7 Parameter Settings

We used a learning rate of $10^{-3}$ and trained with balanced batches of size
3 for 60 epochs based on our empirical observations. The computations were
performed with double precision.

####  VII-A8 Results and Discussion

The average F1 scores and the standard deviations, as well as the number of
each network’s trainable parameters, can be found in Table I. For each sample
size, manifold GCN with the one-hot embedding performs clearly better than its
competitors. This is impressive, especially since it is also the method with
the lowest number of trainable parameters. The large difference in the number
of parameters for our methods is due to the linear layer in the node
embedding. It seems clear that the stronger separation of the embedded nodes
helps the network in the learning process. The fact that it is restricted to
the degree embedding probably also affects the performance of the HGNN model.

We can also observe that our method performs better when there is less
training data. HGNN with a 100-dimensional embedding space starts only to
learn when 1080 graphs are in the training set. When the dimension is reduced
to 6, it has better performance with small data sets, but the results for
bigger sets are worse than with the higher dimension. We think that the better
performance for small data sets comes from the inductive bias that we
introduce through our isometry-invariant network.

To make sure that the results also hold for the large data sets used in [65],
we assembled a set with 6000 graphs and trained on one split. We also included
the SPD manifold as the base manifold there. We obtained the following
classification accuracies:444The difference between our results and those in
[65] is not surprising since, there, graphs with up to 500 nodes were used;
our results show that the task is harder when there are fewer nodes.

  * •

Ours (Hyperbolic one-hot): $0.847$,

  * •

Ours (SPD): $0.818$,

  * •

HGNN ($d=6$): $0.690$,

  * •

HGNN ($d=100$): $0.838$.

These results confirm our findings. They further indicate that, for this task,
the SPD space is not superior to the hyperbolic space.

###  VII-B Alzheimer’s disease

Figure 3: Hippocampus mesh with vertex normals (left) and the same normals,
colored like their footpoints, as points on the 2-sphere (right).

Alzheimer’s disease (AD) is a neurodegenerative disorder that is diagnosed
more and more often worldwide. Studies showed that AD is characterized by an
atrophy pattern of the brain, particularly of the hippocampi; see, for
example, [93]. We derived sphere-valued features from digital representations
of right hippocampi in the form of normals of triangle meshes. Note that the
normals are equivariant under rigid motions of the meshes and invariant under
scalings. Since atrophy affects hippocampal volume, we fed this information
into the manifold GCN network to classify whether the subject had AD or not.

####  VII-B1 Data

For our experiment, we prepared a data set consisting of 60 subjects diagnosed
with AD and 60 cognitive normal (CN) controls based on data from the open
access Alzheimer’s Disease Neuroimaging Initiative555http://adni.loni.usc.edu/
(ADNI) database. Among others, ADNI provides 1632 magnetic resonance images of
brains with segmented hippocampi. In our experiment, we employed the
hippocampus meshes derived in [36]. Each mesh was interpreted as a graph. As
edge weights, we used the edge weights from the mesh’s cotangent Laplacian.
Features were the vertex normals; the latter were computed as the averaged
normals of the faces to which the vertex belonged. Due to the equivariance of
the normals, our network does not need an alignment of the hippocampi.
Nevertheless, the imaging protocol assured a coarse alignment in the
coordinate system of the MRI scanner. A visualization of a hippocampus mesh
with its normals and their location on the sphere is shown in Figure 3.

####  VII-B2 Architecture and Loss

We used a manifold GCN of depth 4 with 16 channels in each layer. Each tMLP
consisted of one layer and used leaky ReLU activations. The last layers are
the same as in the experiment with synthetic graphs (an invariant layer, a
pooling layer, and a final MLP); the only differences are that each has 16
channels and the volume of the hippocampus is appended to the output of the
pooling layer. The network was also trained using the standard cross-entropy
loss.

####  VII-B3 Evaluation and Comparison Methods

We used two different comparison methods. Mesh CNN [94] is a state-of-the-art
deep neural network for learning from polygonal meshes. It takes the original
meshes as input and combines specialized convolution and pooling layers to
predict the class. As the second comparison method, we built a Euclidean GNN
that resembled our proposed network very closely. The manifold GCN block was
mimicked using the convolutions from [84] (instead of the diffusion layers)
and linear layers (instead of the depth-one tMLPs). We measured the accuracy
of the classification, that is, the ratio of correctly classified subjects. We
always split the data set into training, validation, and test sets using a
3:1:1 ratio and used ADAM for training. After each epoch, the accuracy on the
validation set was computed. For each method, the final model was chosen as
the first one with the highest validation score, and its score on the test set
was reported.

####  VII-B4 Software

We used the same software as in the experiment with the synthetic graphs for
our network and its training. Additionally, pyvista 0.42.3 was used to compute
the vertex normals and the volumes of the meshes. For MeshCNN, we used the
code that is available online.666https://ranahanocka.github.io/MeshCNN/ The
Euclidean GNN was built with the graph convolutions provided by jraph 0.0.6.

####  VII-B5 Parameter Settings

We used a learning rate of $10^{-3}$ and trained with batch size 1 for 300
epochs based on our empirical observations. The computations were performed
with double precision.

####  VII-B6 Results and Discussion

The results of our experiment and the number of trainable parameters for each
network are shown in Table II. Our network has the best accuracy among the
tested methods. From the performance of Mesh CNN, it becomes clear that the
model reduction we perform by considering volume and surface normals instead
of the full mesh helps tremendously in this application. Indeed, both methods
that utilize only the former information work considerably better. The fact
that the Euclidean GNN performs worse shows that the classification is
improved by considering the vertex normals as elements of their natural
domain, the sphere, and not as elements of the surrounding space.

Importantly, we did not observe antipodal points that were connected with an
edge, so no problems with cut loci appeared for our method.

TABLE II: Shape classification results. Method | Mean Accuracy | # parameters
---|---|---
GCN | $0.752\pm 0.085$ | 2616
Mesh CNN | $0.592\pm 0.073$ | 1320558
Ours (Sphere) | $\textbf{0.765}\pm 0.075$ | 2554

##  VIII Conclusion

We have presented two new GNN layers for manifold-valued features. They are
equivariant under the symmetries of the domain and the feature space of the
graph and can be combined to form a highly versatile GNN block that can be
used for many deep learning tasks. Contrary to existing GNN methods, our
layers can handle data from a wide range of manifolds, which opens up
possibilities for deep learning applications on manifold-valued data that
could not be tackled with GNNs before.

We applied networks with our new layers to two graph classification tasks,
observing very good performances. Tests with little training data suggest that
architectures that use it need less training data, probably because of the
equivariance properties. This might prove helpful in applications where
training data is hard to come by.

Problems that can arise with GNNs are over-squashing and bottleneck phenomena
[95]. They lead to situations in which it is very hard to transfer information
between distant nodes. A remedy that is presented in [53] is to build layers
from diffusion using implicit instead of explicit numerical schemes. This
helps to overcome, for example, bottlenecks, as global information is used for
an update. This idea can be transferred to our case: future research can focus
on making a diffusion layer based on an implicit (Euler) discretization of the
diffusion equation. It also seems interesting to test the use of
$p$-Laplacians [80] with $p\neq 2$.

## Acknowledgments

Martin Hanik is supported by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy—MATH+: The Berlin
Mathematics Research Center, EXC-2046/1—project ID: 390685689.

##  Appendix A Properties of the Graph Diffusion Equation

In this section, we discuss properties of the diffusion equation (4). In
particular, we prove Theorem IV.1.

For a graph $G=(V,E,w,f)$ be a graph with (ordered) vertices
$V=\\{v_{1},\dots,v_{n}\\}$, we define

$\displaystyle\Omega_{G}:=$
$\displaystyle\\{\boldsymbol{p}=(p_{1},\ldots,p_{n})\in M^{n}\ |\ \exists
f\in\mathcal{H}(V,M)$ $\displaystyle\text{such
that}\;f(v_{i})=p_{i},\;i=1,\ldots,n\\}.$

This set is open, since every maximal neighborhood $D_{p}\subseteq T_{p}M$ in
which $\exp_{p}$ is a diffeomorphism is open. Therefore, $\Omega_{G}$ is a
submanifold of $M^{n}$ (of the same dimension). If $M$ is a Hadamard manifold,
then $\Omega_{G}=M^{n}$.

Let $J_{i}\subseteq\\{1,\dots,n\\}$ denote the index set of the neighbors of
$v_{i}$. We define the following vector field on $\Omega_{G}$:

$\displaystyle\boldsymbol{p}\mapsto\bigg{(}\sum_{j\in
J_{i}}w(v_{i},v_{j})\log_{p_{i}}p_{j}\bigg{)}_{i=1}^{n}.$ (6)

This vector field is smooth in $\Omega_{G}$, since, away from the cut locus,
$\log_{p}(q)$ is smooth both as a function in $p$ and $q$; see, for example,
[78, Prop. 18].

######  Lemma A.1.

i) Equation (4) has a solution if and only if there exists an integral curve
$\gamma:[0,a)\to\Omega_{G}$ with vector field (6) starting in
$(f(v_{1}),\dots,f(v_{n}))$. In this case, we have

$\gamma_{i}(t)=\widetilde{f}(v_{i},t),\quad i=1,\dots,n,\quad t\in[0,a).$ (7)

ii) There exists $a\in\mathbb{R}_{+}$ such that Equation (4) has a unique
solution on $[0,a)$, which is moreover smooth.

###### Proof.

i) Equation (4) is an autonomous system of $n$ ordinary differential equations
on $M$. Indeed, for $\widetilde{f}(v_{i},t):=p_{i}$, $i=1,\dots,n$, we find

$-\Delta\widetilde{f}(v_{i},t)=\sum_{j\in
J_{i}}w(v_{i},v_{j})\log_{p_{i}}(p_{j}).$ (8)

Therefore, if $\widetilde{f}$ solves (4), then we can define the curve
$\gamma$ by (7). Now, Equation (8) implies that the derivative of $\gamma$ at
$\boldsymbol{p}$ is given by (6), so it is an integral curve of the above
vector field. The other direction works by going backwards.
ii) If $\gamma$ is smooth, then also $\widetilde{f}$. The result now follows
from part i) since integral curves of smooth vector fields that start at an
arbitrary point and live for a positive amount of time always exist and are
smooth [79, Prop. 9.2]. ∎

Clearly, if $M$ has a cut locus, then $\Omega_{G}$ is not (geodesically)
complete; when the integral curve hits the boundary, it stops. However, if the
data is sufficiently localized, this cannot happen, as we show in this
appendix.

In the following, a _bounding ball of a graph_ is the smallest closed geodesic
ball that contains the features of the graph. It exists, for example, if the
graph is contained in the maximal neighborhood $D_{p}M$ of some exponential
$\exp_{p}$. A closed geodesic ball is convex if any two of its points can be
connected by a unique geodesic that never leaves the ball. It is well known
that sufficiently small geodesic balls are normal convex neighborhoods [78].

Now we can prove Theorem IV.1, which says that the graph diffusion “lives
forever” if the bounding ball of the initial features is convex.

###### Proof.

(Proof of Theorem IV.1 )

We know from Lemma (A.1) that there is a solution $\widetilde{f}$ that lives
at least for some time $a>0$. We must show that $a=\infty$, and we do so by
contradiction.

Assume therefore that $a<\infty$. Then, the Escape Lemma [79, Lem. 9.19]
implies that the corresponding integral curve $\gamma$ leaves every compact
subset of $\Omega_{g}$ before it stops at time $a$. Let $\overline{B}$ be the
bounding ball of $G$. Clearly, $\overline{B}^{n}$ is a compact subset of
$\Omega_{G}$. We will show that $\gamma$ cannot leave $\overline{B}^{n}$,
leading to the desired contradiction.

To this end, let $v\in V$ be a node whose feature lies on the boundary
$\partial\overline{B}$ of $\overline{B}$. Then, because $\sum_{u\sim
v}w(v,u)\leq 1$ and because geodesics never leave $\overline{B}$, the vector
$-\Delta f(v)$ either points to the interior of $\overline{B}$ or along its
boundary. This is also true for any other $\widehat{f}\in\mathcal{H}(V,M)$ for
which at least one feature lies on $\partial\overline{B}$ and all others in
the interior of $\overline{B}$. Hence, at the boundary
$\partial\overline{B}^{n}$, the vectors (6) point either into the interior of
$\overline{B}^{n}$ or are tangent to $\partial\overline{B}^{n}$. But then, the
integral curve $\gamma$ cannot leave $\overline{B}^{n}$; see [79, Lem. 9.33].
∎

######  Example.

On a $d$-dimensional sphere, any graph that is contained in an open hemisphere
has a convex bounding ball.

As for the corresponding process in the Euclidean space, constant functions
are stable777A graph with features $f$ is stable under (4) if the solution is
independent of $t$, that is, $\widetilde{f}(v,t)=f(v)$ for all $t\in[0,a)$.
under the diffusion (4), and many initial configurations will converge towards
the same one. An example is shown in Figure 1. Indeed, we conjecture that if a
graph’s bounding ball is convex, then its features will always converge to a
constant point. However, in some manifolds, there exist non-constant stable
functions, as the following example shows.

Figure 4: Stable tetrahedron-graph.

######  Example.

In the $2$-sphere, we can inscribe a regular tetrahedron such that the right-
hand side of Equation (3) vanishes at each node, making it stable under
diffusion. (For example, when all weights are set to $1/3$.) The graph is
shown in Figure 4.

The above example represents a special case of a class of diffusion-stable
graphs. Remember that with positive weights $(w_{i})_{i=1}^{n}$ such that
$\sum_{i=1}^{n}w_{i}=1$, the weighted Fréchet mean of $n$ points
$(p_{i})_{i=1}^{n}\in M^{n}$ is defined by

$\textnormal{wFM}\Big{(}(w_{i})_{i=1}^{n},(p_{i})_{i=1}^{n}\Big{)}:=\operatorname*{arg\,min}_{p\in
M}\sum_{i=1}^{n}w_{i}\textnormal{dist}(p,p_{i})^{2}.$ (9)

Taking the gradient of the sum on the right-hand side of (9) yields, with
$\overline{p}:=\textnormal{wFM}\Big{(}(w_{i})_{i=1}^{n},(p_{i})_{i=1}^{n}\Big{)}$,
the optimality condition

$\sum_{i=1}^{n}-2w_{i}\log_{\overline{p}}(p_{i})=0.$ (10)

This implies the following partial characterization of stable configurations:

######  Proposition A.2.

If for each node $v\in V$ with neighbors $u_{1},\dots,v_{m_{v}}\in V$ we have

$f(v)=\textnormal{wFM}\Big{(}\big{(}w(v,v_{i})\big{)}_{i=1}^{m_{v}},\big{(}f(v_{i})\big{)}_{i=1}^{m_{v}}\Big{)},$

then the graph is stable under diffusion.

###### Proof.

The result follows as the right-hand side of the diffusion equation 4 vanishes
using (10). ∎

For positively curved spaces, the proposition yields non-constant, stable
solutions towards which a diffusion process can flow. In particular, as the
hippocampi meshes are topological spheres, the graphs used in this experiment
had global support on the 2-sphere. We thus think that the trained neural
network observed how the graphs diffused towards a globally supported stable
configuration similar to the tetrahedron graph from Figure 4. This makes
contacts with the cut locus even more unlikely, since features do not have to
“cross an equator” in order to diffuse towards a constant function.

##  Appendix B The $\ell$-step map and normal convex neighborhoods

In this appendix, we give some theoretical results on the behavior of the
$\ell$-step map (5) for localized data.

######  Proposition B.1.

Let $M$ be a Riemannian manifold, and let $G=(V,E,w,f)$ be a graph with
positive weights and features $f\in\mathcal{H}(V,M)$ such that the bounding
ball $\overline{B}$ is convex. Assume $\sum_{u\sim v}w(v,u)\leq 1$ for all
$v\in V$. Then, there is a maximal $a>0$ such that for all $t\leq a$

$\textnormal{step}^{\ell}_{t,\alpha}(f)(V)\subset\overline{B}.$

###### Proof.

Let $v\in V$ be a node whose feature $f(v)$ lies on the boundary of
$\overline{B}$. Then, because $\sum_{u\sim v}w(v,u)\leq 1$, the vector
$-\Delta f(v)$ either points to the interior of $\overline{B}$ or along its
boundary. Thus, there is a maximal $t_{v}>0$ such that $\exp_{f(v)}(-t\Delta
f(v))\in\overline{B}$ for all $t\leq t_{v}$. Obviously, such a number also
exists for vertices whose features are in the interior of $\overline{B}$.
Hence, $a=\max_{v\in V}t_{v}$. ∎

We directly get two corollaries.

######  Corollary B.2.

Let the assumptions of Proposition B.1 be fulfilled. Then there exists $a>0$
such that for all $t\leq a$

$\textnormal{step}^{\ell}_{t,\alpha}(f)(v)\in\overline{B}.$

In particular, every step of the $\ell$-step map is well-defined.

By Corollary B.2 a diffusion layer cannot run into trouble with a cut locus,
if the bounding ball of an input graph is convex and the time parameter $t$ is
smaller than some $a>0$. The next corollary tells us even more about the
behavior of the $\ell$-step map.

######  Corollary B.3.

Let the assumptions of Proposition B.1 be fulfilled. Then there exists $a>0$
such that for all $t<a$

$\displaystyle\max_{v,u\in
V}\textnormal{dist}\big{(}\textnormal{step}^{\ell}_{t,\alpha}(f)(v),$
$\displaystyle\textnormal{step}^{\ell}_{t,\alpha}(f)(u)\big{)}$
$\displaystyle<\max_{v,u\in V}\textnormal{dist}\big{(}f(v),f(u)\big{)}.$

In machine learning applications, we expect that it is preferable when the
graphs shrink in diameter (or, at least, do not expand) since most deep
learning layers show such a “contractive” behavior; see [60] for a more in-
depth discussion on this topic. The corollaries suggest that the diffusion
times of a diffusion layer should be initialized close to zero and that the
weights of a graph should be normalized, as proposed in Section VI.

## References

  * [1] W. Torng and R. B. Altman, “Graph convolutional neural networks for predicting drug-target interactions,” _J. Chem. Inf. Model._ , vol. 59, no. 10, pp. 4131–4149, 2019.
  * [2] X.-M. Zhang, L. Liang, L. Liu, and M.-J. Tang, “Graph neural networks and their current applications in bioinformatics,” _Front. Genet._ , vol. 12, 2021.
  * [3] H.-C. Yi, Z.-H. You, D.-S. Huang, and C. K. Kwoh, “Graph representation learning in bioinformatics: trends, methods and applications,” _Brief. Bioinformatics_ , vol. 23, no. 1, p. bbab340, 2022.
  * [4] A. Sanchez-Gonzalez, N. Heess, J. T. Springenberg, J. Merel, M. Riedmiller, R. Hadsell, and P. Battaglia, “Graph networks as learnable physics engines for inference and control,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 4470–4479.
  * [5] A. Sanchez-Gonzalez, J. Godwin, T. Pfaff, R. Ying, J. Leskovec, and P. Battaglia, “Learning to simulate complex physics with graph networks,” in _Proceedings of the 37th International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, vol. 119. PMLR, 2020, pp. 8459–8468.
  * [6] J. Shlomi, P. Battaglia, and J.-R. Vlimant, “Graph neural networks in particle physics,” _Mach. Learn.: Sci. Technol._ , vol. 2, no. 2, p. 021001, 2020\.
  * [7] W. Fan, Y. Ma, Q. Li, Y. He, E. Zhao, J. Tang, and D. Yin, “Graph neural networks for social recommendation,” in _The World Wide Web Conference_ , 2019, pp. 417–426.
  * [8] Y. Wu, D. Lian, Y. Xu, L. Wu, and E. Chen, “Graph convolutional networks with Markov random field reasoning for social spammer detection,” in _Proceedings of the AAAI conference on artificial intelligence_ , vol. 34, no. 01, 2020, pp. 1054–1061.
  * [9] S. Kumar, A. Mallik, A. Khetarpal, and B. Panda, “Influence maximization in social networks using graph embedding and graph neural network,” _Inf. Sci._ , vol. 607, pp. 1617–1636, 2022.
  * [10] L. Wu, Y. Chen, K. Shen, X. Guo, H. Gao, S. Li, J. Pei, and B. Long, “Graph neural networks for natural language processing: A survey,” _Found. Trends Mach. Learn._ , vol. 16, no. 2, pp. 119–328, 2023.
  * [11] N. I. Fisher, T. Lewis, and B. J. Embleton, _Statistical analysis of spherical data_. Cambridge, UK: Cambridge University Press, 1993.
  * [12] V. De Bortoli, E. Mathieu, M. Hutchinson, J. Thornton, Y. W. Teh, and A. Doucet, “Riemannian score-based generative modelling,” _Adv. Neural Inf. Process. Syst._ , vol. 35, pp. 2406–2422, 2022.
  * [13] E. Nava-Yazdani, F. Ambellan, M. Hanik, and C. von Tycowicz, “Sasaki metric for spline models of manifold-valued trajectories,” _Comput. Aided Geom. Des._ , p. 102220, 2023.
  * [14] R. Bergmann, R. H. Chan, R. Hielscher, J. Persch, and G. Steidl, “Restoration of manifold-valued images by half-quadratic minimization,” _Inverse Probl. Imaging_ , vol. 10, no. 2, pp. 281–304, 2016.
  * [15] M. Gräf, S. Neumayer, R. Hielscher, G. Steidl, M. Liesegang, and T. Beck, “An optical flow model in electron backscatter diffraction,” _SIAM J. Imaging Sci._ , vol. 15, no. 1, pp. 228–260, 2022.
  * [16] A. Cherian and S. Sra, “Positive definite matrices: data representation and applications to computer vision,” _Algorithmic Advances in Riemannian Geometry and Applications: For Machine Learning, Computer Vision, Statistics, and Optimization_ , pp. 93–114, 2016.
  * [17] P. T. Fletcher and S. Joshi, “Riemannian geometry for the statistical analysis of diffusion tensor data,” _Signal Process._ , vol. 87, no. 2, pp. 250–262, 2007.
  * [18] K. You and H.-J. Park, “Re-visiting Riemannian geometry of symmetric positive definite matrices for the analysis of functional connectivity,” _NeuroImage_ , vol. 225, p. 117464, 2021.
  * [19] X. Pennec, “Manifold-valued image processing with SPD matrices,” in _Riemannian geometric statistics in medical image analysis_. Elsevier, 2020, pp. 75–134.
  * [20] M. Hanik, M. A. Demirtaş, M. A. Gharsallaoui, and I. Rekik, “Predicting cognitive scores with graph neural networks through sample selection learning,” _Brain Imaging Behav._ , pp. 1–16, 2021.
  * [21] E. Wong, J. S. Anderson, B. A. Zielinski, and P. T. Fletcher, “Riemannian regression and classification models of brain networks applied to autism,” in _Connectomics in NeuroImaging: Second International Workshop, CNI 2018, Held in Conjunction with MICCAI 2018, Granada, Spain, September 20, 2018, Proceedings 2_. Springer, 2018, pp. 78–87.
  * [22] C. Ju and C. Guan, “Graph neural networks on SPD manifolds for motor imagery classification: A perspective from the time-frequency analysis,” _IEEE Trans. Neural Netw. Learn. Syst._ , pp. 1–15, 2023.
  * [23] F. Nerrise, Q. Zhao, K. L. Poston, K. M. Pohl, and E. Adeli, “An explainable geometric-weighted graph attention network for identifying functional networks associated with gait impairment,” in _International Conference on Medical Image Computing and Computer-Assisted Intervention_. Springer, 2023, pp. 723–733.
  * [24] L. Lin, V. Rao, and D. Dunson, “Bayesian nonparametric inference on the Stiefel manifold,” _Stat. Sin._ , pp. 535–553, 2017.
  * [25] R. Chakraborty and B. C. Vemuri, “Statistics on the Stiefel manifold: Theory and applications,” _Ann. Stat._ , vol. 47, no. 1, pp. 415–438, 2019.
  * [26] J. Chen, G. Han, H. Cai, D. Yang, P. J. Laurienti, M. Styner, and G. Wu, “Learning common harmonic waves on Stiefel manifold–a new mathematical approach for brain network analyses,” _IEEE Trans. Med. Imaging_ , vol. 40, no. 1, pp. 419–430, 2020.
  * [27] C. Mantoux, B. Couvy-Duchesne, F. Cacciamani, S. Epelbaum, S. Durrleman, and S. Allassonnière, “Understanding the variability in graph data sets through statistical modeling on the Stiefel manifold,” _Entropy_ , vol. 23, no. 4, p. 490, 2021.
  * [28] P. Turaga, A. Veeraraghavan, A. Srivastava, and R. Chellappa, “Statistical computations on Grassmann and Stiefel manifolds for image and video-based recognition,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 33, no. 11, pp. 2273–2286, 2011.
  * [29] Z. Huang, R. Wang, S. Shan, and X. Chen, “Projection metric learning on Grassmann manifold with application to video based face recognition,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2015, pp. 140–149.
  * [30] Z. Huang, J. Wu, and L. Van Gool, “Building deep networks on Grassmann manifolds,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 32, no. 1, 2018.
  * [31] F. Park, J. Bobrow, and S. Ploen, “A Lie group formulation of robot dynamics,” _Int. J. Rob. Res._ , vol. 14, no. 6, pp. 609–618, 1995.
  * [32] M. Hanik, H.-C. Hege, and C. von Tycowicz, “Bi-invariant dissimilarity measures for sample distributions in Lie groups,” _SIAM J. Math. Data Sci._ , vol. 4, no. 4, pp. 1223–1249, 2022.
  * [33] R. Vemulapalli, F. Arrate, and R. Chellappa, “Human action recognition by representing 3D skeletons as points in a Lie group,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2014, pp. 588–595.
  * [34] R. Vemulapalli and R. Chellapa, “Rolling rotations for recognizing human actions from 3D skeletal data,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2016, pp. 4471–4479.
  * [35] C. von Tycowicz, F. Ambellan, A. Mukhopadhyay, and S. Zachow, “An efficient Riemannian statistical shape model using differential coordinates,” _Med. Image Anal._ , vol. 43, pp. 1–9, 2018.
  * [36] F. Ambellan, S. Zachow, and C. von Tycowicz, “Rigid motion invariant statistical shape modeling based on discrete fundamental forms,” _Med. Image Anal._ , vol. 73, p. 102178, 2021.
  * [37] G. Bachmann, G. Bécigneul, and O. Ganea, “Constant curvature graph convolutional networks,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 486–496.
  * [38] S. Zhu, S. Pan, C. Zhou, J. Wu, Y. Cao, and B. Wang, “Graph geometry interaction learning,” _Adv. Neural Inf. Process. Syst._ , vol. 33, pp. 7548–7558, 2020.
  * [39] L. Sun, Z. Zhang, J. Ye, H. Peng, J. Zhang, S. Su, and S. Y. Philip, “A self-supervised mixed-curvature graph neural network,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 36, no. 4, 2022, pp. 4146–4155.
  * [40] C. Deng, F. Xu, J. Ding, L. Fu, W. Zhang, and X. Wang, “Fmgnn: Fused manifold graph neural network,” 2023.
  * [41] Z. Huang and L. Van Gool, “A Riemannian network for SPD matrix learning,” in _Proceedings of the AAAI conference on artificial intelligence_ , vol. 31, no. 1, 2017, pp. 2036–2042.
  * [42] W. Zhao, F. Lopez, M. J. Riestenberg, M. Strube, D. Taha, and S. Trettel, “Modeling graphs beyond hyperbolic: Graph neural networks in symmetric positive definite matrices,” in _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_. Springer, 2023, pp. 122–139.
  * [43] J. Atwood and D. Towsley, “Diffusion-convolutional neural networks,” _Adv. Neural Inf. Process. Syst._ , vol. 29, pp. 2001–2009, 2016.
  * [44] Y. Li, R. Yu, C. Shahabi, and Y. Liu, “Diffusion convolutional recurrent neural network: Data-driven traffic forecasting,” in _International Conference on Learning Representations_ , 2018.
  * [45] J. Gasteiger, S. Weißenberger, and S. Günnemann, “Diffusion improves graph learning,” _Adv. Neural Inf. Process. Syst._ , vol. 32, pp. 13 366–13 378, 2019.
  * [46] M. M. Bronstein, J. Bruna, T. Cohen, and P. Veličković, “Geometric deep learning: Grids, groups, graphs, geodesics, and gauges,” _arXiv preprint arXiv:2104.13478_ , 2021.
  * [47] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, and S. Y. Philip, “A comprehensive survey on graph neural networks,” _IEEE Trans. Neural Netw. Learn. Syst._ , vol. 32, no. 1, pp. 4–24, 2020.
  * [48] J. Zhao, Y. Dong, M. Ding, E. Kharlamov, and J. Tang, “Adaptive diffusion in graph neural networks,” _Adv. Neural Inf. Process. Syst._ , vol. 34, pp. 23 321–23 333, 2021.
  * [49] Q. Chen, Y. Wang, Y. Wang, J. Yang, and Z. Lin, “Optimization-induced graph implicit nonlinear diffusion,” in _International Conference on Machine Learning_. PMLR, 2022, pp. 3648–3661.
  * [50] R. Liao, Z. Zhao, R. Urtasun, and R. S. Zemel, “Lanczosnet: Multi-scale deep graph convolutional networks,” _arXiv preprint arXiv:1901.01484_ , 2019.
  * [51] B. Chamberlain, J. Rowbottom, M. I. Gorinova, M. Bronstein, S. Webb, and E. Rossi, “Grand: Graph neural diffusion,” in _International Conference on Machine Learning_. PMLR, 2021, pp. 1407–1418.
  * [52] M. Thorpe, T. M. Nguyen, H. Xia, T. Strohmer, A. Bertozzi, S. Osher, and B. Wang, “Grand++: Graph neural diffusion with a source term,” in _International Conference on Learning Representations_ , 2021.
  * [53] B. Chamberlain, J. Rowbottom, D. Eynard, F. Di Giovanni, X. Dong, and M. Bronstein, “Beltrami flow and neural diffusion on graphs,” _Adv. Neural Inf. Process. Syst._ , vol. 34, pp. 1594–1609, 2021.
  * [54] J. Hansen and T. Gebhart, “Sheaf neural networks,” _arXiv preprint arXiv:2012.06333_ , 2020.
  * [55] C. Bodnar, F. Di Giovanni, B. Chamberlain, P. Liò, and M. Bronstein, “Neural sheaf diffusion: A topological perspective on heterophily and oversmoothing in GNNs,” _Adv. Neural Inf. Process. Syst._ , vol. 35, pp. 18 527–18 541, 2022.
  * [56] N. Sharp, S. Attaiki, K. Crane, and M. Ovsjanikov, “Diffusionnet: Discretization agnostic learning on surfaces,” _ACM Trans. Graph._ , vol. 41, no. 3, pp. 1–16, 2022.
  * [57] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, “Neural ordinary differential equations,” _Adv. Neural Inf. Process. Syst._ , vol. 31, 2018\.
  * [58] M. Poli, S. Massaroli, J. Park, A. Yamashita, H. Asama, and J. Park, “Graph neural ordinary differential equations,” _arXiv preprint arXiv:1911.07532_ , 2019.
  * [59] L.-P. Xhonneux, M. Qu, and J. Tang, “Continuous graph neural networks,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 10 432–10 441.
  * [60] R. Chakraborty, J. Bouza, J. Manton, and B. C. Vemuri, “Manifoldnet: A deep neural network for manifold-valued data with applications,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 44, no. 2, pp. 799–810, 2020.
  * [61] S. Sommer and A. Bronstein, “Horizontal flows and manifold stochastics in geometric deep learning,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 44, no. 2, pp. 811–822, 2020.
  * [62] D. Krioukov, F. Papadopoulos, M. Kitsak, A. Vahdat, and M. Boguná, “Hyperbolic geometry of complex networks,” _Phys. Rev. E_ , vol. 82, no. 3, p. 036106, 2010.
  * [63] A. Gu, F. Sala, B. Gunel, and C. Ré, “Learning mixed-curvature representations in product spaces,” in _International Conference on Learning Representations_ , 2019.
  * [64] I. Chami, Z. Ying, C. Ré, and J. Leskovec, “Hyperbolic graph convolutional neural networks,” _Adv. Neural Inf. Process. Syst._ , vol. 32, 2019.
  * [65] Q. Liu, M. Nickel, and D. Kiela, “Hyperbolic graph neural networks,” _Adv. Neural Inf. Process. Syst._ , vol. 32, 2019.
  * [66] J. Dai, Y. Wu, Z. Gao, and Y. Jia, “A hyperbolic-to-hyperbolic graph convolutional network,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 154–163.
  * [67] Y. Zhang, X. Wang, C. Shi, N. Liu, and G. Song, “Lorentzian graph convolutional networks,” in _Proceedings of the Web Conference 2021_. New York: Association for Computing Machinery, 2021, pp. 1249–1261.
  * [68] C. Cruceru, G. Bécigneul, and O.-E. Ganea, “Computationally tractable Riemannian manifolds for graph embeddings,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 35, no. 8, 2021, pp. 7133–7141.
  * [69] B. Xiong, S. Zhu, N. Potyka, S. Pan, C. Zhou, and S. Staab, “Pseudo-Riemannian graph convolutional networks,” _Adv. Neural Inf. Process. Syst._ , vol. 35, pp. 3488–3501, 2022.
  * [70] A. Lou, D. Lim, I. Katsman, L. Huang, Q. Jiang, S. N. Lim, and C. M. De Sa, “Neural manifold ordinary differential equations,” _Adv. Neural Inf. Process. Syst._ , vol. 33, pp. 17 548–17 558, 2020.
  * [71] I. Katsman, E. M. Chen, S. Holalkere, A. Asch, A. Lou, S.-N. Lim, and C. D. Sa, “Riemannian residual neural networks,” in _Thirty-seventh Conference on Neural Information Processing Systems_ , 2023.
  * [72] D. J. Rezende, G. Papamakarios, S. Racaniere, M. Albergo, G. Kanwar, P. Shanahan, and K. Cranmer, “Normalizing flows on tori and spheres,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 8083–8092.
  * [73] E. Mathieu and M. Nickel, “Riemannian continuous normalizing flows,” _Adv. Neural Inf. Process. Syst._ , vol. 33, pp. 2503–2515, 2020.
  * [74] N. Rozen, A. Grover, M. Nickel, and Y. Lipman, “Moser flow: Divergence-based generative modeling on manifolds,” _Adv. Neural Inf. Process. Syst._ , vol. 34, pp. 17 669–17 680, 2021.
  * [75] P. Hagemann, J. Hertrich, and G. Steidl, “Generalized normalizing flows via Markov chains,” in _Non-local data interactions: foundations and applications_. Cambridge, UK: Cambridge University Press, 2023.
  * [76] L. Ruthotto and E. Haber, “An introduction to deep generative modeling,” _DMV Mitteilungen_ , vol. 44, no. 3, pp. 1–24, 2021.
  * [77] J. Thornton, M. Hutchinson, E. Mathieu, V. De Bortoli, Y. W. Teh, and A. Doucet, “Riemannian diffusion Schrödinger bridge,” _arXiv preprint arXiv:2207.03024_ , 2022.
  * [78] P. Petersen, _Riemannian Geometry_. Berlin: Springer, 2006.
  * [79] J. M. Lee, _Introduction to Smooth Manifolds_. Berlin: Springer, 2012.
  * [80] R. Bergmann and D. Tenbrinck, “A graph framework for manifold-valued data,” _SIAM J. Imaging Sci._ , vol. 11, no. 1, pp. 325–360, 2018.
  * [81] G. Gilboa and S. Osher, “Nonlocal operators with applications to image processing,” _SIAM Multiscale. Model. Simul._ , vol. 7(3), pp. 1005–1028, 2008.
  * [82] M. M. Postnikov, _Geometry VI: Riemannian Geometry_. Berlin: Springer, 2013.
  * [83] T. Sakai, _Riemannian geometry_. Providence, USA: American Mathematical Soc., 1996.
  * [84] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” in _International Conference on Learning Representations_ , 2017.
  * [85] C. Deng, O. Litany, Y. Duan, A. Poulenard, A. Tagliasacchi, and L. J. Guibas, “Vector neurons: A general framework for SO(3)-equivariant networks,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 12 200–12 209.
  * [86] X. Pennec, “Intrinsic statistics on Riemannian manifolds: basic tools for geometric measurements,” _J. Math. Imaging Vis._ , vol. 25, pp. 127–154, 2006.
  * [87] P. Erdős and A. Rényi, “On the evolution of random graphs,” _Publication of the Mathematical Institute of the Hungarian Academy of Sciences_ , vol. 5, no. 1, pp. 17–60, 1960.
  * [88] A.-L. Barabási and R. Albert, “Emergence of scaling in random networks,” _Science_ , vol. 286, no. 5439, pp. 509–512, 1999.
  * [89] D. J. Watts and S. H. Strogatz, “Collective dynamics of ‘small-world’networks,” _Nature_ , vol. 393, no. 6684, pp. 440–442, 1998\.
  * [90] Y. Thanwerdas and X. Pennec, “O(n)-invariant Riemannian metrics on SPD matrices,” _Linear Algebra Appl._ , vol. 661, pp. 163–201, 2023.
  * [91] F. Ambellan, M. Hanik, and C. von Tycowicz, “Morphomatics: Geometric morphometrics in non-Euclidean shape spaces,” 2021, https://morphomatics.github.io/.
  * [92] A. A. Hagberg, D. A. Schult, and P. J. Swart, “Exploring network structure, dynamics, and function using networkx,” in _Proceedings of the 7th Python in Science Conference_ , G. Varoquaux, T. Vaught, and J. Millman, Eds., Pasadena, CA USA, 2008, pp. 11–15.
  * [93] S. G. Mueller, N. Schuff, K. Yaffe, C. Madison, B. Miller, and M. W. Weiner, “Hippocampal atrophy patterns in mild cognitive impairment and Alzheimer’s disease,” _Hum. Brain Mapp._ , vol. 31, no. 9, pp. 1339–1347, 2010.
  * [94] R. Hanocka, A. Hertz, N. Fish, R. Giryes, S. Fleishman, and D. Cohen-Or, “MeshCNN: A network with an edge,” _ACM Trans. Graph._ , vol. 38, no. 4, 2019.
  * [95] U. Alon and E. Yahav, “On the bottleneck of graph neural networks and its practical implications,” in _International Conference on Learning Representations_ , 2021.