1 Introduction↩︎

Graph Neural Networks (GNNs) [1]–[3] have become the dominant framework for learning on graph structured data and have achieved widespread success in domains as varied as protein structure prediction [4], material science [5], weather forecasting [6] or catalyst screening [7].

A key theoretical pillar supporting their widespread adoption is the commonly held belief that GNNs define stable and continuous mappings: small changes in the input graph are expected to induce correspondingly small changes in learned representations. Indeed, a substantial body of prior work has established continuity, stability, and robustness for GNNs under a broad range of perturbation modes: Previous investigations e.g. considered GNN-stability under small variations in edge-weights [8], [9] or graph rewiring [10]. Consistency of latent representations has also been established across graphs that faithfully approximate the same metric measure space [11], graphon [12], [13], graphing [14] or graphop [15].¹ Taken together, these results have contributed to a prevailing view in the field that graph neural networks define broadly stable and well-behaved mappings: similar input graphs are mapped to similar latent embeddings.

In this work, we show that this is, however, not always the case: In particular we demonstrate that standard GNN architectures can produce drastically different latent representations for graphs that describe the same underlying object at different resolution scales. As a consequence, models trained on graphs at one resolution, for instance, fail to generalize to equivalent graphs at another resolution.

We trace this failure to the previously understudied fact that graph neural networks are generally not continuous with respect to all natural modes of graph convergence: We observe that there is a natural interpolating sequence of graphs starting at an original graph and terminating at its coarse-grained version. Mathematically, this sequence smoothly interpolates between the original- and coarse heat kernel. This convergence is not respected by GNNs: Even though input graphs converge (in the heat sense) the generated sequence of latent embeddings does not converge.

This observation reveals a fundamental limitation of existing GNN architectures. It also motivates the development of modifications to existing architecures, in order to render them continuous across graph resolution scales.

Concretely, our paper makes the following contributions:

• We identify and empirically demonstrate a previously underappreciated failure mode of standard GNNs, i.e., a lack of continuity across graph resolutions.

• We provide a theoretical explanation linking this failure mode to discontinuity under graph convergence.

• We show that the underlying obstruction is structural, originating in standard propagation-schemes.

• We show how to modify propagation schemes to yield provably scale-continuous GNN architectures.

• We establish theoretical guarantees and present extensive empirical validation demonstrating our modifications facilitate reliable generalization across scales.

2 Graphs and Graph Learning↩︎

We first briefly review graphs and graph neural networks:

2.1 Graphs and their Description↩︎

We consider finite, weighted graphs \(G = (V,E)\) equipped with a node weight matrix \(M = \operatorname{diag}(\mu_v)_{v\in V}\). Each weight \(\mu_v > 0\) reflects node-wise information. Examples are physical mass or charge of atom \(v\) in a molecular graph, number of users represented if \(v\) is a node in a social network, or a volume element if \(v\) is part of a point cloud. Often, all node masses are set to \(\mu_v = 1\) so that \(M = I\). Edges \((u,v)\in E\) of a graph carry symmetric non-negative (potentially binary) weights \(A_{uv}\) which we collect in the adjacency matrix \(A\). The corresponding degree matrix is \(D = \operatorname{diag}(d_v)\) with \(d_v = \sum_{u} A_{uv}\).

The natural self-adjoint Laplacian associated with a node-weighted graph is \(L = M^{-1}(D - A)\). It is self-adjoint with respect to the \(M\)-weighted inner product \(\langle f,g\rangle_M = f^\top M g\). When the node weights encode local volumes or sampling densities it may be thought of as the discrete counterpart of the Laplacian \(L = - \Delta\) on a continuum domain [16]. While the natural Laplacian is central to our work, more commonly used in machine learning applications is instead the normalized Laplacian \(\mathcal{L} = I - D^{-1/2} A\, D^{-1/2}\). It is self-adjoint in the standard Euclidean inner product and its spectrum is contained in the interval \([0,2]\).

2.2 Graph Neural Networks↩︎

Graph neural networks are machine learning methods adapted to handle graph structured data. Corresponding network architectures broadly fall into two categories:

2.2.0.1 Message Passing Networks:

Message Passing Neural Networks (MPNNs) [3], [17] are arguably the most widely used class of GNNs. At each layer \(k\), node features \(h_v^{(k)}\) are updated by aggregating information from neighbouring nodes: \[\begin{align} m_{vu}^{(k)} &= \phi(h_v^{(k)}, h_u^{(k)}, A_{uv}) \label{eq:mf95form} \\ h_v^{(k+1)} &= \text{AGG}_{u \in \mathcal{N}(v)} m_{vu}. \end{align}\tag{1}\] Here \(A_{uv}\) denotes edge weights. Architectutes differ in which message- and aggregation functions are being used [18]–[21], with typcial choices for the neighborhood aggregation function \(\operatorname{AGG}\) being \(\operatorname{sum}\), \(\operatorname{avg}\), or (seldom) \(\operatorname{max}\) [22].

2.2.0.2 Spectral Networks:

Spectral graph concolutional layers [2], [23] apply learnable filters \(g_\theta(\mathcal{L})\) to node-wise information via matrix-multiplication. The filters themselves are composed of learnable scalar functions \(g_{\theta}(\cdot)\) which are applied to the normalized graph Laplacian \(\mathcal{L}\). In practice, the filter \(g_\theta\) is not implemented via eigendecomposition but approximated via polynomials [24]–[27]. The prototypical example is ChebNet, where filters are parametrized as \(g_\theta(\lambda) \;\approx\; \sum_{\ell=0}^{K} \theta_\ell\, T_\ell(\lambda),\) with \(T_\ell\) the \(\ell^{\text{th}}\)-Chebyshev polynomial. With node feature matrix \(X\) and learnable weight matrices \(\{W_\ell\}_{\ell}\) , a spectral layer is then efficiently implemented as \[\begin{align} \label{eq:spectral95prop} X^{(k+1)} = \sum_{\ell=0}^{K} T_\ell(\mathcal{L})\, X^{(k)}\, W_\ell . \end{align}\tag{2}\]

3 The Failure Mode: Lack of Generalization↩︎

To showcase how standard graph neural networks are unable to consistently integrate varying scales experimentally, we make use of the QM\(7\) dataset [28], consisting of organic molecules containing both hydrogen and heavy atoms. Prediction target is the molecular atomization energy. Each molecule is represented by a weighted adjacency matrix whose entries \(A_{ij} = Z_iZ_j\cdot|\vec{x}_i-\vec{x}_j|^{-1}\) correspond to Coulomb energies between atoms \(i,j\), with \(|\vec{x}_i-\vec{x}_j|\) denoting the interatomic distance and \(Z_i\) atomic charge. From a physical perspective, describing a molecule at the level of interacting atoms corresponds to a specific choice of resolution scale: Interactions of individual protons and neutrons inside the various atomic nuclei are discarded.² Instead, only an aggregate description is used wherein each nucleus is only described by a single node.

To test the ability of GNNs to generalize across scales, we additionally also consider a version of QM\(7\) where we lower the resolution scale even further: To this end we aggregate each heavy atomic core together with its surrounding (single-proton) hydrogen atoms into super-nodes. As depicted in Fig. 1, this is done via standard graph coarsification [29], [30]; Appendix 16.1 provides exact details.

This QM\(7_{\text{coarse}}\) dataset models data obtained from a resolution-limited observation process unable to resolve positions of individual (small) hydrogen atoms and only providing information about how many are bound to a given heavy atom. Using the high-resolution graphs \(\{G\}\) of QM\(7\) and the low-resolution graphs \(\{\underline{G}\}\) in QM\(7_{\text{coarse}}\), we then investigate the ability of GNNs³ to consistently handle various scales by confronting models during inference with a resolution-scale different from the one they were trained on.

Table 1 collects results. Mean-absolute-errors (MAEs) during inference increase significantly, when going from a same-resolution setting to a cross-resolution setting. None of the standard architecture types are able to consistently handle more than one scale. Also employing multi-scale propagation schemes (SAG, SAG-M, UFGNet, Lanczos) does not result in consistently incorporating multiple scales: Corresponding cross-resolution MAEs are among the largest.

1.5pt

We can trace the inability of common models to generalize back to the difference in latent embeddings \(\{F\}\) and \(\{\underline{F}\}\), generated for original \(\{G\}\) and coarsified graphs \(\{\underline{G}\}\).

4pt

As evident from Table 2, latent embeddings of graphs describing the same object on varying resolutions are significantly different. This in turn explains the inability to generalize across scales. It should be noted that in practice, this failure cannot be overcome by augmenting the training set, as we have no way of generating faithful high-resolution descriptions given only lower resolution training data.

4 The Underlying Problem: Discontinuity↩︎

Our next goal then is to understand why the latent embeddings \(F, \underline{F}\) generated for original (\(G\)) and coarse (\(\underline{G}\)) graph descriptions of the same underlying object are so vastly different. Counter to conventional wisdom and consensus within the graph learning community we find below that

Hence, they may produce vastly different outputs \((F, \underline{F})\) for highly similar inputs \((G, \underline{G})\).

To see this discontinuity in action, note that in the coarse graphs \(\{\underline{G}\}\) of QM7\(_{\text{coarse}}\), hydrogen atoms have been fused into their nearest heavy atoms. Such a coarse graph hence can be viewed as the limit of a process in which hydrogen atoms are displaced from equilibrium toward their nearest heavy atom. The limit graph is obtained once they are fully merged into the corresponding heavy nodes.

If standard GNN architectures were continuous, the convergence of this graph modification process in Fig. 2 towards a limit graph should be reflected also in the latent space: Latent embeddings of modified graphs should converge to the latent embedding of the limit graph.

Figure 3 compares embeddings \(\{\underline{F}\}\) generated for coarsified graphs \(\{\underline{G}\}\), with embeddings \(\{ F_\omega\}\) of graphs \(\{G_\omega\}\) where hydrogen atoms have been moved to reduce the distance towards their nearest heavy atoms by a factor of \(\omega \geq 1\) (i.e. \(\text{dist}_\text{new} = \text{dist}_\text{equilib.}/\omega\)), but have not yet completely arrived at the positions of nearest heavy atoms.

From Fig. 3, we observe that latent embeddings do not converge (\(\|F_\omega-\underline{F}\|\nrightarrow 0\)). Since the sequence \(\{G_w\}_w\) of input graphs converges⁴ but GNN outputs evidently do not, they indeed cannot be considered continuous. Hence they may map similar graphs to vastly different latent embeddings.

5 The Obstruction: GNN Propagation↩︎

We can understand the underlying reason for this discontinuity by exemplarily investigating the prototypical graph neural network GCN [20]:⁵ Inside a GCN-layer, a node feature matrix \(X\in \mathbb{R}^{N \times F}\) (with number of nodes \(N\) and feature dimension \(F\)) is updated as \[\label{eq:gcn95prop} X\mapsto \hat{A} X W.\tag{3}\] Here \(W \in \mathbb{R}^{F \times F}\) mixes channels. Information flow over the graph is implemented via the so called renormalized adjacency matrix \(\hat{A} \in \mathbb{R}^{N \times N}\), with \(\hat{A}_{ij} \sim A_{ij}/\sqrt{d_id_j}\). When moving hydrogen (H) atoms towards heavy atoms (\(|\vec{x}_{\text{H}}-\vec{x}_{\text{heavy}}| \rightarrow 0\)), corresponding edge weights \(A_{\text{H},\text{heavy}} = 1\cdot Z_{\text{heavy}}\cdot|\vec{x}_{\text{H}}-\vec{x}_{\text{heavy}}|^{-1}\) of the original adjacency matrix \(A\) tend to infinity. Thus also node-degrees associated to heavy atoms tend to infinity. Since distances (and hence weights in the adjacency matrix \(A\)) between heavy atoms remain constant, the “renormalized" entries \(\hat{A}_{\text{heavy},\text{heavy}}\) in \(\hat{A}\), however, tend to zero instead.

Hence, as hydrogen atoms are moved out of equilibrium towards their final positions, the communication between heavy atoms in the modified graphs \(G_\omega\) becomes severely disrupted (\(\hat{A}_{\text{heavy},\text{heavy}} \rightarrow 0\)). Information is only propagated along a severely disconnected effective limit graph (dissected into distinct connected components (Fig 4 (a)); and not along the true coarse limit graph \(\underline{G}\) (Fig. 4 (b)).

As a consequence of the information flows over the graphs \(G_\omega, \underline{G}\) being vastly different, the latent embeddings \(F_\omega, \underline{F}\) that are being generated for the two graphs differ greatly.

At first glance, it may seem that the observations above apply only to the QM7 dataset. There edge weights scale inversely with distance, so sequences of graphs with diverging weights arise naturally. In Section 6 immediately below, however, we discuss that this is merely a proxy for a general fact in graph theory: the natural interpolation between an original graph \(G\) and its coarsified version \(\underline{G}\) proceeds by letting the edge weights within those clusters of \(G\) that collapse to the super-nodes of \(\underline{G}\) tend to infinity [31], [32]. Thus the observations above are indeed applicable to generic datasets.

6 Solution: Laplace Transform Propagation↩︎

Having established that standard GNNs are not continuous across scales due to the way they propagate information over a given graph \(G\), let us now work towards building architectures that instead are continuous.

6.1 From Graphs to Heat Kernels↩︎

To this end, let us begin by formalizing rigorously, in which sense the sequence of graphs \(G_\omega\) considered in Section 5 approaches the coarsified limit graph \(\underline{G}\):

6.1.0.1 Heat Kernel convergence:

When moving hydrogen atoms out of equilibrium, we are significantly increasing certain weights (\(A_{\text{H},\text{heavy}} \sim |\vec{x}_{\text{H}}-\vec{x}_{\text{heavy}}|^{-1} \sim \omega \rightarrow \infty\)). From a heat diffusion perspective, information in a graph equalizes much faster along edges with very large weights. In the limit where edge-weights within certain sub-graphs tend to infinity, information within these clusters equalizes immediately. Each such sub-graph thus effectively behaves as a single node in a coarse grained effective graph \(\underline{G}\).

To quantify this, we recall that the diffusion equation on a graph is given by \(dX(t)/dt = -L \cdot X(t)\) with solution \(X(t) = e^{-L t} \cdot X(0)\). As was established rigorously in [32] (cf. Appendix 17 for a discussion) the convergence \[\label{eq:orig95exp95conv95result95II} \|e^{-tL_\omega} - J^\uparrow e^{- t\underline{L}} J^\downarrow\| \rightarrow 0\tag{4}\] then holds for any fixed \(t>0\) as \(\omega \rightarrow \infty\). Here \(L_\omega, \underline{L}\) are the Laplacians of the respective graphs \(G_\omega, \underline{G}\). The matrices \(J^{\downarrow,\uparrow}\) appear naturally and allow to compare graphs of different sizes: Projection \(J^\downarrow\) assigns the average over strongly connected clusters to the super-node representing this cluster in \(\underline{G}\). Interpolation \(J^\uparrow\) assigns information at a super-node in \(\underline{G}\) to each node in the respective cluster in \(G\).

It is important to note that in order to detect the convergence \(G_\omega \rightarrow \underline{G}\) as in (4 ) one needs to make use of the (non-normalized) Laplacian \(L\). Using e.g. the normalized Laplacian \(\mathcal{L}\) would not lead to convergence (\(e^{-t\mathcal{L}_\omega} \nrightarrow J^\uparrow e^{-t\underline{\mathcal{L}}} J^\downarrow\)). Also directly comparing the Laplacians \(L, L_\omega\) would not detect the convergence: Clearly \(\|L_\omega\| \rightarrow \infty\), so that the sequence of Laplacians \(L_\omega\) alone does not converge.

6.1.0.2 Interpolating between Graphs on different scales:

We can use the convergence guarantee of (4 ) to interpolate ("in edge weight space") between an original graph \(G\) and its coarse grained version \(\underline{G}\) via a family of intermediate graphs \(\{G_\omega\}_\omega\) with ever increasing intra cluster weights.

As we let \(\omega \rightarrow \infty\), the heat kernel on \(G\) then more and more resembles the one on \(\underline{G}\). To visualize this fact, we exemplarily plot in Fig. 6 (c) the difference \(\| e^{-L_\omega t} - J^\uparrow e^{-\underline{L} t}J^\downarrow \|\) for the coarse graining setting of Figure 6 (a, b).

For fixed \(t > 0\) we see that \(\|e^{-tL_\omega} - J^\uparrow e^{- t\underline{L}} J^\downarrow\| \rightarrow 0\) as \(\omega\) increases. Additionally, the decay \(\|e^{-tL_\omega} - J^\uparrow e^{- t\underline{L}} J^\downarrow\| \rightarrow 0\) for increasing \(t\) is faster, the larger \(w\) is chosen. This is congruent with our intuition: The stronger two nodes are connected, the more they act as a single entity.

6.2 From Heat Kernels to Propagation Schemes↩︎

From (4 ) and Fig. 6 we conclude that applying the matrix \(e^{-tL_\omega}\) on \(G\) is more and more the same as projecting to \(\underline G\) via \(J^\downarrow\), applying \(e^{-t\underline{L}}\) there and interpolating back via \(J^\uparrow\).

For a network to be continuous across scales, we need exactly this property: We want the propagation scheme on \(G_\omega\) to more and more resemble the one on \(\underline{G}\) as we increase \(\omega \rightarrow \infty\). As we have seen in Section 4, the propagation rule \(X \mapsto \hat{A} X W\) of GCN does not achieve this, as it leads to disconnected limit graphs instead. However, propagating as \[X \mapsto e^{-tL} X W,\] does facilitate contact and similarity between information flows over \(G_\omega\) and \(\underline{G}\) as \(\omega \rightarrow \infty\).

More generally, suppose we have for each time \(t > 0\) individually that \(\|e^{-L_{\omega}t} - J^\uparrow e^{- t\underline{L}t} J^\downarrow\| < \delta\). If we build up the propagation matrix \(\psi(L_\omega)\) as a weighted sum of such diffusion flows \(e^{-tL_\omega}\) that have progressed to various times (\(\psi(L_\omega) \sim \sum_{k} a_k e^{-t_k L_\omega}\)) and the coefficients \(\{a_k\}_k\) are not too large, then we can estimate \(\| \psi(L_\omega) - J^\uparrow \psi(\underline{L}) J^\downarrow\| \leq \left( \sum_k |a_k|\right)\cdot \delta =: \epsilon\) by a triangle-inequality argument. Thus we can still guarantee that for large \(\omega\) the propagation implemented by the layer-wise update rule \[\label{eq:global95laplacian95proagation} X \mapsto \psi(L_\omega) XW\tag{5}\] over \(G_\omega\) is approximately the same as the effective propagation \(X \mapsto [J^\uparrow \psi(\underline{L}) J^\downarrow] XW\) over \(\underline{G}\).

Generalizing the weighted sum to an integral, we define:

Allowing generalized functions means we e.g. allow:

• Dirac distributions \(\hat{\psi}_{\delta_{s}}(t) := \delta(t-s)\); yielding back heat kernels \(\psi_s(L) = \int_{0}^\infty \delta(t-s) e^{-tL}dt = e^{-s L}\).

• Decaying functions \(\hat{\psi}_k := (-t)^{k-1}e^{-\lambda t}\) yielding \(k^{\text{th}}\)-powers of resolvents \(\psi_k(L) = [(\lambda I+L)^{-1}]^k\).

For such matrices \(\Psi(L)\) we note (details in Appendix 15.2):

7 The Result: Scale-Continuous GNNs↩︎

Lemma 1 above guarantees networks built on the propagation rule (5 ) are continuous across scales. We may think of this propagation rule as a scale continuous modification to the propagation rule (3 ) of GCN. GCN in turn is simultaneously both the simplest spectral- as well as message-passing based graph neural network. Our next goal is therefore to leverage the insights obtained from making GCN scale continuous to extend scale continuity to general spectral- and message-passing based GNNs.

7.1 Laplace Transform Propagation↩︎

7.1.0.1 Spectral Networks:

In order to turn spectral networks continuous across scales, we replace the polynomials \(T_\ell(\mathcal{L})\) of the normalized Laplacian \(\mathcal{L}\) used in their standard layer-wise update (2 ) with functions \(\psi_\ell(L)\) of the graph Laplacian \(L\) as defined in (6 ). The resulting layer wise update reads \[\begin{align} \label{eq:spectral95scale95ct} X^{(k+1)} = \sum_{\ell=0}^{K} \psi_\ell(L)\, X^{(k)}\, W_\ell. \end{align}\tag{7}\]

7.1.0.2 Message Passing Networks:

To determine which modifications are needed in order to turn message passing based networks scale continuous, we gain intuition from the the binary edge setting \((A_{uv} \in \{0,1\})\): There we may write the definition of the message function in (1 ) equivalently simply as \(m_{uv} = A_{uv} \cdot \phi(h_v^{(k)}, h_u^{(k)})\). As we saw in Sections 4 & 5, such adjacency-focused propagation schemes result in networks that are not continuous across scales. We hence make the replacement \(A \mapsto \psi(L)\) in our propagation scheme, and choose \(\operatorname{sum}\) as the aggregation function: \[\begin{align} m_{vu}^{(k)} &= [\psi(L)]_{vu} \cdot \phi(h_v^{(k)}, h_u^{(k)}) \tag{8}\\ h_v^{(k+1)} &= \sum_{u \in G} m_{vu} \tag{9}. \end{align}\] Additionally one initial pass \(h_0 = \Psi(L) X\) on the initial node features \(X\) is performed, before feeding them into our scale-continuous message passing framework (8 ) & (9 ).

7.1.0.3 Scalability:

At first glance, propagation along Laplace-transform matrices as in (6 ) may appear difficult to scale, since matrices such as \(\psi(L) = e^{-tL}\) are typically dense. A naïve implementation would thus require sending \(N^2\) messages over an \(N\)-node graph. This does however not pose a practical bottleneck: most entries are vanishingly small (\([\psi(L)]_{ij} \ll \delta\)), decaying exponentially with (graph- or resistance-) distance \(d(i,j)\) between nodes \(i\) and \(j\) [33]. Hence, propagation remains effectively sparse, and negligible edges can be safely discarded in practice. We explore this further in Appendix 18.

7.2 Theoretical Guarentees↩︎

Modifying spectral- and message-passing based graph neural networks exactly as outlined then indeed yields scale continuous networks, as we show next. Our main Theorem directly below establishes scale-continuity of the modified GNNs outlined in Section 7.1:

The proof of this Theorem is provided in Appendix 15.2. We note that whenever \(G_\omega\) converges to \(\underline{G}\) in the sense of heat kernels, by definition we have \(\|e^{-tL_\omega} - J^\uparrow e^{- t\underline{L}} J^\downarrow\| \rightarrow 0\) for all \(t>0\). Hence the integrand in Theorem 1 converges to zero. Hence also \(\|F_\omega - \underline{F}\|\) converges to zero.

Using the same proof strategy, one may also show:

The proof is provided in Appendix 14.2. Similar guarantees also hold at the node level, (cf. Appendix 14.1).

8 Numerical Investigation of Scale Continuity↩︎

8.0.0.1 Continuity across scales:

Theorem 1 implies that networks employing Laplace-transform propagation schemes are continuous as maps from the space of graphs into their latent spaces. To numerically verify this, we repeat the experiment of Section 4 for models belonging to this category (using resolvent and exponential matrices; cf. Definition 1 in a spectral and a spatial setting). As is evident from Fig. 7, latent embeddings generated by models employing Laplace transform propagation do indeed converge (\(\|F_\omega-\underline{F}\| \rightarrow 0\)). Hence for the models constructed in Section 7.1, we have indeed verified the desired continuity across scales.

8.0.0.2 Resulting Generalization Ability:

In Section 4 we had identified lack of continuity as the obstruction to generalizing across scales. As verified above, graph neural networks based on Laplace-transform propagation are continuous. Hence we expect them to map similar graphs to similar latent embeddings. To verify this, we here repeat the generalization experiment of Tables 1 & 2 in Section 3.

3pt

Comparing with standard GNNs we see that in cross-resolution settings the difference \(\|F - \underline{F}\|\) of latent embeddings generated by methods employing global Laplacian propagation schemes are lower than those of standard graph learning methods by factors of order \(10^2\) to \(10^4\), indicating that the methods developed in Section 7.1 indeed generalize.

The small-to-negligible variations in latent representations then translate to small-to-negligible variations in prediction performance: As we infer from Table 4 below, MAEs generated by Laplace-Transform based GNNs in the cross resolution setting are essentially the same as those corresponding to same-resolution settings. Hence GNNs based on Laplace Transform propagation schemes are indeed able to generalize across scales.

1pt

8.0.0.3 Scale Continuity at the Node Level:

We next consider popular citation networks (cf. Appendix 16.4 for details) where nodes correspond to scientific articles. Labels correspond to the academic discipline of the paper and an edge implies a citation. We then expand individual nodes into connected \(k\)-cliques (cf. Fig. 8). We might interpret this as dissecting each article into subsections, which reference each other. Both typical models and the Laplace-transform propagation based methods introduced in Section 7.1 were then trained on the same (\(k\)-fold expanded) train-set and asked to classify nodes in the (\(k\)-fold expanded) test-partition.

The classification accuracies of methods not employing Laplace-transform propagation decrease significantly with increasing clique size (cf. Fig. 8). We can understand the underlying reason for this using GCN as an Example (cf. Appendix 13 for discussions on other methods): Inside a GCN-layer, a node feature matrix \(X\) is updated as \(X\mapsto \hat{A} X W\), with the renormalized adjacency matrix \(\hat{A}\) given as \(\hat{A}_{ij} \sim A_{ij}/\sqrt{d_id_j}\). As the degree \(d_i\) of each node increases (linearly) with increasing clique-size \(k\), the message-strength \(\hat{A}_{ij}\) between the respective cliques decreases as \(\hat{A}_{ij} \sim 1/k\). Hence information propagation between the cliques becomes disrupted as \(k\) increases: GCNs propagation is more and more similar to a modified version of the original graph graph where edges between cliques are removed (\(G_2\) in Fig. 9 (a) bottom).

This is not the case for Laplace-transform propagation based networks (using either resolvent or exponential propagation matrices): As Corrolary 2 elucidates, similarity of graphs for such models is not determined through the renormalized adjacency as for GCN. Rather it is governed by the similarity the heat kernels corresponding to the graphs, (i.e. the size of \(\|e^{-t L_1} - e^{-t L_2}\|\) for any \(t > 0\)). As we can deduce from Fig. 9 (b) this quantity does not tend to zero with increasing clique size. Hence for the models of Section 7.1, propagation does not happen along a disconnected effective graph. Thus they are able to propagate information also between high connectivity areas and therefore are able to retain a high classification accuracy as the resolution scale is increased.

8.0.0.4 Scale continuity in the binary edge weight setting:

In Section 6.1 we argued that from a diffusion perspective, as edge-weights within certain sub-graphs of a graph \(G\) tend to infinity, each such sub-graph effectively behaves as a single node in a coarse grained effective graph \(\underline{G}\). This allows to continuously interpolate between \(G\) and its coarse grained version \(\underline{G}\) by scaling up edge weights within clusters. In order to be able to generalize across scales, we then modified standard GNNs to turn them continuous with respect to this limiting procedure in Section 7.1.

We next test continuity of GNNs under a related but distinct graph-limiting procedure: Conceptually we can think of increasing the weights of edges within clusters as a proxy for increasing the connectivity within these subgraphs. To increase connectivity, rather than increasing the weights of a fixed number of edges, we might instead also increase the number of edges within a cluster. To do this, we make use of stochastic block models (c.f. Appendix 16.3) (SBM) and vary the intra-cluster probability from \(p = 0\) to \(p = 1\).

We then compare latent embeddings \(F_p\) generated for graphs drawn from an SBM at intra-cluster connectivity \(p\) with the latent embeddings \(\underline{F}\) generated for a coarse grained version of this graph where clusters are aggregated to single nodes. As is evident from Fig. 11 below, latent embeddings generated by Laplace transform based models increasingly approximate the latent embedding \(\underline{F}\) as \(p \rightarrow 1\), while the same is not true for standard GNNs. In latent space, this continuity allows GNNs of Section 7.1 to embedd clique graphs close to graphs where each such clique is represented as a supernode (cf. Table [tab:sbm95lat95diff])).

8.0.0.5 Scale Continuity in the Continuum Limit:

Finally, we consider the setting where graphs discretize a continuous manifold. In this regime, we require that as we increase the resolution scale, model outputs converge to those that would be produced by a model operating directly on the underlying continuous manifold. In particular we want the correlation between two nodes representing (the same two) points on the manifold to stabilize, as the resolution is increased.

To show that this requirement is not fulfilled by standard GNNs, we feed standard models with a dirac impulse at a fixed point (located at \(\phi = 0\), \(\theta = 0\)) on the torus \(\mathbb{T} = [0,2\pi)^2\). We then investigate the output the GNNs generate at a test node, which is positioned at angles \(\phi = \pi/2\), \(\theta = 0\).

As evident from Fig. 13 (a), this requirement is not fulfilled for standard GNNs, as as neithers models’ response stabilizes: For most GNNs the impuls response drops to zero, while for Lanczos it fluctuates as the mesh resolution (i.e. number of nodes) is increased.

In contrast to that, the impulse response for the Laplace-transform based methods of Section 7.1 stays consistent as the mesh resolution is varied. To show that this persists also at graph level, we would ideally want to show that the graph level latent embeddings \(F_N\) generated for a mesh on \(N\) nodes converge to the latent embedding corresponding to the continuous manifold as \(N \rightarrow \infty\). In Appendix 16.5 we detail that we indeed have actual convergence \(F_N \rightarrow F\) towards a limit embedding \(F\) corresponding to the underlying continuous manifold. Since \(F\) is impossible to compute numerically, we verify the condition \(\|F_N - F_{2N}\| \xrightarrow{N \rightarrow \infty} 0\) (necessary for convergence) in Fig. 14 below instead.

9 Discussion↩︎

Our paper analyzed the discontinuity of existing GNNs across scales. We found the underlying obstruction to originate in commonly-used propagation schemes. We derived modifications to turn GNNs continuous, and showed that these modified models based on Laplace-Transform propagation can indeed consistently incorporate varying scales.

Acknowledgements↩︎

This work has received funding from the Swiss State Secretariat for Education, Research, and Innovation (SERI). C.K. thanks Madeline Navarro for helpful discussions during the rebuttal phase. M.B. is partially supported by the EPSRC Turing AI World-Leading Research Fellowship No. EP/X040062/1 and EPSRC AI Hub No. EP/Y028872/1.

Impact Statement↩︎

This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.

10 Discussion of Tasks and Graph Structures for which Scale Continuity is important↩︎

When GNNs are tasked to strictly distinguish non-isomorphic graphs the exact articulation of a graph matters (e.g. program graphs, electrical circuit graphs, etc.). In other settings, multiple graphs may describe the same underlying object (graphs discretizing the same manifold, sampled from the same distribution, or describing the same underlying noisy (social-, or natural-) science dataset. Here we want models to effectively embed the underlying object and not overfit on exact fine-print articulation of graphs describing said object.

A notion of distance respecting this intuition arises from comparing heat kernels as \(\|e^{-t L_1 - e^{-t L_2}}\|\). Using exponentials leads to an (exponential) suppression of large eigenvalues of the Laplacian \(L\). The suppressed information encoded into these large eigenvalues corresponds to fine structure details of the graph (c.f. e.g.[16]) which the model is supposed to discount. When comparing graphs on different numbers of nodes, we have to align them via projection \(J^\downarrow\) and interpolation \(J^\uparrow\) operators determining how to traverse between these graphs.

The quantity \(\|e^{-t L_1} - J^\uparrow e^{-t L_2} J_\downarrow \|\) (c.f. e.g. eq. (4 ) & Fig. 6) measures similarity of coarse structures in \(G_1, G_2\), after alignment. If this difference is small, coarse structures are both aligned and similar. Only then should models generate similar embeddings.

11 Discussion of Related Work↩︎

In this section we provide further details on existing approaches to estimating variatons in latent embeddings generated by graph neural networks:

11.0.0.1 Graphon Neural Networks and the Transferability of Graph Neural Networks [12]:

This seminal work explores the theoretical underpinnings of Graph Neural Networks (GNNs) in the context of graphons, a mathematical generalization of graphs to large-scale, continuous structures. The paper establishes a connection between GNNs and graphons, providing insights into the behavior of GNNs on large, dense graphs (\(|\mathcal{E}|\) is of \(\mathcal{O}(N^2)\), with \(N\) the number of nodes [15]) by modeling these graphs as graphons. This framework helps understand how GNNs operate in the limit of large graphs and illuminates their potential to generalize across different graph structures in this realm. A central focus of the paper is the transferability of GNNs—specifically, their ability to perform well on large graphs that may differ in size or topology from those seen during training. Transferability errors between graphs discretizing the same graphon are established to be of \(\mathcal{O}(N^{-\frac{1}{2}})\), with \(N\) the minimum number of nodes. Since this approach is aplicable only when graphs are large and \(|\mathcal{E}| = \mathcal{O}(N^2)\), it e.g. does not cover the setting of Table 1.

11.0.0.2 Transferability of Graph Neural Networks: an Extended Graphon Approach [13]:

This work is in spirit similar to [12] whose results it extends from considering the adjacency matrix as the graph shift operator to more general graph shift operators and from considering only polynomial filters to allowing for general continuous filter functions. Similar limitations as for [12] apply.

11.0.0.3 Limits, approximation and size transferability for GNNs on sparse graphs via graphops [15]:

In contrast to approaches using graphons, which focus on large dense graphs, this paper instead focuses on transferability on sparse graphs (\(|\mathcal{E}| = \mathcal{O}(N)\)) making use of the concept of Graphops, mathematical operators that can be used to model how GNNs behave on large sparse graphs. The authors also explore how GNNs can transfer learned representations from smaller, sparse graphs to larger ones, and vice versa. By leveraging the Graphop framework, the paper formalizes conditions for successful transferability between graphs of varying sizes. Due to the necessity for graphs to be large (with edges \(|\mathcal{E}| = \mathcal{O}(N)\)) it also does not apply e.g. to the setting considered in Section 3.

11.0.0.4 On Local Distributions in Graph Signal Processing [14]:

This work is rooted in the field of graph signal procesing (GSP) and puts a particular emphasis on the transferability of GSP techniques across different graph structures. The paper focuses on the concept of graphings, which are a probabilistic framework for representing large sparse graphs and their underlying structures. By modeling large graphs through graphings, the authors provide a framework that makes it possible to generalize local distributions and signal processing tasks across different graphs. This work similarly is only applicable to large sparse graphs.

11.0.0.5 Graph Convolutional Neural Networks via Scattering [34]

This work provides a different perspective on spectral GNNs by connecting them to scattering transforms [35]–[37], a concept from signal processing. The authors demonstrate that spectral GNNs can be interpreted as a discrete graph counterpart of scattering transforms, which involve multi-scale wavelet-like operations that capture hierarchical information across different levels of graph structure. A key focus of the paper is the stability of networks when viewed through the scattering framework. The authors argue that scattering transforms offer a more stable approach to graph signal processing compared to traditional networks, especially in the presence of noisy or incomplete graph data. Derived single filter transferability results depend on spectral properties of the utilized Laplacians on the respective graphs. The focus of this work is on the setting where \(\|L - \tilde{L}\|\) is small. This is distinct from the setting of similar heat kernels, which is more general.

11.0.0.6 Stability to Deformations of Manifold Filters and Manifold Neural Networks [38]

:

This work explores the theoretical foundation of manifold filters and manifold neural networks (MNNs). Similarly to the filters analyzed in the present work, manifold filters are defined in terms of Laplace transforms. By framing spectral graph neural networks as discrete approximations of MNNs, the authors analyze conditions under which MNNs remain stable under smooth deformations of the manifold. Stability is shown to depend on specific spectral properties of the filter functions, including Lipschitz continuity and integral Lipschitz continuity, which control the trade-off between robustness and frequency discriminability. The paper establishes that filters meeting these conditions can generalize effectively to new manifolds by adapting to changes in the Laplace-Beltrami operator’s spectrum. More techically, filters are bounded as \(|\psi(L) - \psi(\tilde{L})| \leq K \|L - \tilde{L}\|\). In Theorem 2 absolute perturbations are considered (\(\tilde{L} = L + A\)), in Theorem 3 relative perturbations are considered (\(\tilde{L} = L + EL\)). In both cases the conditions on spectrum and filter functions stem from the fact that Lipschitz-ness does not directly translate to operator Lipschitz-ness when measured in spectral norm (see e.g. [39] for a discussion). In contrast, our work does not primarily focus on (large) graphs approximating a common manifold. We only consider this setting as one example. In comparison to [38], the results we develop for this particular example (c.f. Appendix 16.5) need no assumptions on band-limitedness. Additionally our results also apply to message passing based networks, in contrast to those of [38] which apply to spectral networks.

11.0.0.7 Geometric Graph Filters and Neural Networks: Limit Properties and Discriminability Trade-offs [40]:

Here instead of measuring the linear norm difference \(\|LP - \mathcal{L}P\|\) between a graph Laplacian \(L\) and a manifold Laplacian \(\mathcal{L}\) (which generically would be infinite as \(\mathcal{L}\) is an unbounded operator), the difference of the action of these operators on eigenfunctions (\(\| LP\phi - \mathcal{L}P\phi \|\)) is considered. After a triangle inequality argument, one term that has to be bounded in order to bound the difference in filter outputs is \(\|\phi_i^n - \phi_i\|\) of the \(i^{th}\) eigenfunction and eigenvector respectively. The fidelity of this approximation depends on spectral separation properties (cf. Theorem 4 ibid.), which hence leads to the requirement that the spectrum be \(\alpha\)-separated. This requirement can thus be considered an artifact of considering the linear approximation \(\|\phi_i^n - \phi_i\|\) for each eigenfunction. In contrast, in our approach when applied to the specific example of discretized manifolds (cf. Appendix 16.5) the notion of approximation of the Laplacian on the underlying manifold is different. We bound the quantity \(\| J^\uparrow e^{-tL} J^\downarrow- e^{-t\Delta} \|\) instead. Hence we do not need to bound differences between individual eigenfunctions and eigenvectors and hence avoid dependencies on spectral separations. Additionally, our results also cover message passing based networks, and are not restricted to spectral networks.

11.0.0.8 Transferability of Spectral Graph Convolutional Neural Networks [11]:

As one of the earliest works challenging the then prevailing belief that spectral methods are not transferable, this work was among the first to present theoretical proofs and experimental evidence to demonstrate that these methods can generalize effectively under certain conditions. The key contribution is a theoretical framework in which transferability depends on how well graphs approximate a shared underlying continuous domain, such as a topological space or metric-measure space. The study also develops sufficient conditions for achieving low transferability errors, demonstrating that spectral gnns can perform consistently across graphs with varying sizes, topologies, and dimensions, provided the graphs discretize the same continuous domain. As in our work, filters here are only required to be bounded and Lipschitz continuous (cf. Theorem 17 ibid.). However, signals are assumed to be bandlimted. In the specific example setting of graphs discretizing manifolds, we avoid Levie’s growth of the stability constant with the number of considered eigenvalues (cf. the discussion towards the end of page 12 ibid.) by avoiding approximations of individual eigenfunctions and instead approximating the bounded operator \(e^{-t\Delta}\) directly. Additionally, our results also cover message passing based networks, and are applicable also outside the setting of graphs discretizing ambient spaces.

11.0.0.9 Diffusion Scattering Transforms on Graphs [41]:

This work emphasizes the stability of scattering-based representations against perturbations in graph topology and reindexing. The framework introduces diffusion scattering transforms that leverage diffusion operators to capture multi-scale hierarchical features of graph signals. The authors focus on ensuring that the transforms are robust to changes in graph structure, such as modifications to edge weights or topology. Stability is achieved through the use of diffusion wavelets, which provide a principled way to construct graph filters that are invariant to local perturbations while retaining sensitivity to meaningful global graph features. The stability analysis demonstrates that the scattering transform bounds the impact of graph perturbations in terms of the changes they induce in the graph Laplacian’s spectrum, ensuring reliable performance across varied graph inputs. Results derived in this work do not apply to our settings, as detecting graph variations via the difference \(\|\mathcal{L} - \tilde{\mathcal{L}}\|\) as considered in [41] is insensitive to the notion of variation and convergence we consider here (c.f. Section 6.1). Additionally the use of the normalized Laplacian in [41] precludes detection of the convergence behaviour of Section 6.1.

11.0.0.10 Stability Properties of Graph Neural Networks [9]:

This paper investigates the stability properties of GNNs to perturbations in the underlying graph structure. The authors analyze how small changes in graph topology— such as modifications to edge weights, addition or deletion of edges, or reindexing of nodes—affect the outputs of spectral GNNs. The paper develops a rigorous mathematical framework to assess the stability of GNNs using tools from spectral graph theory. It establishes that GNNs are stable to localized perturbations in the graph topology, with the degree of stability depending on the spectral properties of the graph filters used within the network. Specifically, it is shown that GNNs exhibit a trade-off between stability and discriminability: filters that are more stable to perturbations may sacrifice sensitivity to high-frequency information, which can limit their ability to differentiate fine-grained graph structures. Here as well, Lipschitz type arguments are being used (see e.g. the assumptions of Theorem 1) to establish single filter transferability. Since scalar Lipschitzness does not translate to operator Lipschitzness under spectral norm, additional restrictions on spectrum and filter functions need to be hence imposed. Following this, the authors highlight the importance of filter design in achieving a balance between robustness and expressivity. Filters that adhere to conditions such as Lipschitz continuity or integral Lipschitz continuity are particularly effective in maintaining stability while preserving key graph features. Results derived in this work do not apply to our settings, as detecting graph variations via the difference \(\|L - \tilde{L}\|\) is insensitive to the notion of variation and convergence we consider here (c.f. Section 6.1).

11.0.0.11 Limitless transferability for graph convolutional Networks [42]:

This work studies stability- and transferability properties of spectral graph neural networks, with a particular focus on directed graphs. In spirit, it is the closest to our work here, as one of the main class of filters it investigates is the class of resolvent based spectral filter functions which constitute an example of the more general class of Laplace transform filters considered in this present work. Our work here is however far more general, considering general Laplace transform based propagations schemes, including not only spectral-, but also message passing based networks. Additionally we develop a rich and detailed theory and provide numerical investigations and experimental evaluations in diversex settings.

11.0.0.12 Previous versions of this present work:

Earlier versions and preliminary aspects of this research were presented at various non-archival workshops [43]–[48]. The various aspects discussed there are subsumed and extended by this present work.

12 Laplace Transform Matrices↩︎

In this section we provide an overview of the concept of Laplace transforms. We begin with a recapitulation of complex measures.

12.1 Complex measures on \(\mathbb{R}_{\geq 0}\) and their Theory of Integration↩︎

As reference for this section [49] might serve.

In mathematics, a measure is a formal generalization of concepts such as length, area and volume. We are interested in assigning a generalized notion of length (or mass) to subsets of the real half-line \[\mathbb{R}_{\geq 0} = [0,\infty).\] The set will turn out to be a so called \(\sigma\)-Algebra; i.e. a set \(\Sigma\) of sets for which

• \(\emptyset, \mathbb{R}_{\geq 0} \in \Sigma\)

• \(A,B \in \sigma \Rightarrow A\cap B \in \Sigma\)

• \(A,B\in \Sigma \Rightarrow A \setminus B \in \Sigma\)

• \(A,B\in \Sigma \Rightarrow A \cup B \in \Sigma\).

We now take \(\Sigma_{\mathbb{R}_{\geq 0}}\) to be the smallest such set of sets \(\Sigma\) that contains all open intervals.

A complex measure then is a set-function that assigns to each set in \(\Sigma_{\mathbb{R}_{\geq 0}}\) a complex number in a certain way:

Let us provide some examples:

Any given measure on \(\mathbb{R}_{\geq 0}\) defines a unique way of integrating (known as Lebesgue integration) a function \(f\) defined on \(\mathbb{R}_{\geq 0}\). This proceeds by approximating any function \(f\) via a weighted sequence of indicator functions (with \(A \in \Sigma_{\mathbb{R}_{\geq 0}}\) a set)

\[\chi_A(t) = \begin{cases} 1 &; t \in A\\ 0 &; t \notin A \end{cases}.\]

as \[f(t) \approx f_n(t) := \sum_k a^n_k \chi_{A_k}(t).\] with \(a_k \in \mathbb{C}\). For these functions, one then sets \[\int_{\mathbb{R}_{\geq 0}} f_n d\mu \equiv \sum_k a^n_k\cdot \mu(A_k).\] Since we have \(\lim_{n \rightarrow \infty} f_n = f\), one then simply sets \[\int_{\mathbb{R}_{\geq 0}} f d\mu \equiv \lim_{n\rightarrow\infty} \int_{\mathbb{R}_{\geq 0}} f_n d\mu.\]

12.2 Laplace Transforms↩︎

We say complex valued measure \(\mu\) is finite if we have \[\int_{\mathbb{R}_{\geq 0}} d|\mu|(t) < \infty.\] Here the measure \(|\mu|\) arises from the original measure \(\mu\) via \[|\mu|((a,b)) \equiv |\mu((a,b))|.\] For any such finite measure \(\mu\) we may define its Laplace transform as \[\psi_\mu(z) := \int_{\mathbb{R}_{\geq 0}} e^{-tz} d\mu(t).\] This function \(f_\mu\) is well defined for \(z\) in the right hemisphere \[\mathbb{C}_R := \{z \in \mathbb{C} : \text{Re}(z) \geq 0\}.\] of the complex plane \(\mathbb{C}\), since there we have \[\begin{align} |\psi_\mu(z)| & = \left|\int_{\mathbb{R}_{\geq 0}} e^{-tz} d\mu(t)\right|\\ & \leq \int_{\mathbb{R}_{\geq 0}} |e^{-tz}|d|\mu|(t)\\ & \leq \int_{\mathbb{R}_{\geq 0}} d|\mu|(t) < \infty. \end{align}\]

More specifically, if the integrable function is given as \(\hat{\psi}_k := (-t)^{k-1}e^{-\lambda t}\) (with \(\text{Re}(\lambda) > 0\)), then \(\psi_k(z) = (z+\lambda)^{-k}\):

Using the function \(\psi_k\) of the examples above, a wide class of functions may be parametrized

13 Effective Propagation Schemes↩︎

For definiteness, we here discuss limit-propagation schemes in the setting where edge-weights are large. The discussion for high-connectivity in the Sense of large cliques proceeds analogously.

In this section, we then take up again the setting of Section 5. We reformulate this setting here in a slightly modified language, that is more adapted to discussing effective propagation schemes of standard architectures:
We partition edges on a weighted graph \(G\), into two disjoint sets \(\mathcal{E} = \mathcal{E}_{\text{reg.}} \dot{\cup} \mathcal{E}_{\text{high}}\), where the set of edges with large weights is given by: \[\mathcal{E}_{\text{high}} := \{(i,j)\in \mathcal{E}: w_{ij} \geq S_{\text{high}}\}\] and the set with small weights is given by: \[\mathcal{E}_{\text{reg.}} := \{(i,j)\in \mathcal{E}: w_{ij} \leq S_{\text{reg.}}\}\] for weight scales \(S_\text{high} > S_\text{reg.} > 0\). Without loss of generality, assume \(S_\text{reg.}\) to be as low as possible (i.e. \(S_\text{reg.} = \max_{(i,j)\in \mathcal{E}_\text{reg.}} w_{ij}\)) and \(S_\text{high}\) to be as high as possible (i.e. \(S_\text{large} = \min_{(i,j)\in \mathcal{E}_\text{high}}\)) and no weights in between the scales.

This decomposition induces two graph structures corresponding to the disjoint edge sets on the node set \(\mathcal{G}\): We set \(G_{\text{reg.}} := (\mathcal{G},\mathcal{E}_{\text{reg.}})\) and \(G_{\text{high}}:=(\mathcal{G},\mathcal{E}_{\text{high}})\) c.f. Fig. 15).
We also introduce the set of edges \(\mathcal{E}_{\text{reg., exclusive}} := \{(i,j)\in \mathcal{E}_{\text{reg.}}| \;\forall k \in \mathcal{G}:\;(i,k) \notin \mathcal{E}_{\text{high}}\;\&\;(k,j) \notin \mathcal{E}_{\text{high}} \}\) connecting nodes that do not have an incident edge in \(\mathcal{E}_{\text{high}}\). A corresponding example-graph \(G_{\text{reg., exclusive}}\) is depicted in Fig. 15 (d).

We are now interested in the behaviour of graph convolution schemes if the scales are well separated: \[S_{\text{high}} \gg S_{\text{reg.}}\]

13.1 Spectral Convolutional Filters↩︎

We first discuss resulting limit-propagation schemes for spectral convolutional networks. Such networks implement convolutional filters as a mapping \[x \longmapsto g_\theta(T)x\] for a node feature \(x\), a learnable function \(g_\theta\) and a graph shift operator \(T\).

13.1.1 Need for Normalization↩︎

The graph shift operator \(T\) facilitating the graph convolutions needs to be normalized for established spectral graph convolutional architectures:

For [50], this e.g. arises as a necessity for convergence of the proposed implementation scheme for the rational filters introduced there (c.f. eq. (10) in [50]).

The work [23] needs its graph shift operator to be normalized, as it approximates generic filters via a Chebyshev expansion. As argued in [23], such Chebyshev polynomials form an orthogonal basis for the space \(L^2([-1,1], dx/\sqrt{1-x^2})\). Hence, the spectrum of the operator \(T\) to which the (approximated and learned) function \(g_\theta\) is applied needs to be contained in the interval \([-1,1]\).

In [20], it has been noted that for the architecture proposed there, choosing \(T\) to have eigenvalues in the range \([0,2]\) (as opposed to the normalized ranges \([0,1]\) or \([-1,1]\)) has the potential to lead to vanishing- or exploding gradients as well as numerical instabilities. To alleviate this, [20] introduces a "renormalization trick" (c.f. Section 2.2. of [20] to produce a normalized graph shift operator on which the network is then based.

We can understand the relationship between normalization of graph shift operator \(T\) and the stability of corresponding convolutional filters explicitly: Assume that we have \[\|T\| \gg 1.\] This might e.g. happen when basing networks on the un-normalized graph Laplacian \(\Delta\) or the weight-matrix \(W\) if edge weights are potentially large (such as in the setting \(S_{\text{high}} \gg S_{\text{reg.}}\) that we are considering).

By the spectral mapping theorem (see e.g. [51]), we have \[\label{specctral95mapping} \sigma\left( g_\theta(T)\right) = \left\{g_\theta(\lambda) : \lambda \in \sigma(T)\right\},\tag{10}\] with \(\sigma(T)\) denoting the spectrum (i.e. the set of eigenvalues) of \(T\). For the largest (in absolute value) eigenvalue \(\lambda_{\max}\) of \(T\), we have \[\label{sp95norm95a95l95ev} |\lambda_{\max}| = \|T\|.\tag{11}\] Since learned functions are implemented directly as a polynomial (as e.g. in [23]), we have \[\lim\limits_{\lambda \rightarrow \pm \infty} |g_\theta(\lambda)| = \infty.\] Thus in view of (10 ) and (11 ) we have for \(\|T\|\) sufficiently large, that \[\|g_\theta(T)\| = |g_\theta(\pm\|T\|)|\] with the sign \(\pm\) determined by \(\lambda_{\max}\gtrless 0\). Since non-constant polynomials behave at least linearly for large inputs, there is a constant \(C>0\) such that \[C\cdot\|T\| \leq\|g_\theta(T)\|\] for all sufficiently large \(\|T\|\). We thus have the estimate \[\|x\|\cdot C\cdot \|T\| \leq \|g_\theta(T)x\|\] for at least one input signal \(x\) (more precisely all \(x\) in the eigen-space corresponding to the largest (in absolute value) eigenvalue \(\lambda_{\max}\)). Thus if \(T\) is not normalized (i.e. \(\|T\|\) is not sufficiently bounded), the norm of (hidden) features might increase drastically when moving from one (hidden) layer to the next. This behaviour persists for all input signals \(x\) have components in eigenspaces corresponding to large (in absolute value) eigenvalues of \(T\).

13.1.2 Spectral Normalizations↩︎

with \(\Delta\) the un-normalized graph Laplacian. In the setting \(S_{\text{high}} \gg S_{\text{reg.}}\) we are considering, this leads to an effective feature propagation along \(G_{\text{high}}\) (c.f. also Fig. 17) only, as Theorem 14 below establishes:

It remains to state and establish Lemma 15 referenced at the end of the proof of Theorem 14:

13.1.3 Symmetric Normalizations↩︎

Apart from networks that make use of the symmetrically normalized graph Laplacian \(\mathscr{L}\), some methods, such as most notably [20], instead base their filters on the operator \[T = \tilde{D}^{-\frac{1}{2}}\tilde{W}\tilde{D}^{-\frac{1}{2}},\] with \[\tilde{W} = (W+Id)\] and \[\tilde{D} = D + Id.\] In analogy to Theorem 16, we here establish the limit propagation scheme determined by such operators:

The difference to the result of Theorem 16 is thus that applicability of the limit propagation scheme of Fig. 19 for the GCN [20] is not only contingent upon \(S_{\text{high}}\gg S_{\text{reg.}}\) but also \(S_{\text{high}}\gg1\).

13.2 Spatial Convolutional Filters↩︎

Apart from spectral methods, there of course also exist methods that purely operate in the spatial domain of the graph. Such methods most often fall into the paradigm of message passing neural networks (MPNNs) [3], [53]: With \(X^{\ell}_i\in\mathbb{R}^F\) denoting the features of node \(i\) in layer \(\ell\) and \(w_{ij}\) denoting edge features, a message passing neural network may be described by the update rule (c.f. [3]) \[\label{mpnn95def95eq} X^{\ell+1}_i = \gamma\left(X^{\ell}_i, \coprod\limits_{j \in \mathcal{N}(i)} \phi\left(X_i^\ell, X_j^\ell, w_{ij}\right)\right).\tag{24}\] Here \(\mathcal{N}(i)\) denotes the neighbourhood of node \(i\), \(\coprod\) denotes a differentiable and permutation invariant function (typically "sum", "mean" or "max") while \(\gamma\) and \(\phi\) denote differentiable functions such as multi-layer-perceptrons (MLPs) which might not be the same in each layer. [53].

Before we discuss corresponding limit-propagation schemes, we first establish that MPNNs are not able to reproduce the limit propagation scheme ofFigure 4 (b) and are thus not stable to scale transitions and topological perturbations.

13.2.1 Scale-Sensitivity of Message Passing Neural Networks↩︎

Here we establish that message passing networks (as defined in (24 ) above) are unable to emulate a limit propagation scheme similar to the one in Figure 4 (b). Hence such architectures are also not stable to scale-changing topological perturbations such as coarse-graining procedures.

Given states \(\{X_1^\ell, X_2^\ell, X_3^\ell\}\) in layer \(\ell\), a limit propagation scheme as in Figure 4 (b) would require the updated feature vector of node \(3\) to be given by \[X_{3, \text{desired}}^{\ell +1} := \gamma\left(X^{\ell}_3, \phi\left(X_3^\ell, \frac{ X_1^\ell+X_2^\ell}{2}, (w_{31}+w_{32})\right) \right)\]

However, the actual updated feature at node \(3\) is given as (c.f. (24 )): \[\label{node95three95no95scale95dep95rem} X_{3, \text{actual}}^{\ell +1} := \gamma\left(X^{\ell}_3, \phi\left(X_3^\ell, X_1^\ell, w_{31}\right) \coprod \phi\left(X_3^\ell, X_2^\ell, w_{32}\right)\right)\tag{25}\] Since there is no dependence on \(S_{\text{high}}\) in equation (25 ) – which defines \(X_{3, \text{actual}}^{\ell +1}\) – the desired propagation scheme can not arise, unless it is paradoxically already present at all scales \(S_{\text{high}}\). If it is present at all scales, there is however only propagation along edges in \(\underline{G}\), even if \(S_{\text{high}}\approx S_{\text{reg.}}\), which would imply that the message passing network would not respect the graph structure of \(G\). Hence \(X_{3, \text{actual}}^{\ell +1} \nrightarrow X_{3, \text{desired}}^{\ell +1}\) does not converge as \(S_{\text{high}}\) increases.

13.2.2 Limit Propagation Schemes↩︎

The number of possible choices of message functions \(\phi\), aggregation functions \(\coprod\) and update functions \(\gamma\) is clearly endless. Here we shall exemplarily discuss limit propagation schemes for two popular architectures: We first discuss the most general case where the message function \(\phi\) is given as a learnable perceptron. Subsequently we assume that node features are updated with an attention-type mechanism.

13.2.2.1 Generic message functions:

We first consider the possibility that the message function \(\phi\) in (25 ) is implemented via an MLP using ReLU-activations: Assuming (for simplicity in notation) a one-hidden-layer MLP mapping features \(X_i^\ell \in \mathbb{R}^{F_\ell}\) to features \(X_i^{\ell+1} \in \mathbb{R}^{F_\ell+1}\) we have \[\phi(X^\ell_i,X^\ell_j,w_{ij}) = \text{ReLU}\left(W^\ell_1\cdot X^\ell_i+W^\ell_2\cdot X^\ell_2+ W^\ell_3\cdot w_{ij} +B^\ell \right)\] with bias term \(B^{\ell+1} \in \mathbb{R}^{F_{\ell+1}}\) and weight matrices \(W^{\ell+1}_1,W^{\ell+1}_2 \in \mathbb{R}^{F_{\ell+1}\times F_{\ell}}\) and \(W^\ell_3 \in \mathbb{R}^{F_{\ell+1}}\).

We will assume that the weight-vecor \(W_3^{\ell+1}\) has no-nonzero entries. This is not a severe limitation experimentally and in fact generically justified: The complementary event of at-least one entry of \(W_3\) being assigned precisely zero during training has probability weight zero (assuming an absolutely continuous probability distribtuion according to which weights are learned).

Let us now assume that the edge \((ij)\) belongs to \(\mathcal{E}_{\text{high}}\) and the corresponding weight \(w_{ij}\) is large (\(w_{ij} \gg 1\)). The behaviour of entries \(\phi(X^\ell_i,X^\ell_j,w_{ij})_a\) of the message \(\phi(X^\ell_i,X^\ell_j,w_{ij}) \in \mathbb{R}^{F_{\ell+1}}\) is then determined by the sign of the corresponding entry \(\left(W_3^\ell\right)_a\) of the weight vector \(W_3^\ell\in \mathbb{R}^{F_{\ell+1}}\):

If we have \(\left(W_3^\ell\right)_a<0\), then \(\phi(X^\ell_i,X^\ell_j,w_{ij})_a\) approaches zero for larger edge-weights \(w_{ij}\): \[\label{uninformative95message} \lim\limits_{w_{ij}\rightarrow \infty} \phi(X^\ell_i,X^\ell_j,w_{ij})_a = 0\tag{26}\]

If we have \(\left(W_3^\ell\right)_a>0\), then \(\phi(X^\ell_i,X^\ell_j,w_{ij})_a\) increasingly diverges for larger edge-weights \(w_{ij}\): \[\label{dangerous95message} \lim\limits_{w_{ij}\rightarrow \infty} \phi(X^\ell_i,X^\ell_j,w_{ij})_a = \infty\tag{27}\] For either choice of aggregation function \(\coprod\) in (24 ) among "max", "sum" or "mean" the behaviour in (27 ) leads to unstable networks if the update function \(\gamma\) is also given as an MLP with ReLU activations. Apart from instabilities, we also make the following observation: If \(S_{\text{high}} \gg S_{\text{reg.}}\), then by (27 ) and continuity of \(\phi\) we can conclude that components \(\phi(X^\ell_i,X^\ell_j,w_{ij})_a\) of messages propagated along \(\mathcal{E}_{\text{high}}\) for which \(\left(W_3^\ell\right)_a>0\) dominate over messages propagated along edges in \(\mathcal{E}_{\text{reg.}}\). By (26 ), the former clearly also dominate over components \(\phi(X^\ell_i,X^\ell_j,w_{ij})_a\) of messages propagated along \(\mathcal{E}_{\text{high}}\) for which \(\left(W_3^\ell\right)_a<0\). This behaviour is irrespective of whether "max", "sum" or "mean" aggregations are employed. Hence the limit propagation scheme essentially only takes into account message channels \(\phi(X^\ell_i,X^\ell_j,w_{ij})_a\) for which \((ij)\in \mathcal{E}_{\text{high}}\) and \(\left(W_3^\ell\right)_a>0\).

Similar considerations apply, if non-linearities are chosen as leaky ReLU. If instead of ReLU activations a sigmoid-nonlinearity \(\sigma\) like \(\tanh\) is employed, messages propagated along \(\mathcal{E}_{\text{large}}\) become increasingly uninformative, since they are progressively more independent of features \(X_i^\ell\) and weights \(w_{ij}\). Indeed, for sigmoid activations, the limits (26 ) and (27 ) are given as follows:

If we have \(\left(W_3^\ell\right)_a<0\), then we have for larger edge-weights \(w_{ij}\) that \[\label{uninformative95message95sigma} \lim\limits_{w_{ij}\rightarrow \infty} \phi(X^\ell_i,X^\ell_j,w_{ij})_a = \lim\limits_{y \rightarrow -\infty} \sigma(y).\tag{28}\]

If we have \(\left(W_3^\ell\right)_a>0\), then \[\label{dangerous95message95sigma} \lim\limits_{w_{ij}\rightarrow \infty} \phi(X^\ell_i,X^\ell_j,w_{ij})_a = \lim\limits_{y \rightarrow \infty} \sigma(y).\tag{29}\] In both cases, the messages \(\phi(X^\ell_i,X^\ell_j,w_{ij})\) propagated along \(\mathcal{E}_{\text{large}}\) become increasingly constant as the scale \(S_{\text{high}}\) increases.

13.2.2.2 Attention based messages:

Apart from general learnable message functions as above, we here also discuss an approach where edge weights are re-learned in an attention based manner. For this we modify the method [19] to include edge weights. The resulting propagation scheme – with a single attention head for simplicity and a non-linearity \(\rho\) – is given as \[X_i^{\ell +1} = \rho\left(\sum\limits_{j \in \mathcal{N}(i)}\alpha_{ij}(WX_j^{\ell +1})\right).\] Here we have \(W \in \mathbb{R}^{F_{\ell+1}\times F_\ell}\) and \[\label{att95weight95mat} \alpha_{ij}=\frac{\exp\left(\text{LeakyRelu}\left(\vec{a}^\top\left[WX_i^\ell\mathbin{\|}WX_j^\ell\mathbin{\|}w_{ij}\right]\right)\right)}{\sum\limits_{k \in \mathcal{N}(i)}\exp\left(\text{LeakyRelu}\left(\vec{a}^\top\left[WX_i^\ell\mathbin{\|}WX_k^\ell\mathbin{\|}w_{ik}\right]\right)\right)},\tag{30}\] with \(\mathbin{\|}\) denoting concatenation. The weight vector \(\vec{a}\in \mathbb{R}^{2F_{\ell+1}+1}\) is assumed to have a non zero entry in its last component. Otherwise, this attention mechanism would correspond to the one proposed in [19], which does not take into account edge weights. Let us denote this entry of \(\vec{a}\) ()determining attention on the weight \(w_{ij}\)) by \(a_w\).

If \(a_w <0\), we have for \((i,j)\in \mathcal{E}_{\text{high}}\) that \[\exp\left(\text{LeakyRelu}\left(\vec{a}^\top\left[WX_i^\ell\mathbin{\|}WX_j^\ell\mathbin{\|}w_{ij}\right]\right) \right) \longrightarrow 0\] as the weight \(w_{ij}\) increases. Thus propagation along edges in \(\mathcal{E}_{\text{high}}\) is essentially suppressed in this case.

If \(a_w >0\), we have for \((i,j)\in \mathcal{E}_{\text{high}}\) that \[\exp\left(\text{LeakyRelu}\left(\vec{a}^\top\left[WX_i^\ell\mathbin{\|}WX_j^\ell\mathbin{\|}w_{ij}\right]\right) \right) \longrightarrow \infty\] as the weight \(w_{ij}\) increases. Thus for edges \((i,j) \in \mathcal{E}_{\text{reg.}}\) (i.e. those that are not in \(\mathcal{E}_{\text{high}}\)), we have \[\alpha_{ij} \rightarrow 0,\] since the denominator in (30 ) diverges. Hence in this case, propagation along \(\mathcal{E}_{\text{reg.}}\) is essentially suppressed and features are effectively only propagated along \(\mathcal{E}_{\text{high}}\).

14 Stability of latent representations generated by Laplace-transform based GNNs for graphs on the same node set↩︎

14.1 Node Level Stability↩︎

In this Section, we are concerned with two graphs \(G_1\) and \(G_2\) that share a common node set. We assume that their heat kernels are similar in the sense of \[\begin{align} \|e^{-tL_1} - e^{-t L_2}\| \approx 0. \end{align}\] We are then interested in bounding the variation in latent representations generated by a GNN \(\Phi_{\mathscr{W},\mathscr{B},\Psi}\) (as defined in Section 7) with weights \(\mathscr{W}\), (potential) biases, and Laplace transform propagation matrices \(\Psi\).

14.1.1 Spectral Laplace-transform methods↩︎

We first consider spectral Laplace Transform based spectral GNNs:

Next we note that we can estimate \(\delta\) in terms of \(\max_i \{\int_0^\infty \hat{\psi}_i(t) \|e^{-tL_1} - e^{t L_2}\| dt \}\).

Combining Lemma 2 with Theorem 18 then establishes a bound on node-wise latent embeddings in terms of the heat kernel distance \(\left\| e^{-tL_1} - e^{-tL_2}\right\|\) for spectral Laplace transform based methods.

14.1.2 Message passing based Laplace Transform Methods:↩︎

We next establish similar results for Message passing networks. To this end, we first introduce some notation. We continue to denote the node-feature matrix by \(X\). As discussed in (8 )&(9 ), we have as layer-wise update the implemented the scheme \[\begin{align} m_{vu}^{(\ell)} &= [\psi(L)]_{vu} \cdot \phi(h_v^{(k)}, h_u^{(k)}) \tag{32}\\ h_v^{(\ell+1)} &= \sum_{u \in G} m_{vu} \tag{33}. \end{align}\]

Denote the feature dimension of layer \(\ell\) by \(d_\ell\). Hence we have \(h^{(\ell)}_u \in \mathbb{R}^{\ell}\) and \(\phi(h^{(\ell)}_v, h^{(\ell)}_u) \in \mathbb{R}^{d_{\ell +1}}\).

As shorthand notation, we collect the latent representations in layer \((\ell)\) into a feature matrix \(X^\ell\) whose columns are given as \[\begin{align} [X^\ell]_{u:} = h^{(\ell)}_u. \end{align}\]

We then define the map \(\Phi: \mathbb{R}^{N\times d_\ell} \rightarrow \mathbb{R}^{N\times N \times d_{\ell +1}}\) componentwise as \[\begin{align} [\Phi(X^\ell)]_{vu:} =\phi([X^{\ell}]_{v:},[X^{\ell}]_{u:}) \equiv \phi(h^{\ell}_{v},h^{\ell}_{u}). \end{align}\] With this we find \[\begin{align} [X^{\ell + 1}]_{v:} = \sum_u [\psi(L)]_{vu} \;[\Phi(X^\ell)]_{vu:}. \end{align}\] This now allows us to prove:

Combining the above Theorem with (31 ) then yields

Here the last equality follows since we are taking the maximum over all used propagation matries, and only ever use a single one.

Hence we have established a bound on variations of node-wise latent embeddings in terms of the heat kernel distance \(\left\| e^{-tL_1} - e^{-tL_2}\right\|\) for spectral Laplace transform based methods.

14.2 Graph Level Stability and Proof of Corrolary 2↩︎

Next we promote the node-level stability results we derived in Subsection 14.1 above to graph level stability results:

We may summarize Theorems 18 & 2 as: Given an input node-feature matrix \(X\) on graphs \(G_1, G_2\) sharing a common node set, there exists cosntants \(C\) so that the output variation \(\|X_1^K - X_2^K\|\) of \(K\)-layer deep Laplace transform based GNNs may be bounded as \[\begin{align} \label{eq:fund95ineq} \|X_1^K - X_2^K\| \leq C \max_i \left\{\int_0^\infty \hat{\psi}_i(t) \|e^{-tL_1} - e^{t L_2}\| dt \right\}. \end{align}\tag{34}\]

We are now interested in how these stability results translate to graph level features. To this end we first specify our graph level aggregation method \(\Omega\):

Graph level features \(F_1, F_2\) are then generates as \(F_1 = \Omega(X^K_1)\) and \(F_2 = \Omega(X^K_2)\).

With this we find:

This then proves Corrolary 2.

15 Stability of latent representations generated by Laplace-transform based GNNs for graphs on differing resolution scales↩︎

In this section, we establish bounds on differences in latent representations generated for Laplace-transform based GNNs when confronted with graphs describing the same underlying object at multiple resolutions. We briefly recall the setting:

We have a high resolution graph \(G\) with associated Laplacian \(L\) and node-feature matrix \(X\). We also have a lower resolution graph \(\underline{G}\), with associated Laplacian \(\underline{L}\) and node-feature matrix \(\underline{X} := J^\downarrow X\) arising from the original node feature matrix \(X\) on \(G\) via projection to \(\underline{G}\). Here we have made use of the projection opreator \(J^\downarrow\) introduced in Section 6.1 which averages node information over clusters that are condensed into supernodes. Interpolation \(J^\uparrow\) instead assings information at each supernode in \(\underline{G}\) to every node making up the corresponding cluster in \(G\).

We have [32] (c.f. also Appendix 17) that \(J^\downarrow J^\uparrow = I_{\underline{G}}\) (i.e. equal to the identity matrix on \(\underline{G}\)).

15.1 Node Level Results↩︎

We are then initially interested in the following question: Suppose we have a node feature matrix \(X\) on the graph \(G\). We can generate node-level latent embeddings \(\Xi(X)\) by feeding this node-wise information into a (Laplace-transform based) node-level GNN \(\Xi\). How different is this outcome from projecting \(X\) onto the coarser graph \(\underline{G}\), running the GNN there to generate node-level embeddings \(\Xi(J^\downarrow X)\) and then interpolating the generated embeddings back up to \(G\) via the interpolation operator \(J^\uparrow\). That is to say, we are interested in the difference \(\|\Xi(X) - J^\uparrow \Xi(J^\downarrow X)\|\).

15.1.1 Spectral Laplace-transform methods↩︎

We begin by establishing \(\|\Xi(X) - J^\uparrow \Xi(J^\downarrow X)\| \leq C \cdot \max\limits_{i} \left\{ \int_0^\infty |\hat{\psi}_i(t)| \cdot \|e^{-tL_\omega} - J^\uparrow e^{- t\underline{L}} J^\downarrow\| \right\}\) for spectral methods:

15.1.2 Message passing based Laplace Transform Methods:↩︎

We next establish similar results for message passing networks. Sticking with the notation of the preceding Subsection 15.1.1, we are hence interested in bounding \(X^K - J^\uparrow \underline{X}^K\), with \(\underline{X}\) being the node feature matrix generated after a \(K\)-layer Laplace transform based MPNN.

Finally we provide the Lemma used in the proof above:

15.2 Graph Level Results and Proof of Theorem 1↩︎

It remains to translate the results of Theorems 20 & 4 to the graph level. To this end, we first note the following:

Using the above lemma, we find \[\begin{align} \|\underline{F} - F\| = \| \Omega(\underline{X}^K) - \Omega(X^K)\| = \| \Omega(J^\uparrow \underline{X}^K) - \Omega(X^K)\|. \end{align}\] Using the \(1\)-Lipschitzness of the aggregation function \(\Omega\) then yields \[\begin{align} \| \Omega(J^\uparrow \underline{X}^K) - \Omega(X^K)\|\leq \| J^\uparrow \underline{X}^K - X^K\|. \end{align}\] Together with Theorems 20 & 4 we have thus proved:

16 Additional Experimental Considerations↩︎

16.1 Additional details on initial Coarse Graining Experiment↩︎

16.1.0.1 Dataset:

The dataset we consider is the QM\(\mathbf{7}\) dataset, introduced in [28], [54]. This dataset contains descriptions of \(7165\) organic molecules, each with up to seven heavy atoms, with all non-hydrogen atoms being considered heavy. A molecule is represented by its Coulomb matrix \(C^{\text{Clmb}}\), whose off-diagonal elements \[C^{\text{Clmb}}_{ij} = \frac{Z_iZ_j}{|R_i-R_j|}\] correspond to the Coulomb-repulsion between atoms \(i\) and \(j\). We discard diagonal entries of Coulomb matrices; which would encode a polynomial fit of atomic energies to nuclear charge [28].

For each atom in any given molecular graph, the individual Cartesian coordinates \(R_i\) and the atomic charge \(Z_i\) are (in principle) also accessible individually. To each molecule an atomization energy - calculated via density functional theory - is associated. The objective is to predict this quantity. The performance metric is mean absolute error. Numerically, atomization energies are negative numbers in the range \(-600\) to \(-2200\). The associated unit is \([\textit{kcal/mol}]\).

16.1.0.2 Details on collapsing procedure:

Again, we make use of the QM\(7\) dataset [28] and its Coulomb matrix description \[\label{offdiagII} C^{\text{Clmb}}_{ij} = \frac{Z_iZ_j}{|R_i-R_j|}\tag{37}\] of molecules. We modify (all) molecular graphs in QM\(7\) by deflecting hydrogen atoms (H) out of their equilibrium positions towards the respective nearest heavy atom. This is possible since the QM\(7\) dataset also contains the Cartesian coordinates of individual atoms. Edge weights between heavy atoms then remain the same, while Coulomb repulsions between H-atoms and respective nearest heavy atom increasingly diverge; as is evident from (37 ).

Given an original molecular graph \(G\) with node weights \(\mu_i=Z_i\), the corresponding limit graph \(\underline{G}\) corresponds to a coarse grained description, where heavy atoms and surrounding H-atoms are aggregated into single super-nodes.

Mathematically, \(\underline{G}\) is obtained by removing all nodes corresponding to H-atoms from \(G\), while adding the corresponding charges \(Z_H = 1\) to the node-weights of the respective nearest heavy atom. Charges in (37 ) are modified similarly to generate the weight matrix \(\underline{W}\).

On original molecular graphs, atomic charges are provided via one-hot encodings. For the graph of methane – consisting of one carbon atom with charge \(Z_C =6\) and four hydrogen atoms of charges \(Z_H =1\) – the corresponding node-feature-matrix is e.g. given as \[\begin{align} X = \begin{pmatrix} 0 & 0 &\cdots &0 &1 & 0 \cdots \\ 1 & 0 & \cdots&0 &0 & 0 \cdots \\ 1 & 0 & \cdots&0 &0 & 0 \cdots \\ 1 & 0 & \cdots&0 &0 & 0 \cdots \\ 1 & 0 & \cdots&0 &0 & 0 \cdots \end{pmatrix} \end{align}\] with the non-zero entry in the first row being in the \(6^{\text{th}}\) column, in order to encode the charge \(Z_C=6\) for carbon.

The feature vector of an aggregated node represents charges of the heavy atom and its neighbouring H-atoms jointly.

Node feature matrices are translated as \(\underline{X} = J^\downarrow X\). Applying \(J^\downarrow\) to one-hot encoded atomic charges yields (normalized) bag-of-word embeddings on \(\underline{G}\): Individual entries of feature vectors encode how much of the total charge of the super-node is contributed by individual atom-types. In the example of methane, the limit graph \(\underline{G}\) consists of a single node with node-weight \[\mu= 6 + 1 + 1 + 1 + 1 = 10.\] The feature matrix \[\underline{X} = J^\downarrow X\] is a single row-vector given as \[\begin{align} \underline{X} = \left(\frac{4}{10},0,\cdots,0,\frac{6}{10},0,\cdots\right). \end{align}\]

16.1.0.3 Experimental Setup:

We randomly select \(1500\) molecules for testing and train on the remaining graphs. On QM\(7\) we run experiments for \(5\) different random random seeds and report mean and standard deviation. All experiments were performed on a single NVIDIA Quadro RTX 8000 graphics card.

16.1.0.4 Additional details on training and models:

Typical GNN models are divided into standard architectures (GCN [20], ChebNet [23], GATv2 [55]) GIN ([18] (together with its edge-weight enabled version [56] where applicable)) and multi- scale architectures (UFGNet [57], Lanczos [58]). Apart from UFGNet (already acting as a pooling layer) we also consider self-attention-pooling [59]; both acting on the final layer (SAG) and as acting on the output of each indivifual layer, with resulting layer-wise features concatenated to produce the final embedding (SAG-M). All considered convolutional layers are incorporated into a two layer deep and fully connected graph convolutional architecture. In each hidden layer, we set the width (i.e. the hidden feature dimension) to \(F_1 = F_2 = 64\). For ChebNet, we set the polynomial order to \(K=2\) as is customary. Lanczos uses \(20\) Lanczos iterations, as proposed in the original paper [58]. UFGNet uses Haar wavelets. For all baselines, the standard mean-aggregation scheme is employed after the graph-convolutional layers to generate graph level features. Finally, predictions are generated via an MLP.

For the Spectral\(_{\text{Res.}}\)-architecture, we set \(\lambda = 1\) and and build filters using the \(k=1\) and \(k=2\) atoms in \(\Psi^{\text{Res}} = \{(z+\lambda)^{-k}\}_{k\in\mathbb{N}}\). For MPNN\(_{\text{Res.}}\) we use the \(k=1\) atom as \(\psi\).

For the Spectral\(_{\text{Exp.}}\)-, we set \(t=1\) and and build filters using the \(k=1\) and \(=2\) atoms in \(\Psi^{\text{Exp}} = \{e^{-(kt_0)z}\}_{k\in\mathbb{N}}\). For MPNN\(_{\text{Res.}}\) we use the \(k=1\) atom as \(\psi\).

As aggregation, we employ the graph level feature aggregation scheme discuseed in Definition 19, with node weights set to atomic charges of individual atoms. Predictions are then generated via a final MLP with the same specifications as the one used for baselines.

16.2 Further discussions on stability results of Table 3↩︎

We can further understand the small difference of latent embeddings using Theorem 7.1. Let use use the shorthand notation \[\begin{align} \eta(t) = \|e^{-Lt} - J^\uparrow e^{-t\underline{L}} J^\downarrow\|. \end{align}\]

Exemplarily considering exponential propagation matrices (cf.Section 6) we have (with \(t_k = k\)) that \(\int_0^\infty |\hat{\psi}_k(t)| \|e^{-tL} - J^\uparrow e^{- t\underline{L}} J^\downarrow\| dt = \|e^{-kL} - J^\uparrow e^{- k\underline{L}} J^\downarrow\|\), we thus have \(\|F - \underline{F}\| \lesssim \max_{k \geq 1} |\eta(k)|\). Investigating the differences \(\eta(t) = \|e^{-tL} - J^\uparrow e^{- t\underline{L}} J^\downarrow\|\) of diffusion flows, we find that \(\eta(t)\) drops to zero fast, as exemplarily plotted in Fig. 23 for the first few molecules of QM\(7\). In particular \(|\eta(k)|_{k \geq 1} \approx 0\).

Using this as an upper bound in Theorem 7.1 shows that embeddings \(F, \underline{F}\) of graphs describing the same molecule at different resolution scales are similar. This explains the ability to generalize between scales.

16.3 Transferability on Graphs generated via Stochastic Block Models↩︎

16.3.0.1 Stochastic Block Models:

Stochastic block models [60] are generative models for random graphs that produce graphs containing strongly connected communities. In our experiments in this section, we consider a stochastic block model whose distributions is characterized by four parameters: The number of communities \(c_{\text{number}}\) determine how many (strongly connected) communities are present in the graph that is to be generated. The community size \(c_{\text{size}}\) determines the number of nodes belonging to each (strongly connected) community. The probability \(p_{\text{connect}}\) determines the probability that two nodes within the same community are connected by an edge. The probability \(p_{\text{inter}}\) determines the probabilities that two nodes in different communities are connected by an edge.

16.3.0.2 Experimental Setup:

Since stochastic block models do not generate node-features, we equip each node with a randomly-generated unit-norm feature vector. Given such a graph \(G\) drawn from a stochastic block model, we then compute a version \(\underline{G}\) of this graph, where all communities are collapsed to single nodes. We then compare the feature vectors generated for \(G\) and \(\underline{G}\). All experiments were performed on a single NVIDIA Quadro RTX 8000 graphics card.

16.3.0.3 Experiment: Varying the Connectivity within the Communities:

We fix the parameters \(c_{\text{number}}, c_{\text{size}}\) and \(p_{\text{inter}}\) in our stochastic block model. We then vary the probability \(p_{\text{connect}}\) that two nodes within the same community are connected by an edge from \(p_{\text{connect}} = 0\) to \(p_{\text{connect}} = 1\). This corresponds to varying the connectivity within the communities from very sparse (or in fact no connectivity) to full connectivity (i.e. the community being a clique). In Figure 24 below, we then plot the difference of feature vectors generated by models for \(G\) and \(\underline{G}\) respectively. For each \(p_{\text{connect}} \in [0,1]\), results are averaged over \(100\) graphs randomly drawn from the same stochastic block model.

We have chosen \(p_{\text{inter}} = 2/c_{\text{size}}^2\) so that – on average – clusters are connected by two edges. The choice of two edges (as opposed to \(1, 3,4,5,...\)) between clusters is not important.

16.4 Node Level Transferability and Graphs with Varying Connectivity↩︎

Here, we duplicated individual nodes on popular node-classification datasets (Citeseer & Cora [61], [62]) \(k\)-times to form (fully connected) \(k\)-cliques, while keeping the train-val-test partition constant.

Models were then trained on the same (\(k\)-fold expanded) train-set and asked to classify nodes on the (\(k\)-fold expanded) test-partition. Baselines were chosen to form a representative selection of common information-propagation methods.

16.4.0.1 Additional details on training and models:

All experiments were performed on a single NVIDIA Quadro RTX 8000 graphics card. We closely follow the experimental setup of [63] on which our codebase for this experiment builds: All models are trained for a fixed maximum (and unreachably high) number of \(n=10000\) epochs. Early stopping is performed when the validation performance has not improved for \(100\) epochs. Test-results for the parameter set achieving the highest validation-accuracy are then reported. Ties are broken by selecting the lowest loss (c.f. [19]). Confidence intervals are calculated over multiple splits and random seeds at the \(95\%\) confidence level via bootstrapping.

We train all models on a fixed learning rate of \(\text{lr} = 0.1\). Global dropout probability \(p\) of all models is optimized individually over \(p \in \{0.3, 0.35, 0.4, 0.45, 0.5\}\). We use \(\ell^2\) weight decay and optimize the weight decay parameter \(\lambda\) for all models over \(\lambda \in \{0.0001, 0.0005\}.\) We choose a two-layer deep convolutional architecture with the dimensions of hidden features optimized over \[\label{hidden95dims} K_\ell \in \{32,64,128\}.\tag{38}\]

In addition to the hyperparemeters specified above, some baselines have additional hyperparameters, which we detail here: ChebNet uses \(K=2\) to avoid the known over-fitting issue [20] for higher polynomial orders. The graph attention network [19] uses \(8\) attention heads, as suggested in [19].

For the LTF-models, we optimize depth over \(K=1,2\) with hidden feature dimension optimized over the values in (38 ) as for baselines. We empirically observed in the setting of unweighted graphs, that rescaling the Laplacian as \(L_{nf} := \frac{1}{c_{nf}}L\) with a normalizing factor \(c_{nf}\) on which we base LTF architectures improved performance. We express this normalizing factor in terms of the largest singular value \(\|\Delta\|\) of the (non-normalized) graph Laplacian. It is then selected among \(c_{nf}/\|\Delta\| \in \{0.001,0.01,0.1,2\}\). The value \(\lambda\) for the resolvent is selected among \(\lambda \in \{0.14,0.15, 0.2, 0.25\}\).

16.5 The Setting og Graphs discretizing a common ambient Manifold↩︎

The concept of operators capturing the geometry of underlying spaces also applies to manifolds \(\mathcal{M}\), where the Laplace-Beltrami operator \(\Delta_{\mathcal{M}}\) can be thought of as a continuous analogue of the Graph Laplacian [64]. This is hence is a prime setting for studying transferability. Counter to previous works [11], [65], our framework here allows to derive transferability guarantees beyond the settings of bandlimited signals and probabalistic guarantees.

We consider the setting of two graphs \(G_1, G_2\) discretely approximating the same manifold (c.f. e.g. Fig. 26). This can be made mathematically precise using the concept of generalized norm resolvent convergence (c.f. e.g. [66] for a discussion). Here we note the following: Given projection operators \(J^{\downarrow}_i\) mapping from \(\mathcal{M}\) to \(G_i\) and interpolation operators \(J^{\uparrow}_i\) mapping from \(G_i\) to \(\mathcal{M}\), we may measure the difference \(\|e^{-t\Delta_{\mathcal{M}}} - J_i^{\uparrow}e^{-tL_i}J^{\downarrow}_i\| \leq \delta_i\) in heat diffusion flows on the respective spaces. The fidelity of the discrete approximation is then essentially determined by the size of \(\delta_i \ll 1\). As discussed in detail below, we have in this setting: \[\label{torus95triangle} \|e^{-tL_1}- (J_1^\downarrow J_2^{\uparrow})e^{-tL_2} (J^{\downarrow}_2 J_1^\uparrow)\| \lesssim (\delta_1 + \delta_2)\tag{39}\]

Together with Theorem 1, this then also implies \[\begin{align} \label{eq:cauchy} \|F_1 - F_2\| \lesssim (\delta_1 + \delta_2) \end{align}\tag{40}\] for the generated latent embeddings of the two graph approximations \(G_1\) and \(G_2\). Below we show using the explicit example of the torus, that as the approximating mesh becomes finer and finer (i.e. the number of nodes \(N\) increases), the corresponding approximation error \(\delta\) tends to zero. Equation (40 ) then guarantees that the corresponding graph level latent embeddings form a Cauchy sequence. Since the latent space \(\mathbb{R}^d\) is complete, we thus have convergence of the latent embeddings \(F_N\) (indexed by the mesh size \(N\)) towards a limit embedding (\(F_N \rightarrow F_\infty\)), with the limit embedding \(F_\infty\) corresponding to the latent embedding our Laplace-transform based model would generate if we could to deploy it on the continuous manifold \(\mathcal{M}\).

Hence let us further discuss the setting of two graphs \(G_i\) discretizing the same ambient manifold \(\mathcal{M}\) in the sense of \[\begin{align} \|J^{\uparrow }_ie^{-tL_i}J^{\downarrow}_i- e^{-t\Delta_\mathcal{M}}\| \leq \delta. \end{align}\]

We will assume \(J^{\downarrow}_iJ_i^{\uparrow} = Id_{G_i}\), which is a justified assumption, as out torus example below elucidates. In this setting, we then have

\[\begin{align} &\|e^{-tL_1} - (J_1^\downarrow J_2^{\uparrow})e^{-tL_2}(J^{\downarrow}_2 J_1^\uparrow)\| \\ =& \|e^{-tL_1} - J_1^\downarrow e^{-t\Delta_\mathcal{M}}J_1^{\uparrow} + J_1^\downarrow(\Delta_\mathcal{M} + Id)^{-1}J_1^{\uparrow} - (J_1^\downarrow J_2^{\uparrow})e^{-tL_2}(J^{\downarrow}_2 J_1^\uparrow)\|\\ \leq& \|e^{-tL_1} - J_1^\downarrow e^{-t\Delta_\mathcal{M}}J_1^{\uparrow} \| + \| J_1^\downarrow e^{-t\Delta_\mathcal{M}}J_1^{\uparrow} - (J_1^\downarrow J_2^{\uparrow})e^{-tL_2}(J^{\downarrow}_2 J_1^\uparrow)\|\\ \end{align}\]

We note \[\begin{align} &\|e^{-tL_1} - J_1^\downarrow e^{-t\Delta_\mathcal{M}}J_1^{\uparrow} \| \\ =& \|J^{\downarrow}_1J_1^{\uparrow} e^{-tL_1 } J^{\downarrow}_1J_1^{\uparrow} - J_1^\downarrow e^{-t\Delta_\mathcal{M} }J_1^{\uparrow} \|\\ \leq&\|J_1^\downarrow \|\|J_1^{\uparrow} \| \cdot \|e^{-tL_1} - J_1^{\uparrow}e^{-t\Delta_\mathcal{M}}J^{\downarrow}_1\| \lesssim \delta. \end{align}\]

We consider:

\[\begin{align} &\| e^{-t\Delta_\mathcal{M}} - (J_1^\downarrow J_2^{\uparrow})e^{-tL_2}(J^{\downarrow}_2 J_1^\uparrow)\|\\ \leq&\|J_1^\downarrow \|\|J_1^{\uparrow} \| \cdot \| e^{-t\Delta_\mathcal{M}} -J_2^{\uparrow}e^{-tL_2}J^{\downarrow}_2 \|\\ \lesssim& \| e^{-t\Delta_\mathcal{M} } - J_2^{\uparrow}e^{-tL_2}J^{\downarrow}_2 \| \leq \delta.\\ \end{align}\]

Hence we have indeed established

\[\|e^{-tL_1} - (J_1^\downarrow J_2^{\uparrow})e^{-tL_2}(J^{\downarrow}_2 J_1^\uparrow)\| \lesssim 2 \delta.\]

Next let us consider an explicit example.

To this end, let us revisit the torus-setting introduced in Fig. 12. To see the origin of the estimate (39 ) in the torus setting, we begin by recalling that the standard torus \(\mathbb{T}\) arises as the cartesian product of two circles \(S_1\) of circumference \(2\pi\): \[\mathbb{T} = S^1 \times S^1.\] Let us parametrize these circles via angles \(0 \leq \theta_1, \theta_1 \leq 2 \pi\). The Laplacian on \(\mathbb{T}\) can then be written as \[\Delta_{\mathbb{T}} = -\partial_{\theta_1}^2 -\partial_{\theta_2}^2.\] A set of corresponding normalized eigenfunctions are given as \[\phi_{k_1,k_2} = \frac{1}{2 \pi} e^{-ik_1\theta_1} e^{-ik_2\theta_2}\] with corresponding eigenvalues \[\lambda_{k_1,k_2} = k_1^2 + k_2^2\] and \(k_1, k_2 \in \mathbb{Z}\).

We now consider a regular discretization of \(\mathbb{T}\) using \(N^2\) nodes. This mesh can be thought of as arising from regular discretizations of each \(S^1\) factor; with a node being placed at angles \(\phi = \frac{2\pi}{N}k\) with \(0 \leq k \leq N\). The individual node weight of each node in the mesh discretization of \(\mathbb{T}\) is set to \(\mu = \frac{(2\pi)^2}{N^2}\). We might think of this discretization \(\mathbb{T}_N\) pf \(\mathbb{T}\) as arising via a cartesian product of the group \(\mathbb{Z}/N\mathbb{Z}\) (i.e. the group of integers modulo \(N\)) with itself. Each node of \(\mathbb{T}_N = \mathbb{Z}/N\mathbb{Z} \times \mathbb{Z}/N\mathbb{Z}\) is then specified by a tuple \((a,b) \in \mathbb{T}_N\), with \(a\in \mathbb{Z}/N\mathbb{Z}\) and \(b\in \mathbb{Z}/N\mathbb{Z}\).

The graph Laplacian \(L_N\) on \(\mathbb{T}_N = \mathbb{Z}/N\mathbb{Z} \times \mathbb{Z}/N\mathbb{Z}\) then acts on a scalar node signal \(x_{ab}\) as \[(L_N x)_{ab} = \frac{N^2}{(2\pi)^2}\left(4x_{ab} - x_{(a+1)b}- x_{(a-1)b} - x_{a(b+1)} - x_{a(b-1)} \right).\] Henceforth we will adopt the notation \(x(a,b) \equiv x_{ab}\).
Normalized eigenvectors for this Laplacian \(L_N\) on \(\mathbb{T}_N\) are given as \[\phi_{k_1,k_2}^N = \frac{1}{2\pi}e^{-i\frac{2\pi k_1}{N}a}e^{-i\frac{2\pi k_1}{N}b}\] with \(0 \leq k_1,k_2 \leq (N-1)\). Corresponding eigenvalues are found to be \[\lambda^N_{k_1,k_2} = \frac{N^2}{\pi^2}\left[\sin^2\left(\frac{\pi}{N}\cdot k_1\right) + \sin^2\left(\frac{\pi}{N}\cdot k_2\right)\right].\]

To facilitate contact between \(\mathbb{T}\) and its graph approximation \(\mathbb{T}_N\), we define an interpolation operator \(J^\uparrow_N\) that maps a graph signal \(f(a,b)\) defined on \(\mathbb{T} = \mathbb{Z}/N\mathbb{Z}\times \mathbb{Z}/N\mathbb{Z}\) to a function \(\overline{f}\) defined on \(\mathbb{T}\) by defining \[\overline{f}(\theta_1,\theta_2) = f(a,b)\] whenever \(\frac{2\pi}{N}(a-1) \leq\theta_1 \leq \frac{2\pi}{N}a\) and \(\frac{2\pi}{N}(b-1) \leq \theta_2 \leq \frac{2\pi}{N}b\).
We then take \(J^\downarrow\) to be the adjoint of \(J^\uparrow\) (i.e. \(J^\downarrow = (J^\uparrow)^*\). It is not hard to see that \(J^\downarrow J^\uparrow = Id_{\mathbb{T}_N}\).
We now want to show that (for \(t > 0\)) \[\begin{align} \label{torus95convergence} \|e^{-t\Delta_\mathbb{T}} - J^\uparrow e^{-tL_N}J^\downarrow\| \rightarrow 0 \end{align}\tag{41}\] as \(N \rightarrow \infty\). To this end, denote by \(P_{k_1,K_2}\) the orthogonal projection onto \(\phi_{k_1,k_2}\). Denote by \(P^N_{k_1,K_2}\) the orthogonal projection onto \(\overline{\phi_{k_1,k_2}^N}\). We note \[\begin{align} \|e^{-t\Delta_\mathbb{T}} - J^\uparrow e^{-tL_N}J^\downarrow\| = \left\| \sum_{k_1,k_2 \in \mathbb{Z}} e^{-\lambda_{k_1,k_2}t} P_{k_1,k_2} - \sum\limits_{- -\frac{N-1}{2} \leq p_1, p_2 \leq \frac{N-1}{2}} e^{-\lambda_{k_1,k_2}t} P^N_{p_1,p_2} \right\|. \end{align}\] From this we observe \[\begin{align} &\|e^{-t\Delta_\mathbb{T}} - J^\uparrow e^{-tL_N}J^\downarrow\| = \left\| \sum_{k_1,k_2 \in \mathbb{Z}} e^{-\lambda_{k_1,k_2}t} P_{k_1,k_2} - \sum\limits_{- -\frac{N-1}{2} \leq p_1, p_2 \leq \frac{N-1}{2}} e^{-\lambda^N_{p_1,p_2}t} P^N_{p_1,p_2} \right\|\\ \leq & \left\|\sum_{\frac{N-1}{2} < |k_1|,|k_2| } e^{-\lambda_{k_1,k_2}t} P_{k_1,k_2} \right\| + \left\| \sum\limits_{- -\frac{N-1}{2} \leq k_1, k_2 \leq \frac{N-1}{2}} \left( e^{-\lambda_{k_1,k_2}t} P_{k_1,k_2} - e^{-\lambda^N_{k_1,k_2}t} P^N_{k_1,k_2} \right) \right\| \end{align}\] For the first summand, we already have \[\begin{align} \left\|\sum_{\frac{N-1}{2} < |k_1|,|k_2| } e^{-\lambda_{k_1,k_2}t} P_{k_1,k_2} \right\| \leq e^{-t\frac{(N-1)^2}{2}}. \end{align}\] Hence let us investigate the second summand. We note \[\begin{align} \label{hardstyle} &\left\| \sum\limits_{-\frac{N-1}{2} \leq k_1, k_2 \leq \frac{N-1}{2}} \left( e^{-\lambda_{k_1,k_2}t} P_{k_1,k_2} - e^{-\lambda^N_{k_1,k_2}t} P^N_{k_1,k_2} \right) \right\| \\ \leq &\left\| \sum\limits_{- \frac{N-1}{2} \leq k_1, k_2 \leq \frac{N-1}{2}} \left( e^{-\lambda_{k_1,k_2}t}- e^{-\lambda^N_{k_1,k_2}t} \right) P^N_{k_1,k_2} \right\| + \left\| \sum\limits_{ -\frac{N-1}{2} \leq k_1, k_2 \leq \frac{N-1}{2}} e^{-\lambda_{k_1,k_2}t} ( P_{k_1,k_2} - P^N_{k_1,k_2} ) \right\| \end{align}\tag{42}\] For the first summand we note \[\begin{align} &\left\| \sum\limits_{- \frac{N-1}{2} \leq k_1, k_2 \leq \frac{N-1}{2}} \left( e^{-\lambda_{k_1,k_2}t}- e^{-\lambda^N_{k_1,k_2}t} \right) P^N_{k_1,k_2} \right\| \\ = &\sup\limits_{- \frac{N-1}{2} \leq k_1, k_2 \leq \frac{N-1}{2}} \left|e^{-\lambda_{k_1,k_2}t}- e^{-\lambda^N_{k_1,k_2}t} \right|\\ = &\sup\limits_{- \frac{N-1}{2} \leq k_1, k_2 \leq \frac{N-1}{2}} e^{-t(k_1^2 + k_2^2)} \left|1 - e^{-t\left( \frac{N^2}{\pi^2}\sin^2\left(\frac{\pi}{N}k_1\right) - k_1^2\right)} e^{-t\left( \frac{N^2}{\pi^2}\sin^2\left(\frac{\pi}{N}k_2\right) - k_2^2\right)}\right| \end{align}\] We note \[\begin{align} \left( \frac{N^2}{\pi^2}\sin^2\left(\frac{\pi}{N}k\right)-k^2\right) = \mathcal{O}\left(\frac{k^4}{N^2}\right). \end{align}\] Using \[\frac{N^2}{\pi^2}\sin^2\left(\frac{\pi}{N}N^{\frac{1}{3}}\right) \lesssim N^\frac{2}{3}\] we note \[\begin{align} &\sup\limits_{- \frac{N-1}{2} \leq k_1, k_2 \leq \frac{N-1}{2}} e^{-t(k_1^2 + k_2^2)} \left|1 - e^{-t\left( \frac{N^2}{\pi^2}\sin^2\left(\frac{\pi}{N}k_1\right) - k_1^2\right)} e^{-t\left( \frac{N^2}{\pi^2}\sin^2\left(\frac{\pi}{N}k_2\right) - k_2^2\right)}\right|\\ \leq & \sup\limits_{|k_1|, |k_2| \leq N^{\frac{1}{3}}} e^{-t(k_1^2 + k_2^2)} \left|1 - e^{-t\left( \frac{N^2}{\pi^2}\sin^2\left(\frac{\pi}{N}k_1\right) - k_1^2\right)} e^{-t\left( \frac{N^2}{\pi^2}\sin^2\left(\frac{\pi}{N}k_2\right) - k_2^2\right)}\right|\\ + & \sup\limits_{|k_1|, |k_2| > N^{\frac{1}{3}}} e^{-t(k_1^2 + k_2^2)} \left|1 - e^{-t\left( \frac{N^2}{\pi^2}\sin^2\left(\frac{\pi}{N}k_1\right) - k_1^2\right)} e^{-t\left( \frac{N^2}{\pi^2}\sin^2\left(\frac{\pi}{N}k_2\right) - k_2^2\right)}\right|\\ &\leq e^{-t(2N^{\frac{2}{3}})} + e^{-t(2N^{\frac{2}{3}})} + e^{-t(N^{\frac{2}{3}})} . \end{align}\] Hence it remains to bound the second summand in (42 ). We note

\[\begin{align} &\left\| \sum\limits_{ -\frac{N-1}{2} \leq k_1, k_2 \leq \frac{N-1}{2}} e^{-\lambda_{k_1,k_2}t} ( P_{k_1,k_2} - P^N_{k_1,k_2} ) \right\|\\ \leq & \sum\limits_{ |k_1|, |k_2| \leq \frac{N-1}{2}} e^{-(k_1^2 + k_2^2)t} \|P_{k_1,k_2} - P^N_{k_1,k_2} \|. \end{align}\]

Next we note \[\|P_{k_1,k_2} - P^N_{k_1,k_2} \| \leq 2 \left\|\phi_{k_1,k_2} - \phi_{k_1,k_2}\right\|.\]

It is not hard to see that \[\begin{align} \left\|\phi_{k_1,k_2} - \overline{\phi^N_{k_1,k_2}}\right\| \leq 2C(|k_1| + |k|_2) \frac{2\pi}{N} \end{align}\] for some appropriately chosen \(C > 0\). Hence we have

\[\begin{align} &\left\| \sum\limits_{ -\frac{N-1}{2} \leq k_1, k_2 \leq \frac{N-1}{2}} e^{-\lambda_{k_1,k_2}t} ( P_{k_1,k_2} - P^N_{k_1,k_2} ) \right\|\\ \leq& \sum\limits_{ |k_1|, |k_2| \leq \frac{N-1}{2}} e^{-(k_1^2 + k_2^2)t} \cdot 2C(|k_1| + |k|_2) \frac{2\pi}{N}\\ =& \mathcal{O}(1/N). \end{align}\]

Where the lass claim follows from summability in \(k_1,k_2\). Thus we have in total indeed established that (41 ) holds.

In our experiments we then make use of the operators \(J^{\uparrow\downarrow}_i\) defined above. The function \(f\in L^2(\mathcal{M})\) on the torus from which node features are derived is chosen as \(f =\frac{1}{4 \pi^2} \sin(\phi)\cos(\theta)\). All networks are asked to predict a scalar signal on the respective graphs.

17 Further Discussion of the setting of coarse-graining Graphs↩︎

Let us consider graphs that contain clusters of nodes which are connected by significantly larger edge weights than those of edges outside of these clusters. From a diffusion perspective, information in a graph equalizes faster along edges with large weights.

In the limit where edge-weights within certain sub-graphs tend to infinity, information within these clusters equalizes immediately. Such clusters thus effectively behave as single nodes. We might thus consider a coarse grained graph \(\underline{G}\) where strongly connected clusters are fused together and represented only via single nodes. This naturally leads to the notion of graph coarsification, as first formalized and studied in [30], [67].

In our case at hand the node set \(\underline{\mathcal{G}}\) of the coarse grained graph \(\underline{G}\) is then given by the set of connected components in \(G_{\text{cluster}}\) (c.f. Fig 29). Edges \(\underline{\mathcal{E}}\) are given by elements \((R,P) \in \underline{\mathcal{G}} \times \underline{\mathcal{G}}\) with non-zero accumulated edge weight \(\underline{W}_{RP} = \sum_{r\in R}\sum_{p\in P} W_{rp}\). Node weights in \(\underline{G}\) are defined accordingly by

interpolation is defined as \(J^\uparrow u = \sum_{R \in \underline{\mathcal{G}}} u_R \cdot\mathbb{1}_R\). In this setting, we have (c.f. the discussion below) that

\[\label{orig95exp95conv95result} \|e^{-tL} - J^\uparrow e^{- t\underline{L}} J^\downarrow\| \lesssim 1/w_{\text{high}}^{\text{min}} \;\;\text{for any} \;\;t>0.\tag{43}\]

Here \(w_{\text{high}}^{\text{min}} \gg 1\) denotes the minimal edge weight inside the strongly connected clusters in \(G\). As the strength of the edge-weights in \(G_{\text{cluster}}\) tends to infinity, we have by (43 ) that also \(\eta(t) = \|e^{-Lt} - J^\uparrow e^{- \underline{L}t} J^\downarrow\| \rightarrow 0\) for any \(t >0\). Thus (for \(t>0\)) the diffusion process \(e^{-Lt}\) on \(G\) acts essentially as first projecting the input-signal to \(\underline{G}\) via \(J^\downarrow\), then diffusing information on the coarse grained graph \(\underline{G}\) via \(e^{-\underline{L}t}\) and finally interpolating back to the original graph \(G\) via \(J^\uparrow\).

Following the proof in [32], we here illustrate \[\begin{align} \|(L + Id)^{-1} - J^\uparrow (\underline{L} + Id)^{-1} J^\downarrow\| \lesssim 1/\lambda_{1}(L_{\text{high}}). \end{align}\] in this Appendix. Using [32], this then also implies \[\label{orig95exp95conv95result95app} \|e^{-tL} - J^\uparrow e^{- t\underline{L}} J^\downarrow\| \lesssim 1/w_{\text{high}}^{\text{min}} \;\;\text{for any} \;\;t>0,\tag{44}\] after noting the linear relation in scaling behaviour \(\lambda_1(L_{\text{cluster}}) \sim w_{\text{high}}^{\text{min}}\).

For convenience, we here explicitly restate the definitions leading up to this setting in full clarity and generality:

In order to translate signals between the original graph \(G\) and the limit description \(\underline{G}\), we need translation operators mapping signals from one graph to the other:

Slightly extending established results in [32], we then have the following result:

\[\lambda_{\max}(L_{\text{reg.}}) = \|L_{\text{reg.}}\|.\] We here inclue the proof of this fact for convenience and self-containedness.

As a corollary, one finds

18 Computational cost analysis↩︎

To evaluate the practical efficiency of our proposed framework, we conducted a thorough empirical investigation of the computational costs associated with constructing and operating on a representative propagation matrix—specifically, the resolvent formulation introduced in Def, 1.

18.1 Computational Efficiency on Small-Scale Graphs (QM7)↩︎

We first established a baseline for resource consumption during model training on the QM7 molecular dataset. To quantify computational overhead, we monitored key performance metrics throughout training:

• Peak GPU Memory Allocation was tracked using torch.cuda.memory.max_memory_allocated.

• Training Time per Epoch was measured as the wall-clock time required to complete a single forward and backward pass across all training samples.

As evident from Table 5, our method has significantly less overhead than standard GNN baselines.

18.2 Scalability to Large-Scale Graphs↩︎

To assess how our approach scales to larger, more complex topologies, we extended our evaluation to datasets with significantly higher node counts, namely ms_academic (\(>18\),\(000\) nodes) and pubmed (\(>19\),\(000\) nodes).

18.2.0.1 Verification of effective sparsity:

We first experimentally verify the claim made towards the end of Section 7.1 that the propagation matrix is effectively sparse.

As is evident from Figure 30, the overwhelming majority of entries \(R_{ij}\) in the resolvent propagation matrix \(R\) are indeed effectively zero. Indeed, we list results for small cutoff values \(\epsilon \geq |R_{ij}|\) in Table 6 below.

18.2.0.2 Computational overhead analysis

Using the identical logging infrastructure and metrics as in the QM7 experiments, we tracked the scaling behavior of memory allocation and epoch duration.

Our empirical results validate the theoretical assertions presented in our scalability discussion at the end of Section 7.1.

References↩︎

---

1. Appendix 11 further summarizes these existing perturbation models, also discussing their inapplicability to our setting.↩︎

2. At an even higher resolution scale, also interactions of quarks and gluons (in turn making up protons and neutrons) are neglegted.↩︎

3. Appendix 16.1 contains further details on models.↩︎

4. Section 6.1 discusses the precise sense in which this convergence is formalized mathematically.↩︎

5. Appendix 13 contains derivations for all standard GNN types.↩︎

6. We are using the standard notation ‘\(\lesssim\)’ to denote smaller than, up to a fixed ‘universal’ constant.↩︎