Paper deep dive
TwinS-GCN: Spectral conjugate for Spectral Graph Convolutional Networks
Chun Hei Michael Chan, Flavia Petruso, Dimitri Van De Ville
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Graph convolutional networks propagate information by repeated local aggregation through a graph shift operator; i.e., a $K$-layer network reaches $K$ hops neighborhood. On the one hand, such spreading can lead to oversmoothing. On the other hand, long-range dependencies demand the depth. Transporting information on long distances and without attenuation requires the shift to distinguish a direction of flow, which a symmetric operator cannot perform but a directed one can fulfill. A natural way to extract pure-directionality is to take the skew-symmetric part of the shift operator through the Cartesian split, which, however, generally does not commute with the shift itself, meaning that the filters built on it are not shift-invariant. We instead use the spectral conjugate; i.e., the image of the operator under $\tau:z\mapsto \bar{z}$, which commutes with the shift and splits it into a dissipative and a non-dissipative part. Two filter families follow: a sum filter, whose non-dissipative component transports signal without energy loss, and a ratio filter, ratio in the pair of components rather than polynomial in the shift. Both arise from non-holomorphic kernels, placing them outside the holomorphic class underlying classical spectral convolution. On the directed cycle, the ratio filter becomes an IIR filter with global impulse response, for which we prove a long-range reach gap against every degree-$K$ polynomial filter. Chebyshev reparameterization gives stable vertex-domain layers with real coefficients, yielding TwinS-GCN, which solves graph transfer tasks at reduced depth and is competitive with state-of-the-art graph convolutional networks on node classification benchmarks.
Tags
Links
- Source: https://arxiv.org/abs/2609.32940v1
- Canonical: https://arxiv.org/abs/2609.32940v1
Trouble viewing inline? Open PDF directly →
Full Text
175,546 characters extracted from source content.
Expand or collapse full text
TwinS-GCN: Spectral conjugate for Spectral Graph Convolutional Networks Chun Hei Michael Chan Affiliation: Neuro-X Institute and Department of Electrical Engineering Flavia Petruso Affiliation: École polytechnique fédérale de Lausanne (EPFL) Dimitri Van De Ville Affiliation: Geneva, 1202, Switzerland Abstract Graph convolutional networks propagate information by repeated local aggregation through a graph shift operator; i.e., a K-layer network reaches K hops neighborhood. On the one hand, such spreading can lead to oversmoothing. On the other hand, long-range dependencies demand the depth. Transporting information on long distances and without attenuation requires the shift to distinguish a direction of flow, which a symmetric operator cannot perform but a directed one can fulfill. A natural way to extract pure-directionality is to take the skew-symmetric part of the shift operator through the Cartesian split, which, however, generally does not commute with the shift itself, meaning that the filters built on it are not shift-invariant. We instead use the spectral conjugate; i.e., the image of the operator under τ:z↦z¯τ:z z, which commutes with the shift and splits it into a dissipative and a non-dissipative part. Two filter families follow: a sum filter, whose non-dissipative component transports signal without energy loss, and a ratio filter, ratio in the pair of components rather than polynomial in the shift. Both arise from non-holomorphic kernels, placing them outside the holomorphic class underlying classical spectral convolution. On the directed cycle, the ratio filter becomes an IIR filter with global impulse response, for which we prove a long-range reach gap against every degree-K polynomial filter. Chebyshev reparameterization gives stable vertex-domain layers with real coefficients, yielding TwinS-GCN, which solves graph transfer tasks at reduced depth and is competitive with state-of-the-art graph convolutional networks on node classification benchmarks. 1 Introduction Graph convolutional networks are among the most widely used architectures for learning on relational data, owing to their efficiency, being light weight and to their increased interpretability given the spectral theory that underpins them (Bruna et al., 2014; Defferrard et al., 2016; Kipf and Welling, 2017). However, two limitations are well documented: node representations tending to converge toward a constant as depth increases, known as oversmoothing (Li et al., 2018; Oono and Suzuki, 2020; Cai and Wang, 2020; Rusch et al., 2023) and a K-layer network aggregating over at most K hops neighborhood, which requires sufficient depth to capture long-range dependencies (Topping et al., 2022; Di Giovanni et al., 2023). Residual connections, normalization, and rewiring (Li et al., 2019; Alon and Yahav, 2021) serve as effective solutions for both obstacles, yet, because our goal is to address their common cause, we do not compare against this line of work. The reasons for these limitations are largely related to the nature of the operator used in the convolutional networks. For a symmetric graph operator, high frequencies tend to be attenuated due to its real spectrum, leading to progressively less distinctive features since spectral convolution strengthens this attenuation when the number of convolution layers (depth) increases (Oono and Suzuki, 2020; Cai and Wang, 2020). Escaping this requires eigenvalues on the imaginary axis, an established criterion for stable non-dissipative dynamics (Haber and Ruthotto, 2018; Chang et al., 2019), which A-DGN and its successors obtain by imposing skew-symmetry on the learnable weights (Gravina et al., 2023; Heilig et al., 2025; Hariri et al., 2026). For directed graphs, the graph shift already possesses non-zero imaginary part as it is non-symmetric (Singh et al., 2016; Marques et al., 2020a). A natural way to extract skew-symmetry and thus imaginary spectrum is the Cartesian split, which separates the shift operator into its symmetric and skew-symmetric parts. We show, however, that neither part commutes with the shift unless the graph is normal—which directed graphs generally are not (Trefethen, 2005; Asllani et al., 2018)—so the induced filters are not shift-invariant. To overcome this limitation while still obtaining imaginary spectrum, we instead use the spectral conjugate (Nevanlinna, 2018), which commutes with the shift on any graph, coincides with the transpose when the shift is normal, and remains defined for non-diagonalizable shifts. On this decomposition, we build learnable graph convolutional layers from filters in (Chan et al., 2026a). Composing the resulting layers gives the TwinS-GCN framework. Importantly, since the underlying scalar functions (kernels) are not complex-differentiable, the construction falls outside the holomorphic functional calculus classical spectral filters rely on (Higham, 2008), introducing a new class of spectral graph convolutions. Contributions. (1) We propose the spectral conjugate, rather than the transpose, to split the directed Laplacian Singh et al. (2016) into dissipative and non-dissipative components while preserving shift invariance. (2) We construct the sum and ratio filters from the resulting components, prove they are both non-holomorphic, and derive their vertex-domain forms via the Cauchy-Pompeiu formula. (3) We prove that on the directed cycle the ratio filter is an IIR filter with global impulse response, and establish a reach gap against every degree-K polynomial filter in the shift. (4) We reparameterize both filters in Chebyshev bases, yielding stable vertex-domain layers with real coefficients. (5) TwinS-GCN solves graph transfer tasks at a fraction of the depth conventional polynomial filters require, preserves feature distinctiveness with depth, and is competitive with state-of-the-art GCNs on general benchmarks for node classification. 2 Dissipative, Asymmetric, and Spectral Conjugate Filters 2.1 Why polynomial spectral filters attenuate Let G=(,ℰ,)G=(V,E,W)11 1 Notations are provided in Appendix A be a weighted graph with N=||N=|V| nodes and weight matrix ∈ℝ≥0N×NW ^N× N_≥ 0, and let ∈ℝN×NL ^N× N be the corresponding graph shift operator (Ortega, 2022) with spectrum σ()=λn⊂ℂσ(L)=\ _n\ . Denote its Jordan decomposition as =(+)−1L=P( +N)P^-1, where diag()=σ()diag( )=σ(L) and N is the nilpotent part. The standard spectral filter is the matrix h()h(L) induced by a scalar kernel h that is analytic on a neighborhood ⊂ℂD of σ()σ(L), on the Jordan block nJ_n of size m attached to λn _n, h(n)=∑k=0m−1h(k)(λn)k!nk,h(J_n)= _k=0^m-1 h^(k)( _n)k!\,N_n^k, (1) placing h(λn)h( _n) on the diagonal and its complex derivatives h(k)(λn)h^(k)( _n) on the successive superdiagonals (Sandryhaila and Moura, 2014). Filtering a signal reads as =h()y=h(L)x, and l layers apply h()lh(L)^l. Since nN_n is nilpotent, the sum terminates, so the gain of each spectral component is |h(λn)||h( _n)| up to a factor polynomial in m. The standard GCN layer can be read as the degree one polynomial h(λ)=1−tλh(λ)=1-tλ (Chamberlain et al., 2021), the first-order truncation of the heat kernel h(λ)=e−tλh(λ)=e^-tλ whose gain is upper bounded by e−tℜ(λn)e^-t ( _n) (Appendix C.5). Components with ℜ(λn)>0 ( _n)>0 thus decay geometrically with depth and only those with ℜ(λn)=0 ( _n)=0 survive: as l→∞l→∞ the representation collapses onto the ℜ(λ)=0 (λ)=0 invariant subspace, which is the spectral form of oversmoothing. The following observation shows that under a symmetric shift this subspace is uninformative. Symmetric shifts and non-dissipative propagation. A real symmetric S has σ()⊂ℝσ(S) , so ℜ(λn)=λn ( _n)= _n for every n; and for a Laplacian on a connected graph the eigenvalue 00 is simple with eigenvector 1 (Fiedler, 1973; Chung, 1997). Under a symmetric shift, the constant signal is therefore the only mode with ℜ(λ)=0 (λ)=0, and the only one preserved under arbitrary depth. Real symmetric operators therefore do not allow for non-dissipative propagation. For directed graphs, the operator is asymmetric, enables purely imaginary eigenvalues, and consequently allows non-dissipative propagation. 2.2 What about the Cartesian split Specify G to be a directed graph, and consider the directed Laplacian =−L=D-W with real weight matrix W and =diag()D=diag(W1) the in-degree matrix (Singh et al., 2016). L is non-symmetric, and its spectrum is complex with ℑ(λ)≠0 (λ)≠ 0 in general. To isolate the asymmetric structure, a natural approach is the Cartesian split =12(+T)+12(−T)L= 12(L+L^T)+ 12(L-L^T), whose second term is skew-symmetric and thus has purely imaginary spectrum. Besides pure imaginary spectrum, it is crucial for both terms to be linear shift-invariant (LSI) — to commute with the shift operator (Sandryhaila and Moura, 2014) — to enable the analysis of graph signals in the same frequency domain, and thereby the interpretation of the features. Additionally filters based on either components are LSI w.r.t. to L only when they the components themselves are. Proposition 1 (Cartesian split breaks shift invariance). The Cartesian components 12(±T) 12(L±L^T) commute with L if and only if L is normal. Even when diagonalizable, directed Laplacians are usually non-normal (Trefethen, 2005; Asllani et al., 2018). Hence the Cartesian split yields components that are not LSI. Thereby, a different split satisfying commutation with L with summands respectively carrying spectrum on real σR()=ℜ(λn) _R(L)=\ ( _n)\ and imaginary axis σI()=jℑ(λn) _I(L)=\j ( _n)\, is desirable. This leads us to naturally consider the spectral conjugate (Nevanlinna, 2018). 2.3 Non-holomorphic functional and the spectral conjugate The spectral conjugate of L is the matrix induced by the non-analytic kernel τ:z↦z¯τ:z z (Nevanlinna, 2018), non-analytic since its Wirtinger derivative ∂¯τ:=∂τ/∂z¯=1 ∂τ:=∂τ/∂ z=1 (Appendix C.4). Yet in equation 1, h is required to be analytic, equivalently holomorphic. Therefore to evaluate τ, we leverage non-holomorphic functional calculus. First recall that holomorphic functional calculus (Auscher et al., 1997) presents the evaluation of an analytic function h at λ∈Ωλ∈ via the Cauchy integral h(λ)=12πi∮Γh(z)z−λz,h(λ)= 12π i _ h(z)z-λ\,dz, (2) where Γ is any contour enclosing λ in its interior. Lifting this to the matrix setting results in the standard graph filters (Appendix C.3) and classical spectral graph filters relying on canonical, Chebyshev, or Cayley polynomials falls within this framework (Appendix C.4). Evaluating τ instead follows from the Cauchy-Pompeiu formula (Bell, 2015), which extends the Cauchy integral to smooth but non-holomorphic functions: h(λ)=12πi∮Γh(z)z−λz+12πi∬Ω∂¯h(z)z−λz∧z¯,h(λ)= 12π i _ h(z)z-λ\,dz+ 12π i _ ∂h(z)z-λ\,dz d z, (3) where Ω is the region enclosed by Γ . Evaluating τ at L through equation 3, as in Nevanlinna (2018), yields the subsequent definition of spectral conjugate. Definition 2.1. (Nevanlinna, 2018, Def. 2.5 – Spectral conjugate) The spectral conjugate of L under the non-holomorphic function τ:z↦z¯τ:z z is given by ♯:=τ()=τ()−1=¯−1,L :=τ(L)=Pτ(J)P^-1=P P^-1, (4) where τ()τ(J) is composed of the Jordan blocks τ(n)=λ¯nτ(J_n)= λ_nI. Theorem 1 (Conjugate properties). ♯L commutes with L, is real, and coincides with TL^T if and only if L is normal. Corollary 1 (Shift-invariant spectral split). The decomposition =12(+♯)⏟∘+12(−♯)⏟↑L= 12 (L+L )_L + 12 (L-L )_L (5) is exact, both components commute with L and with each other, are real, and σ(∘)=σR()σ(L )= _R(L), σ(↑)=σI()σ(L )= _I(L). For normal L, it reduces to the Cartesian split. Thereby we term ∘L and ↑L dissipative and asymmetric operators, respectively. In practice, ♯L is computed with the Schur–Parlett algorithm (Parlett, 1976; Davies and Higham, 2003), which evaluates a matrix function on the Schur form instead of the (numerically unstable) Jordan form. Its cost is dominated by the Schur decomposition which scales as O(N3)O(N^3). 2.4 Spectral filters: the sum and ratio kernels Recalling equation 5, one can design kernels that operate separately on σR() _R(L) and σI() _I(L). We consider the two kernels of (Chan et al., 2026a), Sum kernel. Separate polynomial expansions in the real and imaginary parts, hsum(z)=∑k=0K1[k](z+z¯2)k+∑k=0K2[k](z−z¯2)k,[k],[k]∈ℂ.h_sum(z)= _k=0^K_1a[k]\, ( z+ z2 )^k\;+\; _k=0^K_2b[k]\, ( z- z2 )^k, a[k],b[k] . (6) Ratio kernel. Consider the polynomial in the ratio (z−z¯)/(z+z¯)(z- z)/(z+ z), hrat(z)=[0]+∑k=1K[k](z−z¯z+z¯)k,ℜ(z)≠0,0,ℜ(z)=0,[k]∈ℂ,h_rat(z)=c[0]+ cases _k=1^Kc[k] ( z- zz+ z )^\!k,& (z)≠ 0,\\[10.0pt] 0,& (z)=0, cases c[k] , (7) The case ℜ(z)=0 (z)=0 arises only at the direct current (DC) eigenvalue λ=0λ=0. Setting hrat(0)=[0]h_rat(0)=c[0] removes the division-by-zero instability, and allows filtering of the DC component. We next show that the sum and ratio (equation 6, 7) kernels span a new function class, being non-holomorphic. Theorem 2 (Vertex-domain sum and ratio filters). The kernels hsumh_sum and hrath_rat are non-holomorphic. Their corresponding graph filters admit the vertex-domain representations hsum()=∑k=0K1[k](∘)k+∑k=0K2[k](↑)k,hrat()=∑k=0K[k]( )k, h_sum(L)= _k=0^K_1a[k]\, (L )^k+ _k=0^K_2b[k]\, (L )^k, h_rat(L)= _k=0^Kc[k]\, (L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @setlinewidth @inpgf@ignorespaces @setdash0.0pt @inpgf@ignorespaces @roundcap @inpgf@ignorespaces @roundjoin @inpgf@ignorespaces @stroke @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope )^k, (8) where =↑(∘)DL to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope =L (L )^D and (∘)D(L )^D the Drazin inverse (Appendix C.2) of ∘L . Theorem 2 has two consequences. First, it grounds the vertex-domain forms equation 8 for any L, diagonalizable or not, whereas Chan et al. (2026a) define them only in the diagonalizable case. Second, no polynomial or convergent power series in z represents either kernel. On a fixed graph, h()h(L) still equals a polynomial in L interpolating h on σ()σ(L), but one that is graph-specific, of degree up to N, and ill-conditioned (Appendix C.10). The new class is thus distinguished at the kernel level, and at the matrix level by degree and conditioning. 3 Spectral graph convolution networks 3.1 Graph filter design Canonical polynomial implementations of spectral filters are ill-conditioned in their coefficients, and Chebyshev bases are the typical choice to tackle this limitation (Boyd, 2001); however, they apply to real spectra. We show that equation 8 admits a Chebyshev reparameterization nonetheless. Throughout, TkT_k denotes the degree k Chebyshev polynomial, and is recursively defined as Tk(x)=2xTk−1(x)−Tk−2(x),k≥2,T0(x)=1,T1(x)=x.T_k(x)=2x\,T_k-1(x)-T_k-2(x), k≥ 2, T_0(x)=1,\ T_1(x)=x. (9) Sum filter. The operator ∘L has real eigenvalues λR[n]\ _R[n]\, which we rescale to [−1,1][-1,1]: ~∘=∘−cRR,cR=12(maxnλR[n]+minnλR[n]),R=maxn|λR[n]−cR|+ε. L = L -c_RIR_R, c_R= 12 ( _n _R[n]+ _n _R[n] ), R_R= _n| _R[n]-c_R|+ . (10) The operator ↑L has purely imaginary eigenvalues jλI[n]\j _I[n]\, so −j↑-jL has real spectrum λI[n]\ _I[n]\. Because ↑L is a real-valued matrix, its eigenvalues occur in conjugate pairs, giving cI=0c_I=0; no recentering is needed and ~↑=−j↑/RI,RI=maxn|λI[n]|+ε. L =-j\,L /R_I, R_I= _n| _I[n]|+ . (11) Combining degree-K1K_1 and degree-K2K_2 expansions in the two rescaled operators yields hsum()=∑k=0K1[k]Tk(~∘)+∑k=0K2[k]Tk(~↑),h_sum(L)= _k=0^K_1 α[k]\,T_k( L )+ _k=0^K_2 β[k]\,T_k( L ), (12) with learnable ∈ℂK1+1 α ^K_1+1, ∈ℂK2+1 β ^K_2+1. Ratio filter. The quotient operator =↑(∘)DL to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope =L (L )^D likewise has purely imaginary eigenvalues jλI[n]/λR[n]j _I[n]/ _R[n] for λR[n]≠0 _R[n]≠ 0, so the same substitution applies. Setting =−j Z=-jL to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope gives the real spectrum [n]=λI[n]/λR[n] z[n]= _I[n]/ _R[n], occurring in conjugate pairs, hence cQ=12(maxn[n]+minn[n])=0c_Q= 12 ( _n z[n]+ _n z[n] )=0 and ~=/RQ,RQ=maxn|[n]|+ε, Z=Z/R_Q, R_Q= _n| z[n]|+ , (13) so that hrat()=∑k=0K[k]Tk(~)h_rat(L)= _k=0^K η[k]\,T_k( Z) with learnable ∈ℂK+1 η ^K+1. To observe the effect of reparameterization, we provide the frequency response for each filter hsum(),hrat()h_sum(L),h_rat(L) in Appendix B.4. Theorem 3 (Real-valued twins filters). Real-valued filters expressible by hsum()h_sum(L) are exactly those obtained with [k]∈ℝ α[k] , [2k]∈ℝ β[2k] and [2k+1]∈jℝ β[2k+1]∈ jR for all k. Likewise, real-valued filters expressible by hrat()h_rat(L) are exactly those obtained with [2k]∈ℝ η[2k] and [2k+1]∈jℝ η[2k+1]∈ jR for all k. Despite the odd-order coefficients being purely imaginary, factoring j out for these terms brings back the problem to a real-valued optimization. This effectively handles concerns of memory and compute overhead when optimizing over complex-valued parameter space (Krämer, 2024). 3.2 The TwinS-GCN architecture In what follows, we describe the layer design for both the sum and ratio filters. Sum filter based layer design. We propose a simple layer design that leverages the sum graph filter proposed in equation 8. Denote (l)X^(l) the node features of dimension N×CN× C (nodes times channels) at the l-th layer. The (l+1)(l+1) layer update reads as: (l+1)=ψ(ν∑k=0K1Tk(~∘)(l)k(l,∘)+(1−ν)∑k=0K2jkTk(~↑)(l)k(l,↑)+(l)),X^(l+1)=ψ (ν _k=0^K_1T_k( L )X^(l) _k^(l, )+(1-ν) _k=0^K_2j^k\,T_k( L )X^(l) _k^(l, )\,+B^(l) ), (14) where ψ is a pointwise nonlinear function, k(l,∘) _k^(l, ) and k(l,↑) _k^(l, ) are learnable real weight matrices, (l)B^(l) is a real learnable bias term and ν∈[0,1]ν∈[0,1] is real learnable parameter as introduced in Rossi et al. (2024), controlling importance of dissipative and non-dissipative parts. Composing sum filter based layers, we build the TwinS (S) model. Ratio filter based layer design. Similarly, the update rule for the (l+1)(l+1)-th layer, based on the ratio graph filter proposed in equation 8 is given by: (l+1)=ψ(∑k=0KjkTk(~)(l)k(l)+(l)),X^(l+1)=ψ ( _k=0^Kj^k\,T_k( Z)X^(l) _k^(l)\,+B^(l) ), (15) where k(l) _k^(l) is a learnable real weight matrix. Consecutive ratio filters constitute the TwinS (R) model. We also propose to combine both models leading to TwinS (C) (Appendix C.8). 3.3 Theoretical analysis 3.3.1 Mitigating oversmoothing Oversmoothing arises when repeated graph convolution drives node features toward a constant vector, and is diagnosed by the decay across layers of the Dirichlet energy (Rusch et al., 2023) ((l))=12∑(n,m)∈ℰ¯‖n(l)−m(l)‖2,ℰ¯=ℰ∪(m,n):(n,m)∈ℰ.S(X^(l))= 12\!\! _(n,m)∈ E\!\! \|X^(l)_n-X^(l)_m \|^2, E=E∪\(m,n):(n,m) \. (16) As shown by the sum in equation 16, the Dirichlet energy sees only the symmetrized edge set ℰ¯ E, or in another words the symmetric part of the shift, making it a measure of feature distinctiveness rather than of directed smoothness (Marques et al., 2020b). We retain it in that sense, following Choi et al. (2023), and verify the mitigation empirically in Section 5.1 by tracking ((l))S(X^(l)) across depth. We analyze this in particular on the directed cycle NC_N, where L is normal and the graph Fourier basis is unitary, so for a scalar signal (l)x^(l) Parseval identity gives ((l))=((l))T∘(l)=∑n=0N−1λR[n]|^(l)[n]|2.S(x^(l))= (x^(l) )^\!T\!L \,x^(l)= _n=0^N-1 _R[n]\, | x^(l)[n] |^2. (17) Since ↑L has purely imaginary spectrum σI() _I(L), any spectral function of unit modulus applied to it leaves every |^[n]|| x[n]| unchanged and hence preserves equation 17 exactly. The filter eβ↑e^βL relies on such function, which it acts as ^[n]↦ejβλI[n]^[n] x[n] e^jβ _I[n] x[n], a pure phase rotation with (eβ↑)=()S(e^βL x)=S(x). A degree one sum filter with =(0,0)a=(0,0) and =(1,β)b=(1,β) realizes its first-order truncation, sum(l+1)=(+β↑)(l)⟹(sum(l+1))=∑nλR[n](1+β2λI[n]2)|^(l)[n]|2,x^(l+1)_sum= (I+βL )x^(l)\; \;S (x^(l+1)_sum )= _n _R[n] (1+β^2 _I[n]^2 ) | x^(l)[n] |^2, (18) so the truncation never contracts energy. Figure 1: Isolating ↑L allows hsumh_sum to translate the signal while the polynomial kernel h inseparably moves and flattens it. The mechanism is the complementary spectral action of the two components: ∘L has real non-negative spectrum and attenuates high-frequency modes, whereas ↑L transports features by phase alone and does not attenuate. Therefore increasing the depth does not necessarily imply smoothing. The argument extends beyond NC_N. On any directed graph, the eigenvalues of eβ↑e^βL are ejβλI[n]e^jβ _I[n], all of unit modulus, so no mode is attenuated exponentially: for diagonalizable L, the filter is a pure phase rotation of each Fourier coefficient (no attenuation), and for defective L, the nilpotent part adds at most a factor polynomial in β (Appendix C.5). 3.3.2 Expanding receptive field reach Message-passing networks built on a local shift operator share a structural ceiling: a K-layer network moves information at most K hops (Theorem 5 (i)), so long-range aggregation often demands large depth. We show that the quotient operator L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope alleviates this ceiling, taking NC_N with a pulse-translation task as a setting that isolates receptive field reach. Graph ratio and IIR filters. The operator L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope shares the eigenbasis of the shift but has eigenvalues [n]=jλI[n]/λR[n] μ[n]=j _I[n]/ _R[n]. On NC_N these are a Möbius transform of the shift eigenvalues [n]=ej[n] z[n]=e^j θ[n], which turns the graph ratio filter into a ratio function of the shift. Theorem 4 (IIR correspondence). On NC_N with shift eigenvalues [n]=ej[n] z[n]=e^j θ[n], the eigenvalues [n] μ[n] satisfy the relation [n]=([n]+1)/([n]−1) μ[n]=( z[n]+1)/( z[n]-1). Consequently, a degree-K polynomial filter in L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope has spectral response ∑k=0K[k][n]k=∑k=0K[k]([n]+1)k([n]−1)K−k([n]−1)K=QK([n])([n]−1)K, _k=0^Kc[k]\, μ[n]^k\;=\; _k=0^Kc[k]\,( z[n]+1)^k( z[n]-1)^K-k( z[n]-1)^K\;=\; Q_K( z[n])( z[n]-1)^K, (19) i.e. a ratio or infinite impulse response (IIR) filter of type [K/K][K/K] in the shift, with all K poles located at the DC mode z=1z=1. Its vertex-domain impulse response therefore has global support. Although IIR for NC_N, on general graphs the ratio filter is no longer a ratio function of L. Different from graph IIR filters (Isufi et al., 2016; Liu et al., 2018), it needs no iterative solver to perform filtering and inserts directly into any convolutional layer, retaining finite impulse response (FIR) efficiency while achieving IIR reach. The reach of the ratio filter. We measure how far a filter propagates information. For =h()y=h(L)x with x an impulse at node 00, the influence of node 00 on node d (at distance d ) is h[d,0]:=h()[d,0]=∂[d]/∂[0]h[d,0]:=h(L)[d,0]=∂y[d]/∂x[0] (the sensitivity of Di Giovanni et al. (2023)), and we define the reach as the relative sensitivity Sh(d)=|h[d,0]|maxv∈|h[v,0]|∈[0,1].S_h(d)= |h[d,0]| _v |h[v,0]|∈[0,1]. (20) Theorem 5 (Reach of shift, heat, and ratio filters). On NC_N with nodes =0,…,N−1V=\0,…,N-1\, let d∈∖0d \0\, and assume N≥3K+2N≥ 3K+2. Then: (i) (Shift filter) If h=∑m=0Kγmmh= _m=0^K _mW^m, then h[d,0]=γdh[d,0]= _d for d≤Kd≤ K and Sh(d)=0,for every d>K.S_h(d)=0, every d>K. (21) (i) (Heat filter) If h=e−t∘h=e^-tL with t>0t>0, then h(d)≤(t/2)d!.S_h(d)≤\ (t/2)^dd!. (22) (i) (Ratio filter) If h=hrat:=∑k=0K[k]( )kh=h_rat:= _k=0^Kc[k](L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope )^k then for K<d<N−K<d<N-K Sh(d)≈|p(d/N)|maxv∈|p(v/N)|,p(s)=∑k=1K[k](−1)k−12k!Nk−1Bk(s),S_h(d)≈ |p(d/N)| _v |p(v/N)|\,, p(s)= _k=1^Kc[k]\, (-1)^k-12^kk!\,N^k-1B_k(s), (23) where Bk(s)B_k(s) is the Bernoulli polynomial of degree k. Figure 2: FIR models are obstructed by FIR wall whereas graph ratio filter overcomes it on NC_N. The three regimes are qualitatively distinct. Reach of the shift polynomial (FIR) is identically zero beyond K hops (what we refer to as the FIR wall shown in Figure 2), the heat kernel (IIR) decays super-exponentially, and the ratio filter only polynomially. At K=1K=1, equation 23 gives Sh(d)=|1−2d/N|S_h(d)=|1-2d/N|, nonzero at every distance except for d=N/2d=N/2. Considering (K)=maxhratmind∈Sh(d)G(K)= _h_rat _d S_h(d), i.e., the worst-case reach at order K, we obtain (K)>0.95G(K)>0.95 for K≥4K≥ 4 (Figure 4). While the [K/K][K/K] characterization and bound on sensitivity are specific to NC_N, the escape from the K-hop class remains valid for general directed graph as the filters are not FIR. We empirically show this on general graphs through the graph transfer task (Bodnar et al., 2021; Di Giovanni et al., 2023) in Section 5.1, and further experiments reported in Appendix F. 4 Related Works Spectral graph convolutional networks. FaberNet (HoloNet) (Koke and Cremers, 2024) performs spectral filtering on directed graphs by approximating holomorphic functions through Faber polynomials (Coleman and Smith, 1987; Ellacott, 1983), outperforming MagNet (Zhang et al., 2021), a directed GCN built on the Magnetic Laplacian (Furutani et al., 2020). Both constructs filter directly from a graph operator without separating the mechanisms governing propagation, and in particular without the dissipative and asymmetric split of equation 5. This is reflected in the function class they span: being holomorphic, neither can represent the kernels of Theorem 2. Dir-GNN (Rossi et al., 2024) aggregates separately over in- and out-neighborhoods with distinct weight matrices. Since (,T)↦(12(+T),12(−T))(L,L^T) ( 12(L+L^T), 12(L-L^T) ) is bijective, this realises a Cartesian split, and Proposition 1 applies: the components do not commute with the shift, so the filters are not LSI. Advection-diffusion graph neural networks. For diagonalizable shifts, equation 5 coincides with the graph diffusion and advection operators of Chan et al. (2026a) who study them as signal-processing tools, without learning or propagation analysis. This places our construction alongside architectures of similar concept (Eliasof et al., 2024; Wu et al., 2023), which we include as baselines. Further conceptual comparisons with these networks are provided in Appendix C.9. 5 Experiments Baselines include methods discussed in Section 4: MagNet, FaberNet, Dir-GNN, ADR-GNN and AdvDIFFormer. Methods that instead constrain the weight parametrization (Gravina et al., 2023; Heilig et al., 2025) (discussed in Section 1) act on a complementary axis and can be combined with ours (Appendix C.11), so we treat them as orthogonal rather than competing. 5.1 Receptive field reach and Oversmoothing (a) Accuracy as a function of the number of layers (b) Feature energy at each layer Figure 3: a) Accuracy of different models on the graph transfer task. In the box are models with overlapping accuracy of 11. Models with performance lower than chance level (20%20\%) at all depths are displayed in gray. b) Dirichlet energy computed after each convolution layer before activation, averaged within batch. Only models (1010 layers) with accuracy beyond chance level are shown. Input and target signal’s energy are also shown for each graph. In both case, error bars are over trials. Settings. We evaluate reach and oversmoothing jointly on the graph transfer task (Bodnar et al., 2021; Di Giovanni et al., 2023), adapted to directed graphs (Appendix D). The task involves the propagation of a feature vector from a source node to a target node, which are separated by a certain number of hops in the graph. Following Choi et al. (2023), we track the Dirichlet energy (equation 16) of the node features across layers as a measure of feature distinctiveness. Results. From Figure 3a, TwinS (R) and TwinS (C) transfer the feature accurately with only a few layers on all three graphs, and TwinS (S) does so on the crossed cycle and clique path, demonstrating TwinS-GCN’s ability to capture long-range dependencies at minimal depth. AdvDIFFormer also solves the task at low depth, but at a higher cost in parameters and runtime (Table 7). The remaining polynomial methods (MagNet, FaberNet, Dir-GNN, GCN) and ADR-GNN instead require a depth comparable to the distance to target, or does not solve the task. Figure 3b shows that on each graph at least one TwinS-GCN variant maintains its Dirichlet energy across layers, preserving feature distinctness while solving the task. These results are consistent with reach and non-attenuation properties of Section 3.3 beyond directed cycle. 5.2 Benchmark datasets Settings. We evaluate performances on directed node-classification benchmarks: Chameleon, Squirrel (filtered version as in Platonov et al. (2023)) representing links between Wikipedia pages related to chameleons and squirrels (Rozemberczki et al., 2021), Texas, Wisconsin, and Cornell which are datasets showing links between websites of different universities (Pei et al., 2020). We also compare our model against MLP (Rumelhart et al., 1986) and undirected graphs baselines — GCN (Kipf and Welling, 2017) and GAT (Brody et al., 2021). Results. The results of the benchmark experiments are summarized in Table 1. Our proposed TwinS-GCN model matches the baseline methods across the datasets with best overall rank despite lower accuracy on Squirrel-F, demonstrating its effectiveness in capturing the underlying directed structures and improving node classification accuracy. Reasons for the increased performance could be attributed to the properties of long-range reach and mitigation of oversmoothing, as discussed in Section 5.1, along with the class of function our model leverages. Table 1: Node classification accuracy (in percentage) across benchmark datasets. Best result per dataset in bold, second best underlined. Additional information in Appendix D.2. Model Chameleon-F Squirrel-F Texas Wisconsin Cornell Mean Rank MLP 36.73±4.7136.73± 4.71 36.55±1.8236.55± 1.82 82.16±4.3482.16± 4.34 85.16±4.9585.16± 4.95 75.86¯±3.94¯ 75.86± 3.94 6.706.70 Undir. GCN 40.89±4.1240.89± 4.12 39.47±1.4739.47± 1.47 52.16±5.1652.16± 5.16 48.92±3.0648.92± 3.06 52.70±5.3052.70± 5.30 10.2010.20 GAT 39.21±3.0839.21± 3.08 35.62±2.0635.62± 2.06 58.38±6.6358.38± 6.63 49.41±4.0949.41± 4.09 54.32±5.0554.32± 5.05 10.2010.20 Dir. MagNet 42.36±3.5342.36± 3.53 44.19±1.3944.19± 1.39 85.31¯±4.56¯ 85.31± 4.56 83.40±3.9883.40± 3.98 73.33±4.1773.33± 4.17 4.20¯ 4.20 Dir-GNN 44.13±2.9244.13± 2.92 44.28¯±2.17¯ 44.28± 2.17 70.54±6.5270.54± 6.52 67.39±6.0067.39± 6.00 59.82±6.9059.82± 6.90 6.406.40 FaberNet 42.99±3.3342.99± 3.33 45.13±2.0645.13 2.06 83.15±5.3083.15± 5.30 79.54±5.1879.54± 5.18 72.61±4.7272.61± 4.72 4.604.60 AdvDIFFormer 41.51±3.6341.51± 3.63 41.41±1.5941.41± 1.59 81.80±6.3681.80± 6.36 76.54±4.8476.54± 4.84 72.16±5.5072.16± 5.50 7.207.20 ADR-GNN 46.07±4.0746.07 4.07 44.09±1.9744.09± 1.97 82.16±5.1682.16± 5.16 84.51±4.5484.51± 4.54 71.35±4.3971.35± 4.39 4.904.90 Ours TwinS (S) 42.33±2.8042.33± 2.80 39.61±1.5539.61± 1.55 84.59±3.4384.59± 3.43 87.12¯±4.13¯ 87.12± 4.13 75.14±3.1575.14± 3.15 4.604.60 TwinS (R) 44.24¯±3.23¯ 44.24± 3.23 40.15±1.8340.15± 1.83 84.05±3.9184.05± 3.91 87.39±3.9887.39 3.98 76.76±3.0176.76 3.01 2.80 TwinS (C) 41.64±3.5941.64± 3.59 40.10±1.5840.10± 1.58 85.41±4.7185.41 4.71 86.80±3.7986.80± 3.79 75.68±4.0175.68± 4.01 4.20¯ 4.20 6 Conclusion We introduced TwinS-GCN, a model built upon the spectral conjugate split of the shift. We show that the resulting components and induced filters span a novel class of filters with kernels outside ones reachable by the holomorphic class underlying classical spectral convolutions. We illustrate TwinS-GCN’s capability to preserve feature distinctiveness with far-reaching receptive fields through examples and the graph transfer task. TwinS-GCN also shows on real benchmarks competitive accuracy for node classification. Outlooks and limitations. Although the transpose is not suitable for a shift-invariant split, is inexpensive, well conditioned, and its Cartesian components have spectra distinct from L, which may add expressivity. Performing our splits on both L and TL^T could thus improve accuracy at the cost of interpretability. Likewise, the antisymmetric weight parameterizations act on the channel-mixing matrices and are tangent to our operator-level construction. They can nevertheless be combined with our approach. Computationally, the spectral conjugate requires a Schur decomposition, a one-time cost per graph that nonetheless restricts the present implementation to moderate-size graphs. Krylov-Schur (Stewart, 2002) approximation of the leading modes is a route to scale. Ethics statement No conflict of interest to report nor ethical issues are foreseen in this work. Acknowledgments This work was supported by Swiss National Science Foundation, Sinergia project “Precision mapping of electrical brain network dynamics with application to epilepsy”, under Grant 209470. Reproducibility statement Here are listed the information to ensure reproducibility of our work: • Mathematical concepts are either defined or referenced. Proofs of all Theorems are provided in Appendix B. • Accompanying information for statements and concepts in the main text are provided in Appendix C. • We exactly detail our newly introduced (TwinS-GCN) framework in the main body of our paper Section 3.2. Variants of the framework are discussed in Section 3.2 • The experimental setups for our experiments are described in Appendix D. • The datasets used in our experiments are described in Appendix D.2. • Hyperparameter settings of the methods used in the experiments are detailed in the Appendix D.1, D.2. • Our code will be made publicly available conditionally on acceptance of the paper. References Akiba et al. (2019) T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama Optuna: a next-generation hyperparameter optimization framework. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, Cited by: §D.1. Alon and Yahav (2021) U. Alon and E. Yahav On the bottleneck of graph neural networks and its practical implications. In International Conference on Learning Representations, Cited by: §1. Asllani et al. (2018) M. Asllani, R. Lambiotte, and T. Carletti Structure and dynamical behavior of non-normal networks. Science advances 4 (12), p. eaau9403. Cited by: §1, §2.2. Auscher et al. (1997) P. Auscher, A. McIntosh, and A. Nahmod Holomorphic functional calculi of operators, quadratic estimates and interpolation. Indiana University Mathematics Journal, p. 375–403. Cited by: §C.3, §2.3. Bell (2015) S. R. Bell The cauchy transform, potential theory and conformal mapping. CRC press. Cited by: §2.3. Berlureau et al. (2026) P. Berlureau, A. Hariri, V. Kawasaki-Borruat, M. Zosso, and P. Vandergheynst Advection-diffusion on graphs: a bakry-Émery laplacian for spectral graph neural networks. arXiv preprint arXiv:2602.18141. Cited by: §C.9. Bodnar et al. (2021) C. Bodnar, F. Frasca, N. Otter, Y. Wang, P. Lio, G. F. Montufar, and M. Bronstein Weisfeiler and lehman go cellular: cw networks. Advances in neural information processing systems 34, p. 2625–2640. Cited by: §D.1, §3.3.2, §5.1. Boyd (2001) J. P. Boyd Chebyshev and fourier spectral methods. Courier Corporation. Cited by: §3.1. Brody et al. (2021) S. Brody, U. Alon, and E. Yahav How attentive are graph attention networks?. arXiv preprint arXiv:2105.14491. Cited by: §5.2. Bruna et al. (2014) J. Bruna, W. Zaremba, A. Szlam, and Y. LeCun Spectral networks and locally connected networks on graphs. In International Conference on Learning Representations, Cited by: §1. Cai and Wang (2020) C. Cai and Y. Wang A note on over-smoothing for graph neural networks. In International Conference on Machine Learning, Cited by: §1, §1. Campbell and Meyer (2009) S. L. Campbell and C. D. Meyer Generalized inverses of linear transformations. Classics in Applied Mathematics, SIAM. External Links: Document Cited by: §B.3, §C.2. Chamberlain et al. (2021) B. Chamberlain, J. Rowbottom, M. I. Gorinova, M. Bronstein, S. Webb, and E. Rossi Grand: graph neural diffusion. In International conference on machine learning, p. 1407–1418. Cited by: §2.1. Chan et al. (2026a) C. Chan, A. Cionca, V. Škultéty, and D. Van De Ville Graph diffusion-advection operator for directed graph signal processing. arXiv preprint arXiv:2606.16306. Cited by: §B.4, §C.6, §C.7, §C.9, §C.9, §F.2, §1, §2.4, §2.4, §4. Chan et al. (2025) C. H. M. Chan, A. Cionca, and D. Van De Ville Hilbert transform on graphs: let there be phase. IEEE Signal Processing Letters 32, p. 1725–1729. Cited by: §C.7. Chan et al. (2026b) C. H. M. Chan, A. Cionca, and D. V. D. Ville Statistical testing on directed graphs by surrogate data generation. IEEE Transactions on Signal and Information Processing over Networks (), p. 1–12. External Links: Document Cited by: §C.7. Chang et al. (2019) B. Chang, M. Chen, E. Haber, and E. H. Chi AntisymmetricRNN: a dynamical system view on recurrent neural networks. In International Conference on Learning Representations, Cited by: §1. Chapman (2015) A. Chapman Advection on graphs. In Semi-autonomous networks: Effective control of networked systems through protocols, design, and modeling, p. 3–16. Cited by: §C.9. Choi et al. (2023) J. Choi, S. Hong, N. Park, and S. Cho Gread: graph neural reaction-diffusion networks. In International conference on machine learning, p. 5722–5747. Cited by: §3.3.1, §5.1. Chung (1997) F. R. Chung Spectral graph theory. Vol. 92, American Mathematical Soc.. Cited by: §2.1. Coleman and Smith (1987) J. P. Coleman and R. A. Smith The faber polynomials for circular sectors. mathematics of computation 49 (179), p. 231–241. Cited by: §4. Davies and Higham (2003) P. I. Davies and N. J. Higham A Schur–Parlett algorithm for computing matrix functions. SIAM J. Matrix Anal. Appl. 25 (2), p. 464–485. Cited by: §C.1, §2.3. Defferrard et al. (2016) M. Defferrard, X. Bresson, and P. Vandergheynst Convolutional neural networks on graphs with fast localized spectral filtering. Advances in neural information processing systems 29. Cited by: §1. Di Giovanni et al. (2023) F. Di Giovanni, L. Giusti, F. Barbero, G. Luise, P. Lio, and M. M. Bronstein On over-squashing in message passing neural networks: the impact of width, depth, and topology. In International conference on machine learning, p. 7865–7885. Cited by: §1, §3.3.2, §3.3.2, §5.1. Eliasof et al. (2024) M. Eliasof, E. Haber, and E. Treister Feature transportation improves graph neural networks. In Proceedings of the AAAI conference on artificial intelligence, Vol. 38, p. 11874–11882. Cited by: §C.9, 1st item, §D.2, §4. Ellacott (1983) S. Ellacott Computation of faber series with application to numerical polynomial approximation in the complex plane. Mathematics of Computation 40 (162), p. 575–587. Cited by: §4. Fiedler (1973) M. Fiedler Algebraic connectivity of graphs. Czechoslovak mathematical journal 23 (2), p. 298–305. Cited by: §2.1. Furutani et al. (2020) S. Furutani, T. Shibahara, M. Akiyama, K. Hato, and M. Aida Graph signal processing for directed graphs based on the hermitian laplacian. In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2019, Würzburg, Germany, September 16–20, 2019, Proceedings, Part I, p. 447–463. Cited by: §4. Gravina et al. (2023) A. Gravina, D. Bacciu, and C. Gallicchio Anti-symmetric DGN: a stable architecture for deep graph networks. In International Conference on Learning Representations, Cited by: §C.11, §1, §5. Gravina et al. (2025) A. Gravina, M. Eliasof, C. Gallicchio, D. Bacciu, and C. Schönlieb On oversquashing in graph neural networks through the lens of dynamical systems. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 39, p. 16906–16914. Cited by: Figure 5. Haber and Ruthotto (2018) E. Haber and L. Ruthotto Stable architectures for deep neural networks. Inverse problems 34 (1), p. 014004. Cited by: §1. Hariri et al. (2026) A. Hariri, A. Arroyo, A. Gravina, M. Eliasof, C. Schönlieb, D. Bacciu, X. Dong, K. Azizzadenesheli, and P. Vandergheynst Return of chebnet: understanding and improving an overlooked gnn on long range tasks. Advances in Neural Information Processing Systems 38, p. 136166–136196. Cited by: §1. Heilig et al. (2025) S. Heilig, A. Gravina, A. Trenta, C. Gallicchio, and D. Bacciu Port-hamiltonian architectural bias for long-range propagation in deep graph networks. In International Conference on Learning Representations, Cited by: §1, §5. Higham (2008) N. J. Higham Functions of matrices: theory and computation. SIAM. Cited by: §C.1, §1. Isufi et al. (2016) E. Isufi, A. Loukas, A. Simonetto, and G. Leus Autoregressive moving average graph filtering. IEEE Transactions on Signal Processing 65 (2), p. 274–288. Cited by: §3.3.2. Kipf and Welling (2017) T. N. Kipf and M. Welling Semi-supervised classification with graph convolutional networks. In International Conference on Learning Representations, External Links: Link Cited by: §1, §5.2. Koke and Cremers (2024) C. Koke and D. Cremers Holonets: spectral convolutions do extend to directed graphs. In International Conference on Learning Representations, Vol. 2024, p. 49121–49154. Cited by: §4. Krämer (2024) N. Krämer A tutorial on automatic differentiation with complex numbers. arXiv preprint arXiv:2409.06752. Cited by: §3.1. Levie et al. (2018) R. Levie, F. Monti, X. Bresson, and M. M. Bronstein Cayleynets: graph convolutional neural networks with complex rational spectral filters. IEEE Transactions on Signal Processing 67 (1), p. 97–109. Cited by: §C.4. Li et al. (2019) G. Li, M. Muller, A. Thabet, and B. Ghanem DeepGCNs: can gcns go as deep as cnns?. In Proceedings of the IEEE/CVF international conference on computer vision, p. 9267–9276. Cited by: §1. Li et al. (2018) Q. Li, Z. Han, and X. Wu Deeper insights into graph convolutional networks for semi-supervised learning. In Proceedings of the AAAI conference on artificial intelligence, Vol. 32. Cited by: §1. Liu et al. (2018) J. Liu, E. Isufi, and G. Leus Filter design for autoregressive moving average graph filters. IEEE Transactions on Signal and Information Processing over Networks 5 (1), p. 47–60. Cited by: §3.3.2. Marques et al. (2020a) A. G. Marques, S. Segarra, and G. Mateos Signal processing on directed graphs: the role of edge directionality when processing and learning from network data. IEEE Signal Processing Magazine 37 (6), p. 99–116. Cited by: §1. Marques et al. (2020b) A. G. Marques, S. Segarra, and G. Mateos Signal Processing on Directed Graphs: The Role of Edge Directionality When Processing and Learning From Network Data. IEEE Signal Processing Magazine 37 (6), p. 99–116 (en). External Links: ISSN 1053-5888, 1558-0792, Document Cited by: §3.3.1. Nevanlinna (2018) O. Nevanlinna A non-holomorphic functional calculus and the complex conjugate of a matrix. Linear Algebra and its Applications 537, p. 191–220. Cited by: §B.2, §B.3, §C.10, §1, §2.2, §2.3, §2.3, Definition 2.1, Lemma 7, Lemma 9. Oono and Suzuki (2020) K. Oono and T. Suzuki Graph neural networks exponentially lose expressive power for node classification. In International Conference on Learning Representations, Cited by: §1, §1. Ortega (2022) A. Ortega Introduction to Graph Signal Processing. 1 edition, Cambridge University Press (en). External Links: ISBN 978-1-108-55234-9 978-1-108-42813-2, Document Cited by: §2.1. Parlett (1976) B. N. Parlett A recurrence among the elements of functions of triangular matrices. Linear Algebra and its Applications 14 (2), p. 117–121. Cited by: §C.1, §2.3. Pei et al. (2020) H. Pei, B. Wei, K. C. Chang, Y. Lei, and B. Yang Geom-GCN: geometric graph convolutional networks. In International Conference on Learning Representations (ICLR), External Links: 2002.05287 Cited by: §D.2, §5.2. Pei et al. (2022) H. Pei, B. Wei, K. C. Chang, Y. Lei, and B. Yang Geom-gcn: data release. Note: https://github.com/bingzhewei/geom-gcnCornell node and edge files revised in commit 3c5710 (17 July 2022); accessed 19 September 2026 Cited by: §D.2. Platonov et al. (2023) O. Platonov, D. Kuznedelev, M. Diskin, A. Babenko, and L. Prokhorenkova A critical look at the evaluation of gnns under heterophily: are we really making progress?. arXiv preprint arXiv:2302.11640. Cited by: 1st item, 2nd item, Table 5, §5.2. Rossi et al. (2024) E. Rossi, B. Charpentier, F. Di Giovanni, F. Frasca, S. Günnemann, and M. M. Bronstein Edge directionality improves learning on heterophilic graphs. In Learning on graphs conference, p. 25–1. Cited by: §3.2, §4. Rozemberczki et al. (2021) B. Rozemberczki, C. Allen, and R. Sarkar Multi-scale attributed node embedding. Journal of Complex Networks 9 (2), p. cnab014. Cited by: §5.2. Rumelhart et al. (1986) D. E. Rumelhart, G. E. Hinton, and R. J. Williams Learning representations by back-propagating errors. nature 323 (6088), p. 533–536. Cited by: §5.2. Rusch et al. (2023) T. K. Rusch, M. M. Bronstein, and S. Mishra A survey on oversmoothing in graph neural networks. arXiv preprint arXiv:2303.10993. Cited by: §1, §3.3.1. Sandryhaila and Moura (2014) A. Sandryhaila and J. M. F. Moura Discrete Signal Processing on Graphs: Frequency Analysis. IEEE Transactions on Signal Processing 62 (12), p. 3042–3054 (en). External Links: ISSN 1053-587X, 1941-0476, Document Cited by: §2.1, §2.2. Singh et al. (2016) R. Singh, A. Chakraborty, and B. Manoj Graph fourier transform based on directed laplacian. In 2016 International Conference on Signal Processing and Communications (SPCOM), p. 1–5. Cited by: §1, §1, §2.2. Stanimirović et al. (2015) P. S. Stanimirović, D. Pappas, and V. N. Katsikis Computing the drazin inverse using the schur decomposition. Applied Mathematics and Computation 250, p. 366–381. Cited by: §C.2. Stewart (2002) G. W. Stewart A krylov–schur algorithm for large eigenproblems. SIAM Journal on Matrix Analysis and Applications 23 (3), p. 601–614. Cited by: §6. Topping et al. (2022) J. Topping, F. D. Giovanni, B. P. Chamberlain, X. Dong, and M. M. Bronstein Understanding over-squashing and bottlenecks on graphs via curvature. In International Conference on Learning Representations, Cited by: §1. Trefethen (2005) L. N. Trefethen Spectra and pseudospectra: the behaviour of non-normal matrices and operators. In The graduate student’s guide to numerical analysis’ 98: Lecture notes from the VIII EPSRC Summer School in Numerical Analysis, p. 217–250. Cited by: §1, §2.2. Watson (1944) G. N. Watson A treatise on the theory of bessel functions. 2nd edition, Cambridge University Press. Cited by: §B.6. Wu et al. (2023) Q. Wu, C. Yang, K. Zeng, and M. Bronstein Supercharging graph transformers with advective diffusion. arXiv preprint arXiv:2310.06417. Cited by: §C.9, §4. Zhang et al. (2021) X. Zhang, Y. He, N. Brugnone, M. Perlmutter, and M. Hirn Magnet: a neural network for directed graphs. Advances in neural information processing systems 34, p. 27003–27015. Cited by: §D.2, §4. Appendix A Notations Table 2: Summary of notations used in the main text. Notation Description Graph and signals G=(,ℰ,)G=(V,E,W) Directed weighted graph with nodes V, edges ℰE, weights W N=||N=|V| Number of nodes C Number of feature channels ℰ¯ E Symmetrized edge set, ℰ∪(m,n):(n,m)∈ℰE∪\(m,n):(n,m) \ ∈ℝ≥0N×NW ^N× N_≥ 0 Non-negative weight (adjacency) matrix =diag()D=diag(W1) Diagonal in-degree matrix ∈ℝNx ^N Graph signal (single channel) (l)∈ℝN×CX^(l) ^N× C Node feature matrix at layer l ^[k] x[k] k-th graph Fourier coefficient of x 1 All-ones vector (constant signal) Operators and spectral quantities =−L=D-W Directed graph Laplacian =−1,=−1L=U U^-1,L=PJP^-1 Eigendecomposition and Jordan decomposition of L (⋅)D(·)^D Drazin inverse ,P,\ N Generalized eigenvectors and nilpotent part (Jordan form) U Eigenvector matrix (non-unitary for directed L) Diagonal matrix of complex eigenvalues λk _k λR[n],λI[n] _R[n],\ _I[n] Real and imaginary parts of the n-th eigenvalue ∘L Dissipative operator, 12(+♯) 12(L+L ) ↑L Asymmetric operator, 12(−♯) 12(L-L ) L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope Quotient operator, ↑(∘)DL (L )^D ♯L Spectral conjugate, ¯−1P P^-1 σ()σ(L) Spectrum of L σR() _R(L) Real axis spectrum of L σI() _I(L) Imaginary axis spectrum of L (⋅)S(·) Smoothness measure Kernels, filters and reparameterization τ:ℂ→ℂτ:C Conjugate kernel h:⊂ℂ→ℂh:D Spectral kernel hsum,hrath_sum,\ h_rat Sum and ratio kernels K,K1,K2K,\ K_1,\ K_2 Polynomial degrees ([k],[k])(a[k],\ b[k]), ([k],[k])( α[k],\ β[k]) Sum filter coefficients (canonical and reparameterized) [k],[k]c[k], η[k] Ratio filter coefficients (canonical and reparameterized) [k] θ[k] Standard polynomial filter coefficients Tk(⋅)T_k(·) Chebyshev polynomial of the first kind, degree k ~∘,~↑ L ,\ L Spectrally rescaled ∘L and ↑L =−j ,~Z=-j\,L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope ,\ Z Real-spectrum quotient operator and its rescaling cR,cQc_R,\ c_Q Spectral centering constants R,RI,RQR_R,\ R_I,\ R_Q Spectral scaling radii Architecture and analysis k(l,∘),k(l,↑) _k^(l, ),\ _k^(l, ) Learnable weight matrices at layer l (l)B^(l) Learnable bias at layer l ν,ρν,ρ Learnable parameter shared across layers ψ(⋅)ψ(·) Pointwise non-linear activation μk,zk _k,\ z_k Argument-ordering and shift eigenvalues on NC_N NC_N Directed cycle graph on N nodes Complex analysis Γ,Ω ,\ Integration contour and enclosed region ∂¯=∂/∂z¯ ∂=∂/∂ z Wirtinger derivative (antiholomorphic) ∂=∂/∂z∂=∂/∂ z Wirtinger derivative (holomorphic) ¯ D Closure of domain of definition D j Imaginary unit, j2=−1j^2=-1 Appendix B Proofs of Theorems and Propositions B.1 Proof of Proposition 1 (Cartesian split breaks shift invariance) Lemma 6. Denoting the commutator [,]=− [A,B ]=AB-BA, we have that for any ∈ℝN×NL ^N× N, [,12(±T)]=±12[,T]. [L,\; 12(L±L^T) ]\;=\;± 12 [L,\;L^T ]. (24) Proof. For any ∈ℝN×NL ^N× N, [,12(±T)]=12((±T)−(±T))=±12(T−T)=±12[,T]. [L,\; 12(L±L^T) ]= 12 (L(L±L^T)-(L±L^T)L )=± 12(LL^T-L^TL)=± 12 [L,\;L^T ]. (25) ∎ If L is normal, then T=TLL^T=L^TL and hence [,12(±T)]=0 [L,\; 12(L±L^T) ]=0, so the Cartesian components commute with L. Conversely, if the Cartesian components commute with L, then [,T]=0[L,\;L^T]=0 and L is normal. Hence the Cartesian components commute with L if and only if L is normal. ∎ B.2 Proof of Theorem 1 (Conjugate properties) Proof of main theorem. We refer to Nevanlinna (2018)[Proposition 5.1] for the commutativity of L and ♯L . If L is normal, then ♯=¯H=TL =U U^H=L^T. Conversely, if ♯=TL =L^T, then as ♯L commutes with L, L is normal. What remains is the proof of the realness of ♯L when L is real. Lemma 7. (Nevanlinna, 2018, Eq. 1.9 – Block action) Let h∈N−1(¯)h ^N-1( D) and let n=λn+nJ_n= _nI+N_n be a Jordan block of size m. Then h(n)=∑k=0m−1(∂kh)(λn)k!nk,∂=∂z,h(J_n)= _k=0^m-1 (∂^kh)( _n)k!\,N_n^k, ∂= ∂ z, (26) where ∂ denotes the Wirtinger derivative instead of the complex derivative. This is a generalization of equation 1 to non-holomorphic functions. Only ∂ appears in equation 26 and given that the conjugate kernel τ(z)=z¯τ(z)= z satisfies ∂τ≡0∂τ≡ 0, so every term with k≥1k≥ 1 vanishes and τ(n)=λn¯,τ(J_n)= _n\,I, (27) as defined in Definition 2.1. Lemma 8 (Conjugation). Let h∈N−1(¯)h ^N-1( D). Setting h∗(z):=h(z¯)¯h^*(z):= h( z), we have h∗∈N−1(¯)h^* ^N-1( D) and h()¯=h∗(¯)for every ∈ℂN×N. h(A)=h^*( A) every A ^N× N. (28) Proof. Conjugation exchanges the Wirtinger derivatives, ∂F¯=∂¯F¯∂\, F= ∂F, so (∂h∗)(z)=(∂h)(z¯)¯(∂ h^*)(z)= (∂ h)( z) and, inductively, (∂kh∗)(λ¯)=(∂kh)(λ)¯(∂^kh^*)( λ)= (∂^kh)(λ). If =−1A=PJP^-1 then ¯=¯¯¯−1 A= P\, J\, P^-1, and ¯ J is a Jordan form whose blocks are λ¯n+n λ_nI+N_n, the nilpotent parts being real. Applying equation 26 blockwise, h(n)¯=∑k=0m−1(∂kh)(λn)¯k!nk=∑k=0m−1(∂kh∗)(λn¯)k!nk=h∗(n¯), h(J_n)= _k=0^m-1 (∂^kh)( _n)k!\,N_n^k= _k=0^m-1 (∂^kh^*)( _n)k!\,N_n^k=h^*( J_n), (29) and conjugating h()−1P\,h(J)\,P^-1 gives equation 28. ∎ Applying Lemma 8 to the conjugate kernel: τ∗(z)=τ(z¯)¯=z¯=z¯=τ(z)τ^*(z)= τ( z)= z= z=τ(z). Since L is real, ¯= L=L, and equation 28 gives ♯¯=τ()¯=τ∗(¯)=τ()=♯. L = τ(L)=τ^*( L)=τ(L)=L . (30) A matrix equal to its own conjugate is real. ∎ Proof of corollary. The decomposition =12(+♯)+12(−♯)L= 12 (L+L )+ 12 (L-L ) is exact. Both components commute with L and with each other since ♯L commutes with L. As ♯L is real, both components of the decomposition are also real. Also by definition, σ(∘)=σR()σ(L )= _R(L), σ(↑)=σI()σ(L )= _I(L). Following, when L is normal, this decomposition reduces to the Cartesian split. ∎ B.3 Proof of Theorem 2 (Vertex-domain sum and ratio filters) Non-holomorphicity. Decomposing ℜ(z)=(z+z¯)/2 (z)=(z+ z)/2 and ℑ(z)=(z−z¯)/(2j) (z)=(z- z)/(2j), the Wirtinger derivative ∂¯hsum,∂¯hrat ∂h_sum, ∂h_rat simplifies to: ∂¯hsum(z)=∑k=1K1k[k]2k(z+z¯)k−1+∑k=1K2−k[k](2j)k(z−z¯)k−1≠0, ∂h_sum(z)= _k=1^K_1 k\,a[k]2^k(z+ z)^k-1+ _k=1^K_2 -k\,b[k](2j)^k(z- z)^k-1≠ 0, (31) ∂¯hrat(z)=∑k=1K[k]−2kz(z−z¯)k−1jk(z+z¯)k+1≠0, ∂h_rat(z)= _k=1^Kc[k] -2kz(z- z)^k-1j^k(z+ z)^k+1≠ 0, (32) for non-zero ,,a,b,c. This shows that both kernels are not holomorphic, as their Wirtinger derivative is non-zero. Vertex-domain expressions. Both filters follow from the functional calculus of Nevanlinna (2018) applied to kernels built from the identity and the conjugation τ(z)=z¯τ(z)= z. Lemma 9 ((Nevanlinna, 2018, Cor. 3.5)). Let ∈ℂN×NL ^N× N and ⊂ℂD be open with ⊃σ()D⊃σ(L). The map h↦h()h h(L) from N−1(¯)C^N-1( D) to ℂN×NC^N× N is a homomorphism. Write κ∘=12(id+τ)κ = 12(id+τ) and κ↑=12(id−τ)κ = 12(id-τ), so that linearity of Lemma 9 together with τ()=♯τ(L)=L gives κ∘()=∘κ (L)=L and κ↑()=↑κ (L)=L , and multiplicativity gives (κ∘)k()=(∘)k(κ )^k(L)=(L )^k and (κ↑)k()=(↑)k(κ )^k(L)=(L )^k for every k∈ℕk . Since hsum=∑k=0K1[k](κ∘)k+∑k=0K2[k](κ↑)kh_sum= _k=0^K_1a[k](κ )^k+ _k=0^K_2b[k](κ )^k, linearity yields equation 8 directly: hsum()=∑k=0K1[k](∘)k+∑k=0K2[k](↑)k.h_sum(L)= _k=0^K_1a[k]\,(L )^k+ _k=0^K_2b[k]\,(L )^k. (33) The ratio kernel requires more care, since it involves a quotient and is discontinuous where ℜ(z)=0 (z)=0. We assume throughout that λ=0λ=0 is the only eigenvalue of L on the imaginary axis, which holds for the Laplacian of any connected graph. As σ()σ(L) is finite, fix ε>0 >0 small enough that the discs B(λn,ε)B( _n, ) are pairwise disjoint and κ∘κ does not vanish on any B(λn,ε)B( _n, ) with λn≠0 _n≠ 0, and set =⋃nB(λn,ε)D= _nB( _n, ). Being open and disconnected, D admits kernels prescribed independently on each set, which lets us define ξ(z)=1/κ∘(z),z∈B(λn,ε),λn≠0,0,z∈B(0,ε).ξ(z)= cases1/κ (z),&z∈ B( _n, ),\ _n≠ 0,\\[2.0pt] 0,&z∈ B(0, ). cases (34) On each set ξ is either constant or the reciprocal of a non-vanishing polynomial in (z,z¯)(z, z), so ξ∈N−1(¯)ξ ^N-1( D) and Lemma 9 applies. By construction ξ inverts the nonzero part of σ(∘)σ(L ) and sets the rest to 00, which is the spectral characterisation of the Drazin inverse (Campbell and Meyer, 2009), thence ξ()=(∘)Dξ(L)=(L )^D. Multiplicativity then gives (κ↑ξ)()=κ↑()ξ()=↑(∘)D= ,(κ ξ)(L)=κ (L)\,ξ(L)=L (L )^D=L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope , (35) and, applied repeatedly, (κ↑ξ)k()=( )k(κ ξ)^k(L)=(L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope )^k for all k≥0k≥ 0. Finally hrat=[0]+∑k=1K[k](κ↑ξ)kh_rat=c[0]+ _k=1^Kc[k](κ ξ)^k on D, where on the points with λn≠0 _n≠ 0 we have κ↑ξ=κ↑/κ∘κ ξ=κ /κ , recovering the quotient of equation 7, while on B(0,ε)B(0, ) every term with k≥1k≥ 1 vanishes and the kernel reduces to [0]c[0], matching hrath_rat at ℜ(z)=0 (z)=0. Linearity then gives hrat()=∑k=0K[k]( )k.h_rat(L)= _k=0^Kc[k]\, (L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope )^k. (36) ∎ B.4 Proof of Theorem 3 (Real-valued twins filters) As in Chan et al. (2026a), leveraging the Vandermonde matrix, any real-valued filters expressible by hsum()h_sum(L) and hrat()h_rat(L) i.e. elements of respectively: hsum():(,)∈ℂK1+1×ℂK2+1∩ℝN×Nandhrat():∈ℂK+1∩ℝN×N\h_sum(L):(a,b) ^K_1+1×C^K_2+1\ ^N× N \h_rat(L):c ^K+1\ ^N× N (37) as defined in equation 8, can be expressed as elements of hsum():(,)∈ℝK1+1×ℝK2+1\h_sum(L):(a,b) ^K_1+1×R^K_2+1\ and hrat():∈ℝK+1\h_rat(L):c ^K+1\ guaranteeing their equality — this holds for both diagonalizable and non-diagonalizable L thanks to conjugate eigenvalues occurrence in the Jordan decomposition. Since we need hsum(),hrat()h_sum(L),h_rat(L) to be real, we start with ,,∈ℝa,b,c . By identifying the unique coefficients for the Chebyshev expansions of ~∘,~↑,~ L , L , Z, starting with maximal order K [K]TK(~∘) α[K]T_K( L ) =[K]∑k=0K[k](~∘)k=[K]∑k=0K(K[k]R∑m=0k(km)(∘)m(cR)k−m), = α[K] _k=0^Ka_K[k]( L )^k= α[K] _k=0^K ( a_K[k]R_R _m=0^k km(L )^m(c_RI)^k-m ), (38) [K]TK(~↑) β[K]T_K( L ) =[K]∑k=0K[k](~↑)k=[K]∑k=0K(−j)kK[k]RI(↑)k, = β[K] _k=0^Kb_K[k]( L )^k= β[K] _k=0^K(-j)^k b_K[k]R_I(L )^k, (39) [K]TK(~) η[K]T_K( Z) =[K]∑k=0K[k](~)k=[K]∑k=0K(−j)kK[k]RQ( )k. = η[K] _k=0^Kc_K[k]( Z)^k= η[K] _k=0^K(-j)^k c_K[k]R_Q(L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope )^k. (40) Given TKT_K is of maximal order K, no other terms except for [K] α[K] contribute to the expansion in [K]a[K], therefore we have exactly [K]K[K]R=[K] α[K] a_K[K]R_R=a[K]. Similarly, considering the two remaining expressions [K](−j)KK[K]RI=[K] β[K](-j)^K b_K[K]R_I=b[K], and [K](−j)KK[K]RQ=[K] η[K](-j)^K c_K[K]R_Q=c[K]. Hence the following holds: • [K] α[K] is real, • [K] β[K] is real for even K and purely imaginary for odd K, • [K] η[K] is real for even K and purely imaginary for odd K. For lower order coefficients, we can apply the same argument recursively (while subtracting the previous orders for which we have already determined the coefficients) to show that [k] α[k] is real for all k and [k],[k] β[k], η[k] are real for even k and purely imaginary for odd k. As this procedure can be performed for all ,,∈ℝa,b,c , hsum():(,)∈ℂK1+1×ℂK2+1∩ℝN×N \h_sum(L):( α, β) ^K_1+1×C^K_2+1\ ^N× N =hsum():[k]∈ℝ,[2k]∈ℝ,[2k+1]∈jℝ, =\h_sum(L): α[k] ,\ β[2k] ,\ β[2k+1]∈ jR\, (41) hrat():∈ℂK+1∩ℝN×N \h_rat(L): η ^K+1\ ^N× N =hrat():[2k]∈ℝ,[2k+1]∈jℝ. =\h_rat(L): η[2k] ,\ η[2k+1]∈ jR\. (42) This completes the proof of Theorem 3. ∎ Further intuition on reparameterized matrices. Additional intuition on the nature of the real-matrix condition is provided in the following lemma and expressions, which in one hand show how the Chebyshev polynomial parity interacts with the matrix parity of ~↑ L and ~ Z to yield a real-matrix filter and in the other hand show how the spectral response is decomposed into real and imaginary parts that are even and odd in λI _I respectively. Lemma 10 (Parity of Chebyshev polynomials). Let x∈[−1,1]x∈[-1,1], the k-th degree Chebyshev polynomials TkT_k satisfies Tk(−x)=(−1)kTk(x)⇒even functionk even,odd functionk odd.T_k(-x)=(-1)^kT_k(x) caseseven function&k even,\\ odd function&k odd. cases (43) Proof. For x∈[−1,1]x∈[-1,1], Tk(−x)=cos(karccos(−x))=cos(k(π−arccos(x)))=(−1)kTk(x)T_k(-x)= (k (-x))= (k(π- (x)))=(-1)^kT_k(x). ∎ Lemma 11 (Matrix parity of Chebyshev polynomials). Let Q be a matrix with purely imaginary entries, giving i∈ℝN×Ni even,jℝN×Ni odd.Q^i\;∈\; casesR^N× N&i even,\\ jR^N× N&i odd. cases (44) Parallely, from Lemma 10 the Chebyshev polynomial TkT_k has the same parity as k, so it contains only monomials xix^i with i≡k(mod2)i≡ k 2. Combining with equation 44: Ak=Tk()∈ℝN×Nk even,jℝN×Nk odd.A_k=T_k(Q)\;∈\; casesR^N× N&k even,\\ j\,R^N× N&k odd. cases (45) Proof. Write Tk(x)=∑i≡k(2)aixiT_k(x)= _i≡ k\,(2)a_ix^i with ai∈ℝa_i . Then Ak=∑i≡k(2)aiiA_k= _i≡ k\,(2)a_iQ^i. Each summand aiia_iQ^i lies in ℝN×NR^N× N if i even, in jℝN×NjR^N× N if i odd, and all summands share the same parity as k, so the whole sum is either in ℝN×NR^N× N or jℝN×NjR^N× N. ∎ Frequency response as a function of Chebyshev polynomials. The spectral response at eigenvalue λ[n]=λR[n]+jλI[n]λ[n]= _R[n]+j _I[n] of the sum filter hsum(λ[n])h_sum(λ[n]) can be expressed in terms of Chebyshev polynomials as ℜ(hsum(λ[n])) (h_sum(λ[n]) ) =∑k[k]Tk(xR,n)+∑i[2i]T2i(xI,n), = _k α[k]T_k(x_R,n)+ _i β[2i]\,T_2i(x_I,n), (46) ℑ(hsum(λ[n])) (h_sum(λ[n]) ) =∑i[2i+1]T2i+1(xI,n), = _i β[2i+1]\,T_2i+1(x_I,n), (47) where xR,n=(λR[n]−cR)/RRx_R,n=( _R[n]-c_R)/R_R and xI,n=λI[n]/RIx_I,n= _I[n]/R_I. For the ratio filter hrat(λ[n])h_rat(λ[n]) we have ℜ(hrat(λ[n])) (h_rat(λ[n]) ) =∑i[2i]T2i(xQ,n), = _i η[2i]\,T_2i(x_Q,n), (48) ℑ(hrat(λ[n])) (h_rat(λ[n]) ) =∑i[2i+1]T2i+1(xQ,n), = _i η[2i+1]\,T_2i+1(x_Q,n), (49) where xQ,n=(−jjλI[n]λR[n])/RQx_Q,n= (-j j _I[n] _R[n] )/R_Q. B.5 Proof of Theorem 4 (IIR correspondence) Möbius identity. On the cycle, λR[n]=1−cos[n] _R[n]=1- θ[n] and λI[n]=−sin[n] _I[n]=- θ[n], so [n]=jλI[n]λR[n]=−jsin[n]1−cos[n]=−2jsin([n]/2)cos([n]/2)2sin2([n]/2)=−jcot([n]2). μ[n]= j _I[n] _R[n]= -j θ[n]1- θ[n]= -2j ( θ[n]/2) ( θ[n]/2)2 ^2( θ[n]/2)=-j \! ( θ[n]2 ). (50) Writing [k]=ej[n] z[k]=e^j θ[n] and using cot(θ/2)=j(z+1)/(z−1) (θ/2)=j(z+1)/(z-1), [n]=−j⋅j[n]+1[n]−1=[n]+1[n]−1, μ[n]=-j· j\, z[n]+1 z[n]-1= z[n]+1 z[n]-1, (51) which is equation 50. This is a Möbius transform, projecting from the unit circle onto the imaginary axis. Rational structure equation 19. Substituting equation 50 into a degree-K polynomial, ∑k=0K[k][n]k=∑k=0K[k]([n]+1)k([n]−1)K−k([n]−1)K=QK([n])([n]−1)K, _k=0^Kc[k] μ[n]^k= _k=0^Kc[k]( z[n]+1)^k( z[n]-1)^K-k( z[n]-1)^K= Q_K( z[n])( z[n]-1)^K, (52) where QKQ_K is a polynomial of degree ≤K≤ K obtained by multiplying the denominator. This is a rational function of [n] z[n] of type [K/K][K/K] whose denominator ([n]−1)K( z[n]-1)^K places all K poles at =1 z=1, the DC mode. A rational spectral response corresponds to an IIR (ARMA) vertex-domain filter and hence has non-compact support. ∎ B.6 Proof of Theorem 5 (Reachability gap) Starting notation. Δu[ℓ]=u[ℓ+1]−2u[ℓ]+u[ℓ−1] u[ ]=u[ +1]-2u[ ]+u[ -1] is the discrete Laplacian and k:=( )kδ0g_k:=(L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope )^kδ_0 the k-th kernel, so that k[d]g_k[d] is the impulse response of ( )k(L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope )^k at distance d. Since L is normal on CNC_N, the spectral split gives ∘=−12(+T),↑=12(T−),L =I- 12(W+W^T), L = 12(W^T-W), (53) so that, acting on a signal u∈ℝNu ^N, (∘u)[ℓ]=−12(Δu)[ℓ],(↑u)[ℓ]=12(u[ℓ+1]−u[ℓ−1]):(L u)[ ]=- 12( u)[ ], (L u)[ ]= 12 (u[ +1]-u[ -1] ): (54) a second difference and a centered first difference. Proposition 2 (Identity on ratio operator). ∘ =↑L L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope =L , and =L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope 1=0, T =T1^TL to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope =0^T. Proof. All operators commute, so it suffices to observe their spectral definition. ∎ Lemma 12 (Linear form of 1g_1). 1[ℓ]=2ℓ/N−1g_1[ ]=2 /N-1 for 1≤ℓ≤N−11≤ ≤ N-1, and 1[0]=0g_1[0]=0. Proof. Apply Proposition 2 to δ0δ_0 and use equation 54: ∘1=↑δ0=12(δ−1−δ1)⟺Δ1=δ1−δ−1.L g_1=L δ_0= 12(δ_-1-δ_1) g_1=δ_1-δ_-1. (55) By equation 55 the second difference of 1g_1 vanishes at every vertex other than ℓ=±1 =± 1. A sequence with vanishing second difference on a discrete interval is affine on that interval, so 1[ℓ]=a+bℓg_1[ ]=a+b for 2≤ℓ≤N−22≤ ≤ N-2. Evaluating Δ1 g_1 at ℓ=0 =0 and inserting the affine form gives 1[0]=a+bN2g_1[0]=a+ bN2. An additional condition comes from the fact that L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope sets the zero-th frequency to 00, i.e., ^1[0]=μ[0]=0 g_1[0]=μ[0]=0 and therefore ∑ℓ1[ℓ]=0 _ g_1[ ]=0. Combining these two conditions gives a=−bN2a=- bN2 and 1[0]=0g_1[0]=0. Once again evaluating Δ1 g_1 at ℓ=1 =1 gives 1[2]−21[1]+1[0]=1g_1[2]-2g_1[1]+g_1[0]=1, which with the affine form gives b=2/Nb=2/N. Hence 1[ℓ]=2ℓ/N−1g_1[ ]=2 /N-1 for 1≤ℓ≤N−11≤ ≤ N-1. ∎ Lemma 13 (Degree and leading coefficient). For N≥3k+2N≥ 3k+2 there is a polynomial PkP_k of degree exactly k with k[d]=Pk(d)g_k[d]=P_k(d) on the interval k<d<N−k<d<N-k, whose leading coefficient obeys ak=−2ak−1k,a1=2N,⟹ak=(−1)k−12kNk!.a_k=- 2a_k-1k, a_1= 2N, a_k= (-1)^k-12^kN\,k!. (56) Proof. We proceed by induction on k, Base case k=1k=1: Lemma 12 gives P1(d)=2d/N−1P_1(d)=2d/N-1, of degree 11 with a1=2/Na_1=2/N. Inductive hypothesis: Suppose the statement holds at order k−1k-1, so that k−1=Pk−1g_k-1=P_k-1 on k−1<d<N−k+1k-1<d<N-k+1 with degPk−1=k−1 P_k-1=k-1 and leading coefficient ak−1≠0a_k-1≠ 0. Applying Proposition 2 to k= k−1g_k=L to3.81pt to8.37pt -2.47626pt to0.0pt @begingroup@_scopebegin=1 @begingroup@stroke=#000000 @begingroup@fill=#000000 @setlinewidth @begingroup@stroke-width=0.4pt @inpgf@ignorespaces to0.0pt @begingroup@_scopebegin=1 @stroke @drawpath@unclippedM 0 -3.15 L 0 7.24fill:none @inpgf@ignorespaces @inpgf@ignorespaces @begingroup@_scopebegin=1 @transformcm0.01.0-1.00.00.0pt5.23047pt @begingroup@transform=matrix(0.0 1.0 -1.0 0.0 0 7.24) @begingroup@_scopebegin=1 @begingroup@stroke-width=0.32pt @begingroup@stroke-dasharray=none,stroke-dashoffset=0.0pt @begingroup@stroke-linecap=round @begingroup@stroke-linejoin=round @drawpath@unclippedM -1.66 2.21 C -1.52 1.38 0 0.14 0.42 0 C 0 -0.14 -1.52 -1.38 -1.66 -2.21fill:none @closescope @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @stroke @drawpath@unclippedM 0 1.97 M 2.36 1.97 C 2.36 3.27 1.3 4.33 0 4.33 C -1.3 4.33 -2.36 3.27 -2.36 1.97 C -2.36 0.66 -1.3 -0.39 0 -0.39 C 1.3 -0.39 2.36 0.66 2.36 1.97 Z M 0 1.97fill:none @inpgf@ignorespaces @closescope @inpgf@ignorespaces @inpgf@ignorespaces @inpgf@ignorespaces @discardpath @closescope @closescope g_k-1 and using equation 54, −12Δk=↑k−1.- 12\, g_k=L g_k-1. (57) Let d lie in the interval k<d<N−k<d<N-k. Then both d−1d-1 and d+1d+1 lie in k−1<d<N−k+1k-1<d<N-k+1, where the inductive hypothesis applies, so by equation 54 the right-hand side of equation 57 equals Q(d):=12(Pk−1(d+1)−Pk−1(d−1)).Q(d):= 12 (P_k-1(d+1)-P_k-1(d-1) ). (58) Expanding, the terms of degree k−1k-1 cancel and Q has degree exactly k−2k-2 with leading coefficient (k−1)ak−1(k-1)a_k-1. On the interval, equation 57 reads Δk=−2Q g_k=-2Q and to it we seek a particular solution that is a polynomial of degree k with leading coefficient b. Since Δ(bdk)=bk(k−1)dk−2+O(dk−3) (b\,d^k)=b\,k(k-1)\,d^k-2+O(d^k-3), matching the coefficients of dk−2d^k-2 gives bk(k−1)=−2(k−1)ak−1⟹b=−2ak−1k≠0.b\,k(k-1)=-2(k-1)a_k-1 b=- 2a_k-1k≠ 0. (59) The homogeneous solutions of Δu=0 u=0 on a discrete interval are the affine sequences, which for k≥2k≥ 2 cannot cancel the term of degree k, so ak≠0a_k≠ 0. Hence kg_k agrees on the interval with a polynomial of degree exactly k, and ak=ba_k=b, which is equation 56. The hypothesis N≥3k+2N≥ 3k+2 guarantees that the interval contains more than k+1k+1 vertices, so the polynomial is not degenerate. Unrolling the recursion leads to equation 56. ∎ Lemma 14 (Bernoulli form). Let s=d/Ns=d/N and ck=(−1)k−12k/k!c_k=(-1)^k-12^k/k!. Then, uniformly on k<d<N−k<d<N-k, k[d]=Nk−1ckBk(s)+O(Nk−3),g_k[d]=N^k-1c_k\,B_k(s)+O\! (N^k-3 ), (60) with equality for k=1k=1. Proof. Set P~k(s):=Pk(Ns)/Nk−1 P_k(s):=P_k(Ns)/N^k-1, the k-th kernel expressed in the fraction distance s=d/Ns=d/N. Its leading term is ak(Ns)k/Nk−1=Nakska_k(Ns)^k/N^k-1=Na_ks^k, so by equation 56 its leading coefficient is Nak=(−1)k−12k/k!=:ckNa_k=(-1)^k-12^k/k!=:c_k, which is independent of N. This normalization is what allows a limit to be taken when computing derivatives. Under d=Nsd=Ns, the second difference acts as N−2d2/ds2N^-2\,d^2/ds^2 and the centered first difference as N−1d/dsN^-1\,d/ds, each with relative error O(N−2)O(N^-2) on polynomials of fixed degree. Hence differentiating over s, equation 57 becomes Pk′=−2NPk−1′P_k =-2N\,P_k-1 , and substituting Pk=Nk−1P~kP_k=N^k-1 P_k cancels every power of N, leaving P~k′=−2P~k−1′,P~1(s)=2s−1=2B1(s). P_k =-2\, P_k-1 , P_1(s)=2s-1=2B_1(s). (61) Try P~k=ckBk P_k=c_kB_k. The Bernoulli polynomials satisfy Bk′=kBk−1B_k =k\,B_k-1, hence Bk′=k(k−1)Bk−2B_k =k(k-1)B_k-2, so equation 61 holds if and only if ckk(k−1)=−2ck−1(k−1)⟺ck=−2ck−1k,c_k\,k(k-1)=-2\,c_k-1(k-1) c_k=- 2c_k-1k, (62) which is the recursion equation 56 of Lemma 13. With c1=2c_1=2 from equation 61, this gives ck=(−1)k−12k/k!c_k=(-1)^k-12^k/k!. To fully determine solution of the second-order equation 61 it first satisfies, ∫01P~k=0 _0^1 P_k=0, which follows from ^k[0]=μ[0]k=0 g_k[0]=μ[0]^k=0. The second is P~k(1−s)=(−1)kP~k(s) P_k(1-s)=(-1)^k P_k(s), which follows from k[N−ℓ]=(−1)kk[ℓ]g_k[N- ]=(-1)^kg_k[ ]. This equality is surmised from the expression in Lemma 12 where 1g_1 is an odd sequence, furthermore k=1∗kg_k=g_1^*k, and the convolution of two odd sequences is even, so parity alternates with k. Both conditions are satisfied for BkB_k, by its vanishing mean and by the reflection formula Bk(1−s)=(−1)kBk(s)B_k(1-s)=(-1)^kB_k(s). The ansatz of Bernoulli polynomials therefore holds. The only approximation is the replacement of differences by derivatives when considering polynomials, leading to O(N−2)O(N^-2) term. At k steps relative error remains O(N−2)O(N^-2), hence the recursion gives equation 60. For k=1k=1 no approximation is made, since Lemma 12 is exact. ∎ Proof of Theorem 5. (i): Since mδ0=δmW^mδ_0=δ_m, the kernel of h=∑m≤Kγmmh= _m≤ K _mW^m is ∑m≤Kγmδm _m≤ K _mδ_m. It therefore equals γd _d at d≤Kd≤ K and vanishes identically for d>Kd>K, hence Sh(d)=0S_h(d)=0 there. The coefficients γm _m are unconstrained, so the response on the first K hops is arbitrary. (i): By equation 53, e−t∘=e−texp(t2(+T))e^-tL =e^-t ( t2(W+W^T) ). The generating identity et2(z+z−1)=∑dId(t)zde t2(z+z^-1)= _dI_d(t)z^d, where IdI_d is the modified Bessel function of first kind (Watson, 1944), leads to the impules response h[d,0]h[d,0] as e−tId(t)e^-tI_d(t), and the peak is at d=0d=0 because I0≥IdI_0≥ I_d for all d. For the bound of Sh(d)=Id(t)/I0(t)S_h(d)=I_d(t)/I_0(t), compare the defining series termwise: the m-th term of Id(t)I_d(t) is (t/2)d+2m/(m!(d+m)!)(t/2)^d+2m/ (m!\,(d+m)! ) and that of (t/2)d!I0(t) (t/2)^dd!I_0(t) is (t/2)d+2m/(d!(m!)2)(t/2)^d+2m/ (d!\,(m!)^2 ), so their ratio is d!m!/(d+m)!≤1d!\,m!/(d+m)!≤ 1. Summing over m gives Id(t)≤(t/2)d!I0(t)I_d(t)≤ (t/2)^dd!I_0(t), and Sh(d)≤(t/2)d!S_h(d)≤ (t/2)^dd!, which is the claim. (i): The term [0]N η[0]I_N has kernel [0]δ0 η[0]δ_0, which vanishes on K<d<N−K<d<N-K; only the terms with k≥1k≥ 1 contribute there. On the same interval K<d<N−K<d<N-K, Lemma 14 applies to each kg_k with k≤Kk≤ K. A weighted sum kg_k with weights [k] η[k] gives p(s)p(s) in equation 23. Substituting p(s)≈h[d,0]p(s)≈ h[d,0] into Sh(d)S_h(d), first yields |p(s)|/maxq∈[0,1]|p(q)||p(s)|/ _q∈[0,1]|p(q)|. The remaining discrepancy is between the maximum over the vertices v∈v and the supremum over the continuum [0,1][0,1], which for polynomials of fixed degree is O(N−2)O(N^-2). (a) Sensitivity profiles (b) Best filter sensitivity across distances (c) (K)G(K) function Figure 4: a) Sensitivity profiles computed for shift, heat (t=5t=5), heat (t=20t=20), and ratio (K=1K=1) filters. b) Plot of the optimal filter sensitivity across distances. c) Plot of the (K)G(K) function. Throughout N is set to 200200. Appendix C Accompanying information C.1 Spectral conjugate definition Definition C.1 (Spectral conjugate). Let ∈ℂn×nA ^n× n and consider its Jordan canonical form =(+)−1A=P( +N)P^-1, where is diagonal and N is nilpotent, i.e., there exists k∈ℕk s.t. k=N^k=0. Then spectral conjugate ♯A is defined as ♯=¯−1,A =P P^-1, (63) where ¯ is the complex conjugate of . Interestingly enough, the spectral conjugate ♯A is different from the Hermitian conjugate H=T¯A^H= A^T, which is defined as the complex conjugate of the transpose of A. The spectral conjugate ♯A preserves the eigenvectors and generalized eigenvectors of A, while only taking the complex conjugate of its eigenvalues. In contrast, the Hermitian conjugate HA^H does not preserve the eigenvectors of A in general. Satisfying properties. Let τ:ℂ→ℂτ:C be the map τ(z)=z¯τ(z)= z, and functional extension to ∈ℂn×nA ^n× n as τ()=♯τ(A)=A . Then the spectral conjugate satisfies the following properties: • τ(+)=τ()+τ()τ(A+B)=τ(A)+τ(B), if A and B commute, • τ()B=τ()τ(A)B=B\,τ(A), if A and B commute, • τ(α)=α¯τ()τ(αA)= α\,τ(A), for any α∈ℂα , • τ(−1)=τ()−1τ(PAP^-1)=P\,τ(A)\,P^-1, for any invertible matrix P. Computation. In practice ♯A is computed from the complex Schur decomposition =HA=QTQ^H, with T upper triangular. The diagonal is set by conjugating the eigenvalues, ii=ii¯F_i= T_i, and the strictly upper-triangular entries are filled in super-diagonal by super-diagonal using the Parlett recurrence, which enforces the commutation =FT=TF (Parlett, 1976). For clustered or repeated eigenvalues the scalar recurrence is ill-conditioned and the blocked variant is used instead, with each diagonal block set to μ¯ μI for the cluster mean μ and the off-diagonal blocks recovered from Sylvester equations (Davies and Higham, 2003; Higham, 2008). Recomposing, ♯=HA =QFQ^H, yields the spectral conjugate without forming the Jordan decomposition explicitly. C.2 Drazin inverse definition Definition C.2 (Drazin inverse). Let ∈ℂn×nA ^n× n and let k=ind()k=ind(A) denote its index, i.e. the smallest non-negative integer such that rank(k)=rank(k+1)rank(A^k)=rank(A^k+1), equivalently the size of the largest Jordan block associated with the eigenvalue 00. The Drazin inverse DA^D is the unique matrix satisfying DD=D,D=D,k+1D=k.A^DA\,A^D=A^D, A\,A^D=A^DA, A^k+1A^D=A^k. (64) The commutativity condition D=DAA^D=A^DA is what distinguishes the Drazin inverse from the Moore–Penrose pseudoinverse †A , which satisfies †=AA A=A but does not in general commute with A. Spectral characterisation. Writing the Jordan decomposition of A with its invertible and nilpotent parts separated, =()−1, invertible,k=,A=P pmatrixC&0\\ 0&N pmatrixP^-1, C invertible, N^k=0, (65) the Drazin inverse admits the closed form D=(−1)−1.A^D=P pmatrixC^-1&0\\ 0&0 pmatrixP^-1. (66) Equivalently, DA^D acts on the spectrum of A through λ⟼λ−1,λ≠0,0,λ=0,λ\; \; casesλ^-1,&λ≠ 0,\\[2.0pt] 0,&λ=0, cases (67) inverting the non-zero part of the spectrum while annihilating the nilpotent part. In particular DA^D shares the eigenvectors and generalized eigenvectors of A, which is precisely why the two operators commute. Computation. In practice DA^D is obtained from a Schur decomposition =HA=QTQ^H whose eigenvalues are reordered so that the non-zero ones occupy the leading diagonal block: inverting that block, zeroing the trailing nilpotent block, and recovering the coupling term from a Sylvester equation yields DA^D without ever forming the Jordan decomposition explicitly (Stanimirović et al., 2015; Campbell and Meyer, 2009). When combined with the computation of the spectral conjugate (Appendix C.1), the Schur decomposition is already available and hence only the pointwise scalar inversion and computation of coupling terms is required. Ultimately the computation of ∘,↑,↑(∘)DL ,L ,L (L )^D is dominated by the Schur decomposition. C.3 Lifting holomorphic function to matrix The evaluation of a holomorphic function h at a matrix L follows from equation 2 where 1z−λ 1z-λ is replaced by the resolvent (z−)−1(zI-L)^-1 (Auscher et al., 1997) h()=12πi∮Γh(z)(z−)−1z,h(L)= 12π i _ h(z)\,(zI-L)^-1\,dz, (68) where Γ encloses the spectrum σ()σ(L) and (z−)(zI-L) is invertible for all z∈Γz∈ . C.4 Wirtinger derivative and holomorphicity of polynomial and rational kernels Wirtinger derivative. Consider a complex number z=x+jyz=x+jy and its conjugate z¯=x−jy z=x-jy, the Wirtinger derivatives are defined as follows: ∂z=12(∂x−j∂y),∂z¯=12(∂x+j∂y). _z= 12( _x-j _y), _ z= 12( _x+j _y). (69) For notation simplicity, we write ∂:=∂z∂:= _z and ∂¯:=∂z¯ ∂:= _ z. In particular the derivatives satisfy ∂¯z=0 ∂z=0 and ∂z¯=0∂ z=0. Holomorphicity. Define a degree K polynomial kernel as: pK:σ()→ℂ,z↦∑k=0Kckzk,p_K:σ(L) , z _k=0^Kc_kz^k, (70) where ck∈ℂc_k are the coefficients. Define a [K1,K2][K_1,K_2] rational kernel as: rK1,K2:σ()→ℂ,z↦pK1(z)pK2(z).r_K_1,K_2:σ(L) , z p_K_1(z)p_K_2(z). (71) Since no instance of z¯ z is in the expressions of pKp_K and rK1,K2r_K_1,K_2, the Wirtinger derivatives ∂¯pK,∂¯rK1,K2 ∂p_K, ∂r_K_1,K_2 are null. Thereby pKp_K and rK1,K2r_K_1,K_2 are holomorphic. Consequently, as canonical and Chebyshev polynomials can be expressed in the form pKp_K, they are holomorphic. Similarly, Cayley polynomials Levie et al. (2018) are special case of the rational kernels, hence they also are holomorphic. C.5 Kernels gain for non-diagonalizable operator Upper bound on gain of heat kernel Consider the heat kernel h(λ)=e−tλh(λ)=e^-tλ, which is holomorphic since ∂¯h=0 ∂h=0, so that ∂kh=h(k)∂^kh=h^(k) and h(k)(λ)=(−t)ke−tλh^(k)(λ)=(-t)^ke^-tλ. Substituting into equation 1 for a Jordan block nJ_n of size m gives h(n)=e−tλn∑k=0m−1(−t)k!nk=e−tλn(1−t22!⋯(−t)m−1(m−1)!01−t⋯(−t)m−2(m−2)!001⋯(−t)m−3(m−3)!⋱000⋯1).h(J_n)=e^-t _n _k=0^m-1 (-t)^kk!\,N_n^k=e^-t _n pmatrix1&-t& t^22!&·s& (-t)^m-1(m-1)!\\ 0&1&-t&·s& (-t)^m-2(m-2)!\\ 0&0&1&·s& (-t)^m-3(m-3)!\\ & & & & \\ 0&0&0&·s&1 pmatrix. (72) Since ‖nk‖2=1\|N_n^k\|_2=1 (with ∥⋅∥2=σmax(⋅)\|·\|_2= _max(·) the operator norm) for k≤m−1k≤ m-1, the gain of the block obeys ‖h(n)‖2≤e−tℜ(λn)∑k=0m−1tkk!=e−tℜ(λn)O(tm−1), \|h(J_n) \|_2\;≤\;e^-t ( _n) _k=0^m-1 t^kk!\;=\;e^-t ( _n)\,O\! (t^m-1 ), (73) so defectiveness contributes at most a polynomial factor in t, dominated by the exponential whenever ℜ(λn)>0 ( _n)>0. ∎ Lower bound on the gain of eβ↑e^βL . In the same way as for the heat kernel, on a Jordan block nJ_n of size m, ↑L acts as jλI[n]+12nj _I[n]I+ 12N_n, so eβ↑|n=ejβλI[n]eβn/2e^βL |_J_n=e^jβ _I[n]\,e^βN_n/2, with inverse e−jβλI[n]e−βn/2e^-jβ _I[n]\,e^-βN_n/2. Since ‖nk‖2=1\|N_n^k\|_2=1 for k≤m−1k≤ m-1, its smallest singular value satisfies min‖^‖2=1∥eβ↑|n^∥2=∥e−βn/2∥2−1≥(∑k=0m−1(|β|/2)k!)−1≥(1+|β|2)−(m−1). _\| x\|_2=1 \|e^βL |_J_n x \|_2= \|e^-βN_n/2 \|_2^-1\;≥\; ( _k=0^m-1 (|β|/2)^kk! )^-1\;≥\; (1+ |β|2 )^-(m-1). (74) The gain is therefore lower bounded by a factor in β. For frequencies (eigenvalues) in Jordan blocks of size m=1m=1, the gain is unitary. For diagonalizable L (i.e. m=1m=1 for all blocks), the gain is unitary for all frequencies. ∎ C.6 Minimal polynomial This subsection recalls few propositions and corollaries of Chan et al. (2026a) applicable for diagonalizable L on the minimal polynomial of the dissipative (diffusion) and asymmetric (advection) operators. This determines the maximal layer depth, providing additional guidance for the design of the TwinS-GCN architecture. Proposition 3. Let NCN_C be the number of pairs of conjugate eigenvectors, the degree of the minimal polynomial of ∘L is deg(m∘(λ))=N−NCdeg(m_L (λ))=N-N_C. Corollary 2. There exists a set of coefficients ck∈ℂc_k such that the matrix polynomial m∘(∘)=(∘)N−NC+∑k=0N−NC−1ck(∘)k=0.m_L (L )=(L )^N-N_C+ _k=0^N-N_C-1c_k(L )^k=0. (75) Proposition 4. Let NZN_Z be the number of eigenvalues with zero imaginary part, then it holds that the degree of the minimal polynomial of ↑L is deg(m↑(λ))=N−NZ+1deg(m_L (λ))=N-N_Z+1. Corollary 3. There exists a set of coefficients ck∈ℂc_k such that the matrix polynomial m↑(↑)=(↑)N−NZ+1+∑k=0N−NZck(↑)k=0.m_L (L )=(L )^N-N_Z+1+ _k=0^N-N_Zc_k(L )^k=0. (76) C.7 Justification on having both sum and ratio Graph Filters In Chan et al. (2026a), the sum and ratio graph kernels are designed to perform different types of filtering. It is important to first note that the image of hsumh_sum and hrath_rat are not mutually inclusive, and span different subspaces. Additionally, authors of Chan et al. (2026a) have shown that hsumh_sum are better suited for diffusive low-pass filtering, while hrath_rat are better suited for advective low-pass filtering and phase manipulation Chan et al. (2025); Chan et al. (2026b). Finally they can be composed disregarding the order of composition, as the applications of both kernels to L commute. C.8 TwinS-GCN Combined architecture: TwinS (C) Following our sum and ratio graph filters proposed in equation 8, one can simply combine both graph filters as they span different subspaces. Denote (l)X^(l) the node features of dimension N×CN× C (nodes times channels) at the l-th layer. The (l+1)(l+1) layer update reads as: (l+1)=ψ(CLOSEν∑k=0K1Tk(~∘)(l)k(l,∘)+(1−ν)∑k=0K2jkTk(~↑)(l)k(l,↑)+ρ∑k=0KjkTk(~)(l)k(l)+(l)), splitX^(l+1)=ψ (&ν _k=0^K_1T_k( L )X^(l) _k^(l, )+(1-ν) _k=0^K_2j^k\,T_k( L )X^(l) _k^(l, )\\ &+ρ _k=0^Kj^k\,T_k( Z)X^(l) _k^(l)+B^(l) ), split (77) where ψ is a nonlinear function, k(l,∘),k(l,↑) _k^(l, ), _k^(l, ) and k(l) _k^(l) are respectively learnable real weight matrices, (l)B^(l) is a learnable bias term and learnable parameters ν,ρ∈[0,1]ν,ρ∈[0,1]. C.9 Conceptual comparison with advection diffusion graph neural networks As discussed in Section 4, for diagonalizable shifts, equation 5 coincides with the graph diffusion and advection operators of Chan et al. (2026a) (separate discussion in following paragraph). Hence, we compare conceptually the graph diffusion and advection operators obtained through our construct for diagonalizable operators with three related works. Berlureau et al. (2026) reweight the undirected Laplacian by a learned potential field, which remains symmetric and therefore anisotropically diffusive. ADR-GNN (Eliasof et al., 2024) learns the edges of its advection operator, following Chapman (2015), whose construction conserves mass rather than emulating transport and is neither skew-symmetric nor read off the input graph. AdvDIFFormer (Wu et al., 2023) models diffusion by global attention and advection by message passing along directed edges. In each case the directional term is an architectural addition to a dissipative model and carries no purely imaginary spectrum. While similar to the work from Chan et al. (2026a), crucial differences include that our framework applies to non-diagonalizable operators, provide foundations in non-holomorphic calculus, and address problems of oversmoothing and receptive field reach. It is also in practice more stable to use thanks to our contribution in its Chebyshev reparameterization. Finally Chan et al. (2026a) belongs to directed GSP literature and does not tackle topics of neural networks. C.10 On polynomial interpolation of the spectral conjugate On a finite and fixed graph, any spectral window is exactly realizable by a polynomial in the shift operator. If L has distinct eigenvalues λ1,…,λK _1,…, _K and minimal polynomial of degree m, then h()=p()h(L)=p_L(L) for the unique interpolant of degree ≤m−1≤ m-1 matching the data ∂mh(λj)∂^mh( _j) (Nevanlinna, 2018, Prop. 5.1, 5.3). Non-holomorphicity of h supposedly could be overcome. We show by explicit computation that the resulting coefficients are nonetheless unusable, whereas for holomorphic they remain stable. In both experiments below we fit a quadratic p(z)=a0+a1z+a2z2p(z)=a_0+a_1z+a_2z^2 through three eigenvalues by solving the Vandermonde system, and compare two target kernels: the holomorphic h(z)=z2+1h(z)=z^2+1 and the conjugation map τ(z)=z¯τ(z)= z, for which ∂¯τ=1≠0 ∂τ=1≠ 0. We write z0=1+iz_0=1+i and let ε>0 >0. Experiment 1: two close eigenvalues. Take the eigenvalues λ1=0 _1=0, λ2=z0 _2=z_0 and λ3=z0+εd _3=z_0+ d with |d|=1|d|=1, so that λ3 _3 lies at distance ε from λ2 _2 in the direction d. As ε→0 → 0 the node set converges to 0,z0,z0\0,z_0,z_0\ for every choice of d. For the holomorphic target kernel the interpolant is exact and independent of both ε and d, p(z)=z2+1,(a0,a1,a2)=(1,0,1).p(z)=z^2+1, (a_0,a_1,a_2)=(1,0,1). (78) For τ one obtains a0=0a_0=0 and d=1:a1=i(−ε−3−i)ε+1+i,a2=1+iε+1+i; d=1: a_1= i(- -3-i) +1+i, a_2= 1+i +1+i; (79) d=i:a1=i(−ε−3+i)ε+1−i,a2=1+iε+1−i. d=i: a_1= i(- -3+i) +1-i, a_2= 1+i +1-i. (80) These coefficients are bounded for small ε , but the denominators differ depending on the direction d limε→0(a1,a2)=(−1−2i, 1),d=1,( 1−2i,i),d=i. _ → 0(a_1,a_2)= cases(-1-2i,\;1),&d=1,\\ (\;\;1-2i,\;i),&d=i. cases (81) From this example, we see that small numerical inaccuracies from eigensolvers may lead to two completely different filters, caused by the direction of approach for which, smaller distance is ineffective. Experiment 2: three close eigenvalues. Take the eigenvalues λ1=z0 _1=z_0, λ2=z0+ε _2=z_0+ , λ3=z0+iε _3=z_0+i , so all three nodes lie within ε of z0z_0. Again the holomorphic target gives (a0,a1,a2)=(1,0,1)(a_0,a_1,a_2)=(1,0,1) for every ε , while for τ a0(ε)=2(−1+i)ε,a1(ε)=−i(ε+4)ε,a2(ε)=1+iε,a_0( )= 2(-1+i) , a_1( )=- i( +4) , a_2( )= 1+i , (82) so that maxk|ak|=(ε+4)/ε∼4ε−1 _k|a_k|=( +4)/ 4 ^-1. ε=10−1 =10^-1 10−210^-2 10−310^-3 10−410^-4 maxk|ak| _k|a_k|, target z2+1z^2+1 1.01.0 1.01.0 1.01.0 1.01.0 maxk|ak| _k|a_k|, target z¯ z 4.1×1014.1× 10^1 4.0×1024.0× 10^2 4.0×1034.0× 10^3 4.0×1044.0× 10^4 Table 3: Coefficient magnitude under a cluster of three. Beyond the non-uniquely determined interpolation from Experiment.1, the interpolating coefficients may also suffer from arbitrarily large values, making it numerically unstable. While we show ill-defined interpolation and numerical instability for the case of τ (spectral conjugate), it surmised to be extended to general non-holomorphic kernels. Either of these examples are realistic especially for non-diagonalizable operators where there necessarily exist an eigenvalue with algebraic multiplicity greater than 11. C.11 Combine with anti-symmetric weight A-DGN (Gravina et al., 2023) models each node feature n(t)∈ℝCx_n(t) ^C as the state of an ODE whose channel-mixing matrix is constrained to be anti-symmetric (their equation 4): ∂tn(t)=σ((−⊤)n(t)+Φ((t),n)+), _tx_n(t)=σ ((V-V )\,x_n(t)+ (X(t),N_n)+b ), (83) where ∈ℝC×CV ^C× C is learnable and Φ aggregates over the neighborhood nN_n. Anti-symmetry places the eigenvalues of −⊤V-V on the imaginary axis, which is the source of A-DGN’s non-dissipative behavior. A forward-Euler discretization with step ϵε and a small damping γ>0γ>0 for numerical stability (their equation 5) gives, in matrix form, (l+1)=(l)+ϵψ((l)(−⊤−γ)+Φ((l))+).X^(l+1)=X^(l)+ε\,ψ (X^(l)(V-V -γI)+ (X^(l))+B ). (84) The two constructions act on different axes: A-DGN constrains how channels are mixed, whereas TwinS (R) constrains how information moves across nodes. They therefore combine by using the ratio filter of equation 15 as the aggregation Φ : (l+1)=(l)+ϵψ((l)(−⊤−γ)+∑k=0KjkTk(~)(l)k(l)+(l)),X^(l+1)=X^(l)+ε\,ψ (X^(l)(V-V -γI)+ _k=0^Kj^k\,T_k( Z)\,X^(l) _k^(l)+B^(l) ), (85) where V is the anti-symmetric channel-mixing weight and k(l) _k^(l) are the filter weights of TwinS (R). TwinS (S) and TwinS (C) combine in the same way, by substituting equation 14 or equation 77 for Φ . Whether the non-dissipativity guarantees of both components hold jointly is not immediate. We leave this analysis, and the empirical evaluation of the combination, to future work. Appendix D Experimental details D.1 Graph transfer task General settings. The task requires the model to learn a mapping from a source node to a target node on a graph of size N. Denote [n]∈ℝCx[n] ^C to be the feature vector at node n. Following Bodnar et al. (2021) the node features are set as follows: [n]=k,if n is source node,,if n is target node,+ϵ,otherwise.x[n]= casese_k,&if n is source node,\\ 0,&if n is target node,\\ 1+ ε,&otherwise. cases (86) where ke_k is set to be a one-hot vector such that k[n]=δk=n e_k[n]= _k=n with integer k∈1,…,Ck∈\1,…,C\ being randomly selected and ϵ∼(,0.1N) ε (0,0.1I_N). The source node and target node being selected such that the distance between them is d(nsource,ntarget)=rd(n_source,n_target)=r. The ground truth for the target node’s feature is set to be the same as the source node’s. The model is evaluated through mean squared error (MSE) between the predicted target node feature and the ground truth target node feature. The predicted class is the index with highest value (since all others should be receiving less energy) and it is compared with the ground truth class k, leading to an accuracy for each sample. The graphs used for the task are shown in Figure 5. Without loss of generality, in our experiments we select node feature of dimension C=5C=5, consider 200200 samples, split into (60%) training, (20%) validation and (20%) testing. The source node is randomly selected, while the target node is selected such that the distance between them is r=10r=10 which is the graph diameter. Figure 5: Graphs used for the task (same as in Gravina et al. (2025)). Details on Section 5.1. For depth experiment: source and target are held a fixed distance of 1010 apart, and the sweep varies only network depth. For every topology, model and depth ∈1,…,10K∈\1,…,10\ we run an independent hyperparameter search (Appendix D.1), so no depth is handicapped by settings tuned at another. Each configuration is trained on 200200 sampled graphs split 60%/20%/20%60\%/20\%/20\% into 120120 training, 4040 validation and 4040 test samples under a fixed data seed, with the one-hot label placed on the source node and all other nodes carrying an uninformative constant perturbed by Gaussian noise (σ=0.1σ=0.1). Training minimizes the mean squared error between the target node’s output and the one-hot label, plus a suppression term penalizing the squared output of every non-target node, which prevents the network from spreading the label rather than transporting it. We report accuracy: a prediction counts as correct when its largest coordinate matches the label’s. Runs are capped at 300300 epochs and evaluated every 88 epochs, stopping early after 3030 consecutive evaluations (i.e. 240240 epochs) without an improvement in validation MSE. Each configuration is scored at its lowest-validation-MSE epoch and the test accuracy is read at that same epoch, so the test split never informs selection. Every reported figure is the mean and standard deviation over 1010 random initializations. Hyperparameter search. The search uses Optuna’s Tree-structured Parzen Estimator (Akiba et al., 2019), run independently for each (topology, model, depth) triple at 4848 trials, and selects on validation MSE. Depth is swept exhaustively over 1,…,10\1,…,10\ rather than sampled, so each model is tuned 480480 times per topology. The searched ranges are: hidden dimension ∈4,8,16,32∈\4,8,16,32\, non-linearity ∈ReLU,None∈\ReLU,None\, normalization ∈True,False∈\True,False\, bias ∈True,False∈\True,False\, and learning rate ∈0.001,0.005∈\0.001,0.005\. Dropout and weight decay are held at 00 throughout: on a task this small, regularization only slows transport without improving generalization. Three further axes are searched only by the architectures that read them — whether the branch-mixing weights α and β are learned or fixed at 0.50.5 (employed in Dir-GNN, FaberNet, TwinS-GCN), and MagNet’s charge q∈0,0.05,0.10,0.15,0.20,0.25q∈\0,0.05,0.10,0.15,0.20,0.25\. While introduced as the standard combinatorial directed Laplacian, options for operator include left,combinatorial,right,symmetric\left,combinatorial,right,symmetric\ normalized directed Laplacian. Cluster tolerance (Appendix C.1) was set to 1e−91e^-9, as we expect no distinct eigenvalue clusters. A dash in the table (Table 4) means that the hyperparameter is not relevant to the current architecture. Accuracies are the mean and standard deviation over random initializations. For the same model, when several depths achieve perfect accuracy, the table reports the shallowest, since on a task that many models saturate the meaningful winner is the cheapest one. Corresponding run-times for each of these models are shown in Table 7. For Dir-GNN, FaberNet, MagNet TwinS-GCN, we consider only the degree of each layer to be 11 (to fully isolate the effect of depth), meaning that we set K=1K=1 for Dir-GNN, FaberNet, MagNet, while for TwinS-GCN, which start at index K=1K=1 (programming convention) for the identity matrix, we set K=2K=2 to represent degree 11 polynomial. Table 4: Best hyperparameter configuration per topology for the graph transfer task, selected over depth and all searched axes by validation MSE. A dash marks an axis the architecture does not take. α/β report whether the branch-mixing weights are learned; q is MagNet’s charge. Model depth #hidden non-lin. norm. bias lr α β q Accuracy Directed Cycle GCN 11 88 ReLU No No 0.0050.005 – – – 0.325±0.0540.325± 0.054 MagNet 66 44 – – Yes 0.0050.005 – – 0.050.05 0.267±0.0310.267± 0.031 FaberNet 1010 3232 – – Yes 0.0050.005 Yes – – 1.000±0.0001.000± 0.000 Dir-GNN 1010 1616 None Yes No 0.0050.005 Yes – – 1.000±0.0001.000± 0.000 ADR-GNN 22 44 – – – 0.0050.005 – – – 0.283±0.0240.283± 0.024 AdvDIFFormer 11 1616 – – – 0.0050.005 – – – 1.000±0.0001.000± 0.000 TwinS (S) 22 44 None Yes No 0.0010.001 Yes – – 0.258±0.0420.258± 0.042 TwinS (R) 33 88 None Yes No 0.0010.001 – – – 1.000±0.0001.000± 0.000 TwinS (C) 33 1616 None Yes Yes 0.0050.005 No Yes – 1.000±0.0001.000± 0.000 Directed Crossed Cycle GCN 1010 3232 ReLU No Yes 0.0050.005 – – – 0.467±0.0660.467± 0.066 MagNet 33 88 – – No 0.0010.001 – – 0.10.1 0.233±0.0120.233± 0.012 FaberNet 1010 1616 – – Yes 0.0050.005 Yes – – 1.000±0.0001.000± 0.000 Dir-GNN 1010 3232 ReLU No No 0.0050.005 Yes – – 1.000±0.0001.000± 0.000 ADR-GNN 66 88 – – – 0.0050.005 – – – 0.300±0.0200.300± 0.020 AdvDIFFormer 11 3232 – – – 0.0050.005 – – – 1.000±0.0001.000± 0.000 TwinS (S) 22 3232 ReLU Yes Yes 0.0050.005 Yes – – 1.000±0.0001.000± 0.000 TwinS (R) 22 1616 None Yes No 0.0010.001 – – – 1.000±0.0001.000± 0.000 TwinS (C) 22 3232 None Yes Yes 0.0010.001 No No – 1.000±0.0001.000± 0.000 Directed Clique Path GCN 33 1616 None Yes Yes 0.0010.001 – – – 0.275±0.0200.275± 0.020 MagNet 33 44 – – Yes 0.0050.005 – – 0.00.0 0.250±0.0200.250± 0.020 FaberNet 1010 1616 – – No 0.0010.001 Yes – – 0.258±0.0310.258± 0.031 Dir-GNN 1010 1616 None Yes No 0.0050.005 Yes – – 0.767±0.1120.767± 0.112 ADR-GNN 1010 88 – – – 0.0050.005 – – – 0.275±0.0200.275± 0.020 AdvDIFFormer 11 3232 – – – 0.0050.005 – – – 1.000±0.0001.000± 0.000 TwinS (S) 22 3232 None Yes No 0.0050.005 Yes – – 1.000±0.0001.000± 0.000 TwinS (R) 22 1616 ReLU Yes Yes 0.0050.005 – – – 1.000±0.0001.000± 0.000 TwinS (C) 22 1616 None Yes Yes 0.0010.001 Yes Yes – 1.000±0.0001.000± 0.000 D.2 Datasets benchmarking General settings. The benchmarking datasets used in the experiments are summarized in Table 5. We evaluate on transductive node classification over directed graphs. Cornell, Texas and Wisconsin are small web-page networks that are strongly heterophilic (edge homophily 0.12, 0.06 and 0.17 respectively), so neighboring nodes usually carry different labels and edge direction is informative. We use the provided ten train/validation/test splits distributed with the datasets. Every configuration is trained on the training nodes for at most 300 epochs, evaluated every 5 epochs, with early stopping after 20 evaluations without validation improvement. Following standard practice, we score each run at the epoch of highest validation accuracy and report the test accuracy at that same epoch, so the test nodes never influence model selection. Each reported figure is the mean and standard deviation of classification accuracy over 30 runs — the ten splits crossed with three random initializations. Table 5: Statistics of the node classification benchmarks. Dataset #Nodes #Edges #Feat. #Cls. Cornell 183 295 1,703 5 Texas 183 309 1,703 5 Wisconsin 251 499 1,703 5 Chameleon† 2,277 36,101 2,325 5 Squirrel† 5,201 217,073 2,089 5 † Original versions. Filtered variants (Platonov et al., 2023): Chameleon (890 / 13,584), Squirrel (2,223 / 46,998). Among the results provided in Table 1 here are the ones retrieved from literature: • Results of GCN, GAT on Wisconsin, Cornell, Texas are taken from Eliasof et al. (2024) and on Chameleon-F, Squirrel-F, comes from Platonov et al. (2023). • Results of MLP on Chameleon-F, Squirrel-F, comes from Platonov et al. (2023) (ResNet). All remaining entries are computed from re-implementation following original repository. Disclaimer: the results reproduced for MLP, MagNet, ADR-GNN on Wisconsin, Cornell, Texas may differ from accuracies found in their respective original paper Zhang et al. (2021), Eliasof et al. (2024). We have noticed in their open repository that the employed Wisconsin, Cornell, Texas datasets were from the pre-2022 geom-gcn (Pei et al., 2020) release while in our experiments we employ the latest release (Pei et al., 2022), hence we reproduce the results on these datasets. Hyperparameter search. Selected hyperparameters for the benchmarking datasets are reported in Table 6. The search uses Optuna’s TPE sampler (multivariate), and configurations are selected on mean validation accuracy. The shared search space is: depth ∈2,3,4,5∈\2,3,4,5\, hidden dimension ∈32,64,128,256∈\32,64,128,256\, dropout ∈0,0.2,0.4,0.5,0.6∈\0,0.2,0.4,0.5,0.6\, normalization and bias ∈True,False∈\True,False\, learning rate ∈0.001,0.005,0.01,0.05∈\0.001,0.005,0.01,0.05\, weight decay ∈0,5.10−4,5.10−3∈\0,5.10^-4,5.10^-3\, input feature scaling ∈raw,row,std∈\raw,row,std\, jumping knowledge ∈None,max,cat∈\None,max,cat\ and, for spectral models, order K∈1,2,3,4K∈\1,2,3,4\. The non-linearity is fixed to ReLU. Each model also searches its own parameters: MagNet’s charge q∈0,0.05,…,0.25q∈\0,0.05,…,0.25\; the Dir-GNN and FaberNet mixing weight α∈0,0.5,1α∈\0,0.5,1\ and whether it is learned; FaberNet’s zero-order term; the ADR-GNN step size ∈0.02,0.05,0.1,0.2∈\0.02,0.05,0.1,0.2\; and the AdvDIFFormer heads ∈1,2,4∈\1,2,4\ and order Kord∈1,2,3K_ord∈\1,2,3\. For TwinS-GCN, the search covers the branch-mixing weights α,βα,β, the Laplacian normalization ∈left,combinatorial,right,symmetric∈\left,combinatorial,right,symmetric\, and the eigenvalue-clustering tolerance ∈default,10−6,10−3∈\default,10^-6,10^-3\. Operators are rescaled by their maximum absolute entry. Next we discuss few differences in protocol: on WebKB, TwinS-GCN uses 7272 trials, against 9696 for MLP, MagNet, Dir-GNN and FaberNet and 7272 for ADR-GNN and AdvDIFFormer to find hyperparameters. On Chameleon-F and Squirrel-F, TwinS-GCN uses 3232 trials, against 248248 and 4848 respectively. TwinS-GCN thus has the smallest budget on every dataset. Its validation score, however, averages over all 1010 splits, whereas the baselines average over 33 (55 for the attention-based baselines on Chameleon-F and Squirrel-F). Finally, for Texas and Squirrel-F the TwinS-GCN search was widened to depth 11, K∈5,…,8K∈\5,…,8\, weight decay 10−210^-2 and patience 5050, which accounts for the settings in Table 6 that are not within the shared search space. Table 6: Hyperparameter search results for benchmarking datasets. A dash marks an axis the architecture does not take. α/β report whether the branch-mixing weights are learned; q is MagNet’s charge. Cluster tol are 1e−91e^-9 unless specified. Model depth #hidden dropout norm. bias lr w.d. feat. JK K learn α learn β model-specific Accuracy Chameleon-F MagNet 33 128128 0.20.2 – Yes 0.0050.005 0.0050.005 raw – 11 – – q=0q=0 0.424±0.0350.424± 0.035 Dir-GNN 55 128128 0.50.5 No Yes 0.0050.005 0.00050.0005 row cat – No – α=0.5α=0.5 0.441±0.0290.441± 0.029 FaberNet 44 128128 0.60.6 – Yes 0.010.01 0.00050.0005 row max 11 No – α=0.5α=0.5, zero-ord.: No 0.430±0.0330.430± 0.033 AdvDIFFormer 44 128128 0.40.4 – – 0.050.05 00 std – – – – heads 22, Kord=1K_ord=1 0.415±0.0360.415± 0.036 ADR-GNN 44 6464 0.60.6 – – 0.010.01 0.0050.005 std – – – – step 0.050.05 0.461±0.0410.461± 0.041 TwinS (S) 22 256256 0.60.6 Yes No 0.0010.001 00 std cat 22 Yes – right 0.423±0.0280.423± 0.028 TwinS (R) 22 256256 0.50.5 No Yes 0.0010.001 0.0050.005 std cat 44 – – none, tol 10−610^-6 0.442±0.0320.442± 0.032 TwinS (C) 22 128128 0.40.4 Yes No 0.010.01 00 std none 44 No Yes none, tol 10−310^-3 0.416±0.0360.416± 0.036 Squirrel-F MagNet 55 128128 0.20.2 – Yes 0.050.05 00 raw – 22 – – q=0.05q=0.05 0.442±0.0140.442± 0.014 Dir-GNN 44 6464 0.50.5 No Yes 0.0050.005 0.00050.0005 row cat – No – α=0.5α=0.5 0.443±0.0220.443± 0.022 FaberNet 33 6464 0.60.6 – Yes 0.0010.001 0.0050.005 raw None 11 No – α=0.5α=0.5, zero-ord.: No 0.451±0.0210.451± 0.021 AdvDIFFormer 22 6464 0.60.6 – – 0.0050.005 00 std – – – – heads 11, Kord=1K_ord=1 0.414±0.0160.414± 0.016 ADR-GNN 55 6464 0.40.4 – – 0.050.05 0.0050.005 std – – – – step 0.20.2 0.441±0.0200.441± 0.020 TwinS (S) 33 6464 0.60.6 Yes No 0.0010.001 0.00050.0005 std cat 44 Yes – left, tol 10−310^-3 0.396±0.0150.396± 0.015 TwinS (R) 22 3232 0.60.6 Yes No 0.0050.005 0.0050.005 std None 22 – – left, tol 10−310^-3 0.402±0.0180.402± 0.018 TwinS (C) 22 128128 0.60.6 Yes Yes 0.0050.005 00 std none 55 Yes Yes left 0.401±0.0160.401± 0.016 Texas MLP 44 3232 0.40.4 – Yes 0.010.01 0.0050.005 raw – – – – – 0.822±0.0430.822± 0.043 MagNet 22 256256 0.20.2 – Yes 0.010.01 0.0050.005 raw – 11 – – q=0.1q=0.1 0.853±0.0460.853± 0.046 Dir-GNN 22 6464 0.60.6 Yes No 0.050.05 0.0050.005 raw cat – Yes – α=0α=0 0.705±0.0650.705± 0.065 FaberNet 22 6464 0.50.5 – Yes 0.0050.005 0.0050.005 raw cat 22 No – α=1α=1, zero-ord.: Yes 0.832±0.0530.832± 0.053 AdvDIFFormer 33 3232 0.20.2 – – 0.050.05 0.00050.0005 raw – – – – heads 44, Kord=1K_ord=1 0.818±0.0640.818± 0.064 ADR-GNN 22 6464 0.20.2 – – 0.050.05 0.0050.005 raw – – – – step 0.020.02 0.822±0.0520.822± 0.052 TwinS (S) 22 6464 0.50.5 Yes Yes 0.010.01 0.00050.0005 raw max 33 No – left 0.846±0.0340.846± 0.034 TwinS (R) 11 128128 0.50.5 Yes No 0.010.01 0.00050.0005 row cat 88 – – comb. 0.841±0.0390.841± 0.039 TwinS (C) 22 256256 0.20.2 Yes Yes 0.010.01 0.010.01 raw cat 88 Yes No left, tol 10−310^-3 0.854±0.0470.854± 0.047 Wisconsin MLP 44 6464 0.40.4 – Yes 0.010.01 0.0050.005 raw – – – – – 0.852±0.0500.852± 0.050 MagNet 22 3232 0.40.4 – Yes 0.010.01 0.0050.005 raw – 11 – – q=0.15q=0.15 0.834±0.0400.834± 0.040 Dir-GNN 55 6464 0.40.4 Yes No 0.050.05 0.00050.0005 raw max – No – α=0.5α=0.5 0.674±0.0600.674± 0.060 FaberNet 22 256256 0.50.5 – Yes 0.0010.001 0.0050.005 raw cat 22 No – α=1α=1, zero-ord.: Yes 0.795±0.0520.795± 0.052 AdvDIFFormer 22 256256 0.40.4 – – 0.010.01 0.00050.0005 raw – – – – heads 11, Kord=2K_ord=2 0.765±0.0480.765± 0.048 ADR-GNN 33 6464 0.20.2 – – 0.050.05 0.0050.005 raw – – – – step 0.020.02 0.845±0.0450.845± 0.045 TwinS (S) 22 6464 0.20.2 Yes No 0.0050.005 0.00050.0005 row cat 22 Yes – right 0.871±0.0410.871± 0.041 TwinS (R) 22 256256 0.20.2 Yes Yes 0.0050.005 0.0050.005 raw cat 33 – – right 0.874±0.0400.874± 0.040 TwinS (C) 22 6464 00 Yes No 0.010.01 0.0050.005 raw max 22 Yes No right 0.868±0.0380.868± 0.038 Cornell MLP 33 6464 0.20.2 – Yes 0.010.01 0.0050.005 raw – – – – – 0.759±0.0390.759± 0.039 MagNet 22 3232 0.20.2 – Yes 0.010.01 0.0050.005 raw – 11 – – q=0.25q=0.25 0.733±0.0420.733± 0.042 Dir-GNN 22 128128 0.20.2 Yes No 0.050.05 0.0050.005 std None – No – α=0.5α=0.5 0.598±0.0690.598± 0.069 FaberNet 33 3232 0.20.2 – No 0.010.01 0.0050.005 raw cat 44 No – α=0α=0, zero-ord.: Yes 0.726±0.0470.726± 0.047 AdvDIFFormer 22 128128 0.40.4 – – 0.010.01 0.00050.0005 raw – – – – heads 22, Kord=1K_ord=1 0.722±0.0550.722± 0.055 ADR-GNN 22 128128 00 – – 0.0050.005 0.0050.005 raw – – – – step 0.020.02 0.714±0.0440.714± 0.044 TwinS (S) 22 6464 0.50.5 Yes Yes 0.010.01 0.00050.0005 raw max 33 No – sym. 0.751±0.0320.751± 0.032 TwinS (R) 22 6464 0.50.5 Yes Yes 0.010.01 0.00050.0005 raw max 33 – – right 0.768±0.0300.768± 0.030 TwinS (C) 22 6464 0.50.5 Yes Yes 0.010.01 0.00050.0005 raw max 33 No No right 0.757±0.0400.757± 0.040 Appendix E Models Runtime Graph transfer task. We report the runtime of the models used in the graph transfer task. The runtime are provided for all different graphs. The runtime are measured on Apple M1 Max. Table 7 shows the minimal needed runtime/parameters of each model to solve the graph transfer task. Table 8 shows for comparable configuration, the models runtime/parameters. Table 7: Runtime and parameter count per topology for the graph transfer task as specified in Appendix D.1, for the best configuration selected in Table 4. Train time is per-epoch, averaged over 100100 epochs (train repeats is 55) and inference time (infer repeats is 5050) is per forward pass over the full graph (4040 graphs/batch). Build time is provided specifically to TwinS-GCN models which require a one-time precomputing, per graph, of the operators. Accuracy is reported for reference. Model #Params Train (ms/epoch) Inference (ms) Build (ms) Accuracy Directed Cycle GCN 2525 0.67±0.340.67± 0.34 0.09±0.000.09± 0.00 – 0.325±0.0540.325± 0.054 MagNet 269269 9.12±0.289.12± 0.28 1.80±0.011.80± 0.01 – 0.267±0.0310.267± 0.031 FaberNet 39,11039,110 26.94±0.7826.94± 0.78 4.44±0.004.44± 0.00 – 1.000±0.0001.000± 0.000 Dir-GNN 4,4174,417 10.58±0.0610.58± 0.06 1.96±0.001.96± 0.00 – 1.000±0.0001.000± 0.000 ADR-GNN 345345 11.67±0.5211.67± 0.52 1.61±0.011.61± 0.01 – 0.283±0.0240.283± 0.024 AdvDIFFormer 1,3171,317 3.66±0.233.66± 0.23 0.48±0.140.48± 0.14 – 1.000±0.0001.000± 0.000 TwinS (S) 161161 2.88±0.432.88± 0.43 0.43±0.040.43± 0.04 0.00.0 0.258±0.0420.258± 0.042 TwinS (R) 288288 2.60±0.142.60± 0.14 0.34±0.000.34± 0.00 0.00.0 1.000±0.0001.000± 0.000 TwinS (C) 2,7192,719 6.40±0.716.40± 0.71 1.17±0.071.17± 0.07 0.00.0 1.000±0.0001.000± 0.000 Directed Crossed Cycle GCN 8,8058,805 14.49±0.5414.49± 0.54 2.89±0.162.89± 0.16 – 0.467±0.0660.467± 0.066 MagNet 421421 4.92±0.094.92± 0.09 0.98±0.000.98± 0.00 – 0.233±0.0120.233± 0.012 FaberNet 10,34210,342 23.55±0.1723.55± 0.17 4.39±0.144.39± 0.14 – 1.000±0.0001.000± 0.000 Dir-GNN 17,02517,025 10.47±0.2610.47± 0.26 1.80±0.141.80± 0.14 – 1.000±0.0001.000± 0.000 ADR-GNN 3,2133,213 53.38±0.6953.38± 0.69 9.21±0.209.21± 0.20 – 0.300±0.0200.300± 0.020 AdvDIFFormer 4,6774,677 4.45±0.084.45± 0.08 0.63±0.000.63± 0.00 – 1.000±0.0001.000± 0.000 TwinS (S) 1,4291,429 2.98±0.012.98± 0.01 0.52±0.060.52± 0.06 0.50.5 1.000±0.0001.000± 0.000 TwinS (R) 320320 1.73±0.041.73± 0.04 0.23±0.000.23± 0.00 0.90.9 1.000±0.0001.000± 0.000 TwinS (C) 2,1422,142 3.88±0.053.88± 0.05 0.79±0.060.79± 0.06 0.90.9 1.000±0.0001.000± 0.000 Directed Clique Path GCN 453453 3.79±0.093.79± 0.09 0.69±0.000.69± 0.00 – 0.275±0.0200.275± 0.020 MagNet 161161 4.84±0.184.84± 0.18 0.96±0.000.96± 0.00 – 0.250±0.0200.250± 0.020 FaberNet 9,7029,702 21.28±0.1521.28± 0.15 4.01±0.004.01± 0.00 – 0.258±0.0310.258± 0.031 Dir-GNN 4,4174,417 10.55±0.0810.55± 0.08 1.95±0.001.95± 0.00 – 0.767±0.1120.767± 0.112 ADR-GNN 5,2935,293 113.09±0.73113.09± 0.73 21.52±0.8521.52± 0.85 – 0.275±0.0200.275± 0.020 AdvDIFFormer 4,6774,677 4.92±0.114.92± 0.11 0.65±0.010.65± 0.01 – 1.000±0.0001.000± 0.000 TwinS (S) 1,2811,281 2.81±0.042.81± 0.04 0.52±0.060.52± 0.06 0.20.2 1.000±0.0001.000± 0.000 TwinS (R) 362362 1.71±0.031.71± 0.03 0.27±0.030.27± 0.03 0.40.4 1.000±0.0001.000± 0.000 TwinS (C) 1,0881,088 4.30±0.044.30± 0.04 0.78±0.060.78± 0.06 0.40.4 1.000±0.0001.000± 0.000 Table 8: Runtime and parameter count per topology for the graph transfer task as specified in Appendix D.1. The count is under an equivalent configuration (hidden=32, depth=4, polynomial degree (K)=1, dropout=0, bias=True, non-linearity=ReLU). Remaining information is as in Table 7. Model #Params Train (ms/epoch) Inference (ms) Directed Cycle GCN 2,4692,469 8.93±1.808.93± 1.80 1.02±0.241.02± 0.24 MagNet 6,9176,917 17.34±1.3017.34± 1.30 3.01±0.683.01± 0.68 FaberNet 27,20527,205 38.95±0.3438.95± 0.34 6.77±0.606.77± 0.60 Dir-GNN 4,9384,938 9.33±0.219.33± 0.21 1.40±0.121.40± 0.12 ADR-GNN 30,18130,181 134.74±3.63134.74± 3.63 22.85±4.6322.85± 4.63 AdvDIFFormer 17,44517,445 34.32±0.3934.32± 0.39 5.04±1.115.04± 1.11 TwinS (S) 9,8769,876 11.96±0.2511.96± 0.25 2.05±0.722.05± 0.72 TwinS (R) 4,9384,938 10.03±4.0010.03± 4.00 1.61±0.971.61± 0.97 TwinS (C) 14,81414,814 22.52±2.2122.52± 2.21 3.22±1.563.22± 1.56 Directed Crossed Cycle GCN 2,4692,469 12.64±0.4312.64± 0.43 2.20±0.892.20± 0.89 MagNet 6,9176,917 19.03±1.8519.03± 1.85 4.64±1.574.64± 1.57 FaberNet 27,20527,205 49.51±3.2449.51± 3.24 7.91±1.937.91± 1.93 Dir-GNN 4,9384,938 10.00±0.3310.00± 0.33 1.79±0.811.79± 0.81 ADR-GNN 30,18130,181 193.44±16.22193.44± 16.22 35.31±7.3835.31± 7.38 AdvDIFFormer 17,44517,445 41.97±2.8441.97± 2.84 5.76±2.455.76± 2.45 TwinS (S) 9,8769,876 13.75±0.8313.75± 0.83 2.21±0.372.21± 0.37 TwinS (R) 4,9384,938 8.47±0.478.47± 0.47 1.11±0.341.11± 0.34 TwinS (C) 14,81414,814 20.97±1.7420.97± 1.74 3.28±0.913.28± 0.91 Directed Clique Path GCN 2,4692,469 17.46±0.9917.46± 0.99 2.81±1.182.81± 1.18 MagNet 6,9176,917 20.83±3.7420.83± 3.74 4.27±1.244.27± 1.24 FaberNet 27,20527,205 47.63±6.6147.63± 6.61 7.50±1.427.50± 1.42 Dir-GNN 4,9384,938 9.92±0.809.92± 0.80 1.56±0.411.56± 0.41 ADR-GNN 30,18130,181 303.24±15.23303.24± 15.23 53.27±23.5053.27± 23.50 AdvDIFFormer 17,44517,445 33.62±1.1633.62± 1.16 4.62±1.474.62± 1.47 TwinS (S) 9,8769,876 14.47±0.3614.47± 0.36 2.58±0.952.58± 0.95 TwinS (R) 4,9384,938 9.20±0.439.20± 0.43 1.71±0.341.71± 0.34 TwinS (C) 14,81414,814 18.87±0.2418.87± 0.24 3.49±0.943.49± 0.94 Benchmarking datasets. We report the runtime of the models used in the benchmarking of datasets. The runtime are measured on Apple M1 Max. Table 9 shows the runtime/parameters used by each model to solve the benchmarking task. Table 11 shows for comparable configuration, the models runtime/parameters. Table 9: Runtime and parameter count, illustrated with Wisconsin dataset benchmarking of datasets, for the best configuration selected in Table 6. Train time is per-epoch, averaged over 100100 epochs (train repeats is 55) and inference time (infer repeats is 5050) is per forward pass over the full graph. Accuracy is reported to distinguish relevance across models. Model #Params Train (ms/epoch) Inference (ms) Accuracy MLP 117,701117,701 1.44±0.261.44± 0.26 0.19±0.050.19± 0.05 0.852±0.0500.852± 0.050 MagNet 111,429111,429 6.78±0.256.78± 0.25 3.76±0.743.76± 0.74 0.834±0.0400.834± 0.040 Dir-GNN 251,072251,072 4.15±0.124.15± 0.12 1.49±0.201.49± 0.20 0.674±0.0600.674± 0.060 FaberNet 6,029,3176,029,317 55.83±2.1355.83± 2.13 23.56±6.7923.56± 6.79 0.795±0.0520.795± 0.052 AdvDIFFormer 1,095,9411,095,941 13.45±0.5013.45± 0.50 2.67±0.362.67± 0.36 0.765±0.0480.765± 0.048 ADR-GNN 197,125197,125 18.80±0.5818.80± 0.58 8.50±0.298.50± 0.29 0.845±0.0450.845± 0.045 TwinS (S) 452,993452,993 4.19±0.814.19± 0.81 1.17±0.371.17± 0.37 0.871±0.0410.871± 0.041 TwinS (R) 1,508,6131,508,613 16.27±0.1616.27± 0.16 8.10±1.728.10± 1.72 0.874±0.0400.874± 0.040 TwinS (C) 678,849678,849 5.23±0.115.23± 0.11 1.80±0.361.80± 0.36 0.868±0.0380.868± 0.038 Table 10: Build time (ms) per TwinS-GCN model, per graph. Model Cornell Texas Wisconsin Chameleon-F Squirrel-F TwinS (S) 2020 1010 50.050.0 1,2301,230 22,54022,540 TwinS (R) 5050 4040 90.090.0 3,8103,810 58,54058,540 TwinS (C) 5050 4040 90.090.0 3,8103,810 58,54058,540 Table 11: Runtime and parameter count, illustrated with Wisconsin dataset benchmarking of datasets, for equivalent configuration (hidden=32, depth=4, dropout=0, K=2, bias=True). Train time is per-epoch, averaged over 5050 epochs (train repeats is 33) and inference time (infer repeats is 5050) is per forward pass over the full graph. Model #Params Train (ms/epoch) Inference (ms) MLP 56,80556,805 1.27±0.441.27± 0.44 0.17±0.040.17± 0.04 Dir-GNN 113,610113,610 1.72±0.091.72± 0.09 0.39±0.000.39± 0.00 MagNet 173,157173,157 10.67±0.2810.67± 0.28 5.85±0.305.85± 0.30 FaberNet 692,677692,677 14.61±0.2714.61± 0.27 6.29±0.126.29± 0.12 ADR-GNN 84,51784,517 15.03±0.1415.03± 0.14 6.03±0.676.03± 0.67 AdvDIFFormer 79,97379,973 6.93±0.366.93± 0.36 1.54±0.001.54± 0.00 TwinS (S) 227,220227,220 2.97±0.072.97± 0.07 0.68±0.000.68± 0.00 TwinS (R) 113,610113,610 1.58±0.081.58± 0.08 0.41±0.030.41± 0.03 TwinS (C) 340,830340,830 4.16±0.024.16± 0.02 1.07±0.031.07± 0.03 Appendix F Additional experiment for receptive field reach F.1 Graph signal transfer on cycle graph Settings. To gain further insight into the results of Section 5.1, we conduct a long-range dependency experiment on the directed cycle graph NC_N with N=200N=200 nodes. The goal of this task is to find the filter coefficients to transfer a localized signal from a source node to a target node at various distances. In particular, we study the target reconstruction error of the different filters (equation 8) underlying TwinS-GCN, as their behavior determines the receptive field of the model. We start with a localized signal defined as a gaussian centered at an arbitrarily chosen source node, and we select a target node such that the number of hops (shortest path distance) between them is r. The starting signal and the ground truth target signal are thus respectively defined as source[n]=e−(n−nsource)22σ2,target[n]=e−(n−ntarget)22σ2,with d(nsource,ntarget)=r,x_source[n]=e^- (n-n_source)^22σ^2, x_target[n]=e^- (n-n_target)^22σ^2, d(n_source,n_target)=r, (87) where σ is the bump’s width, and d shortest path distance from nsourcen_source to ntargetn_target. First, we generate different pairs of source and target nodes, spanning all distances up to the graph radius r=100r=100 (half the graph), and construct signals with σ=2σ=2. Then, for each fixed order k∈10,20,50k∈\10,20,50\, we estimate the coefficients of the sum and ratio filters hsumh_sum and hrath_rat introduced in Section 2.4 and the polynomial filter h of the same order. The coefficients are optimized for the normalized MSE (NMSE) between the filtered source signal and the ground truth target signal, where NMSE(target,^source,k)=‖target−^source,k‖22‖target‖22NMSE(x_target, x_source,k)= \|x_target- x_source,k \|_2^2 \|x_target \|_2^2 (88) where ^source,k=h^i,k()source x_source,k= h_i,k(L)\,x_source, and h^i,k h_i,k is the estimated k-th order filter of class hi∈h,hsum,hrath_i∈\h,h_sum,h_rat\. The results are shown in Figure 6. Figure 6: Long-range dependency experiment on cycle graph with increasing orders. Results. We clearly observe the effect of the FIR wall in Figure 6, where the kernel h generally struggles to reach 00 NMSE (perfect transfer) even when its order (degree) is higher than the distance between the source and the target, which in the Figure is shown as percentage of the graph radius). The kernel hsumh_sum also performs marginally better than h, consistent with the capacity of leveraging ↑L stated in Section 3.3.1, preventing oversmoothing. Hence in contrast to h, hsumh_sum keeps the initial signal localized, reducing NMSE. Finally, hrath_rat overcomes this wall as even when its order is largely lower than the distance, the signal transfer still happen (NMSE<1NMSE<1). It is however performing slightly worse than polynomial kernels when its order is larger than the distance, as expected from Theorem 5. F.2 Graph signal transfer on temperature graph Settings. Similarly to the experiment in Section F.1, we conduct a long-range dependency experiment on a temperature graph with N=212N=212 nodes, built as described in Figure 7b (Chan et al., 2026a). This time around, we set the source node randomly on the east side of the graph (source regions as the wind goes from east to west) and a target node randomly on the west side. The starting signal and the ground truth target signal are thus respectively defined as source[n]=e−([n]−[nsource])22σ2,target[n]=e−([n]−[ntarget])22σ2,with d(nsource,ntarget)=r,x_source[n]=e^- (c[n]-c[n_source])^22σ^2, x_target[n]=e^- (c[n]-c[n_target])^22σ^2, d(n_source,n_target)=r, (89) with [n]c[n] being the coordinates of node n, σ the bump’s width, and d the shortest path distance from nsourcen_source to ntargetn_target. We generate 1010 different pairs of source and target nodes, leading to 1010 pairs of signals with σ=1σ=1. An example of source-target pair is shown in Figure 7c. The coefficients are optimized for the NMSE between the filtered source signal and the ground truth target signal (equation 88). The results are shown in Figure 8a and 8b. We observe in particular that hrath_rat breaks the FIR wall and is able to better reach the target node despite a similar order than other kernels. With increasing order, the performance of hrath_rat improves, but is surpassed by hsumh_sum. (a) Measured wind velocity field (b) Derived temperature graph (c) Signal centered on source and target nodes used for the temperature displacement task Figure 7: a) Measured wind velocity field. b) Derived temperature graph. c) Signal centered on source and target nodes used for the temperature displacement task. (a) Estimated target signal (filter order: K=10K=10) (b) NMSE for target reconstruction as a function of filter order Figure 8: a) Estimated target signal using kernels h,hsum,hrath,h_sum,h_rat with K=10K=10. b) Mean NMSE across the 10 (source,target) pairs for estimated target signal as a function of filter order. Results. In this experiment we show similar takeaway as in Appendix F.1 but instead for a real graph, that is a temperature graph with direction given by the wind. We see once more a better NMSE performance for hrath_rat on low filter order regime w.r.t. h and hsumh_sum. On the other hand, hsumh_sum is able to achieve an overall better NMSE than h.