HARMONIA: Interpretable Graph Learning through Mixtures of Neural Bases

1 INTRODUCTION

• Post-hoc explanations approximate or summarize a trained GNN rather than determine its prediction mechanism, motivating interpretable-by-design models with explicit feature and structural effects.Their faithfulness depends on how well the auxiliary explanation captures the underlying GNN's behavior.Source: §1 · p. 1
• Separate networks per feature can scale poorly, while one shared basis can underfit heterogeneous feature responses or require expanding capacity for every feature.Graph additive models retain feature-wise response functions while aggregating information from the graph.Source: §1 · pp. 1–2
• Shortest-path summaries omit transition-probability differences across walk lengths and paths, while storing RRWP for all node pairs requires quadratic memory.These limitations make richer structural modeling impractical for large graphs without a more efficient computation.Source: §1 · p. 2
• HARMONIA combines specialized feature bases and sparse RRWP aggregation to scale across feature dimension and graph size while preserving an explicit additive prediction structure.MoNB routes features to specialized functional spaces; SRA computes linear RRWP aggregation through sparse propagation without dense pairwise representations.Source: §1 · p. 2
• The paper analyzes when specialized bases overcome routing error and when random-walk spectral properties make hop effects identifiable or increasingly redundant.It also proves SRA reduces computation from O(Tn^2d) to O(T|E|d), and evaluates predictive performance, explanation faithfulness, and scalability.Source: §1 · p. 2

2 PRELIMINARIES AND PROBLEM SETUP

• Neural Basis Models represent each feature response by combining shared learnable basis functions with feature-specific coefficients.This reduces the number of nonlinear functions learned, but constrains all responses to one common functional space that may not fit heterogeneous features.Source: §2 · pp. 2–3
• RRWP describes node-pair relations through random-walk transition probabilities at multiple walk lengths.A direct implementation requires O(Tn^2) memory and can require O(Tn^2d) computation.Source: §2 · p. 3
• The predictor seeks an explicit decomposition over original features, with each feature response at a source node modulated by a feature-specific topology function for its target node.The resulting logit has an exact feature–node decomposition, but flexible responses and multi-hop influence without dense computation remain challenges.Source: §2 · p. 3

3.1 FEATURE SPECIALIZATION THROUGH MIXTURES OF NEURAL BASES

• MoNB replaces one global basis space with multiple shared experts, routing features to specialized functional spaces while allowing features with similar responses to reuse bases.Feature embeddings generate sparse expert scores, and routing changes each feature's effective basis directions without increasing its B-dimensional coefficient vector.Source: §3.1 · pp. 3–4
• Specialization improves over the best globally shared B-dimensional space exactly when routing error is smaller than the specialization gain.A larger specialization gain allows more routing error before specialization loses its advantage.Source: §3.1 · p. 5

3.2 NORMALIZED RANDOM-WALK STRUCTURAL MODELING

• Normalized RRWP lets features emphasize different propagation depths, but hop weights have unique interpretations only when the corresponding walk operators are linearly independent.Under the stated reversible-walk assumptions, the effective structural horizon identifies when additional hops become structurally redundant; it does not determine the optimal walk length.Source: §3.2 · pp. 5–6
• The effective structural horizon is graph-dependent and can be estimated before training from the transition operator's nonstationary spectral radius.For a chosen resolution, the theorem gives a walk-length budget; its spectral analysis assumes an ergodic reversible random walk, not required by the architecture.Source: §3.2 · p. 6

3.3 SPARSE RRWP AGGREGATION

• SRA avoids quadratic pairwise materialization by repeatedly applying the original sparse transition matrix, even when multi-step transition powers become dense.This reordering computes the same RRWP aggregation exactly and preserves feature-, source-node-, and walk-length-specific prediction decompositions.Source: §3.3 · pp. 6–7
• SRA uses O(T|E|d) computation, compared with O(Tn^2) memory and O(Tn^2d) computation for explicit pairwise materialization.For ergodic reversible random walks, the graph-dependent horizon can set T; for general directed graphs that spectral bound need not hold, while SRA remains exact.Source: §3.3 · p. 7

4.1 REAL-WORLD BENCHMARKS ACCURACY AND SCALABILITY

• HARMONIA outperforms interpretable baselines on five of six datasets, scales to large graphs where some baselines run out of memory, and is best overall on Tolokers and ogbn-arxiv.It remains competitive on Cora and ogbn-products.Source: Table 1 · p. 7, §4.1 · p. 8
• Figure 2 supports choosing T = Hε + 1, beyond which consecutive propagation operators differ by at most ε, avoiding an unnecessarily large propagation budget.The associated cost grows with T as O(T|E|d).Source: §4.1 · p. 8, Figure 2 · p. 8

4.2 EXPLANATION FAITHFULNESS AND COMPLETENESS

• HARMONIA matches GNAN’s perfect feature selection while reducing response error (0.193 versus 0.435), achieving the highest edge NDCG and precision (0.9425 and 0.8162), and leading node recovery in ranking and signed-effect accuracy.Node recovery reaches 0.9658 NDCG@4 and 0.2992 Signed NRMSE.Source: §4.2 · p. 8, Table 2 · p. 9
• HARMONIA’s centered response curves closely track the four ground-truth linear, quadratic, sinusoidal, and saturating functions.The synthetic graphs use four signal features with feature-specific hop profiles and 16 non-contributing features.Source: §4.2 · p. 8, Figure 3 · p. 8
• UMAP visualizations show that too few experts constrain specialization, while too many fragment feature groups and weaken separation and identifiability.The visualizations show Cora features, with colors indicating expert assignment.Source: §4.2 · p. 8, Figure 4 · p. 9

Interpretable additive models

• GAMs, NAMs, and NBMs model feature-wise responses, with NBMs sharing neural bases to reduce parameter costs; GrAMs extend additive structure to relational data.GNAN uses shortest-path distance for structural influence, while G-NAMRFF uses random Fourier features for parameter-efficient feature modeling.Source: §5 · p. 9
• GMAN and SUPERMAN extend GNAN to sparse temporal or heterogeneous graph settings and allow feature or graph grouping, trading fine-grained interpretability for greater expressivity.HARMONIA instead combines specialized shared neural bases with multi-scale RRWP structural modeling.Source: §5 · p. 9

Post-hoc GNN explainers

• Post-hoc methods explain trained GNNs through auxiliary mechanisms, whereas HARMONIA exposes feature and structural contributions directly through its predictive function.Examples include subgraph and feature identification, local feature attribution, and human-readable rules.Source: §5 · p. 9

Conclusion

• HARMONIA scales with feature dimension and graph size while using specialized bases for feature responses and sparse propagation without dense RRWP tensors.The framework preserves fine-grained explanations.Source: §6 · p. 9

Future work

• The conclusion identifies extending HARMONIA to hypergraphs and improving feature and structural component identifiability and stability as future work.It also proposes principled expert-count selection to avoid overspecified configurations.Source: §6 · p. 9

• Generative AI tools were used only for code debugging and improving manuscript grammar.Source: AI USE STATEMENT · p. 10
• The authors developed the experimental design, analysis, results, and conclusions independently of generative AI assistance.Source: AI USE STATEMENT · p. 10
• The authors retain responsibility for the scientific content and interpretations.Source: AI USE STATEMENT · p. 10

B.1 NOTATION AND TECHNICAL SETUP

• MoNB relaxes a shared basis space by allowing features with different response shapes to use specialized experts, which sparse top-m routing can combine.The comparison isolates the potential benefit of specialization from errors due to a particular basis choice or poor training.Source: §B.1 · p. 16

B.2 IRREDUCIBLE APPROXIMATION ERROR BY A GLOBAL BASIS SPACE

• The optimal B-dimensional global basis spans the top B eigenfunctions, while the remaining spectral tail is irreducible approximation error independent of optimization or training.A rapidly declining spectrum makes global sharing more adequate; substantial variation beyond the first B directions indicates a limitation of the shared space.Source: §B.2 · pp. 17, 20

B.3 FUNCTIONAL SPECIALIZATION AND THE SPECIALIZATION GAP

• Specialization replaces the global spectral tail with group-wise local tails, helping when the full collection is complex but feature groups have simpler functional structure.The specialization gap is the global error minus the best oracle mixture error; a near-zero gap means multiple experts offer little representational benefit.Source: §B.3 · pp. 21–22
• When each of C groups fits a subspace of dimension at most B but the full collection spans more than B directions, mixture error is zero while global-basis error is positive.This yields a positive specialization gap.Source: §B.3 · p. 22
• With B = 2, two groups spanning {1, x} and {sin x, cos x} are individually exact, but their four independent functions cannot fit one two-dimensional global basis.Sharing a low-dimensional representation across the groups therefore creates approximation error.Source: §B.3 · p. 22
• The specialization gain assumes an appropriate hard-routing assignment; assigning a feature to an unsuitable expert increases its approximation error and can reduce the oracle gain.The analysis identifies learned routing errors as an additional cost beyond the representational benefit of specialization.Source: §B.3 · pp. 22, 24
• Hard-routed MoNB is equivalent to the specialized-subspace approximation problem: routing selects a feature’s functional space, and its coefficients choose a function within that space.Thus, the analyzed specialization error applies to the actual MoNB parameterization.Source: §B.3 · pp. 23–24
• Across experts, MoNB can span up to CB functional directions while each feature activates at most B; a sufficiently wide global basis could span the same space.The distinction is how capacity is organized, not an absolute expressivity advantage over an arbitrarily wide global model.Source: §B.3 · p. 24

B.4 ROUTING ERROR AND THE SPECIALIZATION–ROUTING TRADE-OFF

• Specialization beats global sharing exactly when routing penalty Rroute is below specialization gap ΓB,C.The benefit can persist despite imperfect routing when the specialization gap exceeds the routing penalty.Source: §B.4 · p. 25
• Routing quality depends on functional loss, not assignment accuracy: few errors can be costly across dissimilar expert spaces, while similar spaces make confusion cheap.The routing penalty reflects both how often harmful assignments occur and their average approximation cost.Source: §B.4 · p. 26
• Adding experts helps only when their added oracle approximation gain exceeds added routing cost, making C a model-selection parameter.Oracle specialization gain is nondecreasing with the number of experts, but learned-model approximation error need not be.Source: §B.4 · p. 27

B.5 STATISTICAL EFFICIENCY OF SPECIALIZED BASIS SHARING

• Coefficient-estimation error scales as σ^2rk/N for a specialized feature versus σ^2R/N for a dense global basis.This comparison assumes fixed basis functions and routing, least-squares estimation, and the stated model assumptions.Source: §B.5 · p. 28

C IDENTIFIABILITY OF RANDOM-WALK STRUCTURAL EXPLANATIONS

• Hop-level explanations require unique coefficients and stability: distinct coefficient vectors may yield the same or nearly the same structural operator.The analysis separates exact identifiability from conditioning under small perturbations.Source: C IDENTIFIABILITY OF RANDOM-WALK… · p. 29

C.1 NORMALIZATION AND SCALE IDENTIFIABILITY

• Simplex-normalized RRWP operators are row-stochastic, removing global scale ambiguity while leaving hop-level uniqueness unresolved.Normalization makes θt,k interpretable as the fraction of structural influence assigned to walk length t.Source: §C.1 · pp. 30–31

C.2 IDENTIFIABILITY OF RANDOM-WALK HOP COEFFICIENTS

• Hop coefficients are identifiable exactly when I, M, ..., M^(T−1) are linearly independent; for diagonalizable M, the limit equals its number of distinct eigenvalues.Equivalently, T cannot exceed the degree qM of M’s minimal polynomial.Source: §C.2 · pp. 31–33

C.3 APPROXIMATE IDENTIFIABILITY AND SPECTRAL CONDITIONING

• The conditioning constant γT measures how much the structural operator changes under normalized coefficient changes; near zero, recovered hop weights can be unstable.A larger γT limits the change in recovered weights caused by a small operator perturbation.Source: §C.3 · pp. 33–34

D.1 DECAY OF HOP DISTINGUISHABILITY

• For an ergodic reversible random walk, the maximum distinction between consecutive walk lengths decays at least geometrically with nontrivial spectral radius ρ.The decay bounds how much one additional walk step can change the propagated structural signal.Source: §D.1 · pp. 35–36

D EFFECTIVE STRUCTURAL HORIZON OF RANDOM-WALK EXPLANATIONS

• The effective structural horizon depends on graph mixing: when ρ is small, hop distinctions fade sooner; when ρ is near one, they persist to larger t.Beyond this graph-dependent range, separating individual hop effects becomes increasingly unstable and less informative.Source: D EFFECTIVE STRUCTURAL HORIZON OF… · p. 35, §D.1 · p. 37

D.2 EFFECTIVE STRUCTURAL HORIZON AND SPECTRAL DEPENDENCE

• The effective horizon Hε is a tolerance-dependent scale beyond which one more random-walk step changes the propagation operator by at most ε.Before the horizon, at least some consecutive walk lengths differ by more than ε.Source: §D.2 · p. 38
• For an ergodic reversible walk, the horizon grows as the nontrivial spectral radius ρ approaches one, so slowly mixing graphs retain distinguishable hop behavior longer.A small spectral gap corresponds to a longer horizon; a large gap leads to faster loss of hop-specific information.Source: §D.2 · pp. 38–39
• At the same tolerance, a graph with ρ = 0.9 can retain useful long-range walk distinctions after a graph with ρ = 0.5 has passed its useful range.The text gives roughly five steps as the bound for ρ = 0.5 and notes that T = 10 may exceed that graph’s informative range while remaining within the other’s.Source: §D.2 · p. 39
• Choosing T far beyond Hε includes neighboring walk operators that differ only weakly, making late-hop weight interpretations increasingly sensitive to finite data and optimization noise.Choosing T within the graph’s effective range focuses the decomposition on meaningfully distinguishable walk lengths.Source: §D.2 · p. 39

E INTERPRETABILITY AND EXPRESSIVITY OF NORMALIZED RRWP

• Simplex-normalized RRWP makes hop weights interpretable as nonnegative shares, but favors diffusion-like aggregation and excludes signed filters such as I−M.This is an expressivity trade-off: normalization makes the hop decomposition easier to interpret while restricting the structural operators represented.Source: E INTERPRETABILITY AND EXPRESSIVITY OF… · p. 39

E.1 SPECTRAL VIEW OF SIMPLEX-NORMALIZED RRWP

• For a lazy random walk, normalized RRWP has a nondecreasing spectral response on [0, 1], preserves the mode at λ = 1, and tends to attenuate smaller-eigenvalue modes.More generally, the response is bounded in magnitude by one for λ in [−1, 1].Source: §E.1 · pp. 40–41

E.2 STABILITY AND THE INTERPRETABILITY–EXPRESSIVITY TRADE-OFF

• Simplex normalization makes propagation non-expansive in ℓ∞, but excludes signed combinations such as I−M and therefore represents a narrower family than unrestricted polynomial filters.The restriction preserves interpretable nonnegative hop shares and prevents the structural operator from amplifying the maximum node-wise perturbation.Source: §E.2 · pp. 42–43

E.3 EXACT ADDITIVE EXPLANATIONS

• HARMONIA’s feature–node contributions sum exactly to the class logit, and regrouping them yields feature-, source-node-, and hop-level explanations from the predictor’s own terms.The exact completeness applies to pre-activation logits; sigmoid or softmax probabilities do not retain the same additive decomposition.Source: §E.3 · pp. 43–44

F SPARSE RRWP AGGREGATION

• Sparse RRWP Aggregation replaces dense pairwise RRWP materialization with repeated sparse propagation, avoiding the pairwise-node cost while computing the same linear operator.The dense descriptor collection requires O(Tn^2) memory; SRA evaluates the aggregation through propagation with the transition matrix.Source: F SPARSE RRWP AGGREGATION · pp. 44–45

F.1 EXACT EQUIVALENCE OF DENSE RRWP AND SPARSE AGGREGATION

• SRA and dense RRWP produce exactly the same representation, so the reformulation changes computation rather than learned hop coefficients, structural contributions, or predictions.The equivalence follows by exchanging finite sums in the linear RRWP aggregation.Source: §F.1 · pp. 45–46

F.2 TIME AND MEMORY COMPLEXITY

• Dense RRWP materialization requires O(Tn2) memory and direct application requires O(Tn2d) arithmetic, whereas SRA propagates feature-response matrices sparsely.The dense higher-order walk operators may become dense even when the transition matrix is sparse.Source: §F.2 · pp. 46–47
• SRA computes the same linear operator in O(T(EM + n)d) time, simplifying to O(T|E|d) when EM = Θ(|E|).For fixed T and d, structural time and memory are linear in n whenever EM = O(n); retaining all propagated states increases memory to O(EM + Tnd), plus coefficient storage.Source: §F.2 · p. 47

G.1 HARDWARE SETTINGS / G.2 DATASETS AND EVALUATION SPLITS

• Node-classification results use five independent seeds; graph classification uses 10-fold cross-validation repeated with three seeds, and the evaluation covers twelve real-world datasets.Experiments used a Linux server with a single RTX 4090 GPU with 24 GB VRAM; OOM denotes failure to complete within that budget.Source: §G.2 · p. 48, §G.1 · p. 48

G.5 MIXTURE-OF-NEURAL-BASIS ROUTING

• HARMONIA routes each feature to selected neural-basis experts, so each feature estimates B basis coefficients while the full expert dictionary contains up to CB basis functions.Feature-identity-based routing keeps each response globally inspectable as a univariate function.Source: §G.5 · pp. 49–50

G.6 PREDICTION HEADS AND GRAPH READOUT

• Node logits retain feature-wise contributions, and sum readout preserves the additive decomposition for graph logits up to the intercept.The task-specific loss is not fully documented in the manuscript, and implementation settings are represented by placeholders rather than reported values.Source: §G.6 · p. 50, §G.7 · p. 50

G.8 COMPUTE AND REPRODUCIBILITY DETAILS

• The large-graph results establish completion within the reported hardware budget, not controlled peak-memory or training-time measurements.The supplied passage ends before specifying the further boundary of the scalability evidence.Source: §G.8 · p. 50
• SRA recursively propagates the feature-response matrix and can stream its forward pass, but ordinary reverse-mode training may retain O(TNd) activations.The streaming memory bound excludes feature-network activations and optimizer state.Source: §G.8 · p. 51

G.9 GRAPH-CLASSIFICATION RESULTS

• HARMONIA has the highest reported mean accuracy among interpretable methods on all six graph-classification datasets, while black-box comparisons are mixed and are not significance tests.It is close to the strongest black-box result on PROTEINS and Mutagenicity, with larger gaps on several other datasets.Source: §G.9 · pp. 51–52

G.10 SYNTHETIC DATA CONSTRUCTION AND EVALUATION PROTOCOL

• Each synthetic data seed defines a 64-node undirected graph with 78 or 79 edges, 20 raw features, and four informative columns; instances share topology within each seed.There are 384 feature/label draws per graph, split into 256 training, 64 validation, and 64 test instances; informative feature locations vary across seeds.Source: §G.10 · p. 52
• The generator combines four response templates with four feature-specific structural profiles, assigned by a predeclared permutation, so structural recovery is evaluated separately for each feature.The functions, informative columns, and hop profiles are hidden during fitting.Source: §G.10 · p. 52, Table 9 · p. 53
• Each data seed is paired with three model seeds, yielding 15 runs per method and 128 deterministic test-node targets per run for matched explanation comparisons.Reported feature–hop means use the 15 runs; edge summaries average targets within runs and then model seeds within each data seed.Source: §G.10 · pp. 52–53

G.11 FEATURE–HOP RECOVERY: EXPLANATIONS, METRICS, AND RESULTS

• HARMONIA and GNAN both reach Precision@4 of 1.000, while HARMONIA’s NDCG@4 is 1.000 versus 0.996 and its Effective NRMSE is 0.193 versus 0.435.HARMONIA’s top-hop accuracy is 0.950 with coverage 1.000; its nonzero Hop L1 of 0.468 indicates that identifying the dominant hop does not recover the full profile exactly.Source: §G.11 · p. 54
• On shared GCN checkpoints, GNNExplainer exceeds the feature-only UCExplainer adapter on all three feature-ranking metrics; the adapter changes predictions for 1,588 of 1,920 targets but changes all 20 features on average.The results distinguish counterfactual validity from feature sparsity.Source: §G.11 · p. 54

G.12 INSTANCE-LEVEL STRUCTURAL EXPLANATIONS

• HARMONIA achieves the highest mean on all three edge-recovery metrics: Spearman 0.7822, NDCG 0.9425, and precision 0.8162.G-NAMRFF is the strongest baseline, with corresponding values of 0.6006, 0.7589, and 0.6430; these aggregate comparisons do not establish statistical significance.Source: Table 11 · p. 55, §G.12 · p. 56
• The recovery benchmark tests a known additive random-walk generator matching HARMONIA’s structural representation, and does not establish recovery under other graph generators.It also does not isolate MoNB’s contribution or establish coefficient identifiability from observed samples.Source: §G.12 · p. 56

G.13 SOURCE-NODE CONTRIBUTION RECOVERY

• HARMONIA achieved the best mean on all five source-node recovery metrics, including NDCG@4 of 0.9658, signed NRMSE of 0.2992, and top-four sign agreement of 0.9792.These results concern training-mean feature-replacement effects under the controlled additive random-walk generator; they do not establish a unique attribution under other baselines or statistical significance across three datasets.Source: §G.13 · pp. 58–59, Table 13 · p. 58
• For the frozen GCN, GNNExplainer raised Spearman to 0.8893 from 0.8327 but lowered Precision@4 to 0.5312 from 0.6771 and NDCG@4 to 0.6914 from 0.8280.The comparison is between GNNExplainer and direct counterfactual effects on the same cohort.Source: §G.13 · p. 59
• UCExplainer achieved counterfactual validity of 0.8993, but its ranking had Spearman −0.0264 and NDCG@4 of 0.2076, while changing all 64 nodes.The reported node changes describe the full explanation before selecting its top-four nodes.Source: §G.13 · p. 59, Table 14 · p. 59

G.14 ADDITIONAL UNIVARIATE RESPONSE RECOVERY

• Mean curve errors were 0.023 for linear, 0.033 for sigmoid, 0.034 for hyperbolic tangent and arctangent, and 0.364 for sine responses.The twelve-function study is separate from the four-informative-feature benchmark; Figure 5 reports the curve NRMSE summaries.Source: §G.14 · p. 59, Figure 5 · p. 60

H.1 FUNCTIONAL SPECIALIZATION AND HETEROGENEITY

• Global NBM and MoNB differed by ∆=0.0019 at low heterogeneity, while MoNB’s approximation gain rose to ∆=0.0366 and 0.0418 at medium and high heterogeneity.Independent per-feature networks remain more flexible, while MoNB retains parameter sharing and features concentrate gating weights on different experts.Source: §H.1 · pp. 60–61, Figure 6 · p. 61

H.2 EXPERT COUNT AND BASIS CAPACITY

• At B=16, increasing experts C from 1 to 5 reduced centered NRMSE from 0.086 to 0.052 under medium heterogeneity and from 0.080 to 0.041 under high heterogeneity, but gains saturated beyond C=5.At C=5 under high heterogeneity, increasing B from 4 to 8 to 16 reduced error from 0.056 to 0.045 to 0.041; C controls available functional spaces, while B controls each expert’s capacity.Source: §H.2 · pp. 61–62, Figure 7 · p. 62

H.3 ROUTING QUALITY AND SPECIALIZATION GAIN

• Specialization remains beneficial only while routing penalty Rroute is below gain ΓB,C; routing quality depends on functional cost, not merely the fraction of incorrect assignments.Across controlled settings, Eroute−Eglobal agrees with Rroute−ΓB,C up to numerical precision, and different corruption policies follow nearly the same damage relationship when measured by Rroute.Source: §H.3 · pp. 62–63

H.4 HOP-LEVEL EXPLANATION RELIABILITY

• Under ϵ=0.01, as γT decreased in representative ill-conditioned settings, operator error remained 0.007–0.018 while hop-allocation coefficient error rose to approximately 0.26–0.54, showing operator recovery did not ensure reliable allocation.Source: §H.4 · p. 64
• Even when operator recovery is accurate and unique, hop-level explanations may be unstable: smaller γT permits larger deviations in recovered hop coefficients.Figure 9 uses only non-singular settings, so the observed instability occurs before exact non-identifiability.Source: §H.4 · p. 64

H.5 EMPIRICAL EFFECTIVE STRUCTURAL HORIZON

• Graph mixing determines the effective structural horizon: complete graphs lose hop distinction rapidly, rings preserve it longer, and block-structured graphs fall between.Every measured Hϵ is below its spectral upper bound bHϵ, which is conservative rather than an exact predictor; truncation length T exposes depths, while Hϵ estimates how many remain distinguishable at the chosen resolution.Source: §H.5 · pp. 65–66, Figure 10 · p. 65

H.6 RRWP NORMALIZATION AND LINEAR–NONLINEAR ABLATIONS

• The DNN RRWP topology function learns nonlinear relationships among channels, while the linear variant enables sparse computation and explicit walk-length attribution.The nonlinear channel mixing increases expressiveness but costs scalability and prevents decomposition into individual t-hop terms.Source: §H.6 · p. 66

H.7.1 REMOVE-AND-RETRAIN EVALUATION

• Top-ranked feature removal causes a much larger mean performance decrease than bottom-ranked or random removal, indicating predictive relevance but not causality.The small positive mean after bottom-feature removal does not establish a statistically significant improvement.Source: §H.7.1 · p. 66, Figure 11 · p. 67

H.7.2 IMPORTANCE STABILITY ACROSS SEEDS

• The five-seed importance scatter is descriptive, not an independent reliability estimate; no numerical Spearman correlation can be inferred from it.A rank-based summary requires run outputs, and constant importance vectors need explicit handling.Source: §H.7.2 · p. 67

H.8 PARAMETER COMPLEXITY

• At D = 5,000, HARMONIA requires only on the order of 10^5 learnable parameters, versus 10^7 for GNAN and G-NAMRFF, a gap of approximately one to two orders of magnitude.Figure 13 reports the comparison in the synthetic multi-class setting with K = 47 output classes; HARMONIA shares nonlinear expert parameters across features, adding feature-specific parameters as D increases.Source: §H.8 · p. 68, Figure 13 · p. 68

H.9 LEARNED FEATURE SHAPE FUNCTIONS

• On Mutagenicity, HARMONIA learns distinct feature-specific response functions with different signs, magnitudes, and nonlinear trends, rather than a common shape.These curves describe the trained predictive model conditional on its learned representation and observed data support; they do not establish causal chemical effects.Source: §H.9 · pp. 68–69, Figure 14 · p. 69

H.10 ABLATION STUDY

• Without normalization, feature-response and structural coefficients have a multiplicative scale ambiguity; the normalized operator is a convex combination of row-stochastic walk operators.The resulting operator is non-expansive in ℓ∞.Source: §H.10 · p. 69
• The ablation holds the MoNB model, routing, walk orders, optimization, and data splits fixed, but removes both simplex nonnegativity and unit-sum constraints.Source: §H.10 · p. 69
• Feature-level Effective NRMSE decreases from 0.2279 to 0.2049 with simplex normalization, a 10.1% relative reduction.Source: §H.10 · p. 70
• Signed source-node contribution NRMSE decreases from 0.3252 to 0.3125, while source node Precision@4 rises from 0.8409 to 0.8509.Source: §H.10 · p. 70
• Hop L1 falls from 0.5491 to 0.4409 with simplex normalization, a 19.7% relative reduction.The comparison uses recovered and oracle hop profiles for four informative features.Source: §H.10 · p. 70, Figure 15 · p. 71
• Edge-ranking Spearman correlation increases from 0.7556 to 0.7823, and Precision@20% rises from 0.8053 to 0.8329.Figure 15b reports consistency across three matched data seeds.Source: §H.10 · p. 70, Figure 15 · p. 71