Exact Functional ANOVA Decomposition for Categorical Inputs

Summary

The paper proposes a closed-form generalized functional ANOVA decomposition for models with categorical inputs, filling a well-defined gap: prior work (Hoeffding, Stone, Hooker, Chastaing et al., Il Idrissi et al.) had established existence and uniqueness of the decomposition under regularity conditions but left practitioners either with sampling-based approximations or with algorithms tied to tree ensembles under full-support density assumptions. The conceptual move is to extend the Walsh–Hadamard/parity basis of Boolean analysis to arbitrary categorical alphabets via a "signed inverse-likelihood" family $\phi^{(A\mathbf{z})}$, and to show that this family (i) spans $L^2$ of any categorical distribution, (ii) satisfies Hooker-style hierarchical orthogonality by construction, and (iii) yields a linear system whose Gram matrix is block-diagonal under independence — collapsing to the standard orthogonal ANOVA and to SHAP-via-Harsanyi-dividends in that regime. The decomposition is then computed by a greedy rank-completion procedure (Algorithm 1) that adaptively picks a basis over the observed support, naturally handling the sparse tabular regime $|X| \ll |E|$ where density-based methods fail. This is a real theoretical advance that unifies three previously separate strands (SHAP under independence, Fourier/parity analysis, dependent functional ANOVA), which is the main reason to take the contribution seriously. Executed replication corroborates the parts of the empirical story that fit in a short compute window: the analytical decomposition is exact, the independent-case KernelSHAP tables reproduce cell-for-cell, and the closed-form vs TreeHFD orthogonality gap (~$10^{-16}$ vs ~$10^{-3}$) is confirmed. The central conceptual limitation, which the authors flag but do not resolve, is that in the very regime the method targets, the decomposition is non-unique: the greedy basis $S_r$ depends on enumeration order, and so do the coefficients, local Shapley values, and reported rankings; no sensitivity study accompanies the identifiability discussion. The paper otherwise relies on independent-case sanity checks, a synthetic Appendix C benchmark, and a POKER HAND comparison against attention rollout, which is a weak baseline. The methodology is sound but the interpretive claims currently outrun the empirical evidence.

Strengths

• Conceptual contribution. The extension of Boolean parity functions to arbitrary categorical distributions via the signed-inverse-likelihood family $\phi^{(A\mathbf{z})}$ is elegant and directly yields hierarchical orthogonality without a case analysis on the dependence structure. Theorem 3.2 plus Corollaries 3.4/3.5 give a clean picture: existence in general, uniqueness under full support, and recovery of both classical orthogonal ANOVA and standard SHAP under independence.
• Unifying framing. The paper stitches together three previously separate lines — Hooker/Chastaing generalized ANOVA, Boolean Fourier analysis (O'Donnell, Gorji et al.), and Shapley-via-Harsanyi-dividends — into a single decomposition, and makes the SHAP correspondence a special case of a more general local-attribution rule.
• Sparse-regime handling. Working on the observed support $X$ rather than the full hypergrid $E$, plus the block-wise vectorization of the Gram system, gives a genuinely usable algorithm for datasets like MUSHROOMS ($|E| \approx 10^{14}$, $|X| = 8124$) where density-based methods such as TreeHFD do not apply.
• Independent-case sanity checks reproduce cleanly. Executed replication reproduces the CAR EVALUATION ISE table to two decimals across 24 cells, the NURSERY ISE table exactly across all 40 cells, and the analytical Table 1 norms to three decimals with $R^2 = 1$ and MSE $\sim 10^{-21}$ — showing the framework really does behave as advertised in the independent / full-support setting.
• TreeHFD orthogonality gap reproduces. Appendix C's headline advantage — closed-form orthogonality at machine precision ($\sim 10^{-19}$ to $10^{-16}$) versus TreeHFD at $10^{-5}$ to $10^{-3}$ — is confirmed by direct execution on the four completed dimensions.

Weaknesses

• Non-uniqueness in the sparse regime is a first-order limitation, not a footnote. In the setting the paper motivates most strongly, the greedy basis $S_r$ depends on the canonical enumeration order, and so do the coefficients, local Shapley values, and global feature-importance rankings. The paper acknowledges this ("basis support $S_r$ is not unique") but never quantifies how much the reported explanations move under a different order. Without a sensitivity study, the sharp POKER HAND "ranks $\gg$ suits" picture and the MUSHROOMS top-3 could partly be artefacts of the enumeration.
• Attention rollout is a weak baseline for POKER HAND. Attention rollout is not a Shapley-style attribution; beating it on ranks-vs-suits does not stress the paper's methodological claim. The natural baselines here are Interventional (Marginal) SHAP, DeepSHAP, or KernelSHAP with a background of poker hands. The sparse-regime experiment currently has no Shapley-family comparator at all.
• Overclaims relative to the theory. The abstract's "without any assumption" hides that the closed-form is unique only under full support (Corollary 3.4). The "decay of information per set cardinality" statement is not implied by hierarchical orthogonality (it forbids leakage upward, not decay), and the paper's own POKER main-effects row ($R^2 \approx 0.003$) is a counterexample. The KernelSHAP "must coincide" wording clashes with the reported ISE values. The main-effects "maximize information" phrasing has no accompanying optimization problem.
• Explanations are evaluated only on the training support. All experiments compute attributions on the same support used to build the decomposition. It is unclear what the framework says about attributions for an unseen categorical configuration; the uniform residual-splitting rule sends part of the attribution to features that never appear in $S_{r_\mathrm{low}}$, purely by convention.
• Residual-splitting convention affects rankings. In the POKER HAND setup ($R^2 = 0.79$ at $r_\mathrm{low} = 5000$), roughly a fifth of the target's variance sits in $f_\mathrm{res}$, which is then distributed uniformly across the 10 features. On DOTA2 the residual is 60% of the variance. The paper does not check that its interpretive conclusions are robust to this choice or report per-feature attribution error under low-rank truncation.
• Notation errors and rendering issues in tables and proofs. Table 5's $|E|$ column ($525, 342, 3115$) is inconsistent with the surrounding $r \ll |E|$ narrative — the intended $52^5, 3^{42}, 2^{113}$ appear to have lost their superscripts, and the DOTA2 base is 2 rather than 3 (as confirmed by the per-column cardinalities emitted at replication time). Table 6's bottom caption says $r_\mathrm{low}$ but the panel is the $r_\mathrm{high}$ setting. Rel. MSE is used in three tables but never defined (it is not $1 - R^2$). The proof of Remark 4.3 concludes with $U \subseteq A$ where the correct condition (matching the body) is $U \cap A \ne \emptyset$. The proof of Corollary 3.5 opens with a non-sequitur uniqueness step. Table 8's "Orthogonality" metric is never defined in prose.
• No error bars. Every quantitative row is a single-run point value even though the torch-based target models are not seeded in the released code, so the target function itself drifts between runs. The claim "instantaneously" for per-sample explanation also overstates cost at $r_\mathrm{low} = 10\,000$.
• Algorithm 1 pseudocode is out of sync with the described method. The while-loop uses rank(B) < r with no reference to $r_\mathrm{low}$, but every experiment runs the low-rank variant.
• LASSO/LARS future-work paragraph clashes with the identifiability discussion. In the sparse regime the overcomplete dictionary is by construction collinear, so calling the LASSO-pruned system "full-rank" without further argument is inconsistent with the paper's own analysis.
• Universal identifiability claim overstates the paper's analysis. The claim that "any ANOVA or GAM-like decomposition" on sparse categorical inputs faces the same non-uniqueness is asserted without proof; regularized/effect-coded alternatives escape in a well-defined sense.

Reproducibility & code

• Analytical, CAR, NURSERY, TreeHFD (headline aggregates). These reproduce well. The analytical case (Table 1) recovers $R^2 = 1$, MSE $\sim 10^{-21}$ with the basis correctly restricted to $\{X_1, X_2, X_4\}$. CAR EVALUATION ISE (24 cells) matches to two decimals; NURSERY ISE (40 cells) matches exactly. Appendix C's Table 8 reproduces on 4 of 5 dimensions with the headline closed-form-vs-TreeHFD orthogonality gap of ~12 orders of magnitude intact; per-cell drift is driven by unseeded $W_1, W_2$ in the synthetic $g(x)$.
• Table 6 top panel (main effects) reproduces; bottom panel does not. Main-effects $r_\mathrm{main}$ counts hit exactly (76 / 83 / 112) and $R^2$/MSE/Rel. MSE match on all three datasets. But the entire $r_\mathrm{high}$ panel was killed by the wall-clock cutoff (high_rank: null in all three emitted JSONs), so the paper's central "high rank helps" numbers ($R^2 = 0.79 / 0.70 / 0.41$) are not independently verified here.
• MUSHROOMS "$10^{-15}$ MSE" does not reproduce. The claim "$R^2 \approx 1$ with a negligible MSE of $10^{-15}$" reproduces on $R^2 = 0.999992$ but the MSE reproduces at $1.97 \times 10^{-6}$ — nine orders of magnitude above the paper. The qualitative "main effects fully reconstruct the signal" claim survives; the "machine-precision" claim as stated does not.
• Timings do not reproduce. CAR EVALUATION takes 21.4 s vs the paper's 0.5 s (~40x slower even accounting for shared-node contention, which the CPU/wall ratio of 7x cannot fully explain); NURSERY takes 194.7 s vs 54 s. The paper's ~100x scaling gap between the two datasets is not reproduced (replicated ratio ~9x). The 0.5 s / 54 s figures may refer to a narrower operation than the decomposition-plus-Shapley cell, but this is not stated.
• Table 7 (MNIST) is only half reproducible. The shipped notebook hardcodes n_set=2000 and does not sweep the four budgets in Table 7. When a driver loop was added, rows 1-2 ($r_\mathrm{low} = 674, 2500$) matched tightly; rows 3-4 ($5000, 10\,000$) were not reached before the cutoff.
• Figure 3 (POKER HAND) is not reproducible from the released code. There is no attention-rollout implementation anywhere in the repository, and the Eq. (23) importance-aggregation step (residual splitting + $\ell_1$ normalization) exists nowhere in the shipped code — the replicator transcribed it from the manuscript. Together with the missing $r_\mathrm{high}$ POKER state, neither series in Figure 3 can be regenerated end-to-end from the bundle.
• Figure 2 (MUSHROOMS) TreeSHAP baseline is not shipped. The notebook plots only the ANOVA main-effects L1 norms; there is no shap.TreeExplainer(model).shap_values(...) call and no side-by-side plotting. Once the replicator wrote the Observational (tree_path_dependent) TreeSHAP path from scratch, the numbers matched the paper to two decimals for all ten features — so the claim holds, but the released code does not verify it.
• Reproducibility hygiene. Torch target models (NURSERY MLP, POKER Transformer, DOTA2 MLP, MNIST MLP) are not torch.manual_seed-ed, so Tables 4/6/7 drift run-to-run. The scikit-learn/XGBoost experiments (CAR, CONNECT-4, MUSHROOMS) are seeded and deterministic. Requirements.txt is unpinned; on Python 3.12 the resolver picks numba 0.53.1 which refuses to build (a numba>=0.61 floor fixes it).
• MUSHROOMS $r_\mathrm{low} = 86$ comes from an undocumented magic number (tol_percent = 1.07, chosen so the loop stops at the last singleton set); this dataset-specific tuning is not surfaced in the paper.

Recommended Changes

Essential.

• Add a basis-order sensitivity experiment. For at least MUSHROOMS and POKER HAND, re-run Algorithm 1 under several random enumeration orders and report the spread of the top-$k$ feature rankings and of the local Shapley values. This addresses the non-uniqueness weakness that is currently the biggest gap between the theory and the interpretive claims.
• Add a Shapley-family baseline on POKER HAND (and preferably other sparse datasets). Replace or augment attention rollout with Interventional/Marginal SHAP, DeepSHAP, or KernelSHAP; this maps directly to the "attention rollout is a weak baseline" and "missing sparse-regime Shapley baseline" weaknesses.
• Fix Table 5 rendering. Rebuild the table so $|E|$ shows the intended $52^5$, $3^{42}$, and $2^{113}$ (verify the DOTA2 base from the per-feature cardinalities emitted by the code — it is 2, not 3), and add a one-line note stating $|E| = \prod_i N_i$.
• Ship the missing analysis code. Add shap.TreeExplainer(...) and the ANOVA-vs-TreeSHAP grouped bar plot to Exp-3 (mushrooms).ipynb; add attention-rollout extraction and the Eq. (23) aggregation step to Exp-4 (poker).ipynb; add a driver loop over $r_\mathrm{low} \in \{674, 2500, 5000, 10\,000\}$ to Exp-5 (mnist).ipynb that emits Table 7. This closes the largest gaps in Reproducibility & code.
• Seed torch training and report per-seed variance. Add torch.manual_seed(...) in the four torch notebooks and report mean $\pm$ std over 3-5 seeds for each row of Tables 4, 6, and 7, so the "no error bars" concern is addressed.
• Reconcile the MUSHROOMS "$10^{-15}$" claim with what the shipped code produces. Either state which error norm the $10^{-15}$ refers to (it is not A.get_L2_Error() as-run) or investigate the observed $10^{-6}$ residual on modern numpy/scipy and adjust the claim.
• Report on the high-rank Table 6 rows with documented compute. State explicitly the wall-clock budget needed to reproduce the $r_\mathrm{high}$ panel (a reproducer will not see 10 min / 37 min / 39 min on a shared node), and — better — ship a driver script that produces the panel end-to-end from a fresh clone.
• Fix the proof of Remark 4.3. Correct the concluding index set from $U \subseteq A$ to $U \cap A \ne \emptyset$ so the proof matches the body.
• Soften the abstract and the "must coincide" wording. Rephrase "without any assumption" to distinguish existence (assumption-free) from uniqueness (needs full support), and qualify the KernelSHAP equality as asymptotic/interventional so it is consistent with the nonzero ISE reported in Tables 3-4 and reproduced by replication.
• Update Algorithm 1 to make $r_\mathrm{low}$ an explicit input, or ship a second box for the low-rank variant with rank(B) < r_low as the stopping criterion.

Suggested.

• Define Rel. MSE in the caption of Table 6 (and reuse in Table 7); currently three related error columns are reported side by side with no stated relationship, and Rel. MSE is not $1 - R^2$.
• Fix the Table 6 bottom-panel caption to read $r_\mathrm{high}$.
• Define the Table 8 orthogonality metric in prose (max $|\langle f_A, f_B \rangle|$ over $B \subsetneq A$, with normalization used) and state whether the closed-form entries are exact-by-construction or Cholesky residuals — this matters given the 12-order gap is the headline claim there.
• Justify (or run a sensitivity check on) the residual-splitting convention on POKER HAND, showing the ranks-vs-suits conclusion is robust to alternatives (residual over active features only, or residual proportional to marginal variance). This is especially important for DOTA2, where the residual is 60% of the variance.
• Rephrase local claims that outrun what is shown: "decay of information per set cardinality", "already recovered by the main effects" (dimension-counting artefact), "instantaneously" (per-sample cost is $O(|S_{r_\mathrm{low}}|)$), the "any ANOVA or GAM-like decomposition" identifiability claim (only shown for this dictionary), and the LASSO/LARS "pruned, full-rank system" wording (collinearity is exactly the identifiability issue).
• Expose MUSHROOMS $r_\mathrm{low}$ as a first-class parameter in the interface used in Exp-3, rather than via the tol_percent = 1.07 magic number.
• Pin Requirements.txt (or add a numba>=0.61 floor) so a fresh install on Python 3.12 builds without manual intervention.
• Report attribution error under low-rank truncation on a controllable synthetic (e.g., a scaled-up variant of the Section 5.1 example) so readers can gauge how the truncated ANOVA-based Shapley diverges from a known ground truth as $r_\mathrm{low}$ varies. This is the missing quantitative link that would justify recommending the low-rank variant in practice.