Spatio-Spectral Graph Neural Networks

Authors: simon  geisler
arthur  kosmala
daniel  herbst
stephan  günnemann
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 22:32:08

Abstract: Spatial Message Passing Graph Neural Networks (MPGNNs) are widely used for learning on graph-structured data. However, key limitations of *ℓ*-step MPGNNs are that their "receptive field" is typically limited to the *ℓ*-hop neighborhood of a node and that information exchange between distant nodes is limited by over-squashing. Motivated by these limitations, we propose *Spatio-Spectral Graph Neural Networks (S²GNNs)* – a new modeling paradigm for Graph Neural Networks (GNNs) that synergistically combines spatially and spectrally parametrized graph filters. Parameterizing filters partially in the frequency domain enables global yet efficient information propagation. We show that S²GNNs vanquish over-squashing and yield strictly tighter approximation-theoretic error bounds than MPGNNs. Further, rethinking graph convolutions at a fundamental level unlocks new design spaces. For example, S²GNNs allow for free positional encodings that make them strictly more expressive than the 1-Weisfeiler-Leman (WL) test. Moreover, to obtain general-purpose S²GNNs, we propose spectrally parametrized filters for directed graphs. S²GNNs outperform spatial MPGNNs, graph transformers, and graph rewirings, e.g., on the peptide long-range benchmark tasks, and are competitive with state-of-the-art sequence modeling. On a 40 GB GPU, S²GNNs scale to millions of nodes.

Main paper: 
 introduction: 
 Spatial Message-Passing Graph Neural Networks (MPGNNs) ushered in various recent breakthroughs. For example, MPGNNs are able to predict the weather with unprecedented precision (Lam et al., 2023) , can be composed as a foundation model for a rich set of tasks on knowledge graphs (Galkin et al., 2023) , and are a key component in the discovery of millions of AI-generated crystal structures (Merchant et al., 2023) . Despite this success, MPGNNs produce nodelevel signals solely considering limited-size neighborhoods, effectively bounding their expressivity. Even with a large number of message-passing steps, MPGNNs are limited in their capability of propagating information to distant nodes due to over-squashing. As evident by the success of global models like transformers (Vaswani et al., 2017) , modeling long-range interactions can be pivotal and an important step towards foundation models that understand graphs.We propose Spatio-Spectral Graph Neural Networks (S 2 GNNs), a new modeling paradigm for tackling the aforementioned limitations, that synergistically combine message passing with spectral filters, explicitly parametrized in the spectral domain. Spectral filters are virtually ignored by prior work but go beyond stacks of message-passing layers or polynomial parametrizations. Due to message passing's finite number of propagation steps, it comes with a distance cutoff p cut (# hops, see Fig. 1 ). Conversely, spectral filters act globally (p max ), even on a truncated frequency spectrum λ cut . Truncating the frequency spectrum for spectral filters is required for efficiency, yet message passing has access to the entire spectrum (right plots in Fig. 1 ). The combination of message passing and spectral filters provably leverages the strengths of each parametrization. Utilizing this combination, S 2 GNNs generalize the concept of "virtual nodes" and distill many important properties of hierarchical message-passing schemes, graph-rewirings, and pooling into a single GNN (see Fig. 3 ). Outside of GNNs, a similar composition is at the core of some State Space Models (SSM) models (Poli et al., 2023) , that deliver transformer-like properties with superior scalability on sequences -as do S 2 GNNs on graphs.Our analysis of S 2 GNNs ( § 3.1) validates their capability for modeling long-range interactions. We prove in § 3.1.1 that combining spectral and spatial filters alleviates the over-squashing phenomenon (Alon & Yahav, 2020; Di Giovanni et al., 2023a,b) , a necessity for effective informationexchange among distant nodes. Our approximation-theoretic analysis goes one step further and proves strictly tighter error bounds in terms of approximation of the target idealized GNN ( § 3.1.2).Vanquishes over-squashing ( § 3.1.1) Superior approximation capabilities ( § 3.1.2) Design Space § 3.2 Spatial message passing (see literature) Spectral filter parametrization ( § 3.2.1) Spectral-domain MLP ( § 3.2.2) Spectral filter for directed graphs ( § 3.2.3) Dual use of partial Eigendecomposition (EVD):free-of-cost Positional Encodings (PE) ( § 3.2.4)Figure 2 : S 2 GNN framework with adjacency matrix A, node features X, and Laplacian L (function of A).Design space of S 2 GNNs ( § 3.2).Except for initial works like (Bruna et al., 2014) and in contrast to spatial MPGNNs, the design decisions for spectral filters are virtually unexplored -and so is their composition. The novel aspects of S 2 GNN's design space include the spectral filter parametrization ( § 3.2.1).We propose the first permutationequivariance-preserving neural network in the spectral domain ( § 3.2.2) and generalize spectral filters to directed graphs ( § 3.2.3). The dual use of the partial eigendecomposition, required for spectral filters, allows us to propose "free-of-cost" positional encodings ( § 3.2.4), that are permutationequivariant, stable, and increase expressivity strictly beyond the 1-Weisfeiler-Leman (WL) test.S 2 GNNs are effective and practical. We empirically verify the shortcomings of MPGNNs and how S 2 GNNs overcome them ( § 4). E.g., we set a new state-of-the-art on peptides-func (Dwivedi et al., 2022) with ≈ 35% fewer parameters, outperforming MPGNNs and graph transformers. Although sequences are just a subdomain of (directed) graphs, we also study how S 2 GNNs compare to specialized sequence models like transformers (Vaswani et al., 2017) or Hyena (Poli et al., 2023) . We find that S 2 GNNs are highly competitive even though they operate on a much more general domain (un-/directed graphs). Last, the runtime and space complexity of S 2 GNNs is equivalent to MPGNNs and, with vanilla full-graph training, S 2 GNNs can handle millions of nodes with a 40 GB GPU. 
background: 
 We study graphs G(A, X) with adjacency matrix A ∈ {0, 1} n×n (or A ∈ R n×n ≥0 if weighted), node features X ∈ R n×d and edge count m. A is symmetric for undirected graphs and, thus, has eigendecomposition λ, V = EVD(A) with eigenvalues λ ∈ R n and eigenvectors V ∈ R n×n : A = V ΛV ⊤ using Λ = diag(λ). Instead of A, we decompose the Laplacian L := I -D -1 /2 AD -1 /2 , with diagonal degree matrix D = diag(A ⃗ 1), since its ordered eigenvalues 0 = λ 1 ≤ λ 2 ≤ • • • ≤ λ n ≤ 2 are similar to frequencies (e.g., low eigenvalues relate to low frequencies, see Fig. 4 ). Likewise, one could use, e.g., L = I -D -1 A or more general variants (Yang et al., 2023) ; however, we focus our explanations on the most common choice L := I -D -1 /2 AD -1 /2 . We choose the matrix of eigenvectors V ∈ R n×n to be orthogonal V V ⊤ = I. We refer to V as the Fourier basis of the graph, with Graph Fourier Transformation (GFT) X = V ⊤ X and its inverse X = V X. To provide an overview, Table 5 lists the symbols used in this work.Spectral graph filters. Many GNNs implement a graph convolution, where node signal X ∈ R n×d is convolved g * G X for every d with a scalar filter g ∈ R n . The graph convolution (Hammond et al., 2011) is defined in the spectral domain as g * G X := V ([V ⊤ g] ⊙ [V ⊤ X]), with elementwise product ⊙ and broadcast of V ⊤ g to match shapes. Instead of spatial g, spectral graph filters parametrize ĝ : [0, 2] → R explicitly and yield V ⊤ g := ĝ(λ) ∈ R n as a function of the eigenvalues.Figure 3 : Message-passing interpretation of V (ĝ ϑ (λ) ⊙ [V ⊤ X]) (spectral filter): via the Fourier coefficients they may exchange information globally and allow intraand inter-cluster message passing. Edge width/color denotes the magnitude/sign of V .Message Passing Graph Neural Networks (MPGNNs) circumvent the EVD via polynomial ĝ(λ) u = p j=0 γ j λ j u since V [ p j=0 γ j diag(λ) j ]V ⊤ X = p j=0 γ j L j X. In practice, many MPGNNs use p = 1: H (l) = (γ 0 I + γ 1 L)H (l-1) with H (0) = X, and stack 1 ≤ l ≤ ℓ layers interleaved with node-wise transformations and activations σ. We refer to Balcilar et al. (2021b) for similar interpretations of MPGNNs like GAT (Veličković et al., 2018) or GIN (Xu et al., 2019) .3 Method S 2 GNNs symbiotically pair spatial Spatial(H (l-1) ; A) MPGNNs and Spectral(H (l-1) ; V , λ) filters, using a partial eigendecomposition. Even though the spectral filter operates on a truncated eigendecomposition (spectrally bounded), it is spatially unbounded. Conversely, spatial MPGNNs are spatially bounded yet spectrally unbounded (see Fig. 1 ).A spectrally bounded filter is sensible for modeling global pair-wise interactions, considering its message-passing interpretation of Fig. 3 . Conceptually, a spectral filter consists of three steps: ① Gather: The multiplication of the node signal with the eigenvectors v ⊤ u X (GFT) is a weighted and signed aggregation over all nodes; ② Apply: the "Fourier coefficients" are weighted; and ③ Scatter broadcasts the signal v u X back to the nodes (inverse GFT). The first eigenvector (here for L = D -A) acts like a "virtual node" (Gilmer et al., 2017) (see also § E). That is, it calculates the average embedding and then distributes this information, potentially interlayered with neural networks. Importantly, the other eigenvectors effectively allow messages to be passed within or between clusters. As we show for exemplary graphs in Fig. 4 , low frequencies/eigenvalues capture coarse structures, while high(er) frequencies/eigenvalues capture details. For example, the second eigenvector in Fig. 4b contrasts the inner with the outer rectangle, while the third eigenspace models both symmetries up/down and left/right. In conclusion, S 2 GNNs augment spatial message-passing with a graphadaptive hierarchy (spectral filter).Thus, S 2 GNNs distill many important properties of hierarchical message-passing schemes (Bodnar et al., 2021) , graph-rewirings (Di Giovanni et al., 2023a) , pooling (Lee et al., 2019) etc. See § J.1 for more examples.S 2 GNN's composition. We study (1) an additive combination for its simpler approximationtheoretic interpretation ( § 3.1.2), or (2) an arbitrary sequence of filters due to its flexibility. In both cases, residual connections may be desirable (see § J.2).H (l) = Spectral (l) (H (l-1) ; V , λ) + Spatial (l) (H (l-1) ; A)(1)EQUATIONSpectral Filter. The building block that turns a spatial MPGNN into an S 2 GNN is the spectral filter:Spectral (l) (H (l-1) ; V , λ) = V ĝ(l) ϑ (λ) ⊙ V ⊤ f (l)θ (H (l-1) )(3) with a point-wise transformation f (l)θ : R n×d (l-1) → R n×d (l) , a learnable spectral filter ĝ(l) ϑ (λ) ∈ R k×d (l) parameterized element-wise as ĝ(l) ϑ (λ) u,v := ĝ(l) v (λ u ; ϑ v ) (see § 3.2.1), and truncated V ∈ R n×k , λ ∈ R k . Due to the combination of message passing with spectral filters, S 2 GNNs' hypothesis class goes beyond (finite-order) polynomials of the Laplacian L (or stacks message passing layers), unlocking a larger class of filters. In Algo. 1, we provide pseudo code for S 2 GNNs (Eq. 1).Truncated spectrum. We omit extra notation for the truncated eigendecompositon EVD(L, k) since it is equivalent to define ĝ(λ j ) = 0 for j > k. However, truncating after the k-th eigenvector requires care with the last eigenspance to maintain permutation equivariance. Due to the ambiguity of eigenvectors in the presence of repeated eigenvalues, we must ensure that we only include eigenspaces in their entirety. That is, we only include {λ j | j ≤ k ∧ λ j ̸ = λ k+1 }. Thus, Spectral(H (l-1) ; EVD(L, k)) is permutation equivariant nonetheless. We defer all proofs to § H. Theorem 1. Spectral(H (l-1) ; EVD(L, k)) of Eq. 3 is equivariant to all n × n permutation matrices P ∈ P: Spectral(P H (l-1) ; EVD(P LP ⊤ , k)) = P Spectral(H (l-1) ; EVD(L, k)).Complementary high-resolution filters. Our Spectral filters are highly discriminative between the frequencies and, e.g., can readily access a single eigenspace. Yet, for efficiency, we limit the spectral filter to a specific frequency band. Due to the combination with message passing, this choice of band does not decide on, say, low-pass behavior; it solely determines where to increase the spectral selectivity. While S 2 GNNs with subsequent guarantees adapt to domain-specific choices for the spectral filter's frequency band, a sensible default is to focus on the low frequencies. The two main reasons for this are (see § J.3 for an extensive list): (1) Low frequencies model the smoothest global signals w.r.t. the graph structure (see Fig. 3 & 4 ). ( 2) Under a relative perturbation model (perturbation budget proportional to degree), stability implies C-integral-Lipschitzness (∃C > 0 : |λ dĝ /dλ| ≤ C), i.e., the filter can vary strongly around zero but must level out for large λ (see Gama et al. (2020) ). Figure 5 : Spectral filters do not exhibit oversquashing on "Clique Path" graphs (Di Giovanni et al., 2023a) . 
theoretical analysis: 
 We show that how S 2 GNNs alleviate oversquashing in § 3.1.1. Next, § 3.1.2 makes the approximation-theoretic advantages precise.Alon & Yahav (2020) pointed out that MPGNNs must pass information through bottlenecks that connect different communities using fixedsize embedding vectors. Topping et al. (2022) and Di Giovanni et al. (2023a) formalize this via an L 1 -norm Jacobian sensitivity analysis:EQUATIONu . MPGNNs' Jacobian sensitivity typically decays O (exp (-r)) with node distance r if the number of walks between the two nodes is small. See § F for results of Di Giovanni et al. (2023a) . S 2 GNNs are not prone to such an exponential sensitivity decay due to their global message scheme. We formalize this in Theorem 2, refer to Fig. 4 for intuition and Fig. 5 for empirical verification. All theoretical guarantees hold if a θ exists such that f θ = I.Theorem 2. An ℓ-layer S 2 GNN can be parametrized s.t. output h (ℓ) v has a uniformly lowerbounded Jacobian sensitivity on a connected graph: ℓ) for nodes u, v ∈ G, a parameter-dependent C ϑ , network width d and edge count m.EQUATIONIn contrast to Di Giovanni et al. (2023a) , we prove a lower bound for S 2 GNNs, guaranteeing a minimum "influence" for any u on v. This is true since S 2 GNNs contain a virtual node as a special case with ĝ(l) ϑ (λ) = 1 {0} , with 1 S denoting the indicator function of a set S (see also § E). However, we find that a virtual node is insufficient for some long-range tasks, including our long-range clustering (LR-CLUSTER) of Fig. 10b . Hence, the exponential sensitivity decay of spatial MPGNNs only shows their inadequacy in long-range settings. Proving its absence is not sufficient to quantify long-range modeling capabilities, noting that the lower bound is not tight for S 2 GNNs on many graphs. We close this gap with our subsequent analysis rooted in polynomial approximation theory.To demonstrate how S 2 GNNs can express a more general hypothesis class than MPGNNs, we study how well an "idealized" GNN (IGNN) can be approximated. Each IGNN layer l can express convolution operators g (l) of any spectral form ĝ(l) : [0, 2] → R. We approximate IGNNs with S 2 GNNs from Eq. 1, with a spectral filter as in Eq. 3 and a spatial part parametrized by a polynomial. While we assume here that the S 2 GNN spectral filter is bandlimited to and a universal approximator on the interval [0, λ max ], the findings generalize to, e.g., a high-pass interval. In the main body, we focus on the key insights for architectures without nonlinear activations. Wang & Zhang (2022) prove that even linear IGNNs can produce any one-dimensional output under certain regularity assumptions on the graph and input signal. Thus, we solely need to consider a single layer. In § H.4, we cover the generic setting including nonlinearities, where multiple layers are helpful. Locality relates to spectral smoothness. The locality of the true/ideal filter g is related to the smoothness of its Fourier transform ĝ. For instance, if g is a low-order polynomial of L, it is localized to a few-hop neighborhood, and ĝ is regularized to vary slowly (Fig. 6a w/o discontinuity). The other extreme is a discontinuous spectral filter ĝ, such as the entirely non-local virtual node filter, ĝ = 1 {0} (discontinuity in Fig. 6a , details in § E). This viewpoint of spectral smoothness illuminates the limitations of finite-hop message passing from an angle that complements spatial analyses in the over-squashing picture. It informs a lower bound on the error, which shows that spatial message passing, i.e, order-p polynomial graph filters g γp with p + 1 coefficients γ p ∈ R p+1 , can converge exceedingly slowly -slower than any inverse root (!) of pto a discontinuous ground truth in the Frobenius-induced operator norm: Theorem 3. Let ĝ be a discontinuous spectral filter. For any approximating sequence g γp p∈N of polynomial filters, an adversarial sequence (G p ) p∈N of input graphs exists such that∄α ∈ R >0 : sup 0̸ =X∈R |Gp |×d ∥(g γp -g) * Gp X∥ F ∥X∥ F = O p -αSuperior S 2 GNN error bound. A spatio-spectral convolution wins over a purely spatial filter when the sharpest irregularities of the ground truth ĝ are within reach of its expressive spectral part. The spatial part, which can "focus" on learning the remaining, smoother part outside of this window, now needs much fewer hops to give a faithful approximation. We illustrate this principle in Fig. 6 where we approximate an additive combination of an order-three polynomial filter with discontinuous lowpass. Only the S 2 filter is faithfully approximating this filter. Formally, we find: Theorem 4. Assume ĝ [λcut,2] is r-times continuously differentiable on [λ cut , 2], and a boundK r (ĝ, λ cut ) ≥ 0 such that d r dλ r ĝ(λ) ≤ K r (ĝ, λ cut ) ∀λ ∈ [λ cut , 2]. An approximating S 2 GNN sequence with parameters ϑ * p , γ * p p∈N exists such that, for arbitrary graph sequences(G p ) p∈N , sup 0̸ =X∈R |Gp |×d ∥(g γ * p + g ϑ * p -g) * Gp X∥ F ∥X∥ F = O K r (ĝ, λ cut )p -rwith a scaling constant that depends only on r, not on ĝ or (G p ) p∈N .The above bound extends to purely spatial convolutions in terms of K r (ĝ, 0) if ĝ is r-times continuously differentiable on the full interval [0, 2]. The S 2 GNN bound of Theorem 4 is then still strictly tighter if K r (ĝ, λ cut ) < K r (ĝ, 0). In particular, taking the limit K 1 (ĝ, 0) → ∞ towards discontinuity makes the purely spatial upper bound arbitrarily loose, whereas a benign filter might still admit a small K 1 (ĝ, λ cut ) for some λ cut > 0. Theorem 3 suggests that this is not an artifact of a loose upper bound but that there is an inherent difficulty in approximating unsmooth filters with polynomials.We conclude the analysis by instantiating the bounds: assuming ĝ is C-integral-Lipschitz for stability reasons (see Gama et al. (2020) and the paragraph before § 3.1.1) yields K 1 (ĝ, λ cut ) = C /λcut, whereas for the electrostatics example ĝσ in § G, we find upper bounds K r (ĝ σ , λ cut ) = r! /λ (r+1). In both cases, the pure spatial bound diverges as smoothness around 0 remains unconstrained.As shown in Fig. 2 , we identify three major, yet unexplored, directions in S 2 GNNs' design space. In § 3.2.1, we discuss how we parametrize the spectral filter. In § 3.2.2, we propose the first neural network for the spectral domain. That is, we allow transformations and non-linearities in the "Fourier" domain. In § 3.2.3, we are the first to instantiate spectral filters for directed graphs. Additionally, due to the availability of the partial eigendecomposition, positional encodings may dual use them to improve epxressivity at negligible cost. In § 3.2.4, we propose the first permutation equivariant, stable and efficient positional encodings that provably admit an expressivity beyond 1-WL. § J provides further details and considerations, like some remarks on batching ( § J.7). For the (sub-) design space of spatial message passing (You et al., 2020) , we refer to its rich literature. For spectral filter function ĝϑ (λ) of Eq. 3, we learn a channelwise linear combination of translated Gaussian basis functions (see "Gaussian smearing" used by Schütt et al. ( 2017)), as depicted in Fig. 7 . This choice (1) may represent any possible ĝϑ (λ) with sufficient resolution (assumption in § 3.1.2); (2) avoids overfitting towards numerical inaccuracies of the eigenvalue calculation; (3) limits the discrimination of almost repeated eigenvalues and, in turn, should yield stability (similar to § 3.2.4). Strategies to cope with a variable λ cut and k (e.g., using attention similar to SpecFormer (Bo et al., 2023a )) did usually not yield superior experimental results.Window. We multiply the learned combinations of Gaussians by an envelope function (we choose a Tukey window) that decays smoothly to zero around cutoff λ cut . This counteracts the so-called "Gibbs phenomenon" (aka "ringing"): as visualized for a path graph/sequence of 100 nodes in Fig. 8 , trying to approximate a spatially-discontinuous target signal using an ideal low-pass range of frequency components results in an overshooting oscillatory behavior near the spatial discontinuity. Dampening the frequencies near λ cut by a smooth envelope/window function alleviates this behavior. We note that the learned filter may, in principle, overrule the windowing at the cost of a higher weight decay penalty. See Algo. 2 for ĝϑ (λ)'s algorithmic description. Depth-wise separable convolution (Sifre, 2014; Howard et al., 2017) : Applying different filters for each dimension is computationally convenient for spectral filters. While "full" convolutions are also possible, we find that such a construction is more prone to over-fitting. In practice, we even use parameter sharing and apply fewer filters than dimensions to counteract over-fitting. We argue that sharing filters among dimensions is similar to the heads in a transformer (Vaswani et al., 2017) .θ . As sketched in Fig. 3 & 4, all nodes participate in the global data transfer. While this global message-passing scheme is graph-adaptive, it does not adjust to the inputs. For adaptivity, we typically consider non-linear feature transformations f(l-1) θ (H (l-1) ), like gating mechanism f (l-1) θ (H (l-1) ) = H (l-1) ⊙σ ′ (H (l-1) W (l) G + ⃗ 1b ⊤ ))with element-wise multiplication ⊙, SiLU function σ ′ , learnable weight W , and bias b. A linear transformation f (l) θ (H (l-1) ) = H (l-1) W (l) is another interesting case since we may first apply the GFT and then the transformation: (V ⊤ H (l-1) )W (l) . Next, we extend this linear transformation to a neural network in the spectral domain by adding multiple transformations and nonlinearities.Applying a neural network s ζ in the spectral domain is highly desirable due to its negligible computational cost if k ≪ n. Moreover, s ζ allows the spectral filter to become data-dependent and may mix between channels. Data-dependent filtering is one of the properties that is hypothesized to make transformers powerful Fu et al. (2023) . We propose the first neural network for the spectral domain of graph filters s(l) ζ : R k×d (l) → R k×d (l)that is designed to preserve permutation equivariance.EQUATIONWe achieve permutation equivariance via sign equivariance sζ (S ⊙ X) = S ⊙ s ζ (X) , ∀S ∈ {-1, 1} k×d (l), combined with a permutation equivariance s ζ (P X) = P s ζ (X) , P ∈ P k , where P k is the set of all k × k permutation matrices. Specifically, we stack linear mappingsW s ∈ R d (l) ×d (l) (without bias) with a gated nonlinearity ϕ( Ĥ) = Ĥ ⊙ σ( ⃗ 1 [m ⊤ W a + b ⊤ a ]) with sigmoid σ, column-wise norm m j = ∥ Ĥ:,j ∥, and learnableW a ∈ R d (l) ×d (l)as well as b a ∈ R d (l) .Directed graphs are an important topic that did not discuss so far. For S 2 GNNs to generalize the capabilities of non-local sequence models like transformers (Vaswani et al., 2017) or SSMs (Poli et al., 2023; Gu & Dao, 2023) it is vital to support direction, e.g., for distinguishing source/beginning and sink/end. However, all discussion before assumed the existence of the eigenvalue decomposition of L. This was the case for symmetric L; however, for directed graphs, L may be asymmetric.To guarantee L is diagonalizable with real eigenvalues, we use the Magnetic Laplacian (Forman, 1993; Shubin, 1994; De Verdière, 2013) which is Hermitian and models direction in the complex domain:L q = I -(D -1 /2 s A s D -1 /2 s ) ⊙ exp[i2πq(A -A ⊤ )] with symmetrized adjacency/degrees A s /D s , potential q ∈ [0, 2π], element-wise exponential exp, and imaginary unit i 2 = -1. While other parametrizations of a Hermitian matrix are also possible, with A ∈ {0, 1} n×n and appropriate choice of q, L q : {0, 1} n×n → C n×n is injective. In other words, every possible asymmetric A maps to exactly one L q and, thus, this representation is lossless. Moreover, for sufficiently small potential q, the order of eigenvalues is well-behaved (Furutani et al., 2020) . In contrast to Koke & Cremers (2024) , a Hermitian parametrization of spectral filters does not require a dedicated propagation for forward and backward information flow. For simplicity we choose q < 1 /nmax with maximal number of nodes n max (with binary A). This choice ensures that the first eigenvector suffices to obtain, e.g., the topological sorts of a Directed Acyclic Graph (DAG). Due to the real eigenvalues of a Hermitian matrix, the presented content generalizes with minor adjustments. Most notably, we use a feature transformation f(l) θ : R n×d (l-1) → C n×d (l)and map back into the real domain after the spectral convolution. We give more implementation details in § J.6 and provide additional background on directed graphs in § C. The availability of the partial eigendecomposition allows for their dual use for positional encodings at negligible cost. Motivated by this, we propose the first efficient (O(km)) and (fully) permutation equivariant spectral Positional Encodings PE that provably increase the expressivity strictly beyond the 1-Weisfeiler-Leman (1-WL) test (Xu et al., 2019; Morris et al., 2019) . In contrast to the Laplacian encodings of Dwivedi & Bresson (2021) , our PE do not require augmenting eigenvectors w.r.t. their sign and maintain permutation equivariance also in the presence of repeated eigenvalues. In comparison to Huang et al. (2024) , our PE come with drastically lower computational cost and have no learnable parameters. Due to the absence of learnable parameters, we need to calculate our PE only once.We construct our k-dimensional positional encodings PE(V , λ) ∈ R n×k asPE(V , λ) = || k j=1 [(V ĥj (λ)V ⊤ ) ⊙ A] • ⃗ 1 (5)with concatenation || and binary adjacency A ∈ {0, 1} n×n . We use a Radial Basis Function (RBF) filter with normalization around each eigenvalue ĥj (λ) = softmax( (λj -λ)⊙(λj -λ) /σ 2 ) with small width σ ∈ R >0 . This parametrization is not only permutation equivariant but also stable according to the subsequent definition via the Hölder continuity. Note that C depends on the eigengap between 1 /(λ k+1 -λ k ) at the frequency cutoff (for exact constant C see proof in § H.5).Definition 1 (Stable PE). (Huang et al., 2024) A PE method PE : R n×k × R k → R n×k is called stable, if there exist constants c, C > 0, such that for any Laplacian L, L ′ , and P * = arg min P ∥L -EQUATIONTheorem 5. The Positional Encodings PE in Eq. 5 are stable according to Definition 1.Next to their stability, our PE can discriminate certain degree-regular graphs (e.g., Fig. 9 ). Since degree-regular graphs cannot be distinguished by 1-WL, our PE makes the equipped GNN (as expressive as 1-WL) strictly more expressive than 1-WL. See § I for continued expressivity analyses.Theorem 6. S 2 GNNs are strictly more expressive than 1-WL with the PE of Eq. 5. 
empirical results: 
 With state-of-the-art performance on the peptides-func task of the long-range benchmark (Dwivedi et al., 2022) , plus strong results on further benchmarks, we demonstrate that S 2 GCN, a GCN paired with spectral filters, is highly capable of modeling long-range interactions ( § 4.1). We assess S 2 GNNs' long sequence performance ( § 4.2) (mechanistic in-context learning) and show that S 2 GCN, a graph machine learning method, can achieve competitive results to state-of-the-art sequence models, including H3, Hyena, and transformers. We exemplify S 2 GNNs' practicality and competitiveness at scale on large-scale benchmarks ( § 4.3) like TPUGraphs (Phothilimthana et al., 2023) , PCQM4Mv2 (Hu et al., 2021) , and Open Graph Benchmark (OGB) Products (Hu et al., 2020) . Further, in § M.8, we report state-of-the-art performance on the heterophilic arXiv-year (Lim et al., 2021) and,in § M.4, we study combinations of spatial and spectral filters beyond Eq. 1 & 2.Setup. We pair different MPGNNs with spectral filters and name the composition S 2 <base>. For example, a S 2 GNN with GAT as base will be called S 2 GAT. We typically perform 3 to 10 random reruns and report the mean ± standard deviation. The experiments of § 4.1 require <11 GB (e.g. Nvidia GTX 1080Ti); for the experiments in § 4.2 & 4.3 we use a 40 GB A100. We usually optimize weights with AdamW (Loshchilov & Hutter, 2019) and cosine annealing scheduler (Loshchilov & Hutter, 2017) . We use early stopping based on the validation loss/score. See § M for more details and https://www.cs.cit.tum.de/daml/s2gnn for code as well as supplementary material.Finding (I): S 2 GCN outperforms state-of-the-art graph transformers, MPGNNs, and graph rewirings on the peptides-func long-range benchmarks (Dwivedi et al., 2022 ) by a substantial margin. Simultaneously, we remain approximately 35% below the 500k parameter threshold and. On peptides-struct we are only outperformed by NBA-GIN (Park et al., 2023) . We extend the best configuration for a GCN of Tönshoff et al. (2023) (see GCN in Table 1 ), lower the number of message passing steps from six to three, and interleave spatial and spectral filters (Eq. 2) with λ cut = 0.7. (He et al., 2023) 0.6921 ± 0.0054 0.2475 ± 0.0015 Graph ViT (He et al., 2023) 0.6942 ± 0.0075 0.2449 ± 0.0016 GRIT (Ma et al., 2023) 0.6988 ± 0.0082 0.2460 ± 0.0012 GPS+HDSE (Luo et al., 2024) 0.7156 ± 0.0058 0.2457 ± 0.0013 Rewiring: DRew-GCN (Gutteridge et al., (Giusti et al., 2023) 0.6569 ± 0.0117 0.2523 ± 0.0013 NBA-GIN (Park et al., 2023) 0.7071 ± 0.0067 0.2424 ± 0.0010 GCN (Tönshoff et al., 2023) 0 Dataset contribution: Distance regression is a task with long-range interactions used in prior work (Geisler et al., 2023; Lim et al., 2023) . Here, the regression targets are the shortest path distances to the only root node (in-degree 0). We generate random trees/DAGs with ≈750 # of nodes on average (details are in § M.7). The target distances often exceed 30 hops. We evaluate on similarly sized graphs as in the training data, i.e., in-distribution (ID) samples, and out-of-distribution (OOD) samples that consist of slightly larger graphs. Details on the dataset construction are in § M.7.Finding (II): spatial MPGNNs are less effective as S 2 GNNs, for long-range interactions. This is evident for peptides Table 1 , clustering Fig. 10 , distance regression Fig. 11 , and over-squashing Fig. 12 . Specifically, if the task requires long-range interactions beyond the receptive field of MPGNNs, they return crude estimates. E.g., in Fig. 11 , the MPGNN predicts (approx.) constantly 20 for all distances beyond its receptive field -roughly the mean in the training data. Moreover, A Vdiag(ĝ ϑ (λ) : , 1 )V Vdiag(ĝ ϑ (λ) : , 2 )V Vdiag(ĝ ϑ (λ) : , 3 )V Vdiag(ĝ ϑ (λ) : , 4 )VFigure 13 : 4 filters on LR-CLUSTER. Large/small entries are yellow/blue, white lines mark clusters.S 2 GNNs may converge faster (see Fig. 25 in § M.6.2) and are more parameter-efficient, as we show on PCQM4Mv2 (Hu et al., 2021) in § M.9.Finding (III): virtual nodes are insufficient. We frequently find that including more than a single eigenvector (k > 1) yields substantial gains. We make this explicit in Fig. 10a , where we append a single spectral layer and sweep over the number of eigenvectors k. We complement these findings with an ablation for the frequency cutoff λ cut on peptides-func in § M.5.Finding (IV): our Positional Encodings PE consistently help, when concatenated to the node features. While this finding is true throughout our evaluation, the differences are more pronounced in certain situations. For example, on LR-CLUSTER in Fig. 10 , the PE help with spectral filter and a small k or without spectral filter and many message passing steps. Finding (V): spectral filters align with clusters, as we illustrate in Fig. 13 for four arbitrary spectral filters learned on LR-CLUSTER. We observe that (a) the spectral filters reflect the true clustering structure, (b) some filters are smooth while others contain details, and (c) they model coarser or finer cluster structures (e.g., first vs. third filter).Following the evaluation of Hyena (Poli et al., 2023) and H3 (Fu et al., 2023) , we benchmark S 2 GCN with sequence models on the associative recall in-context learning task, stemming from mechanistic interpretability (Elhage et al., 2021; Power et al., 2022; Zhang et al., 2023; Olsson et al., 2022) . In associative recall, the model is asked to retrieve the value for a key given in a sequence. For example, in the sequence a,0,e,b,z,9,h,2,=>,z, the target is the value for key z, which is 9 since it follows z in its prior occurrences. We create a sequence/path graph with a node for each "token" (separated by "," in the example above) and label the target node with its value. We assess the performance of S 2 GCN on graphs that vary in size by almost two orders of magnitude and follow Poli et al. (2023) with a vocabulary of 30 tokens. Moreover, we finetune our S 2 GCN on up to 30k nodes.Finding (VI): our spectral filter for directed are effective and may improve generalization, as we find in Fig. 14 (and Table 13 of § M.7).Finding (VII): S 2 GCN a state-of-the-art sequence model, as it performs on par with Hyena and, here, outperforms transformers (Table 2 ).Finding (VIII): S 2 GNNs is practical and scalable. We demonstrate this on the OGB Products graph (2.5 mio. nodes, Table 3 ) and the (directed) 10 million graphs dataset TPUGraphs (average number of nodes ≈10,000, Table 4 ). In both cases, we find full-graph training (without segment training (Cao et al., 2023) ) using 3 (Dir-) GCN layers interlayered with spectral filters, a reasonable configuration on a 40 GB A100. However, for OGB Products, we find that batching is superior, presumably because the training nodes are drawn from a "small" region of the graph (see § K). The cost of partial EVD for each dataset (excluding TPUGraphs and distance regression) is between 1 to 30 minutes on CPUs. We report the detailed costs of EVD and experiments in § M.3. 
related work: 
 Combining spatial and spectral filters has recently attracted attention outside of the graph domain in models like Hyena (Poli et al., 2023) Long-range interactions on graphs. Works that model long-range interactions can be categorized into: (a) MPGNNs on rewired graphs (Gasteiger et al., 2019a,b; Gutteridge et al., 2023) ; (b) higherorder GNNs (Fey et al., 2020; Wollschläger et al., 2024) that, e.g., may pass information to distant nodes through hierarchical message passing schemes; and (c) message passing adaptations to facilitate long-range interactions. For example, Park et al. (2023) propose "non-backtracking" message passing, Errica et al. (2024) adaptively choose the numbers of message passing steps, and Ding et al. (2024) use linear RNNs to aggregate over each node's neighborhoods. While approaches (a-c) can increase the receptive field of GNNs, they are typically still spatially bounded. In contrast, (d) alternative architectures, like graph transformers (Ma et al., 2023; Dwivedi & Bresson, 2021; Kreuzer et al., 2021; Rampášek et al., 2022; Geisler et al., 2023; Deng et al., 2024) with global attention, may model all possible n×n interactions. We provide notes on the limitations of graph transformers with absolute positional encodings in § D, which highlights the importance of capturing the relative relationships between nodes, as S 2 GNNs do. Moreover, in a recent/contemporary non-attention model for all pair-wise interactions, Batatia et al. (2024) use a resolvent parametrization of matrix functions relying on the LDL factorization of a matrix, but do not characterize their approximation-theoretic properties, over-squashing, expressivity on graphs, nor how to deal with directed graphs.In § B, we discuss additional related work w.r.t. expressivity and directed graphs. 
discussion: 
 We propose S 2 GNNs, adept at efficiently modeling complex long-range interactions via the synergistic composition of spatially and spectrally parametrized filters ( § 3). We show that S 2 GNNs share many properties with graph rewirings, pooling, and hierarchical message passing schemes (Fig. 3 & 4). S 2 GNNs outperform the aforementioned techniques with a substantial margin on the peptides long-range benchmark ( § 4.1), and we show that S 2 GNNs are also strong sequence models, performing on par or outperforming state-of-the-art like Hyena or H3 in our evaluation ( § 4.2). Even though we find global graph models, like S 2 GNNs, more prone to overfitting (see § K/L for further limitations/impact), moving to global models aligns with the trend for other deep learning domains. (column) Vector B Matrix B Set B ⊤ Transpose of matrix B B H Conjugate transpose of matrix B i, i 2 = -1 Complex number ⊙ Element-wise multiplication • Function composition, i.e., f 2 • f 1 = f 1 (f 2 (•)) ∥•∥ F Frobenius norm O"Big O" notation for asymptotic growth of function Bound on r-th derivative of spectral filter on interval [λ cut , 2] Spectral (l) (H (l-1) ; V , λ) Spectral(ly parametrized) filter Spatial (l) (H (l-1) ; A)P ∈ P Permutation matrix * G Graph convolution (see § 2) (G p ) p∈N Graph sequence (cf. Theorem 3) G(A, X) Graph A ∈ R n×n ≥0 Adjacency matrix X ∈ R n×d Node features n Number of nodes m Number of edges d Number of attributes L q ∈ C n×n Graph/combinatorial/random walk (Magnetic) Laplacian q ∈ R ≥0 Potential (see Eq. 7) L q=0 = L ∈ R n×n Real-valued Laplacian (q = 0) λ, V = EVD(L) Eigendecomposion s.t. L = V diag(λ)L H λ, V = EVD(L, k) Partial eigendecomposion containing λ 1 ≤ λ 2 ≤ • • • ≤ λ k (≤ λ k+1 ) kFilter in graph convolution g γp Order-p polynomial filter in graph convolution, with coefficients γ p ĝFilter in graph convolution, in spectral domain ĝ(λ)as function evaluated at eigenvalues ĝϑ (λ)as parametrized function evaluated at eigenvalues ĝ(l) ϑ (λ)as parametrized function at layer l, evaluated at eigenvaluesPE(V , λ) Positional encodings f θ Feature transformation in spatial domain s ζFeature transformation in spectral domainExpressivity. Laplacian eigenvectors have been used previously to construct positional encodings that improve the expressivity of GNNs or Transformers (Lim et al., 2023; Wang et al., 2022; Geisler et al., 2023; Huang et al., 2024) . Our positional encodings are similar to the preprocessing of Balcilar et al. (2021a) , where the authors design an edge-level encoding/mask to surpass 1-WL. The hierarchy of Weisfeiler-Leman (WL) tests is a common way to categorize the expressivity of GNNs (Grohe et al., 2021) . Xu et al. (2019) showed that most MPGNNs are bound by or as strong as the 1-WL test. Lim et al. (2023) point out that spectral GNNs suffer from similar limitations as MPGNNs w.r.t. their expressivity. Generally, the development of expressive GNNs is an active research direction, and we refer to Li & Leskovec (2022) for a broad overview.Directed graphs. Rossi et al. (2023) also extend the WL test to directed graphs and propose an MPGNN for directed graphs. How to model direction in graphs is also still an open question and various approaches were proposed (Battaglia et al., 2018; Tong et al., 2020; Zhang et al., 2021; Rossi et al., 2023; Koke & Cremers, 2024) . We utilize a Hermitian Laplacian for direction awareness, namely the Magnetic Laplacian, which was also used by Zhang et al. (2021) for an MPGNN and Geisler et al. (2023) for positional encodings.Undirected vs. directed graphs. For spatial filtering, it is straightforward to plausibly extend the message passing (e.g. Battaglia et al. (2018) ; Rossi et al. (2023) ). However, the spectral motivation and spectral filter on directed graphs require more care. The eigendecomposition is guaranteed to exist for real symmetric matrices. Real symmetric matrices are always diagonalizable, and the eigenvectors will then span a complete orthogonal basis to represent all possible signals X ∈ R n×d . Note that some non-symmetric square matrices are also diagonalizable and, thus, also have an eigendecomposition, albeit the eigenvectors may not be orthogonal. Thus, further consideration is required to generalize the graph Laplacian to general directed graphs.Magnetic Laplacian. For the spectral filter on directed graphs, we build upon a direction-aware generalization, called Magnetic Laplacian (Forman, 1993; Shubin, 1994; De Verdière, 2013; Furutani et al., 2020; Geisler et al., 2023 )EQUATIONwhere A s = A ∨ A ⊤ is the symmetrized graph with diagonal degree matrix D s . ⊙ denotes the element-wise product, exp the element-wise exponential, i = √ -1 an imaginary number, and q the potential (hyperparameter). By construction L q is a Hermitian matrix L q = L H q where the conjugate transpose is equal to L q itself. Importantly, Hermitian matrices naturally generalize real symmetric matrices and have a well-defined eigendecomposition L q = V ΛV H with real eigenvalues Λ and unitary eigenvectors V V H = I. For appropriate choices of the potential q, the order of eigenvalues is well-behaved (Furutani et al., 2020) . Recently Geisler et al. (2023) demonstrated the efficacy of these eigenvectors for positional encodings for transformers. Moreover, the Magnetic Laplacian was used for a spectrally designed spatial MPGNN (Zhang et al., 2021) , extending Defferrard et al. (2017) . Due to the real eigenvalues, one could, in principle, also apply a monomial basis (Chien et al., 2021) , or different polynomial bases stemming from approximation theory (He et al., 2021; Wang & Zhang, 2022; He et al., 2022; Guo & Wei, 2023) .To see why Eq. 7 describes an injection for appropriate choices of q, consider that the sparsity pattern of L q matches A up to the main diagonal. If A contains a self-loop the main diagonal will have a 0 instead of 1 entry at the self-loop location. A -A ⊤ can be directly inferred from the phase exp[i2πq(A -A ⊤ )], assuming that q < 1 /(2 maxu,v Au,v). Thus, it is solely left to obtain A s from I -D -1 /2 s A s D -1 /2 s, which is trivial for a binary adjacency but more involved for real-valued weights. Determining if and when L q is injective for real-valued A is left for future work.Properties of the eigendecomposition. The eigendecomposition is not unique, and thus, one should consider the result of the eigensolver arbitrary in that regard. One ambiguity becomes apparent from the definition of an eigenvalue itself Lv = λv since one can multiply both sides of the equation with a scalar c ∈ C \ {0}: L(cv) = λ(cv). We already implicitly normalized the magnitude of the eigenvectors V by choosing them to be orthogonal (V V ⊤ = I) or unitary (V V H = I). Thus, after this normalization, c only represents an arbitrary sign for real-valued eigenvectors or a rotation on the unit circle in the complex case. Another reason for ambiguity occurs in the case of repeated / multiple eigenvalues (e.g., λ u = λ v for u ̸ = v). In this case, the eigensolver may return an arbitrary set of orthogonal eigenvectors chosen from the corresponding eigenspace.Here, we consider a vanilla graph transformer f (X) that solely becomes structure-aware due to the addition (or concatenation) of positional encodings: f (X + PE(A)). The main point we are going to demonstrate is that a vanilla transformer with such absolute positional encodings PE(A) ∈ R n×d will be limited in its expressivity if the positional encodings are permutation equivariant P PE(A) = PE(P AP ⊤ ) w.r.t. any n × n permutation matrix P ∈ P.The limitation particularly arises in the presence of automorphisms P a AP ⊤ a = A with specifically chosen permutation P a . To be more specific, assume that nodes u and v are automorphic to each other. That is, there exists a P a that will swap the order of u and v (among other nodes) s.t. P a AP ⊤ a = A. By permutation equivariance, we knowP a PE(A) = PE(P a AP ⊤ a ) = PE(A) and, hence, PE(A) u = PE(A) v .We have just shown that automorphic nodes will have the same positional encodings PE if the positional encodings are permutation equivariant. This implies that permutation equivariant positional encodings PE are not even able to capture simple neighboring relationships. For example, consider an undirected sequence/path graph o-o-o-o-o with five nodes. Here, the two end nodes, which we also all first and last node, are automorphic. So are the second and second-last nodes. Assuming the second and second last nodes have different node features (e.g., A-B-C-D-A), that breaks the symmetry, it is still not possible for a transformer with absolute positional encodings to tell the first and last node apart. In other words, in the example, the transformer cannot tell apart the end node with neighboring feature B from the end node with neighboring feature D. This shows a severe limitation of architectures without additional components capturing the relative distances (e.g., as S 2 GNNs can). This concern is not as critical for architectures where the positional encodings are not entirely permutation equivariant (Dwivedi & Bresson, 2021; Kreuzer et al., 2021) , with relative positional encodings (Ma et al., 2023) , and might also be of lesser concern for directed graphs (Geisler et al., 2023) .Adding a fully connected virtual node (Gilmer et al., 2017) is among the simplest ways to add the ability for long-range information exchange. An equivalent method was proposed as a simple oversquashing remedy in the seminal work by Alon & Yahav (2020) . A single Spectral layer amounts to a type of virtual nodes in the special case of f θ = I andEQUATIONAssuming a simply-connected graph G, the unique normalized zero-eigenvector v of the symmetrically-normalized graph LaplacianL = I -D -1 /2 AD -1 /2 has components v u = du 2|E|, where d u denotes the degree of node u ∈ G, and |E| the number of edges in the graph. At node u ∈ G, we therefore findSpectral (l) u (H (l-1) ; V , λ) = √ d u 2|E| v∈G d v h (l-1) v (9) with h (l-1) vdenoting the row of H (l-1) corresponding to node v ∈ G. In other words, filtering out the zero-frequency component of the signal means scattering a global, degree-weighted embedding average to all nodes of the graph. For the unnormalized graph Laplacian, Eq. 9 instead becomes an unweighted average, which is consistent with the usual definition of a virtual node. We refer to Fig. 3 & 4 for additional intuition.We restate two key results from Di Giovanni et al. (2023a) using our notation. They imply the existance of a regime in which 1-hop MPNN architectures suffer from exponentially decaying Jacobian sensitivity. Meanwhile, S 2 GNNs can easily learn a signal of constantly lower-bounded sensitivity, as shown by invoking its trivial subcase of a virtual node in Theorem 2.Theorem 7 (Adapted from Di Giovanni et al. (2023a) ). In an ℓ-layer spatial MPGNN with messagepassing matrix S = c r I + c a A (c r , c a ∈ R + ) and a Lipschitz nonlinearity σ,H (l) = Spatial (l) (H (l-1) ; A) = σ SH (l-1) W (l-1) , 1 ≤ l ≤ ℓ (10)the Jacobian sensitivity satisfies the following upper bound:EQUATIONwith hEQUATIONv denoting the rows of H (0) , H (ℓ) corresponding to the nodes v, u ∈ G, c σ the Lipschitz constant of the nonlinearity, w the maximum entry value over all weight matrices W (l) , and d the network width.The dependence of the upper bound on the matrix power S ℓ vu -not generally present for S 2 GNN by Theorem 2 -leads to a topology-dependence which becomes explicit in the following theorem. It concerns the typical shallow-diameter regime, in which the number ℓ of MPGNN layers is comparable to the graph diameter. Theorem 8 (Adapted from Di Giovanni et al. ( 2023a)). Given an MPNN as in Eq. 10, with c a ≤ 1, let v, u ∈ G be at distance r. Let c σ be the Lipschitz constant of σ, w the maximal entry-value overall weight matrices, d min the minimal degree of G, and γ ℓ (v, u) the number of walks from v to u of maximal length ℓ. For any 0 ≤ k < r, there exists C k > 0 independent of r and of the graph, such that ∂h(r+k) v ∂h (0) u L 1 ≤ C k γ r+k (v, u) 2c σ wd d min r .For 1-hop MPGNNs with 2c σ wd < d min , we therefore identify an exponential decay of sensitivity with node distance r in the weak-connectivity limit for which γ r+k (v, u) increases subexponentially with r. As Di Giovanni et al. (2023a) point out, sharper bounds can be derived under graph-specific information about (S r ) vu .We express the electric potential along a periodic sequence of screened 1D charges as a convolution of a corresponding graph signal with a consistently defined graph filter. This closed-form example underscores our default use of a low-pass window for the spectral part of S 2 GNNs by showing how a continuous problem with a convolutional structure and quickly flattening spectral response (typical for pair interactions in physics and chemistry) discretizes into a graph problem with similar features.The approach exploits the surjective mapping of Fourier modes on [0, n] onto the Laplacian eigenvectors of a cycle graph C n . We consider two corresponding representations of the same problem:EQUATIONEQUATION• A continuous representation (Eq. 12) in terms of a 1D distribution ρ of n point charges q 1 , . . . , q n and their periodically repeating image charges, written as a sum of Dirac combs at equidistant offsets l with 0 ≤ l ≤ n -1, interacting via the potential profile ϕ σ obtained from solving in Gauss' law of electrostatics for a 1D point charge screened by a Gaussian cloud of opposite background charge with standard deviation σ. The screening ensures convergence to a finite potential and its exact form is insignificant (we choose the Gaussiantype screening due to its analytical tractability). Note that ϕ σ (x) ≃ const. -|x| for x → 0 (the unscreened 1D potential in the direction normal to an infinitely wide capacitor plate), while the screening guarantees an exponential dropoff to zero as x → ∞,• A graph representation (Eq. 14) by placing the n charges q 1 , . . . , q n onto a cycle graph C n .We derive the graph filter g σ from a consistency condition (Eq. 13) between both representations: the graph convolution (g σ * G q) has to yield the electric potential V sampled at the charge loci if we want g σ to act like the continuous convolution kernel ϕ σ in the discrete graph picture.The Fourier transform of ϕ σ (in the convention without integral prefactor and with a 2π frequency prefactor) reads φσ (κ) = 1 πκ 2 1 -exp -1 2 σ 2 κ 2 . For the density, the Poisson sum formula gives ρ(κ) = n-1 k=0 1 √ n qk ∆ 1 (κ -k n ) with qk = 1 √ n n-1 j=0 q i exp -i2π k n j .The coefficients qk are precisely the components of the graph Fourier transform of q (physically, they amount for the structure factor). By the convolution theorem, V (κ) = φσ (κ)ρ(κ). By noting that all integer-shifted frequencies in the Dirac combs ∆ 1 • -k n (or all Brillouin zones, in physics terminology) yield the same phase exp i2π k n l if we only sample V (x) at the charge locix = l, 0 ≤ l ≤ n -1, we can write V (l) = 1 2π n-1 k=0 qk m∈Z φσ k n + m 1 √ n exp i2π k n l .Through pattern-matching with the consistency condition of Eq. 13, we can therefore identify that the graph filter is a sum over Brillouin zones,(ĝ σ )(λ k ) = 1 2π m∈Z φσ k n + m , where λ k denotes the eigenvalues of the normalized C n graph Laplacian, λ k = 1 -cos 2πkn . To fulfill this relation for all n, k we setĝσ (λ) = 1 2π m∈Z φσ 1 2π arccos(1 -λ) + mWe claim now (and prove in a later paragraph) that for λ > λ 0 > 0 and a sufficiently large choice σ > σ(r, λ 0 ), the absolute r-th derivative satisfies the upper bound| d r dλ r ĝσ (λ)| ≤ | d r dλ r ĝ∞ (λ)|,where we can think of ĝ∞ as the limit of taking σ → ∞ (i.e., a constant background charge):ĝ∞ (λ) = 1 2π m∈Z φ∞ 1 2π arccos(1 -λ) + m , φ∞ (κ) = 1 πκ 2The merit of this is that unlike the screened ĝσ (λ), ĝ∞ (λ) can be solved analytically to find closedform bounds on the absolute derivatives | d r dλ r ĝσ (λ 0 )|. By invoking the sum expansion form of the trigamma functionΨ 1 (z) = ∞ m=0 1 (z+m) 2 , the reflection identity ψ 1 (1 -z) + ψ 1 (z) = π 2sin 2 πz , and the half-angle formula sin 2 x 2 = 1-cos(x), we findĝ∞ (λ) = 1 2π 2 Ψ 1 1 2π arccos(1 -λ) + Ψ 1 1 - 1 2π arccos(1 -λ) = 1 2 sin 2 1 2 arccos(1 -λ) = 1 λ ,a remarkably simple result. We can now readily evaluate | d r dλ r ĝ∞ (λ)| = r! λ r+1 , but it remains to prove that this upper-bounds | d r dλ r ĝσ (λ)| for any λ > λ 0 > 0 and sufficiently large σ > σ(r, λ 0 ). For compactness, define the expressions z(λ) := 1 2π arccos(1 -λ) ∈ 0, 1 2 (strictly increasing in λ), y σ (z) := exp -1 2 σ 2 z 2 , and z = 1 -z ≥ z. Consider the series of "term-by-term" derivativesd dz ĝσ (λ(z)) = - 1 π 2 ∞ m=0 1 (z + n) 3 (1 -y σ (z + m)) - 1 (z + n) 3 (1 -y σ (z + m)) + ∞ m=0 O (y σ (m)) d 2 dz 2 ĝσ (λ(z)) = 3 π 2 ∞ m=0 1 (z + n) 4 (1 -y σ (z + m)) + 1 (z + n) 4 (1 -y σ (z + m)) + ∞ m=0 O (y σ (m)). . . They converge uniformly on 0, 1 2 as they clearly are Cauchy sequences under uniform bound (moreover, well-definedness in z = 0 follows by applying l'Hospital's rule -physically, this is the merit provided by including Gaussian screening in our model). Therefore, they indeed converge to the respective derivatives d r dz r ĝσ (λ(z)) (justifying the above notation). The same holds for the corresponding series for d r dz r ĝ∞ (λ(z)): they are not defined in z = 0, but otherwise still converge as they match the known series expansion of the polygamma function. Given λ 0 > 0 and thus z(λ 0 ) > 0, taking σ larger than some σ(r, λ 0 ) guarantees that d r dz r ĝσ (λ(z)) and d r dz r ĝ∞ (λ(z)) are of the same sign for λ > λ 0 (z(λ) > z(λ 0 )). This holds for all orders r ∈ N since we see by induction that the product rule always yields one term analogous to the first respective terms above, and otherwise only terms of O (y σ (m))). Then, observing that 0≤ 1 -y σ (x) < 1 ∀ x ≥ 0 and z ≥ z implies | d r dz r ĝσ (λ 0 (z 0 ))| ≤ | d r dz r ĝ∞ (λ 0 (z 0 ))|.The same must hold for the λ-derivatives by the chain rule.One interesting question is whether ĝσ is also C-integral-Lipschitz for some constant C > 0. We discuss this stability-informed criterion (Gama et al., 2020) in the main body as a domain-agnostic prior assumption about the "ideal" graph filter if no other ground truth knowledge informing additional smoothness bounds (such as here) is available. While the above bound is too loose to certify this directly (| d dλ ĝσ (λ)| ≤ Cλ -1 would be needed), integral-Lipschitzness under some constant follows from the fact that | d dλ ĝσ (λ)| is bounded on [0, 2]: by the uniform convergence of the term-by-term derivative series, it is continuous. Well-definedness of the product d dz ĝσ dz dλ has to be checked in λ = 0, where it follows by continuous extension using l'Hospital's rule. As a continuous function defined on a compact interval, | d dλ ĝσ | assumes a maximum.We next prove the permutation equivariance of the spectral filter in Eq. 3:Theorem 1. Spectral(H (l-1) ; EVD(L, k)) of Eq. 3 is equivariant to all n × n permutation matrices P ∈ P: Spectral(P H (l-1) ; EVD(P LP ⊤ , k)) = P Spectral(H (l-1) ; EVD(L, k)).for the general case of parametrizing a Hermitian "Laplacian" L ∈ C n×n , L H = L. Note that this proof does not rely in any means on the specifics of L, solely that the eigendecomposition exists L = V ΛV H with unitary eigenvectors V V H = I. For practical reasons, it is suitable to define L(A) as a function of A. A similar proof for real-valued eigenvectors is given by (Lim et al., 2023) .The specific spectral filter we consider isEQUATIONwith arbitrary f : C d1 → C d2 , applied row-wise to H (l-1) ∈ C n×d1 . Analogously, h : C d2 → C d3 is applied row-wise. We choose complex functions to emphasize generality, although we restrict Spectral to real in-and outputs in all experiments. The graph filter is defined as element-wise function ĝu,v (λ) := ĝv (λ u , {λ 1 , λ 2 , . . . , λ k }) that depends on the specific eigenvalue λ and potentially the set of eigenvalues {λ 1 , λ 2 , . . . , λ k } (or its vector representation λ) of the partial eigendecomposition.We need to make sure that the partial decomposition includes all eigenvalues of the same magnitude, i.e., λ u ̸ = λ u ′ , ∀u ∈ {1, 2, . . . , k}, u ′ ∈ {k + 1, k + 2, . . . , n}. In practice, this is achieved by choosing large enough k to accommodate all eigenvalues λ cut < λ k+1 , or by dropping trailing eigenvalues where λ j = λ k+1 for j ∈ {1, 2, . . . , k}. Generally, it is also not important that we consider the k smallest eigenvalues in the spectral filter. We only need to ensure that the spectral filter is either calculated on all or no eigenvalues/-vectors of an eigenspace.Proof. Assuming functions ϕ(X) and ψ(X) are permutation equivariant, then ϕ(ψ(X)) is permutation equivariant ϕ(ψ(P X)) = ϕ(P ψ(X)) = P ϕ(ψ(X)) for any n × n permutation P ∈ P. Thus, it sufficies to prove permutation equivariance for h, f,V (ĝ(λ) ⊙ [V H X]) independently, where X ∈ C n×d2 .Regardless of the complex image and domain of h and f , they are permuation equivariant since they are applied row-wisef (X) = [f (X 1 ) f (X 2 ) . . . f (X n )]H and reordering the rows in X ∈ C n×d1 also reorders the outputs: f (P X) = P f (X).For finalizing the proof of permutation equivariance, we first rearrangeY = V (ĝ(λ) ⊙ [V H X]) = k u=1 v u (ĝ u,: (λ u ) ⊙ [v H u X]) and Y :,v = k u=1 ĝu,v (λ u )v u v H u X :,v .This construction (a) is invariant to the ambiguity that every eigenvector v u can be arbitrarily rotatedc u v u by {c u ∈ C | |c u | = 1}. That is, (c u v u )(c u v u ) H = c u cu v u v H u = v u v H u .Moreover, (b) in the case of j repeated eigenvalues {s + 1, s + 2, . . . , s + j} where λ s+1 = λ s+2 = • • • = λ s+j , we can choose a set of orthogonal eigenvectors arbitrarily rotated/reflected from the j-dimensional eigenspace (basis symmetry). The given set of eigenvectors can be arbitrarily transformed V :,s+1:s+j Γ j by a matrix chosen from the unitary group Γ j ∈ U (j). Sinces+j u=s ĝu,v (λ u )v u v H u X :,v = ĝs,v (λ s ) s+j u=s v u v H u X :,v = ĝs,v (λ s ) V :,s+1:s+j V H :,s+1:s+j X :,vwe simply need to show that the expression is invariant to a transformation by Γ j :V :,s+1:s+j Γ j (V :,s+1:s+j Γ j ) H = V :,s+1:s+j Γ j Γ H j V H :,s+1:s+j = V :,s+1:s+j V H :,s+1:s+jTo see why Γ j ∈ U (j) is a sufficient choice in the light of repeated/multiple eigenvalues, consider the defintion of eigenvalues/vectorsLV :,s+1:s+j = L | | | v s+1 v s+2 . . . v s+j | | | = | | | v s+1 v s+2 . . . v s+j | | |     λ s+1 0 . . . 0 0 λ s+2 . . . 0 . . . . . . . . . . . . 0 0 . . . λ s+j     = λ s+1 | | | v s+1 v s+2 . . . v s+j | | | = λ s+1 V :,s+1:s+jwe can now multiply both sides from the right with an arbitrary matrix B ∈ C j×j . To preserve the unitary property V :,s+1:s+j V H :,s+1:s+j = I, we require (V :,s+1:s+j B)(V :,s+1:s+j B) H = I. Thus, the eigenvectors can be arbitrarily transformed byΓ j ∈ U (j) instead of B ∈ C j×j .This concludes the proof.We restate Theorem 2 in more detail and also considering graphs that contain multiple connected components. The unchanged bottom line is that S 2 GNNs can express signals lower-bounded by a constant that is unaffected by local properties of the graph topology, instead of suffering from exponential sensitivity decay like spatial MPGNNs. Theorem (Theorem 2, formal). Consider an ℓ-layer S 2 GNN of the form Eq. 1. Let ( θ, ϑ, θ) be parameters of the spatial GNN, spectral filters ĝ(l) ϑ , and feature transformation f θ . Assume the existence of parameters θ such that Spatial (l) (H (l-1) ; A, θ) = 0 ∀1 ≤ l ≤ ℓ and θ such that f θ = I. Then, a filter choice ϑ exists such that the ℓ-layer S 2 GNN of the additive form Eq. 1 can express a signal h (ℓ) v (H (0) ; θ, ϑ, θ) with uniformly lower-bounded Jacobian sensitivity,EQUATIONwith hEQUATIONv denoting the rows ofH (0) , H (ℓ) corresponding to the nodes u ̸ = v ∈ G, connected component C ⊂ G containing |E C | edges, network width d and parameter-dependent constant K ϑ .Proof. Choose θ such that Spatial (l) (H (l-1) ; A) = 0 ∀1 ≤ l ≤ ℓ (typically by setting all weights and biases to zero), θ such that f θ = I, and set ϑ such thatEQUATIONfor some K ϑ > 0. This choice of filter parameters ϑ lets Spectral act like a type of virtual node across all hidden dimensions k: In the standard orthonormal basis of the 0-eigenspace given byEQUATIONwhere C enumerates all connected components, and d u denotes the degree of node u, we findEQUATIONwith C (v) denoting the connected component containing v. Particularly, note that applying the spectral layer more than once does not affect the result since the projector onto an eigenvector is idempotent (up to K ϑ ). The result must also hold in any other orthonormal basis of the 0-eigenspace due to the invariance of Spectral under orthogonal eigenbasis transformations. Differentiating with respect to h (0) u , taking the L 1 norm and using√ d u d v ≥ 1 shows the statement.Theorem 3. Let ĝ be a discontinuous spectral filter. For any approximating sequence g γp p∈N of polynomial filters, an adversarial sequence (G p ) p∈N of input graphs exists such that∄α ∈ R >0 : sup 0̸ =X∈R |Gp|×d ∥(g γp -g) * Gp X∥ F ∥X∥ F = O p -αThe proof makes use of a result by S. Bernstein (Natanson, 1964) : Theorem 9 (Bernstein). Let f : [0, 2π] → C be a 2π-periodic function. Then f is α-Hölder continuous for some α ∈ (0, 1) if, for every p ∈ N, there exists a degree-p trigonometric polynomialT p (x) = a 0 + p j=1 a j cos(jx) + p j=1 b j sin(jx) with coefficients a j , b j ∈ C, such that sup 0≤x≤2π |f (x) -T p (x)| ≤ C(f ) p αwhere C(f) is a positive number depending on f .Proof. Given a discontinuous filter ĝ : [0, 2] → R, construct the function f : [0, 2π] → C fulfilling the prerequisites of Theorem 9 by pre-composing f := ĝ • (cos(•) + 1). We proceed via contradiction. Suppose that there is an α ∈ (0, 1) and a sequence of degree-p polynomial filters, ĝγp (λ) = p j=0 γ j λ j , γ = (γ 0 , . . . , γ p ) ⊤ ∈ R p+1 , such that ∥ĝ γp -ĝ∥ ∞ = O(p -α ). Then, the sequence of trigonometric polynomials T p := ĝγp • (cos(•) + 1) fulfills the condition of Theorem 9. This would imply that f = ĝ • (cos(•) + 1) is α-Hölder continuous, meaning that a constant K > 0 exists such that|ĝ(cos(x) + 1) -ĝ(cos(y) + 1)| ≤ K|x -y| α ∀ x, y ∈ [0, 2π]Considering λ 0 ∈ [0, 2], λ → λ 0 and x = arccos(λ 0 -1), y = arccos(λ -1) (using the arccos branch in which both λ 0 , λ eventually end up) shows a contradiction to the assumed discontinuity of ĝ. Therefore, no polynomial filter sequence ĝγp p∈N together with an α ∈ (0, 1) exist such that ∥ĝ γp -ĝ∥ ∞ = O(p -α ). In particular, for any sequence ĝγp p∈N , a sequence of adversarial values(λ p ) p∈N , λ p ∈ [0, 2] exists such that ∄α ∈ (0, 1) : |ĝ γp (λ p ) -ĝ(λ p )| = O(p -α )The proof is finished if we can find a sequence of graphs (G p ) such that the symmetricallynormalized graph Laplacian L p of G p contains λ p as an eigenvalue. In this case, we can construct adversarial input signals X p on the graphs G p by setting the first embedding channel to an eigenvector corresponding to λ p , and the remaining channels to zero, such that g γp -g * Gp X p = |ĝ γp (λ p ) -ĝ(λ p )|X p . In particular, it then holds that∄α ∈ R + : sup 0̸ =X∈R |Gp|×d ∥(g γp -g) * Gp X∥ F ∥X∥ F = O p -αIf we assume only simple graphs, such a construction is unfortunately not possible since the set of all simple graphs and therefore the set of all realizable eigenvalues is countable, whereas the adversarial values λ p could lie anywhere in the uncountable set [0, 2]. We can, however realize arbitrary eigenvalues by using weighted graphs with three nodes. Consider a cyclic graph structure and tune the weight of edge (1, 2) to sin 2 (θ p ) and the weight of edges (2, 3) and (3, 1) to cos 2 (θ p ) with θ p ∈ 0, π 2 . The symmetrically-normalized graph Laplacian, Remark. If one wishes to restrict the set of possible adversarial graph sequences (G p ) p∈N to include only simple graphs, a version of Theorem 3 still holds where we restrict the assumption to filters ĝ which are piecewise-continuous with discontinuities on a finite set of points D ⊂ S, where S ⊂ [0, 2] denotes the countable set of eigenvalues realizable by simple graphs. This still covers a large class of filters to which order-p polynomial filters can provably converge slower than any inverse root of p in the operator norm, and includes the virtual node filter (discontinuous only in λ = 0) presented as an example in the main body. The proof is fully analogous up to the point of constructing λ p . If λ p ∈ D, we can find a graph that realizes it exactly. Now assume λ p / ∈ D. We note that the set S is dense in [0, 2] (clear from considering, e.g., the cyclic graphs C n with symmetricallynormalized Laplacian eigenvalues λ k = 1 -cos 2πk n ). Since we assume that ĝ and therefore also |ĝ γp -ĝ| is piecewise-continuous anywhere but on D ⊂ S and D is finite, we can find an open neighborhood N (λ p ) for any λ p / ∈ D on which ĝ is continuous. Using that S is dense in [0, 2], we find a graph sequence G(l)EQUATIONEQUATIONEQUATIONp l∈N → λ p for which∥ĝ γp ( λ(l) p ) -ĝ( λ(l) p )∥ → ∥ĝ γp (λ p ) -ĝ(λ p )∥.Therefore, by the same reasoning as in the proof of Theorem 3, we find that there can be no α ∈ (0, 1) for which sup 0̸ =X∈R |Gp |×d∥(gγ p -g) * Gp X∥F ∥X∥F is of O (p -α ).We first introduce the setting and notation to state Theorem 4 in its general version. We study how well S 2 GNNs can approximate "idealized" GNNs (IGNNs) containing L graph convolution layers 1 ≤ l ≤ L, each of which can express a convolution operator g with any spectral representation ĝ(l) : [0, 2] → R d (l). An IGNN layer therefore has the structureEQUATIONwithH (l) ∈ R n×D (l) , W (l) ∈ R D (l) ×D (l-1) and V ∈ R n×n .We compare this to S 2 GNNs with ℓ = (m + 1)L layers for m ≥ 1, in the additive form of Eq. 1,EQUATION) Each layer 1 ≤ l ≤ ℓ parametrizes a spatio-spectral convolution. The spectral part satisfies Eq. 3,EQUATION)with embeddingsH (l) ∈ R n×d (l), linear feature transforms f(l) θ := W (l) spec ∈ R d (l) ×d (l-1)and a spectral filter ĝ(l) ϑ : [0, 2] → R that is fully supported and a universal approximator on [0, λ cut ]. Note we assume here that in every layer, there is only one spectral filter which gets reshaped as to act on every hidden component, whereas in practice, we relax this assumption to different filters per component, which can only be more expressive. The spatial part is a polynomial filter of the formSpatial (l) (H (l-1) ; A) = σ     p j=0 γ (l) j L j   H (l-1) W (l) spat   = σ V ĝ(l) γ (λ) ⊙ V ⊤ H (l-1) W (l) spat with W (l) spat ∈ R d (l) ×d (l-1), polynomial order p (fixed across layers), and a spectral representationEQUATIONp ) ⊤ ∈ R p+1 . We note that theorem 4 extends immediately to the case of directed graphs if the spatial part is instead a polynomial of the magnetic Laplacian (see section 3.2.3) over complex-valued embeddings like in Zhang et al. (2021) .Note that the layer-wise hidden dimensions D (l) vs. d (l) of the IGNN vs. S 2 GNN do not have to agree except at the input layer, d (0) = D (0) (of course, both networks receive the same input H (0) = H (0) = X), and at the output layer, d (ℓ) = D (L) . We now state the general version of Theorem 4. Theorem (Theorem 4, general). Assume an L-layer IGNN with filters ĝ(l) such that ĝ(l)[λcut,2] ∈ C r [λ cut , 2] and d r dλ r ĝ(l) [λcut,2] ∞ ≤ K max r (λ cut ) for all 1 ≤ l ≤ L. Let ∥ĝ (l) ∥ ∞ ≤ ∥ĝ∥ max ∞ and ∥W (l) ∥ 2 ≤ ∥W ∥ max 2 for all 1 ≤ l ≤ L. Assume that σ = [ • ] ≥ is the ReLu function. Then,(1) For a fixed polynomial order p ≥ 2, an approximating sequence S 2 GNN m m∈N of [(m + 1)L]layer S 2 GNNs exists such that, for arbitrary graph sequences(G m ) m∈N , sup 0̸ =X∈R |Gp |×d ∥ (S 2 GNN m ) Gm -(IGNN) Gm (X)∥ F ∥X∥ F = O C L (∥ĝ∥ max ∞ , ∥W ∥ max 2 ) K max r (λ cut ) (pm) -r , C L (∥ĝ∥ max ∞ , ∥W ∥ max 2 ) = ∥W ∥ max 2 L-1 l=1 ∥ĝ∥ max ∞ ∥W ∥ max 2 + (∥ĝ∥ max ∞ ∥W ∥ max 2 ) lwith a leading-order scaling constant that depends only on r. Here, ( • ) Gm denotes the instantiation of all model filters on the eigenvalues of an input graph G m , which maps both models onto a G mdependent function R D (0) → R D (L) .(2) For fixed m ≥ 1, an approximating sequence S 2 GNN p p∈N of [(m + 1)L]-layer S 2 GNNs with increasing layer-wise polynomial order p exists such that, for all (G p ) p∈N , the same bound holds.Proof. We first prove the following lemma, narrowing down the previous theorem to a single layer.Lemma 1. Let IGNN (l) denote a single IGNN layer as in Eq. 20, with a filter ĝ(l) such that ĝ(l) [λcut,2] is r-times continuously differentiable on [λ cut , 2] and satisfies a bound K r ĝ(l) , λ cut ≥ 0,d r dλ r ĝ(l) (λ) ≤ K r ĝ(l) , λ cut ∀λ ∈ [λ cut , 2]. Let σ = [ • ]≥ be the ReLu function, and let ∥W (l) ∥ 2 denote the spectral norm of W (l) . Then, (1) For fixed polynomial order p ≥ 2, an approximating sequence S 2 GNN (l) m m∈N of (m+1)-layer S 2 GNNs exists such that, for arbitrary graph sequences (G m ) m∈N , sup0̸ =X∈R |Gp |×d (S 2 GNN (l) m ) Gm -(IGNN (l) ) Gm (X) F ∥X∥ F = O [∥W (l) ∥ 2 K r (ĝ, λ cut )](pm) -rwith a scaling constant that depends only on r. Here, ( • ) Gm denotes the instantiation of all model filters on the eigenvalues of an input graph G m , which maps both models onto a G m -dependent function R D (l-1) → R D (l) .(2) For fixed m ≥ 1, an approximating sequence S 2 GNN (l) p p∈N of (m + 1)-layer S 2 GNNs with increasing layer-wise polynomial order p exists such that, for all (G p ) p∈N , the same bound holds.Remark. The proof of the simplified Theorem 4 used in the main body is analogous to the proof of Lemma 1 just without the nonlinearity, which has the following consequences:• The final layer m + 1 which we only need to apply one last nonlinearity to the output (since the spectral part of all layers, including the previous layer m, has none) becomes obsolete, so the final layer instead becomes m,• The two limits (1) and ( 2) are equivalent by the reduction to an mp-order polynomial filter,• We do not need the dimension-doubling "trick" outlined below to get rid of the nonlinearity in the proof and instead set all feature transform matrices in layers 1 through m -1 to the identity and the final ones to W * (m)spec = W (l) , W * (m) spat = W (l) .Proof of Lemma 1. We first note that m S 2 GNN spatial parts, each of order p, would act like an (mp)-order polynomial filter (factorized into m order-p polynomials), were it not for the nonlinearities in between. However, using the fact that σ is the ReLu function, we can choose intermediate hidden dimensions twice the size of the input dimension and then use the linear transforms to store a positive and a negative copy of the embeddings, add them back together after applying each ReLu, just to split the result back into a positive and negative copy for the next layer. This essentially gets us rid of σ. Throughout the proof, we use a star superscript to denote the specific parameters that will ultimately satisfy our bound, whereas we put no star above parameters that are yet to be fixed in a later part of the proof.For m ≥ 2, the trick discussed above works if we setW * (1) spec = 1 2 (I -I) ∈ R D (l-1) ×2D (l-1) , W * (2) spec , . . . , W * (m-1) spec = 1 2 I -I (I -I) ∈ R 2D (l-1) ×2D (l-1) , W * (m+1) spec = W (l) I -I ∈ R 2D (l-1) ×D (l) , W * (1) spat = (I -I) ∈ R 2D (l-1) ×D (l-1) , W * (2) spat , . . . , W * (m-1) spat = I -I (I -I) ∈ R 2D (l-1) ×2D (l-1) , W * (m+1) spat = W (l) I -I ∈ R 2D (l-1) ×D (l) .In the case m = 1, pick the matrices W * (1) spec , W * (1)spat from above for the first, and the matrices W * (m+1) spec, W * (m+1) spat from above for the second layer.Set ĝ * (m+1) γ (λ) = 1 and ĝ * (m+1) ϑ (λ) = 0. Given these choices and a graph G with eigenvalues λ,(S 2 GNN (l) ) G (X) = σ V ĝspsp (λ) ⊙ V ⊤ H (l-1) W (l) , ĝspsp = m j=1 ĝ(j) ϑ + ĝ(j) γ We see that ĝspsp [λmax,2] = m j=1 ĝ(j) γ since ĝ(j) ϑ [λmax,2] = 0 for 1 ≤ j ≤ m.This can express any polynomial up to order mp on [λ max , 2], since we assumed a layer-wise p ≥ 2 and any polynomial with real coefficients factorizes into real-coefficient polynomials of degree less or equal to 2 by the fundamental theorem of algebra. On the interval [0, λ max ], on the other hand, the filter ĝspsp [0,λmax] can express any IGNN filter ĝ(l) [0,λmax] . For m = 1, this is immediately clear. Else, set ĝ(j)ϑ to constants C j ∈ R ≥ , 1 ≤ j ≤ m -1 large enough that none of the polynomials C j + ĝ(j) γ , 1 ≤ j ≤ m -1, has a zero in [0, λ max ]. Defining ĝ(m) ϑ = ĝ(l) [0,λmax ] m j=1 Cj +ĝ (j) γ - ĝ(m) γ [0,λmax]gives the desired function.We proceed by making use of a result by D. Jackson (Natanson, 1964) , which is essentially a converse to Theorem 9 which we used to prove Theorem 3:Theorem (Jackson's theorem on an interval). Let a < b ∈ R, k, r ∈ N with k ≥ r -1 ≥ 0, f ∈ C r [a, b].Then, a polynomial p k of degree less or equal to k exists such that∥p k -f ∥ ∞ ≤ b -a 2 π 2 r 1 (k + 1)k . . . (k -r + 2) d r dx r f ∞Since ĝspsp can express any polynomial up to order mp on [λ max , 2] and, for any such polynomial, find parameters for the spectral parts that match the ideal filter ĝ(l)[0,λmax] exactly (not contributing to the supremum error), we can directly transfer this theorem to our case. Define S 2 GNN (l) m from the lemma by setting the linear feature transforms and final-layer filters as above. For the filters in layers 1 through m, define γ * (1) , . . . , γ * (m) such that m j=1 ĝ * (j) γ factorizes into into the polynomial from Jackson's theorem on [λ max , 2], and ϑ * (1) , . . . , ϑ * (m) to match ĝ(l) on [0, λ max ]. This defines a filter ĝ(l) spsp . We then find, for mp ≥ r -1 ≥ 0,∥ĝ (l) spsp -ĝ(l) ∥ ∞ ≤ 2 -λ max 2 π 2 r 1 (mp + 1) mp . . . (mp -r + 2) d r dλ r ĝ(l) [0,λmax] ∞ Therefore, ∥ĝ (l) spsp -ĝ(l) ∥ ∞ is of O (K r (ĝ, λ cut )(mp) -r) and we can find a scaling constant that depends only on r. Since the Lipschitz constant of σ is 1, we find for any graph G with eigenvalues λ and any graph signal0 ̸ = X ∈ R |G| , (S 2 GNN (l) m ) G -(IGNN (l) ) G (X) F ∥X∥ F ≤ V ĝ(l) spsp -ĝ(l) (λ) ⊙ V ⊤ XW (l) F ∥X∥ F ≤ ∥ĝ (l) spsp -ĝ(l) ∥ ∞ (V V ⊤ )XW (l) F ∥X∥ F ≤ ∥ĝ (l) spsp -ĝ(l) ∥ ∞ ∥W (l) ∥ 2 = O [∥W (l) ∥ 2 K r (ĝ, λ cut )](mp) -rwith a scaling constant that depends only on r. Exactly the same procedure and bounds hold if we instead keep m fixed and increase p. This finishes the proof of Lemma 1.We can now prove the main theorem by induction. Lemma 1 gives the initial step. Now, assume the theorem holds for L IGNN layers. We can then choose S 2 GNN m m∈N =S 2 GNN (L+1) m • S 2 GNN (L•••••1) m m∈N, where S 2 GNN (L+1) m are the approximating models fulfilling Lemma 1, while S 2 GNN(L•••••1) mfulfill the induction assumption. We assume fixed p and increasing m, but the proof is fully analogous in the other case. Applying the same decomposition to (IGNN m ) m∈N lets us express the error on a graph sequence (G m ) m∈N as(S 2 GNN m ) Gm -(IGNN) Gm (X) F ∥X∥ F = (S 2 GNN (L+1) m • S 2 GNN (L•••••1) m ) Gm -(IGNN (L+1) m • IGNN (L•••••1) m ) Gm (X) F ∥X∥ F ≤ (∥X∥ F ) -1 (S 2 GNN (L+1) m • S 2 GNN (L•••••1) m ) Gm -(S 2 GNN (L+1) m • IGNN (L•••••1) m ) Gm (X) F + (∥X∥ F ) -1 (S 2 GNN (L+1) m • IGNN (L•••••1) m ) Gm -(IGNN (L+1) m • IGNN (L•••••1) m ) Gm (X) F ≤ ∥ĝ∥ max ∞ ∥W ∥ max 2 + O(K max r (λ cut ) (pm) -r ) O C L (∥ĝ∥ max ∞ , ∥W ∥ max 2 ) K max r (λ cut ) (pm) -r + O(K max r (λ cut ) (pm) -r )(∥ĝ∥ max ∞ ∥W ∥ max 2 ) L = O C L+1 (∥ĝ∥ max ∞ , ∥W ∥ max 2 ) K max r (λ cut ) (pm) -r .We first used the triangle inequality in line 3 to split the difference into two terms. Next, we bound the first term using the induction assumption, as well as the Lipschitz constant of S 2 GNN (L+1) m, which in turn, by Lemma 1, is the Lipschitz constant of IGNN (L+1) m up to a term of O(K max r (λ cut ) (pm) -r ). We moreover bound the second term using the Lipschitz constant ofIGNN (L•••••1) m, as well as Lemma 1 to arrive at the final result.We next prove the stability of our positional encodings: Theorem 5. The Positional Encodings PE in Eq. 5 are stable according to Definition 1.Recall the definition of stability via Hölder continuity: Definition 1 (Stable PE). (Huang et al., 2024) A PE method PE : R n×k × R k → R n×k is called stable, if there exist constants c, C > 0, such that for any Laplacian L, L ′ , and P * = arg min P ∥L -P LEQUATIONFor this proof, we build on the work of Huang et al. (2024) where the authors show that under the assumptions of Definition 2, and some minor adjustments, a positional encoding of the following form Eq. 23 is stable (Theorem 10).SPE(V , λ) = ρ V diag (ϕ 1 (λ)) V ⊤ , V diag (ϕ 2 (λ)) V ⊤ , . . . , V diag (ϕ k (λ)) V ⊤ (23)Definition 2. The key assumptions for SPE are as follows:• ϕ ℓ and ρ are permutation equivariant.• ϕ ℓ is K ℓ -Lipschitz continuous: for any λ, λ ′ ∈ R k , ∥ϕ ℓ (λ) -ϕ ℓ (λ ′ )∥ F ≤ K ℓ ∥λ -λ ′ ∥. • ρ is J-Lipschitz continuous: for any [B 1 , B 2 , . . . , B k ] ∈ R n×n×k and [B ′ 1 , B ′ 2 , . . . , B ′ k ] ∈ R n×n×k , ∥ρ (B 1 , B 2 , . . . , B k ) -ρ (B ′ 1 , B ′ 2 , . . . , B ′ k )∥ F ≤ J k l=1 ∥B ℓ -B ′ ℓ ∥ F .Theorem 10 (Stability of Eq. 23 by Huang et al. (2024) ). Under Definition 2, SPE (Eq. 23) is stable with respect to the input Laplacian: for Laplacians L, L ′ ,EQUATIONwhere the constants areα 1 = 2J k l=1 K ℓ , α 2 = 4 √ 2J k l=1 M ℓ , and α 3 = J k l=1 K ℓ .Here M ℓ = sup λ∈[0,2] k ∥ϕ ℓ (λ)∥ and again P * = arg min P ∈Π(n) L -P * LP ⊤ * F . The eigengap γ = λ k+1 -λ k is the difference between the (k + 1)-th and k-th smallest eigenvalues, and γ = +∞ if k = n.We prove a similar bound for general weighted adjacency matrices A ∈ R n×n ≥0 (note that such a stability result would be trivial if we restrict A ∈ {0, 1} n×n , since any function on a finite set is Lipschitz continuous). To achieve this, we need a technical assumption in order to ensure that the function values do not blow up and degree normalization is indeed a Lipschitz continuous function: We assume that the domain of A is restricted to (symmetric) matrices whose degrees are uniformly bounded by some constants 0 < Dmin < Dmax :EQUATIONTo decompose the proof into smaller pieces we commonly use the well-known fact that the composition of Lipschitz continuous functionsf 1 •f 2 , with constants C 1 and C 2 , is also Lipschitz continuous ∥f 1 (f 2 (y)) -f 1 (f 2 (x))∥ ≤ C 1 C 2 ∥y -x∥ with constant C 1 C 2 .Proof. Our proposed encoding (Eq. 5) matches roughly Eq. 23. Specifically,ϕ ℓ (λ) = softmax( (λj -λ)⊙(λj -λ) /σ 2 ) with σ ∈ R >0 . However, ρ ℓ (B 1 , B 2 , . . . , B k ) does not directly match || k j=1 [B j ⊙ A] • ⃗ 1, since itis also a function of the adjacency A. Nevertheless, we show that ϕ ℓ is K ℓ -Lipschitz continuous and ρ is J-Lipschitz continuous, where we also bound the change of A.We will start with ϕ ℓ (λ) = softmax( (λj -λ)⊙(λj -λ) /σ 2 ). The softmax is well-known to be of Lipschitz constant 1 w.r.t. the L 2 vector norm/Frobenius norm. -x /σ has a Lipschitz constant of 1/σ. This leaves us with the quadratic term ψ u (λ) = (λ u -λ) ⊙ (λ u -λ) where we bound the norm of the Jacobian J ψu =          -2(λ u -λ 1 ) 0 . . . 0 . . . 0 0 -2(λ u -λ 2 ) . . .EQUATIONthat is zero everywhere except for the diagonal entries, excluding its u-th entry. Thus,∥J ψu ∥ F ≤ 2k max v∈{1,2,...,k} (λ v -λ u ) ≤ 2k(λ k -λ 1 ) ≤ 4k, as 0 = λ 1 ≤ λ k ≤ 2.We can therefore use K ℓ := 4k/σ. Now we continue with ρℓ (A, B 1 , B 2 , . . . , B k ). For f (A, B) = (B⊙A)• ⃗ 1 with a general weighted adjacency A ∈ R n×n , we consider∥(B ⊙ A) • ⃗ 1 -(B ′ ⊙ A ′ ) • ⃗ 1∥ F ≤ (A) ∥ ⃗ 1∥ 2 ∥B ⊙ A -B ′ ⊙ A ′ ∥ F = √ n∥B ⊙ A -B ′ ⊙ A ′ ∥ F = √ n∥B ⊙ A -B ′ ⊙ A + B ′ ⊙ A -B ′ ⊙ A ′ ∥ F ≤ (B) √ n∥B ⊙ A -B ′ ⊙ A∥ F + √ n∥B ′ ⊙ A -B ′ ⊙ A ′ ∥ F = √ n∥(B -B ′ ) ⊙ A∥ F + √ n∥B ′ ⊙ (A -A ′ )∥ F ≤ (C) √ n max u,v A u,v ≤ (D) Dmax ∥B -B ′ ∥ F + max u,v B ′ u,v ≤ (E) 1 √ n ∥A -A ′ ∥ F (F ).(27) (A) holds by Cauchy-Schwarz, (B) by triangle inequality, (C) by Cauchy-Schwarz, (D) follows from the domain of A, and (E) is true since the largest eigenvalue ofB = V ϕ ℓ (λ)V ⊤ is 1 because ϕ ℓ (λ) j ≤ 1, ∀1 ≤ j ≤ k.To further bound (F), i.e. ∥A -A ′ ∥ F , note that ∥L -L ′ ∥ F = ∥D -1/2 AD -1/2 - D ′-1/2 A ′ D ′-1/2 ∥ F . For g(A) := D 1/2 AD 1/2 ,our initial assumption from Eq. 25 yields the existence of a Lipschitz constant C Dmin, Dmax for g, which can be verified by computing the partial derivatives of g. Thus, we can bound∥A -A ′ ∥ F = ∥g(D -1/2 AD -1/2 ) -g(D ′-1/2 A ′ D ′-1/2 )∥ F ≤ C Dmin Dmax , Dmax Dmin ∥D -1/2 AD -1/2 -D ′-1/2 A ′ D ′-1/2 ∥ F = C Dmin Dmax , Dmax Dmin ∥L -L ′ ∥ F =: α 4 ∥L -L ′ ∥ F . (28)As concatenation of k vectors || k j=1 x has a Lipschitz constant of 1, we have J = √ n Dmax . Moreover, we have an additional term for the RHS of Eq. 24 with constant α 4 √ nk, coming from (F) and Eq. 28.To finalize the proof, we restate the beginning of the proof of Huang et al. (2024) and incorporate the additionalEQUATION)Including the extra term stemming from our A-dependent ρℓ (A, B 1 , B 2 , . . . , B k ), the stability guarantee reads∥SPE(EVD(L), L) -P * SPE (EVD (L ′ ) , L)∥ F ≤ (α 1 + α 2 ) k 5/4 ∥L -P * LP ⊤ * ∥ F + α 2 k γ + α 3 + α 4 √ nk L -P * LP ⊤ * F(30) with the newly introduced α 4 arising as Lipschitz constant of (inverse) degree normalization. The proof is complete.Windowing for "eigengap" independent bounds. Note that C depends on the eigengap between 1 /λ k+1 -λ k at the frequency cutoff. One should be able to improve upon this bound with windowing (see Fig. 7 )), effectively lowering the Lipschitz constant of ĥj (λ) around λ k . We leave a formal treatment of this insight to future work.We next prove the expressivity of a GNN/S 2 GNN in combination with our positional encodings: Theorem 6. S 2 GNNs are strictly more expressive than 1-WL with the PE of Eq. 5.For this, we assume that the positional encodings are the only node attributes, subsuming a constant feature or that there is a linear transformation on the raw features. We require that the choice of spatial MPGNN / spectral filter is at least as expressive as the 1-WL test, which is the case, e.g., for GIN. Moreover, we assume that the node-level embeddings are aggregated to the graph level using summation.Proof. To show that GNN(PE(V , λ)) is strictly more expressive as 1-WL. For all graphs that 1-WL can distinguish, the GNN may learn to ignore the PE. Thus, we only need to prove that the positional encodings/node features of PE(V , λ) suffice to distinguish some graphs that 1-WL could Figure 15 : Our positional encodings PE Eq. 5 illustrated in the node colors and sizes. We plot all 5 (rows) 3-regular graphs with 8 nodes and all possible dimensions of the encoding (columns). We use σ = 0.001. Color denotes the sign, and size encodes the absolute value. We hypothesize that the visual "smoothness" between graphs and dimensions is due to our PE's stability (Theorem 5). not distinguish. For all graphs that 1-WL can distinguish we know, by assumption, that the GNN can distinguish the graphs.As Li et al. (2020) point out, 1-WL (and MPGNN that are as capable as 1-WL) cannot distinguish degree-regular graphs with the same number of nodes and degrees. A degree regular graph is a graph where each node has the same degree. This is closely related to Theorem 11.We next show that our PE alone distinguishes certain degree-regular graphs. In this construction, we consider all 3-regular graphs with n = 8 nodes for this (see Fig. 15 ). The encodings ⃗ 1 PE result in the following values with σ = 0.001 and rounded to max 2 decimal places:⃗ 1 ⊤ PE(EVD(L 1 )) = [3 1.73 1 0.41 -1 -1 -1.73 -2.41] ⃗ 1 ⊤ PE(EVD(L 2 )) = [3 1.56 0.62 0.62 0 -1.62 -1.62 -2.41] ⃗ 1 ⊤ PE(EVD(L 3 )) = [3 1.73 1 0.41 -1 -1 -1.73 -2.41] ⃗ 1 ⊤ PE(EVD(L 4 )) = [3 1 1 0.41 0.41 -1 -2.41 -2.41] ⃗ 1 ⊤ PE(EVD(L 5 )) = [3 1 1 1 -1 -1 -1 -3] (31)By constructing examples, this shows that our PE can distinguish 4 out of the 5 3-regular graphs with 8 nodes. Thus, our PE may distinguish at least some graphs that 1-WL cannot. This concludes the proof.While it is well-known that common spatial MPGNNs are at most as expressive as 1-WL and that spectrally designed GNNs can be more expressive than 1-WL (Theorem 2 of Balcilar et al. (2021a) ), we show that spectral GNNs are not able to distinguish degree-regular graphs. This upper bound was not known/formalized prior to our work (Bo et al., 2023b) . Fortunately, our PE largely mitigates the limitation. The improved expressivity of our positional encodings, along with their efficiency, stems from the element-wise product with A (see also Geerts (2021) ).Theorem 11. Spectral filters V diag(ĝ(λ))V ⊤ ⃗ 1 are strictly less expressive than 3-WL with LaplacianL = D -A, L = I -D -1 A, or L = I -D -1/2 AD -1/2 .Proof. The proof relies on properties of the eigenvectors for the different choicesL u = D -A, L rw = I -D -1 A, or L s = I -D -1/2 AD -1/2 . For L u ⃗ 1 = λ 0 ⃗ 1 = 0 and L rw ⃗ 1 = λ 0 ⃗ 1 = 0 the first eigenvector is constant.The first eigenvector of L s is D 1/2 ⃗ 1 (ignoring normalization). Thus, for degree-regular graphs, the first eigenvector of L s is also constant.By the orthogonality of eigenvectors, v u ⊥ v v if u ̸ = v, we know that all other eigenvectors are orthogonal to constant node features. Consequently, the "Fourier transformed" node features are V ⊤ ⃗ 1 = √ n 0 . . . 0 for all three choices L u , L rw , and L s . Since this is true for all degreeregular graphs, spectral GNNs cannot distinguish degree-regular graphs with the same number of nodes.Since the 3-WL test can distinguish some degree-regular graphs, 3-WL is strictly more expressive than a spectral GNN.Corollary 1. "Spectrally designed" MPGNNs that use a polynomial parametrization of filter diag(ĝ(λ)) are strictly less expressive than 3-WL with the same choices for L.Proof. With a polynomial parametrization of the spectral filter ĝ(λ), we knowV (ĝ(λ)⊙[V ⊤ x]) = V diag(ĝ(λ))V ⊤ x = ĝ(L)x = p j=0 γ j L j x (see § 2).Due to this equivalence between a spectral and spatial filter and the constant node features x = ⃗ 1, any polynomial filter p j=0 γ j L j ⃗ 1 cannot distinguish degree-regular graphs. This argument also holds if the polynomial filter is normalized by the maximum eigenvalue as done by ChebNet (Defferrard et al., 2017) .We next provide insights, details, and remarks on the details and variants of S 2 GNNs, accompanying the main section § 3. The structure roughly follows the main body.Next to the overview in Fig. 2 and the method description of the main part, we provide pseudocode in Algo. 1 for a Spatio-Spectral Graph Neural Network on a graph with node attributes, and in Algo. 2 for a spectral filter.Algorithm 1 Spatio-Spectral Graph Neural Network (S 2 GNN), implementing Eq. 11: Input: Adjacency A ∈ R n×n ≥0 , , node attributes X ∈ R n×d (0) , number of eigenvectors k 2: V , λ ← EVD(L(A), k) 3: H (0) ← X + PE(V , λ) 4: for l ∈ {1, 2, . . . , ℓ} do 5: H (l) ← Spectral (l) (H (l-1) ; V , λ) + Spatial (l) (H (l-1) ; A) 6: Return H (ℓ)Algorithm 2 Real-valued spectral filter of Eq. 41: Input: Node embeddings H (l-1) ∈ R n×d (l-1) , eigenvalues λ ∈ [0, 2] k , eigenvectors V ∈ R n×k 2: ĝ(l) ϑ (λ) ← Smearing(λ)W ⊙ Window(λ) 3: Ĥ(l-1) ← V ⊤ H (l-1) 4: Ĥ(l) ← s (l) ζ (ĝ (l) ϑ (λ) ⊙ Ĥ(l-1) ) 5: H (l) ← V Ĥ(l) 6: Return H (l) (a) (b) (c) (d) (e)Figure 16 : Sketch of intra-and inter-cluster message passing capabilitiesV (ĝ ϑ (λ) ⊙ [V ⊤ X] = [ k j=1 ĝϑ (λ j )v j v ⊤ j ]X.The "star" node reflects the global Fourier coefficient and colors/widths illustrate its signed and weighted message passing. We show the first four eigenvectors, order nodes left to right in v j , and sum repeated eigenvalues.In Fig. 16 , we provide further examples of hierarchies/eigenspaces spectral filters have access to, complementing Fig. 3 & 4 . Here and in the main part, we use the main diagonal of j s.t. λj =λu v j v ⊤ j for deciding on the edge weights of the graph structures, potentially summing over multiple eigenvectors with identical values λ j = λ u . We take the product j s.t. λj =λu sign(v j ) for visualizing the sign of the n edges for the global aggregation.For all graphs, the first eigenvector denotes the constant signal (for L = D -A). For (a-d), we observe that the second eigenspace roughly describes a half oscillation, i.e., the left vs. right part of the graph. The third eigenspace separates the middle parts. For (a), the fourth eigenspace models the interactions between the extremal nodes. For (b-d), the frequency increments again, effectively clustering the graph in four roughly equal pieces. For (e), the eigenspaces model the interplay between (automorphic) inner and outer structures, as well as the vertical and horizontal symmetry.Composing a residual connection with a graph filterG = diag(ĝ(λ)) ∈ R n×n yields Y = V GV ⊤ H + H = V (G + I)V ⊤ H, chaining multiple filters (without nonlinearities) re- sults in V G 2 V ⊤ V G 1 V ⊤ H = V G 2 G 1 V ⊤ H. Chaining and residual connections resolve to V (G 2 G 1 + G 2 + G 1 + I)V ⊤ H.Hence, an arbitrary sequence of graph filters (Eq. 2) can be more flexible due to the interactions between filters. Note that this composition is only true in the absence of nonlinearities. Nevertheless, the main intuition about how filters interact remains approximately the same also in the light of nonlinearities.A sensible default is to focus on the low frequencies. We specifically identify the following six reasons: (1) Low frequencies model the smoothest global signals w.r.t. the high-level graph structure (see Fig. 3 & 4 ). ( 2) Gama et al. (2020) find that, under a relative perturbation model (perturbation budget proportional to connectivity), stability implies C-integral-Lipschitzness (∃C > 0 : |λ dĝ /dλ| ≤ C), i.e., the filter can vary arbitrarily around zero but must level out towards larger λ. Stability to graph perturbations is a strong domain-agnostic prior. (3) Many physical long-range interactions are power laws with a flattening frequency response. For example, we construct an explicit graph filter modeling the electric potential of charges in a 1D "ion crystal" ( § G) and find that a lowpass window is optimal. (4) Sequence models like Hyena (Poli et al., 2023) apply global low-pass filters through their exponential windows. ( 5) Cai et al. (2023) prove that an MPGNN plus virtual node (see § E) can emulate DeepSets (Zaheer et al., 2017) and, thus, approximate self-attention to any precision. Nonetheless, we find that a virtual node alone does not necessarily yield good generalization ( § 3.1.1 & 4.1). ( 6) Nonlinearities "spill" features between frequency bands (Gama et al., 2020) . This includes spillage from higher frequencies to the band of the spectral filter. Gama et al. (2020) argue that this spillage makes it possible to learn stable yet expressive graph filters and is also a feature of stable message passing models.For scaling a single to graphs of different orders of magnitude, it can be beneficial to rescale the eigenvalues before learning the filter ĝϑ (λ). That is, we use ĝ(l) ϑ ( λ) with rescaled λ. For example, the eigenvalues for a path/sequence are λ j ≈ (1 -cos( πj /n)). Thus, the resolution is poor, especially for the eigenvalues close to zero since cos approaches slope 0. For this reason, we consider rescaling the eigenvalues withEQUATIONor λj = n /π cos -1 (1 -λ j ) (33) The latter is, e.g., convenient for identifying the second lowest eigenvalue regardless of n. Due to the poor numerical properties of these relations, we evaluate cos-1 (1 -λ j ) = tan -1 ( √ 2λj -λ 2 j/1-λj ) instead.While the GFT and its inverse preserve the norm of the input (e.g., ∥ x∥ 2 = ∥V ⊤ x∥ 2 = ∥x∥ 2 ), this is not true if operating on a truncated frequency spectrum or if the filter ĝϑ (λ) suppresses certain frequencies. For example, in the example of a virtual node (for simplicity here with L = D -A), a signal x that is zero at every node but one at a single node, then the signal will be equally scattered to every frequency. Then, suppressing all frequencies but λ = 0, yields ∥V1 {0} V ⊤ x∥ 2 = 1 / √ n.Motivated by this unfortunate scaling, we also consider normalization in the spectral domain. Specifically, we normalizeĤ = ĝϑ (λ)⊙ V ⊤ f θ (H) ∈ R k×d s.t. Ĥj ← (1-a j ) Ĥj +a j Ĥj /∥ Ĥj ∥2with learnable a ∈ [0, 1] d . This allows, e.g., broadcasting a signal from one node without impacting its scale. However, we empirically find that this normalization only helps only marginally in the over-smoothing experiment (Di Giovanni et al., 2023a) and otherwise can destabilize training. We also consider variants where the norm in the spectral domain is scaled with the norm of the signal in the spatial domain with more or less identical results. We hypothesize that such normalization is counter-productive for, e.g., a bandpass filter if the signal does not contain the corresponding frequencies.For the spectral filter of Eq. 3, we use f(l) θ ( Ĥ(l) ) = H (l) ⊙ [σ(H (l) W (l) G,ℜ ) + i • σ(H (l) W (l)G,ℑ )] and subsequently map the result of Spectral back the real domain, e.g., using wEQUATIONℑ ∈ R d and real ℜ(•) as well as imaginary component ℑ(•). For the positional encodings PE(V , λ) of § 3.2.4, we use A s in Eq. 5 and concatenate real as well as imaginary components. The neural network for the spectral domain s ζ of § 3.2.2 generalizes without adjustment. Similar to Koke & Cremers (2024) , one could also employ complex weights; however, we do not.We use readily available eigensolvers (scipy) and, thus, use a fixed number of eigenvectors (typically k ≪ n) instead of determining k based on λ cut . The partial eigendecomposition is of complexity O(km) for m edges, while the spectral filter has complexity O(kdn). On a different remark, we batch multiple graphs using block diagonal matrices (Fig. 17 ). Spectral graph-level readouts. The key insight is that frequencies are a global concept, and hence, the GFT can be used for global readouts in graph-level tasks. With k ≪ n, such a readout is practically free in the presence of intermediate spectral layers and of O(kn) otherwise. Thus, there is the opportunity for a computationally convenient aggregation of global information, including a sort of graph-level "jumping knowledge" (Xu et al., 2018) . The only caveat is that the Fourier coefficients are not unique due to the ambiguity in the eigendecomposition. To maintain permutation equivariance, we take the absolute value and aggregate over dimension k in Eq. 4 instead of the multiplication with V .                                                                    A 1 A 2 . . . A b 0 0 i n i i n i                                                                     V 1 V 2 . . . V b 0 0 i k i i n iWe observe that such intermediate readout can improve performance slightly, e.g., on TPUGraphs. However, we leave a systematic evaluation of its benefits for future work.Linear bottle necks. To circumvent overfitting, we commonly replace the linear transformationsW X in f (l) θ ( Ĥ(l) ) and ĝϑ (λ) with low-rank bottlenecks W 2 W 1 X, s.t. W ∈ R d×d , W 2 ∈ R d×d ′ , W 1 ∈ R d ′ ×d , and d ′ < d.We expect that many common graph benchmarks do not have or only insignificant long-range interactions. We observe that MPGNNs are less likely to overfit, perhaps since locality is a good inductive bias in many circumstances (Bronstein et al., 2021) . Moreover, we observe that the spectral filter ( § 3.2.1) may converge slowly and get stuck in local optima. We find that a sufficient amount of randomly initialized filters mitigates this issue to a large extent. Further, one can introduce inductive biases via windowing functions (Fig. 7 ), like the exponential window used by Hyena (Poli et al., 2023) .Even if the true causal model generating the target consists of long-range interactions, it might be sufficient to model the training data solely using (potentially spurious) local interactions. This might be especially true if the training nodes are samples from a "small" vicinity of the graph (e.g., OGB Products (Hu et al., 2020) ).Closely related to the previous point is the amount of available training data. We hypothesize that S 2 GNNs are more data-hungry than their purely spatial counterpart. That is, to reliably detect (nonspurious) long-range interactions in the training data, a sufficient amount of data is required. Similar findings have been made, e.g., in the image domain (Dosovitskiy et al., 2021) .Except for heterophilic graphs, direction plays a small role in graph machine learning even though many benchmark tasks actually consist of directed graphs (Rossi et al., 2023) . Moreover, there is a lack of benchmarks involving directed graphs, which require long-range interactions. Note that most of the theoretical findings generalize to directed graphs under appropriate modeling decisions/assumptions. However, we do not make this discussion explicit since MPGNNs for directed graphs are still actively researched.Since a lot of the previous points hover around the insufficiency of the available benchmarks, we propose two new tasks § 4.1 and derive further datasets, e.g., for associative recall.While we demonstrate the practicality of S 2 GNNs in § 4.3 on large-scale benchmarks, the partial eigendecomposition EVD starts to become costly on the largest graphs we use for evaluation. Even though we did not experiment with lowering the requested precision, etc., we expect that for scaling further, naïve approaches might not be sufficient. One direction could be to utilize GPUs instead of CPUs or to adapt concepts, e.g., from spectral clustering (von Luxburg, 2007) .Even though there are many important reasons why we should utilize a spectral filter on the low end of the spectrum, there might be tasks for which this choice is suboptimal. One way to estimate the frequency band to which one should apply a spectral filter is via a polynomial regression and then determine where the derivative is maximal. Note that it is efficient to calculate the eigenvectors around an arbitrary location of the spectrum, e.g., with the "shift-invert mode" of scipy/ARPACK (Lehoucq et al., 1998) .Due to the many possible design decisions of spectrally parametrized filters, the neglect of spectral filters in prior work, and the lack of appropriate benchmarks, it was not possible to ablate all the details. We expect that future work will discuss the specific building blocks in greater detail.We expect that S 2 GNNs will have similar societal implications as other model developments like Convolutional Neural Networks (CNNs) (LeCun et al., 1989) , LSTMs (Hochreiter & Urgen Schmidhuber, 1997) , transformers (Vaswani et al., 2017) , or modern Graph Neural Networks (Gilmer et al., 2017) . Since such models may be used as building blocks in architectures for predictive tasks, generative modeling, etc., they have a wide range of positive and negative implications. Nevertheless, we expect that S 2 GNNs will not have more negative implications than other machine learning model innovations.This section provides further details on the experimental setup ( § M.1), the computational cost ( § M.3), and graph constructions with additional experimental results for the clustering tasks ( § M.6); likewise we provide details for the distance regression ( § M.7), arXiv-year ( § M.8), and provide nodes on the graph construction in TPUGraphs ( § M.10). Note that the sections on clustering ( § M.6) and distance regression ( § M.7) also contain ablations and further insights.Implementation. The code base is derived from Cao et al. (2023) , which on the other hand derive the code of Rampášek et al. (2022) . The implementation heavily relies on PyTorch geometric (Fey & Lenssen, 2019) .Datasets. We collect the main statistics, including licenses, for the datasets in Table 6 . The provided code will download all datasets along with the experiment execution, except for TPUGraphs, where one should follow the official instructions. Due to the high variation in results, we merge all "layout" datasets and present the results on this joint dataset. We use the fixed public splits for all experiments and proceed accordingly for our datasets (see § M.6 and § M.7).Hyperparameters. While we provide full parameters for all experiments and models in our code, we gather an overview of the used S 2 GNNs variants here. The parameters were determined through cascades of random search throughout the development of the method. We list the most important parameters in Table 7 . Usage of external results. The performance of baselines is commonly taken from leaderboards and the respective accompanying papers. This specifically includes the results in Table 1 , Table 2 , and Table 11 .Setup. For clustering ( § M.6), distance regression ( § M.7), and arXiv-year ( § M.8) we report the detailed setup in the respective sections. For the other tasks, the relevant details are:• Peptides: We follow the setup and implementation of Rampášek et al. (2022) . That is, we train for 250 epochs with a batch size of 200. We rerun experiments on 10 random seeds.• Over-squashing: We derive the setup from Di Giovanni et al. (2023a). In the main part (Fig. 5 ), for the GCN, we report the numbers of their Figure 3 for a GCN on "Clique Path" graphs. For the spectral filter, we actually consider the more challenging setting where we do not train one model per graph size. Instead, we train one model for all sequence lengths. The task is to retrieve the correct of five possible classes on the other end of the graph. In the extended experiment of Fig. 12 , we compose the dataset of "Clique Path" and "Ring" graphs (see Di Giovanni et al. (2023a) ). To avoid m = O(n 2 ), we limit the fully connected clique to 15 nodes. For training and validation, we enumerate all graphs with even n ∈ {4, 6, . . . , 50} and train for 500 epochs. For test, we enumerate the graphs with even n ∈ {52, 54, . . . , 100}. We rerun experiments on 10 random seeds.• Associative recall: We construct one dataset consisting of key-value sequences of length 20 to 999. As Poli et al. (2023) , we use a vocabulary of 30. We sample 25,000/500 random graphs for train/validation. For the test set, we randomly generate 500 graphs for the sequence lengths of 1,000 to 1,199. We train for 200 epochs. In the experiment with validation/test sequence length 30k (Table 2 ), we generate 10,000 training graphs of length 29,500 to 30,499 and finetune S 2 GNNsGCN from the smaller setup. We rerun experiments on 10 random seeds. • OGB Products: Even though full-graph training with 3 layers GCN plus one spectral layer fits into a 40 GB A100 GPU, we find that batched training works better. We randomly divide the graph during training into 16 parts and train a 6-layer S 2 GAT with spectral layers after the second and last message passing step. Inference is performed on the entire graph at once. We rerun experiments on 5 random seeds.• TPUGraphs: This is the only dataset where we use a transformer to model ĝ instead of the procedure detailed in § 3.2.1. We fix the number of eigenvectors to k = 100 and do not apply any windowing. Due to the large variation of results, we merge all "layout" tasks into a single dataset. Since the default graph construction is not able to express all relevant information, we adapt it as detailed in § M.10, however, the empirical impact was small. TPUGraphs "layout" consists of a few hundred distinct graph structures with a large variation on the node-level configuration/features. We sample 10,000 configurations for each graph structure of each "layout" sub-split. Here, we introduce two batch dimensions:(1) batching over multiple graphs and (2) batching over the configurations. In each training step of the 1,000 epochs, we sample a small subset of configurations per graph structure and apply a pairwise hinge loss to rank the configurations. We do not perform random reruns due to the computational cost.In Fig. 6 and Fig. 8 , we provide qualitative insights about the approximation of filters and ringing.In Fig. 6 , we construct a true filter by adding a discontinuous filter at λ = 0 and a polynomial filter of order 3. For the spectral part, we use the true filter values and fit a Chebyshev polynomial on the remaining part. We then plot the response of the true filter and its approximations on a path graph with 21 nodes and L = I -D -1 A.Similarly, in Fig. 8 , we use a path graph with 100 nodes and L = I -D -1 A. We then construct a perfect low pass (k = 25) and approximate a rectangular wave.We report the computational cost for the experiments in Table 7 for a single random seed. On top of the pure cost of reproducing our numbers, we conducted hyperparameter searches using random search. Partially, we required 100s of runs to determine good parameter ranges. A generally wellworking approach was first to reproduce the results of the best available MPGNN in prior work. Thereafter, we needed to assess how likely additional capacity would lead to overfitting. Usually, we reduced the number of message-passing steps, added the spectral filter, and determined appropriate values for the number of eigenvectors k. In the light of overfitting, it is a good idea to lower the number of Gaussians in the smearing of the filter parametrization ( § 3.2.1), introduce bottle-neck layers ( § J), and use fewer spectral filters than hidden dimensions.Runtime with precalculated eigenvectors. In Fig. 18 , we contrast the runtime cost of a spectral convolution with spatial messages passing on ogb-arXiv (170k nodes) of Hu et al. (2020) , using an Nvidia GTX 1080Ti. This essentially compares a sparse matrix multiplication (adjacency matrix) with matrix multiplications on dense "tall and skinny" matrices (GFT). we find that one GCN-layer here is as costly as a spectral filter with approx. k = 2, 500 eigenvectors.10Large-scale benchmarks. On the large-scale datasets OGB Products and TPUGraphs, we perform full-graph training (without, e.g., segment training (Cao et al., 2023) ) using 3 DirGCN layers interlayered with spectral filters targeting a pair-wise hinge loss. The spectral GNN uses the Magnetic Laplacian to incorporate direction. The spatial MPGNN closely resembles the model of Rossi et al. (2023) , except that we half the dimension for the forward and backward message passing and concatenate the result. We shrink the dimensions to model the direction at a very low cost. We conclude that S 2 GNNs can be very practical even if applied at scale and can effectively model long-range interactions also on large graphs. Eigendecompositon. We show the computational cost for the eigendecomposition of a random Erdős Rényi graph (every edge has equal likelihood to be drawn). We use scipy (CPU) and PyTorch (GPU) with default arguments. For the sparse decomposition with PyTorch, we use the svd_lowrank method. Note that the default parameters for PyTorch are usually leading to large numerical errors. Fig. 19 demonstrates that the cost of the eigendecomposition is manageable. For large graphs like ogbn-products (2.5 mio. nodes), the EVD takes around 30 minutes with k = 100 on 6 CPU cores of an AMD EPYC 7542. Note that the default parameters of the eigensolver allow for 1000s of iterations or until the error in the 32-bit float representation achieves machine precision.In the main body we present two ways to combine a spatial and spectral filter: An additive combination (Eq. 1) and an arbitrary sequence of filters (Eq. 2). In this section, we perform an ablation analysis on the peptides-func benchmark and report the results in Table 8 .Instead of summation of the spatial and spectral parts, concatenation is another possible option. Getting input features H (l-1) ∈ R n×d (l-1) , we design Spectral (l) and Spatial (l) to map to R n×d (l) /2 and update the embeddings asEQUATIONAdditionally, we consider normalization of the addends at the end of each embedding update, dividing by 1/ √ 2 (concatenation w/ residual) and 1/ √ 3 (summation w/ residual) as an attempt to keep the variance constant.Inspired by recent advancements in state space models like Mamba (Gu & Dao, 2023) , we also consider modeling an update step in a similar way, identifying the convolutional part with Spatial (l) and the SSM part with Spectral (l) .The following table shows results for the different design choices to combine the spatial and spectral parts, with all hyperparameters being precisely the ones reported in Table 7 for peptides-func. Nevertheless, we find that the spectral filter is well aligned with the cluster structure in these tasks.We plot this some exemplary filter in Fig. 23 . The findings match the explanations of § 4.1 also for CLUSTER. In the remainder of the section, we provide full details of the experiment setups. Moreover, we provide additional results not presented in the main text, including ablations.Setup. To sample an input graph, we start by generating C = 6 p-dimensional cluster centers µ c ∼ U [0, 10] p for c ∈ {0, . . . , C -1} (we use p = 2). Next, we draw n c ∈ {100, . . . , 199} points x ic ∼ N (µ c , 4I p ) which will represent the nodes of the graph. Subsequently, we update the class memberships such that every point is in its most likely class according to the underlying probabilistic model. Finally, we connect each node v to its e v ∼ U ({1, . . . , 10}) closest neighbors by Euclidean distance ∥ • ∥ 2 . This whole procedure is repeated until the generated graph is connected. We then discard the location information and only keep the graph structure. In this way, we generate graphs of an average diameter of ≈ 33. See being in class c -1 and a feature value of 0 means that the class is unknown and has to be inferred by the model. Only one node v c per class is randomly chosen to be labeled and all remaining node features are set to 0. The output labels are defined as the class labels. We use weighted cross entropy loss for training and class-size-weighted accuracy as a target metric. We generate 10,000 training and 1,000 val/test graphs each and report the average ± standard deviation over 3 random reruns.As an underlying spatial model baseline, we use a vanilla GCN (Kipf & Welling, 2017) .We compare this to S 2 GCN, only applying one spectral convolution immediately before the last spatial layer. We investigate the influence of the number k ∈ {0, 1, . . . , 10} of eigenvectors to be taken into account with 4 spatial layers, with k = 0 indicating the absence of a spectral layer (see Fig. 10a , and Table 9 for the underlying data). We also vary the number of spatial MP layers from 2 to 10 and compare the performance of a purely spatial GCN to the corresponding S 2 GCN with one spectral convolution (see Fig. 10b , and Table 10 for the underlying data).Throughout all evaluations, we maintain a consistent hyperparameter configuration: Specifically, we use an inner dimension of 128, GELU (Hendrycks & Gimpel, 2016) as an activation function, no dropout, and residual connections for all spatial and spectral layers. For the spectral layer, we implement the gating mechanism f (l)θ , but abstain from a neural network in the spectral domain ( § 3.2.2), bottlenecks, or parameter sharing. We train for 50 epochs with a batch size of 50, using the AdamW optimizer (Loshchilov & Hutter, 2019) with a base learning rate of 0.003, a weight decay of 0.0001, a cosine scheduler and 5 warmup epochs.Further discussion. The clustering task comes naturally to S 2 GCN, as a spectral layer can simulate certain variations of spectral clustering (von Luxburg, 2007) : Suppose H (l-1) ∈ R n×C is a one-hot encoding of the cluster labels, i.e. H ≡ id, the spectral layer Spectral (l) from Eq. 3 turns into H(l) = V V ⊤ H (l-1) . Hence, H (l) v,c = V ⊤ v,:V vc,: encodes a notion of similarity between a node v and each labeled node v c . This relates to the Euclidean distance∥V v,: -V vc,: ∥ 2 = ∥V v,: ∥ 2 2 + ∥V vc,: ∥ 2 2 -2V ⊤ v,: V vc,:which is more typically used for spectral clustering.M.6.2 SBM Clustering CLUSTER (Dwivedi et al., 2023) Setup. We conduct an ablation study on the original CLUSTER task (Dwivedi et al., 2023) , which uses a similar setup to our GMM clustering task, however drawing from a SBM instead: For each cluster, n c ∈ {5, . . . , 35} nodes are sampled. Nodes in the same community are connected with a probability of p = 0.55, while nodes in different communities are connected with a probability of q = 0.25. While there is no need for long-range interactions in this task, considering that the average diameter of the graphs is just ≈ 2.17, separating the clusters is much harder than in the GMM clustering task (see Fig. 23 for example adjacency matrices from the SBM and GMM models). We use weighted cross entropy loss for training and class-size-weighted accuracy as a target metric. We report the average ± standard deviation over 3 random reruns.Models. In our ablation study, we consider GCN (Kipf & Welling, 2017) , GAT (Veličković et al., 2018) , and GatedGCN (Bresson & Laurent, 2018) as MPGNN baselines, following a setup similar to Dwivedi et al. (2023) . We consider models with 4 layers (roughly 100k parameters) and 16 layers (roughly 500k parameters), while keeping most hyperparameters the same as in the benchmark, including inner dimension, dropout, and the number of heads for GAT. However, our reported baseline results and parameter counts differ slightly as we are using a different post-MP head, where we maintain a constant dimension until the last layer, in contrast to Dwivedi et al. ( 2023) who progressively shrink the inner dimension. We construct the corresponding S 2 GNNs by modifying each baseline model, replacing the 3 rd and the 5 th /15 th layers with spectral layers, ensuring a roughly equivalent parameter count. Additionally, each model is optionally supplemented by our positional encodings PE ( § 3.2.4).We further conduct a hyperparameter search on the most promising base MPGNN candidate, Gat-edGCN, which leads to an optimized version of S 2 GNN. This optimized model has 18 spatial MPGNN layers, spectral layers between all spatial layers, and additional RWSE encodings. The inner dimension is adjusted to keep the model well below a parameter budget of 500k. Finally, we also evaluate S 2 GCN and S 2 GAT using these hyperparameter settings. Throughout all evaluations, we use GELU (Hendrycks & Gimpel, 2016) as an activation function, residual connections for all spatial and spectral layers, and implement the gating mechanism f (l)θ without employing a neural network in the spectral domain. We use a batch size of 128 for training the 4-layer models and 64 for all other models. For the spectral layer, we use the partial eigendecomposition corresponding to the lowest k = 50 eigenvalues (k = 100 for the optimized S 2 GNN versions), spectral normalization, and λ cut = 1.3. For the optimized models, we employ parameter sharing with 128 heads, and a bottleneck of 0.25 in feature gating. We use 8 attention heads for all GAT versions in accordance with Dwivedi et al. ( 2023) (in- For the purely spatial models, we use p = 0.0 as dropout (similar to Dwivedi et al. ( 2023)). We observe this to lead to overfitting for models with spectral layers, for which we set p ∈ {0.1, 0.2}. Hyperparameters differing between the compared models are listed in Table 12 . We train for 100 epochs using the AdamW optimizer (Loshchilov & Hutter, 2019) with a base learning rate of 0.001, no weight decay, and a cosine scheduler with 5 warmup epochs.Results. Results for the CLUSTER task are presented in Table 12 , Table 11 and Fig. 24 . Introducing a spectral layer significantly enhances performance on the 4-layer architectures, both with and without positional encodings. The effect is most pronounced on GCN, where replacing just a single GCN layer by a spectral layer boosts accuracy from 0.504 to 0.655. Notably, introducing two spectral layers still has a consistent positive effect on all 16-layer architectures.Setup. We generate directed random trees with one source by sampling trees with n ∈ {500, . . . , 999} nodes, picking one node at random to declare as a source and introducing edge directions accordingly. To construct random DAGs with long distances, we start from such directed random trees and proceed by adding ⌊n/10⌋ edges at random, choosing each edge direction such that the resulting graph is still a DAG. Additionally, we mark the source node with a node feature. Besides evaluating all models in an in-distribution regime, we also assess the generalization power of the methods by drawing out-of-distribution val/test splits from slightly larger graphs of n ∈ {1000, . . . , 1099} and n ∈ {1100, . . . , 1199} nodes each. We use L 2 loss for training and R 2 as a target metric. We sample 50,000 training and 2,500 val/test graphs each and report the average ± standard deviation over 3 random reruns. Models. As a MPGNN baseline, we use a five-layer directed version of GCN, DirGCN (Rossi et al., 2023) , with three post-message-passing layers, and concatenating instead of averaging over the source-to-target and target-to-source parts. We compare these baselines to S 2 DirGCN of the form Eq. 2 with four spectral layers, alternating spatial and spectral convolutions and employing residual connections. We benchmark versions of S 2 DirGCN that ignore edge direction in the spectral convolution against directed versions in which we set q = 0.001. In all cases, we use the partial eigendecomposition corresponding to the k = 50 lowest eigenvalues. All models are optionally enriched by the positional encodings from § 3.2.4. Throughout all evaluations, we use an inner dimension of 236, GELU (Hendrycks & Gimpel, 2016) as an activation function, and dropout p = 0.05.For the spectral layers, we utilize the gating mechanism f (l)θ , not employing a neural network in the spectral domain, we use spectral normalization, λ cut = 0.1, and a bottleneck of 0.03 in the spectral layer. We train for 50 epochs, using a batch size of 36 and the AdamW optimizer (Loshchilov & Hutter, 2019) with a base learning rate of 0.001, a weight decay of 0.008, and a cosine scheduler with 5 warmup epochs.Results. In Table 13 , we show the performance of the different models on DAGs and trees. We observe that the simple MPGNNs are notably surpassed by all versions of S 2 DirGCN. While S 2 DirGCN achieves nearly perfect predictions on the tree tasks in both the directed and undirected case, the undirected version is outperformed by the directed version on the DAG tasks. Here, performance also reduces slightly in the out-of-distribution regime. The great performance on the tree M.9 Large-Scale PCQM4Mv2 (Hu et al., 2021) We show that S 2 GNNs are very parameter efficient. Even though we only conduct a very rudimentary hyperparameter search, S 2 GNNs keep up with state of the art approaches. Specifically, we adapt the hyperparameters from peptides-func and achieve comparable performance to the state of the art (excluding external data) with about 3-20% of the number of parameters. The "XLA" collections of TPUGraphs contain many constructs that are most certainly suboptimal for a machine-learning-based runtime prediction. However, in our preliminary experiments, we could not show that our graph construction yielded better results in a statistically significant manner. Nevertheless, we include this discussion since it might be insightful.To understand the challenges with the default graph construction, note that in the TPUGraphs dataset each node represents an operation in the compuational graph of Accelerated Linear Algebra (XLA). Its incoming edges are the respective operands, and the outgoing edges signal where the operation's result is used. Thus, the graph describes how the tensors are being transformed. An (perhaps unnecessary) challenge for machine learning models arises from using tuple, which represents a sequence of tensors of arbitrary shapes. In this case, the model needs to reason how the tuple is constructed, converted, and unpacked again. Moreover, directly adjacent tensors/operations can be very far away in the graphs of TPUGraphs.We identified and manually "fixed" three cases to eliminate this problem largely in the TPUGraphs dataset: Tuple-GetTupleElement, While, and Conditional. Since we could not access the configurations in the HLO protobuf files and C++ XLA extraction code, we decided to perform these optimizations ourselves. However, it might be a better strategy to utilize the existing XLA compiler etc.Additionally, to the subsequently described graph structure changes, we extract the order of operands from the HLO protobud files. Outgoing edges are assumed to be unordered except for the Get-TupleElement operation, where the tuple index is used as order. Moreover, we extracted all features masked in the C++ code and then excluded constant features.The dataset contains aggregations via the XLA Tuple operation that are often directly followed by a GetTupleElement operation. To a large extent, these constructs are required for the subsequently discussed While and Conditional operations. Importantly, the model could not extract the relationships through a tuple aggregation since the tuple_index was not included in the default features. Moreover, the resulting tuple hub nodes severely impact the Fourier basis of the graphs (see § 2). We illustrate the graph simplification in Fig. 27 and denote the edge order of incoming edges from top to bottom. The edge order represents the order of operands.We propose dropping immediate Tuple-GetTupleElement constructs and directly connecting predecessors and successors. For this, we generate a graph solely consisting of direct connections and then resolve multiple consecutive Tuple-GetTupleElement constructs via a graph traversal (depthfirst search).Figure 27 : Tuple-GetTupleElement simplification: the Tuple aggregates the output of multiple predecessors/operations and then the GetTupleElement extracts the tensor according to its index (number in respective nodes). We propose dropping immediate Tuple-GetTupleElement constructs and connecting predecessors and successors.We perform the Tuple-GetTupleElement simplification after dealing with While and Conditionals. However, for the sake of simplicity, we will avoid using tuples in the subsequent explanations for While and Conditional. In other words, the subsequent explanations extend to functions with multiple arguments via the use of tuples.The While operation has the signature While(condition, body, init) where condition is a function given the init or output of the body function. Note that in the TPUGraph construction, body as well as condition only represent the outputs of the respective function and their operands need to be extracted from HLO.To avoid hub nodes and to retain the dataflow between operations (important for decisions about the layout), operands and outputs are connected directly. Technically, we am modeling a do-while construct because the condition is not connected to the inputs. Since the successors of the while are of type GetTupleElement, they relabeled to a new node type, signaling the end of a while loop.To support nested while loops, each node in the body is assigned a new node feature signaling the number of while body statements it is part of.Figure 28 : Instead of aggregating everything into a hub node, we propose to connect respective inputs and outputs.Conditional(branch_index, branch_computations, branch_operands) is the most common signature of the Conditional operation, where the integer-valued branch_index selects which branch_computations is executed with the respective input in branch_operands. Similarly to the While operation, we introduce new node types for the inputs of computations and the successors (they are GetTupleElement operations).Figure 29 : Instead of aggregating everything into a hub node, we propose to connect respective inputs and outputs. Here as an example with two conditional computations.Guidelines: The proposed framework S 2 GNN does not pose an unusual risk on top of the common risks for research on machine learning architectures. Our evaluation does not include human subjects, participants, or data. We discuss the broader impact in § L. Thus, in summary, we conform to the ethics guidelines in every aspect.• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes] Justification: We discuss the broader impact in § L. It should be noted that the proposed framework S 2 GNN does not pose an unusual risk on top of the common risks for research on machine learning architectures. We discuss the broader impact in § L.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The proposed framework S 2 GNN does not pose an unusual risk on top of the common risks for research on machine learning architectures. Our evaluation does not include human subjects, participants, or data.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.• Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: We cite all origins of the datasets, even if we derive the benchmark from someone else. We provide an overview of the datasets, including the licenses in Table 6 . We provide code, including configuration, at https://www.cs.cit.tum.de/daml/s2gnn. It also contains an (anonymized) license for its usage and discloses the provenance of the project setup.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 13. New Assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [Yes] Justification: The code and configuration for the dataset on GMM clustering, distance regression, and associative recall will be provided. The generation of all resources takes around one hour and uses solely CPUs. We provide code, including configuration, at https://www.cs.cit.tum.de/daml/s2gnn. The datasets are described in § M, § M.6, and § M.7. Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?38th Conference on Neural Information Processing Systems (NeurIPS 2024).† equal contribution. 
 
 Appendix:  


Decision: This paper studies an essential and timely topic. It proposes a novel hybrid method of spatio-spectral Graph Neural Networks to combine spatial and spectral graph filters and ensure better efficiency and bounds. After a fruitful discussion between the authors and the reviewers, this paper has an overall positive feeling.

The authors should carefully read and address all remaining concerns in the final version of the paper.

Reviews: review: paper_summary: paper_summary: The proposed method Spatio-Spectral Graph Neural Networks (S2GNNs) combine the spectral GCN and spatial GCNs embeddings linearly. The paper considers a deep dive into spectral filter properties and attempt to motivate the combination with spatial filtering. There are well known results repeated in the method section, which could have been part of background. Also, the contribution towards directed graphs is not clear.
strength_weakness: strength_weakness: 1.	There are well known results repeated in the method section. For example, it is established that the spectral filters are permutation equivariant, Locality relates to spectral smoothness.
2.	The paper sometimes fails to connect the dots. For example, MPNNs (spatial) suffer from oversquashing, but at the same time spectral GNNs can handle oversquashing issue. So why do we need to combine them both? Its not clear in Introduction.
3.	Contribution towards the directed graphs is unclear, any spectral filters can be adapted to directed graphs via magnetic Laplacian. The main point of the paper is to demonstrate the superiority against existing methods for undirected settings first.
questions: questions: GNNs aim to learn data driven filters. “Ideal” filters may NOT be a “good” choice always. Then what is the need of implementing an ideal discontinuous spectral filter discussed in Section 3.2?

In the experiments, did you take the directionality into account in S2GNN? If, yes what is the performance of S2GNN without directionality?
limitations: limitations: n/a
paper_summary: paper_summary: The paper presents Spatio-Spectral Graph Neural Networks (S$^2$GNNs), a novel paradigm that combines spatial and spectral parameters to overcome the limitations of $\ell$-step MPGNNs, notably their restricted receptive fields and over-squashing issues. S$^2$GNNs achieve global information propagation efficiently, surpassing MPGNNs in approximation accuracy, expressivity, and scalability, demonstrated through superior performance on various graph tasks and ability to handle large-scale datasets.
strength_weakness: strength_weakness: 1. **Deficiency in Ablation Studies for Embedding Updates**: Eq.~(1) within the paper that updates to the latent embeddings are achieved through the summation of spectral and spatial layers. This method warrants a broader range of ablation studies, such as exploring concatenation, normalization, or adopting a randomized approach to select between spectral or spatial optimizations like [1], to validate the robustness and efficiency of these embeddings.

2. **Inadequate Comparative Analysis for Over-Squashing and Long-Range Interactions**: The paper's treatment of over-squashing and long-range interactions could be enhanced by incorporating comparative analyses with established methods. Notably, the GCNII[2], which utilizes skip connections and PageRank to facilitate up to 64 layers in a GNN, and methods based on random walks that intensify depth-first search capabilities for extended influence, should be examined. This inclusion would enhance the analysis and provide a more rigorous evaluation of the proposed model’s capabilities in handling long-range interactions.

[1] Chang, Heng, et al. "Not all low-pass filters are robust in graph convolutional networks." Advances in Neural Information Processing Systems 34 (2021): 25058-25071.
[2] Chen, Ming, et al. "Simple and deep graph convolutional networks." International conference on machine learning. PMLR, 2020.
questions: questions: See above
limitations: limitations: NA
paper_summary: paper_summary: This paper proposes Spatio-Spectral Graph Neural Networks (S²GNNs) to address the limitations of ℓ-step Spatial Message Passing Graph Neural Networks (MPGNNs), such as limited receptive fields and over-squashing. S²GNNs combine spatial and spectral graph filters for efficient global information propagation, offering tighter error bounds and greater expressiveness than MPGNNs. They outperform traditional MPGNNs, graph transformers, and graph rewirings on benchmarks and scale efficiently to millions of nodes.
strength_weakness: strength_weakness: 1. This paper's writing and structure can be further improved. Adding descriptions to the notation and including pseudocode would enhance readability.
2. The experiment includes too few baselines. For instance, only GAT and GCN are used as baselines in Tables 3 and 4, while some of the latest methods, such as SpecFormer, which also uses Eigenvector, are not considered.
questions: questions: 1. Please refer to the weaknesses mentioned above.
2. What is the upper bound of the graph size that the EVD used in S²GNN can support? For example, can it handle graphs with billions of nodes, such as papers100m?
3. In real tasks, are the eigenvectors $k$ sensitive to performance? For instance, on products and TPU graphs, would using a larger $k$ improve performance?
limitations: limitations: The authors discussed limitations in the Appendix.
paper_summary: paper_summary: The authors are proposing Spatio-Spectral Graph Neural Networks ($S^2GNNs$), a hybrid model that combines (Spatial) Message Passing Neural Networks (MPNNs) and Spectral Graph Neural Networks. The authors argue that, by combining message-passing with spectral filters, the model can better model both local and global information, alleviating phenomena such as over-squashing and under-reaching. The authors discuss the design space of their proposed method and showcase a technique for obtaining position encodings that can make their model more expressive than the 1-WL test. The authors perform multiple experiments on both synthetic and real-world datasets.
strength_weakness: strength_weakness: - The paper's main weaknesses are that the main text is somewhat unfocused, and some related work is not mentioned. I will detail in the following:
	- While the authors generally do a very good job of providing valuable theoretical insights regarding their design choices, I would argue that some of them are unnecessary and distract the reader from the main points of the paper. For instance, I believe the discussion about "high-resolution filters for low frequencies" from lines 149-167 to be somewhat hand-wavy and distracting, especially the remark that "as graph problems discretize continuous phenomena, stability to the graph perturbations is a strong domain-agnostic prior". While I don't think that the discussion is bad, I feel that the supplementary materials would be a much better place in which to include it. 
	- I believe that Subsection 3.3 (Parametrizing Spectral Filters) is somewhat confusing, especially for a reader who is not very familiar with spectral methods. I did not initially understand the discussion about the window (Lines 247-253), and I believe it should contain more context. The authors should either expand the section or move it to the supplementary materials.
	- The proposed positional encodings are a nice addition, but similar laplacian positional encodings, such as [[1]], have not been discussed as related work. Other works, such as [[2]], are cited in the paper but not in the context of the laplacian positional encodings that they propose.
	- Some other recent methods deal with long-range relationships on graphs and have not been discussed in the related work. Notably, GRED [[3]] uses Linear Recurrent Units to model distant relationships with a mechanism similar to k-hop MPNNs and also draws inspiration from SSMs; AMP [[4]] has a message filtering mechanism and is adaptively increasing the number of layers depending on the dataset; NBA [[5]] is passing messages in a non-backtracking manner, obtaining a slightly better score than $S^2 GCN$ on *peptides-struct* ($0.2424$ v.s. $0.2447$), and PR-MPNN [[6]] rewires graphs adaptively by learning a probability distribution for the adjacency matrices.
- The authors mention on line 69 that the runtime of their method is equivalent to that of an MPNN, but I could not find a table containing a direct comparison. It would be great if this were included somewhere. There are training times on page 35, Table 6, but no comparison with a simple MPNN (such as a GCN or a GIN model).


Overall, I believe that the paper is very strong. The model is conceptually simple but contains novel ideas. The empirical results are very good, and the authors test their method on many different scenarios. I found the mechanistic in-context learning experiments to be particularly interesting. I believe the work could have a considerable impact on how we design GNN architectures going forward.

As stated previously, the main weakness is that the paper feels somewhat unfocused, and some important related work is omitted. This is a somewhat strange complaint to have, but I feel that the scope of the paper is very broad. In my opinion, the overall style of the paper would be better suited for a journal submission if some details from the appendix were moved into the main text. For a conference with more limited space, I would advise the authors to move some of the more in-depth details to the appendix. This is in no way an argument for rejection, I just believe that the paper would be a lot nicer to read if it didn't contain these details in the main text.

I recommend a 7 (Accept), but I will update my score to an 8 (Strong Accept) if the authors expand their related work discussion and provide more details for Subsection 3.3. I would also recommend to the authors that they move some of the details from Section 3 to the Appendix, but not doing so will not affect my final score. 



[1]: https://arxiv.org/abs/2110.07875
[2]: https://arxiv.org/pdf/2012.09699
[3]: https://arxiv.org/pdf/2312.01538
[4]: https://arxiv.org/pdf/2312.16560
[5]: https://arxiv.org/pdf/2310.07430
[6]: https://arxiv.org/abs/2310.02156
questions: questions: - In the graph rewiring literature, the models are often benchmarked on other graph-level datasets (such as QM9, Zinc, or the TUDataset collection). Did the authors try their method on these datasets? I am not requesting these experiments, but they would be a nice addition.
- The Long-Range Graph Benchmark also contains COCO-SP, PascalVOC-SP, and PCQM-Contact. Did the authors try to benchmark their models on these datasets? I'm personally not a fan of COCO-SP and PascalVOC-SP because they're obtained from Computer Vision datasets, so the GNN inductive bias might not be appropriate for them. Still, PCQM-Contact is a molecular dataset, and it would be interesting to see some results on it.
limitations: limitations: - The limitations are addressed in the appendix.

Comments: title: Conclusion of Discussion Period
Official Comment by Reviewer ohZY
comment: Dear AC, dear Reviewers, 

We greatly appreciate the discussion of our work and we cherish the feedback! In our view, the presented paper makes profound technical contributions and accordingly contains an above-average amount of content (see next paragraph). We made the conscious decision not to withhold any significant contribution from the main body and thus make sure that our work is reviewed in its entirety. While this enforces a more compact writing style, we would like to highlight that, indeed, none of the presented aspects was deemed unnecessary by the reviewers. Instead, reviewers EZD5, H69b, and rAa9 uniformly acknowledged our theoretical foundations/analysis of our Spatio-Spectral Graph Neural Networks (S$^2$GNNs), along with our method's general applicability and strong empirical results. We are moreover certain that the 11% extra space in the camera-ready version will suffice for the discussed adjustments to increase explicitness and, hence, address all the changes suggested by the reviewers.

Arguably, in comparison to impactful works like ChebNet [1] or GCN [2], our technical contributions are of much greater scope. In essence, ChebNet [1] made the Chebyshev-polynomial coefficients of [4] learnable, GCN [2] uses a fixed parametrization again but stacks multiple layers interleaved with non-linearities. Instead, we provide a new modeling paradigm rooted in approximation theory and signal processing that goes beyond pure recursive neighborhood aggregations and unlocks new design spaces.

We appreciate further feedback that would benefit our work.

[1] Defferrard et al. "Convolutional Neural Networks on Graphs with Fast Localized Spectral Filtering" (2016)

[2] Kipf & Welling "Semi-Supervised Classification with Graph Convolutional Networks" (2017)

[3] Veličković et al. "Graph Attention Networks" (2017)

[4] Hammond et al. "Wavelets on Graphs via Spectral Graph Theory" (2009)
We thank reviewer ohZY for acknowledging that we resolved his concerns in our rebuttal and reviewer rAa9 for the additional suggestion on including "stronger arguments in the Introduction". We propose to extend the suggested change in our rebuttal to reviewer ohZY:
> Conversely, spectral filters act globally ($p_{\max}$), even with truncation of the frequency spectrum ($\lambda_{\text{cut}}$). The truncation of the frequency spectrum is required for efficiency. Yet, the combination of spatial and spectral filters provably leverages the strengths of each parametrization. **On an intuitive level, the spectral filter extends the entirely global aggregation of a so-called virtual node to a graph-adaptive intra- and inter-cluster message passing (see Fig. 3 & 4). In contrast, message passing recursively aggregates over local neighborhoods.** Utilizing the strengths of both parametrizations ...

We would appreciate clarification if this does not fully address reviewer rAa9's suggestion. Moreover, we would be grateful for further suggestions from reviewer ohZY for improving the paper's content, methodology, or impact s.t. we can strengthen our contribution and score.
I thank the authors for their well put response. I have increased my score.

Also, thank you AC and reviewers for engaging in the discussion.
We appreciate the reviewer's feedback and thank them for updating their score! We would appreciate if the reviewer has any further specific suggestions for improving the paper's content, methodology, or impact, s.t. we can strengthen our contribution and score.
Thank you for the author's response. Considering the comments from other reviewers, I will raise the score to 6. However, the writing of the paper still needs improvement.
I thank the AC for their involvement in the rebuttal discussion period.

Regarding the review:

- I believe weaknesses 1 and 3 are fair concerns. However, in my opinion, the Authors have properly addressed them in their rebuttal.

- Regarding weakness 2 - I believe that the Authors motivate their approach very well throughout the paper. Nevertheless, including stronger arguments in the Introduction section might make the paper easier to follow. The Authors seem to acknowledge this in their rebuttal.

Overall, I believe that the concerns raised by the Reviewer are legitimate but have been addressed by the Authors. Moreover, in my opinion, the raised concerns are minor in the context of the proposed work. The overall work is on combining spectral convolutions and message-passing, analyzing why the combination works well in practice, especially in modeling long-range relationships. Therefore, I would argue that the initial concerns are insufficient to justify a rejection (3).

Finally, while I can't make any assumptions regarding whether the Authors' responses satisfied Reviewer ohZY, I believe that the authors did a good job answering the questions and clarifying the concerns raised by the Reviewer.

Nevertheless, I am aware that my review is the most positive one, and I believe that input from the other Reviewers would be valuable for this discussion.
Dear Reviewer and Authors,

Since this paper has low and high scores, we should discuss the paper in more detail. 

Specifically, can the reviewers please discuss the review of ohZY with the lowest score and the response from the authors? Do you agree about the weaknesses 1-3? Are you satisfied with the response?
Please feel free to share any additional thoughts you might have. 

Thank you
AC
Dear Reviewer ohZY

Could you please read the rebuttal and respond to the authors?
We are approaching the end of the discussion period, and your feedback will be critical.

Thank you! 
AC
We appreciate the reviewer's feedback and thank them for updating their score! Since we resolved all weaknesses/concerns, we would appreciate if the reviewer has any further suggestions for improving the paper's content, methodology, or impact. We would greatly appreciate their input to strengthen our contribution and score.
Thank you for your response. All of my concerns have been addressed. I will arise my score accordingly.
We thank the reviewer for raising the score to 8 (strong accept) and the strong support on the concerns of the other reviewers! We value the suggestions very much and will take meticulous care on addressing them in the camera ready version.
I want to thank the authors for their rebuttal. I believe they have addressed my concerns and other reviewers' concerns well.

As per the initial review, I have updated my score to an 8 (strong accept). I am keeping my confidence to 4 since there might be some details I missed or some parts of the paper that I might have misunderstood.

Thank you for the paper. I enjoyed reading it, and I believe that the work is very solid.