Provably Accurate Shapley Value Estimation via Leverage Score Sampling

## Abstract

Abstract Originally introduced in game theory, Shapley values have emerged as a central tool in explainable machine learning, where they are used to attribute model predictions to specific input features. However, computing Shapley values exactly is expensive: for a general model with n 𝑛 n italic_n features, O ⁢ ( 2 n ) 𝑂 superscript 2 𝑛 O(2^{n}) italic_O ( 2 start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT ) model evaluations are necessary. To address this issue, approximation algorithms are widely used. One of the most popular is the Kernel SHAP algorithm, which is model agnostic and remarkably effective in practice. However, to the best of our knowledge, Kernel SHAP has no strong non-asymptotic complexity guarantees. We address this issue by introducing Leverage SHAP , a light-weight modification of Kernel SHAP that provides provably accurate Shapley value estimates with just O ⁢ ( n ⁢ log ⁡ n ) 𝑂 𝑛 𝑛 O(n\log n) italic_O ( italic_n roman_log italic_n ) model evaluations. Our approach takes advantage of a connection between Shapley value estimation and agnostic active learning by employing leverage score sampling , a powerful regression tool. Beyond theoretical guarantees, we show that Leverage SHAP consistently outperforms even the highly optimized implementation of Kernel SHAP available in the ubiquitous SHAP library [Lundberg & Lee, 2017].

## 1 Introduction

While AI is increasingly deployed in high-stakes domains like education, healthcare, finance, and law, increasingly complicated models often make predictions or decisions in an opaque and uninterpretable way.
In high-stakes domains, transparency in a model is crucial for building trust.
Moreover, for researchers and developers, understanding model behavior is important for identifying areas of improvement and applying appropriate safe guards.
To address these challenges, Shapley values have emerged as a powerful game-theoretic approach for interpreting even opaque models (Shapley, ; Štrumbelj & Kononenko, ; Datta et al., ; Lundberg & Lee, ). These values can be used to effectively quantify the contribution of each input feature to a model’s output, offering at least a partial, principled explanation for why a model made a certain prediction.

Concretely, Shapley values originate from game-theory as a method for determining fair ‘payouts’ for a cooperative game involving $n$ players. The goal is to assign higher payouts to players who contributed more to the cooperative effort. Shapley values quantify the contribution of a player by measuring how its addition to a set of other players changes the value of the game. Formally, let the *value function* $v:2^{[n]}\to\mathbb{R}$ be a function defined on sets $S\subseteq[n]$.
The Shapley value for player $i$ is:

$$ $\displaystyle\phi_{i}$ $\displaystyle=\frac{1}{n}\sum_{S\subseteq[n]\setminus\{i\}}\frac{v(S\cup\{i\}) -v(S)}{\binom{n-1}{|S|}}.$ (1) $$

The denominator weights the marginal contribution of player $i$ to set $S$ by the number of sets of size $|S|$, so that the marginal contribution to sets of each size are equally considered. With this weighting, Shapley values are known to be the unique values that satisfy four desirable game-theoretic properties: Null Player, Symmetry, Additivity, and Efficiency (Shapley, ). For further details on Shapley values and their theoretical motivation, we refer the reader to Molnar ().

A popular way of using Shapley values for explainable AI is to attribute predictions made by a model $f:\mathbb{R}^{n}\to\mathbb{R}$ on a given input $\mathbf{x}\in\mathbb{R}^{n}$ compared to a baseline input $\mathbf{y}\in\mathbb{R}^{n}$ (Lundberg & Lee, ).
The players are the features and $v(S)$ is the prediction of the model when using the features in $S$; i.e., $v(S)=f(\mathbf{x}^{S})$ where $\mathbf{x}^{S}_{i}=x_{i}$ for $i\in S$ and $\mathbf{x}^{S}_{i}=y_{i}$ otherwise.(^1^11There are multiple ways to define the baseline $\mathbf{y}$. The simplest is to consider a fixed vector (Lundberg & Lee, ). Other approaches define $\mathbf{y}$ as a random vector that is drawn from a data distribution (Lundberg et al., ; Janzing et al., ) and take $v(S)=\mathbb{E}[\mathbf{x}^{S}]$, where the expectation is taken over the random choice of $\mathbf{y}$. The focus of our work is estimating Shapley values once the values function $v$ is fixed, so our methods and theoretical analysis are agnostic to the specific approach used.)
When $v$ is defined in this way, Shapley values measure how each feature value in the input contributes to the prediction.

Shapley values also find other applications in machine learning, including in feature and data selection. For feature selection, the value function $v(S)$ is taken to be the model loss when
when using the features in $S$ (Marcílio & Eler, ; Fryer et al., ).
For data selection, $v(S)$ is taken to be the the loss when using the data observations in $S$ (Jia et al., ; Ghorbani & Zou, ).

### 1.1 Efficient Shapley Value Computation

Naively, computing all $n$ Shapley values according to Equation 1 requires $O(2^{n})$ evaluations of $v$ (each of which involves the evaluation of a learned model) and $O(2^{n})$ time. This cost can be reduced in certain special cases, e.g. when computing feature attributions for linear models or decision trees
(Lundberg et al., ; Campbell et al., ; Amoukou et al., ; Chen et al., ).

More often, when $v$ is based on an arbitrary model, like a deep neural network, the exponential cost in $n$ is avoided by turning to approximation algorithms for estimating Shapley values, including sampling, permutation sampling, and Kernel SHAP (Strumbelj & Kononenko, ; Lundberg & Lee, ; Mitchell et al., ). The Kernel SHAP method is especially popular, as it performs well in practice for a variety of models, requiring just a small number of black-box evaluations of $v$ to obtain accurate estimates to $\phi_{1},\ldots,\phi_{n}$. The method is a corner-stone of the ubiquitous SHAP library for explainable AI based on Shapley values (Lundberg & Lee, ).

Kernel SHAP is based on an elegant connection between Shapley values and least squares regression (Charnes et al., ).
Specifically, let $[n]$ denote $\{1,\ldots,n\}$, $\emptyset$ denote the empty set, and $\mathbf{1}$ denote an all $1$’s vector of length $n$. The Shapley values $\bm{\phi}=[\phi_{1},\ldots,\phi_{n}]\in\mathbb{R}^{n}$ are known to satisfy:

$$ $\displaystyle\bm{\phi}=\operatorname*{arg\,min}_{\mathbf{x}:\langle\mathbf{x}, \mathbf{1}\rangle=v([n])-v(\emptyset)}\|\mathbf{Zx-y}\|_{2},$ (2) $$

where $\mathbf{Z}\in\mathbb{R}^{2^{n}-2\times n}$ is a specific structured matrix whose rows correspond to sets $S\subseteq[n]$ with $0<|S|<n$, and $\mathbf{y}\in\mathbb{R}^{2^{n}-2}$ is vector whose entries correspond to values of $v(S)$.
(We precisely define $\mathbf{Z}$ and $\mathbf{y}$ in Section 2.)

Since solving the regression problem in Equation 2 directly would require evaluating $v(S)$ for all $2^{n}-2$ subsets represented in $\mathbf{y}$,
Kernel SHAP solves the problem approximately via *subsampling*. Concretely, for a given number of samples $m$ and a discrete probability distribution $\mathbf{p}\in[0,1]^{2^{n}-2}$ over rows in $\mathbf{Z}$, consider a sampling matrix $\mathbf{S}\in\mathbb{R}^{m\times 2^{n}-2}$, where each row of $\mathbf{S}$ is 0 0 except for a single entry $1/\sqrt{p_{j}}$ in the $j^{\text{th}}$ entry with probability $p_{j}$.
The Kernel SHAP estimate is given by

$$ $\displaystyle\tilde{\bm{\phi}}=\operatorname*{arg\,min}_{\mathbf{x}:\langle \mathbf{x},\mathbf{1}\rangle=v([n])-v(\emptyset)}\|\mathbf{SZx-Sy}\|_{2}.$ (3) $$

Importantly, computing this estimate only requires at most $m$ evaluations of the value function $v$, since $\mathbf{Sy}$ can be constructed from observing at most $m$ entries in $\mathbf{y}$. Value function evaluations typically dominate the computational cost of actually solving the regression problem in Equation 3, so ideally $m$ is chosen as small as possible.
In an effort to reduce sample complexity, Kernel SHAP does not sample rows uniformly. Instead, the row corresponding to subset $S$ is chosen with probability proportional to:

$$ $\displaystyle w(|S|)=\left(\binom{n}{|S|}|S|(n-|S|)\right)^{-1}.$ (4) $$

The specific motivation for this distribution is discussed further in Section 2, but the choice is intuitive: the method is more likely to sample rows corresponding to subsets whose size is close to 0 0 or $n$, which aligns with the fact that these subsets more heavily impact the Shapley values (Equation 1). Ultimately, however, the choice of $w(|S|)$ is heuristic.

In practice, the Kernel SHAP algorithm is implemented with further optimizations (Lundberg & Lee, ; Covert & Lee, ; Jethani et al., ).

- •
Paired Sampling: If a row corresponding to set $S$ is sampled, the row corresponding to its complement $\bar{S}\vcentcolon=[n]\setminus S$ is also sampled. This paired sampling strategy intuitively balances samples so each player $i$ is involved in the exact same number of samples. Paired sampling substantially improves Kernel SHAP (Covert & Lee, ).
- •
Sampling without Replacement: When $n$ is small, sampling with replacement can lead to inaccurate solutions even as the number of samples $m$ approaches and surpasses $2^{n}$, even though $2^{n}$ evaluations of $v$ are sufficient to exactly recover the Shapley values. In the SHAP library implementation of Kernel SHAP, this issue is addressed with a version of sampling without replacement: if there is a sufficient ‘sample budget’ for a given set size, all sets of that size are sampled (Lundberg & Lee, ). Doing so substantially improves Kernel SHAP as $m$ approaches $2^{n}$.

### 1.2 Our Contributions

Despite its ubiquity in practice, to the best of our knowledge, no non-asymptotic theoretical accuracy guarantees are known for Kernel SHAP when implemented with $m<2^{n}-2$ row samples (which corresponds to $m$ evaluations of the value function, $v$). We address this issue by proposing Leverage SHAP, a light-weight modification of Kernel SHAP that 1) enjoys strong theoretical accuracy guarantees and 2) consistently outperforms Kernel SHAP in experiments.

Leverage SHAP begins with the observation that a nearly optimal solution to the regression problem in Equation 2 can be obtained by sampling just $\tilde{O}(n)$ rows with probability proportional to their *statistical leverage scores*, a natural measure for the “importance” or “uniqueness” of a matrix row (Sarlós, ; Rauhut & Ward, ; Hampton & Doostan, ; Cohen & Migliorati, ).(^2^22Technically, this fact is known for *unconstrained* least square regression. Some additional work is needed to handle the linear constraint in Equation 2, but doing so is relatively straightforward.)
This fact immediately implies that, *in principal*, we should be able to provably approximate $\bm{\phi}$ with a nearly-linear number of value function evaluations (one for each sampled row).

However, leverage scores are expensive to compute, naively requiring at least $O(2^{n})$ time to write down for a matrix like $\mathbf{Z}$ with $O(2^{n})$ rows. Our key observation is that this bottleneck can be avoided in the case of Shapley value estimation: we prove that the leverage scores of $\mathbf{Z}$ have a simple closed form that admits efficient sampling without ever writing them all down. Concretely, we show that the leverage score of the row corresponding to any subset $S\subset[n]$ is proportional to $\binom{n}{|S|}^{-1}$.

This suggests a similar, but meaningfully different sampling distribution than the one used by Kernel SHAP (see Equation 4).
Since all subsets of a given size have the same leverage score, we can efficiently sample proportional to the leverage scores by sampling a random size $s$ uniformly from $\{1,\ldots,n-1\}$ then selecting a subset $S$ of size $s$ uniformly at random.

###### Theorem 1.1 .

For any $\epsilon>0$ and constant $\delta>0$, the Leverage SHAP algorithm uses $m=O(n\log n+n/\epsilon)$ evaluations of $v$ in expectation(^3^33Formally, if we include the dependence on $\delta$, we require $m=O(n\log(n/\delta)+n/(\epsilon\delta))$.) and $O(mn^{2})$ additional runtime to return estimated Shapley values $\tilde{\bm{\phi}}$ satisfying $\langle\tilde{\bm{\phi}},\mathbf{1}\rangle=v([n])-v(\emptyset)$ and, with probability $1-\delta$,

$$ $\displaystyle\|\mathbf{Z}\tilde{\bm{\phi}}-\mathbf{y}\|_{2}^{2}\leq(1+\epsilon )\|\mathbf{Z}\bm{\phi}-\mathbf{y}\|_{2}^{2}.$ (5) $$

In words, Theorem 1.1 establishes that, with a near-linear number of function evaluations, we can compute approximate Shapley values whose *objective value* is close to the true Shapley values. We also require $O(mn^{2})$ additional runtime to solve the linearly constrained regression problem in Equation 3, although this cost can be reduced in practice using, e.g. iterative methods.

By leveraging the fact that $\mathbf{Z}$ is a well-conditioned matrix, the bound in Equation 5 also implies a bound on the average squared error, $\|\tilde{\bm{\phi}}-\bm{\phi}\|_{2}^{2}$, which is provided in Section 4.

Beyond our theoretical results, we also show that leverage score sampling can be naturally combined with paired sampling, without sacrificing theoretical guarantees. As for Kernel SHAP, doing so improves its performance in experiments. Moreover, we show that a natural “without replacement” version of leverage score sampling (which is even simpler than the method employed by Kernel SHAP) leads to additional accuracy improvements when $m$ is large. Overall, these simple optimizations lead to an algorithm that consistently outperforms the optimized version of Kernel SHAP available in the SHAP Library. We illustrate this point in Figure 1.

More extensive experiments are included in Section 5. We find that the improvement of Leverage SHAP over Kernel SHAP is especially substantial in settings where $n$ is large in comparison to $m$ and when we only have access to noisy estimates of the Shapley values. This is often the case in applications to explainable AI where, as discussion earlier, $v(S)$ is an expectation over fucntions involving random baselines, which is estimated via a finte sample.

Figure: Figure 1: Predicted versus true Shapley values for all features in 8 datasets (we use $m=5n$ samples for this experiment). Points near the identity line indicate that the estimated Shapley value is close to its true value. The plots suggest that our Leverage SHAP method is more accurate than the baseline Kernel SHAP algorithm, as well as the optimized Kernel SHAP implementation available in the SHAP library. We corroborate these findings with more experiments in Section 5.
Refer to caption: x1.png

###### Theorem 1.1 .

### 1.3 Related Work

#### Shapley Values Estimation.

As discussed, naively computing Shapley values requires an exponential number of evaluations of $v$. While more efficient methods exist for certain structured value functions (see references in Section 1.1), for generic functions, faster algorithms involve some sort of approximation.
The most direct way of obtaining an approximation to is approximate the summation definition of Shapley values (Equation 1), which involves $O(2^{n})$ terms for each subset of $[n]$, with a random subsample of subsets (Castro et al., ; Strumbelj & Kononenko, ; Štrumbelj & Kononenko, ). The first methods for doing so use a different subsample for each player $i$.

A natural alternative is to try selecting subsets in such a way that allows them to be reused across multiple players (Illés & Kerényi, ; Mitchell et al., ). However, since each term in the summation involves both $v(S)$ and $v(S\cup\{i\})$, it is difficult to achieve high levels of sample reuse when working with the summation definition. One option is to split the sum into two, and separately estimate $\sum_{S\subseteq[n]\setminus\{i\}}v(S\cup\{i\})/\binom{n-1}{|S|}$ and $\sum_{S\subseteq[n]\setminus\{i\}}v(S)/\binom{n-1}{|S|}$. This allows substantial subset reuse, but tends to perform poorly in practice due to high variance in the individual sums Wang & Jia ().
By solving a global regression problem that determines the entire $\bm{\phi}$ vector, Kernel SHAP can be viewed as a more effective way of reusing subsets to obtain Shapley values for all players.

In addition to the sampling methods referenced above, Covert & Lee () propose a modified Kernel SHAP algorithm, prove it is unbiased, and compute its asymptotic variance.
However, they find that the modified version performs worse than Kernel SHAP.
Additionally, we note that some recent work in the explainable AI setting takes advantage of the fact that we often wish to evaluate feature impact for a large number of different input vectors $\mathbf{x}\in\mathbb{R}^{n}$, which each induce their own set function $v$. In this setting, it is possible to leverage information gathered for one input vector to more efficiently compute Shapley values for another, further reducing sample cost (Schwab & Karlen, ; Jethani et al., ).
This setting is incomparable to ours; like Kernel SHAP, Leverage SHAP works in the setting where we have black-box access to a *single* value function $v$.

#### Leverage Scores.

As discussed, leverage scores are a natural measure of importance for matrix rows. They are widely used as importance sampling probabilities in randomized matrix algorithms for problems ranging from regression (Cohen et al., ; Alaoui & Mahoney, ), to low-rank approximation (Drineas et al., ; Cohen et al., ; Musco & Musco, ; Rudi et al., ), to graph sparsification (Spielman & Srivastava, ), and beyond (Agarwal et al., ).

More recently, leverage score sampling has been applied extensively in active learning as a method for selecting data examples for training (Cohen & Migliorati, ; Avron et al., ; Chen & Price, ; Erdélyi et al., ; Chen & Derezinski, ; Gajjar et al., ; ; Cardenas et al., ; Shimizu et al., ). The efficient estimation of Shapley values via the regression problem in Equation 4 can be viewed as an active learning problem, since our primary cost is in observing entries in the target vector $\mathbf{y}$, each of which corresponds to an expensive value function evaluation.

For active learning of linear models with $n$ features under the $\ell_{2}$ norm, the $O(n\log n+n/\epsilon)$ sample complexity required by leverge score sampling is known to be near-optimal. More complicated sampling methods can achieve $O(n/\epsilon)$, which is optimal in the worst-case (Batson et al., ; Chen & Price, ). However, it is not clear how to apply such methods efficiently to an exponentially tall matrix, as we manage to do for leverage score sampling.

## 2 Background and Notation

#### Notation.

Lowercase letters represent scalars, bold lowercase letters vectors, and bold uppercase letters matrices.
We use the set notation $[n]=\{1,\dots,n\}$ and $\emptyset=\{\}$.
We let $\mathbf{0}$ denote the all zeros vector, $\mathbf{1}$ the all ones vector, and $\mathbf{I}$ the identity matrix, with dimensions clear from context.
For a vector $\mathbf{x}$, $x_{i}$ is the $i^{\text{th}}$ entry (non-bold to indicate a scalar).
For a matrix $\mathbf{X}\in\mathbb{R}^{\rho\times n}$, $[\mathbf{X}]_{i}\in\mathbb{R}^{1\times n}$ is the $i^{\text{th}}$ row.
For a vector $\mathbf{x}$, $\|\mathbf{x}\|_{2}=(\sum_{i}x_{i}^{2})^{1/2}$ denotes the Euclidean ($\ell_{2}$) norm and $\|\mathbf{x}\|_{1}=\sum_{i}|x_{i}|$ is the $\ell_{1}$ norm.
For a matrix $\mathbf{X}\in\mathbb{R}^{n\times m}$, $\mathbf{X}^{+}$ denotes the Moore–Penrose pseudoinverse.

#### Preliminaries.

Any subset $S\subseteq[n]$ can be represented by a binary indicator vector $\mathbf{z}\in\{0,1\}^{n}$, and we use
use $v(S)$ and $v(\mathbf{z})$ interchangeably. We construct the matrix $\mathbf{Z}$ and target vector $\mathbf{y}$ appearing in Equation 2 by indexing rows by all $\mathbf{z}\in\{0,1\}^{n}$ with $0<\|\mathbf{z}\|_{1}<n$:

- •
Let $\mathbf{Z}\in\mathbb{R}^{2^{n}-2\times n}$ be a matrix with $[\mathbf{Z}]_{\mathbf{z}}=\sqrt{w(\|\mathbf{z}\|_{1})}\mathbf{z}^{\top}$.
- •
Let $\mathbf{y}\in\mathbb{R}^{2^{n}-2}$ be the vector where $[\mathbf{y}]_{\mathbf{z}}=\sqrt{w(\|\mathbf{z}\|_{1})}(v(\mathbf{z})-v(\mathbf
{0}))$.

Above, $w(s)=\left(\binom{n}{s}s(n-s)\right)^{-1}$ is the same weight function defined in Equation 4.

As discussed, the Kernel SHAP method is based on an equaivalence between Shapley values and the solution of a contrained regression problem involving $\mathbf{Z}$ and $\mathbf{y}$. Formally, we have:

###### Lemma 2.1 (Equivalence (Lundberg & Lee, 2017 ; Charnes et al., 1988 ) ) .

$$ $\displaystyle\bm{\phi}$ $\displaystyle=\operatorname*{arg\,min}_{\mathbf{x}:\langle\mathbf{x},\mathbf{1 }\rangle=v(\mathbf{1})-v(\mathbf{0})}\|\mathbf{Zx-y}\|_{2}^{2}$ (6) $\displaystyle=\operatorname*{arg\,min}_{\mathbf{x}:\langle\mathbf{x},\mathbf{1 }\rangle=v(\mathbf{1})-v(\mathbf{0})}\sum_{\mathbf{z}:0<\|\mathbf{z}\|_{1}<n}w (\|\mathbf{z}\|_{1})\cdot\left(\langle\mathbf{z},\mathbf{x}\rangle-(v(\mathbf{ z})-v(\mathbf{0}))\right)^{2}.$ (7) $$

For completeness, we provide a self-contained proof of Lemma 2.1 in Appendix F.
The form in Equation 7 inspires the heuristic choice to sample sets with probabilities proportional to $w(\|\mathbf{z}\|_{1})$ in Kernel SHAP, as larger terms in the sum should intutively be sampled with higher probability to reduce variance of the estimate. However, as discussed, a more principaled way to approximately solve least squares regression problems via subsampling is to use probabilities proportional to the statistical leverage scores. Formally, these scores are defined as follows:

###### Definition 2.2 (Leverage Scores) .

Consider a matrix $\mathbf{X}\in\mathbb{R}^{\rho\times n}$.
For $i\in[\rho]$, the leverage score of the $i^{\text{th}}$ row $[\mathbf{X}]_{i}\in\mathbb{R}^{1\times n}$ is
$\ell_{i}\vcentcolon=[\mathbf{X}]_{i}(\mathbf{X^{\top}X})^{+}[\mathbf{X}]_{i}^{
\top}.$

###### Lemma 2.1 (Equivalence (Lundberg & Lee, 2017 ; Charnes et al., 1988 ) ) .

###### Definition 2.2 (Leverage Scores) .

## 3 Leverage SHAP

Statistical leverage scores are traditionally used to approximately solve unconstrained regression problems. Since
Shapley values are the solution to a linearly constrained problem, we first reformulate this into an unconstrained problem. Ultimately, we sample by leverage scores of the reformulated problem, which we prove have a simple closed form that admits efficient sampling.

Concretely, we have the following equivalence:

###### Lemma 3.1 (Constrained to Unconstrained) .

Let $\mathbf{P}$ be the projection matrix $\mathbf{I}-\frac{1}{n}\mathbf{11}^{\top}$.
Define $\mathbf{A}=\mathbf{ZP}$ and $\mathbf{b}=\mathbf{y}-\mathbf{Z1}\frac{v(\mathbf{1})-v(\mathbf{0})}{n}.$
Then

$$ $\displaystyle\operatorname*{arg\,min}_{\mathbf{x}:\langle\mathbf{x},\mathbf{1} \rangle=v(\mathbf{1})-v(\mathbf{0})}\|\mathbf{Zx-y}\|_{2}^{2}=\operatorname*{ arg\,min}_{\mathbf{x}}\left\|\mathbf{Ax-b}\right\|_{2}^{2}+\mathbf{1}\frac{v( \mathbf{1})-v(\mathbf{0})}{n}.$ (8) $$

Further, we have that $\min_{\mathbf{x}:\langle\mathbf{x},\mathbf{1}\rangle=v(\mathbf{1})-v(\mathbf{0
})}\|\mathbf{Zx-y}\|_{2}^{2}=\min_{\mathbf{x}}\left\|\mathbf{Ax-b}\right\|_{2}
^{2}.$

Lemma 3.1 is proven in Appendix H. Based on the equivalent, we consider methods that construct a sampling matrix $\mathbf{S}$ and return:

$$ $\displaystyle\tilde{\bm{\phi}}=\operatorname*{arg\,min}_{\mathbf{x}}\|\mathbf{ SAx-Sb}\|_{2}^{2}+\mathbf{1}\frac{v(\mathbf{1})-v(\mathbf{0})}{n}.$ (9) $$

The main question is how to build $\mathbf{S}$. Our Leverage SHAP method does so by sampling with probabilities proportional to the leverage scores of A. Since naively these scores would be intractable to compute, requiring $O(\rho n^{2})$ time, where $\rho=2^{n}-2$ is the number of rows in $\mathbf{A}$, our method rests on the derivation of a simple closed form expression for the leverage scores, which we prove below.

###### Lemma 3.1 (Constrained to Unconstrained) .

### 3.1 Analytical Form of Leverage Scores

###### Lemma 3.2 .

Let $\mathbf{A}$ be as defined in Lemma 3.1. The leverage score of the row in $\mathbf{A}$ with index $\mathbf{z}\in\{0,1\}^{n}$, where $0<\|\mathbf{z}\|_{1}<n$, is equal to
$\ell_{\mathbf{z}}=\binom{n}{\|\mathbf{z}\|_{1}}^{-1}$.

The proof of Lemma 3.2 depends on an explicit expression for the matrix $\mathbf{A^{\top}A}$:

###### Lemma 3.3 .

Let $\mathbf{A}$ be as defined in Lemma 3.1. $\mathbf{A^{\top}A}=\frac{1}{n}\mathbf{P}$.

This fact will also be relevant later in our analysis, as it implies that $\mathbf{A^{\top}A}$ is well-conditioned. Perhaps surprisingly, all of its non-zero singular values are equal to $1/n$.

###### Proof of Lemma 3.3 .

Recall that $\mathbf{A}=\mathbf{ZP}$, so $\mathbf{A}^{\top}\mathbf{A}=\mathbf{P^{\top}Z^{\top}}\mathbf{ZP}$. We start by deriving explicit expressions for the $(i,j)$ entry of
$\mathbf{Z}^{\top}\mathbf{Z}$, denoted by $[\mathbf{Z}^{\top}\mathbf{Z}]_{(i,j)}$, for all $i\in[n]$ and $j\in[n]$. In particular, we can check that $[\mathbf{Z}^{\top}\mathbf{Z}]_{(i,j)}=\sum_{\mathbf{z}\in\{0,1\}^{n}}z_{i}z_{j
}w(\|\mathbf{z}\|_{1}).$
So, when $i\neq j$,

$$ $\displaystyle[\mathbf{Z}^{\top}\mathbf{Z}]_{(i,j)}=\sum_{\mathbf{z}:z_{i},z_{j }=1}w(\|\mathbf{z}\|_{1})$ $\displaystyle=\sum_{s=2}^{n-1}\binom{n-2}{s-2}w(s).$ (10) $$

Let $c_{n}$ denote the above quantity, which does not depend on $i$ and $j$. I.e., $[\mathbf{Z}^{\top}\mathbf{Z}]_{(i,j)}=c_{n}$ for $i\neq j.$

For the case $i=j$, let $j^{\prime}$ be any fixed index not equal to $i$. We have that

$$ $\displaystyle[\mathbf{Z}^{\top}\mathbf{Z}]_{(i,i)}$ $\displaystyle=\sum_{\mathbf{z}:z_{i}=1}w(\|\mathbf{z}\|_{1})=\sum_{\mathbf{z}: z_{i}=1,z_{j^{\prime}}=0}w(\|\mathbf{z}\|_{1})+\sum_{\mathbf{z}:z_{i}=1,z_{j^{ \prime}}=1}w(\|\mathbf{z}\|_{1})$ $\displaystyle=\sum_{s=1}^{n-1}\binom{n-2}{s-1}w(s)+\sum_{s=2}^{n-1}\binom{n-2} {s-2}w(s)=\frac{1}{n}+c_{n}.$ (11) $$

The last equality can be verified via a direct calculation.
By Equations 10 and 3.1, we have that $\mathbf{Z^{\top}Z}=\frac{1}{n}\mathbf{I}+c_{n}\mathbf{1}\mathbf{1}^{\top}$.
So we conclude that:

$$ $\displaystyle\mathbf{A}^{\top}\mathbf{A}=\mathbf{P^{\top}Z^{\top}ZP}$ $\displaystyle=\left(\mathbf{I}-\frac{1}{n}\mathbf{11^{\top}}\right)\left(\frac {1}{n}\mathbf{I}+c_{n}\mathbf{11}^{\top}\right)\left(\mathbf{I}-\frac{1}{n} \mathbf{11^{\top}}\right)=\frac{1}{n}\left(\mathbf{I}-\frac{1}{n}\mathbf{11^{ \top}}\right).$ $$

With Lemma 3.3, we are now ready to analytically compute the leverage scores.

###### Proof of Lemma 3.2 .

The leverage scores of $\mathbf{A}$ are given by $\ell_{\mathbf{z}}=[\mathbf{ZP}]_{\mathbf{z}}\left(\mathbf{A^{\top}A}\right)^{+
}[\mathbf{ZP}]_{\mathbf{z}}^{\top}$.
By Lemma 3.3, we have that $\mathbf{A^{\top}A}=\frac{1}{n}\left(\mathbf{I}-\frac{1}{n}\mathbf{11}^{\top}\right)$ and thus $(\mathbf{A^{\top}A})^{+}=n\mathbf{I}-\mathbf{11}^{\top}$.
Recall that $[\mathbf{Z}]_{\mathbf{z}}=\sqrt{w(\|\mathbf{z}\|_{1})}\mathbf{z}^{\top}$ so $[\mathbf{ZP}]_{\mathbf{z}}=\sqrt{w(\|\mathbf{z}\|_{1})}\mathbf{z}^{\top}\left(
\mathbf{I}-\frac{1}{n}\mathbf{11^{\top}}\right)$. We thus have:

$$ $\displaystyle\ell_{\mathbf{z}}$ $\displaystyle=\sqrt{w(\|\mathbf{z}\|_{1})}\left(\mathbf{z}-\frac{\|\mathbf{z} \|_{1}}{n}\mathbf{1}\right)^{\top}(n\mathbf{I}-\mathbf{11}^{\top})\sqrt{w(\| \mathbf{z}\|_{1})}\left(\mathbf{z}-\frac{\|\mathbf{z}\|_{1}}{n}\mathbf{1}\right)$ $\displaystyle=w(\|\mathbf{z}\|_{1})\|\mathbf{z}\|_{1}(n-\|\mathbf{z}\|_{1})= \binom{n}{\|\mathbf{z}\|_{1}}^{-1}.\qed$ (12) $$

Notice that the leverage scores simply correspond to sampling every row proportional to the number of sets of that size.
In retrospect, this is quite intuitive given the original definition of Shapley values in Equation 1.
As shown in Figure 2, unlike the Kernel SHAP sampling distribution, leverage score sampling ensures that all subset sizes are equally represented in our sample from $\mathbf{A}$.

Figure: Figure 2: Let $S$ be a subset sampled with the Kernel SHAP or Leverage SHAP probabilities. The plots show the distribution of the set size $|S|$ for different $n$. As in in the definition of Shapley values, the leverage score distribution places equal *total weight* on each subset size, contrasting with Kernel SHAP which over-weights small and large subsets.
Refer to caption: x2.png

###### Lemma 3.2 .

###### Lemma 3.3 .

###### Proof of Lemma 3.3 .

###### Proof of Lemma 3.2 .

### 3.2 Our Algorithm

With Lemma 3.2 in place, we are ready to present Leverage SHAP, which is given as Algorithm 1. In addition to leverage score sampling, the algorithm incorporates paired sampling and sampling without replacement, both of which are used in optimized implementations of Kernel SHAP.

For paired sampling, the idea is to sample rows with probabilities proportional to the leverage scores, but in a *correlated way*: any time we sample index $\mathbf{z}$, we also select its complement, $\bar{\mathbf{z}}$, where $\bar{z}_{i}=1-z_{i}$ for $i\in[n]$. Note that, by the symmetry of $\mathbf{A}$’s leverage scores, $\ell_{\mathbf{z}}=\ell_{\bar{\mathbf{z}}}$.
Similarly, $w(\|\mathbf{z}\|_{1})=w(\|\bar{\mathbf{z}}\|_{1})$. Moving forward, let $\mathcal{Z}$ denote the set of pairs $(\mathbf{z},\bar{\mathbf{z}})$ where $0<\|\mathbf{z}\|_{1}<n$.

To perform paired sampling without replacement, we select indices $(\mathbf{z},\bar{\mathbf{z}})$ independently at random with probability $\min(1,c(\ell_{\mathbf{z}}+\ell_{\bar{\mathbf{z}}}))=\min(1,2c\ell_{\mathbf{z}})$, where $c>1$ is an oversampling parameter. All rows that are sampled are included in a subsampled matrix $\mathbf{Z^{\prime}P}$ and reweighted by the inverse of the probability with which they were sampled. The expected number of row samples in $\mathbf{Z^{\prime}P}$ is equal to $\sum_{(\mathbf{z},\bar{\mathbf{z}})}\min(1,2c\ell_{\mathbf{z}})$. We choose $c$ via binary search so that this expectation equals $m-2$, where $m$ is our target number of value function evaluations (two evaluations are reserved to compute $v(\mathbf{1})-v(\mathbf{0}))$. I.e., we choose $c$ to solve the equation:

$$ $\displaystyle m-2=\sum_{s=1}^{n-1}\min\left(1,2c\binom{n}{s}^{-1}\right)\binom {n}{s}.$ (13) $$

Note that our procedure is different from the with-replacement Kernel SHAP procedure described in Section 1.1. Sampling without replacement requires more care to avoid iterating over the exponential number of pairs in $\mathcal{Z}$. Fortunately, we can exploit the fact that there are only $n-1$ different set sizes, and all sets of the same size have the same leverage score. This allows us to determine how many sets of a different size should be sampled (by drawing a binomial random variable), and we can then sample the required number of sets of each size uniformly at random. Ultimately, we can collect $m$ samples in $O(mn^{2})$ time. Details are deferred to Appendix G.

Because we sample independently without replacement, it is possible that the number of function evaluations $|\mathcal{Z}^{\prime}|+2$ used by Algorithm 1 exceeds $m$.
However, because this number is itself a sum of independent random variables, it tightly concentrates around its expectation, $m$.
If necessary, Leverage SHAP can be implemented so that the number of evaluations is deterministically $m$ (i.e., instead of sampling from a Binomial on Line 6 in Algorithm 2, set $m_{s}$ to be the expectation of the Binomial).
Since this version of Leverage SHAP no longer has independence, our theoretical analysis no longer applies.
However, we still believe the guarantees hold even without independence, and leave a formal analysis to future work.

Figure: Algorithm 1 Leverage SHAP

## 4 Theoretical Guarantees

As discussed, Leverage SHAP offers strong theoretical approximation guarantees.
Our main result is Theorem 1.1 where we show that, with $m=O(n\log(n/\delta)+\frac{n}{\epsilon\delta})$ samples in expectation, Algorithm 1 returns a solution $\tilde{\bm{\phi}}$ that, with probability $1-\delta$, satisfies

$$ $\displaystyle\|\mathbf{A}\tilde{\bm{\phi}}-\mathbf{b}\|_{2}^{2}\leq(1+\epsilon )\|\mathbf{A}\bm{\phi}-\mathbf{b}\|_{2}^{2}.$ (14) $$

We prove this guarantee in Appendix A.
The proof modifies the standard analysis of leverage scores for active least squares regression, which requires two components: 1) a *subspace embedding guarantee*, proven with a matrix Chernoff bound applied to a sum of rank-$1$ random matrices, and 2) and an *approximate matrix-multiplication guarantee*, proven with a second moment analysis.
We replace these steps with 1) a matrix Bernstein bound applied to a sum of rank-$2$ random matrices and 2) a block-approximate matrix multiplication guarantee.

While well-motivated by the connection between Shapley values and linear regression, the objective function in Theorem 1.1 is not intuitive.
Instead, we may be interested in the $\ell_{2}$-norm error between the true Shapley values and the estimated Shapley values.
Fortunately, we can use special properties of our problem to show that Theorem 1.1 implies the following corollary on the $\ell_{2}$-norm error.

###### Corollary 4.1 .

Suppose $\tilde{\bm{\phi}}$ satisfies
$\|\mathbf{A}\tilde{\bm{\phi}}-\mathbf{b}\|_{2}^{2}\leq(1+\epsilon)\|\mathbf{A
\bm{\phi}-b}\|_{2}^{2}$.
Let $\gamma=\|\mathbf{A}\bm{\phi}-\mathbf{b}\|_{2}^{2}/\|\mathbf{A}\bm{\phi}\|_{2}^
{2}$.
Then

$$ $\displaystyle\|\bm{\phi}-\tilde{\bm{\phi}}\|_{2}^{2}\leq\epsilon\gamma\|\bm{ \phi}\|_{2}^{2}.$ $$

The proof of Corollary 4.1 appears in Appendix B.
The statement of the corollary is in terms of a problem-specific parameter $\gamma$ which intuitively measures how well the optimal coefficients $\bm{\phi}$ can reach the target vector $\mathbf{b}$.
In Appendix C, we explore how the performance varies with $\gamma$ in practice.
We find (see e.g., Figures 5, 6, 7) that performance does erode as $\gamma$ increases for all regression-based algorithms, suggesting that $\gamma$ is not an artifact of our analysis.

###### Corollary 4.1 .

## 5 Experiments

In the experiments, we evaluate Leverage SHAP and Kernel SHAP based on how closely they align with the ground truth Shapley values.
We run our experiments on eight popular datasets from the SHAP library (Lundberg & Lee, ).
We find that Leverage SHAP can even outperform the highly optimized Kernel SHAP in the SHAP library, especially for large $n$ or when access to the set function is noisy.

Figure: Figure 3: The $\ell_{2}$-norm error between the estimated Shapley values and ground truth Shapley values as a function of sample size. The lines report the median error while the shaded regions encompass the first to third quartile over 100 runs. Besides the setting where the exact Shapley values can be recovered, Leverage SHAP reliably outperforms even Optimized Kernel SHAP: the third quartile of the Leverage SHAP error is roughly the median error of Optimized Kernel SHAP.
Refer to caption: x3.png

Implementation Details.
In order to compute the ground truth Shapley values for large values of $n$, we use a tree-based model for which we can compute the exact Shapley values efficiently using Tree SHAP (Lundberg & Lee, ). All of our code is written in Python and can be found in the supplementary material. (We will make a public repository after the review process.) We use the SHAP library for the Optimized Kernel SHAP implementation, Tree SHAP, and the datasets. We use XGBoost for training and evaluating trees, under the default parameters.

Additional experiments. So that we do not overcrowd the plots, we show performance on Kernel SHAP, the Optimized Kernel SHAP implementation, and Leverage SHAP in the main body.
We also run all of the experiments with the ablated estimators in Appendix D.
The ablation experiments (see e.g., Figure 9) suggest that Bernoulli sampling (sampling without replacement) improves Leverage SHAP for small $n$ and the paired sampling improves Leverage SHAP for all settings.
Because it is more interpretable, we report the $\ell_{2}$-norm error metric to compare the estimated to the ground truth Shapley values. We also report the linear objective error for similar experiments in Appendix E.

Figure 3 plots the performance for Kernel SHAP, the Optimized Kernel SHAP, and Leverage SHAP as the number of samples varies (we set $m=5n,10n,20n,\ldots,160n$).
For each estimator, the line is the median error and the shaded region encompasses the first and third quartile.(^4^44We report the median and quartiles instead of the mean and standard deviation since the mean minus standard deviation could be negative for the small error values in our experiments.)
Kernel SHAP gives the highest error in all settings of $m$ and $n$, pointing to the importance of the paired sampling and sampling without replacement optimizations in both Optimized Kernel SHAP and Leverage SHAP.
Because both methods sample without replacement, Optimized Kernel SHAP and Leverage SHAP achieve essentially machine precision as $m$ approaches $2^{n}$.
For all other settings,
Leverage SHAP reliably outperforms Optimized Kernel SHAP; the third quartile of Leverage SHAP is generally at the median of Optimized Kernel SHAP.
The analogous experiment for the linear objective error depicted in Figure 11 indicates the same findings.

**Table 1: Summary statistics of the $\ell_{2}$-norm error for every dataset. We adopt the Olympic medal convention: gold, silver and bronze cells signify first, second and third best performance, respectively. Except for ties, Leverage SHAP gives the best performance across all settings.**
|  | IRIS | California | Diabetes | Adult | Correlated | Independent | NHANES | Communities |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Kernel SHAP |  |  |  |  |  |  |  |  |
| Mean | 0.00261 | 0.0208 | 15.4 | 0.000139 | 0.00298 | 0.00324 | 0.0358 | 130.0 |
| 1st Quartile | 5.69e-07 | 0.0031 | 3.71 | 1.48e-05 | 0.00166 | 0.00163 | 0.0106 | 33.5 |
| 2nd Quartile | 9.52e-06 | 0.0103 | 8.19 | 3.86e-05 | 0.00249 | 0.00254 | 0.0221 | 53.6 |
| 3rd Quartile | 0.00181 | 0.029 | 20.1 | 0.000145 | 0.00354 | 0.00436 | 0.0418 | 132.0 |
| Optimized Kernel SHAP |  |  |  |  |  |  |  |  |
| Mean | 3.28e-14 | 0.00248 | 2.33 | 1.81e-05 | 0.000739 | 0.000649 | 0.00551 | 21.8 |
| 1st Quartile | 2.12e-14 | 0.000279 | 0.549 | 2.16e-06 | 0.00027 | 0.000187 | 0.000707 | 5.85 |
| 2nd Quartile | 3.55e-14 | 0.00138 | 1.26 | 5.43e-06 | 0.000546 | 0.000385 | 0.0024 | 13.0 |
| 3rd Quartile | 4.22e-14 | 0.0036 | 3.03 | 1.63e-05 | 0.00101 | 0.000964 | 0.00665 | 25.1 |
| Leverage SHAP |  |  |  |  |  |  |  |  |
| Mean | 3.28e-14 | 0.000186 | 0.63 | 5.21e-06 | 0.000458 | 0.000359 | 0.00385 | 14.7 |
| 1st Quartile | 2.12e-14 | 1.91e-05 | 0.0631 | 6.3e-07 | 0.000139 | 9.51e-05 | 0.000333 | 3.6 |
| 2nd Quartile | 3.55e-14 | 8.31e-05 | 0.328 | 2.33e-06 | 0.000376 | 0.000235 | 0.00149 | 8.9 |
| 3rd Quartile | 4.22e-14 | 0.000231 | 0.769 | 7.09e-06 | 0.000617 | 0.000556 | 0.00401 | 15.3 |

Table 1 depicts the estimator performance for $m=10n$.
The table shows the mean, first quartile, second quartile, and third quartile of the error for each estimator on every dataset.
Except for ties in the setting where $m\geq 2^{n}$, Leverage SHAP gives the best performance followed by Optimized Kernel SHAP and then Kernel SHAP.

Figure: Figure 4: The $\ell_{2}$-norm error between the estimated and ground truth Shapley values as a function of noise in the set function. Leverage SHAP gives the best performance in all settings; the third quartile of the Leverage SHAP generally matches the median of the Optimized Kernel SHAP error.
Refer to caption: x4.png

Beyond sample size, we evaluate how the estimators perform with noisy access to the set functions.
This setting is particularly relevant for Shapley values in explainable AI, since the set function may be an expectation that is approximated in practice.
Figure 4 shows plots performance as noise is added to the set functions:
Instead of observing $v(S)$, the estimators observe $v(S)+\zeta$ where $\zeta\sim\mathcal{N}(0,\sigma^{2})$ is normally distributed with standard deviation $\sigma$ (we set $\sigma=0,5\times 10^{-3},1\times 10^{-2},5\times 10^{-2},\ldots,1$).
As in the sample size experiments, Leverage SHAP meets or exceeds the other estimators, suggesting the utility of the algorithm as a robust estimator in explainable AI.
The analogous experiment for the linear objective error depicted in Figure 12 indicates the same findings.

## 6 Conclusion

We introduce Leverage SHAP, a principled alternative to Kernel SHAP, designed to provide provable accuracy guarantees with nearly linear model evaluations. The cornerstone of our approach is leverage score sampling, a powerful subsampling technique used in regression. To adapt this method for estimating Shapley values, we reformulate the standard Shapley value regression problem and analytically compute the leverage scores of this reformulation. Leverage SHAP efficiently uses these scores to produce provably accurate estimates of Shapley values. Our method enjoys strong theoretical guarantees, which we prove by modifying the standard leverage score analysis to incorporate the empirically-motivated paired sampling and sampling without replacement optimizations. Through extensive experiments on eight datasets, we demonstrate that our algorithm outperforms even the optimized version of Kernel SHAP, establishing Leverage SHAP as a valuable tool in the explainable AI toolkit.

## 7 Acknowledgements

This work was supported by the National Science Foundation under Grant No. 2045590 and Graduate Research Fellowship Grant No. DGE-2234660.
We are grateful to Yurong Liu for introducing us to the problem of estimating Shapley values.

## Appendix A Proof of Approximation Guarantee

In this section, we prove the following theorem.

###### Theorem A.1 .

Let $m=O(n\log(n/\delta)+n\log(1/\delta)/\epsilon)$.
Algorithm 1 produces an estimate $\tilde{\bm{\phi}}$ such that, with probability $1-\delta$,

$$ $\displaystyle\|\mathbf{A}\tilde{\bm{\phi}}-\mathbf{b}\|_{2}^{2}\leq(1+\epsilon )\|\mathbf{A\bm{\phi}-b}\|_{2}^{2}.$ $$

Because of the connection between the constrained and unconstrained regression problems described in Lemma 3.1, the theorem implies the theoretical guarantees in Theorem 1.1.
The time complexity bound in Theorem 1.1 follows from Lemma G.1 in Appendix G.

Consider the following leverage score sampling scheme where rows are selected in blocks.
(We use the word “block” instead of “pair” in the formal analysis for generality.)
Let $\Theta$ be a partition into blocks of equal size (size 2 in our case) where the leverage scores within a block are equal.
Block $\Theta_{i}$ is independently sampled with probability $p^{+}_{i}\vcentcolon=\min(1,c\sum_{k\in\Theta_{i}}\ell_{k})$ for a constant $c$.
The constant $c$ is chosen so that the expected number of blocks sampled is $m$, i.e., $\sum_{i}p^{+}_{i}=m$.
Let $m^{\prime}$ be the number of blocks sampled.
Let $\mathbf{S}\in\mathbb{R}^{|\Theta_{i}|m^{\prime}\times\rho}$ be the (random) sampling matrix.
Each row of $\mathbf{S}$ corresponds to a row $k$ from a sampled block $i$:
every entry in the row is 0 except for the $k$th entry which is $1/\sqrt{p^{+}_{i}}$.

In order to analyze the solution returned by Algorithm 1, we will prove that the sampling matrix $\mathbf{S}$ preserves the spectral norm and Frobenius norm.

###### Lemma A.2 (Bernoulli Spectral Approximation) .

Let $\mathbf{U}\in\mathbb{R}^{\rho\times n}$ be a matrix with orthonormal columns.
Consider the block random sampling matrix $\mathbf{S}$ described above with rows sampled according to the leverage scores of $\mathbf{U}$.
When $m=\Omega(n\log(n/\delta)/\epsilon^{2})$,

$$ $\displaystyle\|\mathbf{I}-\mathbf{U^{\top}S^{\top}SU}\|_{2}\leq\epsilon$ (15) $$

with probability $1-\delta$.

###### Proof of Lemma A.2 .

We will use the following matrix Bernstein bound (see e.g., Theorem 6.6.1 in Tropp et al. ()).

###### Fact A.3 (Matrix Bernstein) .

Consider a finite sequence $\{\mathbf{X}_{i}\}$ of independent, random, Hermitian matrices with dimension $n$. Assume that $\mathbb{E}[\mathbf{X}_{k}]=\mathbf{0}$ and $\|\mathbf{X}_{i}\|_{2}\leq L$ for all $i$.
Define $\mathbf{X}=\sum_{i}\mathbf{X}_{i}$ and let $V=\|\mathbb{E}[\mathbf{X}^{2}]\|_{2}=\|\sum_{i}\mathbb{E}[\mathbf{X}_{i}^{2}]
\|_{2}$.
Then for any $\epsilon>0$,

$$ $\displaystyle\Pr\left(\|\mathbf{X}\|_{2}\geq\epsilon\right)\leq n\exp\left( \frac{-\epsilon^{2}/2}{V+L\epsilon/3}\right).$ (16) $$

Let $\mathbf{U}_{(k,)}\in\mathbb{R}^{1\times n}$ be the $k$th row of $\mathbf{U}$.
Similarly, let $\mathbf{U}_{\Theta_{i}}\in\mathbb{R}^{|\theta_{i}|\times n}$ be the matrix with rows $\mathbf{U}_{(k,)}$ for $k\in\Theta_{i}$.

We will choose $\mathbf{X}_{i}=\mathbf{U}_{\Theta_{i}}^{\top}\mathbf{U}_{\Theta_{i}}-\frac{1}{
p^{+}_{i}}\mathbf{U}_{\Theta_{i}}^{\top}\mathbf{U}_{\Theta_{i}}\mathbbm{1}[i
\textnormal{ selected}]$.
Then $\mathbb{E}[\mathbf{X}_{i}]=\mathbf{0}$ and $\mathbf{X}=\mathbf{I}-\mathbf{U}^{\top}\mathbf{S}^{\top}\mathbf{SU}$.

First, we will upper bound $\max_{i}\|\mathbf{X}_{i}\|_{2}$.
If $i$ is not selected or $p^{+}_{i}=1$, then $\|\mathbf{X}_{i}\|_{2}=0$ so we only need to consider the case where $i$ is selected and $p^{+}_{i}=c\sum_{k\in\Theta_{i}}\ell_{k}<1$.
Then

$$ $\displaystyle\|\mathbf{X}_{i}\|_{2}\leq\left|1-\frac{1}{p^{+}_{i}}\right|\| \mathbf{U}_{\Theta_{i}}^{\top}\mathbf{U}_{\Theta_{i}}\|_{2}\leq\frac{1}{p^{+}_ {i}}\sum_{k\in\Theta_{i}}\|\mathbf{U}_{(k,)}\|_{2}^{2}=\frac{1}{c}=L$ (17) $$

Next, we will upper bound $\|\sum_{i}\mathbb{E}[\mathbf{X}_{i}^{2}]\|_{2}$.
Again, notice that $\mathbb{E}[\mathbf{X}_{i}^{2}]=\mathbf{0}$ if $p^{+}_{i}=1$ so we only need to consider the case where $p^{+}_{i}<1$.
Then

$$ $\displaystyle\sum_{i}\mathbb{E}[\mathbf{X}_{i}^{2}]=\sum_{i:p^{+}_{i}<1}p^{+}_ {i}\left(1-\frac{1}{p^{+}_{i}}\right)^{2}\left(\mathbf{U}_{\Theta_{i}}^{\top} \mathbf{U}_{\Theta_{i}}\right)^{2}+(1-p^{+}_{i})\left(\mathbf{U}_{\Theta_{i}}^ {\top}\mathbf{U}_{\Theta_{i}}\right)^{2}.$ (18) $$

Notice that $p^{+}_{i}(1-1/p^{+}_{i})^{2}+(1-p^{+}_{i})=p^{+}_{i}-2+1/p^{+}_{i}+1-p^{+}_{i}
=1/p^{+}_{i}-1\leq 1/p^{+}_{i}.$
Therefore

$$ $\displaystyle(\ref{eq:expected_X})\preceq\sum_{i:p^{+}_{i}<1}\frac{1}{p^{+}_{i }}\mathbf{U}_{\Theta_{i}}^{\top}\mathbf{U}_{\Theta_{i}}\mathbf{U}_{\Theta_{i}} ^{\top}\mathbf{U}_{\Theta_{i}}$ (19) $$

Observe that entry $(k,k^{\prime})$ of $\mathbf{U}_{\Theta_{i}}\mathbf{U}_{\Theta_{i}}^{\top}\in\mathbb{R}^{|\Theta_{i
}|\times|\Theta_{i}|}$ is $\mathbf{U}_{k}\mathbf{U}_{k^{\prime}}^{\top}$.
So the absolute value of each entry is $|\mathbf{U}_{k}\mathbf{U}_{k^{\prime}}^{\top}|\leq\|\mathbf{U}_{k}\|_{2}\|
\mathbf{U}_{k^{\prime}}\|_{2}=\ell_{k}^{1/2}\ell_{k^{\prime}}^{1/2}$ by Cauchy-Schwarz.
Define $\ell^{\max}_{i}=\max_{k\in\Theta_{i}}\ell_{k}$ and $\ell^{\min}_{i}=\min_{k\in\Theta_{i}}\ell_{k}$.
By the Gershgorin circle theorem, $\mathbf{U}_{\Theta_{i}}\mathbf{U}_{\Theta_{i}}^{\top}\preceq 2\ell^{\max}_{i}|
\Theta_{i}|^{2}\mathbf{I}$.
Equivalently, $\mathbf{A}\preceq\mathbf{B}$, $\mathbf{x}^{\top}\mathbf{A}\mathbf{x}\leq\mathbf{x}^{\top}\mathbf{B}\mathbf{x}$ for all $\mathbf{x}$.
Consider an arbitrary $\mathbf{z}$.
We have $\mathbf{z^{\top}C^{\top}ACz}\leq\mathbf{z^{\top}C^{\top}BCz}$ since $\mathbf{Cz}$ is some $\mathbf{x}$.
It follows that $\mathbf{U}_{\Theta_{i}}^{\top}\mathbf{U}_{\Theta_{i}}\mathbf{U}_{\Theta_{i}}^{
\top}\mathbf{U}_{\Theta_{i}}\preceq 2\ell^{\max}_{i}|\Theta_{i}|^{2}\mathbf{U}
_{\Theta_{i}}^{\top}\mathbf{U}_{\Theta_{i}}$.
Then

$$ $\displaystyle(\ref{eq:expected_X_bound})\preceq\frac{1}{c}\sum_{i=1}^{|\Theta_ {i}|}\frac{2\ell^{\max}_{i}|\Theta_{i}|^{2}\mathbf{U}_{\Theta_{i}}^{\top} \mathbf{U}_{\Theta_{i}}}{\sum_{k\in\Theta_{i}}\ell_{k}}\preceq\frac{1}{c}\max_ {i}2|\Theta_{i}|\frac{\ell^{\max}_{i}}{\ell^{\min}_{i}}\mathbf{I}.$ (20) $$

Since $|\Theta_{i}|\leq 2$ and the leverage scores in a block are all equal, $\|\mathbb{E}[\mathbf{X}_{j}^{2}]\|_{2}\leq 4/c$.

Recall that $c$ is chosen so that $m=\sum_{i}\min(1,c\sum_{k\in\Theta_{i}}\ell_{k})\leq c\sum_{i}\sum_{k\in\Theta
_{i}}\ell_{k}=cn$, where the last equality follows because the leverage scores sum to $n$.
So $1/c\leq n/m$.

Plugging in to Fact A.3, we have that

$$ $\displaystyle\Pr\left(\|\mathbf{I}-\mathbf{U^{\top}S^{\top}SU}\|_{2}\geq \epsilon\right)\leq n\exp\left(\frac{-m\epsilon^{2}/2n}{4+\epsilon/3}\right)\leq\delta$ (21) $$

provided $m=\Omega(n\log(n/\delta)/\epsilon^{2})$.

∎

Next, we will prove that the sampling matrix preserves the Frobenius norm.

###### Lemma A.4 (Frobenius Approximation) .

Let $\mathbf{U}\in\mathbb{R}^{\rho\times n}$ be a matrix with orthonormal columns.
Consider the block random sampling matrix $\mathbf{S}$ described above with rows sampled according to the leverage scores of $\mathbf{U}$.
Let $\mathbf{V}\in\mathbb{R}^{\rho\times n^{\prime}}$.
As long as $m\geq\frac{1}{\delta\epsilon^{2}}$, then

$$ $\displaystyle\left\|\mathbf{U^{\top}S^{\top}SV}-\mathbf{U^{\top}V}\right\|_{F} \leq\epsilon\|\mathbf{U}\|_{F}\|\mathbf{V}\|_{F}$ (22) $$

with probability $1-\delta$.

###### Proof.

We adapt Proposition 2.2 in Wu () to the setting where blocks are sampled without replacement.

$$ $\displaystyle\mathbb{E}[\|\mathbf{U^{\top}S^{\top}SV}-\mathbf{U^{\top}V}\|_{F} ^{2}]$ $\displaystyle=\sum_{h_{1}=1}^{n}\sum_{h_{2}=1}^{n}\mathbb{E}[(\mathbf{U^{\top} S^{\top}SV}-\mathbf{U^{\top}V})_{(h_{1},h_{2})}^{2}]$ (23) $\displaystyle=\sum_{h_{1}=1}^{n}\sum_{h_{2}=1}^{n}\textrm{Var}[(\mathbf{U^{ \top}S^{\top}SV})_{(h_{1},h_{2})}]$ (24) $$

where the last equality follows because $\mathbb{E}[\mathbf{U^{\top}S^{\top}SV}]=\mathbf{U^{\top}V}$.
We will consider an individual entry in terms of the blocks $i$.
If $p^{+}_{i}=1$, then the entry is 0 so we only need to consider the case where $p^{+}_{i}<1$.
Then, using independence, we have

$$ $\displaystyle\textrm{Var}[(\mathbf{U^{\top}S^{\top}SV})_{(h_{1},h_{2})}]$ $\displaystyle=\sum_{i:p^{+}_{i}<1}\textrm{Var}[\frac{1}{{p^{+}_{i}}}\mathbbm{1 }[i\textrm{ selected}](\mathbf{U}_{(h_{1},\Theta_{i})}^{\top}\mathbf{V}_{( \Theta_{i},h_{2})})]$ (25) $\displaystyle\leq\sum_{i:p^{+}_{i}<1}\frac{1}{p^{+}_{i}}\left(\mathbf{U}_{(h_{ 1},\Theta_{i})}^{\top}\mathbf{V}_{(\Theta_{i},h_{2})})\right)^{2}$ (26) $$

so

$$ $\displaystyle(\ref{eq:variance_entries})$ $\displaystyle\leq\sum_{i:p^{+}_{i}<1}\|\mathbf{U}_{(,\Theta_{i})}^{\top} \mathbf{V}_{(\Theta_{i},)}\|_{F}^{2}\leq\sum_{i:p^{+}_{i}<1}\frac{1}{p^{+}_{i} }\|\mathbf{U}_{(\Theta_{i},)}\|_{F}^{2}\|\mathbf{V}_{(\Theta_{i},)}\|_{F}^{2}$ (27) $\displaystyle=\sum_{i:p^{+}_{i}<1}\frac{\sum_{k\in\Theta_{i}}\|\mathbf{U}_{k} \|_{2}^{2}}{c\sum_{k\in\Theta_{i}}\ell_{k}}\|\mathbf{V}_{(\Theta_{i},)}\|_{F}^ {2}=\frac{1}{c}\|\mathbf{V}\|_{F}^{2}\leq\frac{1}{m}\|\mathbf{U}\|_{F}^{2}\| \mathbf{V}\|_{F}^{2}$ (28) $$

where the last equality follows because $m\leq cn$ and $\|\mathbf{U}\|_{F}^{2}=n$.

By Markov’s inequality,

$$ $\displaystyle\Pr\left(\left\|\mathbf{U^{\top}S^{\top}SV}-\mathbf{U^{\top}V} \right\|_{F}>\epsilon\|\mathbf{U}\|_{F}\|\mathbf{V}\|_{F}\right)$ $\displaystyle\leq\frac{\mathbb{E}\left[\left\|\mathbf{U^{\top}S^{\top}SV}- \mathbf{U^{\top}V}\right\|_{F}^{2}\right]}{\epsilon^{2}\|\mathbf{U}\|_{F}^{2} \|\mathbf{V}\|_{F}^{2}}\leq\frac{1}{m\epsilon^{2}}\leq\delta$ (29) $$

as long as $m\geq\frac{1}{\delta\epsilon^{2}}$.
∎

With the spectral and Frobenius approximation guarantees from Lemmas A.2 and A.4, we are ready to prove Theorem A.1.

###### Proof of Theorem A.1 .

Observe that

$$ $\displaystyle\|\mathbf{A}\tilde{\bm{\phi}}-\mathbf{b}\|_{2}^{2}=\|\mathbf{A} \tilde{\bm{\phi}}-\mathbf{A}\bm{\phi}+\mathbf{A}\bm{\phi}-\mathbf{b}\|_{2}^{2} =\|\mathbf{A}\tilde{\bm{\phi}}-\mathbf{A}\bm{\phi}\|_{2}^{2}+\|\mathbf{A}\bm{ \phi}-\mathbf{b}\|_{2}^{2}$ (30) $$

where the second equality follows because $\mathbf{A}\bm{\phi}-\mathbf{b}$ is orthogonal to any vector in the span of $\mathbf{A}$.
So to prove the theorem, it suffices to show that

$$ $\displaystyle\|\mathbf{A\tilde{\bm{\phi}}-A\bm{\phi}}\|_{2}^{2}\leq\epsilon\| \mathbf{A\bm{\phi}-b}\|_{2}^{2}.$ (31) $$

Let $\mathbf{U}\in\mathbb{R}^{\rho\times n^{\prime}}$ be an orthonormal matrix that spans the columns of $\mathbf{A}$ where $n^{\prime}\leq n$.
There is some $\mathbf{y}$ such that $\mathbf{Uy}=\mathbf{A}\bm{\phi}$ and some $\tilde{\mathbf{y}}$ such that $\mathbf{U}\tilde{\mathbf{y}}=\mathbf{A}\tilde{\bm{\phi}}$.
Observe that $\|\mathbf{A\tilde{\bm{\phi}}-A\bm{\phi}}\|_{2}=\|\mathbf{U}\tilde{\mathbf{y}}-
\mathbf{Uy}\|_{2}=\|\tilde{\mathbf{y}}-\mathbf{y}\|_{2}$ where the last equality follows because $\mathbf{U}^{\top}\mathbf{U}=\mathbf{I}$.

By the reverse triangle inequality and the submultiplicavity of the spectral norm, we have

$$ $\displaystyle\|\tilde{\mathbf{y}}-\mathbf{y}\|_{2}$ $\displaystyle\leq\|\mathbf{U^{\top}S^{\top}SU}(\tilde{\mathbf{y}}-\mathbf{y}) \|_{2}+\|\mathbf{U^{\top}S^{\top}SU}(\tilde{\mathbf{y}}-\mathbf{y})-(\tilde{ \mathbf{y}}-\mathbf{y})\|_{2}$ (32) $\displaystyle\leq\|\mathbf{U^{\top}S^{\top}SU}(\tilde{\mathbf{y}}-\mathbf{y}) \|_{2}+\|\mathbf{U^{\top}S^{\top}SU}-\mathbf{I}\|_{2}\|\tilde{\mathbf{y}}- \mathbf{y}\|_{2}.$ (33) $$

Because $\mathbf{U}$ has the same leverage scores as $\mathbf{A}$ and the number of rows sampled in $\mathbf{S}$ is within a constant factor of $m$, we can apply Lemma A.2: With $m=O(n\log\frac{n}{\delta})$, we have $\|\mathbf{U^{\top}S^{\top}SU}-\mathbf{I}\|_{2}\leq\frac{1}{2}$ with probability $1-\delta/2$.
So, with probability $1-\delta/2$,

$$ $\displaystyle\|\tilde{\mathbf{y}}-\mathbf{y}\|_{2}\leq 2\|\mathbf{U^{\top}S^{ \top}SU}(\tilde{\mathbf{y}}-\mathbf{y})\|_{2}.$ (34) $$

Then

$$ $\displaystyle\|\mathbf{U^{\top}S^{\top}SU}(\tilde{\mathbf{y}}-\mathbf{y})\|_{2}$ $\displaystyle=\left\|\mathbf{U^{\top}S^{\top}}\left(\mathbf{SU}\tilde{\mathbf{ y}}-\mathbf{Sb}+\mathbf{Sb}-\mathbf{SU}\mathbf{y}\right)\right\|_{2}$ (35) $\displaystyle=\left\|\mathbf{U^{\top}S^{\top}S}\left(\mathbf{U}\mathbf{y}- \mathbf{b}\right)\right\|_{2}$ (36) $$

where the second equality follows because $\mathbf{SU}\tilde{\mathbf{y}}-\mathbf{Sb}$ is orthogonal to any vector in the span of $\mathbf{SU}$.
By similar reasoning, notice that $\mathbf{U^{\top}}(\mathbf{Uy-b})=\mathbf{0}$.
Then, as long as $m=O(\frac{n}{\delta\epsilon})$, we have

$$ $\displaystyle\left\|\mathbf{U^{\top}S^{\top}S}\left(\mathbf{U}\mathbf{y}- \mathbf{b}\right)\right\|_{2}\leq\frac{\sqrt{\epsilon}}{2\sqrt{n}}\|\mathbf{U} \|_{F}\|\mathbf{Uy-b}\|_{2}$ (37) $$

with probability $1-\delta/2$ by Lemma A.4.
Since $\mathbf{U}$ has orthonormal columns, $\|\mathbf{U}\|_{F}^{2}\leq n$.
Then, combining inequalities yields

$$ $\displaystyle\|\mathbf{A\tilde{\bm{\phi}}-A\bm{\phi}}\|_{2}^{2}=\|\tilde{ \mathbf{y}}-\mathbf{y}\|_{2}^{2}\leq 2\|\mathbf{U^{\top}S^{\top}SU}(\tilde{ \mathbf{y}}-\mathbf{y})\|_{2}\leq\epsilon\|\mathbf{Uy-b}\|_{2}=\epsilon\| \mathbf{A}\bm{\phi}-\mathbf{b}\|_{2}^{2}$ (38) $$

with probability $1-\delta$.

∎

###### Theorem A.1 .

###### Lemma A.2 (Bernoulli Spectral Approximation) .

###### Proof of Lemma A.2 .

###### Fact A.3 (Matrix Bernstein) .

###### Lemma A.4 (Frobenius Approximation) .

###### Proof.

###### Proof of Theorem A.1 .

## Appendix B Proof of Approximation Corollary

We will establish several properties specific to our problem and then use these properties to prove Corollary 4.1.

###### Lemma B.1 (Properties of $\mathbf{A}$ , $\bm{\phi}$ , and $\bm{\tilde{\phi}}$ ) .

The Shapley values $\bm{\phi}$, the matrix $\mathbf{A}$, and the estimated Shapley values $\bm{\tilde{\phi}}$ produced by Algorithm 1 satisfy

$$ $\displaystyle\|\mathbf{A}(\tilde{\bm{\phi}}-\bm{\phi})\|_{2}^{2}=\frac{1}{n}\| \tilde{\bm{\phi}}-\bm{\phi}\|_{2}^{2}.$ (39) $$

and

$$ $\displaystyle\|\mathbf{A}\bm{\phi}\|_{2}^{2}=\frac{1}{n}\left(\|\bm{\phi}\|_{2 }^{2}-\frac{(v(\mathbf{1})-v(\mathbf{0}))^{2}}{n}\right).$ (40) $$

###### Proof.

Even though $\mathbf{A}$ does not have $n$ singular values, we can exploit the special structure of $\mathbf{A}$, $\bm{\phi}$, and $\tilde{\bm{\phi}}$.
Let $\mathbf{A}=\mathbf{U\Sigma V^{\top}}$ be the singular value decomposition of $\mathbf{A}$.
By Lemma 3.3, we know $\mathbf{\Sigma}\in\mathbb{R}^{(n-1)\times(n-1)}$ is a diagonal matrix with $\frac{1}{\sqrt{n}}$ on the diagonal
and $\mathbf{V}\in\mathbb{R}^{n\times(n-1)}$ has $n-1$ orthonormal columns that are all orthogonal to $\mathbf{1}$.
Let $\mathbf{v}_{1},\ldots,\mathbf{v}_{n-1}\in\mathbb{R}^{n}$ be these vectors.
Then $\mathbf{A}=\sum_{i=1}^{n}\frac{1}{\sqrt{n}}\mathbf{u}_{i}\mathbf{v}_{i}^{\top}$ where $\mathbf{u}_{1},\ldots,\mathbf{u}_{n-1}\in\mathbb{R}^{2^{n-1}}$ are the $n-1$ orthonormal columns in $\mathbf{U}$.

By Lemma 2.1 and 3.1, we can write

$$ $\displaystyle\bm{\phi}=\sum_{i=1}^{n-1}\mathbf{v}_{i}s_{i}+\mathbf{1}\frac{v( \mathbf{1})-v(\mathbf{0})}{n}\qquad\textnormal{and}\qquad\tilde{\bm{\phi}}= \sum_{i=1}^{n-1}\mathbf{v}_{i}\tilde{s}_{i}+\mathbf{1}\frac{v(\mathbf{1})-v( \mathbf{0})}{n}$ (41) $$

for some $s_{i}$ and $\tilde{s}_{i}$ where $i\in[n-1]$.
Then
$(\tilde{\bm{\phi}}-\bm{\phi})=\sum_{i=1}^{n-1}\mathbf{v}_{i}(\tilde{s}_{i}-s_{
i})$
and

$$ $\displaystyle\mathbf{A}(\tilde{\bm{\phi}}-\bm{\phi})=\sum_{i=1}^{n-1}\frac{1}{ \sqrt{n}}\mathbf{u}_{i}\mathbf{v}_{i}^{\top}\sum_{j=1}^{n-1}\mathbf{v}_{j}( \tilde{s}_{j}-s_{j})=\sum_{i=1}^{n-1}\frac{1}{\sqrt{n}}\mathbf{u}_{i}(\tilde{s }_{i}-s_{i})$ (42) $$

so
$\|\mathbf{A}(\tilde{\bm{\phi}}-\bm{\phi})\|_{2}^{2}=\frac{1}{n}\|\tilde{\bm{
\phi}}-\bm{\phi}\|_{2}^{2}.$

Similarly, we can write

$$ $\displaystyle\mathbf{A}\bm{\phi}=\sum_{i=1}^{n-1}\frac{1}{\sqrt{n}}\mathbf{u}_ {i}\mathbf{v}_{i}^{\top}\left(\sum_{j=1}^{n-1}\mathbf{v}_{j}s_{j}+\mathbf{1} \frac{v(\mathbf{1})-v(\mathbf{0})}{n}\right)=\sum_{i=1}^{n-1}\frac{1}{\sqrt{n} }\mathbf{u}_{i}s_{i}$ (43) $$

where the second equality follows because $\mathbf{v}_{1},\ldots,\mathbf{v}_{n-1}$ are orthogonal to $\mathbf{1}$.
Then $\|\mathbf{A}\bm{\phi}\|_{2}^{2}=\frac{1}{n}\left(\|\bm{\phi}\|_{2}^{2}-\frac{(
v(\mathbf{1})-v(\mathbf{0}))^{2}}{n}\right).$

∎

See 4.1

###### Proof of Corollary 4.1 .

We have

$$ $\displaystyle\|\mathbf{A}\tilde{\bm{\phi}}-\mathbf{b}\|_{2}^{2}=\|\mathbf{A} \tilde{\bm{\phi}}-\mathbf{A}\bm{\phi}+\mathbf{A}\bm{\phi}-\mathbf{b}\|_{2}^{2} =\|\mathbf{A}\tilde{\bm{\phi}}-\mathbf{A}\bm{\phi}\|_{2}^{2}+\|\mathbf{A}\bm{ \phi}-\mathbf{b}\|_{2}^{2}$ (44) $$

where the second equality follows because $\mathbf{A}\bm{\phi}-\mathbf{b}$ is orthogonal to any vector in the span of $\mathbf{A}$.

Then, by the assumption, we have

$$ $\displaystyle\|\mathbf{A}\tilde{\bm{\phi}}-\mathbf{A}\bm{\phi}\|_{2}^{2}\leq \epsilon\|\mathbf{A}\bm{\phi}-\mathbf{b}\|_{2}^{2}.$ (45) $$

By the definition of $\gamma=\frac{\|\mathbf{A}\bm{\phi}-\mathbf{b}\|_{2}^{2}}{\|\mathbf{A}\bm{\phi}
\|_{2}^{2}}$ and Lemma B.1,

$$ $\displaystyle\epsilon\|\mathbf{A}\bm{\phi}-\mathbf{b}\|_{2}^{2}=\epsilon\gamma \|\mathbf{A}\bm{\phi}\|_{2}^{2}=\epsilon\gamma\frac{1}{n}\left(\|\bm{\phi}\|_{ 2}^{2}-\frac{(v(\mathbf{1})-v(\mathbf{0}))^{2}}{n}\right).$ (46) $$

Finally, by Lemma B.1 along with Equations 45 and 46, we have

$$ $\displaystyle\frac{1}{n}\|\tilde{\bm{\phi}}-\bm{\phi}\|_{2}^{2}=\|\mathbf{A}( \tilde{\bm{\phi}}-\bm{\phi})\|_{2}^{2}\leq\epsilon\|\mathbf{A}\bm{\phi}- \mathbf{b}\|_{2}^{2}\leq\frac{1}{n}\epsilon\gamma\|\bm{\phi}\|_{2}^{2}.$ (47) $$

The corollary statement then follows after multiplying both sides by $n$.
∎

###### Lemma B.1 (Properties of 𝐀 𝐀 \mathbf{A} bold_A , ϕ bold-italic-ϕ \bm{\phi} bold_italic_ϕ , and ϕ ~ bold-~ bold-italic-ϕ \bm{\tilde{\phi}} overbold_~ start_ARG bold_italic_ϕ end_ARG ) .

###### Proof.

###### Proof of Corollary 4.1 .

## Appendix C Additional Experiments: Performance by γ 𝛾 \gamma italic_γ

The problem-specific term $\gamma=\|\mathbf{A}\bm{\phi}-\mathbf{b}\|_{2}^{2}/\|\mathbf{b}\|_{2}^{2}$ appears in the $\ell_{2}$-norm guarantee of Corollary 4.1.

In Figure 5, we explore experimentally whether $\gamma$ is an artifact of our analysis or a necessary parameter.
We perform the experiment by explicitly building the matrix $\mathbf{A}$ (only possible for small $n$) and then choosing $\mathbf{b}$ as a linear combination of a vector in the column span of $\mathbf{A}$ and a vector not in the column span of $\mathbf{A}$.
Then, we back out the corresponding set function $v$ using the definition of $\mathbf{b}$.
We induce different values of $\gamma$ as we vary how much of $\mathbf{b}$ is a vector in the column span of $\mathbf{A}$.
Because all three of the algorithms we consider perform worse as $\gamma$ increases, the experiment suggests that $\gamma$ is not an artifact of our analysis.
The same trend appears for the ablated estimators, as shown in Figure 6.
Encouragingly, $\gamma$ does not appear to be a relevant factor for the objective error, as shown in Figure 7.

Figure: Figure 5: The $\ell_{2}$-norm error as a function of $\gamma$, the problem-specific parameter in Corollary 4.1. The lines indicate the median and the shaded regions encompass the first and third quartile over 100 runs. All three of the main algorithms we consider have higher $\ell_{2}$-norm error as $\gamma$ grows, suggesting that $\gamma$ is not an artifact of our analysis.
Refer to caption: x5.png

Figure: Figure 6: The $\ell_{2}$-norm error as a function of $\gamma$, the problem-specific parameter in Corollary 4.1. The lines indicate the mean over 100 runs. All algorithms we consider have higher $\ell_{2}$-norm error as $\gamma$ grows, suggesting that $\gamma$ is not an artifact of our analysis.
Refer to caption: x6.png

Figure: Figure 7: The linear objective error as a function of $\gamma$, the problem-specific parameter in Corollary 4.1. The lines indicate the mean over 100 runs. All algorithms we consider have no dependence on $\gamma$ grows, suggesting that $\gamma$ is not a relevant parameter for the linear objective error. This finding aligns with Theorem 1.1, which has no dependence on $\gamma$.
Refer to caption: x7.png

## Appendix D Additional Experiments: ℓ 2 subscript ℓ 2 \ell_{2} roman_ℓ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT -norm error with Ablation

Figure: Figure 8: Predicted versus true Shapley values for all features in 8 datasets (we use $m=5n$ samples for this experiment). Points near the identity line indicate that the estimated Shapley value is close to its true value. The plots suggest that our Leverage SHAP method is more accurate than the baseline Kernel SHAP algorithm, as well as the optimized Kernel SHAP implementation available in the SHAP library and ablated versions of Leverage SHAP.
Refer to caption: x8.png

Figure: Figure 9: The $\ell_{2}$-norm error between the estimated Shapley values and ground truth Shapley values as a function of sample size. The lines report the mean error over 100 runs. Leverage SHAP reliably gives the best performance, exceeding second-best Optimized Kernel SHAP for small $n$ and Leverage SHAP w/o Bernoulli Sampling for large $n$.
Refer to caption: x9.png

Figure: Figure 10: The $\ell_{2}$-norm error between the estimated and ground truth Shapley values as a function of noise in the set function. Leverage SHAP gives the best performance in almost all settings, with the exception of NHANES where Leverage SHAP with replacement (without Bernoulli) gives slightly better performance.
Refer to caption: x10.png

**Table 2: Summary statistics of the $\ell_{2}$-norm error for every dataset and estimator. Except for ties, Leverage SHAP gives the best performance for small $n$. For large $n$, Leverage SHAP with replacement (without Bernoulli) sampling can mildly outperform.**
|  | IRIS | California | Diabetes | Adult | Correlated | Independent | NHANES | Communities |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Kernel SHAP |  |  |  |  |  |  |  |  |
| Mean | 0.00261 | 0.0208 | 15.4 | 0.000139 | 0.00298 | 0.00324 | 0.0358 | 130.0 |
| 1st Quartile | 5.69e-07 | 0.0031 | 3.71 | 1.48e-05 | 0.00166 | 0.00163 | 0.0106 | 33.5 |
| 2nd Quartile | 9.52e-06 | 0.0103 | 8.19 | 3.86e-05 | 0.00249 | 0.00254 | 0.0221 | 53.6 |
| 3rd Quartile | 0.00181 | 0.029 | 20.1 | 0.000145 | 0.00354 | 0.00436 | 0.0418 | 132.0 |
| Optimized Kernel SHAP |  |  |  |  |  |  |  |  |
| Mean | 3.28e-14 | 0.00248 | 2.33 | 1.81e-05 | 0.000739 | 0.000649 | 0.00551 | 21.8 |
| 1st Quartile | 2.12e-14 | 0.000279 | 0.549 | 2.16e-06 | 0.00027 | 0.000187 | 0.000707 | 5.85 |
| 2nd Quartile | 3.55e-14 | 0.00138 | 1.26 | 5.43e-06 | 0.000546 | 0.000385 | 0.0024 | 13.0 |
| 3rd Quartile | 4.22e-14 | 0.0036 | 3.03 | 1.63e-05 | 0.00101 | 0.000964 | 0.00665 | 25.1 |
| Leverage SHAP |  |  |  |  |  |  |  |  |
| Mean | 3.28e-14 | 0.000186 | 0.63 | 5.21e-06 | 0.000458 | 0.000359 | 0.00385 | 14.7 |
| 1st Quartile | 2.12e-14 | 1.91e-05 | 0.0631 | 6.3e-07 | 0.000139 | 9.51e-05 | 0.000333 | 3.6 |
| 2nd Quartile | 3.55e-14 | 8.31e-05 | 0.328 | 2.33e-06 | 0.000376 | 0.000235 | 0.00149 | 8.9 |
| 3rd Quartile | 4.22e-14 | 0.000231 | 0.769 | 7.09e-06 | 0.000617 | 0.000556 | 0.00401 | 15.3 |
| Kernel SHAP Paired |  |  |  |  |  |  |  |  |
| Mean | 0.000117 | 0.000539 | 1.13 | 7.61e-06 | 0.000682 | 0.000621 | 0.00589 | 24.6 |
| 1st Quartile | 4.06e-14 | 4.15e-05 | 0.174 | 9.07e-07 | 0.000209 | 0.000148 | 0.000566 | 6.6 |
| 2nd Quartile | 1.57e-09 | 0.00022 | 0.566 | 2.45e-06 | 0.000539 | 0.000389 | 0.00188 | 13.0 |
| 3rd Quartile | 2.26e-06 | 0.000586 | 1.55 | 8.08e-06 | 0.000957 | 0.00097 | 0.00609 | 25.0 |
| Leverage SHAP wo Bernoulli |  |  |  |  |  |  |  |  |
| Mean | 0.000144 | 0.000478 | 1.05 | 8.07e-06 | 0.000445 | 0.000374 | 0.00293 | 13.6 |
| 1st Quartile | 4.06e-14 | 3.46e-05 | 0.146 | 6.99e-07 | 0.000136 | 9.91e-05 | 0.0003 | 3.38 |
| 2nd Quartile | 3.83e-10 | 0.000206 | 0.433 | 2.84e-06 | 0.000295 | 0.000271 | 0.00109 | 7.62 |
| 3rd Quartile | 3.96e-06 | 0.000492 | 1.29 | 9.13e-06 | 0.00064 | 0.000562 | 0.00329 | 13.3 |
| Leverage SHAP wo Bernoulli, Paired |  |  |  |  |  |  |  |  |
| Mean | 0.0019 | 0.0165 | 13.3 | 0.000106 | 0.00202 | 0.00207 | 0.0207 | 73.1 |
| 1st Quartile | 2.01e-07 | 0.00302 | 2.15 | 1.37e-05 | 0.000999 | 0.00099 | 0.00631 | 17.3 |
| 2nd Quartile | 8.77e-06 | 0.00882 | 5.58 | 4.16e-05 | 0.00166 | 0.00172 | 0.0131 | 30.2 |
| 3rd Quartile | 0.000348 | 0.0226 | 15.9 | 0.000104 | 0.00245 | 0.00278 | 0.0251 | 65.2 |

## Appendix E Additional Experiments: Linear Objective Error with Ablation

Figure: Figure 11: Linear objective error by sample size. The lines indicate the mean error over 100 runs. Leverage SHAP quickly achieves the lowest error in all settings. (For small $n$, numerical instability can lead to error below 1.)
Refer to caption: x11.png

Figure: Figure 12: Linear objective error by standard deviation of noise. The noise grows large as $\sigma$ approaches one for all algorithms. The plots indicate no clear winner across all four datasets.
Refer to caption: x12.png

## Appendix F Equivalence

See 2.1

We provide a cleaner, direct proof using the unconstrained problem and the exact characterization of $\mathbf{Z^{\top}Z}$.

###### Proof of Lemma 2.1 .

We will show that

$$ $\displaystyle\bm{\phi}=\operatorname*{arg\,min}_{\mathbf{x}}\|\mathbf{Ax-b}\|_ {2}^{2}+\mathbf{1}\frac{v(\mathbf{1})-v(\mathbf{0})}{n}$ (48) $$

where
$\mathbf{A}=\mathbf{ZP}$ and $\mathbf{b}=\mathbf{y}-\mathbf{Z1}\frac{v(\mathbf{1})-v(\mathbf{0})}{n}.$
Then Lemma 2.1 follows from Lemma 3.1.

We know

$$ $\displaystyle\operatorname*{arg\,min}_{\mathbf{x}}\|\mathbf{Ax-b}\|_{2}^{2}=( \mathbf{A^{\top}A})^{+}\mathbf{A^{\top}b}.$ (49) $$

By Lemma 3.3, we have $(\mathbf{A^{\top}A})^{+}=(\mathbf{P^{\top}Z^{\top}ZP})^{+}=n\mathbf{I}-\mathbf
{11^{\top}}=n\mathbf{P}$.
Then

$$ $\displaystyle(\mathbf{A^{\top}A})^{+}\mathbf{A}^{\top}=n\mathbf{P}\mathbf{P}^{ \top}\mathbf{Z}^{\top}=n\mathbf{P}\mathbf{Z}^{\top}.$ (50) $$

Let $\mathbf{z}\in\{0,1\}^{n}$ with $0<\|\mathbf{z}\|_{1}<n$.
Then the $\mathbf{z}$th column of $(\mathbf{A^{\top}A})^{+}\mathbf{A}^{\top}$ is given by

$$ $\displaystyle n\sqrt{w(\|\mathbf{z}\|_{1})}\mathbf{z}-\sqrt{w(\|\mathbf{z}\|_{ 1})}\|\mathbf{z}\|_{1}\mathbf{1}.$ (51) $$

Recall $w(s)=\frac{(s-1)!(n-s-1)!}{n!}$.
The weights are carefully designed so that $w(s)s=w(s+1)(n-s-1)$.
Then, for $i\in[n]$, we have the $i$th entry of $(\mathbf{A^{\top}A})^{+}\mathbf{A}^{\top}\mathbf{b}$ given by

$$ $\displaystyle\sum_{\mathbf{z}:0<\|\mathbf{z}\|_{1}<n}w(\|\mathbf{z}\|_{1})(nz_ {i}-\|\mathbf{z}\|_{1})\left(v(\mathbf{z})-v(\mathbf{0})-\|\mathbf{z}\|_{1} \frac{v(\mathbf{1})-v(\mathbf{0})}{n}\right)$ $\displaystyle=\sum_{S\subseteq[n]\setminus\{i\}}\frac{|S|!(n-|S|-1)!}{n!} \biggl{(}v(S\cup\{i\})-v(\emptyset)-(|S|+1)\frac{v([n])-v(\emptyset)}{n}$ $\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\quad-v(S)+v(\emptyset)+(|S|) \frac{v([n])-v(\emptyset)}{n}\biggr{)}$ (52) $$

where the equality follows by converting from binary vectors to sets.
Observe that $z_{i}=1$ if and only if the corresponding set $S\cup\{i\}$ contains $i$.
Continuing, we have

$$ $\displaystyle(\ref{eq:set_conversion})$ $\displaystyle=\sum_{S\subseteq[n]\setminus\{i\}}\frac{|S|!(n-|S|-1)!}{n!} \biggl{(}v(S\cup\{i\})-v(S)-\frac{v([n])-v(\emptyset)}{n}\biggr{)}$ $\displaystyle=\phi_{i}-\frac{v([n])-v(\emptyset)}{n}\sum_{S\subseteq[n] \setminus\{i\}}\frac{|S|!(n-|S|-1)!}{n!}.$ (53) $$

Then Equation 48 follows since
$\sum_{S\subseteq[n]\setminus\{i\}}\frac{|S|!(n-|S|-1)!}{n!}=\frac{1}{n}\sum_{S
\subseteq[n]\setminus\{i\}}\binom{n-1}{|S|-1}^{-1}=1$.

∎

###### Proof of Lemma 2.1 .

## Appendix G Sampling Without Replacement

Since there are an exponential number of rows,
the naive method of considering each row independently will not work.
Our approach is as follows:
For each set size $s\in[\lfloor n/2\rfloor]$, we sample the number of samples that will be chosen.
This is distributed as a Binomial distribution with $\binom{n}{s}$ trials and probability $\min(1,2c\binom{n}{s}^{-1})$.
Then, because the probabilities are the same for all subsets of size $s$, we exploit this symmetry and sample random indices uniformly from $[\binom{n}{s}]$ without replacement.
We index all subsets of size $s$ and then construct the subsets corresponding to the randomly sampled indices without enumerating all subsets.
The Bernoulli sampling code appears in Algorithm 2 and the subset construction code appears in Algorithm 3, both deferred to Appendix G.

Like before, there is one special case when $n$ is even and there is a middle set size $s=\lfloor n/2\rfloor$.
Here, the paired samples $\mathbf{z}$ and $\bar{\mathbf{z}}$ have the same size $s$ so we need to partition the set of all subsets of size $s$ so we do not risk sampling the same pair twice.
We accomplish this by sampling $\mathbf{z}$ from all subsets of $n-1$ items with size $s-1$ then appending $\mathbf{z}_{n}=1$ and computing the complement $\bar{\mathbf{z}}$.

Figure: Algorithm 2 BernoulliSample$(n,c)$: Efficient Bernoulli Sampling

Figure: Algorithm 3 Combo$(n,s,i)$: Compute the $i$th Combination in Lexicographic Order

###### Lemma G.1 (Time Complexity) .

Let $m^{\prime}$ be the number of samples used in Algorithm 1 and define $T_{m^{\prime}}$ as the time complexity of evaluating $v$ in parallel on $m^{\prime}$ inputs.
Algorithm 1 runs in time $O(m^{\prime}n^{2}+T_{m^{\prime}})$.

For most settings, we expect $T_{m^{\prime}}$ to dominate the time complexity.
For example, consider a shallow fully connected neural network with inner dimension just $n$.
Even evaluating the forward pass on $m^{\prime}$ inputs for a single layer of this network takes time $O(m^{\prime}n^{2})$.

###### Proof of Lemma G.1 .

BernoulliSample takes $O(mn^{2})$ time:
For each pair of samples, we call Combo at most twice which takes time at most $O(n^{2})$ since there are only two loops, each over at most $n$ items.
For simplicity of presentation, the algorithm calls $v$ three times: once with input $\mathbf{1}$, once with input $\mathbf{0}$, and once with $m^{\prime}-2$ inputs.
However, by concatenating the inputs, the algorithm can be easily modified to call $v$ once on $m^{\prime}$ inputs for a time complexity of $T_{m^{\prime}}$.
Computing the optimal subsampled solution $\tilde{\bm{\phi}}^{\perp}$ takes $O(n^{3})$ time to compute a factorization of a $n\times n$ matrix and $O(m^{\prime}n^{2})$ time to compute the $n\times n$ matrix.
Since $m^{\prime}\geq n$, the time complexity is $O(m^{\prime}n^{2})$.
∎

###### Lemma G.1 (Time Complexity) .

###### Proof of Lemma G.1 .

## Appendix H Constrained Regression

See 3.1

###### Proof of Lemma 3.1 .

We will decompose $\mathbf{x}$ into one vector that is orthogonal to $\mathbf{1}$ and another vector that is a scaling of $\mathbf{1}.$ That is, $\mathbf{x}=\mathbf{x}^{\prime}+c\mathbf{1}$ where $\langle\mathbf{x^{\prime}},\mathbf{1}\rangle=0$ for some $c$.
In order to satisfy the constraint, it follows that $c=\frac{v(\mathbf{1})-v(\mathbf{0})}{n}$.
Then

$$ $\displaystyle\operatorname*{arg\,min}_{\mathbf{x}:\langle\mathbf{x},\mathbf{1} \rangle=v(\mathbf{1})-v(\mathbf{0})}\|\mathbf{Zx}-\mathbf{y}\|_{2}^{2}$ $\displaystyle=\operatorname*{arg\,min}_{\mathbf{x^{\prime}}+c\mathbf{1}: \langle\mathbf{x^{\prime}},\mathbf{1}\rangle=0}\|\mathbf{Z}(\mathbf{x^{\prime} }+c\mathbf{1})-\mathbf{y}\|_{2}^{2}$ (54) $\displaystyle=\operatorname*{arg\,min}_{\mathbf{x^{\prime}}:\langle\mathbf{x^{ \prime}},\mathbf{1}\rangle=0}\|\mathbf{Zx^{\prime}}-(\mathbf{y}-c\mathbf{Z} \mathbf{1})\|_{2}^{2}+c\mathbf{1}$ (55) $\displaystyle=\operatorname*{arg\,min}_{\mathbf{x}^{\prime}}\|\mathbf{ZPx^{ \prime}}-(\mathbf{y}-c\mathbf{Z}\mathbf{1})\|_{2}^{2}+c\mathbf{1}$ (56) $$

We used $\mathbf{P}$ to project off any component in the direction of $\mathbf{1}$ and thereby remove the constraint.
Plugging in the value of $c$ with the definitions of $\mathbf{A}$ and $\mathbf{b}$ yields the first equation.
The second equation follows by a similar argument.
∎

###### Proof of Lemma 3.1 .

#### Traditional Solution Used in Kernel SHAP

Consider the problem

$$ $\displaystyle\mathbf{x^{*}}=\operatorname*{arg\,min}_{\mathbf{x}:\langle \mathbf{x},\mathbf{1}\rangle=v(\mathbf{1})-v(\mathbf{0})}\|\mathbf{Zx}-\mathbf {y}\|_{2}^{2}.$ (57) $$

For Leverage SHAP, we reformulate the problem to an unconstrained regression problem with Lemma 3.1 and solve it using standard least squares.

In the standard Kernel SHAP implementation, the constrained problem is solved directly using the method of Lagrange multipliers.
We put the analysis in here for completeness.

The Lagrangian function with respect to $\mathbf{x}$ and multiplier $\lambda\in\mathbb{R}$ is

$$ $\displaystyle\mathcal{L}(\mathbf{x},\lambda)$ $\displaystyle=\|\mathbf{Zx}-\mathbf{y}\|_{2}^{2}+\lambda(\langle\mathbf{x}, \mathbf{1}\rangle-v(\mathbf{1})+v(\mathbf{0}))$ (58) $$

At the optimal solution $\mathbf{x}^{*}$ and $\lambda^{*}$, the gradients with respect to $\mathbf{x}$ and $\lambda$ are 0.
This observation implies the following two equations:

$$ $\displaystyle\nabla_{\mathbf{x}}\mathcal{L}(\mathbf{x}^{*},\lambda^{*})=2 \mathbf{Z}^{\top}(\mathbf{Z}\mathbf{x}^{*}-\mathbf{y})+\lambda^{*}\mathbf{1}=0 \quad\implies\quad\mathbf{x}^{*}=(\mathbf{Z^{\top}Z})^{-1}\left(\mathbf{Z}^{ \top}\mathbf{y}-\frac{\lambda^{*}}{2}\mathbf{1}\right)$ (59) $\displaystyle\nabla_{\lambda}\mathcal{L}(\mathbf{x}^{*},\lambda^{*})={\langle \mathbf{x}^{*}},\mathbf{1}\rangle-v(\mathbf{z})+v(\mathbf{0})=0.$ (60) $$

Together, Equations 59 and 60 tell us that

$$ $\displaystyle\langle\mathbf{1},{\mathbf{x}^{*}}\rangle=v(\mathbf{z})-v(\mathbf {0})=\mathbf{1}^{\top}(\mathbf{Z^{\top}Z})^{-1}\left(\mathbf{Z^{\top}b}-\frac{ \lambda^{*}}{2}\mathbf{1}\right)$ (61) $$

Plugging back into Equation 59, we have that

$$ $\displaystyle\mathbf{x}^{*}=(\mathbf{Z^{\top}Z})^{-1}\left(\mathbf{Z^{\top}y}- \frac{\mathbf{1}^{\top}(\mathbf{Z^{\top}Z})^{-1}\mathbf{y}-v(\mathbf{1})+v( \mathbf{0})}{\mathbf{1}^{\top}(\mathbf{Z^{\top}Z})^{-1}\mathbf{1}}\mathbf{1} \right).$ (62) $$

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: 1. It would be nice to have a fuller picture of the Shapley value problem from a broader approximation algorithms standpoint. The linear regression approach is ultimately just one approach to the problem. Do we have a sense for the best possible approximation one can hope for? Are there any computational inapproximability results?
1. Also, it's natural to ask whether additional structure makes the problem easier (e.g. for the specific setting of feature attribution). I realize this is indeed true for decision trees and such. But I am curious if the authors have an idea of a "dream result" under less restrictive but still reasonable structural assumptions, especially for feature attribution. For feature attribution for a general black box model, is the Leverage SHAP approach likely to be the best?
1. One thing that was not totally clear to me from the experiments is whether Leverage SHAP strictly dominates Kernel SHAP at every point in the running time vs accuracy graph. That is, for any fixed running time budget, is it always better to run Leverage SHAP? Theoretically, one concern could be that Leverage SHAP necessarily requires sampling $m = \Omega(n \log n)$ rows (IIUC based on Thm 1.1), whereas I believe Kernel SHAP allows you to pick $m$ arbitrarily (albeit without a guarantee). Thus perhaps for small running time budgets, maybe Kernel SHAP can sometimes be more effective than Leverage SHAP. Or perhaps Kernel SHAP is just better optimized as a practical implementation.
I think Figure 3 (sample size $m$ vs error) nearly answers this question, but my question is whether there is any subtlety in how $m$ translates to actual running time. And in general if there is any catch to the "strictly dominates Kernel SHAP" question.

Minor nit: in two places (lines 132 and 266) there is a typo: "principal" -> "principle".
soundness: 4
strengths: Shapley values are a basic and important topic in interpretable AI and beyond, finding wide application in practice. The problem of efficiently estimating them well is a very well-motivated one. This paper makes a very nice and useful contribution to this problem. The key theoretical insight of analyzing the form of the leverage scores is simple but very clever and elegant, and allows them to make use of a very well-studied toolbox in statistics (although there is still technical work to be done). It immediately feels like the "right" way to solve the problem. The resulting algorithm is theoretically sound, clean, simple, as well as effective in practice. The paper is also very clearly written, with a clear description of all the relevant context as well as clear exposition in general. I did not verify all the proofs in complete detail but they seemed correct to me.
summary: This paper provides a new algorithm for approximating Shapley values with provable guarantees. Shapley values have widespread application in ML as a way of formalizing individual feature contributions to a model's final prediction, although this paper considers the fully general setting in terms of a generic value function $v : 2^{[n]} \to \mathbb{R}$ (where $n$ is the number of features or "players"). 

The algorithm builds on a known way of formulating the Shapley value as the solution of a certain $2^n$-dimensional linear regression problem. The widely used "Kernel SHAP" algorithm approximates the solution by essentially subsampling the rows of the design matrix in a certain way. But a more principled approach is to subsample according to leverage scores, a well-studied concept in statistics. The catch is that naively, leverage score take time polynomial in the size of the design matrix (so $2^{O(n)}$) to compute. The key idea the authors exploit to get around this is that for the specific Shapley design matrix, the leverage scores can actually be written down in a simple closed form.

This allows them to efficiently solve the underlying regression problem and carry over known guarantees from the leverage score toolbox. Specifically, they show that the estimated Shapley values are close to the true Shapley values in a certain sense (both in terms of the optimum achieved and the values themselves). They show further refinements using paired sampling without replacement.

Finally, they show various experiments demonstrating that the resulting algorithm indeed outperforms the previous best Kernel SHAP algorithm in practice.
weaknesses: I do not see any major weaknesses. I do think would be helpful for the authors to discuss the limitations of the Leverage SHAP algorithm a bit more (e.g. does it strictly dominate all prior algorithms?), and provide some context on what still remains open in this space (see below for related questions).

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 2
questions: - Which algorithm in the literature does the Optimized Kernel SHAP in the experiments correspond to?
- Do the theoretical guarantees for leverage SHAP continue to hold when optimizations like paired sampling, and sampling without replacement are not applied? Is there a difference in the guarantees with and without the optimizations?
- In Table 2, Leverage SHAP w/o Bernoulli sampling seems to outperform in some cases. Is there a way to understand this? An analysis or discussion around this would be very useful to better understand the method's behavior.
soundness: 3
strengths: - Estimating Shapley scores accurately and efficiently is an important problem in explainable machine learning. The paper provides a theoretically principled approach for this problem.
- The approach seems to outperform Kernel SHAP and optimized Kernel SHAP baselines in the experiments.
summary: The paper develops a computationally efficient method for approximating Shapely values, which are useful in interpretable machine learning. The proposed method provides a modification of the well-known Kernel SHAP method, with additional non-asymptotic theoretical convergence guarantees and better empirical performance. The authors use a modified sampling without replacement approach to optimize their method and report experiments with ablation studies for the optimizations.
weaknesses: - The main theoretical result (Theorem 1.1) is somewhat unsatisfactory as it does not directly compare the true and estimated Shapely values. The authors address this via Corollary 4.1, but it has a non-intuitive $\gamma$ term which is can be large and makes the approximation guarantees weaker. Are there conditions under which $\gamma$ is guaranteed to be small? This would better help understand the limitations of current theoretical results.
- The experiments could include more baselines like (Jethani et al., 2021), (Mitchell et al., 2022b), and (Yang & Salman, 2019) for a more comprehensive comparison with the state-of-the-art.
- The technical novelty in proving the new theoretical results also appears to be limited. The main result seems similar to the active learning guarantee Theorem 1.1 of Shimizu et al. ICLR 2024, and it is not clear what additional technical insights are needed to develop the current result. A discussion of novel technical insights needed to develop the current result would be helpful.

Typos:
- Line 155, Line192 finte, Line 255 contrained
- References are incorrectly bracketed

---

confidence: 3
contribution: 4
initial_rating: 8
presentation: 4
questions: N.A.
soundness: 4
strengths: This paper is very well written, introduces the context of their work beautifully, and provides both a theoretical and practical contribution to the field.
summary: This paper modifies the Kernel SHAP algorithm, which approximates the Shapley value and model-agnostically quantifies the role of each feature in a model prediction. The Kernel SHAP did not enjoy theoretical convergence guarantees, and the authors propose a modification of this algorithm that offers a theoretical convergence guarantee and similar numerical performance.
weaknesses: A weakness is that it might feel niche, but as a non-specialist in interpretable AI, I cannot judge the importance of the Shapley values. If this information is important, then the author's contribution is quite important because it removes some level of heuristic thanks to their theoretical contribution.

Meta-review: This work considers the important Shapley value estimation problem. Authors show that for specific design matrices, there is a closed-form solution that can be used for leverage score sampling, which reduces the computational cost by orders of magnitude. All reviewers agree that the paper is well written and easy to follow, the proposed algorithm is interesting and significant.

Decision: Accept (Spotlight)