PolySHAP: Polynomial Shapley Explanations

• PolySHAP is an explainable AI method that uses polynomial regression to derive consistent and efficient Shapley value estimates while incorporating higher-order feature interactions.
• The approach models the explanation game with a multilinear polynomial that extends beyond additive effects, improving empirical accuracy and addressing computational barriers.
• Empirical benchmarks confirm that PolySHAP enhances error metrics and precision, with theoretical guarantees ensuring convergence to true Shapley values via a Möbius transformation.

PolySHAP is a method in explainable artificial intelligence (XAI) designed to estimate Shapley values more accurately and efficiently by fitting low-degree polynomials to the model-explanation game, thereby capturing feature interactions absent from traditional KernelSHAP. Shapley values, grounded in cooperative game theory, quantify individual feature contributions by averaging marginal effects over all feature subset combinations, but their direct computation requires $2^d$ game evaluations — an exponential barrier for models with moderately large $d$. KernelSHAP ameliorates this cost by approximating the explanation game as a linear function via weighted least squares over a sampled subset of feature combinations, accounting only for additive effects. PolySHAP generalizes KernelSHAP by fitting higher-order (degree-$k$) polynomial models to efficiently incorporate non-linear feature interactions, providing provably consistent Shapley value estimates with improved empirical accuracy on benchmark datasets (Fumagalli et al., 26 Jan 2026).

1. Limitations of KernelSHAP and Additive Approximation

KernelSHAP estimates Shapley values for a model-agnostic game $\nu:2^D\to \mathbb{R}$ by solving a weighted linear least squares problem. The exact formulation is:

$$\phi^{SV}[\nu] = \underset{\phi\in\mathbb{R}^d,\; \sum_i\phi_i=\nu(D)}{\arg\min} \sum_{S\subseteq D} \mu(S) \left(\nu(S) - \sum_{i\in D}\phi_i 1[i\in S]\right)^2,$$

where the weight $\mu(S) = 1/\binom{d-2}{|S|-1}$ for $0<|S|$m\ll 2^d$ subsets $S_1,\ldots,S_m$ from a distribution $p(S)$ (commonly proportional to $\mu$) and fits the linear model:

$$\hat{\phi} = \underset{\phi,\; \sum\phi_i = \nu(D)}{\arg\min} \sum_{\ell = 1}^m \frac{\mu(S_\ell)}{p(S_\ell)}\left(\nu(S_\ell) - \sum_i \phi_i1[i\in S_\ell]\right)^2.$$

This procedure is computationally tractable, scaling as $O(m d^2)$, but is constrained to additive effects—it cannot capture interactions beyond linearity.

2. PolySHAP Polynomial Regression Formulation

PolySHAP constructs a multilinear polynomial of degree $k$ over the feature set, generalizing the additive model. The "interaction frontier" is given by:

$I_{\leq k} = \{T\subseteq D:\; 1\leq|T|\leq k\},$

with one coefficient $\phi_T$ per subset $T\in I_{\leq k}$. The polynomial approximation is:

$$\nu(S) \approx f_\phi(S) := \sum_{T\in I_{\leq k}} \phi_T \mathbf{1}[T\subseteq S],$$

where $\mathbf{1}[T\subseteq S]=\prod_{j\in T} 1[j\in S]$. The fitting objective is:

$$\hat{\phi} = \underset{\phi\in\mathbb{R}^{d'}:\;\sum_T\phi_T = \nu(D)}{\arg\min} \sum_{\ell=1}^m \frac{\mu(S_\ell)}{p(S_\ell)}\left(\nu(S_\ell) - \sum_{T\in I_{\leq k}} \phi_T \mathbf{1}[T\subseteq S_\ell]\right)^2,$$

where $d' = |I_{\leq k}| = \sum_{i=1}^k \binom{d}{i}$. This is solved via constrained least squares regression, using a design matrix $$\tilde{X}_{\ell,T} = \sqrt{\mu(S_\ell)/p(S_\ell)} \mathbf{1}[T\subseteq S_\ell]$$ and target vector $$\tilde{y}_\ell = \sqrt{\mu(S_\ell)/p(S_\ell)}\nu(S_\ell)$$. The solution is subject to the constraint $\mathbf{1}^T\phi = \nu(D)$, handled by standard projection techniques.

3. Shapley Value Extraction via Möbius–Shapley Transformation

After determining $\hat{\phi}$, the estimated Shapley values are recovered by mapping monomial coefficients to individual feature attributions via:

$$\phi_i^{SV}[\nu] = \hat{\phi}_{\{i\}} + \sum_{T\in I_{\leq k}:\,i\in T} \frac{\hat{\phi}_T}{|T|}.$$

For $k=1$, this reduces to the KernelSHAP additive solution $\phi_i^{SV} = \hat{\phi}_{\{i\}}$. For $k>1$, higher-order coefficients correct the Shapley estimate for non-additive interactions. This transformation, rigorously derived in [(Fumagalli et al., 26 Jan 2026), Theorem 3.1], constitutes a Möbius–Shapley conversion and ensures coherent attribution, including interaction effects.

4. Consistency Guarantees

PolySHAP is consistent in the sense that, as the number of sampled subsets $m\to 2^d$, the estimated polynomial coefficients $\hat{\phi}^{I}$ converge to the population minimizer, and the resulting Shapley values $\phi_i^{SV}[\nu]$ converge to the true Shapley values for the explanation game. The proof strategy employs (i) recasting both degree-1 and degree-$k$ weighted least squares problems into constrained matrix regression, (ii) use of a projection lemma (Lemma A.1) establishing equivalence between degree-$k$ and projected degree-1 solutions in full enumeration, and (iii) demonstration that reconstructed Shapley values via Möbius–Shapley transformation recover the population solution exactly, under full-rank design and unbiased sampling assumptions (Fumagalli et al., 26 Jan 2026).

5. Paired Sampling and Algebraic Equivalence to Quadratic PolySHAP

Paired (antithetic) sampling is a KernelSHAP heuristic involving sampling subsets in complementary pairs ($S$ and $D\setminus S$), empirically observed to improve accuracy by reducing estimator variance. PolySHAP establishes a formal equivalence: KernelSHAP with paired sampling solves the same normal equations as 2-PolySHAP projected to degree-1. Specifically:

$$\text{KernelSHAP}_{\text{paired}} = \text{PolySHAP}_{k=2} \rightarrow \text{Shapley}.$$

The proof rests on symmetries in the cross-moment matrix under paired data, which ensure all quadratic columns collapse to effective additive corrections, yielding identical Shapley estimates [(Fumagalli et al., 26 Jan 2026), Theorem 4.1]. This result provides a theoretical foundation for the practical success of paired sampling and shows it implicitly incorporates all second-order interactions at no additional computational cost.

6. Empirical Performance and Practical Recommendations

Empirical benchmarks conducted in (Fumagalli et al., 26 Jan 2026) demonstrate that PolySHAP improves mean squared error (MSE), top-5 precision, and Spearman correlation over additive KernelSHAP as the polynomial degree $k$ increases ($k=2,3,4$), provided the number of samples $m\gg d'$. Paired sampling enables KernelSHAP ($k=1$) to perform equivalently to quadratic PolySHAP ($k=2$), but higher-order ($k\geq 3$) models yield further improvements when computational budgets allow. Recommended regime for moderate $d$ is $k \approx 2$ or $3$; for high-dimensional settings, employing a "partial frontier" (random subset of degree-$k$ monomials) is suggested. Leverage-score sampling ($p(S)\propto$ row leverage) further reduces variance, guaranteeing $m = O(d'\log d')$ for $\epsilon$-accurate fits. Total computational cost is $O(m d'^2 + d'^3)$ for regression and $O(d d')$ for Shapley extraction, with model evaluation typically the dominant factor except when $d'$ is extremely large.

7. Summary and Scope

PolySHAP generalizes KernelSHAP by fitting interaction-informed low-degree polynomial models to the model-agnostic explanation game, enabling consistent recovery of Shapley values including non-linear interactions. The method is grounded in weighted least squares regression constrained to sum-to-game-total, with Shapley value extraction via Möbius transformation. The use of paired sampling is formally justified as algebraically equivalent to including all second-order interactions. PolySHAP provides improved statistical accuracy and theoretical consistency while remaining computationally feasible for modest degree and dimension. Empirical guidelines favor degree $k=2$ or $3$ and leverage-score or paired sampling strategies for practical deployment (Fumagalli et al., 26 Jan 2026).