Provably Accurate Shapley Value Estimation via Leverage Score Sampling

In plain words

Shapley values explain a model's prediction by sharing out credit among its input features: each feature gets a number for how much it moved the output. Computing them exactly means asking the model about every combination of features, which grows too fast to be practical as the number of features rises.

Kernel SHAP is the popular shortcut. It asks the model about a sample of feature combinations and solves a regression problem to recover the credit numbers. Which combinations it samples matters a lot, and its sampling rule was a heuristic that had no accuracy guarantee at the time.

This paper swaps that rule for leverage-score sampling, a standard way to pick the most informative rows of a regression. The authors show that the scores for this problem have a simple form, so the sample can be drawn without listing every combination. The result, Leverage SHAP, is a lightweight change to Kernel SHAP that comes with accuracy guarantees. In the authors' experiments it has lower error than the optimized Kernel SHAP in the widely used SHAP library.

TL;DR

Leverage SHAP changes only how Kernel SHAP chooses its feature subsets. It samples them with leverage scores, which have a closed form, and then solves a projected weighted least-squares problem. The paper gives a guarantee on the regression objective, not directly on the Shapley values. Author-reported: about 50% lower error than Optimized Kernel SHAP, averaged over budgets from 5n to 160n evaluations. Not independently re-run.

Problem

Exact Shapley values need one model evaluation per feature subset, so O(2^n) for n features §1.1. Kernel SHAP reduces this by estimating the values from m sampled subsets, with rows weighted by w(|S|) = 1 / (C(n, |S|) · |S| · (n − |S|)) eq. 4. That sampling distribution is a heuristic. When the paper was submitted, no non-asymptotic accuracy guarantee was known for Kernel SHAP with fewer than 2^n − 2 samples §1.2. Noisy value estimates, for example expectations over a random baseline, make the problem harder and are the setting the authors target §1.2.

Method

The method rests on a regression view of Shapley values. Write z for a 0/1 indicator of a feature subset, v(z) for the model's value on that subset, and n for the number of features. The Shapley vector φ is the solution of a least-squares problem over the rows z with 0 < |z| < n, weighted by w(|z|), subject to the constraint that the values sum to v(all) − v(∅) §2, eq. 6–7.

1. Remove the constraint. Let P = I − 11ᵀ/n, A = ZP and b = y − Z1 (v(all) − v(∅)) / n. The constrained problem becomes an unconstrained least-squares problem in A and b, and the answer is shifted back by (v(all) − v(∅))/n · 1 §3, eq. 8, eq. 9.
2. Leverage scores have a closed form. Because AᵀA = P/n, the leverage score of the row for subset z is 1 / C(n, |z|) §3.1. All subsets of one size share a score, so sampling proportional to leverage gives every size the same total weight. Equivalently, draw a size s uniformly from 1 to n − 1, then a subset of that size uniformly §1.2, fig. 2.
3. Pair the samples. Each sampled subset is evaluated together with its complement, because the two have equal leverage and equal kernel weight §3.2. The idealised scheme keeps each pair independently with probability min(1, 2c·ℓ_z), where c is a scale chosen so the expected number of evaluations is m − 2 §3.2, eq. 13. Each kept row is weighted by w(|z|) divided by its inclusion probability alg. 1.
4. Sample without replacement, size by size. The sampler needs one binomial or fixed count per size and then a uniform draw of subsets of that size, so it never enumerates the 2^n rows app. J, alg. 2, alg. 3. For even n the middle size is handled separately so that a pair is never drawn twice app. J.
5. Solve. Take the weighted least-squares solution of the sampled, projected system. Use the minimum-norm solution, because AᵀA is singular alg. 1 §3.

Algorithm 1 of the paper is the theoretical version. The released estimator differs in step 3 and step 4 (see Paper vs code).

# Leverage SHAP with the released per-size sampler  [alg. 1] [app. J]
# Input:  n players, v (set function), m (budget; paper requires m >= n)
# Output: phi, shape (n,), with sum(phi) = v(all) - v(empty)

def leverage_shap(n, v, m):
    # Cap the budget at 2^n and round it up to an even number  [alg. 1]
    m = min(m, 2 ** n)
    m = m + m % 2
    # Two evaluations go to the empty and the full coalition  [alg. 1]
    v_empty = v(zeros(n))
    v_all = v(ones(n))
    # Expected rows per size is min(C(n, s), c), with c chosen so the
    # total is m - 2  [eq. 13]  [code: sampling.py#L85-L142]
    expected = expected_rows_per_size(n, m - 2)
    # Integer counts per size, symmetric under s <-> n - s
    # [code: sampling.py#L172-L175]
    counts = symmetric_round(expected)
    pairs = []
    for size in range(1, n // 2 + 1):
        # Middle size of even n: a pair holds two subsets of this size,
        # so half as many pairs are drawn
        # [app. J] [code: sampling.py#L189-L195]
        is_middle = 2 * size == n
        pair_count = counts[size] // 2 if is_middle else counts[size]
        # Draw subsets of this size uniformly, without replacement; each one
        # is paired with its complement  [§3.2]  [app. J]
        for z in draw_subsets(n, size, pair_count, is_middle):
            pairs.append((z, complement(z)))
    rows = [row for pair in pairs for row in pair]
    # One evaluation per row: m - 2 in total  [alg. 1]
    values = [v(row) for row in rows]
    # Target: value minus empty value, minus this row's share of the gap
    # [alg. 1]  [code: leverage_shap.py#L49-L51]
    targets = [
        value - v_empty - (v_all - v_empty) * sum(row) / n
        for row, value in zip(rows, values)
    ]
    # Weight: kernel weight w(|z|) over the realised inclusion probability
    # counts[s] / C(n, s)  [eq. 4]  [code: leverage_shap.py#L52-L53]
    weights = [
        kernel_weight(n, sum(row)) / (counts[sum(row)] / comb(n, sum(row)))
        for row in rows
    ]
    # Projection that removes the all-ones direction, so the efficiency
    # constraint holds by construction  [§3, app. K]
    P = identity(n) - ones(n, n) / n
    design = rows @ P
    # Minimum-norm weighted least squares, not an explicit inverse
    # [code: leverage_shap.py#L63-L66]
    x = weighted_min_norm_lstsq(design, targets, weights)
    # Give back each player an equal share of the efficiency gap  [alg. 1]
    return x + (v_all - v_empty) / n

Cost: the model is called m times, and the extra work is O(m·n² + n³) on top of those calls app. J. The paper notes that the model calls usually dominate app. J.

Theory: the main result bounds the regression objective. With m = O(n log(n/δ) + n/(εδ)) evaluations in expectation, the returned vector has objective at most (1 + ε) times the optimum, with probability 1 − δ, for constant δ §1, eq. 5. A corollary converts this to the squared ℓ2 error of φ, with a problem-specific factor γ = ‖Aφ − b‖² / ‖Aφ‖² §4. The paper reports γ only for n ≤ 16 app. F.

Hyperparameters

Results

Author-reported. Setup: eight datasets from the SHAP library, XGBoost model, exact Tree SHAP ground truth, 100 random runs per setting, error is relative squared ℓ2 error §5. Values at m = 10n tab. 1. Q1–Q3 are the first to third quartiles. Leverage SHAP is the best method in every column except ties tab. 1.

Table 1: relative squared ℓ2 error at m = 10n tab. 1

Table 2: error ratio, Leverage SHAP / Optimized Kernel SHAP, by budget tab. 2

A value of 1 means m ≥ 2^n, where both methods are exact tab. 2.

The paper's headline is 50.0% of the error of Optimized Kernel SHAP on average, computed over the 38 entries that are not 1 tab. 2. See Pitfalls for what this average excludes.

Table 3: quartiles of γ over 100 runs, small-n datasets tab. 3

Other results, described in words:

• Sample size fig. 3: Kernel SHAP has the highest error in all settings. Leverage SHAP is below the median of Optimized Kernel SHAP on the small-n datasets. On Correlated and Independent its third quartile is about the median of Optimized Kernel SHAP. On NHANES and Communities, its upper quartile sits at roughly one and a half times that median fig. 3.
• Noise fig. 4: Leverage SHAP is best at low and moderate noise (σ ≤ 0.1) on every dataset. At σ ≥ 0.5 it is close to Optimized Kernel SHAP.
• Ablations tab. 4 app. G: paired sampling reduces error by a factor of about four or more. Without-replacement sampling helps for small n, by a factor between 1.1 and 10. Unpaired Leverage SHAP is much worse; on California at 10n its mean error is 0.0124 against 0.000249 for Leverage SHAP tab. 4.
• Objective error fig. 10, fig. 11 app. H: the regression objective, not φ, is shown. Leverage SHAP or its Binomial variant gives the best results in most settings with m < 2^n.

Paper vs code

Code: official repository rtealwitter/leverageshap, commit 9702c2cc6386e13fca20b9caeb8747dedf253cef, cloned with --depth 1 [paper footnote, §5]. The repository owner is an author of the paper.

• Sampler. The theorem covers the independent sampler of Algorithm 2. The released estimator draws fixed per-size counts, rounded and symmetric, instead §3.2 sampling.py. The paper says so and writes "we believe the guarantees still hold" §3.2.
• Minimum budget. Algorithm 1 takes m ≥ n alg. 1. The code raises any budget below 6 to 6 with a printed warning leverage_shap.py.
• Scale search. Algorithm 1 uses a binary search for c alg. 1, line 2. The code inverts the same piecewise-linear sum in closed form, with the same target m − 2 sampling.py.
• Reweighting. Algorithm 1 weights each row by w(|z|) / min(1, 2c·ℓ_z) alg. 1, line 499. The code weights by w(|z|) divided by the realised per-size inclusion probability counts[s] / C(n, s) leverage_shap.py sampling.py.
• Middle size. Appendix J and the code agree: for even n, the middle size draws half as many pairs from subsets of size n/2 − 1 with the last player always included app. J sampling.py.
• Appendix G text vs appendix J. Appendix G says the Binomial variant draws middle-layer pairs from Binomial(C(n, n/2)/2, p) "as in Algorithm 2" app. G. Appendix J says the opposite, that the draw is from Binomial(C(n−1, s−1), p) "rather than halving" app. J. The ablation code uses C(n−1, s−1) ablations.py. Treat the appendix G sentence as a typo.
• Model calls. Algorithm 1 uses m − 2 calls on sampled rows plus two on v(∅) and v(all) alg. 1. The code matches this total.

Pitfalls

• Theorem and code differ. The guarantee is for the independent sampler of Algorithm 2, not the deterministic per-size counts used in the released code and every reported experiment §3.2. Expect the released estimator to behave as the paper's experiments show, not as the proof guarantees.
• Realised budget. In the deterministic version the number of evaluations equals the rounded budget. In the Binomial variant it is random; the paper reports a mean within 2.7% of min(m, 2^n) app. G.
• Small budgets. Below n evaluations the projected system is rank-deficient. When m < 6 the code uses 6 evaluations instead, with a printed warning leverage_shap.py. Report the realised count, not the requested one.
• Overflow for large n. The released sampler uses integer combinations. When the count is too large for numpy, it falls back to drawing with replacement, with no warning sampling.py. The paper reports that this fallback affects NHANES and Communities in its ablation runs app. G.
• Noise touches the empty and full coalitions. The noise wrapper perturbs every prediction, including v(∅) and v(all), because the edge-case calls go through the same model benchmark.py. The paper says only that the estimators observe v(S) + ζ app. D.
• Numerical sensitivity. AᵀA is singular by construction, so use the minimum-norm least-squares solution, not an inverse leverage_shap.py §3. For small n the objective error can fall below 1 because of numerical instability fig. 10 caption.
• The "50%" average. The headline average is over the 38 entries of Table 2 that are not 1. Those 38 are the budgets below 2^n tab. 2. If the ten entries equal to 1 are included, the average ratio is about 0.60 (my arithmetic from Table 2), a reduction of about 40%.
• Asymptotic claim. The abstract's O(n log n) holds when ε and δ are treated as constants. The theorem allows constant δ, and its bound also has a 1/ε term §1.
• γ is not controlled by the theorem. The corollary's ℓ2 bound scales with γ. The paper measures γ only for n ≤ 16 app. F, so the bound's size is unknown for larger n.
• Baseline choice is not stated. The paper does not say how the baseline is chosen. The code uses column means and a random explicand. Shapley values depend on the baseline through v(S), so another baseline changes the target datasets.py.
• Names. "Kernel SHAP" in the paper's ablations is the paper's own unpaired regression. "Optimized Kernel SHAP" is the SHAP library baseline used in the main comparisons app. G §5.
• Hardness. A multiplicative guarantee on every φ_i is NP-hard in general, which is why the guarantee is on the objective app. C.

Lineage

• Builds on: Kernel SHAP (Lundberg & Lee, 2017), cited as the baseline in §1. No card.
• Compared against: Optimized Kernel SHAP and Permutation SHAP from the SHAP library, and the paper's own Kernel SHAP variants §5 app. G.
• Uses an idea from active learning: leverage-score sampling for regression §1.3.
• Related follow-up named by the paper: a unified framework that covers Kernel SHAP and Leverage SHAP and gives non-asymptotic guarantees for Kernel SHAP, by Chen et al. (2025) §1.3.