Higher-Order Feature Attribution

• Higher-order feature attribution is a set of techniques that quantify the collective impact of feature groups on a model's output beyond individual contributions.
• It distinguishes main effects from interaction effects by using methods like Integrated Gradients, Möbius transforms, and cross-derivative approaches.
• These methods improve model interpretability by isolating genuine multi-feature interactions, offering scalable and actionable insights into complex predictive models.

Higher-order feature attribution denotes a family of post-training analysis methods that assess how groups of input features, rather than only individual coordinates, contribute to a model output. The topic arises because the interpretation of standard feature attributions is straightforward when features act independently, but becomes less direct when the predictive model involves interactions such as multiplicative relationships or joint feature contributions. In that regime, first-order attributions satisfy completeness but do not, by themselves, distinguish pure main effects from joint, multiplicative, or higher-order interactions. Recent work has therefore developed several formalizations of higher-order attribution, including path-integral constructions built on Integrated Gradients, Möbius-transform and Harsanyi-dividend decompositions, cross-derivative interaction scores, submodular set-function models, and formal explanation-based scores (Butler et al., 7 Oct 2025).

1. Problem setting and conceptual scope

Standard feature-attribution methods assign to each input feature $x_i$ a scalar score $a_i$ intended to measure how much $x_i$ moved the prediction away from a chosen baseline. This is sufficient for additive models, but it is incomplete when the effect of one feature depends on another. The canonical motivating example is

$f(x_1,x_2)=3\,x_1x_2,$

for which the marginal effect of $x_1$ depends on the value of $x_2$. In such cases, higher-order attribution seeks to isolate and quantify contributions associated to subsets of features rather than to singleton features alone (Butler et al., 7 Oct 2025).

Different research traditions formalize this objective in different spaces. In the Integrated Gradients line, higher-order terms are defined by iterating IG operators and using higher mixed partials. In the set-function line, one studies the contribution of subsets $S\subseteq D$ through ablation-based function values $f(S)$, Möbius transforms, or Harsanyi dividends. In cross-derivative approaches, the interaction object is the local mixed partial $$\partial^{\,k}F/\partial x_{i_1}\cdots\partial x_{i_k}$$. In formal XAI, higher-order attribution appears as a score assigned to subsets $S$ according to how often they occur inside minimal formal explanations. A plausible implication is that “higher-order feature attribution” is less a single algorithm than a common explanatory objective realized by several mathematically distinct decompositions (Jiang et al., 2023).

A recurring distinction is between overall effect and pure interaction effect. First-order scores often remain useful as overall main-effect scores, but higher-order terms are introduced to identify what is not already explained by lower-order subsets. This distinction is explicit in Möbius-transform methods, in Harsanyi-dividend aggregation, and in the cumulant analogy developed for higher-order IG (Jiang et al., 2023).

2. Integrated-Gradients-based higher-order attribution

A general theory based on Integrated Gradients fixes a model $a_i$0 and a baseline $a_i$1. The first-order IG attribution to feature $a_i$2 is

$a_i$3

The $a_i$4th-order attribution term for the feature-index tuple $a_i$5 is then defined by composing IG operators in analogy with higher derivatives: $a_i$6 For $a_i$7 this recovers usual IG; for $a_i$8 it yields the Integrated Hessians formulas of Janizek et al. (Butler et al., 7 Oct 2025).

This construction supports several structural results. For second order, $a_i$9, and each attribution operator $x_i$0 is linear. A marginalization identity states that, for each $x_i$1,

$x_i$2

Completeness extends across orders: $x_i$3 and more generally

$x_i$4

If $x_i$5, then all mixed attributions vanish: $x_i$6 for $x_i$7. The stated proof ingredients are repeated application of the fundamental theorem of calculus along the baseline-to-input path, the Fubini theorem, and Clairaut’s theorem (Butler et al., 7 Oct 2025).

The simple bilinear example illustrates how the hierarchy behaves. For $x_i$8 with baseline $x_i$9, the first-order attributions satisfy $f(x_1,x_2)=3\,x_1x_2,$0 and $f(x_1,x_2)=3\,x_1x_2,$1, while the second-order term is $f(x_1,x_2)=3\,x_1x_2,$2. One checks that $f(x_1,x_2)=3\,x_1x_2,$3 and $f(x_1,x_2)=3\,x_1x_2,$4. Here the nonzero $f(x_1,x_2)=3\,x_1x_2,$5 is the explicit interaction signal. In the synthetic polynomial

$f(x_1,x_2)=3\,x_1x_2,$6

a third-order run recovers nonzero third-order scores exactly on triples $f(x_1,x_2)=3\,x_1x_2,$7 and $f(x_1,x_2)=3\,x_1x_2,$8, while second order picks up the $f(x_1,x_2)=3\,x_1x_2,$9 term but splits the triple interactions into pairwise footprints (Butler et al., 7 Oct 2025).

3. Set-function, Möbius, and formal explanation formulations

A second major line of work treats a model as a set function over feature subsets. Let $x_1$0, and write $x_1$1 for the model output on the counterfactual input in which only features in $x_1$2 are present. The Möbius transform

$x_1$3

isolates the pure interaction effect of $x_1$4, and inclusion–exclusion gives

$x_1$5

A Weighted Möbius Score is then defined by

$x_1$6

with the faithfulness constraint $x_1$7 whenever $x_1$8. In this framework, Shapley values, the Grabisch–Roubens Shapley interaction index, Shapley–Taylor interaction indices, pure indirect effect, total indirect effect, and ArchAttribute all arise as special choices of the weights (Jiang et al., 2023).

The conceptual role of the Möbius transform is close to the role played by Harsanyi dividends in cooperative games. H-Sets uses the dividend

$x_1$9

where $x_2$0, to isolate the unique contribution of exactly set $x_2$1, removing lower-order overlaps. The paper states that $x_2$2 captures precisely the interaction effect of $x_2$3 not accounted for by any proper subset, yields an orthogonal decomposition over all $x_2$4 coalitions, and is the unique additive Möbius transform (Mehrotra et al., 23 Apr 2026).

Formal Feature Attribution provides a logically different higher-order semantics. For a classifier $x_2$5, a formal explanation of why $x_2$6 is a subset-minimal sufficient set $x_2$7 such that every input agreeing with $x_2$8 on $x_2$9 is classified as $S\subseteq D$0. If $S\subseteq D$1 denotes the set of all such formal explanations, then for any subset $S\subseteq D$2,

$S\subseteq D$3

Thus higher-order attribution becomes the fraction of all minimal explanations containing $S\subseteq D$4. A weighted version $S\subseteq D$5 additionally favors short explanations by summing $S\subseteq D$6 over explanations containing the feature of interest (Yu et al., 2023).

These formulations differ in primitive object—path integrals, counterfactual subset values, or logical explanations—but they converge on a common aim: to separate pure interaction structure from lower-order accumulation. This suggests that higher-order attribution is fundamentally a decomposition problem, with the choice of decomposition determined by the semantics of “feature presence,” “interaction,” and “explanation.”

4. Statistical, differential, and topological interpretations

The IG-based theory makes explicit connections to classical statistics. The $S\subseteq D$7th-order IG term is compared to ANOVA and Sobol indices because it isolates the contribution of exactly $S\subseteq D$8 features to the output change. It is also compared to the $S\subseteq D$9th-order joint cumulant, which measures pure $f(S)$0-way dependence among random variables. Under that analogy, a large second-order attribution $f(S)$1 indicates a strong two-feature interaction just as a large covariance signals linear dependence (Butler et al., 7 Oct 2025).

Cross-derivative methods make the differential viewpoint explicit. For a differentiable $f(S)$2 and subset $f(S)$3, the local $f(S)$4-way interaction effect is

$f(S)$5

When $f(S)$6, this is the Hessian cross term; for $f(S)$7, it is the natural generalization to higher mixed partials. The T-NID procedure aggregates these quantities over a small representative set of samples rather than over an entire dataset, and Taylor-CAM extends the same principle to visual explanations by replacing first-order Grad-CAM weights with higher-order cross derivatives (Lerman et al., 2020).

The topological interpretation advanced in the IG-based theory organizes attributions by the combinatorial topology of the complete simplex on $f(S)$8 vertices. Nodes ($f(S)$9-simplices) carry first-order scores $$\partial^{\,k}F/\partial x_{i_1}\cdots\partial x_{i_k}$$0, edges ($$\partial^{\,k}F/\partial x_{i_1}\cdots\partial x_{i_k}$$1-simplices) carry second-order scores $$\partial^{\,k}F/\partial x_{i_1}\cdots\partial x_{i_k}$$2, triangles ($$\partial^{\,k}F/\partial x_{i_1}\cdots\partial x_{i_k}$$3-simplices) carry third-order scores $$\partial^{\,k}F/\partial x_{i_1}\cdots\partial x_{i_k}$$4, and so on. The collection of all $$\partial^{\,k}F/\partial x_{i_1}\cdots\partial x_{i_k}$$5th-order attributions can therefore be viewed as a signal supported on the set of $$\partial^{\,k}F/\partial x_{i_1}\cdots\partial x_{i_k}$$6-simplices of the full simplex. In this representation, completeness and marginalization become boundary-operator relationships in the simplicial complex, which opens the door to topological graph signal processing tools such as Hodge Laplacians and simplicial Fourier bases for analysis and visualization (Butler et al., 7 Oct 2025).

A plausible implication is that higher-order attribution admits not only scalar ranking interpretations but also structured-signal interpretations. In that view, a model explanation is not merely a list of scores but a hierarchy of signals over a simplicial object, enabling localization of interaction patterns by order and by combinatorial support.

5. Algorithms, scalability, and computational regimes

The principal computational difficulty is combinatorial growth with interaction order. For higher-order IG, a naive $$\partial^{\,k}F/\partial x_{i_1}\cdots\partial x_{i_k}$$7th-order computation requires a $$\partial^{\,k}F/\partial x_{i_1}\cdots\partial x_{i_k}$$8-fold integral over $$\partial^{\,k}F/\partial x_{i_1}\cdots\partial x_{i_k}$$9; with an $S$0-point quadrature in each dimension, the cost is $S$1. The paper states that second order is usually feasible, third order is on the edge, and beyond that one must be selective. Suggested strategies are low-order Monte Carlo, path-bundling, and sparsity thresholding, where first-order scores are estimated first and higher-order scores are computed only for the most influential features (Butler et al., 7 Oct 2025).

Within the Weighted Möbius framework, computing $S$2 for all $S$3 requires

$S$4

model evaluations in the naïve case. The paper describes exact enumeration as feasible for $S$5 up to $S$6 when $S$7, a recursive dynamic-programming relation $S$8, sampling and Monte Carlo approximations for Shapley-type scores, and greedy subset selection or orthogonal matching pursuit to identify a small number of high-magnitude interactions without full enumeration (Jiang et al., 2023).

Submodular formulations address a different computational problem: how to learn a context-dependent scoring function over subsets. A deep submodular set function $S$9 is learned so that its continuous extension scores given attribution heatmaps, and final feature attribution is defined by the marginal gain

$a_i$00

during a greedy traversal. The stated effect is diminishing returns: once correlated or redundant features are already in $a_i$01, adding one of them yields a smaller marginal gain, thereby discounting redundancy and increasing specificity (Manupriya et al., 2021).

In graphs, higher-order attribution has been developed for GNN-LRP. Subgraph relevance $a_i$02 is defined as a sum over all walks that stay inside a node subset $a_i$03, which naively costs $a_i$04. The message-passing reformulation sGNN-LRP replaces the exponential walk sum with a dynamic layer-wise update of pooled relevance messages and computes the same subgraph score in $a_i$05. In practice, with the forward-hook implementation, the cost is the same as two standard GNN forward/backward passes, with only constant overhead for masking (Xiong et al., 21 May 2026).

Image-domain discovery of higher-order sets introduces yet another regime. H-Sets first detects locally interacting pairs through input Hessians $a_i$06, recursively merges them into sets using a threshold $a_i$07 and maximum size $a_i$08, and then attributes the resulting sets using IDG-Vis, a set-level extension of Integrated Directional Gradients. With $a_i$09 sets and $a_i$10, the total detection cost is $a_i$11 Hessian terms, described as tractable with GPU-accelerated second derivatives (Mehrotra et al., 23 Apr 2026).

6. Empirical behavior, applications, and interpretive issues

Empirical studies span synthetic functions, NLP, images, graphs, and formal verification settings. In sentiment analysis with BERT on SST-2, the Weighted Möbius framework treats each token’s hidden representation as a mediator, uses $a_i$12 for first-order attribution and $a_i$13 for pairwise interactions, and reports that lower Transformer layers exhibit stronger pairwise interaction signals, while higher layers concentrate attribution on single tokens, especially CLS. The example “not”$a_i$14“bad” has a positive second-order score despite each token alone being negative (Jiang et al., 2023).

In chain-of-thought prompt engineering on GPT-3.5 for last-letter concatenation, sentences or phrases are treated as features. The reported findings distinguish methods sharply: Möbius scores capture saturation through negative pairwise scores when both individual sentences are positive; TIE simply sums all interactions and is uninformative; Shapley-value and Shapley-interaction scores split interaction credit evenly across members, masking which specific combination is most critical; ArchAttribute recovers exactly the model’s accuracy lift from including particular sentences. Phrase-level scores further suggest that CoT prompts often work by highlighting key phrases rather than by the human-like logical structure (Jiang et al., 2023).

In image classification, H-Sets combines Hessian-based interaction detection with IDG-Vis and Harsanyi dividends. Across VGG16, ResNet101, DenseNet121, and MobileNetV3 on ImageNet and CUB, it is evaluated against IG, GradCAM, CAFO/CASO, Archipelago, and MOXI using sparsity via Gini index and faithfulness via ROAD_AOPC. Representative ImageNet results list, for example, ResNet sparsity values IG $a_i$15, Arch $a_i$16, MoXI $a_i$17, H-Sets $a_i$18, and faithfulness values IG $a_i$19, Arch $a_i$20, MoXI $a_i$21, H-Sets $a_i$22. The paper states that H-Sets yields saliency maps that are both sparser and more faithful, and that qualitative maps emphasize compact, semantically coherent object parts rather than diffuse gradients or coarse superpixels (Mehrotra et al., 23 Apr 2026).

Cross-derivative approaches also report quantitative advantages. On synthetic interaction benchmarks, T-NID with cross derivatives attains average pairwise AUC $a_i$23, compared with NID $a_i$24 and GLIDER $a_i$25. On Sort-Of-CLEVR relational questions, Taylor-CAM reports top-1 region-pair accuracy $a_i$26 for the “furthest” question and $a_i$27 for the “same shape” question, outperforming Grad-CAM* and GLIDER. A user study with $a_i$28 further reports higher object-guess and question-guess accuracy for Taylor-CAM than for Grad-CAM and GLIDER (Lerman et al., 2020).

For GNN explanations, sGNN-LRP addresses the scalability bottleneck directly. On BA-2motif with GIN depths $a_i$29, naïve GNN-LRP takes $a_i$30 ms / $a_i$31 ms / $a_i$32 ms, while sGNN-LRP takes $a_i$33 ms / $a_i$34 ms / $a_i$35 ms, corresponding to $a_i$36 speedup. On real datasets, AUAC and AUPC are used to assess activation and pruning behavior, and tuning the discount factor $a_i$37 often improves AUAC/AUPC over $a_i$38 (Xiong et al., 21 May 2026).

Formal Feature Attribution highlights a different interpretive issue: soundness. The framework is presented as an alternative to LIME and SHAP that avoids explanation unsoundness and out-of-distribution sampling, but exact FFA is expensive, with decision problems at the second level of the polynomial hierarchy. The paper therefore proposes an anytime approximation based on enumeration of formal explanations, and reports that even after $a_i$39 seconds of enumeration, approximate FFA outperforms LIME and SHAP in both error and ranking metrics on tabular and image settings (Yu et al., 2023).

Across these literatures, a recurring misconception is that higher-order attribution is synonymous with pairwise heatmaps. The surveyed work instead treats order as a hierarchy: edges are only one level of a broader decomposition that can extend to triples, larger subsets, subgraphs, or minimal explanation sets. Another misconception is that interaction attribution is necessarily incompatible with attribution axioms. The record is mixed: some image methods are described as coarse or as failing core interpretability axioms, whereas H-Sets states that IDG-Vis plus Harsanyi dividends satisfies non-negativity, normality, monotonicity, superadditivity, approximate completeness, sensitivity, implementation invariance, linearity, and symmetry-preserving (Mehrotra et al., 23 Apr 2026).

Higher-order feature attribution therefore occupies a junction of explainable AI, statistics, cooperative game theory, formal methods, and topological signal processing. Its central problem is stable across formulations: decomposing prediction changes into main effects and genuine multi-feature interactions. What varies is the semantics of interaction, the representation space in which it is measured, and the computational compromises required to expose it at useful scales.