Bundle Neural Network for message diffusion on graphs

## Abstract

Abstract The dominant paradigm for learning on graph-structured data is message passing. Despite being a strong inductive bias, the local message passing mechanism suffers from pathological issues such as over-smoothing, over-squashing, and limited node-level expressivity.
To address these limitations we propose Bundle Neural Networks (BuNN), a new type of GNN that operates via message diffusion over flat vector bundles – structures analogous to connections on Riemannian manifolds that augment the graph by assigning to each node a vector space and an orthogonal map. A BuNN layer evolves the features according to a diffusion-type partial differential equation. When discretized, BuNNs are a special case of Sheaf Neural Networks (SNNs), a recently proposed MPNN capable of mitigating over-smoothing. The continuous nature of message diffusion enables BuNNs to operate on larger scales of the graph and, therefore, to mitigate over-squashing. Finally, we prove that BuNN can approximate any feature transformation over nodes on any (potentially infinite) family of graphs given injective positional encodings, resulting in universal node-level expressivity. We support our theory via synthetic experiments, and showcase the strong empirical performance of BuNNs over a range of real-world tasks, achieving state-of-the-art results on several standard benchmarks in transductive and inductive settings.

## 1 Introduction

Graph Neural Networks (GNNs) [, ] are widely adopted machine learning models designed to operate over graph structures. A majority of GNN models function locally, respecting the topology by having nodes exchange messages directly with neighboring nodes. Such models are referred to as *Message Passing Neural Networks (MPNNs) []*. While message passing has repeatedly proven to be a useful inductive bias in many tasks, the local mechanism was shown to suffer from certain pathological issues, namely over-smoothing [, , ], over-squashing [, , ], and limited expressivity [, ].

To address these issues, we propose Bundle Neural Networks (BuNNs), a new type of GNN that operates over *flat vector bundles* – a structure that augments the graph by assigning to each node a vector space and an orthogonal map. BuNNs *do not perform ‘explicit’ message passing* through multiple steps of information exchange between neighboring nodes, but instead operate via what we call message diffusion. Each layer involves a node update step, and a diffusion step evolving the features according to a vector diffusion PDE as in Singer and Wu []. The resulting architecture enjoys the desirable properties of Sheaf Neural Networks [], in that they can avoid over-smoothing, but are more scalable and achieve better performance on a range of benchmark datasets. We also show how the continuous nature of message diffusion together with the bundle structure allows to mitigate issues such as over-squashing. Additionally, we prove that under the right conditions, BuNNs are compact uniform approximators, a type of universality result for node-level function approximation.

In summary, our contributions are the following:

- •
We derive BuNNs from heat equations over flat vector bundles, and show that flat vector bundles are more amenable to computation than general vector bundles (Section 3).
- •
We study the diffusion process and show that it can mitigate over-smoothing and over-squashing (Section 4). We support these claims with synthetic experiments (Section 6.1).
- •
We prove that, under certain conditions, BuNNs are compact uniform universal approximators, and this is, to our knowledge, the first of such results (Section 5).
- •
We show that BuNNs perform well on a range of tasks, for instance achieving a new state-of-the-art result on the Peptides-func dataset (Section 6.2).

## 2 Background

Let $\mathsf{G}=(\mathsf{V},\mathsf{E})$ be an undirected graph on $n=|\mathsf{V}|$ nodes with edges $\mathsf{E}$. We represent the edges via an adjacency matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$ where the entry $\mathbf{A}_{uv}$ for $u,v\in\mathsf{V}$ is $1$ if the edge $(u,v)\in\mathsf{E}$ and 0 0 otherwise. Let $\mathbf{D}$ be the diagonal degree matrix with entry $\mathbf{D}_{vv}=d_{v}$ equal to the degree of $v$. The graph Laplacian is defined as $\mathbf{L}:=\mathbf{D}-\mathbf{A}$ and the random walk normalized graph Laplacian is defined as $\mathbf{\mathcal{L}}:=\mathbf{I}-\mathbf{D^{-1}A}$. We assume that at each node $v\in\mathsf{V}$ we are given a $c$-dimensional signal $\mathbf{x}_{v}$ and group such signals into a matrix $\mathbf{X}\in\mathbb{R}^{n\times c}$. A *feature transformation* on a graph $\mathsf{G}$ is a permutation-equivariant map $f_{\mathsf{G}}:\mathbb{R}^{n\times c}\rightarrow\mathbb{R}^{n\times c^{\prime}}$ that transforms the node signals.

GNNs and feature transformations. A $\operatorname{GNN}_{\boldsymbol{\Theta}}$ is a (continuous) map parameterized by $\boldsymbol{\Theta}$ that takes as input a graph alongside node signals $\mathsf{G}=(\mathsf{V},\mathsf{E},\mathbf{X})$ and outputs a transformed signal $(\mathsf{V},\mathsf{E},\mathbf{X^{\prime}})$. A GNN on a graph $\mathsf{G}$ is therefore a feature transformation $\operatorname{GNN}_{\Theta}:\mathbb{R}^{n\times c}\rightarrow\mathbb{R}^{n
\times c^{\prime}}$. Given a collection of graphs $\mathcal{G}$, with feature space $\mathbb{R}$, a feature transformation $F$ of $\mathcal{G}$ is an assignment of every graph $\mathsf{G}\in\mathcal{G}$ to a feature transformation $F_{G}:\mathbb{R}^{n_{G}\times c}\rightarrow\mathbb{R}^{n_{G}\times c^{\prime}}$. The set of (continuous) feature transformations over a collection of graphs in $\mathcal{G}$ is denoted $\mathcal{C}\left(\mathcal{G},\mathbb{R}^{c},\mathbb{R}^{c^{\prime}}\right)$.

Cellular sheaves. A *cellular sheaf* [] $\left(\mathcal{F},\mathsf{G}\right)$ over an undirected graph $\mathsf{G}=\left(\mathsf{V},\mathsf{E}\right)$ augments $\mathsf{G}$ by attaching to each node $v$ and edge $e$ a vector space space called stalks and denoted by $\mathcal{F}(v)$ and $\mathcal{F}(e)$, usually the stalks are copies of $\mathbb{R}^{d}$ for some $d$. Additionally, every incident node-edge pair $v\trianglelefteq e$ gets assigned a linear map between stalks called restriction maps and denoted $\mathcal{F}_{v\trianglelefteq e}:\mathcal{F}(v)\rightarrow\mathcal{F}(e)$. Given two nodes $v$ and $u$ connected by an edge $(v,\ u)$, we can transport a vector $\mathbf{x}_{v}\in\mathcal{F}(v)$ from $v$ to $u$ by first mapping it to the stalk at $e=(v,u)$ using $\mathcal{F}_{v\trianglelefteq e}$, and mapping it to $\mathcal{F}(u)$ using the transpose $\mathcal{F}_{u\trianglelefteq e}^{T}$.
As a generalization of the graph adjacency matrix, the sheaf adjacency matrix $\mathbf{A}_{\mathcal{F}}\in\mathbb{R}^{nd\times nd}$ is defined as a block matrix in which each $d\times d$ block $(\mathbf{A}_{\mathcal{F}})_{uv}$ is $\mathcal{F}_{u\trianglelefteq e}^{T}\mathcal{F}_{v\trianglelefteq e}$ if there is an edge between $u$ and $v$ and $\mathbf{0}_{d\times d}$ if there is not. We similarly define the block diagonal * degree matrix* $\mathbf{D}_{\mathcal{F}}\in\mathbb{R}^{nd\times nd}$ as $\left(\mathbf{D}_{\mathcal{F}}\right)_{vv}:=d_{v}\mathbf{I}_{d\times d}$. Finally, we define the sheaf Laplacian as $\mathbf{L}_{\mathcal{F}}=\mathbf{D_{\mathcal{F}}}-\mathbf{A}_{\mathcal{F}}$. These matrices act as bundle generalizations of their well-known standard graph counterparts and we recover such matrices when $\mathcal{F}(v)\cong\mathbb{R}$ and $\mathcal{F}_{v\trianglelefteq e}{v}=1$ for all $v\in\mathsf{V}$ and $e\in\mathsf{E}$.

Vector bundles. When restriction maps are orthogonal, we call the sheaf a vector bundle, a structure analogous to connections on Riemannian manifolds. For this reason, the sheaf Laplacian also takes the name connection Laplacian []. In this case the product $\mathcal{F}_{u\trianglelefteq e}^{T}\mathcal{F}_{v\trianglelefteq e}$ is also orthogonal and is denoted $\mathbf{O}_{uv}$ referring to the transformation a vector undergoes when moved across a manifold via parallel transport. In this case we denote the node-stalk at $v$ by $\mathcal{B}(v)$, the bundle-adjacency by $\mathbf{A}_{\mathcal{B}}$ and the bundle Laplacian $\mathbf{L}_{\mathcal{B}}$, and its random-walk normalized version $\mathbf{\mathcal{L}}_{\mathcal{B}}:=\mathbf{I}_{dn\times dn}-\mathbf{D}_{
\mathcal{B}}^{-1}\mathbf{A}_{\mathcal{B}}$.

Consider a $d$-dimensional vector field over the graph, i.e. a $d$-dimensional feature vector at each node denoted $\mathbf{X}\in\mathbb{R}^{nd}$ in which the signals are stacked column-wise. Similarly to the graph case, the operation $\mathbf{D}_{\mathcal{B}}^{-1}\mathbf{A}_{\mathcal{B}}\mathbf{X}$ is an averaging over the vector field, and $\mathbf{\mathcal{L}}_{\mathcal{B}}$ a measure of smoothness, since:

$$ $\left(\mathbf{D}_{\mathcal{B}}^{-1}\mathbf{A}_{\mathcal{B}}\mathbf{X}\right)_{ u}=\frac{1}{d_{u}}\sum_{u:\left(v,u\right)\in\mathsf{E}}\mathbf{O}_{uv}\mathbf {x}_{v}\ \in\mathbb{R}^{d},\ \text{and}\ \left(\mathbf{\mathcal{L}}_{\mathcal{ B}}\mathbf{X}\right)_{u}=\frac{1}{d_{u}}\sum_{u:\left(v,u\right)\in\mathsf{E}} \left(\mathbf{x}_{u}-\mathbf{O}_{uv}\mathbf{x}_{v}\right)\ \in\mathbb{R}^{d}.$ $$

## 3 Bundle Neural Networks

In this Section, we derive the Bundle Neural Network model from heat diffusion equations over flat vector bundles. We then discuss its relationship to GCN [] and NSD []. We provide algorithmic implementations and additional practical details in the Appendix (Section E).

Heat diffusion over bundles. The *bundle Dirichlet energy* $\mathcal{E}_{\mathcal{B}}(\mathbf{X})$ of a vector field $\mathbf{X}\in\mathbb{R}^{nd}$ is:

$$ $\mathcal{E}_{\mathcal{B}}(\mathbf{X}):=\mathbf{X}^{T}\mathbf{\mathcal{L}}_{ \mathcal{B}}\mathbf{X}=\frac{1}{2}\sum_{(v,u)\in\mathsf{E}}\frac{1}{d_{u}} \left\|\mathbf{x}_{u}-\mathbf{O}_{uv}\mathbf{x}_{v}\right\|_{2}^{2}.$ $$

The gradient of the Dirichlet energy $\nabla_{\mathbf{X}}\mathcal{E}_{\mathcal{B}}(\mathbf{X})$ is the random-walk Laplacian $\mathbf{\mathcal{L}}_{\mathcal{B}}$. We can write down a heat diffusion equation (or simply heat equation) over a vector bundle as a gradient flow, whose evolution equation with initial condition $\mathbf{X}(0)=\mathbf{X}$ is:

$$ $\partial_{t}\mathbf{X}(t)=-\mathbf{\mathcal{L}}_{\mathcal{B}}\mathbf{X}(t)\ , \ \mathbf{X}(0):=\mathbf{X}\ .$ $$

It is well known (see Hansen and Gebhart []) that the solution to this type of equation can be written using matrix exponentiation as:

$$ $\mathbf{X}(t)=\exp(-t\mathbf{\mathcal{L}}_{\mathcal{B}})\mathbf{X}(0)$ $$

We call the operator $\mathcal{H}_{\mathcal{B}}(t):=\exp(-t\mathbf{\mathcal{L}}_{\mathcal{B}})\in
\mathbb{R}^{nd\times nd}$ the *bundle heat kernel*, which is defined as

$$ $\mathcal{H}_{\mathcal{B}}(t):=\lim_{K\to\infty}\sum_{k=0}^{K}\frac{\left(-t \mathcal{L}_{\mathcal{B}}\right)^{k}}{k!}.$ $$

Computing the heat kernel is necessary to solve the heat equation. An exact solution can be computed using spectral methods. For small $t$, one can instead consider the truncated Taylor expansion centered at 0 0, which amounts to fixing a $K$ in the above equation. Bodnar et al. [] instead approximates the solution using the Euler discretization with unit time step of the heat equation.

However, all these methods pose a challenge when the bundle structure is learned as in Bodnar et al. [], since the heat kernel has to be recomputed after every gradient update of the bundle structure. This high computational overhead limits the usability of Sheaf Neural Networks in applications.

Figure: Figure 2: Comparison of different Laplacians, from left to right: standard Laplacian, general bundle Laplacian, flat bundle Laplacian. Flat bundle Laplacian is amenable since it can be factored into a synchronization step, a standard Laplacian, and a desynchronization step.
Refer to caption: x1.png

Flat vector bundles. To address the scalability issues of general sheaves and vector bundles, we consider the special case of flat vector bundles in which every node $u$ gets assigned an orthogonal map $\mathbf{O}_{u}$, and every connection factorizes as $\mathbf{O}_{vu}=\mathbf{O}_{v}^{T}\mathbf{O}_{u}$. Consequently, the bundle Laplacian factors:

$$ $\mathbf{\mathcal{L}}_{\mathcal{B}}=\mathbf{O}^{T}\left(\mathbf{\mathcal{L}} \otimes\mathbf{I}_{d}\right)\mathbf{O},$ $$

where $\mathbf{O}\in\mathbb{R}^{nd\times nd}$ is block diagonal with $v$-th block being $\mathbf{O}_{v}$. We call the matrices $\mathbf{O}$ and $\mathbf{O}^{T}$ synchronization and desynchronization, respectively. We compare different Laplacians in Figure 2. This factorization avoids the $\mathcal{O}\left(d^{3}|\mathsf{E}|\right)$ cost of computing the restriction map over each edge. Additionally, it allows to cast the bundle heat equation into a standard graph heat equation:

###### Lemma 3.1 .

For every node $v$, the solution at time $t$ of the heat equation on a connected bundle $\mathsf{G}=\left(\mathsf{V},\ \mathsf{E},\ \mathbf{O}\right)$ with input node features $\mathbf{X}$ satisfies:

$$ $\left(\mathcal{H}_{\mathcal{B}}(t)\mathbf{X}\right)_{v}=\sum_{u\in\mathsf{V}} \mathcal{H}(t,v,u)\mathbf{O}_{v}^{T}\mathbf{O}_{u}\mathbf{x}_{u},$ $$

where $\mathcal{H}(t)$ is the standard graph heat kernel, and $\mathcal{H}(t,\ v,\ u)\in\mathbb{R}$ its the entry at $(v,\ u)$.

This reduces the computation of the bundle heat kernel to that of the cheaper graph heat kernel, an operator that does not change depending on the bundle and can, therefore, also be pre-computed.

Our model. The BuNN layer occurs in four steps, as illustrated in Figure 1 in which steps 1, 2, 3, and 4 correspond to the Equations below of the same number. First, the bundle maps $\mathbf{O}_{v}$ are computed using a neural network $\boldsymbol{\phi}$ and the use of Householder reflections [] or direct parameterization when $d=2$. Second, an encoder step updates the initial node signals via a matrix $\mathbf{W}\in\mathbb{R}^{d\times d}$, and bias $\mathbf{b}\in\mathbb{R}^{d}$. Next, the features are diffused over the *learned* vector bundle. Finally, a non-linearity $\sigma$ is applied. We summarize the steps in the following equations:

$$ $\displaystyle\mathbf{O}^{(\ell)}_{v}$ $\displaystyle:=\boldsymbol{\phi}^{\left(\ell\right)}(\mathsf{G},\mathbf{X^{( \ell)}},v)\ \ \ \forall v\in\mathsf{V}$ (1) $\displaystyle\mathbf{h}^{(\ell)}_{v}$ $\displaystyle:={\mathbf{O}_{v}^{\left(\ell\right)}}^{T}\mathbf{W}^{(\ell)}{ \mathbf{O}_{v}^{\left(\ell\right)}}\mathbf{x}^{(\ell)}_{v}+\mathbf{b}^{(\ell)} \ \ \ \forall v\in\mathsf{V}$ (2) $\displaystyle\mathbf{Z}^{(\ell+1)}$ $\displaystyle:=\mathcal{H}_{\mathcal{B}}(t)\mathbf{H}^{(\ell)}$ (3) $\displaystyle\mathbf{X}^{(\ell+1)}$ $\displaystyle:=\sigma\left(\mathbf{Z}^{(\ell+1)}\right)$ (4) $$

The diffusion time $t$ in Equation 3 is a hyperparameter determining the scale at which messages are diffused. For the case of small $t$, we approximate the heat kernel via its truncated Taylor series of degree $K$, and for large $t$, we use the exact solution using spectral methods. For simplicity of exposition, the steps above describe an update given a single bundle (i.e., $c=d$), meaning that $\mathbf{x}_{v}^{(\ell)}\in\mathbb{R}^{d}$. In general, we allow for multiple bundles and vector field channels, which we detail in Appendix E.
We note that both Equations 2 and 3 are linear (or affine), and the non-linearities lie in 1 and 4. Equation 2 may be interpreted as a bundle-aware encoder, while Equation 3 is the *message diffusion* step guided by the heat kernel.

We note that without the bundle structure, Equation 3 would converge exponentially fast to constant node representations over the graph (e.g., Theorem $1$ in []), potentially leading to over-smoothing. Accordingly, the bundle is crucial in this formulation to prevent node features from converging to constant. We note that this is a limitation of existing diffusion-based GNNs [, ].

Link to Graph Convolutional Networks. It is possible to derive Graph Convolutional Networks (GCNs) [] as an approximation of BuNNs operating over a trivial bundle. Setting $t=1$ Equation 3 becomes $\mathbf{Z}^{(l+1)}=\exp(-\mathcal{L}_{\mathsf{G}})\mathbf{H}^{(l)}$. The approximation $\exp(-\mathcal{L}_{\mathsf{G}})\approx 1-\mathcal{L}_{\mathsf{G}}$ gives the update $\mathbf{Z}^{(l+1)}=(1-\mathcal{L}_{\mathsf{G}})\mathbf{H}^{(l)}=\mathbf{A}_{
\mathsf{G}}\mathbf{H}^{(l)}$ recovering the GCN update. Section 4.1 shows that BuNNs can avoid over-smoothing, highlighting an important difference between GCNs and BuNNs.

Comparison with Sheaf Neural Networks. Flat vector bundles are a special case of *cellular sheaves* [, ], meaning that our model has close connections to Sheaf Neural Networks (SNNs) [, , , , ], with the closest relative being the Neural Sheaf Diffusion (NSD) model from Bodnar et al. []. BuNNs distinguish themselves from NSD in several ways. First, NSD approximates the heat equation using a time-discretized solution to the heat equation, which results in a standard message passing algorithm. In contrast, the direct use of the heat kernel allows BuNNs to break away from the explicit message-passing paradigm. Secondly, the use of flat bundles increase scalability since the bundle maps are computed at the node level. Additionally, flat bundles guarantees *path independence*, a requirement for the theory in Bodnar et al. [] to hold, often not satisfied for general sheaf constructions such as ones used in NSD. Thirdly, we allow $\phi$ to be any GNN while NSD restricts $\phi$ to be an MLP. We found that incorporating the graph structure in $\phi$ improved the experimental results. The update in Equation 2 is also different to NSD, and is necessary to prove our main theoretical result, Theorem 5.3. We provide experimental comparisons to NSDs in the experimental section and show that BuNNs significantly outperform their sheaf counterparts. We provide a summarized comparison between GCNs, SNNs, and BuNNs in Table 1.

**Table 1: Comparison between GCN, SNN, and BuNN in terms of propagation and message type, as well as their capability of mitigating issues related to GNNs.**
|  | GCN | SNN | BuNN |
| --- | --- | --- | --- |
| Propagation | $\mathbf{A}_{\mathsf{G}}\mathbf{X}$ | $\mathbf{\mathcal{L}}_{\mathcal{F}}\mathbf{X}$ | $\mathcal{H}_{\mathcal{B}}(t)\mathbf{X}$ |
| Message type | Standard<br>(Local) | Standard<br>(Local) | Diffusive<br>(Global) |
| No under-reaching | ✗ | ✗ | ✓ |
| No over-squashing | ✗ | ✗ | ✓ |
| No over-smoothing | ✗ | ✓ | ✓ |

###### Lemma 3.1 .

## 4 Properties of Bundle Neural Networks

In this Section, we give a formal analysis of the BuNN model. In Section 4.1, we derive the fixed points of the bundle heat diffusion, which is the subspace of signals towards which solutions converges exponentially fast, and show that despite converging to fixed points, BuNNs can retain information at the node level and therefore avoid over-smoothing. Section 4.2 discusses how our model can capture long-range interactions and mitigate over-squashing.

### 4.1 Fixed points and over-smoothing.

A significant limitation of MPNNs is the problem of over-smoothing, the phenomenon of node features becoming increasingly indistinguishable with the increase of the MPNN’s depth. Over-smoothing has been repeatedly highlighted as one of the critical problems in GNNs, making training deep GNNs much more challenging. Such behavior is related to the diffusion used in MPNNs being similar to that of heat diffusion on graphs [], which converges to uninformative fixed points(^1^11Up to degree scaling if using the symmetric-normalized Laplacian, see Li et al. [].) resulting in a loss of information at the node level.
Bodnar et al. [] show that the richer sheaf (or bundle) structure, however, gives rise to more interesting limiting behavior. Indeed, by Lemma 3.1, since $\lim_{t\rightarrow\infty}\mathcal{H}(t,\ v,\ u)=\frac{d_{u}}{2|\mathsf{E}|}$, the limit over time of a solution is $\frac{1}{2|E|}\sum_{u\in\mathsf{V}}d_{u}\mathbf{O}_{v}^{T}\mathbf{O}_{u}
\mathbf{x}_{u}$.
To better understand the space of fixed points, notice that any signal $\mathbf{Y}\in\mathbb{R}^{n\times d}$ expressible as $\mathbf{y}_{v}=\frac{1}{2|E|}\sum_{u\in\mathsf{V}}d_{u}\mathbf{O}_{v}^{T}
\mathbf{O}_{u}\mathbf{x}_{u}$ for some $\mathbf{X}$ is a fixed point. For any two nodes $u,\ v$ we get:

$$ $\mathbf{O}_{v}\mathbf{y}_{v}=\frac{1}{2|E|}\sum_{w\in\mathsf{V}}d_{w}\mathbf{O }_{w}\mathbf{x}_{w}=\mathbf{O}_{u}\mathbf{y}_{u}.$ (5) $$

Consequently, the fixed points of vector diffusion have a global geometric dependency where all nodes relate to each other by some orthogonal transformation, e.g. the output in Figure 1.

Equation 5 provides insight into the smoothing properties of BuNNs. When the bundle is trivial, the equation reduces to $\mathbf{y}_{v}=\mathbf{y}_{u}$, with no information at the node level, i.e. over-smoothing. When it is non-trivial, the output signal can vary across the graph (e.g., $\mathbf{y}_{v}\neq\mathbf{y}_{u}$ for two nodes $u$ and $v$). We summarize this in Proposition 4.1.

###### Proposition 4.1 .

Let $\mathbf{Y}$ be the output of a BuNN layer with $t=\infty$, where $G$ is a connected graph, and the bundle maps are not all equal. Then, there exists $u,v\in\mathsf{V}$ connected such that $\mathbf{y}_{v}\neq\mathbf{y}_{u}$.

###### Proposition 4.1 .

### 4.2 Over-squashing and long range interactions.

While message-passing in MPNNs constitutes a strong inductive bias, this local mechanism is problematic when the task requires the MPNN to capture interactions between distant nodes. These issues have been attributed mainly to the *over-squashing* of information [].
Topping et al. [] and Di Giovanni et al. [] formalize over-squashing through a sensitivity analysis, giving upper bounds on the partial derivatives of the outputs of the MPNN with respect to the input. In particular, they show that under weak assumptions on the message function, the Jacobian of MPNNs satisfies

$$ $|\partial\left(\mathrm{MPNN}_{\Theta}(\mathbf{X})\right)_{u}/\partial\mathbf{x }_{v}|\leq c^{r}\left(\mathbf{A}^{r}\right)_{uv},$ (6) $$

for any nodes $u,\ v$, where $r$ is the depth of the network and $c$ is a constant. Given two nodes $u,v$ at a distance $r$, message-passing will require at least $r$ layers for the two nodes to communicate and overcome *under-reaching*. If $r$ is large, distant nodes can communicate, but Di Giovanni et al. [] show that over-squashing becomes dominated by vanishing gradients. Built on top of such sensitivity analysis, we compute the Jacobian for BuNN:

###### Lemma 4.2 .

Let BuNN be a linear layer defined by Equations 2 & 3 with hyperparameter $t$. Then, for any connected graph and nodes $u,\ v$, we have

$$ $\frac{\partial\left(\operatorname{BuNN}\left(\mathbf{X}\right)\right)_{u}}{ \partial\mathbf{x}_{v}}=\mathcal{H}(t,u,v)\mathbf{O}_{u}^{T}\mathbf{W}\mathbf{ O}_{v},$ $$

The form of the Jacobian in Lemma 4.2 differs significantly from the usual form in Equation 6. First, all nodes communicate in a single BuNN layer since $\mathcal{H}(t,u,v)>0$ for all $u,v$, and $t$, allowing for *direct pair-wise communication between nodes*, making a BuNN layer operate globally similarly to the attention mechanism in the Transformer model, and therefore overcome under-reaching. Secondly, taking $t$ to be large allows BuNNs to operate on a larger scale, allowing stronger communication between distant nodes and overcoming over-squashing without the vanishing gradient problem. Indeed, $\mathcal{H}(t,u,v)$ converges to $\frac{d_{u}}{2|\mathsf{E}|}$ exponentially fast with $t\rightarrow\infty$, showing that for larger scales the sensitivity between two nodes does not depend on how they are connected.

Further, Lemma 4.2 gives a finer picture of the capabilities of BuNNs. For example, to mitigate over-squashing a node may decide to ignore information received from certain nodes while keeping information received from others. This allows the model to reduce the receptive field of certain nodes. We formalize such a behavior in Corollary 4.3.

###### Corollary 4.3 .

Consider n nodes $u,\ v,\ \text{and}\ w_{i},\text{ for}\ i=1,\ldots n-2$, of a connected graph with $2$ dimensional bundle such that $\mathbf{x}_{v}\neq\mathbf{x}_{w_{i}}$ and $\mathbf{x}_{u}\neq\mathbf{x}_{w_{i}}\ \forall i$. Then in a BuNN layer, at a given channel, the node $v$ can learn to ignore the information from all $w_{i}s$ while keeping information from $u$.

###### Lemma 4.2 .

###### Corollary 4.3 .

## 5 Expressivity of the model

In this Section, we study the expressive power of BuNNs at the *node level* from a function approximation perspective. We note that such expressive power is strictly stronger than the more common graph-level expressivity, which usually relates GNNs to the Weisfeiler-Leman hierarchy. In node-level approximation, we desire our model to be able to approximate *arbitrary feature transformations over the nodes*. Further, existing work generally focuses on a fixed graph or graphs of bounded size [, , ]. In contrast, we focus on the more challenging case of potentially infinite families of graphs, known as *uniform approximation* []. Such analysis studies the class of feature transformation that GNN architectures can approximate over a potentially *infinite collection* of graphs $\mathcal{G}$. Given input and output dimensions $c$ and $c^{\prime}$, a family of GNNs will parameterize a subset of continuous feature transformation $\mathcal{C}(\mathcal{G},\mathbb{R}^{c},\mathbb{R}^{c^{\prime}})$ over $\mathcal{G}$. It is natural to ask which feature transformations a specific family of GNNs can or cannot parameterize. To this end, we introduce the notion of *compact uniform approximation* as an amendment to the definition of uniform approximation from Rosenbluth et al. [] and discuss our specific choice in the Appendix (Section  D).

###### Definition 5.1 .

Let $\mathcal{F}\subseteq\mathcal{C}(\mathcal{G},\mathbb{R}^{c},\mathbb{R}^{c^{
\prime}})$ be a set of feature transformations over a family of graphs $\mathcal{G}$, and let $H\in\mathcal{C}(\mathcal{G},\mathbb{R}^{c},\mathbb{R}^{c^{\prime}})$ a feature transformation over $\mathcal{G}$. We say that $\mathcal{F}$ compactly uniformly approximates $H$, if for all finite subsets $\mathcal{K}\subseteq\mathcal{G}$, for all compact $K\subset\mathbb{R}^{c}$, and for all $\epsilon>0$, there exists an $F\in\mathcal{F}$ such that for all $\mathsf{G}\in\mathcal{K}$ and $\mathbf{X}\in K^{n_{\mathsf{G}}}$, we have that $||F_{\mathsf{G}}(\mathbf{X})-H_{\mathsf{G}}(\mathbf{X})||_{\infty}\leq\epsilon.$

As we show in Proposition 5.2, bounded-depth MPNNs fail to have this property for any arbitrary $\mathcal{G}$, since fixing the depth of the MPNNs to $\ell$, we can take a single compactly featured graph $\mathsf{G}\in\mathcal{G}$ with diameter larger than $\ell$. As there is a node whose receptive field does not include all nodes in such a $\mathsf{G}$, the architecture cannot uniformly approximate every function on $\mathcal{G}$.(^2^22This phenomenon in MPNNs is often called ‘under-reaching’.)

###### Proposition 5.2 .

There exists a family $\mathcal{G}$ consisting of connected graphs such that bounded-depth MPNNs are not compact uniform approximators, even if enriched with unique positional encoding.

It was shown in Rosenbluth et al. [] that *even with injective Positional Encodings (PEs)*, there are functions that MPNNs with virtual-nodes cannot approximate uniformly while transformers can, and vice-versa. Instead, we prove that our model is compactly uniformly universal given injective PEs. To the best of our knowledge, this is the first positive uniform approximation result for a GNN, demonstrating the remarkable modeling capabilities of BuNNs.

###### Theorem 5.3 .

Let $\mathcal{G}$ be a set of connected graphs with injective positional encodings. Then $2$-layer deep BuNNs with encoder/decoder have compact uniform approximation over $\mathcal{G}$.

In particular, let $\epsilon>0$ and $H$ be a feature transformation on a finite subset $\mathcal{K}\subseteq\mathcal{G}$ and $K\subseteq\mathbb{R}^{d}$ a compact set, then there is a $2$-layer deep BuNN with width of order $\mathcal{O}\left(\sum_{\mathsf{G}\in\mathcal{K}}|\mathsf{V}_{\mathsf{G}}|\right)$ that approximates $H$ over $\bigsqcup_{\mathsf{G}\in\mathcal{K}}K^{n_{\mathsf{G}}}\subseteq\bigsqcup_{
\mathsf{G}\in\mathcal{G}}\mathbb{R}^{n_{\mathsf{G}}d}$. In other words, the required hidden dimension of BuNN is only linearly dependent on the number of nodes in the family of graphs.

###### Definition 5.1 .

###### Proposition 5.2 .

###### Theorem 5.3 .

## 6 Experiments

In this Section, we evaluate our model through a range of synthetic and real-world experiments. We first validate the theory from Section 4 with two synthetic tasks. We then evaluate BuNNs on popular real-world benchmarks. We use truncated Taylor approximations for small values of $t$ ($\leq 1$), while for larger $t$, we use the exact spectral solution. We provide supplementary information on the implementation and precise experimental details in the Appendix (Sections E and F respectively).

### 6.1 Synthetic experiments: over-squashing and over-smoothing

Tasks. In this experiment, we compare BuNNs to MPNNs on two new datasets. The tasks are node-regression tasks, in which nodes must average the input features of a subset of nodes as illustrated in Figure 3. The input graph contains two types of nodes, denoted $A$ and $B$, whose features are sampled from disjoint distributions. The target for nodes in $A$ is to output the average input feature over nodes in $B$ and vice-versa.

Figure: Figure 3: Synthetic over-squashing (left) and over-smoothing (right) tasks on barbell and clique graphs, respectively. In both cases, the blue nodes output the average over the red nodes and vice-versa.
Refer to caption: x2.png

As a first experiment, to evaluate the capacity to mitigate over-squashing, we consider the case where the underlying graph is a barbell graph consisting of two fully connected graphs - each being a type and bridged by a single edge. This bottleneck makes it hard to transfer information from one cluster to another. As a second experiment, we test the capacity to mitigate over-smoothing by considering the fully connected graph, in which all nodes are connected. The fully connected graph is a worst-case scenario for over-smoothing since after one step of message passing, the features are already fully averaged over the graph. This over-smoothing makes it challenging to predict different outputs over the different clusters.

Setup. As a first baseline, we consider a constant predictor always predicting 0 0, the expected mean over the whole graph, making it unaware of any cluster-specific information. As a second baseline, we consider the cluster-specific constant predictor predicting the expected mean over the opposite cluster, that is, $\frac{\pm\sqrt{3}}{2}$, depending on the cluster. Additionally, we consider GNN baselines to be a node-level MLP, GCN [], GraphSAGE [], GAT [], and NSD []. The width of MPNNs is fixed to the minimal depth to avoid under-reaching, namely $3$ for the Barbell and $1$ for the fully connected graph, and ensure the width is large ($>128$) considering the task. We compare these to a BuNN with an MLP learning the bundle maps of a comparable number of parameters and the same depth. We use Adam optimizer with $10^{-3}$ learning rate, batch size $1$, and train for $500$ epochs. We use $100$ samples for training and $100$ samples for testing.

Results. The results for $N=10$ are reported in Table 2. All MPNNs perform poorly on the over-squashing task on the barbell graph. All MPNNs perform comparably to baseline $2$, which is aware of nodes in the same clusters but not the other, showing their incapability to transfer information between clusters. This is explained by the fact that the resistance distance between nodes from different clusters is $1.0$ while it is $0.2$ for nodes in the same cluster, making communication between nodes from different clusters hard. Indeed, resistance distance tightly connects to over-squashing, and the two nodes’ communication strength is inversely proportional to this metric Di Giovanni et al. []. On the other hand, BuNN achieves almost perfect accuracy on this task, which supports the claim that BuNN mitigates over-squashing. We note, however, that to solve the task perfectly, we need $t\geq 10$, allowing BuNN to operate on a larger scale more adapted to the task. To solve this task, BuNN can assign the orthogonal maps to separate the two types of nodes, making each node listen only to the nodes of the other type, a behavior proved to be possible in Corollary 4.3.

**Table 2: Mean squared error (MSE) of different models for the over-squashing and over-smoothing experiments on the barbell and clique graphs, respectively.**
| Task | Barbell<br>(over-squashing) | Clique<br>(over-smoothing) |
| --- | --- | --- |
| Base. 1 | $30.97\pm 0.42$ | $30.94\pm 0.42$ |
| Base. 2 | $1.00\pm 0.07$ | $0.99\pm 0.08$ |
| MLP | $1.08\pm 0.07$ | $1.10\pm 0.08$ |
| GCN | $1.05\pm 0.08$ | $29.65\pm 0.34$ |
| SAGE | $0.90\pm 0.29$ | $0.86\pm 0.10$ |
| GAT | $1.07\pm 0.09$ | $20.97\pm 0.40$ |
| NSD | $1.09\pm 0.15$ | $0.08\pm 0.02$ |
| BuNN | $\mathbf{0.01\pm 0.07}$ | $\mathbf{0.03\pm 0.01}$ |

Similarly, the over-smoothing task on the clique graph is also challenging. Indeed, GCN and GAT perform exceptionally poorly, comparably to Baseline $1$, which only has access to global information and not cluster-specific information. To the best of our knowledge, these are the only models with formal proofs of over-smoothing [, ] respectively. GraphSAGE performs slightly better because it processes neighbors differently than it processes nodes themselves, but it is still close to node-level MLPs and Baseline 2. Moreover, NSD solves the task as expected due to their capability to mitigate over-smoothing. Lastly, BuNN also learns to solve this task due to its capacity to mitigate over-smoothing. Indeed, it can learn to ignore nodes from a given cluster, a behavior proved to be possible in Corollary 4.3.

### 6.2 Real-world tasks

We evaluate BuNNs on the Long Range Graph Benchmark [] and the heterophilic tasks from Platonov et al. [] and. In our experiments, we report results both for the BuNN model described in Section 3 and a slight variation inspired by MixHop GNNs [], which we call BuNN-Hop. We found this variation to perform exceptionally well on some real-world tasks. We provide the implementation details in the Appendix (Section E).

**Table 3: Results for the heterophilic tasks. Accuracy is reported for roman-empire and amazon-ratings, and ROC AUC is reported for minesweeper, tolokers, and questions. Best results are denoted by First and Second. Asterisk^∗ denotes that some runs ran out of memory on an NVIDIA A10 GPU (24 GB).**
|  | roman-empire | amazon-ratings | minesweeper | tolokers | questions |
| --- | --- | --- | --- | --- | --- |
| GCN | $73.69\pm 0.74$ | $48.70\pm 0.63$ | $89.75\pm 0.52$ | $83.64\pm 0.67$ | $76.09\pm 1.27$ |
| SAGE | $85.74\pm 0.67$ | $53.63\pm 0.39$ | $93.51\pm 0.57$ | $82.43\pm 0.44$ | $76.44\pm 0.62$ |
| GAT | $80.87\pm 0.30$ | $49.09\pm 0.63$ | $92.01\pm 0.68$ | $83.70\pm 0.47$ | $77.43\pm 1.20$ |
| GAT-sep | $88.75\pm 0.41$ | $52.70\pm 0.62$ | $93.91\pm 0.35$ | $83.78\pm 0.43$ | $76.79\pm 0.71$ |
| GT | $86.51\pm 0.73$ | $51.17\pm 0.66$ | $91.85\pm 0.76$ | $83.23\pm 0.64$ | $77.95\pm 0.68$ |
| GT-sep | $87.32\pm 0.39$ | $52.18\pm 0.80$ | $92.29\pm 0.47$ | $82.52\pm 0.92$ | $78.05\pm 0.93$ |
| NSD | $80.41\pm 0.72$ | $42.76\pm 0.54$ | $92.15\pm 0.84$ | $78.83\pm 0.76^{*}$ | $69.69\pm 1.46^{*}$ |
| BuNN | $91.75\pm 0.39$ | $53.74\pm 0.51$ | $98.99\pm 0.16$ | $84.78\pm 0.80$ | $78.75\pm 1.09$ |
| BuNN-Hop | $89.63\pm 0.59$ | $52.01\pm 0.45$ | $97.91\pm 0.30$ | $84.33\pm 0.77$ | $78.68\pm 1.04$ |

Heterophilic datasets. As we have shown in Section 4, BuNNs are provably capable of avoiding over-smoothing. It is, therefore, natural to test how BuNN performs on heterophilic graphs where over-smoothing is recognized as an important limitation (see, e.g., Yan et al. []). We follow their methodology to evaluate BuNN on the $5$ heterophilic tasks proposed in Platonov et al. []. We run the models with $10$ different seeds and report the mean and standard deviation of the test accuracy for roman-empire and amazon-ratings, and mean and standard deviation test ROC AUC for minesweeper, tolokers, and questions. We use the classical baselines from Platonov et al. [] and provide the hyper-parameters in the Appendix (Section F.2).

We report the results in Table 3. BuNN achieves the best score on all tasks, with an average relative improvement of $4.4\%$. Its score on minesweeper is particularly impressive,
which is significantly ahead of the rest and for which BuNN solves the task perfectly. We found that the optimal value of $t$ over our grid search varies across datasets, being $1$ for amazon-ratings and $100$ for roman-empire. BuNN-Hop also performs similarly but tends slightly inferior to BuNN on $3$ of the datasets and inferior on $2$. BuNN consistently outperforms the sheaf-based model NSD. Such results showcase the strong modeling capacity of BuNN in heterophilic settings.

Long Range Graph Benchmark. In Section 4.2, we show that BuNNs have desirable properties when it comes to over-squashing and modeling long-range interactions. To verify such claims empirically, we evaluate BuNN on tasks from the Long Range Graph Benchmark (LRGB) []. We consider the Peptides dataset consisting of $15\,535$ graphs which come with two associated graph-level tasks, Peptides-func and Peptides-struct.

**Table 4: Results for the Peptides-struct and Peptides-func tasks from the Long Range Graph Benchmark (results are $\times 100$ for clarity). The best result is bold.**
|  | Peptides-func | Struct |
| --- | --- | --- |
| Model | Test AP $\uparrow$ | Test MAE $\downarrow$ |
| GCN | $68.60\pm 0.50$ | $24.60\pm 0.07$ |
| GINE | $66.21\pm 0.67$ | $24.73\pm 0.17$ |
| GatedGCN | $67.65\pm 0.47$ | $24.77\pm 0.09$ |
| GPS | $65.34\pm 0.91$ | $25.09\pm 0.14$ |
| DReW | $71.50\pm 0.44$ | $25.36\pm 0.15$ |
| GAPH ViT | $69.42\pm 0.75$ | $\mathbf{24.49}\pm 0.16$ |
| BuNN | $71.03\pm 0.22$ | $24.98\pm 0.12$ |
| BuNN-Hop | $\mathbf{71.92}\pm 0.22$ | $24.63\pm 0.25$ |

The graph classification task in Peptides-func is to predict the function of the peptide from $10$ classes, while the regression task in Peptides-struct is inferring the 3D properties of the peptides. In both cases, we follow the standard experimental setup detailed by Dwivedi et al. [] alongside the updated suggestions from Tönshoff et al. []. The performance metric is Average Precision (AP) (higher is better) for Peptides-func and Mean Absolute Error (MAE) (lower is better) for Peptides-struct. We run each experiment on $4$ distinct seeds and report mean and standard deviation over these runs. Baseline models are taken from Tönshoff et al. [] and include MPNNs, transformer models, and the current SOTA models [, ].

We report the results in Table 4. BuNNs achieve, to the best of our knowledge, a new state-of-the-art result on Peptides-func. The results on Peptides-struct are also competitive, with BuNN-Hop remaining within a standard deviation from the best result. The strong performance provides further evidence of the long-range capabilities of BuNNs.

## 7 Conclusion

In this work, we proposed Bundle Neural Networks – a new type of GNN that operates via message diffusion on graphs. We gave a formal analysis of BuNNs showing that message diffusion can mitigate issues such as over-smoothing - since the heat equation over vector bundles admits a richer set of fixed points - and over-squashing - since BuNNs can operate at a larger scale than standard MPNNs. We also prove compact uniform approximation of BuNNs, characterizing their expressive power and establishing their superiority over MPNNs. Finally, we showed that BuNNs perform well on a range of synthetic and real-world benchmarks including both inductive and transductive tasks.

## Appendix A Limitation and Future work

A limitation of our framework is that while message diffusion allows to operate on different scales of the graph, the computation of the heat kernel for large $t$ requires spectral methods and is therefore computationally expensive. An exciting research direction consists of using existing computational methods to approximate the heat kernel efficiently for large values of $t$. A limitation of our experiments is that we consider inductive graph regression/classification and transductive node regression/classification tasks but no link prediction task. A limitation in our theory is that Theorem  ef thm:unifapp assumes injective positional encodings at the node level, which might only sometimes be available; future work could characterize the expressiveness when these are unavailable.

## Appendix B Societal Impact

This paper presents work whose goal is to advance the field of Machine Learning. Our work proposes a novel framework and model that can potentially be used for positive and negative societal impacts. However, due to the technical nature of our work, we are not aware of any particular potential positive or negative societal impact that should be highlighted here.

## Appendix C Proofs

In this Section, we provide proof of the theoretical results from the main text. Namely Lemma 3.1, Proposition 4.1, Lemma 4.2, Corollary 4.3, Proposition 5.2, and finally Theorem 5.3.

###### Lemma 3.1 .

For every node $v$, the solution at time $t$ of heat diffusion on a connected bundle $\mathsf{G}=\left(\mathsf{V},\ \mathsf{E},\ \mathbf{O}\right)$ with input node features $\mathbf{X}$ satisfies:

$$ $\left(\mathcal{H}_{\mathcal{B}}(t)\mathbf{X}\right)_{v}=\sum_{u\in\mathsf{V}} \mathcal{H}(t,v,u)\mathbf{O}_{v}^{T}\mathbf{O}_{u}\mathbf{x}_{u},$ (7) $$

where $\mathcal{H}(t)$ is the standard graph heat kernel, and $\mathcal{H}(t,\ v,\ u)\in\mathbb{R}$ its the entry at $(v,\ u)$.

###### Proof.

Since $\mathbf{\mathcal{L}}_{\mathcal{B}}=\mathbf{O}_{\mathcal{B}}^{T}\mathbf{
\mathcal{L}}\mathbf{O}_{\mathcal{B}}$ we get $\mathcal{H}_{\mathcal{B}}(t,u,v)=\mathbf{O}_{\mathcal{B}}^{T}\mathcal{H}(t,u,v
)\mathbf{O}_{\mathcal{B}}$.
∎

###### Lemma 3.1 .

###### Proof.

### C.1 Over-smoothing: proofs of Section 4.1 .

In this section, we prove the result of over-smoothing. This result follows straightforwardly from Lemma 3.1.

###### Proposition 4.1 .

Let $\mathbf{Y}$ be the output of a BuNN layer with $t=\infty$, where $G$ is a connected graph, and the bundle maps are not all equal. Then, $u,v\in\mathsf{V}$ is connected such that $\mathbf{y}_{v}\neq\mathbf{y}_{u}$ almost always.

###### Proof.

Consider a stable signal $\mathbf{X}\in\ker\mathcal{L}_{\mathcal{B}}$ and pick $u,v\in\mathsf{V}$ such that $\mathbf{O}_{u}\neq\mathbf{O}_{v}$. As $\mathbf{X}$ is stable, it must have 0 0 Dirichlet energy, so we must have that $\mathbf{O}_{u}\mathbf{y}_{u}=\mathbf{O}_{v}\mathbf{y}_{v}$, but as $\mathbf{O}_{u}\neq\mathbf{O}_{v}$ we have that $\mathbf{y}_{u}\neq\mathbf{y}_{v}$, which holds except in degenerate cases such as when the matrices $\mathbf{O}_{u}$ are reducible, or when the original signal is the zero vector.
∎

###### Proposition 4.1 .

###### Proof.

### C.2 Over-squashing: proofs of Section 4.2 .

In this Section we prove our results on over-squashing and long-range interactions. The Jacobian result follows straightforwardly from Lemma 3.1 and the definition of a BuNN layer, and the Corollary follows from the Jacobian result.

###### Lemma 4.2 .

Let BuNN be a linear layer defined by Equations 2 & 3. Then, for any connected graph and nodes $u,\ v$, we have

$$ $\frac{\partial\left(\operatorname{BuNN}\left(\mathbf{X}\right)\right)_{u}}{ \partial\mathbf{x}_{v}}=\mathcal{H}(t,u,v)\mathbf{O}_{u}^{T}\mathbf{W}\mathbf{ O}_{v},$ $$

and therefore

$$ $\lim_{t\rightarrow\infty}\frac{\partial\left(\operatorname{BuNN}\left(\mathbf{ X}\right)\right)_{u}}{\partial\mathbf{x}_{v}}=\frac{d_{v}}{2|\mathsf{E}|} \mathbf{O}_{u}^{T}\mathbf{W}\mathbf{O}_{v}.$ $$

###### Proof.

The result follows from the closed-form solution of the heat kernel from Lemma 3.1. We start by applying the bundle encoder from Equation 2 that updates each node representation as $\mathbf{h}_{v}=\mathbf{O}_{v}^{T}\mathbf{W}\mathbf{O}_{v}\mathbf{x}_{v}+
\mathbf{b}$. We therefore compute

$$ $\displaystyle\frac{\partial\left(\operatorname{BuNN}\left(\mathbf{X}\right) \right)_{u}}{\partial\mathbf{x}_{v}}$ $\displaystyle=\frac{\partial}{\partial\mathbf{x}_{v}}\left[\sum_{v\in\mathsf{V }}\mathcal{H}(t,u,v)\mathbf{O}_{u}^{T}\mathbf{O}_{v}\left(\mathbf{O}_{v}^{T} \mathbf{W}\mathbf{O}_{v}\mathbf{x}_{v}+\mathbf{b}\right)\right]$ (8) $\displaystyle=\mathcal{H}(t,u,v)\mathbf{O}_{u}^{T}\mathbf{W}\mathbf{O}_{v}.$ (9) $$

The second statement follows from the fact that $\mathcal{H}(t,u,v)\rightarrow\frac{d_{u}}{2|E|}$
∎

To illustrate the flexibility of such a result, we examine a setting in which we want nodes to select which nodes they receive information from, therefore ‘reducing’ their receptive field.

###### Corollary 4.3 .

Consider n nodes $u,\ v,\ \text{and}\ w_{i},\ \text{for}\ i=1,\ldots n-2$ of a connected graph $\mathsf{G}$ with $2$ dimensional features such that $\mathbf{x}_{v}\neq\mathbf{x}_{w_{i}}$ and $\mathbf{x}_{u}\neq\mathbf{x}_{w_{i}}\ \forall i$. Assume the bundle structure is learned with a node-level MLP. Then, in a given channel, node $v$ can learn to ignore the information from all $w_{i}s$ while keeping information from $u$.

###### Proof.

We denote $\mathbf{y}$ the output of the layer, and index the two dimensions by super-scripts, i.e. $\mathbf{y}=\begin{pmatrix}\mathbf{y}^{(1)}\\
\mathbf{y}^{(2)}\end{pmatrix}$. Our goal is to have $\frac{\partial\mathbf{y}_{v}^{(1)}}{\partial\mathbf{y}_{u}}\neq\begin{pmatrix}
0&0\end{pmatrix}$, while $\frac{\partial\mathbf{y}_{v}^{(1)}}{\partial\mathbf{y}_{w_{i}}}=\begin{pmatrix
}0&0\end{pmatrix}$ for all $i$. This would make the first channel of the output at $v$ insensitive to the input at all $w_{i}$s while being sensitive to the input at node $u$.

Fix $\mathbf{O}_{v}=\mathbf{O}_{u}=\begin{pmatrix}0&1\\
1&0\end{pmatrix}$ and $\mathbf{O}_{w_{i}}=\begin{pmatrix}1&0\\
0&1\end{pmatrix}$. Such maps can always be learned by an MLP, by the assumptions on $\mathbf{x}_{v}$, $\mathbf{x}_{u}$, and $\mathbf{x}_{w_{i}}$ and by the universality of MLPs. Let the weight matrix be $\mathbf{W}=\begin{pmatrix}w_{11}&w_{12}\\
w_{21}&w_{22}\end{pmatrix}$. By Lemma 4.2 we get $\frac{\partial\mathbf{y}_{v}}{\partial\mathbf{x}_{u}}=\mathcal{H}(t,v,u)
\mathbf{O}_{v}^{T}\mathbf{W}\mathbf{O}_{u}=\mathcal{H}(t,v,u)\begin{pmatrix}w_
{22}&w_{12}\\
w_{21}&w_{11}\end{pmatrix}$ and $\frac{\partial\mathbf{y}_{v}}{\partial\mathbf{x}_{w_{i}}}=\mathcal{H}(t,v,w_{i
})\mathbf{O}_{v}^{T}\mathbf{W}\mathbf{O}_{w_{i}}=\mathcal{H}(t,v,w_{i})\begin{
pmatrix}w_{21}&w_{22}\\
w_{11}&w_{12}\end{pmatrix}$. Setting $w_{21}$ and $w_{22}$ to 0 0 gives $\frac{\partial\mathbf{y}^{(1)}_{v}}{\partial\mathbf{x}_{w_{i}}}=\begin{pmatrix
}0&0\end{pmatrix}$ and $\frac{\partial\mathbf{y}^{(1)}_{v}}{\partial\mathbf{x}_{u}}=\mathcal{H}(t,v,u)
\begin{pmatrix}0&w_{12}\end{pmatrix}\neq\mathbf{0}$, as desired.
∎

###### Lemma 4.2 .

###### Proof.

###### Corollary 4.3 .

###### Proof.

### C.3 Expressivity of BuNNs: proofs of Section 5 .

We now turn to BuNN’s expressivity. Before proving that BuNNs have compact uniform approximation, we prove that MPNNs fail to have this property. This proves BuNNs’ superiority and shows that uniform expressivity is a good theoretical framework for comparing GNN architectures.

###### Proposition 5.2 .

There exists a family $\mathcal{G}$ consisting of connected graphs such that bounded-depth MPNNs are not compact uniform approximators, even if enriched with unique positional encoding.

###### Proof.

Let $\mathcal{G}$ be any family of connected graphs with an unbounded diameter (for example, the $n\times n$ grids with $n\rightarrow\infty$). Let the depth of the MPNN be $L$. Let $\mathsf{G}\in\mathcal{G}$ be a graph with diameter $>L$, and let $u$ and $v$ be two nodes in $\mathsf{V}_{\mathsf{G}}$ at distance $>L$. Note that the output at $v$ will be insensitive to the input at $u$, and therefore, the MPNN cannot capture feature transformations where the output at $v$ depends on the input at $u$. This argument holds even when nodes are given unique positional encodings.
∎

We now turn to our main theoretical contribution. The proof of Theorem 5.3 is split into two parts. The first proves that $1$-layer BuNNs have compact uniform approximation over *linear feature transformations*. The second part is extending to continuous feature transformation, which is an application of classical results.

We start by recalling what a linear feature transformation over a family of graphs $\mathcal{G}$ is:

###### Definition C.1 .

A linear feature transformation $L\in\mathcal{C}(\mathcal{G},\mathbb{R}^{c},\mathbb{R}^{c^{\prime}})$ over a family of graphs $\mathcal{G}$ is an assignment of each graph $\mathsf{G}\in\mathcal{G}$ to a linear map $L_{\mathsf{G}}:\mathbb{R}^{n_{\mathsf{G}}c}\rightarrow\mathbb{R}^{n_{\mathsf{G
}}c^{\prime}}$. Here, linearity means that for any two node-signals $\mathbf{X}_{1}\in\mathbb{R}^{nc}$ and $\mathbf{X}_{2}\in\mathbb{R}^{nc^{\prime}}$, and any real number $\alpha\in\mathbb{R}$, it holds that $L_{\mathsf{G}}\left(\alpha\mathbf{X}_{1}\right)=\alpha L_{\mathsf{G}}\left(
\mathbf{X}_{1}\right)$, and $L_{\mathsf{G}}\left(\mathbf{X}_{1}+\mathbf{X}_{2}\right)=L_{\mathsf{G}}\left(
\mathbf{X}_{1}\right)+L_{\mathsf{G}}\left(\mathbf{X}_{2}\right)$.

In the proof, we use the following Theorem, which adapts classical results about the universality of MLPs.

###### Theorem C.2 .

If a class of neural networks has compact uniform approximation over $\mathcal{G}$ with respect to linear functions and contains non-polynomial activations, then it has compact universal approximation over $\mathcal{G}$ with respect to continuous functions.

###### Proof.

Classical theorems such as Theorem 1 in [] allow us to approximate any continuous function over a compact set by composing a linear map $C$, an activation $\sigma$, and an affine map $A\cdot+b$. By assumption, we can implement the linear map, the activation, and the affine map; hence, by composing them, we can approximate any continuous function over the compact set.
∎

We are now ready to prove the paper’s main result: that, given injective positional encodings, BuNNs are compact universal approximators.

###### Theorem 5.3 .

Let $\mathcal{G}$ be a set of connected graphs with injective positional encodings. Then there is a $2$-layer deep BuNN with compact uniform approximation over $\mathcal{G}$.

In particular, let $\epsilon>0$ and $h$ be a feature transformation supported on $\bigsqcup_{\mathsf{G}\in\mathcal{K}}K^{n_{\mathsf{G}}}\subseteq\bigsqcup_{
\mathsf{G}\in\mathcal{G}}\mathbb{R}^{n_{\mathsf{G}}d}$ with $\mathcal{K}\subseteq\mathcal{G}$ finite and $K\subseteq\mathbb{R}^{d}$ a compact set, then there is a $2$-layer deep BuNN with width $\mathcal{O}\left(\sum_{\mathsf{G}\in\mathcal{K}}|\mathsf{V}_{\mathsf{G}}|\right)$ that approximates $h$ with uniform error $<\epsilon$.

###### Proof.

Reducing to linear approximation. It suffices to show that a BuNN layer can approximate any linear feature transformation $L$ because we can apply classical results such as Theorem C.2 to get universal approximation of $2$-layer deep networks with activation. Following Definition 5.1, we aim to show that we can approximate a linear feature transformation $L$ on any compact subset. For this, we fix $\epsilon>0$, the finite subset $\mathcal{K}\subseteq\mathcal{G}$, and compact feature space $K\subseteq\mathbb{R}^{c}$. In fact, we assume that $K=\mathbb{R}^{c}$ since approximating a linear map on any compact feature space is equivalent to approximating it on the whole space because a linear map defined on a neighborhood of the 0 0 vector can be extended uniquely to the whole vector space. Our goal is therefore to find a parameterization of a single BuNN layer such that for any graph $\mathsf{G}\in\mathcal{K}$ and for any input feature $\mathbf{X}\in\mathbb{R}^{n_{\mathsf{G}}c}$, we have $\|L_{\mathsf{G}}\left(X\right)-\operatorname{BuNN}_{\mathsf{G}}\left(\mathbf{X
}\right)\|_{\infty}<\epsilon$. We will show that $L$ can be parameterized exactly. Since $L$ is linear, it suffices to find a linear BuNN layer that satisfies for any $\mathsf{G}\in\mathcal{K}$ and any $\mathbf{X}\in\mathbb{R}^{n_{\mathsf{G}}c}$, $\frac{\partial\left(L\left(\mathbf{X}\right)\right)_{u}}{\partial\mathbf{x}_{v
}}=\frac{\partial\left(\operatorname{BuNN}\mathbf{X}\right)_{u}}{\partial
\mathbf{x}_{v}}$. By Lemma 4.2, we have $\frac{\partial\operatorname{BuNN}(\mathbf{X})_{u}}{\partial\mathbf{x}_{v}}=
\mathcal{H}(t,u,v)\mathbf{O}_{u}\mathbf{W}\mathbf{O}_{v}^{T}$. Hence, since MLPs are universal and the positional encodings are injective, it suffices to find bundle maps $\mathbf{O}:\mathsf{V}\rightarrow O\left(k\right)$ and $\mathbf{W}$ such that $\frac{1}{n_{\mathsf{G}}d_{u}}\sum_{v\in\mathsf{V}}\mathbf{O}_{u}^{T}\mathbf{W}
\mathbf{O}_{v}=\frac{\partial\left(LX\right)_{u}}{\partial X_{v}}$ for every $u,v\in\mathsf{V}$.

Defining the encoder and decoder: In order to find such a BuNN, we first need a linear encoder $\operatorname{lift}:\mathbb{R}^{c}\rightarrow\mathbb{R}^{2ck}$ which will be applied at every node before applying a $2ck$ dimensional BuNN layer. The lifting transformation maps each node vector $\mathbf{X}_{u}$ to the concatenation of $k$ vectors $\mathbf{X}_{u}$ interleaved with $k$ vectors $\vec{\mathbf{0}}\in\mathbb{R}^{c}$. This is equivalent to the linear transformation given by left multiplication by $\left(\mathbf{I}_{c\times c},\mathbf{0},\ldots,\mathbf{I}_{c\times c},\mathbf{
0}\right)^{T}\in\mathbb{R}^{2ck\times c}$. After the $2ck$ dimensional BuNN network, we will also need a linear decoder $\operatorname{pool}:\mathbb{R}^{2ck}\rightarrow\mathbb{R}^{c}$ applied to every node individually, which is the sum of the $k$ different $c$-dimensional vectors that are at even indices. This is equivalent to left multiplication by the matrix $\left(\mathbf{I}_{c\times c},\mathbf{0}_{c\times c},\ldots,\mathbf{I}_{c\times
c
},\mathbf{0}_{c\times c}\right)\in\mathbb{R}^{c\times 2ck}$. These two can be seen as a linear encoder and linear decoder, often used in practical GNN implementations. We prove the result by adding the lifting and pooling layers and using the higher dimensional $\widehat{BuNN}$ layer, so $\operatorname{BuNN}=\operatorname{pool}\circ\widehat{\operatorname{BuNN}}\circ
\operatorname{lift}$.

Defining the ‘universal bundle’: We fix $k=\sum_{\mathsf{G}\in\mathcal{K}}\left|\mathsf{V}_{\mathsf{G}}\right|$, so we can interpret our embedding space as a lookup table where each index corresponds to a node $v\in\bigsqcup_{\mathsf{G}\in\mathcal{K}}V_{\mathsf{G}}$. In turn, we can think of the parameter matrix $\mathbf{W}\in\mathbb{R}^{\left(\sum_{G\in\mathcal{K}}\left|\mathsf{V}_{\mathsf
{G}}\right|\right)\times\left(\sum_{G\in\mathcal{K}}\left|\mathsf{V}_{\mathsf{
G}}\right|\right)}$ as a lookup table where each entry corresponds to a pair of nodes in our dataset $\mathcal{K}$. Still thinking of the indices of the $2ck$ dimensions as $2c$-dimensional vectors indexed by the $k$ nodes in our dataset, we define $\mathbf{O}_{u}\in O\left(2ck\right)$ as a block diagonal matrix with $k$ different $2c$-dimensional blocks. These are all set to the identity except for the block at the index corresponding to node $u$, which is defined as
$\begin{pmatrix}\mathbf{\mathbf{0}_{c\times c}}&\mathbf{I}_{c\times c}\\
\mathbf{I}_{c\times c}&\mathbf{\mathbf{0}_{c\times c}}\end{pmatrix}$
which is a $2c\times 2c$ matrix that acts by permuting the first $c$ dimensions with the second $c$ dimensions.

Computing the partial derivatives. Since our model $\operatorname{BuNN}$ is a composition of linear maps, and since the maps $pool$ and $lift$ are applied node-wise, we get

$$ $\displaystyle\begin{split}&\frac{\partial\left(\operatorname{BuNN}\left( \mathbf{X}\right)\right)_{u}}{\partial\mathbf{x}_{v}}\\ &=\operatorname{pool}\frac{\partial\left(\widehat{\operatorname{BuNN}}\left( \operatorname{lift}\left(\mathbf{X}\right)\right)\right)_{u}}{\partial \operatorname{lift}\left(\mathbf{X}_{v}\right)}\operatorname{lift}\\ &=\left(\mathbf{I}_{c\times c},\mathbf{0}_{c\times c},\ldots,\mathbf{I}_{c \times c},\mathbf{0}_{c\times c}\right)\mathcal{H}(t,u,v)\mathbf{O}_{u}^{T} \mathbf{W}\mathbf{O}_{v}\left(\mathbf{I}_{c\times c},\mathbf{0}_{c\times c}, \ldots,\mathbf{I}_{c\times c},\mathbf{0}_{c\times c}\right)^{T}\\ &=\mathcal{H}(t,u,v)\sum_{1\leq k_{1},\ k_{2}\leq k}\left(\mathbf{I}_{c\times c },\mathbf{0}_{c\times c}\right)\mathbf{O}_{u}^{T}\mathbf{W}\mathbf{O}_{v}\left (\mathbf{I}_{c\times c},\mathbf{0}_{c\times c}\right)^{T}\end{split}$ $$

We proceed by partitioning the indexing by $(k_{1},\ k_{2})$ into four cases. The first case is $C_{1}=\left\{(k_{1},k_{2})\text{ such that }\left(k_{1}\neq u,v\text{ and }k_{
2}\neq u,v\right)\text{or}(k_{1}=k_{2}=u\text{ or }k_{1}=k_{2}=v)\right\}$ for which both $\mathbf{O}_{u}$ and $\mathbf{O}_{v}$ act like the identity. The second case is $C_{2}=\left\{(k_{1},k_{2})\text{ such that }k_{1}=u\text{ and }k_{2}\neq u,v\right\}$ where $\mathbf{O}_{u}$ flips the first $c$ rows with the second $c$ rows and $\mathbf{O}_{v}$ acts like the identity. $C_{3}=\left\{(k_{1},k_{2})\text{ such that }k_{2}=j\text{ and }k_{1}\neq u,v\right\}$ where $\mathbf{O}_{v}$ flips the first $c$ columns with the second $c$ columns, and $\mathbf{O}_{u}$ acts like the identity on the rows. Finally, the last case is when $k_{1}=u$ and $k_{v}=v$ in which $\mathbf{O}_{v}$ flips the rows and $\mathbf{O}_{v}$ flips the columns.

$$ $\displaystyle\begin{split}\ldots&=\mathcal{H}(t,u,v)\sum_{1\leq k_{1},\ k_{2} \leq k}\left(\mathbf{I}_{c\times c},\mathbf{0}_{c\times c}\right)\mathbf{O}_{u }^{T}\mathbf{W}\mathbf{O}_{v}\left(\mathbf{I}_{c\times c},\mathbf{0}_{c\times c }\right)^{T}\\ &=\mathcal{H}(t,u,v)\left(\mathbf{I}_{c\times c},\mathbf{0}_{c\times c}\right) \left[\sum_{(k_{1},k_{2})\in C_{1}}\begin{pmatrix}\mathbf{W}_{00}^{k_{1}k_{2}} &\mathbf{W}_{01}^{k_{1}k_{2}}\\ \mathbf{W}_{10}^{k_{1}k_{2}}&\mathbf{W}_{11}^{k_{1}k_{2}}\\ \end{pmatrix}+\sum_{(k_{1},k_{2})\in C_{2}}\begin{pmatrix}\mathbf{W}_{10}^{k_{ 1}k_{2}}&\mathbf{W}_{11}^{k_{1}k_{2}}\\ \mathbf{W}_{00}^{k_{1}k_{2}}&\mathbf{W}_{01}^{k_{1}k_{2}}\\ \end{pmatrix}\right.\\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left.+\sum_{ (k_{1},k_{2})\in C_{3}}\begin{pmatrix}\mathbf{W}_{01}^{k_{1}k_{2}}&\mathbf{W}_ {00}^{k_{1}k_{2}}\\ \mathbf{W}_{11}^{k_{1}k_{2}}&\mathbf{W}_{10}^{k_{1}k_{2}}\\ \end{pmatrix}+\begin{pmatrix}\mathbf{W}_{11}^{ij}&\mathbf{W}_{10}^{ij}\\ \mathbf{W}_{01}^{ij}&\mathbf{W}_{00}^{ij}\\ \end{pmatrix}\right]\left(\mathbf{I}_{c\times c},\mathbf{0}_{c\times c}\right) ^{T}\end{split}$ $$

Finally, after applying $\left(\mathbf{I}_{c\times c},\mathbf{0}_{c\times c}\right)$ on the left and $\left(\mathbf{I}_{c\times c},\mathbf{0}_{c\times c}\right)^{T}$ on the right (an operation that selects the upper left $c\times c$ block), we observe that setting all $\mathbf{W}_{00}^{k_{1}k_{2}}=\mathbf{W}_{01}^{k_{1}k_{2}}=\mathbf{W}_{10}^{k_{
1}k_{2}}$ to $\mathbf{0}_{c\times c}$ and setting $\mathbf{W}_{11}^{uv}:=\frac{1}{\mathcal{H}(t,u,v)}\frac{\partial{\left(L
\mathbf{X}\right)_{u}}}{\partial\mathbf{x}_{v}}$ if the nodes corresponding to $k_{1}$ and $k_{2}$ lie in the same graph and $\mathbf{0}_{c\times c}$ otherwise. This allows us to conclude that any linear layer can be parameterized.
∎

###### Proposition 5.2 .

###### Proof.

###### Definition C.1 .

###### Theorem C.2 .

###### Proof.

###### Theorem 5.3 .

###### Proof.

## Appendix D Discussion on compact uniform approximation versus uniform approximation

A strong definition of expressivity that deals with infinite collections of graphs was proposed in Rosenbluth et al. []. This definition subsumes graph-isomorphism testing (where the input feature on graphs is constant) and is finer than graph-level function approximation since it works at the node level. Furthermore, it also deals with infinite families of graphs, as opposed to most mainstream theorems of GNN expressivity, which are proved for graphs of bounded size (e.g. Azizian and marc lelarge [], Geerts and Reutter []) (see Section 2 for the notation and definition of features transformations).

###### Definition D.1 (From Rosenbluth et al. [ 27 ] ) .

Let $c,\ c^{\prime}\in\mathbb{N}$ and take $\mathbb{R}$ as feature space. Consider a collection of graphs $\mathcal{G}$. Let $\Omega\subseteq\mathcal{C}\left(\mathcal{G},\mathbb{R}^{c},\mathbb{R}^{c^{
\prime}}\right)$ be a set of feature transformations over $\mathcal{G}$, and let $H\in\mathcal{C}\left(\mathcal{G},\mathbb{R}^{c},\mathbb{R}^{c^{\prime}}\right)$ a feature transformation over $\mathcal{G}$. We say that $\Omega$ uniformly additively approximates $H$, notated $\Omega\approx H$ if $\forall\epsilon>0\ \forall\text{ compact }K^{n}\subset\mathbb{R}^{nc}\ \exists
F
\in\Omega\ \text{such that:},\ \forall\mathsf{G}\in\mathcal{G}\ \forall X\in K
^{n_{\mathsf{G}}c}\ \left\|F_{\mathsf{G}}\left(X\right)-H_{\mathsf{G}}\left(X
\right)\right\|_{\infty}\leq\epsilon$
where the sup norm $\|\cdot\|_{\infty}$ is taken over all nodes and dimensions of $n_{\mathsf{G}}c^{\prime}$ dimensional output.

Note that this definition differs from our Definition 5.1 in that it requires uniform approximation over all graphs in $\mathcal{G}$ simultaneously, while we allow the width to vary with the finite subset $\mathcal{K}\subseteq\mathcal{G}$, similar to how classical results allow the width to vary with the compact set over which to approximate the function. Such a definition has proven useful in Rosenbluth et al. [] to distinguish different aggregation functions and in Rosenbluth et al. [] to distinguish MPNNs with virtual nodes from Transformers. However, we argue that the definition above is too strong for a finite parameter GNN. This is because it requires uniform approximation over a *non-compact set*, which contrasts with traditional work on expressivity and is generally unfeasible and impractical. Indeed, finite-parameters MLPs are not universal over the whole domain $\mathbb{R}$ under the $\ell_{\infty}$-norm. On an infinite collection of featured graphs, the topology is the disjoint union topology on $\bigsqcup_{\mathsf{G}\in\mathcal{G}}\mathbb{R}^{n_{\mathsf{G}}d}$, a compact subset consists of a finite set of graphs, and for each graph $\mathsf{G}$ only non-zero on a compact subset of $\mathbb{R}^{nd}$. For these reasons, we introduce Definition 5.1, which is still rich enough to distinguish between BuNNs and MPNNs.

###### Definition D.1 (From Rosenbluth et al. [ 27 ] ) .

## Appendix E Algorithmic and Implementation details

In this Section, we provide more details on the implementation of BuNNs. We start by discussing how to use several vector-field channels when the input dimension is greater than the bundle dimension. We then discuss how to use several bundles at once when a single bundle is insufficient. We then combine both views, namely having several vector-field channels on several bundles at once. Finally, we describe how we compute our bundle maps in the experiments.

Extending to several vector-field channels. When the signal dimension exceeds the bundle dimension, i.e. $c>d$, we cannot directly apply BuNNs to the input signal. In that case, we first transform the signal into a hidden dimension, a multiple of the bundle dimension, i.e.$c=dp$. We reshape the input signal into $p$ channels of $d$-dimensional vector fields, where we apply the diffusion step (Equation 3) on each $p$ channels simultaneously, and we apply the weight matrix $\mathbf{W}\in\mathbb{R}^{dp\times dp}$ by first flattening the node signals into $dp$ dimensions, then multiplying by $\mathbf{W}$, and then reshaping it into $p$ channels of $d$ dimensional vector fields.

Extending to several bundles. Learning a high dimensional orthogonal matrix O($d$) becomes expensive since the manifold of orthogonal matrices is $\frac{d(d-1)}{2}$ dimensional. However, we can compute many low-dimensional bundles in parallel. In practice, we found that using several $2$-dimensional bundles was enough. Computing $b$ different $2$-dimensional bundles requires only $b$-parameters since the manifold O($2$) is $1$-dimensional. We, therefore, also use different ‘bundle channels’ given by an additional hyper-parameter – the number of bundles, which we denote $b$. Given an input signal of dimension $c=db$, we can decompose the signal into $b$ bundle channels of dimension $d$. We can compute the diffusion step (Equation 3) for each bundle in parallel. For the update step (Equation 2), we apply the weight matrix $\mathbf{W}\in\mathbb{R}^{bd\times bd}$ by first flattening the node signals into $bd$ dimensions, then multiplying by $\mathbf{W}$, and then reshaping it into $b$ bundle channels of $d$ dimensional vector fields over $b$ different bundle structures.

###### Remark E.1 .

We note that using $b$ different $d$ dimensional bundles is equivalent to parameterizing a subset of one $bd$-dimensional structure, consisting of the orthogonal map $\mathbf{O}\in O(bd)\subset\mathbb{R}^{bd\times bd}$ that are block diagonal matrices $\mathbf{O}=\bigoplus_{i=1\ldots b}\mathbf{O}_{i}$,
with each $\mathbf{O}_{i}\in O(d)$.

Extending to several bundles and vector-field channels. We can combine the above two observations. Given an input signal of dimension $c=bdp$, we can subdivide this into $b$ different bundle structures of dimension $d$ and $p$ channels for each bundle. We diffuse on the appropriate bundle structure and flatten the vector fields into a $c\times c$ vector before applying the learnable parameters.

Computing the bundle maps. In our experiments, we noticed that having several bundles of dimension $2$ was more efficient than one bundle of large dimensions, while there was no clear performance gain when using higher dimensional bundles. To compute the $b$ bundle maps $\mathbf{O}_{v}$ we therefore only need $b$ rotation angles $\theta_{v}$, one per bundle. In our experiments, we use Housholder reflections using the python package noa [] or direct parameterization. For direct parameterization, we do the following: since the matrix group $O(2)$ is disconnected, we always take $b$ to be even and parameterize half the bundles as rotation matrices $r(\theta)=\begin{pmatrix}\cos\left(\theta\right)&\sin\left(\theta\right)\\
-\sin\left(\theta\right)&\cos\left(\theta\right)\end{pmatrix}$ and the other half to correspond to matrices with determinant $-1$, which can be parameterized by $r^{*}(\theta)=\begin{pmatrix}\cos\left(\theta\right)&\sin\left(\theta\right)\\
\sin\left(\theta\right)&-\cos\left(\theta\right)\end{pmatrix}$. We compute the angles $\theta$ as in Equation 1 where the network $\boldsymbol{\phi}^{\left(\ell\right)}$ is either an MLP or a GNN. The network $\boldsymbol{\phi}$ is either shared across layers or differing at every layer.

Taylor approximation algorithm. We now provide pseudo-code on how we implement Equations 2, and  3. We then proceed with a complexity analysis. The key idea of the algorithm is that the bundle heat kernel can be approximated efficiently using the standard graph heat kernel.

Figure: Algorithm 1 Taylor expansion implementation of a BuNN layer

The complexity of the algorithms is as follows. There are $3$ matrix-vector multiplications done at each node in lines $1$, $2$, and $8$, which are done in $\mathcal{O}\left(3d^{2}|\mathsf{V}|\right)$. The for loops consist of matrix-matrix multiplications, which are done in $\mathcal{O}\left(d|\mathsf{E}|\right)$ with sparse matrix-vector multiplication. The memory complexity is $\mathcal{O}\left((d+d^{2})|\mathsf{V}|\right)$ since we need to store $d$ dimensional vectors and the orthogonal maps for each node. The exact implementation is described in Algorithm 1

Spectral method. We now describe how to implement a BuNN layer using the eigenvectors and eigenvalues of the Laplacian.

Figure: Algorithm 2 Spectral implementation of a BuNN layer

###### Remark E.1 .

### E.1 Householder reflections.

Many different parameterizations of the group $O(n)$ exist. While direct parameterizations are possible for $n=2,\ 3$ it becomes increasingly complex to do so for larger $n$, and a general method working for all $n$ is a desirata. While there are several methods to do so, we use Householder reflection since it is used in related methods []. We use the Pytorch package from [obukhov2021torchhouseholder]. Given given $k$ vectors $v_{i}\in\mathbb{R}^{d}$, define the Householder matrices as $H_{i}=I-2\frac{v_{i}v_{i}^{T}}{\|v_{i}\|_{2}^{2}}$, and define $U=\prod_{i=1}^{k}H_{i}$. All orthogonal matrices may be obtained using the product of $d$ such matrices. Hence the map $\mathbf{R}^{d\times d}\rightarrow O(d)$ mapping $V=(v_{i})$ to $U$ is a parameterization of the orthogonal group. We use pytorch implementations allowing autograd provided in [obukhov2021torchhouseholder].

### E.2 BuNN-Hop extension.

In practice we also explored a variation of BuNNs inspired by MixHop [] which we call BuNN-Hop (described in Algorithm 3). We found that a slightly different model performed better on the Peptides benchmark datasets. One layer of the model can be derived by running a discrete time diffusion process with update $\mathbf{X}(t)=\mathbf{A}_{\mathcal{B}}\mathbf{X}(t-1)$ and integrating over it using attention term $\boldsymbol{\alpha}(t)$. The resulting update equation is $\mathbf{Y}=\sum_{k=1}^{K}\boldsymbol{\alpha}_{k}\mathbf{A}^{k}_{\mathcal{B}}
\mathbf{X}$ where $\boldsymbol{\alpha}$ satisfy $\sum_{k=1}^{K}\boldsymbol{\alpha}_{k}=1$.

Figure: Algorithm 3 BuNN-Hop layer implementation.

## Appendix F Experiment details

In this Section we provide additional information about the experiments on the heterophilic graph benchmarks, the LRGB benchmarks, and the synthetic experiments. All experiments were ran on a cluster using NVIDIA A10 (24 GB) GPUs, each experiment using at most $1$ GPU. Each machine in the cluster has 64 cores of Intel(R) Xeon(R) Gold 6326 CPU at 2.90GHz, and $\sim$500GB of RAM available. The synthetic experiments from Section 6.1 were run on CPU and each run took roughly $20$ minutes. The heterophilic experiments from Section 6 were run GPU and varied between $5$ minutes to $1.5$ hours. The LRGB experiments were run on GPU and varied between $0.5$ hours and $4$ hours.

### F.1 LRGB: training and tuning.

For peptides-func and peptides-struct we use a fixed parameter budget of $\sim 500k$ as in Dwivedi et al. []. We fix hyper-parameters to be the best GCN hyper-parameters from Tönshoff et al. [], and tune only BuNN-specific parameters as well as the use of BatchNorm. In Table 5, we report the grid of hyper-parameters that we searched, and denote in bold the best combinations of hyper-parameters. The parameters fixed from Tönshoff et al. [] are the following:

- •
Dropout $0.1$
- •
Learning rate $0.001$
- •
Head depth $3$
- •
Positional Encoding: LapPE for struct and RWSE for func
- •
Optimizer: AdamW with a cosine annealing learning rate schedule and linear warmup.
- •
Batch size $200$
- •
Weight decay 0 0
- •
We use skip connection as implemented in Dwivedi et al. [] and not in Tönshoff et al. []. That is, the skip connection does not skip the non-linearity.

For the BuNN specific parameters, we use $2$ dimensional bundles, whose angles $\theta$ we compute with the help of a small SumGNN architecture using a sum aggregation as defined by $\mathbf{\theta}_{v}^{(\ell)}=\sigma\left(\mathbf{W}_{\text{s}}\mathbf{x}^{(
\ell)}_{v}+\mathbf{W}_{\text{n}}\sum_{u\in\mathcal{N}(v)}\mathbf{x}^{(\ell)}_{
u}\right)$ where the input dimension is the hidden dimension, the hidden dimension is twice the number of bundles and the output is the number of bundles. The number of SumGNN layers is a hyper-parameter we tune. When it is 0 0 we use a $2$ layer MLP with hidden dimension also twice number of bundles. For each hyper-parameter configuration, we set the hidden dimension in order to respect to the parameter budget. We use truncated Taylor series for this implementation.

**Table 5: Grid of hyper-parameters for peptides-func and peptides-struct.**
| Parameters | All Values | Best Values |  |  |  |
| --- | --- | --- | --- | --- | --- |
|  | BuNN | BuNN-Hop |  |  |  |
|  |  | func | struct | func | struct |
| Num bundles b | $4$, $8$, $16$ | $16$ | $16$ | $16$ | $16$ |
| Number of BuNN layers | $1\ -\ 6$ | $6$ | $4$ | $6$ | $4$ |
| Number of SumGNN layer | $0\ -\ 3$ | $1$ | 0 0 | $6$ | $4$ |
| Max expansion degree $K$ | $1$  $2$  $4$  $8$  $16$  $32$  $64$  $128$ | $128$ | $1$ | $64$ | $64$ |
| Time $t$ | $0.1$, $1$, $5$, $6$, $7$, $8$, $9$, $10$, $\infty$ | $1$ | $6$ | N/A | N/A |

### F.2 Heterophilic graphs: training and tuning.

For the heterophilic graphs we use the source code from Platonov et al. [] in which we add our layer definition. We report all training parameters that we have tuned. Namely, we use GELU activation functions, the Adam optimizer with learning rate $3\times 10^{-5}$, and train all models for $2000$ epochs and select the best epoch based on the validation set. To compute the bundle maps, we compute the parameters $\theta$ with a GraphSAGE architecture shared across layers ($\phi$ method = shared) or different at each layer ($\phi$ method = not shared), with hidden dimension dimension the number of bundles. The number of layers of this GNN is a hyper-parameter we tuned, which when set to 0 0 we use a $2$ layer MLP. For each task we manually tuned parameters, which are subsets of the combinations of parameters in the grid from Table 6. The implementation of the heat kernel used is either truncated Taylor series with degree $8$, or the Spectral implementation. We report the best performing combination of parameters in Table 7.

**Table 6: Parameters searched when tuning on the heterophilic graph benchmark datasets.**
| Parameters | All Values |
| --- | --- |
| Hidden dim | $256,512$ |
| Num bundles $b$ | $8,16,32,64,128,256$ |
| Bundle dimension $d$ | $2$ |
| Number of BuNN layers | $1\ -\ 8$ |
| Number of GNN layer $\phi$ | $0\ -\ 8$ |
| Taylor expansion degree $K$ | $8$ |
| Time $t$ | $0.1,1,10,100$ |
| $\phi$ method | shared, not shared |
| Dropout | $0.0,\ 0.2$ |
| Learning rate | $3\times 10^{-4},3\times 10^{-5}$ |

**Table 7: Best parameter for each dataset in the heterophilic graph benchmarks.**
| Parameters | Best Values |  |  |  |  |
| --- | --- | --- | --- | --- | --- |
|  | roman-empire | amazon-ratings | minesweeper | tolokers | questions |
| Hidden dim | $512$ | $512$ | $512$ | $512$ | $256$ |
| Num bundles b | $64$ | $64$ | $256$ | $256$ | $128$ |
| Bundle dim | $2$ | $2$ | $2$ | $2$ | $2$ |
| Number of BuNN layers | $6$ | $2$ | $8$ | $6$ | $6$ |
| Number of GNN layer | $8$ | 0 0 | $8$ | $7$ | $6$ |
| Taylor expansion degree $K$ | $N/A$ | $N/A$ | $8$ | $8$ | $8$ |
| Time $t$ | $100$ | $1$ | $1$ | $1$ | $1$ |
| $\phi$ method | not shared | not shared | not shared | not shared | shared |
| Dropout | $0.2$ | $0.2$ | $0.2$ | $0.2$ | $0.2$ |
| Learning rate | $3\times 10^{-4}$ | $3\times 10^{-4}$ | $3\times 10^{-5}$ | $3\times 10^{-5}$ | $3\times 10^{-5}$ |

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 2
questions: What price is paid for going rom the full cellular sheaves to the flat vector bundle in expressivity? In practice, in what settings is the flat vector bundle approach good enough? For example, it seems like that for sparse graphs the orthogonal maps associated with each node would be good enough.

Similarly, can the authors provide a comparison between the structure of the fixed points for cellular sheaves and flat vector bundles? It seems like some of the "richness" of the long time limit of diffusion is lost with the flat vector bundles. How does this reflect on expressivity of flat vector bundles compared with cellular sheaves?

Does the transformation learned for each node need to be orthogonal in practice? It seems that the approach is somewhat limited because the learned matrices must be orthogonal. It is difficult to learn matrices that are constrained to be orthogonal. Can the authors relax this requirement in practice and use invertible matrices or a fully unconstrained matrix?

It seems like a key aspect of the proposed approach is how the orthogonal transformations at each node are learned? The authors mention MLPs or even GNNs that take in the graph structure and positional encodings of the nodes. What type of positional encodings are best for learning the orthogonal transformations? Why are the node features not taken into account when computing the orthogonal map for each node? Are the authors trying to keep the graph geometry distinct from the feature information? 

In Table 3, how does the BuNN approach outperform the NSD for the data sets in which memory is not a limit for NSD? Is the depth of the models in BuNN larger helped by spectral methods? It seems like depth was 1 for the amazon-ratings data set. Does NSD overfit the data because of its use of cellular sheaves?
soundness: 3
strengths: This reviewer generally liked the paper. The idea of simplifying cellular sheaves using flat vector bundles is clever. This idea seems to solve the key problem with the computational complexity challenges faced by neural sheaf diffusion. Another strength of the paper is the clever synthetic toy models which demonstrate the gain from the proposed bundle architecture very clearly, especially when it comes to the over-squashing problem. The empirical performance of the proposed method is also impressive, especially on heterophilic graphs like minesweeper. The novel notion of expressivity introduced by the authors is also intriguing.
summary: This paper extends previous work on graph neural networks equipped with cellular sheaves by introducing a more efficient way of computing the heat diffusion over bundles. The key idea of the authors is to assume a simpler form for the sheaves, namely a flat vector bundle, which can be thought of as an orthogonal map associated with each node (as opposed to one associated with each node-edge pair). This assumption significantly simplifies the computations necessary to learn the bundles from data. In particular, the heat diffusion kernel can now be computed efficiently during the learning of the bundles. This enables the authors to use spectral methods to run the diffusion process for long periods of time addressing the over-squashing problem on graphs. This flat vector bundles also bring the strength of the cellular sheaves in that the fixed point of diffusion is no longer a trivial state where all nodes share the same features. Rather, the long time limit now has richer structure where variation is allowed across the nodes. This addresses the over-smoothing problem. The authors propose GNN architecture using flat vector bundles and heat kernels (computed both using Taylor expansions and spectral methods) and demonstrate improved performance on synthetic data sets and real data sets.
weaknesses: A weakness of this paper is that the central idea is arguably a simple extension of the neural sheaf diffusion paper. It is not clear that this constitute significant advance. This is balanced by the impressive empirical results in the paper.

The authors highlight well the advantages of replacing the cellular sheaves with flat vector bundles. However, it is not clear what is the cost of doing so. How does expressivity suffer with this assumption? Are there circumstances under which this simplification works better? For example, does the sparsity of the graph play role. More discussion on this would have helped the reader. The authors nicely show the existence of fixed points for the diffusion process under flat vector bundles and that features across nodes are related by orthogonal transformations. How does this compare with the structure of fixed points under general cellular sheaves? How does this impact the expressivity of the bundle neural networks compared with their sheave counterpart?

This reviewer struggled at times to understand how much the principled geometric approach brought over from topology and differential geometry contributed to formulation of the proposed GNN. Although interesting, can the author's approach be though of as learning a transformation on each node. With the formality relaxed, it is not clear if this transformation needs to be orthogonal. The authors seems to go to quite a bit of trouble to ensure that the learned matrices at each node are orthogonal. Of course, Riemannian manifolds require this, but is it really needed? Would invertible matrices done just as good of a job? It seems like in practice they are limited to using only 2-dimensional vector spaces on each node, with the transformation matrices parameterized explicitly to be orthogonal.

This reviewer was bothered at times by the strong use of language by the authors. The claim that the proposed bundle approach solves the over-squashing problem is over-stated. Any spectral method will address over-squashing. The authors have proposed an efficient way to compute the heat kernel more efficiently by replacing cellular sheaves with flat vector bundles. This allows them to run diffusion for longer times. There is no conceptual new solution to over-squashing.

---

confidence: 2
contribution: 2
initial_rating: 6
presentation: 2
questions: Please see weaknesses
soundness: 2
strengths: 1. The paper addresses an important problem of oversmoothing and oversquashing seen in GNNs via vector bundles, building on cellular sheafs in GNNs.
1. The use of flat vector bundles with orthogonal maps to reduce the computational complexity is interesting.
1. The paper is well-organized.
summary: The paper proposes a new graph neural network architecture called Bundle Neural Networks (BuNNs). The architecture is based on message diffusion on flat vector bundles which are topological structures that assign to each node a vector space and an orthogonal map. They seem to be inspired from Sheaf neural networks and, in different forms, are claimed to address the problems of overs-moothing and over-squashing. Theoretical analysis on the feature transformation with injective positional encodings shows uniform expressivity. Further analysis on the properties of BuNNs and experimental results on both synthetic and real-world datasets, are provided for validating the proposed model.
weaknesses: 1. The experiments to validate oversmoothing and oversquashing are not comprehensive. The authors can test the oversquashing on Neighbours match dataset [1], which has been used well in multiple papers. Further the authors can explore using effective resistance metric to measure oversquashing [3] and [4].
1. Related works which are designed via dirichlet energy optimization like [2] are not discussed. May help to discuss the difference with [2] in terms of results on oversmoothing.

### References:
[1]$~$ Uri Alon and Eran Yahav. On the bottleneck of graph neural networks and its practical implications. arXiv preprint arXiv:2006.05205, 2020.

[2]$~$ Fu, Guoji, et al. "Implicit graph neural diffusion based on constrained Dirichlet energy minimization." arXiv preprint arXiv:2308.03306 (2023).

[3]$~$ Mitchell Black, Zhengchao Wan, Amir Nayyeri, and Yusu Wang. Understanding oversquashing in gnns through the lens of effective resistance. In International Conference on Machine Learning, pages 2528–2547. PMLR, 2023.

[4]$~$ Dong, Yanfei, et al. "Differentiable Cluster Graph Neural Network." arXiv preprint arXiv:2405.16185 (2024).

---

confidence: 4
contribution: 2
initial_rating: 6
presentation: 2
questions: * See weaknesses
soundness: 3
strengths: The proposed method is sound and very general. BuNNs have many favorable and interesting properties. I particularly like the property regarding fixed points. The paper is nicely written and mostly understandable.
summary: This paper proposes to perform message-passing over a diffused graph using flat vector bundles. This is achieved by transforming each node state within its own vector space using a learned and distinct orthogonal map. This is the message that a node sends to each neighbor, where the inverse of the orthongonal map from the receiving node is used as a decoding mapping. The authors show that this can avoid over-smoothing, over-squashing, and can be a compact uniform approximation over a set of graphs. Experiments show the effectiveness of the proposed approach.
weaknesses: * Eq. 2: From the theoretical derivation, this step does not seem needed as $\mathcal{H}_B$ already applies $O_v$ and $O_u$. The only provided reason is that it helps to prove Theorem 5.3. Are there other reasons why $W$ and $b$ are introduced? What happens when these are removed?
* BuNNs and their properties are mostly compared to MPNNs applied to the original. It seems fairer to me to also compare its properties to Graph Transformers and infinite-depth MPNNs (e.g. [1]) (Table 1).
* Theorem 5.3: 
    * This statement is quite interesting but I do not fully understand it. In which cases are injective positional encodings achievable? Would it be satisfied when $\mathcal{G}$ is a set of WL distinguishable graphs and we use WL colors are positional encodings?
    * How does it deal with two nodes $u,v$ in the same graph with the same WL unfolding tree? By my understanding of the injective positional encodings, they could still have the same initial state. The MLP (Eq. 1) would produce the same $\mathbf{O}_v = \mathbf{O}_u$ for both nodes.
    * Which changes would be needed for graph transformers to also satisfy this universality?
* Table 2: A comparison with Graph Transformers or infinite-depth MPNNs would be insightful
* Table 3: It seems like the authors performed experiments for NSD and BuNN, but hyperparameter ranges and optimal hyperparameters are only presented for BuNN.
* Details about runtimes are missing. It would be important to have execution times for at least one experiment for BuNN for a comprehensive evaluation.

---
[1] Gu et al., Implicit graph neural networks, NeurIPS (2020).

Meta-review: The manuscript proposes a novel graphical neural network architecture that sends messages on flat vector bundles. The reviewers unaminously supports the manuscript for acceptance. During the discussion phase, the authors also address most of the concerns by the reviewers. Therefore, the metareviewer would recommend acceptance of the paper.

Decision: Accept (Spotlight)