Gauge-Invariant Learnable Spectral Positional Encodings for Directed Graphs via Hermitian Block Krylov Subspaces
July 08, 2026
Abstract
Spectral positional encodings (PEs) for directed graphs face two obstacles: magnetic Laplacians require an \(O(n^3)\) Hermitian eigendecomposition per potential, and their complex eigenvectors are defined only up to unitary gauge, which prior work handles with basis-invariant architectures. We propose learnable spectral PEs of the form \(h_\theta(A_q)\,R\), where \(A_q\) is a normalized magnetic operator, \(h_\theta\) a learnable scalar spectral response, and \(R\) a block of random probes. Because the PE is a matrix function of the operator, it is gauge-invariant by construction. We compute it in a Hermitian block Krylov subspace from sparse matrix–vector products only, prove that \(k = O(\log(1/\varepsilon))\) block steps suffice uniformly over heat–resolvent response families, and give a covering-number argument for why low-dimensional structured families generalize where free per-eigenvalue weights overfit. On a directed SBM whose symmetrization is uninformative by construction, direction-blind PEs stay at chance while magnetic Krylov PEs converge to the exact-eigendecomposition oracle as the depth grows. The same probes yield gauge-invariant pairwise features with \(1/\sqrt{s}\) Monte-Carlo error, and the undirected \(q{=}0\) case improves heterophilous benchmarks over no-PE and polynomial baselines.
1 Introduction↩︎
Positional encodings are a key ingredient of graph transformers and a practical route to lifting the expressive power of message-passing GNNs beyond 1-WL [1], [2]. The dominant spectral construction, Laplacian PE, uses eigenvectors of the symmetric graph Laplacian and inherits three well-known problems: an \(O(n^3)\) eigendecomposition, instability under spectral perturbation, and sign/basis ambiguity of eigenvectors [3]–[5]. Learnable spectral PEs (LLPE) replace raw eigenvectors with a trainable filter of the spectrum [6], which improves task adaptivity but worsens the scalability problem, since the filter is applied in the eigenbasis.
Directed graphs sharpen every one of these issues. Direction matters in citation, program, and circuit graphs, yet the adjacency matrix is no longer symmetric, and the standard fix of symmetrization destroys the information of interest. The magnetic Laplacian [7] restores symmetry in the form of a complex Hermitian operator \(L_q\) whose phases encode edge directions, and Multi-\(q\) magnetic PEs [8] show that a set of potentials \(q\) is provably necessary to express directed walk profiles. But eigenvector-based magnetic PEs require a dense Hermitian eigendecomposition per potential, and complex eigenvectors are defined only up to per-eigenspace unitary gauge, which [8] address with dedicated basis-invariant architectures.
1.0.0.1 This paper.
We take a different route. We define the PE directly as a learnable matrix function of the magnetic operator applied to a block of random probes, \[\label{eq:pe} Z_q(\theta) \;=\; h_\theta(A_q)\,R \in \mathbb{C}^{n \times s}, \qquad R \sim \mathcal{CN}(0, \tfrac{1}{s} I),\tag{1}\] and approximate it in a Hermitian block Krylov subspace \(\mathcal{K}_k(A_q, R)\) using sparse matrix–vector products only. This design has three consequences. The probes are part of the randomized encoding: for a fixed draw, equivariance is conditional on permuting \(R\) with the nodes; if probes are freshly sampled after relabeling, the encoding is permutation equivariant in distribution because the Gaussian law is exchangeable (Proposition 2).
(1) Gauge invariance without invariant architectures. A matrix function \(h(A_q)\) is defined independently of any eigendecomposition; eigenvectors, their phases, and their basis choices never appear. The invariance that prior work recovers with special networks [3], [8] holds here by construction (Proposition 1).
(2) Approximation is provably cheap, uniformly over the learnable family. For structured response families built from heat kernels and resolvents, we show the block Krylov approximation error decays exponentially in the number of steps \(k\), uniformly over the family parameters (Theorem 1). The learnable PE family computed in \(k\) Krylov steps therefore \(\varepsilon\)-covers the exact LLPE family with \(k = O(\log(1/\varepsilon))\).
(3) Structure controls estimation. Free per-eigenvalue spectral weights [6] create an \(n\)-parameter hypothesis class per potential; heat–resolvent mixtures use \(O(m)\) parameters with \(m \ll n\). We give a covering-number argument (Proposition 3) and matching experiments showing free weights overfit at chance level in the low-label regime where structured families succeed.
Empirically, we construct a cyclic directed SBM whose symmetrized graph is (approximately) a homogeneous Erdős–Rényi graph, so all class signal lives in edge directions. Direction-blind PEs (\(q{=}0\)) are at chance; magnetic Krylov PEs recover the structure and converge monotonically to the exact-eigendecomposition oracle as \(k\) grows. On real directed node-classification benchmarks (WebKB, Wikipedia networks) we find that no spectral PE beats random probes. We report this negative result and take it as evidence that these standard benchmarks do not test directional structure, which further motivates controlled diagnostics.
2 Related Work↩︎
2.0.0.1 Spectral PEs and their invariance problem.
Laplacian PEs [1] require eigendecomposition and suffer sign/basis ambiguity, addressed architecturally by SignNet/BasisNet [3] and stable variants [4], [5]. PEARL [9] approximates equivariant functions of eigenvectors with message-passing and random or basis inputs, achieving linear complexity. Random feature propagation [10] also propagates random vectors but targets fixed diffusion, not a learnable spectral family with guarantees.
2.0.0.2 Learnable spectral filters.
Spectral GNNs learn polynomial [11]–[13] or set-to-set [14] filters; LLPE [6] learns per-eigenvalue weights for PEs. LanczosNet [15] used a Lanczos basis inside a GNN, without random probes, directed operators, or approximation guarantees for a learnable family.
2.0.0.3 Directed graph PEs.
MagNet [16] and magnetic PEs [7] introduced magnetic Laplacians to graph learning; Multi-\(q\) Mag-PE [8] showed multiple potentials are needed to express walk profiles and built basis-invariant networks over complex eigenvectors. All require dense Hermitian eigendecompositions per potential. Our PE inherits the multi-\(q\) insight but replaces eigenvectors with gauge-invariant matrix functions computed by sparse matvecs.
2.0.0.4 Krylov methods for matrix functions.
The convergence of (block) Lanczos approximations to \(f(A)b\) is classical [17], [18]: the error is governed by the best uniform polynomial approximation of \(f\) on the spectral interval. Our Theorem 1 is an application of this machinery made uniform over a learnable response family, as required for a learning guarantee.
3 Method↩︎
3.0.0.1 Setup.
Let \(G = (V, E)\) be a directed graph with \(n = |V|\), adjacency \(a_{uv} \in \{0,1\}\) and no self-loops. For a potential \(q \in [0, \tfrac12]\) define the symmetrized weights and phases \[\begin{align} a^{\mathrm{sym}}_{uv} &= \tfrac{1}{2}(a_{uv} + a_{vu}), & \Theta^{(q)}_{uv} &= 2\pi q\,(a_{uv} - a_{vu}), \end{align}\] the Hermitian matrix \(H_q(u,v) = a^{\mathrm{sym}}_{uv} e^{i\Theta^{(q)}_{uv}}\), the degree \(d_u = \sum_v a^{\mathrm{sym}}_{uv}\), and the normalized operator \[A_q \;=\; -\,D^{-1/2} H_q D^{-1/2}, \qquad \operatorname{spec}(A_q) \subseteq [-1, 1].\] \(A_q = L_q - I\) where \(L_q\) is the normalized magnetic Laplacian; \(q = 0\) recovers the symmetrized (direction-blind) operator.
3.0.0.2 Learnable spectral responses.
A response family is a set \(\mathcal{H} = \{h_\theta : [-1,1] \to \mathbb{R},\;\theta \in \Theta\}\). We study: (a) Chebyshev: \(h_\theta(\xi) = \sum_{m=0}^{M} c_m T_m(\xi)\); (b) heat–resolvent mixtures on \(\mu = \xi + 1 \in [0,2]\): \[\label{eq:hr} h_\theta(\xi) = \beta + \sum_{j=1}^{m} \alpha_j\, e^{-t_j \mu} + \sum_{j=1}^{m} \gamma_j\, (\mu + \tau_j)^{-1},\tag{2}\] with \(t_j \in [t_{\min}, t_{\max}]\), \(\tau_j \ge \tau_{\min} > 0\); (c) MLP: \(h_\theta(\xi) = \mathrm{MLP}(\phi(\xi))\) with fixed Fourier features \(\phi\); and (d) free weights (LLPE oracle): one parameter per eigenvalue. Per-head RMS normalization of \(h_\theta\) over the spectrum, with a learnable gain, puts all families on a common scale; we found this essential for fair comparisons (Section 5).
3.0.0.3 Hermitian block Krylov approximation.
Given complex Gaussian probes \(R \in \mathbb{C}^{n \times s}\), \(k\) steps of QR-stabilized block Lanczos with full reorthogonalization produce an orthonormal basis \(Q \in \mathbb{C}^{n \times r}\) (\(r \le ks\)) of \(\mathcal{K}_k(A_q, R) = \mathrm{span}\{R, A_q R, \dots, A_q^{k-1} R\}\), the projected operator \(T = Q^{\mathsf H} A_q Q\), and \(G = Q^{\mathsf H} R\). The PE for potential \(q\) is \[\label{eq:krylovpe} \widehat{Z}_q(\theta) = Q\, h_\theta(T)\, G = Q V h_\theta(\Xi) V^{\mathsf H} G,\tag{3}\] where \(T = V \Xi V^{\mathsf H}\) is the \(r \times r\) Hermitian eigendecomposition, computed once in \(O(r^3) = O(k^3 s^3)\), independent of \(n\). Multiple potentials \(\vec{q} = (q_1, \dots, q_Q)\) are handled with one Krylov cache each; features \([\mathrm{Re}\,\widehat Z_{q_i};\, \mathrm{Im}\,\widehat Z_{q_i}]\) are concatenated over potentials and heads and projected to the PE dimension. Training only differentiates through \(h_\theta\) on the fixed cache, so each epoch costs \(O(Q\,(r^2 s + n r s))\). The total precompute is \(O(Q\,k\,(\mathrm{nnz}(A) \,s + n k s^2))\) sparse work, and no dense \(n \times n\) object is ever formed.
4 Theory↩︎
Throughout, \(A \in \mathbb{C}^{n\times n}\) is Hermitian with \(\operatorname{spec}(A) \subseteq [-1,1]\), \(R \in \mathbb{C}^{n \times s}\), and \(Q, T, G\) are produced by \(k\) exact block Lanczos steps. Proofs are in Appendix 8.
Proposition 1 (Gauge invariance and equivariance). \(Z_q(\theta) = h_\theta(A_q) R\) depends on \(A_q\) only through the operator itself: it is invariant to any choice of eigenbasis (including per-eigenspace unitary gauge transformations of degenerate eigenspaces and global phase). Moreover, for any permutation matrix \(P\), \(h_\theta(P A_q P^\top) (PR) = P\, Z_q(\theta)\), and the same statements hold for the Krylov approximation 3 .
Unlike eigenvector PEs, there is no sign, phase, or basis ambiguity to repair, because no eigenvector is ever selected. The random probes should not be confused with deterministic node features: Proposition 2 states the two precise equivariance notions we use in experiments.
Proposition 2 (Random probes and pairwise entries). Let \(R = [r_1,\ldots,r_s]\) have iid complex Gaussian columns with \(\mathbb{E}r_t r_t^{\mathsf H}=\sigma^2 I\), and let \(F_a=h_a(A)\) and \(F_b=h_b(A)\) be bounded matrix functions of the same Hermitian operator.
For a fixed probe draw, \(F_a(PAP^\top)(PR)=P F_a(A)R\). If instead \(R'\) is freshly sampled after relabeling, \(F_a(PAP^\top)R' \stackrel{d}{=} P F_a(A)R\).
For any node pair \((i,j)\), \[\widehat K_{ab}(i,j)=\frac{1}{s\sigma^2} \sum_{t=1}^{s} (F_a r_t)_i\,\overline{(F_b r_t)_j}\] is an unbiased estimator of \((F_aF_b^{\mathsf H})_{ij}\) with variance \(O(\left\lVert e_i^{\mathsf T}F_a \right\rVert_2^2\left\lVert e_j^{\mathsf T}F_b \right\rVert_2^2/s)\).
Lemma 1 (Block polynomial exactness). For every polynomial \(p\) of degree at most \(k-1\), \(\;Q\,p(T)\,Q^{\mathsf H} R = p(A)\,R\).
Theorem 1 (Uniform Krylov approximation). Let \(\mathcal{H}\) be any response family and let \(E_{k-1}(h) = \inf_{\deg p \le k-1} \sup_{x \in [-1,1]} |h(x) - p(x)|\). Then \[\sup_{\theta \in \Theta} \left\lVert h_\theta(A) R - Q h_\theta(T) G \right\rVert_F \;\le\; 2 \left\lVert R \right\rVert_F\, \sup_{\theta \in \Theta} E_{k-1}(h_\theta).\] For the heat–resolvent family 2 with \(t_j \le t_{\max}\), \(\tau_j \ge \tau_{\min}\), and mixture coefficients bounded by \(B\) in \(\ell_1\), \[\sup_{\theta} E_{k-1}(h_\theta) \;\le\; C(B, t_{\max}, \tau_{\min})\; \rho^{\,k}, \qquad \rho < 1,\] where \(\rho = \big(1 + \tau_{\min} + \sqrt{\tau_{\min}^2 + 2\tau_{\min}} \big)^{-1}\) for the resolvent part and the heat part decays super-geometrically. Consequently \(k = O\!\big(\log(1/\varepsilon)\big)\) block steps suffice for a uniform \(\varepsilon \left\lVert R \right\rVert_F\) approximation of the entire learnable PE family.
Remark 1. Chebyshev responses of degree \(M \le k-1\) are reproduced exactly* (Lemma 1), so the Krylov PE family contains the polynomial LLPE family of degree \(k-1\) and \(\varepsilon\)-covers the analytic families. In floating point, full reorthogonalization keeps the bounds meaningful [18]; our implementation tracks orthogonality and projection residuals explicitly.*
Proposition 3 (Capacity of structured responses). Fix compact parameter sets and suppose the downstream loss is bounded and Lipschitz in the PE. For heat–resolvent responses with \(m\) components per potential, the induced predictor class has covering entropy \(O(Qm\log(C/\varepsilon))\). Consequently its empirical-process deviation on \(\ell\) labeled nodes scales as \(\widetilde{O}(\sqrt{Qm/\ell})\). A free per-eigenvalue response over rank \(r\) has covering entropy \(\Omega(Qr\log(C/\varepsilon))\) and admits the corresponding worst-case rate \(\widetilde{O}(\sqrt{Qr/\ell})\); in exact LLPE, \(r\) can be as large as \(n\).
Proposition 3 separates the two error sources: the Krylov depth \(k\) controls the approximation bias (Theorem 1) while the response family controls the estimation variance, and the two can be tuned independently.
5 Experiments↩︎
5.0.0.1 Naming convention.
Variants are named operator-solver-family. Operators: Mag (magnetic, grid \(\vec{q}\)) and Sym (\(q{=}0\), direction-blind); the prefix is dropped in the undirected experiments of Section 5.1, where the operator is the real symmetric \(A = L - I\). Solvers: Direct (explicit matrix recursion on \(A\)), Krylov (Eq. 3 ), and Exact (dense eigendecomposition oracle). Families: Cheb, Heat, HR (heat–resolvent), MLP, and Free (per-eigenvalue weights, the LLPE oracle); see Section 3. For example, Mag-Krylov-HR is the heat–resolvent family on magnetic operators computed in the Krylov subspace.
5.0.0.2 A controlled directed benchmark.
We build a cyclic directed SBM (dsbm): \(C{=}3\) classes on a directed cycle; forward edges (\(c \to c{+}1\)) appear with probability \(p_f\), backward with \(p_b \ll p_f\), and intra-class with \((p_f + p_b)/2\). By construction every unordered pair is connected in the symmetrized graph with probability \(\approx p_f + p_b\) regardless of classes: the symmetrized graph carries no class signal. Nodes have no features; the classifier sees the PE only. We use \(n = 600\), \(s = 32\) probes, potentials \(\vec{q} = (0, \tfrac16, \tfrac13)\), \(k = 10\) steps, 5 seeds, and report test accuracy at best validation.
| PE variant | Test accuracy | Precompute (s) |
|---|---|---|
| Random probes | 31.9 \(\pm\) 1.1 | 0.00 |
| LapPE (sym.eigvecs) | 32.6 \(\pm\) 1.5 | 2.04 |
| RWSE (sym.walks) | 33.4 \(\pm\) 0.2 | 0.21 |
| Mag-PE (eigvecs) | 51.1 \(\pm\) 2.3 | 6.50 |
| Sym-Krylov-Cheb (\(q{=}0\)) | 33.0 \(\pm\) 2.0 | 0.04 |
| Sym-Krylov-HR (\(q{=}0\)) | 34.5 \(\pm\) 2.5 | 0.04 |
| Sym-Krylov-MLP (\(q{=}0\)) | 31.1 \(\pm\) 2.2 | 0.04 |
| Mag-Krylov-Cheb | 83.2 \(\pm\) 7.8 | 0.12 |
| Mag-Krylov-Heat | 98.8 \(\pm\) 0.7 | 0.12 |
| Mag-Krylov-HR | 97.6 \(\pm\) 0.8 | 0.12 |
| Mag-Krylov-MLP | 98.4 \(\pm\) 1.0 | 0.12 |
| Mag-Exact-HR (oracle) | 98.2 \(\pm\) 0.9 | 0.06 |
| Mag-Exact-MLP (oracle) | 99.3 \(\pm\) 0.2 | 0.06 |
| Mag-Exact-Free (LLPE oracle) | 34.3 \(\pm\) 3.3 | 0.06 |
5.0.0.3 Main comparison (Table 2).
Random probes, the standard undirected baselines (LapPE on the symmetrized Laplacian, RWSE on symmetrized random walks), and all \(q{=}0\) variants sit at chance, validating the dataset design: the class signal lives entirely in edge directions. A fixed magnetic eigenvector PE (Mag-PE with the same potentials \(q\) and a deterministic gauge fix, in the spirit of [7], [8]) does see direction but reaches only \(51.1\%\): without a learnable spectral response, the raw eigenvector coordinates are a poor interface for the downstream MLP. All learnable magnetic Krylov families recover the directional structure (\(83\)–\(99\%\)). The free-weight LLPE oracle collapses to chance despite having the most capacity, consistent with the overfitting regime described by Proposition 3.
Figure 1: DSBM diagnostics, left to right: (a) test accuracy vs. block Krylov depth \(k\) (dashed lines: exact eigendecomposition oracles), the experimental counterpart of Theorem 1; (b) label efficiency in exact mode; free per-eigenvalue weights need far more labels than the \(O(m)\)-parameter heat–resolvent family, consistent with Proposition 3; (c) single-potential accuracy vs.\(q\), peaking near the cycle-aligned potential \(q = 1/C\) [8]; (d) PE precompute time vs.\(n\) at constant expected degree, measured on the same machine as Tables 12–13; the Krylov PE scales with sparse matvecs up to \(n{=}50\)k while dense Hermitian eigendecomposition is cubic and is skipped beyond \(n = 3000\)..
5.0.0.4 Depth sweep (Figure 1a).
Accuracy increases monotonically in \(k\): the heat–resolvent variant moves from \(49.5\%\) at \(k{=}2\) to \(85.9\%\) at \(k{=}3\), \(97.7\%\) at \(k{=}8\), and \(98.4\%\) at \(k{=}20\), against an exact oracle at \(98.2\%\). Most of the gap is closed by \(k \approx 8\), i.e.a basis of \(r \le ks \ll n\) dimensions, consistent with the \(O(\log 1/\varepsilon)\) rate of Theorem 1.
5.0.0.5 Label efficiency (Figure 1b).
With approximation error removed (exact mode), capacity ordering determines label efficiency in the direction suggested by Proposition 3: the \(O(m)\)-parameter heat–resolvent family reaches \(91.1\%\) with only \(2\%\) labeled nodes, the MLP response needs \({\sim}5\%\), and free per-eigenvalue weights stay near chance even with \(20\%\) labels (\(34.3\%\)) and reach only \(48.7\%\) at \(40\%\) labels. A separate probe-count sweep (Table 11) confirms the Monte-Carlo side of Proposition 2: increasing the number of probes from \(s{=}8\) to \(s{=}32\) roughly halves the standard deviation and improves accuracy.
5.0.0.6 Potential alignment (Figure 1c).
With a single potential the response is sharply resonant: accuracy is at chance for \(q \le 0.1\), rises through \(72.9\%\) at \(q{=}0.15\) to a \(98\)–\(99\%\) plateau for \(q \in [0.2, 1/3]\) (the band aligned with the \(C{=}3\) cyclic structure), and collapses back to chance as \(q \to 1/2\). Since the informative band is task-dependent and unknown a priori, a small grid of potentials, each with its own cheap Krylov cache, inherits the best scale, supporting the multi-\(q\) design of [8] within a learnable-response framework.
5.0.0.7 Scaling (Figure 1d).
At constant expected degree the Krylov precompute is essentially flat (\(0.09\)s at \(n{=}600\) to \(0.17\)s at \(n{=}50\mathrm{k}\), all three potentials included), while dense Hermitian eigendecomposition grows cubically (\(0.05\)s at \(n{=}600\), \(0.22\)s at \(n{=}3000\)) and is infeasible at the \(50\)k scale, where the dense eigenvectors of a single potential alone would occupy tens of gigabytes. This is the regime in which the small-graph SPE comparison no longer applies: exact magnetic spectra become the bottleneck, while Krylov still uses sparse matvecs and stays below \(0.2\) seconds in our diagnostic. Accuracy at fixed \(k{=}10\) drifts only mildly with \(n\) (\(98.1\%\) at \(n{=}600\) vs.\(95.6\%\) at \(n{=}6000\)), recoverable by scaling \(k\) (or, at larger scales, \(s\)) with the spectral resolution required. Table 12 extends this diagnostic to \(50\)k-node directed graphs: exact magnetic PEs are skipped as infeasible, while Mag-Krylov remains well above random and direction-blind baselines. A 50k-node ablation (Table 13) further separates the two approximation parameters: increasing \(k\) from \(10\) to \(20\) at fixed \(s\) does not improve accuracy, whereas increasing \(s\) from \(32\) to \(64\) recovers most of the large-graph drop, raising accuracy from \(60.7\%\) to \(69.8\)–\(70.6\%\) while keeping preprocessing below one second. This suggests that, in this regime, the dominant residual error is random-probe resolution rather than Krylov depth.
| PE variant | Chameleon | Cornell | Texas | Wisconsin |
|---|---|---|---|---|
| Random probes | 41.5 \(\pm\) 1.9 | 61.3 \(\pm\) 2.4 | 63.3 \(\pm\) 1.8 | 73.6 \(\pm\) 1.3 |
| Sym-Krylov-HR (\(q{=}0\)) | 38.1 \(\pm\) 1.6 | 61.3 \(\pm\) 1.7 | 61.9 \(\pm\) 2.4 | 71.8 \(\pm\) 0.9 |
| Sym-Krylov-MLP (\(q{=}0\)) | 30.2 \(\pm\) 1.4 | 63.6 \(\pm\) 3.1 | 65.0 \(\pm\) 1.6 | 71.5 \(\pm\) 1.5 |
| Mag-Krylov-HR | 30.8 \(\pm\) 1.8 | 61.1 \(\pm\) 5.4 | 60.8 \(\pm\) 2.1 | 72.3 \(\pm\) 4.0 |
| Mag-Krylov-MLP | 30.9 \(\pm\) 3.9 | 57.2 \(\pm\) 3.7 | 61.6 \(\pm\) 3.9 | 70.3 \(\pm\) 1.7 |
5.0.0.8 Real directed graphs (Table 3).
On WebKB and Wikipedia networks (loaded without symmetrization), all spectral PEs, directional or not, perform within noise of the random-probe baseline, both with raw features and in PE-only mode. We report this as a negative result: the node labels of these standard heterophilic benchmarks are not encoded in global directed structure, so they cannot distinguish between directional PEs. This mirrors why [8] evaluate on direction-critical tasks (distance prediction, sorting networks, circuits) rather than node classification. We therefore complement this node-level diagnostic with the directed pairwise regression suite in Section 5.2.
5.1 The Undirected Special Case (\(q = 0\))↩︎
At \(q = 0\) the framework reduces to the real symmetric operator \(A = L - I\), so every framework-level claim can be validated independently of direction. We summarize experiments run with the undirected instantiation (synthetic heterophilous graphs, MLP backbone, setup in Appendix 9).
| Method | Acc.(%) | \(\Delta\) |
|---|---|---|
| Random probes | 32.83 \(\pm\) 1.72 | \(-3.24\) |
| Direct-Cheb | 36.07 \(\pm\) 4.18 | – |
| Krylov-Cheb | 36.07 \(\pm\) 4.18 | \(+0.00\) |
| Krylov-MLP | 33.48 \(\pm\) 1.49 | \(-2.60\) |
| Krylov-HR | 47.25 \(\pm\) 6.87 | \(+11.18\) |
| Exact-HR (oracle) | 65.60 \(\pm\) 9.53 | \(+29.52\) |
| Krylov-Heat | 72.00 \(\pm\) 6.01 | \(+35.93\) |
| Exact-Heat (oracle) | 74.77 \(\pm\) 7.29 | \(+38.70\) |
| \(M\) | \(k\) | Direct-Cheb | Krylov-Cheb | Krylov-Heat |
|---|---|---|---|---|
| 2 | 6 | 33.62 \(\pm\) 1.70 | 33.65 \(\pm\) 1.71 | 68.01 \(\pm\) 7.26 |
| 4 | 8 | 36.48 \(\pm\) 3.75 | 36.48 \(\pm\) 3.75 | 69.85 \(\pm\) 6.57 |
| 8 | 12 | 36.07 \(\pm\) 4.18 | 36.07 \(\pm\) 4.18 | 72.00 \(\pm\) 6.01 |
| 16 | 20 | 34.61 \(\pm\) 3.61 | 34.61 \(\pm\) 3.61 | 73.51 \(\pm\) 6.22 |
| Method | Minesweeper | Tolokers | Questions |
|---|---|---|---|
| No PE / MLP | 50.26 \(\pm\) 0.98 | 73.49 \(\pm\) 0.84 | 66.89 \(\pm\) 1.42 |
| LapPE | 52.46 \(\pm\) 1.02 | 79.92 \(\pm\) 0.99 | 71.48 \(\pm\) 1.19 |
| RWSE | 52.42 \(\pm\) 1.11 | 77.43 \(\pm\) 0.60 | – |
| RFP | 50.39 \(\pm\) 0.16 | 76.47 \(\pm\) 0.49 | 67.88 \(\pm\) 0.81 |
| Direct-Cheb | 50.29 \(\pm\) 1.65 | 72.52 \(\pm\) 0.86 | 68.07 \(\pm\) 1.13 |
| Krylov-Heat | 51.19 \(\pm\) 1.48 | 76.04 \(\pm\) 1.05 | 69.24 \(\pm\) 1.61 |
Three observations transfer directly. (i) Exactness in practice: Krylov-Cheb and Direct-Cheb agree to all reported digits at every order (Tables 4 and 5), so the
Krylov compression loses nothing on the polynomial family while costing a small dense problem instead of \(M{+}1\) stored \(n \times s\) matrices. (ii) Structured families under a fixed
budget: the heat family gains \(+36\) points over the best polynomial at identical Krylov cost and sits within \(3\) points of its exact oracle. (iii) Real data: on
heterophilous benchmarks (Table 6) Krylov-Heat improves over no-PE and polynomial baselines on all three datasets, while LapPE is the strongest method overall. We read this as a property of the benchmarks rather than
of the encodings: these undirected graphs have rich node features, and a fixed low-pass eigenvector summary suffices; there is no directional signal for a learnable response to exploit (cf.the chance performance of all symmetric PEs in Table 2). The Krylov PE remains an order of magnitude cheaper to precompute (Questions: \(1.2\)s vs.\(14.5\)s for LapPE’s sparse eigsh) and carries no
sign/basis-ambiguity caveats. RWSE’s sparse-power computation does not complete on the largest graph. Undirected scaling experiments (up to \(n{=}9000\): Krylov \(0.19\)s vs.exact \(30.9\)s precompute, with a \(4\)-point accuracy gap at fixed \(k\)) appear in Appendix 9.
5.2 Directed Pairwise Regression↩︎
| Method | Precompute | Stored PE state | Pairwise readout |
|---|---|---|---|
| Exact Mag-PE + SPE | \(O(Q n^3)\) dense eigensolve | \(O(Q n k_{\mathrm{eig}})\)–\(O(Qn^2)\) | \(O(P Q k_{\mathrm{eig}}^2)\) |
| Mag-Krylov + probes | \(O(Qk(\mathrm{nnz}\,s + nks^2))\) | \(O(Qns + Qk^2s^2)\) | \(O(P Q H^2 s)\) |
| Measured scaling | exact skipped at \(n{>}3000\); Krylov \(0.17\)–\(0.76\)s at \(50\)k | dense eigenvectors vs.probe blocks | same cached probes |
We evaluate on the directed shortest-path-distance (SPD) benchmark of [8]: random connected DAGs (\(16\)–\(63\) nodes for training, \(64\)–\(71\) for testing), where the model regresses \(\mathrm{spd}(i,j)\) for all reachable pairs. We integrate our encoder into their public framework and keep dataset, predictor head, batch size, and epochs identical across methods (\(40\)k training graphs; \(3\) seeds). Their headline method couples multi-\(q\) magnetic eigenvectors with the SPE invariant network [5], which consumes the full eigendecomposition and per-pair products of eigenvector entries, which is exact but costs \(O(n^3)\) per graph.
The matrix-function view offers a readout that eigenvector PEs do not: for any pair \((i,j)\), the normalized probe inner products \(\frac{1}{s\sigma^2}\sum_{t=1}^{s}(h_a(\mathbf{L}_q)\mathbf{z}_t)_i \overline{(h_b(\mathbf{L}_q)\mathbf{z}_t)_j}\) are unbiased Hutchinson estimates of the gauge-invariant matrix entries \((h_a h_b^{*})(\mathbf{L}_q)_{ij}\), computed from the same precomputed probes as the node-level PE at no extra preprocessing cost. Table 8 isolates the contribution of each design axis: the pairwise readout (vs.node-level features alone), the response-family capacity \(H\), the cross-head products (\(H^2\) effective responses), and the probe count \(s\) controlling the \(O(1/\sqrt{s})\) Monte-Carlo error of Proposition 2.
| Configuration | Test MSE |
|---|---|
| Mag-PE eigvecs + MLP encoder | \(0.314 \pm 0.003\) |
| Mag-Krylov, node-level only (\(H{=}4\), \(s{=}8\)) | \(0.272 \pm 0.008\) |
| + pairwise readout | \(0.127 \pm 0.003\) |
| + capacity \(H{=}16\) | \(0.069 \pm 0.002\) |
| + probes \(s{=}32\) | \(0.051 \pm 0.001\) |
| + cross-head products (\(H^2\)) | \(0.042 \pm 0.001\) |
| Mag-PE + SPE (exact eigendecomposition) | \(0.0037 \pm 0.0004\) |
Each row corresponds to a quantity in our analysis, and each improves the error in the predicted direction. The pairwise readout, which eigenvector-PE baselines do not support, halves the error twice over; enlarging the response family helps until the Monte-Carlo noise floor binds; and quadrupling \(s\) buys the \(\approx 2\times\) noise reduction the \(1/\sqrt{s}\) rate predicts, with diminishing returns thereafter (at \(s{=}128\) the error improves only to \(\approx 0.033\) at \(4\times\) the memory). The exact-spectrum SPE remains an order of magnitude better on these small graphs: SPD demands near-exact pairwise resolvent information, and a method that pays \(O(n^3)\) for the exact spectrum will win whenever that price is payable. The comparison quantifies the trade-off: \(O(\mathrm{nnz}\cdot M\cdot r)\) precompute with a controlled \(1/\sqrt{s}\) error, versus exact spectral readouts that do not scale past small graphs.
We next keep the same protocol and evaluate the full directed pairwise suite: shortest-path distance (SPD), longest-path distance (LPD), and the step-4 walk-profile vector (WP). These targets stress different global properties of the directed DAGs, while reusing the same encoder implementation and the same pairwise cross-head readout. Table 9 shows that the Krylov matrix-function PE remains far stronger than the standard Mag-PE encoder, improving the error by \(7.5\times\) on SPD, \(2.1\times\) on LPD, and \(7.4\times\) on WP. SPE remains the best method on average, but its variance is much larger on the harder LPD/WP targets: in both cases one seed is close to, or worse than, the structured Krylov model.
| Method | SPD | LPD | WP |
|---|---|---|---|
| Mag-PE eigvecs + MLP encoder | \(0.314 \pm 0.003\) | \(0.446 \pm 0.004\) | \(2.978 \pm 0.006\) |
| Mag-Krylov + pairwise cross | \(0.042 \pm 0.001\) | \(0.210 \pm 0.004\) | \(0.401 \pm 0.005\) |
| Mag-PE + SPE (exact eigendecomposition) | \(0.0037 \pm 0.0004\) | \(0.102 \pm 0.058\) | \(0.306 \pm 0.132\) |
Thus we do not interpret Tables 8–9 as a small-graph leaderboard against SPE. Exact spectra are the right tool when the graphs are tiny enough and the task demands near-exact pairwise spectral entries. The contribution is instead a controlled approximation path: the same gauge-invariant matrix-function PE supports node features and pairwise readouts, exposes a probe count \(s\) with \(1/\sqrt{s}\) error, and trades some accuracy for sparse, memory-light preprocessing that remains usable when dense eigendecomposition is no longer an option.
6 Discussion and Limitations↩︎
Our diagnostic benchmark is intentionally adversarial to symmetrization; real graphs mix directional and undirected signal, where the margin between \(q{=}0\) and magnetic PEs naturally narrows. The potentials \(\vec{q}\) are currently a fixed grid: a learnable \(q\) requires differentiating through the operator and rebuilding the Krylov cache, an interesting direction we leave open. Our theory covers approximation and a covering-number capacity bound, not the optimization dynamics of the response families; the RMS normalization we found necessary in practice (Section 5) suggests conditioning of the spectral parameterization deserves study in its own right. Finally, extending beyond the pairwise regression benchmarks of Section 5.2 to the circuit and sorting-network tasks of [8], which require transformer backbones, is the natural next experimental step.
7 Conclusion↩︎
Learnable spectral PEs for directed graphs do not require eigenvectors. A matrix-function formulation over magnetic operators is gauge-invariant by construction, admits uniform Krylov approximation guarantees with \(O(\log 1/\varepsilon)\) sparse steps, and separates the design into three independent parameters: the Krylov depth \(k\) (approximation bias), the response-family capacity, and the probe count \(s\) (Monte-Carlo resolution). This decomposition, validated here on directed graphs, is not specific to the directed setting.
8 Proofs↩︎
8.1 Proposition 1↩︎
A matrix function of a Hermitian matrix is defined by the spectral calculus \(h(A) = \sum_i h(\lambda_i) \Pi_i\) where \(\Pi_i\) are the (unique) orthogonal spectral projectors; any eigenbasis choice within an eigenspace yields the same \(\Pi_i\), hence the same \(h(A)\) and the same \(Z = h(A)R\). Permutation equivariance follows from \(h(PAP^\top) = P h(A) P^\top\) for unitary \(P\). The Krylov statement follows because \(Q\) is constructed from \(\{A^j R\}_{j<k}\): replacing \((A, R) \mapsto (PAP^\top, PR)\) maps \(Q \mapsto PQ\) (up to the same internal QR conventions), leaving \(T = Q^{\mathsf H} A Q\) and \(G = Q^{\mathsf H} R\) unchanged and mapping \(\widehat Z \mapsto P \widehat Z\). 0◻
8.2 Proposition 2↩︎
For the conditional equivariance claim, use the identity \(h(PAP^\top)=P h(A)P^\top\) and multiply by \(PR\). For distributional equivariance, an iid complex Gaussian matrix satisfies \(P^\top R' \stackrel{d}{=} R\), hence \(h(PAP^\top)R' = P h(A) P^\top R' \stackrel{d}{=} P h(A)R\).
For the pairwise estimator, write the \(i\)th row of \(F_a\) as \(u^{\mathsf H}=e_i^{\mathsf T}F_a\) and the \(j\)th row of \(F_b\) as \(v^{\mathsf H}=e_j^{\mathsf T}F_b\). A single probe contributes \(X_t=(u^{\mathsf H}r_t)\overline{(v^{\mathsf H}r_t)}.\) Since \(\mathbb{E}r_t r_t^{\mathsf H}=\sigma^2 I\), \[\mathbb{E}X_t = u^{\mathsf H}(\mathbb{E}r_t r_t^{\mathsf H})v = \sigma^2 u^{\mathsf H}v = \sigma^2 (F_aF_b^{\mathsf H})_{ij}.\] Thus \((s\sigma^2)^{-1}\sum_t X_t\) is unbiased. For the variance, Isserlis’ formula for circular complex Gaussians gives \[\begin{align} \mathbb{E}|X_t|^2 &= \mathbb{E}|u^{\mathsf H}r_t|^2 |v^{\mathsf H}r_t|^2 \\ &= \sigma^4\big(\left\lVert u \right\rVert_2^2\left\lVert v \right\rVert_2^2 + |u^{\mathsf H}v|^2\big) \\ &\le 2\sigma^4\left\lVert u \right\rVert_2^2\left\lVert v \right\rVert_2^2 . \end{align}\] After subtracting \(|\mathbb{E}X_t|^2\) and averaging independent probes, \[\mathrm{Var}(\widehat K_{ab}(i,j)) \le \frac{2}{s}\, \left\lVert e_i^{\mathsf T}F_a \right\rVert_2^2\left\lVert e_j^{\mathsf T}F_b \right\rVert_2^2 .\] 0◻
8.3 Lemma 1↩︎
By induction: \(A^j R \in \mathcal{K}_k\) for \(j \le k-1\), and on \(\mathcal{K}_k\) the compression \(T = Q^{\mathsf H} A Q\) satisfies \(Q T^j Q^{\mathsf H} R = A^j R\) for \(j \le k-1\) (standard block Lanczos argument; rank deflation only shrinks the space when it is already invariant). Linearity extends to polynomials. 0◻
8.4 Theorem 1↩︎
Fix \(\theta\) and any polynomial \(p\) with \(\deg p \le k-1\). By Lemma 1, \[\begin{align} &\left\lVert h_\theta(A)R - Q h_\theta(T) G \right\rVert_F \\ &\quad\le \left\lVert (h_\theta - p)(A) R \right\rVert_F + \left\lVert Q (h_\theta - p)(T) G \right\rVert_F \\ &\quad\le \big( \left\lVert (h_\theta - p)(A) \right\rVert_2 + \left\lVert (h_\theta - p)(T) \right\rVert_2 \big) \left\lVert R \right\rVert_F \\ &\quad\le 2 \sup_{x \in [-1,1]} |h_\theta(x) - p(x)| \cdot \left\lVert R \right\rVert_F, \end{align}\] using \(\operatorname{spec}(T) \subseteq [\lambda_{\min}(A), \lambda_{\max}(A)] \subseteq [-1,1]\) (Cauchy interlacing for compressions) and the spectral mapping theorem. Taking the infimum over \(p\) and the supremum over \(\theta\) gives the first claim. For the rates: each resolvent \(x \mapsto (x + 1 + \tau)^{-1}\) is analytic on a Bernstein ellipse with parameter \(\rho^{-1} = 1 + \tau + \sqrt{\tau^2 + 2\tau}\), giving \(E_{k-1} \le C_\tau \rho^k\) [19]; each heat term \(x \mapsto e^{-t(x+1)}\) is entire, with \(E_{k-1}(e^{-t(\cdot+1)}) \le C e^{-k^2/(2t)}\)-type super-geometric decay for \(k \gtrsim t\). Sub-additivity of \(E_{k-1}\) over the mixture with \(\ell_1\) coefficient bound \(B\) and the compactness of \([t_{\min}, t_{\max}] \times [\tau_{\min}, \infty)\) in the relevant constants yield the uniform bound. 0◻
8.5 Proposition 3↩︎
Let \(\Theta \subset \mathbb{R}^{d}\) be a compact parameter set for a response family, and suppose \(\theta\mapsto h_\theta(\xi)\) is \(L_h\)-Lipschitz uniformly over \(\xi\in[-1,1]\). Because the Krylov cache is fixed during training, \[\widehat Z(\theta)-\widehat Z(\theta') =Q\big(h_\theta(T)-h_{\theta'}(T)\big)G\] is Lipschitz in \(\theta\) with constant at most \(L_h\left\lVert Q \right\rVert_2\left\lVert G \right\rVert_F\). Composing with a downstream predictor and a bounded \(L_\ell\)-Lipschitz loss preserves Lipschitzness up to constants. Therefore an \(\varepsilon\)-net of \(\Theta\) induces a \(C\varepsilon\)-net of the loss class. For a bounded \(d\)-dimensional compact parameter set, the covering number is \[\mathcal{N}(\varepsilon,\Theta,\left\lVert \cdot \right\rVert_2)\le (C/\varepsilon)^d, \qquad \log\mathcal{N}=O(d\log(C/\varepsilon)).\] Standard covering-number bounds for bounded losses then yield uniform deviation \[\sup_{f\in\mathcal{F}} |R(f)-\widehat R_\ell(f)| =\widetilde{O}\!\left(\sqrt{d/\ell}\right).\] For heat–resolvent mixtures, \(d=O(Qm)\) (a constant number of parameters per component and per potential). For free per-eigenvalue weights over a rank-\(r\) spectral representation, \(d=\Omega(Qr)\), and in exact LLPE one may have \(r=n\). This gives the two rates in the proposition and explains the observed low-label overfitting of the free family. 0◻
9 Undirected experimental details↩︎
The undirected experiments of Section 5.1 use the real symmetric operator \(A = L - I = -D^{-1/2} \mathrm{Adj}\, D^{-1/2}\) (the \(q{=}0\) instantiation) and random Gaussian probes shared across all variants. Synthetic graphs are heterophilous SBMs with \(n = 1200\), \(3\) classes, \(p_{\mathrm{in}} = 0.008\), \(p_{\mathrm{out}} = 0.035\) (probabilities scaled \(O(1/n)\) in the scaling study), no node features (PE-only classification); \(s = 32\) probes, PE dimension \(64\), \(8\) heads, \(8\) heat/resolvent components, \(2\)-layer MLP backbone with hidden width \(128\), dropout \(0.1\), AdamW with learning rate \(10^{-3}\) and weight decay \(10^{-5}\), up to \(500\) epochs with patience \(150\) and best-validation restoration. Controlled-comparison results average 10 seeds, scaling results 3 seeds; heterophilous benchmarks use the standard splits of the respective datasets with original features and ROC-AUC for the binary tasks. All synthetic numbers were re-verified by rerunning the released configuration end to end. “Krylov-Heat” denotes the heat-only structured family; “Krylov-HR” the heat–resolvent mixture; oracles use a full eigendecomposition. The Chebyshev-order sweep ties the Krylov budget to the polynomial order via \(k = M + 4\) so that polynomial exactness (Lemma 1) is guaranteed with margin.
| \(n\) | Method | Acc.(%) | Precompute | Train time |
|---|---|---|---|---|
| 1200 | Direct-Cheb | 36.11 \(\pm\) 1.73 | 0.006s | 2.01s |
| 1200 | Krylov-Heat | 69.44 \(\pm\) 0.66 | 0.034s | 2.07s |
| 1200 | Exact-Heat | 71.75 \(\pm\) 1.25 | 0.043s | 4.14s |
| 3000 | Direct-Cheb | 54.89 \(\pm\) 8.71 | 0.014s | 6.29s |
| 3000 | Krylov-Heat | 81.08 \(\pm\) 1.46 | 0.077s | 5.18s |
| 3000 | Exact-Heat | 81.14 \(\pm\) 0.59 | 0.528s | 33.98s |
| 5000 | Direct-Cheb | 66.74 \(\pm\) 8.13 | 0.019s | 8.36s |
| 5000 | Krylov-Heat | 86.11 \(\pm\) 1.78 | 0.102s | 6.67s |
| 5000 | Exact-Heat | 89.39 \(\pm\) 0.79 | 3.441s | 96.14s |
| 9000 | Direct-Cheb | 73.61 \(\pm\) 4.39 | 0.033s | 12.55s |
| 9000 | Krylov-Heat | 87.04 \(\pm\) 0.36 | 0.192s | 12.81s |
| 9000 | Exact-Heat | 91.01 \(\pm\) 0.54 | 30.926s | 310.16s |
| Probe count | Accuracy |
|---|---|
| \(s=8\) | \(73.31 \pm 7.01\) |
| \(s=32\) | \(86.86 \pm 3.23\) |
| \(n\) | Random | Sym-Krylov-HR | Mag-Krylov-HR |
|---|---|---|---|
| 10k | \(33.56 \pm 0.39\) | \(33.60 \pm 0.21\) | \(\mathbf{79.02 \pm 2.46}\) |
| 20k | \(33.64 \pm 0.05\) | \(33.61 \pm 0.33\) | \(\mathbf{73.88 \pm 4.24}\) |
| 50k | \(33.24 \pm 0.21\) | \(33.26 \pm 0.05\) | \(\mathbf{60.72 \pm 1.59}\) |
| Setting | Mag-Krylov-HR | Precompute |
|---|---|---|
| \(k{=}10,\;s{=}32\) | \(60.72 \pm 1.59\) | \(0.17\)s |
| \(k{=}20,\;s{=}32\) | \(60.52 \pm 1.66\) | \(0.38\)s |
| \(k{=}10,\;s{=}64\) | \(69.78 \pm 2.70\) | \(0.33\)s |
| \(k{=}20,\;s{=}64\) | \(\mathbf{70.62 \pm 2.73}\) | \(0.76\)s |
10 Implementation details↩︎
Block Lanczos uses QR with rank revealing via \(|R_{ii}|\) thresholding (tolerance \(10^{-7}\)) and two-pass reorthogonalization; we log orthogonality error \(\left\lVert Q^{\mathsf H} Q - I \right\rVert_F / \sqrt{r}\) (\({\sim}10^{-7}\) in float32), Hermiticity error of \(T\), and the projection residual \(\left\lVert AQ - QT \right\rVert_F / \left\lVert AQ \right\rVert_F\).
10.0.0.1 Directed experiments.
Cyclic DSBM with \(n = 600\) (except the scaling study), \(C = 3\) classes, \(p_f = 0.05\), \(p_b = 0.005\), no node features; potentials \(\vec{q} = (0, \tfrac16, \tfrac13)\), \(s = 32\) complex Gaussian probes, \(k = 10\) block steps (except the depth sweep), PE dimension \(32\), \(4\) heads, \(6\) heat/resolvent components, MLP responses with \(8\) fixed Fourier frequencies; stratified splits with \(10\%\) train / \(20\%\) validation. Training: AdamW, lr \(2\cdot10^{-3}\), weight decay \(10^{-4}\), dropout \(0.5\), up to 300 epochs with patience 50 on validation accuracy, model restored to the best epoch. Per-head RMS normalization with learnable gain is applied to all response families (Section 3); disabling it degrades Cheb and MLP families drastically (e.g.Mag-Krylov-Cheb \(34\% \to 93\%\) with normalization) while leaving the internally normalized heat–resolvent family unchanged, which is why we treat it as part of the parameterization rather than a tunable trick. The backbone is a 2-layer MLP on \([X; \mathrm{PE}]\); all variants share probes, splits, seeds, and backbone. Code is in the supplement. Large-scale directed runs use the same directed SBM with edge probabilities scaled as \(1/n\) to keep expected degree constant: \((p_f,p_b)=(0.0024,0.00030)\) for \(n{=}10\)k, \((0.0012,0.00015)\) for \(n{=}20\)k, and \((0.00048,0.00006)\) for \(n{=}50\)k. The large-scale table uses \(80\) epochs with patience \(20\) and skips exact variants above \(n=3000\). The probe-stability table fixes \(n{=}3000\) and compares \(s{=}8\) vs.\(s{=}32\) over 10 seeds.