Permutation SHAP: Sampling-Based Shapley Estimation

• Permutation SHAP is a model-agnostic method that estimates Shapley values by averaging marginal contributions over random feature orderings.
• It utilizes a paired-sampling trick by coupling each permutation with its reverse to reduce variance, achieve asymptotic normality, and exactly recover block-additive attributions.
• Extended approaches include sampling from the product of marginals and applying kernel herding or sequential Bayesian quadrature to enhance convergence in high-dimensional settings.

Searching arXiv for recent and foundational papers relevant to Permutation SHAP. Permutation SHAP denotes a family of Shapley-value procedures in which permutations play the central computational role. In its standard usage, PermutationSHAP is a simple, model-agnostic way to compute Shapley attributions by averaging marginal contributions over random feature orderings; in its paired-sampling form, each sampled ordering is coupled with its reverse, yielding asymptotic normality, exactness for purely second-order interactions, and an additive-recovery guarantee for block-additive value functions (Mayer et al., 18 Aug 2025). In a separate recent usage, “permutation SHAP” refers to computing SHAP or KernelSHAP on data obtained by independently permuting each column, thereby sampling from the product of marginals $\mu^*$ and enabling theoretically justified feature discarding from aggregate scores (Bhattacharjee et al., 29 Mar 2025).

1. Definition, notation, and scope

Let $f:\mathbb{R}^p\to\mathbb{R}$ be a predictor and let $x\in\mathbb{R}^p$ be the instance to explain. For any subset $S\subseteq\{1,\dots,p\}$, write $x_S$ for the input where only features in $S$ are retained and the rest are masked, intervened, or set to some baseline. If $\pi$ is a permutation of $\{1,\dots,p\}$ and $S_i(\pi)$ is the set of features preceding $i$ in $f:\mathbb{R}^p\to\mathbb{R}$0, then the classical Shapley value is

$f:\mathbb{R}^p\to\mathbb{R}$1

Equivalently, with $f:\mathbb{R}^p\to\mathbb{R}$2,

$f:\mathbb{R}^p\to\mathbb{R}$3

Because $f:\mathbb{R}^p\to\mathbb{R}$4 becomes infeasible once $f:\mathbb{R}^p\to\mathbb{R}$5, PermutationSHAP replaces the full average over $f:\mathbb{R}^p\to\mathbb{R}$6 by sampling random permutations (Mayer et al., 18 Aug 2025).

The permutation formulation is also written as an expectation over $f:\mathbb{R}^p\to\mathbb{R}$7 in the sampling literature. With $f:\mathbb{R}^p\to\mathbb{R}$8, one has

$f:\mathbb{R}^p\to\mathbb{R}$9

which makes the approximation problem a numerical integration problem on permutation space. In that setting, a simple unbiased Monte Carlo estimator converges at $x\in\mathbb{R}^p$0 and requires $x\in\mathbb{R}^p$1 calls to $x\in\mathbb{R}^p$2 if one exploits the sequential structure of $x\in\mathbb{R}^p$3 (Mitchell et al., 2021).

The cited literature uses the phrase “permutation SHAP” in two related but distinct senses:

Usage	What is permuted	Primary purpose
PermutationSHAP	Feature orderings $x\in\mathbb{R}^p$4	Local Shapley estimation by averaging marginal gains
Extended-support permutation SHAP	Columns of the data matrix	Global aggregation over $x\in\mathbb{R}^p$5 for safe feature discarding

The first sense concerns local attribution for a fixed instance $x\in\mathbb{R}^p$6. The second concerns aggregate scores over an extended distribution and is therefore conceptually closer to global feature selection than to pointwise explanation.

2. Paired-sampling PermutationSHAP estimator

A paired-sampling trick reduces variance by coupling each sampled permutation $x\in\mathbb{R}^p$7 with its reverse $x\in\mathbb{R}^p$8. If $x\in\mathbb{R}^p$9 are i.i.d., define

$S\subseteq\{1,\dots,p\}$0

and

$S\subseteq\{1,\dots,p\}$1

The paired-sampling PermutationSHAP estimate is

$S\subseteq\{1,\dots,p\}$2

Its implementation is direct: initialize $S\subseteq\{1,\dots,p\}$3 for all features, repeatedly draw a uniform permutation, compute the marginal gains for that ordering and for its reverse, and accumulate their average. Because each iteration uses exactly $S\subseteq\{1,\dots,p\}$4 model-evaluations, the total cost is approximately $S\subseteq\{1,\dots,p\}$5 (Mayer et al., 18 Aug 2025).

This estimator is model-agnostic in the sense that it only requires the ability to evaluate the value function $S\subseteq\{1,\dots,p\}$6 for selected subsets. Its computational profile is therefore dominated by the number of model evaluations induced by the sampled orderings rather than by model-specific structure. In the summary provided for the method, the cost per sample is $S\subseteq\{1,\dots,p\}$7 model-evaluations (Mayer et al., 18 Aug 2025).

The ordering-based construction is important. Unlike regression-based approximations, it accumulates incremental gains along a path through feature space, and the reverse permutation introduces a symmetric counterpart to that path. This structure is what underlies the variance reduction and the exact recovery results stated below.

3. Asymptotic normality, exactness, and additive recovery

The paired estimator admits a central limit theorem. Define the random-vector contribution $S\subseteq\{1,\dots,p\}$8 by

$S\subseteq\{1,\dots,p\}$9

and let

$x_S$0

Then, as $x_S$1,

$x_S$2

In particular, each component $x_S$3 converges at rate $x_S$4, and $x_S$5 can be estimated empirically from the sampled permutations (Mayer et al., 18 Aug 2025).

A stronger result holds for purely second-order interactions. If the value function is a bilinear form,

$x_S$6

so that only feature-pairs interact, then even $x_S$7 suffices: for any single permutation $x_S$8,

$x_S$9

exactly. The proof sketch given for the result shows that, in the purely quadratic case, the paired sum becomes independent of $S$0 (Mayer et al., 18 Aug 2025).

The distinctive structural guarantee is the additive-recovery property. Suppose

$S$1

where $S$2 is a partition into disjoint additive blocks. Then for any single permutation $S$3, the paired-sampling estimator with one draw satisfies

$S$4

The proof sketch proceeds by telescoping the incremental gains along $S$5 within each block while features from other blocks act as dummy variables (Mayer et al., 18 Aug 2025).

These theorems sharply distinguish errors internal to a block from errors across blocks. The block totals are exact after one paired sample, whereas finer-grain errors within a block shrink at rate $S$6. This suggests that paired PermutationSHAP is particularly well matched to models with low-order interactions or a block-additive structure.

4. Comparison with KernelSHAP

KernelSHAP is presented in the cited material as a global weighted-least-squares procedure. In the paired-sampling comparison, it solves

$S$7

and its sampling version is a regression of $S$8 onto $S$9 (Mayer et al., 18 Aug 2025).

The two estimators share some asymptotic features. KernelSHAP is unbiased and converges at $\pi$0, while paired PermutationSHAP satisfies a CLT at rate $\pi$1 with explicitly computable covariance. The difference emphasized in the comparison is not consistency but structure: PermutationSHAP’s incremental-gain sampling respects the ordering-based telescoping that recovers block sums exactly, whereas the regression-based KernelSHAP “mixes” contributions across blocks (Mayer et al., 18 Aug 2025).

As a consequence, KernelSHAP does not, in general, return the correct aggregate $\pi$2 over an additive partition, even though it is unbiased in the large-sample limit. The additive-recovery guarantee is therefore exclusive to paired-sampling PermutationSHAP in the formulation summarized by the paper. In practice, the same summary states that paired-sampling halves variance, expressed as $\pi$3 versus $\pi$4, and that for purely second-order interactions even one pair yields exact $\pi$5 (Mayer et al., 18 Aug 2025).

The comparison is methodologically important because it identifies a property that is not a generic consequence of Shapley-value axioms alone. It arises from the specific pathwise estimator defined over permutations.

5. Sampling design on permutation space

Once Shapley estimation is written as an expectation over $\pi$6, the quality of the permutation samples becomes a numerical integration problem. Standard Monte Carlo sampling is unbiased but can exhibit slow convergence, and variance can be high for models with strong feature interactions. Variance-reduction schemes such as antithetic, stratified, and Owen’s multilinear improve constants but not the $\pi$7 rate (Mitchell et al., 2021).

One line of work formulates the problem in an RKHS over permutations. With a positive-definite kernel $\pi$8, permutations are embedded into an RKHS $\pi$9. The cited paper lists three kernels: the Kendall kernel, the universal Mallows kernel

$\{1,\dots,p\}$0

and the Spearman inner-product kernel. In all three cases, $\{1,\dots,p\}$1 can be computed in closed form (Mitchell et al., 2021).

Kernel herding greedily builds an unweighted set of permutations to reduce RKHS discrepancy, and for universal kernels such as Mallows it achieves $\{1,\dots,p\}$2 convergence in RKHS norm. Sequential Bayesian quadrature models the integrand on permutations as a zero-mean Gaussian process and selects permutations by minimizing posterior variance of the integral estimate. The worst-case quadrature error satisfies

$\{1,\dots,p\}$3

where $\{1,\dots,p\}$4 is the kernel discrepancy; herding drives $\{1,\dots,p\}$5 down for equal weights, and SBQ optimally chooses weights to minimize it (Mitchell et al., 2021).

A separate construction exploits the geometry of the permutohedron and the hypersphere $\{1,\dots,p\}$6. After centering and scaling, permutations lie on a translated copy of $\{1,\dots,p\}$7, and a point on the sphere can be mapped back to a permutation via rank ordering. The paper gives two fast schemes: orthogonal spherical codes, which produce $\{1,\dots,p\}$8 permutations with negative correlation and good coverage, and Sobol-on-the-sphere sampling, which uses a low-discrepancy sequence, inverse-CDF transforms for hyperspherical coordinates, and an $\{1,\dots,p\}$9 map back to permutation space. Both yield uniform $S_i(\pi)$0 in expectation and run in $S_i(\pi)$1 time (Mitchell et al., 2021).

Empirically, the same study reports that on small–medium $S_i(\pi)$2 with $S_i(\pi)$3, kernel herding and SBQ with Mallows $S_i(\pi)$4 achieve $S_i(\pi)$5 error decay and outperform MC’s $S_i(\pi)$6, while on larger $S_i(\pi)$7 or very expensive models the linear-time sphere methods often match or beat antithetic MC by stronger negative correlation or low discrepancy. The experiments show convergence to a smaller RMSE in the same number of model evaluations (Mitchell et al., 2021).

6. Extended-support permutation SHAP and feature discarding

A separate line of work studies aggregate SHAP values for feature removal and shows that small aggregate SHAP values on the original data support do not necessarily imply that the corresponding feature does not affect the function. The issue is that SHAP evaluates $S_i(\pi)$8 on points outside the data support, where dependence on a feature can be masked. The proposed remedy is to aggregate SHAP values over the extended distribution

$S_i(\pi)$9

whose support is $i$0 (Bhattacharjee et al., 29 Mar 2025).

Sampling from $i$1 can be implemented by independently permuting each column of the observed data matrix. If $i$2 is the original data and $i$3 is an independent random permutation of $i$4 for each column $i$5, define

$i$6

Each row of $i$7 is then an approximate sample from $i$8 because coordinates are independent across columns. The paper defines the extended-support aggregate

$i$9

Its main theorem states that

$f:\mathbb{R}^p\to\mathbb{R}$00

A robust version states that if $f:\mathbb{R}^p\to\mathbb{R}$01 and $f:\mathbb{R}^p\to\mathbb{R}$02, then there exists $f:\mathbb{R}^p\to\mathbb{R}$03 not depending on $f:\mathbb{R}^p\to\mathbb{R}$04 such that

$f:\mathbb{R}^p\to\mathbb{R}$05

Thus, a small average absolute SHAP value over the extended support guarantees that one can replace $f:\mathbb{R}^p\to\mathbb{R}$06 by a function independent of feature $f:\mathbb{R}^p\to\mathbb{R}$07 with small squared error under $f:\mathbb{R}^p\to\mathbb{R}$08 (Bhattacharjee et al., 29 Mar 2025).

The same framework extends to KernelSHAP. KernelSHAP fits a weighted linear model on subsets sampled from the Shapley-kernel distribution

$f:\mathbb{R}^p\to\mathbb{R}$09

To obtain “Permutation KernelSHAP” on $f:\mathbb{R}^p\to\mathbb{R}$10, one runs KernelSHAP on the scrambled data $f:\mathbb{R}^p\to\mathbb{R}$11 so that all value-function evaluations use samples from $f:\mathbb{R}^p\to\mathbb{R}$12 instead of $f:\mathbb{R}^p\to\mathbb{R}$13. If the average of $f:\mathbb{R}^p\to\mathbb{R}$14 over $f:\mathbb{R}^p\to\mathbb{R}$15 is at most $f:\mathbb{R}^p\to\mathbb{R}$16, then there exists a function $f:\mathbb{R}^p\to\mathbb{R}$17 independent of $f:\mathbb{R}^p\to\mathbb{R}$18 with

$f:\mathbb{R}^p\to\mathbb{R}$19

where $f:\mathbb{R}^p\to\mathbb{R}$20 is the KernelSHAP estimation error (Bhattacharjee et al., 29 Mar 2025).

The operator-theoretic explanation is expressed through the Shapley Lie algebra. The linear operators

$f:\mathbb{R}^p\to\mathbb{R}$21

generate a solvable Lie algebra, and on finite-dimensional invariant subspaces they can be simultaneously triangularized. In the summary provided, this yields invertibility of the operator $f:\mathbb{R}^p\to\mathbb{R}$22 and identifies the kernel of the SHAP operator $f:\mathbb{R}^p\to\mathbb{R}$23 with the subspace of $f:\mathbb{R}^p\to\mathbb{R}$24-independent functions. Within that framework, column-wise permutation is not merely a heuristic data-augmentation step; it is the mechanism by which one samples from $f:\mathbb{R}^p\to\mathbb{R}$25 and obtains the safe-discarding guarantees (Bhattacharjee et al., 29 Mar 2025).