Message passing neural networks (MPNNs) have been shown to have limitations in terms of expressivity and modeling long-range interactions. Augmenting MPNNs with a virtual node (VN) removes the locality constraint of the layer aggregation and has been found to improve performance on a range of benchmarks. We provide a comprehensive theoretical analysis of the role of VNs and benefits thereof, through the lenses of oversmoothing, oversquashing, and sensitivity analysis. First, in contrast to prior belief, we find that VNs typically avoid replicating anti-smoothing approaches to maintain expressive power. Second, we characterize, precisely, how the improvement afforded by VNs on the mixing abilities of the network and hence in mitigating oversquashing, depends on the underlying topology. Finally, we highlight that, unlike Graph-Transformers (GT), classical instantiations of the VN are often constrained to assign uniform importance to different nodes. Consequently, we propose a variant of VN with the same computational complexity, which can have different sensitivity to nodes based on the graph structure. We show that this is an extremely effective and computationally efficient baseline on graph-level tasks.

§ INTRODUCTION

Graph Neural Networks (GNNs) [Scarselli et al., 2009, Gori et al., 2005, Micheli, 2009] are a popular framework for learning on graphs. GNNs typically adopt a message passing paradigm [Gilmer et al., 2017], in which node features are aggregated over their local neighborhood recursively, resulting in architectures known as Message Passing Neural Networks (MPNNs). Whilst MPNNs have linear complexity in the number of edges and an often beneficial locality bias, they have been shown to have several limitations. In terms of expressivity, they are at most as powerful as the Weisfeiler-Leman test [Weisfeiler and Leman, 1968, Morris et al., 2019, Xu et al., 2018] and cannot count certain substructures [Chen et al., 2020, Bouritsas et al., 2023]. Besides, the repeated aggregation of local messages can either make  node representations indistinguishable—a process known as `oversmoothing' [Nt and Maehara, 2019, Oono and Suzuki, 2020]—or limit the models' ability to capture long-range dependencies—a phenomenon called `oversquashing' [Alon and Yahav, 2021].

In particular, oversquashing strongly depends on the graph-topology, since it describes the inability of MPNNs to reliably exchange information between pairs of nodes with large commute time [Di Giovanni et al., 2023].
To overcome this limitation, methods have moved away from only performing aggregation over the local neighborhood and instead have altered the graph connectivity to reduce its commute time. Common examples include variations of multi-hop message-passing <cit.>, adding global descriptors [Horn et al., 2022], or rewiring operations based on spatial or spectral quantities [Topping et al., 2021, Arnaiz-Rodríguez et al., 2022]. A notable case is given by Graph-Transformers (GTs) [Kreuzer et al., 2021, Ying et al., 2021], where global attention is used to weight the messages of each node pair.

Virtual Nodes. While GTs incur a large memory cost, a significantly more efficient method that incorporates global information in each layer to combat oversquashing, is adding a virtual node (VN) connected to all nodes in the graph [Pham et al., 2017]. This is widely used in practice [Gilmer et al., 2017, Battaglia et al., 2018] with empirical successes [Hu et al., 2020, Hwang et al., 2022]. However, besides recent work connecting VNs to GTs [Cai et al., 2023, Rosenbluth et al., 2024], few efforts have been made to analyze in what way and to which extent VN improves the underlying model.

Contributions. We propose a theoretical study of VN through the lenses of oversmoothing, oversquashing and sensitivity analysis. Our analysis sheds light on the connection between VN and existing anti-smoothing techniques, the impact VN has on reducing oversquashing by affecting the commute time, and the gap between VN and GTs in terms of the sensitivity to other nodes' features. Specifically, our main contributions are the following:

In <Ref> we contest prior belief that a reason for the success of MPNN + VN is their ability to simulate methods which prevent oversmoothing. We tackle this point by comparing the expressivity of MPNN + VN and anti-smoothing techniques through polynomial filters.

* In <Ref> we provide the first, systematic study of the impact of VN on the oversquashing phenomenon. We show that, differently from more general rewiring techniques such as GTs, the improvement brought by VN to the mixing abilities of the network, can be characterized and bounded in terms of the spectrum of the input graph.

* Relying on sensitivity analysis of node features, in <Ref> we find a gap between VN and GT in their ability to capture heterogeneous node importance. Our theory motivates a new formulation of MPNN + VN, simply denoted by MPNN + VN$_G$, that better leverages the graph structure to learn a heterogeneous measure of node relevance, at no additional cost.

* Finally, in <Ref>, we first validate our theoretical analysis through extensive ablations and experiments. Next, we evaluate MPNN + VN$_G$ and show that it consistently surpasses the baseline MPNN + VN, precisely on those tasks where node heterogeneity matters.

§ PRELIMINARIES AND RELATED WORK

Message Passing Neural Networks
Let $G = (V, E)$ be an unweighted, undirected and connected graph with node set $V$ and edge set $E$. The connectivity of the graph is encoded in the adjacency matrix $\textbf{A} \in \mathbb{R}^{n \times n}$, where $n = |V|$, and the 1-hop neighborhood of a node is $N_i = \{j \in V : (j, i) \in E\}$. Furthermore, node features $\{h_i:i\in V\}\subset\R^d$ are provided as input. Message passing layers are defined using learnable update and aggregation functions $\mathsf{up}$ and $\mathsf{agg}$ as follows

\begin{align}\label{eq:MPNN}
h^{(\ell+1)}_i = \mathsf{up}^{(\ell)} ( h^{(\ell)}_i, ~\mathsf{agg}^{(\ell)} ( \{h^{(\ell)}_j : j \in N_i\})),
\end{align}

where $h^{(\ell)}_i$ is the feature of node $i$ at layer $\ell$, $\mathsf{up}$ is a learnable update function (typically an MLP) and $\mathsf{agg}$ is invariant to node permutations. Below, we refer to models described as in (<ref>) as MPNNs.

Limitations of MPNNs
Aggregating information over the 1-hop of each node iteratively, leads to systematic limitations. MPNNs without node identifiers [Loukas, 2019] are not universal [Maron et al., 2018] and are at most as powerful as the 1-WL test in distinguishing graphs [Morris et al., 2019, Xu et al., 2018].
Additionally, when the task depends on features located at distant nodes (i.e. long-range interactions), MPNNs suffer from oversquashing, defined as the compression of information from exponentially growing receptive fields into fixed-size vectors [Alon and Yahav, 2021]. In fact, oversquashing limits the MPNN's ability to solve tasks with strong interactions among nodes at large commute time [Di Giovanni et al., 2023, Black et al., 2023, Di Giovanni et al., 2023]. Finally, MPNNs may incur oversmoothing, whereby node features attain the same value as the number of layers increases [Nt and Maehara, 2019, Oono and Suzuki, 2020, Cai and Wang, 2020].

Graph Transformers (GTs) An extension of the MPNN-class is given by Graph-Transformers (GTs) [Kreuzer et al., 2021, Rampášek et al., 2022, Ying et al., 2021, Hussain et al., 2022, Wu et al., 2021, Chen et al., 2022], where a global attention mechanism connecting all pairs of nodes augments (or entirely replaces [Ma et al., 2023]) the 1-hop aggregation in MPNNs. GTs increase the expressive power of MPNNs by heavily relying on structural and positional encodings [Dwivedi et al., 2021], and entirely bypass the issue of oversquashing since the global attention breaks all potential bottlenecks arising from the graph topology. This comes at a price though, with GTs incurring quadratic memory cost and being sensitive to the choice of positional (structural) encodings.

MPNNs With Virtual Nodes An alternative approach to enhancing MPNNs without leading to the significant memory costs of GTs, consists in introducing a virtual node (VN) connected to all other nodes in the graph [Gilmer et al., 2017, Battaglia et al., 2018]. A standard formulation of message passing with virtual nodes [Cai et al., 2023], referred to as MPNN + VN in the following, is:

\begin{align}\label{eq:mpnn+VN}
h^{(\ell+1)}_{\mathsf{vn}} &= \mathsf{up}^{(\ell)} ( h^{(\ell)}_{\mathsf{vn}}, ~\mathsf{agg}_{\mathsf{vn}}^{(\ell)} ( \{h^{(\ell)}_j : j \in V\})), \\
h^{(\ell+1)}_i &= \mathsf{up}^{(\ell)} ( h^{(\ell)}_i, ~\mathsf{agg}^{(\ell)} ( \{h^{(\ell)}_j : j \in N_i\}), h^{(\ell)}_{\mathsf{vn}}). \notag
\end{align}

The class MPNN + VN extends MPNNs in (<ref>) to a multi-relational graph $G_{\mathsf{vn}}$ obtained by adding a node VN connected to all $i\in V$, where any new edge is distinguished from the ones in the input graph.

Understanding VN: Open Questions and Challenges Augmenting MPNNs with VN has often been shown to improve performance [Hwang et al., 2022, Sestak et al., 2024]. Nonetheless, a careful analysis of how and when VN helps is still lacking. Key related work is Cai et al., 2023, which relied on the universal approximation property of DeepSets [Segol and Lipman, 2019] to show that MPNN + VN with hidden state dimensions of $O(n^d)$ and $O(1)$ layers can approximate a self-attention layer. Their argument though discards the topology and holds in a non-uniform regime, raising questions about its practical implications. Additionally, Cai et al., 2023 conjectured that MPNN + VN benefits from being able to simulate strategies for mitigating oversmoothing such as PairNorm [Zhao and Akoglu, 2019]. Crucially though, VN is mainly added to networks so to model non-local interactions, as follows from (<ref>) where any node pair is separated by at most 2 hops through VN. Nonetheless, neither a formal analysis on the extent to which MPNN + VN mitigates oversquashing, nor a fine-grained comparison between GTs and MPNN + VN from the point of view of sensitivity analysis (i.e. interactions among node features in a layer), have been conducted so far.

Outline of This Work We provide a systematic study of MPNN + VN through the lenses of oversmoothing, oversquashing, and sensitivity analysis. We show that: In contrast to prior belief, the ability of VN to emulate anti-smoothing approaches is likely not beneficial (<ref>); The improvement brought by VN in mitigating oversquashing depends on the spectrum of the graph (<ref>); Differently from GTs, the sensitivity of VN to distinct node features is often uniform, and hence we propose VN$_G$, a new formulation of VN that can learn heterogeneous node importance at no additional cost (<ref>); The framework VN$_G$ consistently improves over VN, closing the empirical gap with GTs (<ref>).

Graph-Level Tasks In this work, we focus on graph-level tasks; this way, we can compare the mixing abilities of MPNN + VN with that of MPNN, as per the graph-level analysis in Di Giovanni et al., 2023 (<ref>). Graph-level benchmarks (both classification and regression) offer a comprehensive test bed to compare VN and GTs, which is aligned with similar evaluations [Gutteridge et al., 2023, Barbero et al., 2023]. Nonetheless, we emphasize that the statements in <Ref> hold layerwise, and hence do not depend on whether the task is graph-level or not; additionally, we can extend the results in <Ref> to node-level tasks following <cit.>.

§ MPNN + VN AND SMOOTHING

Despite the empirical success of MPNN + VN, the reasons behind these improvements still lack a thorough understanding. To this aim, Cai et al., 2023 postulated that MPNN + VN is also beneficial due to its ability to replicate PairNorm [Zhao and Akoglu, 2019], a framework designed to mitigate oversmoothing. In this section we provide theoretical arguments against this conjecture, proving that MPNN + VN loses expressive power when it emulates PairNorm.

Comparison to PairNorm Methods such as PairNorm [Zhao and Akoglu, 2019], aim at mitigating oversmoothing, which occurs when the node representations converge to the same vector in the limit of many layers, independent of the task [Nt and Maehara, 2019, Oono and Suzuki, 2020]. To address the claims in Cai et al., 2023 and compare MPNN + VN and PairNorm, in this section we adopt the assumptions in [Zhao and Akoglu, 2019] and other theoretical works on oversmoothing [Di Giovanni et al., 2023], and consider networks whose layers are linear.
In particular, to assess advantages of MPNN + VN over PairNorm, we study their expressivity. Since aggregating $m$ linear layers yields a polynomial filter of the form $p(\mathbf{A})\mathbf{H}$, where $\mathbf{H}$ are the input node features, we choose as expressivity metric precisely the ability of a model to learn such polynomial filters.

In fact, under mild assumptions, polynomial filters are comparable to the 1-WL test [Wang and Zhang, 2022]. Additionally, when the input features $\mathbf{H}$ are constant, polynomial filters (<ref>) coincide with graph moments, a metric of expressivity studied in Dehmamy et al., 2019, further validating our approach.

Measuring Expressivity With Polynomial Filters    Since in this section we restrict the analysis to models that have linear layers, we introduce a simple MPNN baseline, which is obtained from (<ref>) by choosing a sum aggregation and a linear update with a residual connection:

\begin{equation}\label{eq:MPNN_lin}
\mathbf{H}^{(\ell + 1)} = \mathbf{H}^{(\ell)} + \mathbf{A}\mathbf{H}^{(\ell)}\mathbf{W}^{(\ell)},
\end{equation}

where $\mathbf{W}^{(\ell)}$ is a $d\times d$ weight matrix and $\mathbf{H}^{(\ell)}$ is the $n\times d$ feature matrix at layer $\ell$. Similarly, for MPNN + VN, we consider a layer update with mean aggregation for VN, akin to Cai et al., 2023:

\begin{equation}\label{eq:MPNN_VN_lin}
\mathbf{H}^{(\ell + 1)} = \mathbf{H}^{(\ell)} + \mathbf{A}\mathbf{H}^{(\ell)}\mathbf{W}^{(\ell)} + \frac{1}{n}\mathbf{1}\mathbf{1}^\top\mathbf{H}^{(\ell-1)}\mathbf{Q}^{(\ell)},
\end{equation}

with $\mathbf{Q}^{(\ell)}$ a weight matrix and $\mathbf{1}\in\R^n$ the vector of ones. For simplicity, we assume that the output after $m$ layers $\mathbf{Y}_m$ also has dimension $d$ and choose a mean pooling operation in the final layer, so that given a weight matrix $\boldsymbol{\Theta}$ we have $ \mathbf{Y}_m = \boldsymbol{\Theta}\mathsf{Mean}(\mathbf{H}^{(m)})$, with $\mathbf{H}^{(m)}$ computed by either (<ref>) or (<ref>). We can express $\mathbf{Y}_m$ as a graph-level polynomial filter in terms of $d\times d$ weight matrices $\{\boldsymbol{\Theta}_k\}$:

\begin{equation}\label{eq:poly_mean}
\mathbf{Y}_m  = \mathsf{Mean}\big(\sum_{k = 0}^{m}\mathbf{A}^k \mathbf{H}\boldsymbol{\Theta}_k\big)^\top.
\end{equation}

We note now that MPNN + VN (<ref>) recovers PairNorm when $\mathbf{Q}^{(\ell)} = -\mathbf{I}$, with $\mathbf{I}$ the $d\times d$-identity matrix. We denote this special choice as MPNN - VN. The proof of the following Theorem and all subsequent results, can be found in Appendices <ref>, <ref> and <ref>.

There are polynomial filters (<ref>) that can be learned by MPNN but not by MPNN - VN. Conversely, any polynomial filter learned by MPNN (<ref>) can also be learned by MPNN + VN (<ref>). Furthermore, there exist polynomial filters that can be learned by MPNN + VN but not by MPNN.

Theorem <ref> implies that the ability of an MPNN - VN to replicate PairNorm comes at the cost of expressivity. Conversely, MPNN + VN (when they do not replicate PairNorm) are more expressive than MPNN when learning polynomial filters. In light of Theorem <ref>, we argue that, contrary to the claim in Cai et al., 2023, MPNN + VN generally avoids replicating PairNorm to maintain its greater expressivity. We empirically validate this statement in Section <ref> and report the smoothing effects of VN on common benchmarks in Appendix <ref>.

In summary, our analysis suggests that on graph-level tasks, benefits of MPNN + VN are independent of their ability to replicate methods such as PairNorm. In fact, for a graph-level prediction, whose output is computed using a global mean, some `smoothing' is actually unavoidable. In particular, the mean pooling in (<ref>) corresponds to projecting the final features precisely onto the subspace spanned by $\mathbf{1}\in \mathbb{R}^{n}$, i.e., the kernel of the (unnormalized) graph Laplacian associated with oversmoothing [Cai and Wang, 2020, Di Giovanni et al., 2023]. This allows us to clarify a key difference when studying oversmoothing for node-level and graph-level tasks. In the former case, if two node features become indistinguishable, then they will be assigned the same label, which is often independent of the task, thereby resulting in harmful oversmoothing. Conversely, in the latter case, even if the node features converge to the same representation, this behavior may result in beneficial oversmoothing if such common representation is controlled to match the desired, global output. We argue that for this reason, oversmoothing has hardly been observed as an issue on graph-level tasks. We refer to Proposition <ref> in Appendix <ref> for results relating VN to the choice of pooling and to ablations in Appendix <ref>.

§ MPNN + VN AND OVERSQUASHING

In Section <ref> we used polynomial filters to show that, on graph-level tasks, MPNN + VN loses expressive power when it replicates anti-oversmoothing methods like PairNorm. Nonetheless, VNs are often added to alleviate oversquashing, by reducing the diameter of the graph to 2. However, no formal analysis on the improvement brought by VN to oversquashing has been derived thus far. In this section, we fill this gap in literature and characterize how the spectrum of the input graph affects the impact of VN to commute time and hence the mixing abilities of the network. To study oversquashing, we analyze general realizations of (<ref>), whose layer updates have the form:

\begin{align}\label{eq:mpnn+vn_explicit}
h_{\mathsf{vn}}^{(\ell+1)} &= \sigma\Big(\boldsymbol{\Omega}_{\mathsf{vn}}^{(\ell)}h_{\mathsf{vn}}^{(\ell)} + \frac{1}{\tilde{n}}\sum_{j=1}^n\phi^{(\ell)}_{\mathsf{vn}}(h_{\mathsf{vn}}^{(\ell)},h_j^{(\ell)})\Big), \notag \\
h_i^{(\ell+1)}&= \sigma\Big(\boldsymbol{\Omega}^{(\ell)}h_i^{(\ell)} + \sum_{j=1}^n\mathsf{A}_{ij}\psi^{(\ell)}(h_i^{(\ell)},h_j^{(\ell)}) + \psi^{(\ell)}_{\mathsf{vn}}(h_i^{(\ell)}, h_{\mathsf{vn}}^{(\ell)})\Big),
\end{align}

where $\sigma$ is a pointwise nonlinear map, $\mathbf{\mathsf{A}}$ is a (potentially) normalized adjacency matrix, $\boldsymbol{\Omega}$ is a weight matrix, $\psi,\psi_{\mathsf{vn}},\phi_{\mathsf{vn}}$ are (learnable) message functions, and $\tilde{n}$ depends on our normalization choice (e.g. $\tilde{n} = 1$ for sum, or $\tilde{n} = n$ for mean). A key observation is that (<ref>) coincide with the MPNNs analyzed in Di Giovanni et al., 2023, but operating on the multi-relational graph $G_{\mathsf{vn}}$ with adjacency

\begin{equation}\label{eq:change_adjacency_VN}
\mathbf{A}_{\mathsf{vn}} = \begin{pmatrix} \mathbf{A} & \mathbf{1} \\ \mathbf{1}^\top & 0\end{pmatrix}, \quad \mathbf{1} = (1,\ldots,1)^\top \in \R^n,
\end{equation}

where edges connecting VN to nodes in $G$ have different type.

Di Giovanni et al., 2023, Black et al., 2023, Di Giovanni et al., 2023 have shown that oversquashing prevents the underlying model from exchanging information between nodes at large commute time $\tau$, where $\tau(i,j)$ measures the expected number of steps for a random walk to commute between $i$ and $j$. Accordingly, to assess if and how VN helps mitigate oversquashing, we need to determine whether the commute time $\tau_{\mathsf{vn}}$ of $G_{\mathsf{vn}}$ is smaller than the commute time $\tau$ of the original graph $G$. Below, we let $v_\ell$ be an orthonormal basis of eigenvectors for the graph Laplacian $\mathbf{L} = \mathbf{D} - \mathbf{A}$, with associated eigenvalues  $0 = \lambda_0 < \lambda_1 \leq \ldots \leq \lambda_{n-1}$.

The commute time between nodes $i,j$ after adding a VN changes as

\begin{equation}\label{eq:pair_change_main}
\tau_{\mathsf{vn}}(i,j) - \tau(i,j) = 2|E| \sum_{\ell = 1}^{n-1}\frac{1}{\lambda_\ell(\lambda_\ell +1)}\Big(\frac{n}{|E|}\lambda_\ell - 1\Big)(v_\ell(i) - v_\ell(j))^2.
\end{equation}

In particular, the average change in commute time is:

\begin{equation}\label{eq:averageCT}
\frac{1}{n^2}\sum_{i,j = 1}^{n}(\tau_{\mathsf{vn}}(i,j) - \tau(i,j)) = \frac{4|E|}{n}\sum_{\ell = 1}^{n-1}\frac{1}{\lambda_\ell(\lambda_\ell +1)}\Big(\frac{n}{|E|}\lambda_\ell - 1\Big).
\end{equation}

The result above highlights how a VN's impact can be determined in terms of the spectrum of the input graph. While there are cases, e.g. when the graph is complete, where (<ref>) is positive, for many real-world graphs adding VN reduces the overall commute time: We validate this claim in Section <ref>.

Di Giovanni et al., 2023 introduced the notion of mixing (see Appendix <ref>), to formally assess the amount of interactions between $i$ and $j$ required by the underlying task, and compared it to the mixing induced by the MPNN-prediction after $m$ layers.  In particular, they proved that oversquashing prevents MPNNs of bounded depth from solving tasks with required strong mixing between nodes at large commute time. Below, `bounded weights' means that weight matrices and message functions derivatives have bounded spectral norms—see Appendix <ref> for details.
There are graph-functions with mixing between nodes $i\neq j$ larger than some constant independent of $i,j$, such that for an MPNN of bounded weights to learn these functions, the  number of layers $m$ must satisfy $m\geq \tau(i,j)/8$.

For real-world graphs like peptides, small molecules, or images, (<ref>) is negative (see <ref>). In light of Theorem <ref>, this means that adding a VN should reduce the minimal number of layers required to learn functions with strong mixing between nodes. We make such connection explicit below:

Given $G$ and nodes $i,j$ for which (<ref>) is negative, there are graph-functions with mixing between $i\neq j$ larger than some constant independent of $i,j$, such that for an MPNN + VN of bounded weights to learn these functions, the minimal number of layers $m$ becomes

\begin{equation}
m\geq \frac{\tau_{\mathsf{vn}}(i,j)}{8} - \frac{|E|}{8}\cdot\sum_{\ell = 1}^{n-1}\frac{1}{\lambda_\ell(\lambda_\ell +1)}\Big(\frac{n}{|E|}\lambda_\ell - 1\Big)(v_\ell(i) - v_\ell(j))^2.
\end{equation}

The result in <Ref> shows that for real-world graphs where adding a VN significantly reduces the commute time—i.e., (<ref>) is negative—the minimal number of layers required by MPNN + VN to learn graph functions with strong mixing among $i,j$, is smaller than that of MPNN. This is the first result showing the extent to which MPNN + VN alleviates oversquashing by increasing the mixing abilities of the network. Crucially, MPNN + VN cannot do better than Theorem <ref>, meaning that their efficacy is a function of the spectrum.
While this implies that MPNN + VN may be sub-optimal, when compared to graph rewiring techniques such as GTs which can modify the commute time arbitrarily, the overall benefits of VN stem from their ability to reduce the commute time and hence mitigate oversquashing at very limited memory cost. We discuss next how we can further refine the formulation of VN to further close their gap with GTs, without paying any additional computational price.

§ COMPARING MPNN + VN AND GTS THROUGH SENSITIVITY ANALYSIS

In Section <ref> we showed that MPNN + VN can mitigate oversquashing and the extent to which this is possible. Another successful approach for combating oversquashing is given by Graph Transformers (GTs), that can entirely rewire the graph through the attention module. In fact, MPNN + VN are often used as a more efficient alternative to GTs [Cai et al., 2023]. In this section we further compare MPNN + VN and GTs through sensitivity analysis and show that the single layer of a GT can learn a heterogeneous node scoring, while MPNN + VN generally cannot. In light of this analysis, we propose a simple variation of MPNN + VN, called $\acro$, which better uses the graph to learn a heterogeneous measure of node relevance while retaining the same efficiency as MPNN + VN.

MPNN + VN Layers Are Homogeneous We start by reporting the layer update of the GPS architecture [Rampášek et al., 2022], one instance of the GT class that encodes both local and global information:

\begin{equation}\label{eq:GPS}
h_i^{(\ell+1)} = f^{(\ell)}(h_{i,{\rm loc}}^{(\ell+1)}
+\mathbf{Q}^{(\ell)}\sum_{k=1}^n a(h_i^{(\ell)}, h_k^{(\ell)})h_k^{(\ell)}),
\end{equation}

where $f$ is an MLP, $h_{\rm loc}$ is the local update given by the choice of MPNN, and $a$ is the attention module. Depending on data and positional encoding, GTs can capture heterogeneous relations. We can quantify such heterogeneity via sensitivity analysis by deriving that the Jacobian $\partial h_i^{(\ell+1)}/\partial h_k^{(\ell)}$ is, in general, a function that depends on $k$, meaning that the state of node $i$ at layer $\ell + 1$ is affected by the state of a node $k$ at the previous layer, as a function varying with $k$. Conversely, MPNN + VN depends on different nodes more uniformly. Explicitly, consider an MPNN + VN as in (<ref>), where the VN update is

\begin{align}
h_{\mathsf{vn}}^{(\ell+1)} &= \sigma(\boldsymbol{\Omega}^{(\ell)}_{\mathsf{vn}}h_{\mathsf{vn}}^{(\ell)} + \mathbf{W}_{\mathsf{vn}}^{(\ell)}\mathsf{Mean}(\{h_{j}^{(\ell)}\})). \label{eq:mpnn_vn_jacobian}
\end{align}

Given a node $k$ at a distance larger than $2$ from node $i$, any message from node $k$ at layer $\ell-1$ is first received by node $i$ at layer $\ell+1$ through the VN.

For MPNN + VN whose VN update is (<ref>), the Jacobian $\partial h_i^{(\ell+1)}/\partial h_k^{(\ell-1)}$ is independent of $k$ whenever $k$ and $i$ are separated by more than 2 hops.

We see that node $i$ receives a global yet homogeneous update from VN, where
each node $k$ at distance $>2$ contributes the same to its state (after two layers). This in stark contrast to the case of GTs.

$\acro$: A New Formulation Of VN Inspired by frameworks such as (<ref>),
to allow for more heterogeneous sensitivity through the VN, we propose $\acro$, a simple variation to MPNN + VN with the same computational complexity, whose layer updates have the form:

\begin{align}\label{eq:mpnn+VN_G}
h^{(\ell+1)}_{i,{\rm loc}} &= \mathsf{up}^{(\ell)}( h^{(\ell)}_i, \mathsf{agg}^{(\ell)} ( \{h^{(\ell)}_j : j \in N_i\})),  \notag \\
h^{(\ell+1)}_{\mathsf{vn}} &= \mathsf{up}_{\mathsf{vn}}^{(\ell)}( h^{(\ell)}_{\mathsf{vn}}, \mathsf{agg}_{\mathsf{vn}}^{(\ell)} ( \{h^{(\ell+1)}_{j,{\rm loc}} : j \in V\})), \notag \\
h^{(\ell+1)}_i &= \widetilde{\mathsf{up}}^{(\ell)} ( h^{(\ell+1)}_{i,{\rm loc}},  h^{(\ell+1)}_{\mathsf{vn}}).
\end{align}

Differently from MPNN + VN (<ref>), we first compute local updates based on the choice of aggregation (i.e. the underlying MPNN model), then we compute a global update through VN using such local representations (as illustrated in Appendix <ref>), and finally we combine the local and non-local representations. To justify how $\acro$ can assign more heterogeneous values to the global updates associated with the VN, we study the sensitivity of node features. Similarly to (<ref>) and (<ref>), we consider an instance of $\acro$ of the form:

\begin{equation}\label{eq:mpnn+vn_G_jacobian}
h_i^{(\ell+1)} = h_{i,{\rm loc}}^{(\ell+1)} + \mathsf{Mean}(\{h_{j,{\rm loc}}^{(\ell+1)}\}), \,\,\, h_{j,{\rm loc}}^{(\ell+1)} = \sigma(\boldsymbol{\Omega}^{(\ell)}h_i^{(\ell)} + \sum_{j\in N_i} \psi_{ij}^{(\ell)}(h_i^{(\ell)},h_j^{(\ell)})).
\end{equation}

Differently from (<ref>), two nodes now can exchange information after a single layer. Below, we let $z_i^{(\ell)}$ be the argument of $\sigma$ in (<ref>), for each $i\in V$, and $\nabla_s$ be the Jacobian with respect to variable $s$.

Given $i\in V$ and $k\in V\setminus N_i$, the Jacobian $\partial h_i^{(\ell+1)}/\partial h_k^{(\ell)}$ computed using (<ref>) is

\begin{equation*}
\frac{\partial h_i^{(\ell+1)}}{\partial h_k^{(\ell)}} = \frac{1}{n}\Big(\hspace{-1.5pt}{\rm diag}(\sigma'(z_k^{(\ell)}))(\boldsymbol{\Omega}^{(\ell)} + \hspace{-3pt}\sum_{u\in N_k}\hspace{-4pt}\nabla_1 \psi^{(\ell)}_{ku}(h_k^{(\ell)},h_u^{(\ell)})) + \hspace{-3pt}\sum_{u\in N_k}\hspace{-4pt}{\rm diag}(\sigma'(z_u^{(\ell)}))\nabla_2\psi_{uk}^{(\ell)}(h_u^{(\ell)},h_k^{(\ell)})\hspace{-1.5pt} \Big).
\end{equation*}

We see that the sensitivity of node $i$ to the feature of $k$ is a function of the graph-topology, i.e., the message functions evaluated over the neighbors of $k$. In particular, our proposed variation $\acro$ uses the same message functions $\psi_{ij}$ to learn a heterogeneous global update. Related to recent advances in graph pooling [Ranjan et al., 2020, Lee et al., 2019, Grattarola et al., 2022, Bianchi and Lachi, 2023], we utilize the graph structure in the global aggregator; nodes that are more relevant to each of their neighbors, are more likely to contribute more to the global update, instead of weighting each node the same as for MPNN + VN. We note that advantages of this approach depend on the choice of the underlying
message function $\psi$ in (<ref>) (see
Table <ref> in Appendix <ref>).

In light of our sensitivity analysis, we argue that for those tasks where an empirical gap between MPNN + VN and GT can be found, this is due to the latter leveraging its ability to assign heterogeneous node scores in its global update. Conversely, when no significant gap arises, we claim this is due to the task not requiring heterogeneous node scores. We validate these points in Section <ref>.

§ EXPERIMENTS

The purpose of this section is to first empirically verify the theoretical claims made in Sections <ref>, <ref> and <ref>. Then, we evaluate our proposed MPNN + VN$_G$ on a diverse set of benchmarks (<ref>).

§.§ Ablations

Our ablations aim to validate the previous theoretical results, by answering the following questions:

(Q1) To what extent does MPNN + VN replicate PairNorm on graph-level tasks? [<ref>]

(Q2) How does adding a VN affect the average commute time on real-world graphs? [<ref>]

(Q3) Which tasks show a gap between MPNN + VN and GT due to node heterogeneity? [<ref>]

The effect of PairNorm and MPNN + VN on LRGB tasks.

Method     Peptides-Func ($\uparrow$)     Peptides-Struct ($\downarrow$)

GCN                 0.5930 $\pm 0.0023$         0.3496 $\pm 0.0013$

GCN + PairNorm                 0.4980  $\pm 0.0031$           0.3471   $\pm 0.0016$

GCN + VN                 0.6623  $\pm 0.0038$           0.2488   $\pm 0.0021$

GINE                 0.5498 $\pm 0.0079$             0.3547 $\pm 0.0045$

GINE + PairNorm                 0.4698 $\pm 0.0053$            0.3562 $\pm 0.0007$

GINE + VN                 0.6346  $\pm 0.0071$           0.2584   $\pm 0.0011$

GatedGCN                 0.5864 $\pm 0.0077$            0.3420 $\pm 0.0013$

GatedGCN  + PairNorm                 0.4674 $\pm 0.0040$            0.3551 $\pm 0.0008$

GatedGCN + VN                 0.6477  $\pm 0.0039$           0.2523   $\pm 0.0016$

Effect of adding a virtual node on the average commute time of four graph datasets.

To answer (Q1) and assess whether MPNN + VN behaves differently to PairNorm in practice, we evaluate both on two graph-level tasks from the Long-Range Benchmarks [Dwivedi et al., 2022] using three different MPNN architectures. It is clear from Table <ref> that whilst PairNorm has a damaging impact on the results as hypothesized in <ref>, MPNN + VN achieves significant gains over the standard MPNN architecture. Further datasets and results using a GatedGCN can be found in Table <ref>. One explanation for the observed phenomenon is that, on these datasets, MPNN + VN favors alignment between the layerwise representations and the final output representation (a smoothing process). This phenomenon can be observed when comparing the cosine similarity of the pooled representations at each layer to the final pooled representation. Plots of this cosine similarity over the layers can be seen in Appendix <ref> and clearly show that the VN causes a smoothing effect in contrast to PairNorm.

For (Q2) we calculated Equation (<ref>) on 100 graphs from four commonly used real-world datasets covering peptides, small molecules and images. The distribution of the change of commute time across these four datasets is shown in Figure <ref> and highlights the tendency for this to be negative. This demonstrates that adding a VN is highly beneficial on real-world graphs since it significantly reduces their commute time, thereby increasing their mixing abilities as per Theorem <ref>.

Here, we aim to provide an empirical justification for the similar performances of MPNN and GT [Cai et al., 2023] and answer (Q3). To measure the heterogeneity of node relevance in the global update, we look at the similarity between the rows of the attention matrix in a GT. In fact, similar rows in $a(h_i^{(\ell)}, h_k^{(\ell)})$ in (<ref>) correspond to assigning similar weights to each node representation, and hence less heterogeneity. To measure this, we calculate the mean standard deviation (std) of each column. We calculated this std on the first layer attention matrix, after training GPS with a GatedGCN as the base MPNN, on four real-world benchmarks from [Dwivedi et al., 2022] and [Dwivedi et al., 2023]. In Table <ref>, we can see that for datasets where the std is small (Peptides-func, Peptides-struct, MNIST), using an MPNN + VN can match or even improve over the GPS implementation. On the other hand, CIFAR10 has a larger std, and we see that GPS outperforms the VN architecture. These results show that the performance gap between MPNN + VN and GTs is related to whether node heterogeneity is used by the GT on the benchmark, considering that MPNN + VN assigns equal sensitivity to different node features (Proposition <ref>). To further show the differences in the attention matrices and the amount of homogeneity on these tasks, we have visualized them in Appendix <ref>. Additionally, we can measure the importance of heterogeneity for the task by removing the heterogeneity in the GT and seeing if we lose performance. We can do this in the GPS framework by setting all the rows of the attention matrix to be equal to their mean. We call this GPS + projection, and it forces the GPS model to become completely homogeneous. Again, we see that on datasets where the homogeneous GatedGCN+VN does well and the std of the attention matrix columns is low, we do not lose much performance with this mean projection. Whereas on CIFAR10, we find that the heterogeneity of nodes is important.

Effects of projecting the non-local part of GPS onto the mean and its comparison to using a VN. We also report the standard deviation of the column sums in the first attention layer.

Method     Pept-Func ($\uparrow$)     Pept-Struct ($\downarrow$)      MNIST ($\uparrow$)     CIFAR10 ($\uparrow$)

GPS                 0.6534 $\pm .0091$         0.2509 $\pm .0014$        98.051 $\pm .126$         72.298 $\pm .356$

GPS + projection               0.6498 $\pm .0054$            0.2487 $\pm .0011$       98.176  $\pm .120$         71.455  $\pm .513$

GatedGCN+PE+VN      0.6712 $\pm .0066$           0.2481 $\pm .0015$      98.122  $\pm .102$         70.280  $\pm .380$

std attention layer     0.0011     0.0007     0.0006     0.0038

§.§ Evaluating Our Proposed $\acro$

We first evaluate $\acro$ on a diverse set of 5 graph-level datasets, comparing against classical MPNN benchmarks, Graph-Transformers [Ma et al., 2023, Shirzad et al., 2023, Rampášek et al., 2022, Kreuzer et al., 2021, Hussain et al., 2022], graph rewiring [Gutteridge et al., 2023], a random-walk based method [Tönshoff et al., 2023], and a generalization of ViT/MLP-Mixer [He et al., 2023]. Crucially, we also compare VN$_G$ with the standard implementation of VN [Cai et al., 2023].

Test performance on two LRGB datasets [Dwivedi et al., 2022] and three other benchmarks from [Dwivedi et al., 2023]. For Peptides-Func and Peptides-Struct, $\pm$ std is shown over 4 runs whilst the remaining datasets are over 10 runs (missing values from literature are indicated by `-'). The first, second and third best results for each task are color-coded.

Method     Pept-Func ($\uparrow$)     Pept-Struct ($\downarrow$)      MNIST ($\uparrow$)     CIFAR10 ($\uparrow$)     MalNet-Tiny ($\uparrow$)

GCN                 0.5930 $\pm 0.0023$         0.3496 $\pm 0.0013$        90.705 $\pm 0.218$         55.710 $\pm 0.381$     81.0

GINE                 0.5498 $\pm 0.0079$         0.3547 $\pm 0.0045$        96.485 $\pm 0.252$         55.255 $\pm 1.527$     88.98 $\pm 0.56$

GatedGCN                 0.5864 $\pm 0.0077$         0.3420 $\pm 0.0013$        97.340 $\pm 0.143$         67.312 $\pm 0.311$     92.23 $\pm 0.65$

GatedGCN+PE                 0.6765 $\pm 0.0047$         0.2477 $\pm 0.0009$        -         69.948 $\pm 0.499$     -

GatedGCN+PE-ViT                 0.6942 $\pm 0.0075$         0.2465 $\pm 0.0015$        98.460 $\pm 0.090$         71.580 $\pm 0.090$     -

GatedGCN+PE-Mixer                 0.6932 $\pm 0.0017$         0.2508 $\pm 0.0007$        98.320 $\pm 0.040$         70.600 $\pm 0.220$     -

CRaWl                 0.7074 $\pm 0.0032$         0.2506 $\pm 0.0022$        97.940 $\pm 0.050$         69.010 $\pm 0.259$     -

DRew                 0.7150 $\pm 0.0044$         0.2536 $\pm 0.0015$        -         -     -

SAN+RWSE                 0.6439 $\pm 0.0075$         0.2545 $\pm 0.0012$        -         -     -

EGT                  -         -        98.173 $\pm 0.087$         68.702 $\pm 0.409$     -

GRIT                 0.6988 $\pm 0.0082$         0.2460 $\pm 0.0012$        98.108 $\pm 0.111$         76.468 $\pm 0.881$     -

GPS                 0.6534 $\pm 0.0091$         0.2509 $\pm 0.0014$        98.051 $\pm 0.126$         72.298 $\pm 0.356$     93.50 $\pm 0.41$

Exphormer               0.6527 $\pm 0.0043$            0.2481 $\pm 0.0007$         98.414 $\pm 0.035$        74.690 $\pm 0.125$       94.02 $\pm 0.21$

GatedGCN+PE+VN              0.6712 $\pm 0.0066$            0.2481 $\pm 0.0015$       98.122 $\pm 0.102$         70.280 $\pm 0.380$     92.62 $\pm 0.57$

GatedGCN+PE+$\text{VN}_{G}$                0.6822 $\pm 0.0052$            0.2458  $\pm 0.0006$      98.626 $\pm 0.100$         76.080 $\pm 0.330$     93.67 $\pm 0.37$

Experimental Details. We evaluated $\acro$ on the Long-Range Graph Benchmark using a fixed 500k parameter budget and averaging over four runs. These molecular datasets (Peptides-Func,
Peptides-Struct) have been proposed to test a method’s ability to capture long-range dependencies. Additionally, we used two graph-level image-based datasets from Benchmarking GNNs [Dwivedi et al., 2023], where we run our model over 10 seeds. We also used a code dataset MalNet-Tiny [Freitas et al., 2020]. This dataset consists of function call graphs with up to 5,000 nodes. These graphs are considerably larger than previously considered graph-level benchmarks and showcase our model's ability to improve over MPNN baselines on a large dataset. We then evaluated our approach on three graph-level tasks from the Open Graph Benchmark [Hu et al., 2020], namely molhiv, molpcba and ppa. Details on datasets and training parameters used, can be found in Appendix <ref>.

We can see from Tables <ref> and <ref> that our method performs well across a variety of datasets and achieves the highest performance on Peptides-Struct, MNIST, ogbg-molhiv and ogbg-ppa. Furthermore, $\acro$ improves over MPNN + VN on all tasks, with the largest percentage improvement being shown on CIFAR10 (8.25%). The latter is perfectly aligned with the discussion in Section <ref> and <ref>, where we have described the importance of heterogeneity to the VN and justified how heterogeneity is particularly required on CIFAR10. Additionally, on the peptides datasets, which have small degrees but large diameters, we find that using $\acro$ can perform better than a variety of transformer-based approaches. Not only does this suggest that MPNN + $\hetvn$  works well in practice but that complex long-range dependencies might not play a primary role for these tasks. This was also supported by Tönshoff et al., 2023 where the lower performance of MPNNs compared to GTs was found to be caused by insufficient hyperparameter tuning. We also show strong performance on the larger MalNet-Tiny dataset. Previous work has shown that using a full transformer outperforms linear attention-based methods [Rampášek et al., 2022] on this task. However, we show that we can match the performance of the full transformer and still be more computationally efficient.

Test performance on graph-level OGB benchmarks [Hu et al., 2020]. Shown is the mean % $\pm$ std of 10 runs (missing values from literature are
indicated with `-'). The first, second and third best results for each task are color-coded.

Method     MolHIV ($\uparrow$)      MolPCBA ($\uparrow$)      PPA ($\uparrow$)

GCN                 76.06 $\pm 0.97$       20.20 $\pm 0.24$        68.39 $\pm 0.84$

GIN                 75.58 $\pm 1.40$       22.66 $\pm 0.28$        68.92 $\pm 1.00$

DeeperGCN                 78.58 $\pm 1.17$       27.81 $\pm 0.38$        77.12 $\pm 0.71$

CRaWL                 -       29.86 $\pm 0.25$        -

SAN                 77.85 $\pm 0.28$       27.65 $\pm 0.42$        -

GPS                 78.80 $\pm 1.01$        29.07 $\pm 0.28$        80.15 $\pm 0.33$

GCN+VN     75.99 $\pm 1.19$       24.24 $\pm 0.34$        68.57 $\pm 0.61$

GIN+VN     77.07 $\pm 1.49$       27.03 $\pm 0.23$        70.37 $\pm 1.07$
GatedGCN+PE+VN              76.87 $\pm
1.36$            28.48 $\pm 0.26$        80.27 $\pm 0.26$

GatedGCN+PE+$\text{VN}_{G}$              79.10  $\pm 0.86$           29.17  $\pm
0.27$        81.17 $\pm 0.30$

§ CONCLUSION

In <ref> we have highlighted how MPNN + VN may not emulate anti-smoothing approaches such as PairNorm to maintain expressive power and have shown in <ref> that this holds in practice. In <ref>, we characterize how the spectrum of the input graph affects the impact of VN on the commute time and hence to oversquashing. We have also evaluated this on real-world benchmarks in <ref>. In <ref>, we show that standard VN implementations cannot assign different scores to nodes in the global update. Consequently, we propose $\acro$, a variation of MPNN + VN sharing the same computational complexity, that uses the graph to learn heterogeneous node importance. We show that this outperforms MPNN + VN on a range of benchmarks, precisely corresponding to those where GTs outperform MPNN + VN by learning heterogenos node scores.

Limitations and Future Work. Our expressivity analysis in <ref> relies on layers which are linear. Although our experiments in <ref> show that it holds in practical, non-linear scenarios, understanding VN and oversmoothing in a non-linear setting will be interesting future work. Additionally, we have not empirically explored the change in commute time afforded by GT and rewiring-based approaches. Future work can compare these approaches and VN, to assess their mixing abilities.

[Alon and Yahav, 2021]
U. Alon and E. Yahav.
On the bottleneck of graph neural networks and its practical
In International Conference on Learning Representations, 2021.
URL <https://openreview.net/forum?id=i80OPhOCVH2>.

[Arnaiz-Rodríguez et al., 2022]
A. Arnaiz-Rodríguez, A. Begga, F. Escolano, and N. M. Oliver.
Diffwire: Inductive graph rewiring via the lovász bound.
In The First Learning on Graphs Conference, 2022.
URL <https://openreview.net/forum?id=IXvfIex0mX6f>.

[Barbero et al., 2023]
F. Barbero, A. Velingker, A. Saberi, M. Bronstein, and F. Di Giovanni.
Locality-aware graph-rewiring in gnns.
arXiv preprint arXiv:2310.01668, 2023.

[Battaglia et al., 2018]
P. W. Battaglia, J. B. Hamrick, V. Bapst, A. Sanchez-Gonzalez, V. Zambaldi,
M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, et al.
Relational inductive biases, deep learning, and graph networks.
arXiv preprint arXiv:1806.01261, 2018.

[Bianchi and Lachi, 2023]
F. M. Bianchi and V. Lachi.
The expressive power of pooling in graph neural networks.
In Thirty-seventh Conference on Neural Information Processing
Systems, 2023.

[Black et al., 2023]
M. Black, Z. Wan, A. Nayyeri, and Y. Wang.
Understanding oversquashing in gnns through the lens of effective
In International Conference on Machine Learning, pages
2528–2547. PMLR, 2023.

[Bouritsas et al., 2023]
G. Bouritsas, F. Frasca, S. Zafeiriou, and M. M. Bronstein.
Improving graph neural network expressivity via subgraph isomorphism
IEEE Transactions on Pattern Analysis and Machine
Intelligence, 450 (1):0 657–668, Jan. 2023.
ISSN 1939-3539.
URL <http://dx.doi.org/10.1109/TPAMI.2022.3154319>.

[Cai and Wang, 2020]
C. Cai and Y. Wang.
A note on over-smoothing for graph neural networks.
arXiv preprint arXiv:2006.13318, 2020.

[Cai et al., 2023]
C. Cai, T. S. Hy, R. Yu, and Y. Wang.
On the connection between mpnn and graph transformer.
arXiv preprint arXiv:2301.11956, 2023.

[Chen et al., 2022]
D. Chen, L. O’Bray, and K. Borgwardt.
Structure-aware transformer for graph representation learning.
In International Conference on Machine Learning, pages
3469–3489. PMLR, 2022.

[Chen et al., 2020]
Z. Chen, L. Chen, S. Villar, and J. Bruna.
Can graph neural networks count substructures?
Advances in neural information processing systems,
33:0 10383–10395, 2020.

[Dehmamy et al., 2019]
N. Dehmamy, A.-L. Barabási, and R. Yu.
Understanding the representation power of graph neural networks in
learning graph topology.
Advances in Neural Information Processing Systems, 32, 2019.

[Di Giovanni et al., 2023]
F. Di Giovanni, L. Giusti, F. Barbero, G. Luise, P. Lio, and M. M. Bronstein.
On over-squashing in message passing neural networks: The impact of
width, depth, and topology.
In International Conference on Machine Learning, pages
7865–7885. PMLR, 2023a.

[Di Giovanni et al., 2023]
F. Di Giovanni, J. Rowbottom, B. P. Chamberlain, T. Markovich, and M. M.
Understanding convolution on graphs via energies.
Transactions on Machine Learning Research, 2023b.

[Di Giovanni et al., 2023]
F. Di Giovanni, T. K. Rusch, M. M. Bronstein, A. Deac, M. Lackenby, S. Mishra,
and P. Veličković.
How does over-squashing affect the power of gnns?
arXiv preprint arXiv:2306.03589, 2023c.

[Dwivedi et al., 2021]
V. P. Dwivedi, A. T. Luu, T. Laurent, Y. Bengio, and X. Bresson.
Graph neural networks with learnable structural and positional
In International Conference on Learning Representations, 2021.

[Dwivedi et al., 2022]
V. P. Dwivedi, L. Rampášek, M. Galkin, A. Parviz, G. Wolf, A. T. Luu,
and D. Beaini.
Long range graph benchmark.
Advances in Neural Information Processing Systems,
35:0 22326–22340, 2022.

[Dwivedi et al., 2023]
V. P. Dwivedi, C. K. Joshi, A. T. Luu, T. Laurent, Y. Bengio, and X. Bresson.
Benchmarking graph neural networks.
Journal of Machine Learning Research, 240
(43):0 1–48, 2023.

[Freitas et al., 2020]
S. Freitas, Y. Dong, J. Neil, and D. H. Chau.
A large-scale database for graph representation learning.
arXiv preprint arXiv:2011.07682, 2020.

[Ghosh et al., 2008]
A. Ghosh, S. Boyd, and A. Saberi.
Minimizing effective resistance of a graph.
SIAM review, 500 (1):0 37–66, 2008.

[Gilmer et al., 2017]
J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl.
Neural message passing for quantum chemistry.
In International conference on machine learning, pages
1263–1272. PMLR, 2017.

[Gori et al., 2005]
M. Gori, G. Monfardini, and F. Scarselli.
A new model for learning in graph domains.
In Proceedings. 2005 IEEE International Joint Conference on
Neural Networks, 2005., volume 2, pages 729–734 vol. 2, 2005.

[Grattarola et al., 2022]
D. Grattarola, D. Zambon, F. M. Bianchi, and C. Alippi.
Understanding pooling in graph neural networks.
IEEE Transactions on Neural Networks and Learning Systems,
page 1–11, 2022.
ISSN 2162-2388.
URL <http://dx.doi.org/10.1109/TNNLS.2022.3190922>.

[Gutteridge et al., 2023]
B. Gutteridge, X. Dong, M. M. Bronstein, and F. Di Giovanni.
Drew: Dynamically rewired message passing with delay.
In International Conference on Machine Learning, pages
12252–12267. PMLR, 2023.

[He et al., 2023]
X. He, B. Hooi, T. Laurent, A. Perold, Y. LeCun, and X. Bresson.
A generalization of vit/mlp-mixer to graphs.
In International Conference on Machine Learning, pages
12724–12745. PMLR, 2023.

[Horn et al., 2022]
M. Horn, E. D. Brouwer, M. Moor, Y. Moreau, B. Rieck, and K. Borgwardt.
Topological graph neural networks.
In International Conference on Learning Representations, 2022.
URL <https://openreview.net/forum?id=oxxUMeFwEHd>.

[Hu et al., 2020]
W. Hu, M. Fey, M. Zitnik, Y. Dong, H. Ren, B. Liu, M. Catasta, and J. Leskovec.
Open graph benchmark: Datasets for machine learning on graphs.
Advances in neural information processing systems,
33:0 22118–22133, 2020.

[Hussain et al., 2022]
M. S. Hussain, M. J. Zaki, and D. Subramanian.
Global self-attention as a replacement for graph convolution.
In Proceedings of the 28th ACM SIGKDD Conference on Knowledge
Discovery and Data Mining, KDD ’22. ACM, Aug. 2022.
URL <http://dx.doi.org/10.1145/3534678.3539296>.

[Hwang et al., 2022]
E. Hwang, V. Thost, S. S. Dasgupta, and T. Ma.
An analysis of virtual nodes in graph neural networks for link
In The First Learning on Graphs Conference, 2022.

[Kingma and Ba, 2014]
D. P. Kingma and J. Ba.
Adam: A method for stochastic optimization.
arXiv preprint arXiv:1412.6980, 2014.

[Kreuzer et al., 2021]
D. Kreuzer, D. Beaini, W. Hamilton, V. Létourneau, and P. Tossou.
Rethinking graph transformers with spectral attention.
Advances in Neural Information Processing Systems,
34:0 21618–21629, 2021.

[Lee et al., 2019]
J. Lee, I. Lee, and J. Kang.
Self-attention graph pooling.
In International conference on machine learning, pages
3734–3743. PMLR, 2019.

[Loukas, 2019]
A. Loukas.
What graph neural networks cannot learn: depth vs width.
In International Conference on Learning Representations, 2019.

[Lovász, 1993]
L. Lovász.
Random walks on graphs.
Combinatorics, Paul erdos is eighty, 20
(1-46):0 4, 1993.

[Ma et al., 2023]
L. Ma, C. Lin, D. Lim, A. Romero-Soriano, P. K. Dokania, M. Coates, P. Torr,
and S.-N. Lim.
Graph inductive biases in transformers without message passing.
arXiv preprint arXiv:2305.17589, 2023.

[Maron et al., 2018]
H. Maron, H. Ben-Hamu, N. Shamir, and Y. Lipman.
Invariant and equivariant graph networks.
In International Conference on Learning Representations, 2018.

[Micheli, 2009]
A. Micheli.
Neural network for graphs: A contextual constructive approach.
IEEE Transactions on Neural Networks, 200
(3):0 498–511, 2009.

[Morris et al., 2019]
C. Morris, M. Ritzert, M. Fey, W. L. Hamilton, J. E. Lenssen, G. Rattan, and
M. Grohe.
Weisfeiler and leman go neural: Higher-order graph neural networks.
In Proceedings of the AAAI conference on artificial
intelligence, volume 33, pages 4602–4609, 2019.

[Nikolentzos et al., 2020]
G. Nikolentzos, G. Dasoulas, and M. Vazirgiannis.
K-hop graph neural networks.
Neural Networks, 130:0 195–205, 2020.

[Nt and Maehara, 2019]
H. Nt and T. Maehara.
Revisiting graph neural networks: All we have is low-pass filters.
arXiv preprint arXiv:1905.09550, 2019.

[Oono and Suzuki, 2020]
K. Oono and T. Suzuki.
Graph neural networks exponentially lose expressive power for node
In International Conference on Learning Representations, 2020.
URL <https://openreview.net/forum?id=S1ldO2EFPr>.

[Pham et al., 2017]
T. Pham, T. Tran, H. Dam, and S. Venkatesh.
Graph classification via deep learning with virtual nodes.
arXiv preprint arXiv:1708.04357, 2017.

[Rampášek et al., 2022]
L. Rampášek, M. Galkin, V. P. Dwivedi, A. T. Luu, G. Wolf, and
D. Beaini.
Recipe for a general, powerful, scalable graph transformer.
Advances in Neural Information Processing Systems,
35:0 14501–14515, 2022.

[Ranjan et al., 2020]
E. Ranjan, S. Sanyal, and P. Talukdar.
Asap: Adaptive structure aware pooling for learning hierarchical
graph representations.
In Proceedings of the AAAI Conference on Artificial
Intelligence, volume 34, pages 5470–5477, 2020.

[Rosenbluth et al., 2024]
E. Rosenbluth, J. Tönshoff, M. Ritzert, B. Kisin, and M. Grohe.
Distinguished in uniform: Self-attention vs. virtual nodes.
In The Twelfth International Conference on Learning
Representations, 2024.
URL <https://openreview.net/forum?id=AcSChDWL6V>.

[Scarselli et al., 2009]
F. Scarselli, M. Gori, A. C. Tsoi, M. Hagenbuchner, and G. Monfardini.
The graph neural network model.
IEEE Transactions on Neural Networks, 200
(1):0 61–80, 2009.

[Segol and Lipman, 2019]
N. Segol and Y. Lipman.
On universal equivariant set networks.
arXiv preprint arXiv:1910.02421, 2019.

[Sestak et al., 2024]
F. Sestak, L. Schneckenreiter, J. Brandstetter, S. Hochreiter, A. Mayr, and
G. Klambauer.
Vn-egnn: E (3)-equivariant graph neural networks with virtual nodes
enhance protein binding site identification.
arXiv preprint arXiv:2404.07194, 2024.

[Shirzad et al., 2023]
H. Shirzad, A. Velingker, B. Venkatachalam, D. J. Sutherland, and A. K. Sinop.
Exphormer: Sparse transformers for graphs.
arXiv preprint arXiv:2303.06147, 2023.

[Tönshoff et al., 2023]
J. Tönshoff, M. Ritzert, E. Rosenbluth, and M. Grohe.
Where did the gap go? reassessing the long-range graph benchmark.
arXiv preprint arXiv:2309.00367, 2023a.

[Tönshoff et al., 2023]
J. Tönshoff, M. Ritzert, H. Wolf, and M. Grohe.
Walking out of the weisfeiler leman hierarchy: Graph learning beyond
message passing.
Transactions on Machine Learning Research, 2023b.

[Topping et al., 2021]
J. Topping, F. Di Giovanni, B. P. Chamberlain, X. Dong, and M. M. Bronstein.
Understanding over-squashing and bottlenecks on graphs via curvature.
arXiv preprint arXiv:2111.14522, 2021.

[Wang et al., 2020]
G. Wang, R. Ying, J. Huang, and J. Leskovec.
Multi-hop attention graph neural network.
arXiv preprint arXiv:2009.14332, 2020.

[Wang and Zhang, 2022]
X. Wang and M. Zhang.
How powerful are spectral graph neural networks.
In International Conference on Machine Learning, pages
23341–23362. PMLR, 2022.

[Weisfeiler and Leman, 1968]
B. Weisfeiler and A. Leman.
The reduction of a graph to canonical form and the algebra which
appears therein.
nti, Series, 20 (9):0 12–16, 1968.

[Wu et al., 2021]
Z. Wu, P. Jain, M. Wright, A. Mirhoseini, J. E. Gonzalez, and I. Stoica.
Representing long-range context for graph neural networks with global
Advances in Neural Information Processing Systems,
34:0 13266–13279, 2021.

[Xu et al., 2018]
K. Xu, W. Hu, J. Leskovec, and S. Jegelka.
How powerful are graph neural networks?
In International Conference on Learning Representations, 2018.

[Ying et al., 2021]
C. Ying, T. Cai, S. Luo, S. Zheng, G. Ke, D. He, Y. Shen, and T.-Y. Liu.
Do transformers really perform badly for graph representation?
Advances in Neural Information Processing Systems,
34:0 28877–28888, 2021.

[Zhao and Akoglu, 2019]
L. Zhao and L. Akoglu.
Pairnorm: Tackling oversmoothing in gnns.
arXiv preprint arXiv:1909.12223, 2019.

§ OUTLINE OF THE APPENDIX

We report proofs and additional details pertaining to <Ref> in <Ref>. In particular, we derive further theoretical results connecting VN and the choice of global pooling.
Next, we report the analysis on commute time and oversquashing when adding a VN in <Ref>. Proofs for the sensitivity analysis results in <Ref> are finally reported in <Ref>.

Regarding additional ablations and experimental details, we comment on training details, hyperparameters, datasets, and hardware in <Ref>. Next, we describe additional experiments confirming that VN introduces smoothing on graph-level tasks in <Ref>, focusing on the impact of the final pooling. We further validate the role played by the underlying MPNN model in MPNN + VN$_G$ in <Ref>, propose a visualization of the framework in <Ref>, and finally report attention maps learned by Graph-Transformer on benchmarks—used as a proxy-measure for the heterogeneity required by the task—in <Ref>.

§ PROOFS OF RESULTS IN SECTION <REF>

First, we note that similarly to Cai et al., 2023, the VN state is initialized to be zero, i.e., the first layer of <ref> is

\begin{equation}
\mathbf{H}^{(1)} = \mathbf{H} + \mathbf{A}\mathbf{H}\mathbf{W}^{(0)}
\end{equation}

where $\mathbf{H}$ is matrix of input node features.
It is trivial to check that any polynomial filter that can be learned by an MPNN can also be learned by MPNN + VN, considering that if we set the additional weight matrices $\{\mathbf{Q}_\ell\}$ to be zero in (<ref>), then we recover the MPNN case in (<ref>). Accordingly, to prove that MPNN + VN is strictly more expressive than MPNN when learning polynomial filters, we only need to show that there exist $m\in\mathbb{N}$ and polynomial filters of degree $m$ that can be learned by an MPNN + VN of $m$ layers, but not by MPNN. We show that such a polynomial can already be found when $m=2$. Namely, given a linear 2-layer MPNN + VN, we can write the 2-layer representation $\mathbf{H}^{(2)}$ explicitly as

\begin{align}
\mathbf{H}^{(2)} &= \mathbf{H} + \mathbf{A}\mathbf{H}(\mathbf{W}^{(0)} + \mathbf{W}^{(1)}) + \mathbf{A}^2\mathbf{H}\mathbf{W}^{(0)}\mathbf{W}^{(1)} + \frac{1}{n}\mathbf{1}\mathbf{1}^\top\mathbf{H}\mathbf{Q}^{(1)}. \label{eq:app_2_layer_expansion}
% + \mathbf{Q}^{(1)} + \mathbf{Q}^{(0)}\mathbf{Q}^{(1)}) + \frac{1}{n}\mathbf{1}\mathbf{1}^\top\mathbf{A}\mathbf{H}\mathbf{W}^{(0)}\mathbf{Q}^{(1)} + \frac{1}{n}\mathbf{A}\mathbf{1}\mathbf{1}^\top\mathbf{H}\mathbf{Q}^{(0)}\mathbf{W}^{(1)}.
\end{align}

Recall that given the choice of global mean pooling, the final output is

\[
\mathbf{Y}_2 = \tilde{\boldsymbol{\Theta}}\mathsf{Mean}(\mathbf{H}^{(2)}).
\]

Since $\mathsf{Mean}(\mathbf{H}^{(2)}) = (\mathbf{H}^{(2)})^\top\mathbf{1}/n$, we find

\[
\mathbf{Y}_2 = \boldsymbol{\Theta}_0 \mathsf{Mean}(\mathbf{H}) + \boldsymbol{\Theta}_1 \mathsf{Mean}(\mathbf{A}\mathbf{H}) + \boldsymbol{\Theta}_2 \mathsf{Mean}(\mathbf{A}^2\mathbf{H}),
\]

where the weights characterizing the family of quadratic polynomial filters are defined by the system:

\begin{align}\label{eq:quadratic_constraints}
\boldsymbol{\Theta}_0 &= \tilde{\boldsymbol{\Theta}}\Big(\mathbf{I} + (\mathbf{Q}^{(1)})^\top\Big), \\
\boldsymbol{\Theta}_1 &=\tilde{\boldsymbol{\Theta}}\Big((\mathbf{W}^{(0)})^\top + (\mathbf{W}^{(1)})^\top\Big), \\
\boldsymbol{\Theta}_2 &= \tilde{\boldsymbol{\Theta}}\Big((\mathbf{W}^{(1)})^\top(\mathbf{W}^{(0)})^\top \Big).
\end{align}

Assume now that the class of polynomial filters we wish to learn, have a vanishing zeroth-order term, i.e., $\boldsymbol{\Theta}_0 = \mathbf{0}$. For the MPNN not containing a VN, corresponding to $\mathbf{Q}^{(1)} =  \mathbf{0}$ in (<ref>), then we must have $\tilde{\boldsymbol{\Theta}} = \mathbf{0}$. Therefore, the only quadratic polynomial without a zeroth-order term that can be learned by a linear MPNN as in (<ref>) of degree 2 is the trivial polynomial. Conversely, an MPNN + VN can, for example, learn the weights $\mathbf{Q}^{(1)} = -\mathbf{I}$ to ensure $\boldsymbol{\Theta}_0 = \mathbf{0}$, and have non-trivial first and second-order terms using the $\mathbf{W}$-weight matrices.
This proves that on any given graph there exist polynomial filters that can be learned by MPNN + VN but not by MPNN.

To conclude the proof, we now need to show that MPNN + VN loses expressive power when it emulates PairNorm. We observe that for the case of MPNN - VN, i.e., $\mathbf{Q}^{(\ell)} = -\mathbf{I}$ for $\ell \in \mathbb{N}$, after two layers, the term $\boldsymbol{\Theta}_0$ is always zero. Since MPNN after two layers can learn quadratic polynomial filters with non-trivial zero-th order terms, this completes the proof.

§.§ Additional Results on the Role of Pooling

In this section we discuss how the choice of the global pooling affects the benefits of VN. Consider an alternative global pooling, which instead of projecting the final features onto $\mathbf{1}/n$ (i.e. computing the mean), projects the final-layer features onto $\mathbf{v}\in\mathbb{R}^n$ satisfying $\mathbf{v}^\top\mathbf{1} = 0,$

\begin{equation}\label{eq:pooling_v}
\mathbf{Y}_m^\mathbf{v} = \boldsymbol{\tilde{\Theta}}(\mathbf{H}^{(m)})^\top \mathbf{v} = \sum_{k=0}^{m}\boldsymbol{\Theta}_k(\mathbf{A}^k\mathbf{H})^\top \mathbf{v}.
\end{equation}

We can then show that there exist polynomial filters that MPNN + VN is no longer able to express when we replace mean pooling with $\mathbf{v}$-pooling.

Given $\mathbf{v}$ orthogonal to $\mathbf{1}$, a linear MPNN + VN of $m$ layers can generate weights $\{\boldsymbol{\Theta}_k\}$ in the output $\mathbf{Y}_m$ in (<ref>) that cannot be generated by a linear MPNN + VN using a global pooling induced by $\mathbf{v}$, as in (<ref>).

Once again, we show that it suffices to consider the case of quadratic polynomial filters, i.e., $m=2$ given in the system of equations in (<ref>). However, since $\mathbf{1}$ and $\mathbf{v}$ are orthogonal, then there is no $\mathbf{Q}^{(1)}$-term in <ref>. As such, we see that the quadratic polynomial $\mathbf{Y}_2^{\mathbf{v}}$ has zero-th order coefficient vanishing if and only if $\tilde{\boldsymbol{\Theta}} = \mathbf{0}$, i.e., the quadratic polynomial is trivial. Conversely, MPNN + VN can instead learn non-trivial quadratic polynomial without zero-th order term by simply setting $\mathbf{Q}^{(1)} = -\mathbf{I}$.

A Fourier Perspective on Pooling. Proposition <ref> highlights that MPNN + VN may lose the ability to express polynomial filters when we no longer have alignment between VN and the global pooling. This result raises a an important question: What is the bias resulting from a certain choice of global pooling? Consider linear MPNNs as above where the global pooling is obtained by projecting onto some $\mathbf{v}\in\mathbb{R}^n$, with the case $\mathbf{v} = \mathbf{1}$ recovering the mean (sum) pooling. Given $\{\psi_i\}\subset \mathbb{R}^{n}$, a basis of orthonormal eigenvectors of $\mathbf{A}$, with eigenvalues $\{\lambda_i\}$, we write $\mathbf{v} = \sum_i c_i \psi_i$, where $\{c_i\}$ are the graph Fourier coefficients of $\mathbf{v}$. Given $1\leq r \leq d$, the channel $r$ of the final output as per (<ref>) takes on the following form in the Fourier domain

\begin{equation}
y_r = \sum_{k=0}^{m}\sum_{p=1}^{d}(\Theta_k)_{rp}\sum_{i=0}^{n-1}c_i\lambda_i^k \langle h^p, \psi_i\rangle,
\end{equation}

where $h^p\in\mathbb{R}^n$ is the channel $p$ of the input features. Accordingly, the Fourier coefficients $\{c_i\}$ of $\mathbf{v}$ and the associated frequencies $\{\lambda_i\}$, determine which projection $\langle h^p, \psi_i\rangle$ affects the output the most.
The fact that on many datasets mean (sum) pooling is beneficial, and hence adding VN is mostly advantageous as per Proposition <ref>, highlights that the underlying tasks depend on the projection of the features onto those eigenvectors $\{\psi_i\}$ aligned with the vector $\mathbf{1}$.

§ PROOFS OF RESULTS IN SECTION <REF>

In this section, we report the proofs of the theoretical results in Section <ref> along with additional details. We start by proving <Ref>.

We start by recalling that the effective resistance between two nodes $i$ and $j$, denoted by $R(i,j)$, can be computed as [Ghosh et al., 2008]:

\begin{equation}\label{eq:formula_ER_ij}
R(i,j)   = L_{ii}^\dagger + L_{jj}^\dagger - 2L_{ij}^\dagger,
\end{equation}

where $\mathbf{L}^\dagger$ is the pseudo-inverse of the graph Laplacian, i.e.,

\begin{equation}
\mathbf{L}^\dagger = \sum_{\ell = 1}^{n-1}\frac{1}{\lambda_\ell}v_\ell v_\ell^\top.
\end{equation}

Accordingly, we can write the effective resistance between two nodes $i$ and $j$ explicitly as

\begin{equation}
R(i,j) = \sum_{\ell = 1}^{n-1}\frac{1}{\lambda_\ell}(v_\ell(i) - v_\ell(j))^2.
\end{equation}

Naturally, for the graph $G_{\mathsf{vn}}$ we obtain from adding VN, the formula changes as

\begin{equation}
R_{\mathsf{vn}}(i,j) = \sum_{\ell = 1}^{n}\frac{1}{\lambda^{\mathsf{vn}}_\ell}(v^{\mathsf{vn}}_\ell(i) - v^{\mathsf{vn}}_\ell(j))^2,
\end{equation}

where $\mathbf{L}_{\mathsf{vn}}v^{\mathsf{vn}}_\ell = \lambda^{\mathsf{vn}}_\ell v^{\mathsf{vn}}_\ell$. When we add the virtual node though, the adjacency matrix changes as in <ref>, meaning that the graph Laplacian can be written as

\begin{equation}
\mathbf{L}_{\mathsf{vn}} = \mathbf{D}_{\mathsf{vn}} - \mathbf{A}_{\mathsf{vn}} = \begin{pmatrix}
\mathbf{D} + \mathbf{I} & \mathbf{0} \\ \mathbf{0}^\top & n
\end{pmatrix} - \begin{pmatrix}
\mathbf{A} & \mathbf{1} \\ \mathbf{1}^\top & 0
\end{pmatrix} = \begin{pmatrix}
\mathbf{L} + \mathbf{I} & -\mathbf{1} \\ -\mathbf{1}^\top & n
\end{pmatrix}.
\end{equation}

We recall now that $v_0 = \frac{1}{\sqrt{n}}\mathbf{1}$, i.e., $\mathbf{L}\mathbf{1} = \mathbf{0}$. Accordingly, we find that $v^{\mathsf{vn}}_0 = \frac{1}{\sqrt{n+1}}\begin{pmatrix}\mathbf{1} \\ 1\end{pmatrix}$, i.e., $\mathbf{L}_{\mathsf{vn}}v^{\mathsf{vn}}_0 = \mathbf{0}$. Similarly, from the block decomposition of $\mathbf{L}_{\mathsf{vn}}$ we derive that

\begin{equation}
\mathbf{L}_{\mathsf{vn}}\begin{pmatrix}
v_\ell \\ 0
\end{pmatrix} = \begin{pmatrix}
(\mathbf{L} + \mathbf{I})v_\ell \\ -\mathbf{1}^\top v_\ell
\end{pmatrix} = (\lambda_\ell + 1)\begin{pmatrix}
v_\ell \\ 0
\end{pmatrix},
\end{equation}

where we have used that $v_\ell$ are orthogonal to $v_0 = \frac{1}{\sqrt{n}}\mathbf{1}$ for all $\ell \geq 1$. Accordingly $\begin{pmatrix}
v_\ell \\ 0
\end{pmatrix}$ are eigenvectors of $\mathbf{L}_{\mathsf{vn}}$ for all $1\leq \ell \leq n-1$ with eigenvalue $\lambda_\ell + 1$. Finally we note that

\begin{equation}
\mathbf{L}_{\mathsf{vn}}\frac{1}{\sqrt{n^2 + n}}\begin{pmatrix}
\mathbf{1} \\ -n
\end{pmatrix} = \frac{1}{\sqrt{n^2 + n}}\begin{pmatrix}
(\mathbf{L} + \mathbf{I})\mathbf{1} +n\mathbf{1} \\ -\mathbf{1}^\top\mathbf{1} -n^2
\end{pmatrix} =(1+n) \frac{1}{\sqrt{n^2 + n}} \begin{pmatrix}
\mathbf{1} \\ -n
\end{pmatrix},
\end{equation}

which means that $\frac{1}{\sqrt{n^2 + n}}\begin{pmatrix}
\mathbf{1} \\ -n
\end{pmatrix}$ is the final eigenvector of $\mathbf{L}_{\mathsf{vn}}$ with eigenvalue $n+1$.
We can then use the spectral decomposition of $\mathbf{L}_{\mathsf{vn}}$ to write the effective resistance between nodes $i$ and $j$ after adding VN as

\begin{align*}
R_{\mathsf{vn}}(i,j) &= \sum_{\ell = 1}^{n}\frac{1}{\lambda^{\mathsf{vn}}_\ell}(v^{\mathsf{vn}}_\ell(i) - v^{\mathsf{vn}}_\ell(j))^2 = \sum_{\ell = 1}^{n-1}\frac{1}{\lambda_\ell + 1}(v_\ell(i) - v_\ell(j))^2 + \frac{1}{n+1}\frac{1}{n^2 + n}(1 - 1)^2 \notag \\
&= \sum_{\ell = 1}^{n-1}\frac{1}{\lambda_\ell + 1}(v_\ell(i) - v_\ell(j))^2.
\end{align*}

Finally, we recall that on a graph, the commute time $\tau$ is proportional to the effective resistance [Lovász, 1993], i.e.,

\begin{equation*}
\tau(i,j) = 2|E|R(i,j).
\end{equation*}

Accordingly, we can use that on $G_{\mathsf{vn}}$ the number of edges is equal to $|E| + n$ and compare the commute time after adding VN with the baseline case:

\begin{align*}
\tau_{\mathsf{vn}}(i,j) - \tau(i,j) &=  \sum_{\ell = 1}^{n-1}\Big(\frac{2(|E| + n)}{\lambda_\ell + 1} - \frac{2|E|}{\lambda_\ell}\Big)(v_\ell(i) - v_\ell(j))^2 \\
&= 2|E|\sum_{\ell = 1}^{n-1}\frac{1}{\lambda_\ell(\lambda_\ell + 1)}\Big( \frac{n}{|E|}\lambda_\ell - 1\Big)(v_\ell(i) - v_\ell(j))^2,
\end{align*}

which is precisely <ref>. We can the sum over all pairs of nodes to obtain

\begin{align*}
\sum_{i,j=1}^n(\tau_{\mathsf{vn}}(i,j) - \tau(i,j)) &= 2|E|\sum_{i,j=1}^n\sum_{\ell = 1}^{n-1}\frac{1}{\lambda_\ell(\lambda_\ell+1)}\Big(\frac{n}{|E|}\lambda_\ell - 1\Big)(v_\ell(i)-v_\ell(j))^2 \\
&= 2|E|\sum_{\ell = 1}^{n-1}\frac{1}{\lambda_\ell(\lambda_\ell+1)}\Big(\frac{n}{|E|}\lambda_\ell - 1\Big)\sum_{i,j=1}^n(v_\ell(i)-v_\ell(j))^2 \\
&= 2|E|\sum_{\ell = 1}^{n-1}\frac{1}{\lambda_\ell(\lambda_\ell+1)}\Big(\frac{n}{|E|}\lambda_\ell - 1\Big)\sum_{i,j=1}^n\Big(v^2_\ell(i) + v^2_\ell(j) -2v_\ell(i)v_\ell(j)\Big).
\end{align*}

We now use that $v_\ell$ form an orthonormal basis, accordingly: $\sum_{i=1}^n v_\ell(i)^2 = 1$. Similarly, since $v_\ell$ is orthogonal to $v_0$ and hence to $\mathbf{1}$, we get $v_\ell^\top \mathbf{1} = \sum_{i=1}^n v_\ell(i) = 0$. Therefore,

\begin{equation*}
\sum_{i,j=1}^n \Big(v^2_\ell(i) + v^2_\ell(j) -2v_\ell(i)v_\ell(j)\Big) = 2n -2\sum_{i=1}^n v_\ell(i) \sum_{j=1}^n v_\ell(j) = 2n.
\end{equation*}

We finally  obtain:

\begin{equation*}
\frac{1}{n^2} \sum_{i,j=1}^n (\tau_{\mathsf{vn}}(i,j) - \tau(i,j)) = \frac{4|E|}{n}\sum_{\ell = 1}^{n-1}\frac{1}{\lambda_\ell(\lambda_\ell+1)}\Big(\frac{n}{|E|}\lambda_\ell - 1\Big),
\end{equation*}

which completes the proof.

Remark: We note that according to (<ref>), we expect $\tau_{\mathsf{vn}}(i,j) < \tau(i,j)$ whenever the two nodes have large commute time on the input graph, i.e. due to the graph having a bottleneck where $\lambda_1$ is small and the Fiedler eigenvector is such that $v_1(i)\cdot v_i(j) < 0$ since $i,j$ belong to different communities. As such, adding a VN should indeed beneficial for graphs with bottlenecks, while it may slightly increase the commute time between nodes that were well connected in the graph due to the addition of a VN spreading information out from the pair. This intuition is validated in <Ref>, where we show that on real-world graphs adding a VN decreases the average commute time significantly.

If we let $\alpha:= 1 + n/|E|$, then

\begin{equation}
-\frac{4|E|\alpha}{\lambda_1(\lambda_1 + 1)} \leq \tau_{\mathsf{vn}}(i,j) -\alpha\tau(i,j) \leq -\frac{4|E|\alpha}{\lambda_{n-1}(\lambda_{n-1}+1)}.
\end{equation}

From (<ref>) and the proof above, we derive

\begin{align*}
\tau_{\mathsf{vn}}(i,j) - \alpha\tau(i,j) &= \sum_{\ell=1}^{n-1}\Big(\frac{2(|E| + n)}{\lambda_\ell + 1} - \frac{2(|E|+n)}{\lambda_\ell}\Big)(v_\ell(i) - v_\ell(j))^2  \\
&= -2(|E| + n)\sum_{\ell = 1}^{n-1}\frac{1}{\lambda_\ell(\lambda_\ell + 1)}(v_\ell(i) - v_\ell(j))^2 \\
&\geq -\frac{2|E|\alpha}{\lambda_1(\lambda_1 + 1)}\sum_{\ell=1}^{n-1}(v_\ell(i) - v_\ell(j))^2 = -\frac{4|E|\alpha}{\lambda_1(\lambda_1 + 1)}
\end{align*}

where the last equality follows from

\begin{align*}
\sum_{\ell = 1}^{n-1}(v_\ell(i) - v_\ell(j))^2 &= \sum_{\ell = 1}^{n-1}(v^2_\ell(i) + v^2_\ell(j) -2v_\ell(i)v_\ell(j))
\\
&= 2\Big(1 - \frac{1}{n}\Big) - 2\sum_{\ell = 0}^{n-1}v_\ell(i)v_\ell(j) + 2v_0(i)v_0(j) = 2\Big(1 - \frac{1}{n}\Big) + \frac{2}{n} = 2,
\end{align*}

where we have used the orthonormality of the eigenvectors for $i\neq j$. The upper bound can be proved in the same way using that $\lambda_\ell \leq \lambda_{n-1}$ for each $\ell$.

Consider an MPNN + VN. We recall that $\mathsf{mix}_y(i,j)$ is the mixing induced by a graph-level function among the features of nodes $i$ and $j$, and is defined as [Di Giovanni et al., 2023]

\[
\mathsf{mix}_y(i,j) = \max_{\mathbf{H}}\max_{1\leq\alpha,\beta\leq d}\left| \frac{\partial^2 y(\mathbf{H})}{\partial h_i^\alpha \partial h_j^\beta} \right|.
\]

According to Di Giovanni et al., 2023, when the mixing induced by the MPNN after $m$ layers, denoted here as $\mathsf{mix}^{(m)}(i,j)$, is smaller than the one required by the downstream task, then we have an instance of harmful oversquashing which prevents the network from learning the right interactions necessary to solve the problem at hand. Since we can write the system in (<ref>) as the system of MPNN analyzed in Di Giovanni et al., 2023, operating on the graph $G_{\mathsf{vn}}$, we can extend their conclusions in Theorem 4.4. Namely, using their notations, we find that given bounded weights and derivatives for their network, the minimal number of layers required to learn a graph-level function $y$ with mixing $\mathsf{mix}_y(i,j)$ is

\[
m \geq \frac{\tau_{\mathsf{vn}}(i,j)}{8} + \alpha_{\mathsf{vn}}\mathsf{mix}_{y}(i,j) + \beta_{\mathsf{vn}},
\]

where the constants $\alpha_{\mathsf{vn}}$ and $\beta_{\mathsf{vn}}$ are independent of nodes $i,j$, and without loss of generality we have taken $c_2 = 1/2$ in Theorem 4.4 of Di Giovanni et al., 2023— clearly, the result generalize to any $c_2\leq 1$. As such, we see that a necessary condition for the family of graph-level functions with mixing $\mathsf{mix}_y(i,j) > -\beta_{\mathsf{vn}}/\alpha_{\mathsf{vn}}$ to be learned by MPNN + VN is that

\[
m \geq \frac{\tau(i,j)}{8} \geq \Big(1 + \frac{n}{|E|}\Big)\Big(\frac{\tau(i,j)}{8} - \frac{|E|}{2\lambda_1(\lambda_1 + 1)}\Big),
\]

where we have used <ref> in the last inequality. This completes the proof.

Accordingly, MPNN + VN can improve upon the minimal number of layers required to learn functions with strong mixing between nodes $(i,j)$, when compared to the baseline MPNN. In fact, such improvement can be characterized precisely in terms of spectral properties of the graph Laplacian, and in light of <Ref> of $\lambda_1$.

§ PROOFS OF RESULTS IN SECTION <REF>

Recall that we consider an MPNN + VN whose layer update is of the form

\begin{align*}
h_{\mathsf{vn}}^{(\ell+1)} &= \sigma\Big(\boldsymbol{\Omega}_{\mathsf{vn}}h_{\mathsf{vn}}^{(\ell)} + \mathbf{W}_{\mathsf{vn}}^{(\ell)}\frac{1}{n}\sum_{j=1}^{n}h_j^{(\ell)}\Big) \notag \\
h_i^{(\ell+1)}&= \sigma\Big(\boldsymbol{\Omega}^{(\ell)}h_i^{(\ell)} + \sum_{j=1}^n\mathsf{A}_{ij}\psi^{(\ell)}(h_i^{(\ell)},h_j^{(\ell)}) + \psi^{(\ell)}_{\mathsf{vn}}(h_i^{(\ell)}, h_{\mathsf{vn}}^{(\ell)})\Big).
\end{align*}

If $i,k\in V$ with $k$ outside the 2-hop neighborhood of $i$, then any message sent from $k$ to $i$ will arrive after 2 layers through the VN. Since the VN update is homogeneous though, node $i$ cannot distinguish which node $k$ sent the message. Explicitly, we can compute the derivatives and obtain

\begin{align*}
\frac{\partial h_i^{(\ell+1)}}{\partial h_k^{(\ell-1)}} = \frac{1}{n}\sigma'(z_i^{(\ell)})\nabla_{2}\psi_{\mathsf{vn}}^{(\ell)}(h_i^{(\ell)},h_\mathsf{vn}^{(\ell)})\sigma'(z_{\mathsf{vn}}^{(\ell-1)})\mathbf{W}_{\mathsf{vn}}^{(\ell)},
\end{align*}

where $z_i^{(\ell)}$ is the argument of the nonlinear activation $\sigma$, $\sigma'(z)$ is the diagonal derivative matrix computed at $z$ and $\nabla_2$ denotes the Jacobian w.r.t. the second variable. As we can see, the derivative is independent of $k$ and is the same for each node $k$ outside the 2-hop neighborhood of $i$, which shows that the layer of MPNN + VN is typically homogeneous and hence fails to assign different relevance (sensitivity) to nodes. This completes the proof.

Consider an instance of $\acro$ whose layer updates have the form

\begin{align*}
h_i^{(\ell+1)} &= h_{i,{\rm loc}}^{(\ell+1)} + \mathsf{Mean}(\{h_{j,{\rm loc}}^{(\ell+1)}\}), \\
h_{j,{\rm loc}}^{(\ell+1)} &= \sigma(\boldsymbol{\Omega}^{(\ell)}h_i^{(\ell)} + \sum_{j\in N_i} \psi_{ij}^{(\ell)}(h_i^{(\ell)},h_j^{(\ell)})).
\end{align*}

Given $i,k\in V$ where node $k$ is outside the 1-hop neighborhood of $i$, any message sent from $k$ to $i$ will arrive after 1 layer through the VN—in other words, $\partial h_{i,{\rm loc}}^{(\ell+1)}/\partial h_k^{(\ell)} = 0$.  We can compute the derivatives and obtain

\begin{align*}
\frac{\partial h_i^{(\ell+1)}}{\partial h_k^{(\ell)}} &= \frac{1}{n}\sum_{j=1}^n \frac{\partial{h_{j,{\rm loc}}^{(\ell+1)}}}{{\partial h_k^{(\ell)}}} \\
&=\frac{1}{n}\sum_{j=1}^n \sigma'(z_j^{(\ell)})\Big(\boldsymbol{\Omega}^{(\ell)}\delta_{jk} + \sum_{u\in N_j}\nabla_1\psi_{ju}(h^{(\ell)}_j,h^{(\ell)}_u)\delta_{jk} + \nabla_2\psi_{ju}(h^{(\ell)}_j,h^{(\ell)}_u)\delta_{uk} \Big) \\
&=\frac{1}{n}\Big(\sigma'(z_k^{(\ell)})\Big(\boldsymbol{\Omega}^{(\ell)} + \sum_{u\in N_k}\nabla_1 \psi^{(\ell)}_{ku}(h_k^{(\ell)},h_u^{(\ell)})\Big) + \sum_{u\in N_k}\sigma'(z_u^{(\ell)})\nabla_2\psi_{uk}^{(\ell)}(h_u^{(\ell)},h_k^{(\ell)}) \Big),
\end{align*}

where $\delta_{pq}$ is the Kronecker delta with indices $p,q$ and in the last sum we have replaced the dumb index $j$ with $u$. This completes the proof.

§ EXPERIMENTAL DETAILS

In this section, we provide further details about our experiments.

§.§ Hardware

All experiments were run on a single V100 GPU.

§.§ Description of Datasets

Below, we provide descriptions of the datasets on which we conduct experiments.

MNIST and CIFAR10 (CC BY-SA 3.0 and MIT License) [Dwivedi et al., 2023]. These datasets test graph classification on the popular image classification datasets. The original images are converted to graphs using super-pixels which represent small regions of homogeneous intensity in the images. They are both 10-class classification tasks and follow the original standard (train/validation/test) splits; 55K/5K/10K for MNIST and 45K/5K/10K for CIFAR10.

ogbg-molhiv and ogbg-molpcba (MIT License) . These are molecular property prediction datasets which use a common node and edge featurization that represents chemophysical properties. Ogbg-molhiv is a binary classification task, predicting the molecule's ability to inhibit HIV replication and ogbg-molpcba is a multi-task binary classification where 128 bioassays are predicted.

ogbg-ppa (CC-0 License) [Hu et al., 2020]. This dataset consists of protein-protein interaction networks derived from 1581 species and categorized into 37 taxonomic groups. Nodes represent proteins and edges encode the normalized level of 7 different associations between proteins. The task is to classify which of the 37 groups the network belongs to.

MalNet-Tiny (CC-BY License) [Freitas et al., 2020]. This is a subset of MalNet which contains function call graphs (FCGs) derived from Android APKs. There are 5,000 graphs with up to 5,000 nodes with each graph coming from a benign software. The goal is to predict the type of software from the structure of the FCG. The benchmarking version of this dataset typically uses Local Degree Profile as the set of node features.

Peptides-func and Peptides (CC-BY-NC 4.0) [Dwivedi et al., 2022]. These datasets are composed of atomic peptides. Ppetides-func is a multi-label graph classification task where there are 10 nonexclusive peptide functional classes. Peptides-struct is a regression task involving 11 3D structural properties of the peptides.

§.§ Dataset Splits and Random Seeds

All the benchmarks follow the standard train/validation/test splits. The test performance at the epoch with the best validation performance is reported and is averaged over multiple runs with different random seeds. All the benchmarking results, including the extra ablations, are based on 10 executed runs, except for Peptides-func and Peptides-struct which are based on the output of four runs.

§.§ Hyperparameters

Considering the large number of hyperparameters and datasets, it was not possible to do an exhaustive grid search and to find the optimal parameters. Here we describe how the final hyperparameters shown in Tables <ref> and <ref> were obtained.

OGB datasets. For both ogbg-molpcba and ogbg-ppa we used the same hyperparameters as used in Rampášek et al., 2022 but the hidden dimension of the MPNN was adjusted so that we had a similar parameter budget. For ogbg-molhiv, we found it to be beneficial to reduce the number of layers to 4 to align with the number of layers used in Bouritsas et al., 2023 but kept the same parameters for the positional encoding and the optimization process.

Peptides-Func and Peptides-Struct.
For these datasets we optimized the hyperparameters over the following ranges:

* Dropout [0, 0.1, 0.2],

* FeedForward Block [True, False],

* Depth [4, 6, 8, 10],

* Positional Encoding [none, LapPE, RWSE],

* Layers Post Message-Passing [1, 2, 3],

and we used the optimization parameters recently proposed by [Tönshoff et al., 2023] where they train for 250 epochs with an AdamW optimizer [Kingma and Ba, 2014], and a cosine annealing learning rate scheduler with a base learning rate of 0.001. When optimizing over the number of layers of message-passing, we changed the hidden dimension to ensure that the parameter budget was around 500K.

MNIST, CIFAR10, MalNet-Tiny. On these benchmarks, we used the same dropout, positional encondings and optimization parameters as used by Shirzad et al., 2023. The only parameter we optimised for was the number of layers in the range [3, 5, 7]. Additionally, we changed the hidden dimension in accordance with the number of layers to match the number of parameters that were used in Shirzad et al., 2023.

Best performing hyperparameters for GatedGCN+PE+$\hetvn$ in Table <ref>.

Hyperparameter     Peptides-Func     Peptides-Struct      MNIST      CIFAR10     MalNet-Tiny

#Layers        4       4         5       7       5

Hidden dim       136         136        46         40       72

Dropout        0.0        0.2         0.1       0.1         0.0

Graph pooling       mean       mean           mean       mean       mean

FeedForward Block       False          False        True     True     True

Positional Encoding      RWSE-20        LapPE-16          LapPE-8          LapPE-8       None

PE dim         16         16         8        8      -

PE encoder     Linear          Linear           Linear      DeepSet     -

Batch size      200        200          16         16        16

Learning Rate     0.001        0.001            0.001      0.001      0.0005

#Epochs        250        250         100       100           150

#Parameters      492,217       492,217           110,148        117,066       281,453

Best performing hyperparameters for GatedGCN+PE+$\hetvn$ in Table <ref>.

Hyperparameter     ogbg-molhiv     ogbg-molpcba      ogbg-ppa

#Layers         4      5          3

Hidden dim       90       384        284

Dropout        0.0      0.2        0.1

Graph pooling      mean       mean        mean

FeedForward Block      True       True       True

Positional Encoding     RWSE-16       RWSE-16        None

PE dim         16      20       -

PE encoder     Linear       Linear        -

Batch size      32      512        32

Learning Rate     0.0001       0.0005       0.0003

#Epochs        100      100        200

#Parameters     360,221       7,519,580      2,526,501

§ SMOOTHING AND THE IMPORTANCE OF GRAPH STRUCTURE

In this section, we explore the performance on various benchmarks of adding or removing the current pooled representation at each layer. We take the pooling function to be the mean and at each layer we either subtract or add this from the output of the message passing layer. Adding the mean increases the smoothness of the representations and we see that this leads to performance improvements on the LRGB tasks in Table <ref>. Again, this highlights the benefits of smoothing on some graph-level tasks. On the other hand, we see that subtracting the mean is beneficial for ogbg-molpcba. Subtracting the mean at each layer before using the mean pooling function in the final layer, means that we are ignoring the output of the message-passing layer which is aligned with the final representation. We can see this as a good measure of the importance of the graph structure as we are removing any locally aggregated updates. This suggests that it is beneficial to ignore the graph topology for ogbg-molpcba.

Analyzing the performance of GatedGCN where we subtract or add the mean of the message-passing output at each layer.

Method     Peptides-Func ($\uparrow$)     Peptides-Struct ($\downarrow$)     ogbg-molhiv ($\uparrow$)      ogbg-molpcba ($\uparrow$)

GatedGCN                 0.5864 $\pm 0.0077$            0.3420 $\pm 0.0013$       0.7827 $\pm 0.0111$      0.2714   $\pm 0.0014$

GatedGCN  + mean                 0.6692 $\pm 0.0042$            0.2522 $\pm 0.0012$       0.7677  $\pm 0.0138$       0.2569  $\pm 0.0034$

GatedGCN  - mean                 0.4675 $\pm 0.0040$            0.3566  $\pm 0.0007$       0.7594 $\pm 0.0096$     0.2866  $\pm 0.0016$

§.§ Cosine Distance Plots

We analyzed the layerwise smoothing on the graph-level LRGB datasets [Dwivedi et al., 2022] using a GatedGCN with and without a VN/attention layer. To do this, we measured the cosine distance between the representations at each layer and the final pooled graph representation. We found that using a VN or an attention layer reduced the distance between the earlier layer representations and the final pooled representation and that this lead to an improvement in performance. This suggests that these approaches can cause a beneficial smoothing towards the final representation.

The cosine distance between the embedding at each layer and the final graph representation after training a GatedGCN, as well as with a VN and a transformer layer on Peptides-func.

The cosine distance between the embedding at each layer and the final graph representation after training a GatedGCN, as well as with a VN and a transformer layer on Peptides-struct.

§ ADDITIONAL ABLATION STUDIES

In this section, we explore additional ablations to improve our understanding of virtual nodes and our heterogeneous extension. In Table <ref>, we look at the performance improvement of our heterogeneous $\text{VN}_G$ when the base MPNN is a GCN or a GatedGCN. Whilst we generally see an improvement with our implementation, this improvement is much larger when we use a GatedGCN. As previously argued, this is due to the fact that in a GatedGCN, the aggregation weights each neighbor through a learnable gate and we can thus have a learnable node importance based on the graph topology. This is in contrast to a GCN where we have a homogeneous node update.

Performance of $\acro$ in comparison to virtual node with positional encodings.

Method     ogbg-molhiv ($\uparrow$)      ogbg-molpcba ($\uparrow$)     peptides-func ($\uparrow$)     peptides-struct ($\downarrow$)     CIFAR10 ($\uparrow$)

GCN+PE+VN                 0.7599 $\pm 0.0119 $        0.2456 $\pm 0.0034$        0.6562 $\pm 0.0070$       0.2540 $\pm 0.0008$      68.957 $\pm 0.381$

GCN+PE+$\hetvn$              0.7678 $\pm
0.0111$            0.2481 $\pm 0.0032$        0.6537 $\pm 0.0026$            0.2497 $\pm 0.0008$     68.756 $\pm 0.172$

GatedGCN+PE+VN              0.7687 $\pm
0.0136$            0.2848 $\pm 0.0026$        0.6712 $\pm 0.0066$            0.2481 $\pm 0.0015$     70.280 $\pm 0.380$

GatedGCN+PE+$\hetvn$             0.7884 $\pm
0.0099$            0.2917 $\pm 0.0027$        0.6822 $\pm 0.0052$            0.2458 $\pm 0.0006$     76.080 $\pm 0.301$

GCN % Increase              +1.04            +1.02         -0.38           +1.69     -0.29

GatedGCN % Increase             +2.56            +2.42        +1.64            +0.89     +8.25

As an extension to Section <ref>, we compared the performance of GatedGCN with augmentations which involve adding a VN
in Table <ref>.
Additionally, we look at the effect of applying PairNorm on the datasets. We see that applying PairNorm to the GatedGCN, a common technique to mitigate oversmoothing, actually reduces performance on all of these datasets. Moreover, using a virtual node always outperforms using PairNorm. This further implies that a VN is not recreating PairNorm and suggests that, on these graph-level tasks, smoothing may be beneficial.

Performance of GatedGCN and its extensions on four benchmark tasks.

Method     Peptides-Func ($\uparrow$)     Peptides-Struct ($\downarrow$)     ogbg-molhiv ($\uparrow$)      ogbg-molpcba ($\uparrow$)

GatedGCN                 0.5864 $\pm 0.0077$            0.3420 $\pm 0.0013$       0.7827 $\pm 0.0111$      0.2714   $\pm 0.0014$

GatedGCN  with PairNorm                 0.4674 $\pm 0.0040$            0.3551 $\pm 0.0008$      0.7645 $\pm 0.0128$     0.2621 $\pm 0.0026$

GatedGCN + VN                 0.6477 $\pm 0.0039$            0.2523 $\pm 0.0016$      0.7676  $\pm
0.0172$            0.2823 $\pm 0.0026$

§ VISUALIZATION OF $\ACRO$

Here we visualize the differences between our implementation of $\acro$ and a standard MPNN+VN. Our model uses the local MPNN acting on the original graph to weight the nodes based on their importance. This means that the VN can perform a heterogeneous global update.

Comparing MPNN+VN with our proposed $\acro$.

§ ATTENTION MAPS

In order to get a better understanding of the homogeneity of the attention layer in the GPS framework, we visualized the first layer attention matrices for various datasets. This complements our analysis in Section <ref> where we relate the gap in performance between a MPNN + VN and GPS to the level of homogeneity of these attention patterns. The attention pattern for CIFAR10 is less homogeneous, as measured by the standard deviation of the columns sums.

First layer attention maps of the self-attention matrix in the GPS framework for different datasets.