Generalized Functional ANOVA in Closed-Form: A Unified View of Additive Explanations

Abstract: The functional ANOVA, or Hoeffding decomposition, provides a principled framework for interpretability by decomposing a model prediction into main effects and higher-order interactions. For independent inputs, this classical decomposition is explicit. It is closely connected to SHAP values, generalized additive models, and orthogonal polynomial expansions, and therefore constitutes a fundamental tool for additive explainability. In the more general and realistic dependent setting, however, obtaining a tractable representation and estimating the decomposition from data remain challenging. In this work, we address this problem for continuous inputs. By combining Hilbert space methods with the generalized functional ANOVA, we build an explicit decomposition Riesz Basis allowing to easily compute the decomposition. Our formulation recovers the classical independent case and its associated orthogonal decomposition. Building on this representation, we propose a simple but mighty algorithm to estimate the decomposition from a data sample in a model-agnostic setting and we compare it empirically with several state-of-the-art explanation methods, demonstrating the power of the approach.

• The paper derives a closed‐form functional ANOVA decomposition that generalizes classical and Shapley-based methods under dependent features.
• It introduces a model-agnostic estimator that scales polynomially, enforces hierarchical orthogonality, and accurately attributes feature effects.
• Empirical results demonstrate near-complete variance recovery (R²≈0.85–1.00) on real-world datasets, outperforming established explainable AI techniques.

Generalized Functional ANOVA in Closed-Form: A Unified View of Additive Explanations

Introduction and Motivation

The opaque nature of modern nonlinear predictors, especially deep neural networks and tree ensembles in tabular data settings, creates significant challenges for interpretability and trustworthy machine learning. Most approaches for explainable AI (XAI) either constrain the model class to be interpretable (e.g., Generalized Additive Models, EBMs, NAMs), or provide post-hoc explanations by decomposing model outputs into feature attributions (SHAP, LIME, DeepSHAP). However, these approaches suffer from either computational intractability, use of heuristic attributions, or the lack of a principled, unified framework for additive explanations under dependence among input features.

The paper "Generalized Functional ANOVA in Closed-Form: A Unified View of Additive Explanations" (2605.18422) systematically addresses these limitations. It derives a closed-form, explicit, functional ANOVA decomposition for predictor models with continuous, potentially dependent, bounded inputs, employing Hilbert space machinery. The work generalizes classical ANOVA/Hoeffding decompositions by constructing a Riesz basis satisfying hierarchical orthogonality—a nontrivial extension under arbitrary dependence. The authors then propose a practical, ultra-fast estimator of the decomposition, operating purely in a model-agnostic regime, and demonstrate its empirical advantage and unifying capacity relative to state-of-the-art explanation and interpretable model classes.

Generalized Functional ANOVA: Theory and Representation

Functional ANOVA (Hoeffding decomposition) expresses any square-integrable function $\nu(\mathbf{X})$ as

$$\nu(\mathbf{X}) = \nu_\emptyset + \sum_{i=1}^p \nu_i(\mathbf{X}_i) + \sum_{i < j} \nu_{i,j}(\mathbf{X}_i, \mathbf{X}_j) + ...$$

where $\nu_S$ is a function of only the subset $\mathbf{X}_S$ and lower-order terms are constrained to capture only irreducible effects (hierarchical orthogonality).

When features are independent, component expressions are explicit via the Möbius transform of conditional expectations, allowing orthogonal decomposition, direct variance allocation, and connection to Shapley values (via Harsanyi dividends). For dependent features, prior approaches either require explicit knowledge of joint/marginal densities or yield computationally intractable estimands.

This paper’s main theoretical breakthrough is the characterization and explicit construction of a Riesz basis for $L^2$ under arbitrary continuous distributions with bounded density, utilizing scaled Legendre polynomials and an inverse-likelihood correction mechanism to ensure hierarchical orthogonality. The resulting components $\nu_S$ resolve the constrained optimization posed by Hooker [2007], and yield closed-form expressions generalizing both classical ANOVA and Shapley-based decompositions.

Figure 1: Main effects for California Housing; the method’s component estimates (black) align against TreeHFD and TreeSHAP attributions on an XGB model, capturing additive structure robustly under dependence.

Algorithmic Framework and Estimation

Operationalizing the functional ANOVA decomposition, the estimator is constructed by minimizing empirical least squares over the Riesz basis functions, with coefficients indexed by subset and polynomial order. To ensure tractability, an interaction order ($K$) and polynomial degree ($d$) truncation is adopted, leveraging the "sparsity-of-effects" and "reluctance" principles established in the literature on high-dimensional interaction modeling.

This estimator is fully model-agnostic: it requires only black-box access to the prediction function and samples from the input distribution (or observed data points). Marginal densities needed in the denominator of basis functions are also estimated from the data, typically by low-degree tensorized Legendre expansions.

The computational pipeline comprises:

1. Feature normalization to $[-1,1]^p$ (componentwise $\tanh$ scaling)
2. Marginal density estimation (degree-$$\nu(\mathbf{X}) = \nu_\emptyset + \sum_{i=1}^p \nu_i(\mathbf{X}_i) + \sum_{i < j} \nu_{i,j}(\mathbf{X}_i, \mathbf{X}_j) + ...$$0 Legendre expansion)
3. Design matrix construction for all relevant basis functions
4. LARS-based model selection + SVD-based sparse linear solving for coefficients

This practical procedure scales polynomially in sample size and feature number for low-order $$\nu(\mathbf{X}) = \nu_\emptyset + \sum_{i=1}^p \nu_i(\mathbf{X}_i) + \sum_{i < j} \nu_{i,j}(\mathbf{X}_i, \mathbf{X}_j) + ...$$1 (usually $$\nu(\mathbf{X}) = \nu_\emptyset + \sum_{i=1}^p \nu_i(\mathbf{X}_i) + \sum_{i < j} \nu_{i,j}(\mathbf{X}_i, \mathbf{X}_j) + ...$$2 suffices in tabular data), and enforces empirical hierarchical orthogonality among components.

Figure 2: Decomposition for a trained MLP on Bike Sharing: local attribution network plot (left); comparison of estimated main effects for hour and atemp by this method (black) versus KernelSHAP and DeepSHAP (center, right).

Empirical Evaluation and Comparative Analysis

Extensive experiments span synthetic (analytical) as well as diverse real-world tabular datasets (California Housing, Census Income, Bike Sharing, etc.), evaluating various target models (XGB, MLP, EBM, NAM). The key points are:

• The proposed estimator matches or surpasses SHAP-based and HFD-based methods in reconstructing theoretical main effects, especially in the presence of dependent or correlated features.
• Consistent $$\nu(\mathbf{X}) = \nu_\emptyset + \sum_{i=1}^p \nu_i(\mathbf{X}_i) + \sum_{i < j} \nu_{i,j}(\mathbf{X}_i, \mathbf{X}_j) + ...$$3 values in the $$\nu(\mathbf{X}) = \nu_\emptyset + \sum_{i=1}^p \nu_i(\mathbf{X}_i) + \sum_{i < j} \nu_{i,j}(\mathbf{X}_i, \mathbf{X}_j) + ...$$4 range demonstrate near-complete variance recovery by main and pairwise terms, confirming the adequacy of low-order truncated decomposition.
• The computational cost remains minimal even for large datasets (e.g., $$\nu(\mathbf{X}) = \nu_\emptyset + \sum_{i=1}^p \nu_i(\mathbf{X}_i) + \sum_{i < j} \nu_{i,j}(\mathbf{X}_i, \mathbf{X}_j) + ...$$5 to $$\nu(\mathbf{X}) = \nu_\emptyset + \sum_{i=1}^p \nu_i(\mathbf{X}_i) + \sum_{i < j} \nu_{i,j}(\mathbf{X}_i, \mathbf{X}_j) + ...$$6 data points), validating the estimator’s practical merit.
• Hierarchical orthogonality is empirically maintained, as quantified by maximum component-wise cosines well below $$\nu(\mathbf{X}) = \nu_\emptyset + \sum_{i=1}^p \nu_i(\mathbf{X}_i) + \sum_{i < j} \nu_{i,j}(\mathbf{X}_i, \mathbf{X}_j) + ...$$7.
• Existing explanation methods (KernelSHAP, TreeSHAP, DeepSHAP, TreeHFD), and interpretable models (NAM, EBM), are shown to implicitly or explicitly approximate subcomponents of the generalized functional ANOVA; only this method guarantees uniqueness and hierarchical orthogonality universally.

Figure 3: Main effect estimation for Age on Census Income—method’s output (black) versus KernelSHAP/DeepSHAP (MLP, left) and TreeHFD/TreeSHAP (XGB, right); correct attribution and denoising of irrelevant variable.

Figure 4: Agreement between natively-learned main effects in EBM and NAM and those reconstructed post-hoc using the proposed ANOVA estimator on Diabetes; the decomposition aligns precisely, supporting the "unifying" claim.

Theoretical and Practical Implications

The ability to construct closed-form functional ANOVA decompositions for continuous, dependent features establishes a rigorous foundation for additive explainability. This result mathematically subsumes and clarifies the relationships between ANOVA, Shapley value attributions, and model-based explanations (GAM, EBM/NAM), while removing the need for internal model access or explicit density knowledge.

• Theoretical implications: The unique basis with hierarchical orthogonality allows well-defined variance allocation, interaction strength measurement, and unbiased interpretation of feature contributions in black-box regression/classification.
• Practical implications: The estimator operates on large datasets in polynomial time and is agnostic to the model class, enabling plug-and-play explainability for any black-box predictor in tabular analysis.

Limitations and Future Directions

Several limitations are acknowledged:

• The method is currently limited to compact (bounded) continuous domains due to the use of Legendre polynomials; extensions to unbounded, hybrid, or categorical features require new orthogonalization schemes.
• The curse of dimensionality surfaces theoretically via exponential subset enumeration; reliance on effect-sparsity and low-order truncation addresses this in practice, but more sophisticated penalization/selection strategies are desirable.
• The density estimation module could be improved with advanced generative or kernel approaches, and the linear solver could benefit from hardware acceleration.

Future research should investigate general penalized estimators (sparse/max-norm), basis construction for mixed-type or unbounded supports, and functional alignment diagnostics for model auditing.

Conclusion

This work delivers a mathematically principled, computationally efficient, and unifying framework for additive explanation of black-box predictors via a closed-form, generalized functional ANOVA decomposition under dependent, continuous features. The connection to existing explanation paradigms is clarified, and the practical estimator provides highly accurate, hierarchically orthogonal attributions with minimal computational burden. This marks a theoretical synthesis and practical consolidation of feature attribution in complex machine learning models.