We thank the reviewer for their careful review, highlighting typos, and giving us the opportunity to clarify.

We will fix all the typos pointed out by the reviewer in the camera-ready version of the manuscript. We will also include a more formal statement of Theorem 2.2 and improve the readability of Table 1 and Algorithm 1. We thank the reviewer for spotting the typo in Algorithm 1; maxval should have only be defined in the context of Algorithm 2 (sampling without replacement), which is the maximum number of samples before which we introduce the Poisson approximation for the sampling distribution.

(1) Proof Sketch.

In response to the reviewer's suggestion (as well as the suggestion of reviewer 7LtQ), we will now be including a proof sketch in the camera-ready version of our paper. It will be included in the main text if space permits, or if not, at the beginning of Appendix A. We hope this would be helpful in providing an overview of the key ideas used in deriving the main theoretical results of our study.

(2) Regression Reformulation.

The main value in the regression formulation to our paper, is that it turns Shapley value computation into a standard linear algebra problem. This is what enables the efficient implementation of KernelSHAP in the SHAP library. In the present work, it also opens a path that allows us to apply existing theory from randomized numerical linear algebra in order to analyze existing methods like KernelSHAP and unbiased KernelSHAP.

(3) Experimental Section.

The experimental section showcases the ability to estimate Shapley values in high dimensions and it compares the different estimators used in practice. For the datasets and models that we tested, we find that the performance of different methods are comparable, with a slight advantage noticed towards the leverage scores sampling, without replacement, least squares estimator.

(4) Benefits of $\lambda$. The main reason we have introduced $\lambda$ is so that our framework encompasses key past estimators like KernelSHAP, LeverageSHAP, and unbiased KernelSHAP. However, it does raise some questions about the "best choice" of $\lambda$, which we discuss in our response to Reviewer TpwW.

We thank the reviewer for their time in reviewing the paper and for stating their main concern clearly. In the rebuttal below, we will address your questions, while highlighting prescriptions for practical applications.

Weaknesses.

(1) Prescriptions for practical applications.

Based on our theoretical and experimental results, the following points are of relevance to practical applications:

a. We can use the knowledge of how the error scales with the number of samples (given by our theoretical analysis) to get a fairly good estimate of the error (as elaborated in our response to reviewer TpwW). This in turn can be used to prescribe the number of samples needed for a specified error.

b. In practice, the best estimators tend to be those employing leverage score sampling without replacement and using least squares estimator; however, this remains dataset and model dependent.

Questions.

(1) Choice of $\lambda$.

Unfortunately, while our analysis reveals that $\lambda$ should be chosen to try to minimize the size of $b_\lambda$ or $P_{U} b_\lambda$, it is currently unclear to us how $\lambda$ should be chosen a priori, particularly in a model agnostic way. Our numerical experiment suggests that $\lambda = \alpha$ (what is used in KernelSHAP and LeverageSHAP) may be better than $\lambda = 0$ (what is used in unbiased KernelSHAP), but we do not view this as a comprehensive experiment.

More broadly, variance reduction techniques, such as the one proposed in the supplementary material and concurrent work like "Regression-adjusted Monte Carlo Estimators for Shapley Values and Probabilistic Values" seem to be worthy of future study. We expect that the best techniques will be model dependent, but that they may substantially improve the performance of Shapley value estimators.

(2) Adaptive Estimators / Optimizing Bounds with Respect to a Sampling Scheme.

Yes, this is a good idea, and is something we have been working on. As our analysis reveals, the algorithms converge in a predictable way with $m$ (i.e. the convergence is the Monte--Carlo rate $\sim 1/\sqrt{m}$). We have found it works quite well to use the estimate from a larger value of $m$ to approximate the error at some $m_0 \ll m$. As long as $m\gg m_0$, the estimate using $m$ samples is a good proxy for the true solution, relative to the size of the error using $m_0$ samples. We can then estimate the error at step $m$ from the estimate of the error at $m_0$ and the known rate of convergence.

In addition to improving the bounds for a specific distribution, it is possible to use prior information to obtain better/adaptive estimators. Specifically, we can use the previous samples to estimate the vector $b_\lambda$ or $P_{U} b_\lambda$, and adapt our sampling distribution accordingly. The algorithm and analysis for such an approach would be more involved since one would ideally like to reuse the previous samples to estimate the Shapley values. While it may not be possible to improve the worst-case performance with such a strategy, it could potentially lead to better performance in practice.

We thank the reviewer for carefully reviewing our paper and for also highlighting typos. We will fix the typos that the reviewer has pointed out in the camera-ready version. We address the questions/comments of the reviewer below.

Questions.

(1) Proof Sketch.

We agree with the reviewer that it would be nice to include a proof sketch in order to present a high-level picture of how our sample complexity bounds are derived. We will include such a proof sketch in the camera-ready version of our manuscript. It will be included in the main text if space permits, or if not, at the beginning of Appendix A.

(2) Leverage Scores and Modified Row-Norm Sampling Comparison.

As the reviewer pointed out, modified row-norm sampling performs almost as good as leverage scores in the worst case, while they can perform better a factor of $\sqrt{d}$ in some cases. Whether this advantage in the performance is seen in practice depends on the model at hand.

For the models that we tested in the experiments, we generally find that leverage scores do slightly better than modified row-norm sampling. On the other hand, to show that we can, in principle, obtain a large advantage for modified row-norm sampling, we have constructed a toy model in Appendix E, where we explicitly demonstrate a factor of $\sqrt{d}$ advantage of modified row-norm sampling over leverage scores.

The key idea used for constructing the model in Appendix E is that modified row-norm sampling and leverage scores place different importance on subsets $S \subset [d]$ of features, and is determined solely by the size of the subsets. Specifically, the weights on subsets of small and large size (compared to $d$) is smaller for the modified row-norm sampling distribution as compared to leverage scores (see Theorem A.6). If we have models for which the entries of the vector $b_\lambda$ or $P_{U} b_\lambda$ corresponding to subsets of moderate size is small, we can expect to see an advantage of modified row-norm sampling over leverage scores. As to whether there are models that are used in practice (such as neural networks) where such a condition holds, and consequently, gives an advantage for modified row-norm sampling over leverage scores, needs to be studied in more detail in the future. We note that a similar discussion holds for comparison of kernel weights with leverage scores, and can be found in Appendix E. We can elaborate on these details in the camera-ready version of the manuscript, and also highlight the open question of whether an advantage of modified row-norm sampling over leverage scores can be obtained for models used in practice.

(3) Experiments Beyond Tree Models.

In the present version, estimates are limited to XGBoost since it is possible to compute exactly the true Shapley values for the model using TreeSHAP algorithm.

It is possible to compute results for neural networks too, but only when the true Shapley values can be computed efficiently ($d\approx 15$ and smaller). Below, we present experimental results for a neural network trained on the Adult dataset. We measure the error for both least squares and matrix-vector estimators, and samples drawn with and without replacement using the leverage score distribution.

We find the results to be consistent with our experiments on trees, with least squares approximation generally performing better than matrix vector estimator.

We thank the reviewer for their time in carefully reviewing our paper. We answer the reviewer's questions below:

(1) Implications.

The KernelSHAP method is widely used to approximate Shapley values (e.g. the SHAP library has 24k stars on Github). However, until our work, there were no strong theoretical bounds guaranteeing the accuracy of KernelSHAP. While the method was observed to work well in practice, the lack of theoretical guarantees represented a major gap in the explainability pipeline. By providing theoretical guarantees for the method, our work further justifies the use of such libraries for providing model explainability in critical fields.

Moreover, by approximating values of $||b_\lambda||$ and $||Ub_\lambda||$, it is possible to estimate the number of samples that will provide a specific $\epsilon, \delta$ error bound. In future work, these bounds could be made precise by analyzing the behavior of $||b_\lambda||$ as the number of samples varies. This also opens up the possibility of adaptively choosing the distribution based on the samples that we have seen so far so as to improve the performance (which was also highlighted by reviewer TpwW), and can be explored in future work.

(2) Comparison Between Estimator Performance (sampling with/without replacement).

Based on our theoretical analysis for sampling with and without replacement in the unpaired case, we expect sampling without replacement to perform as well as or better than sampling with replacement. Experimentally, for the case of paired sampling (where we always sample a subset along with its complement), we find that the sampling with and without replacement perform very similarly for the least squares estimator. However, sampling with replacement outperforms sampling without replacement for matrix multiplication estimator in the paired case: we leave the problem developing bounds for the paired case for future research.

(3) Fine-Grained Results.

We believe the information that the reviewer is looking for can be found in Appendix H (Figures 6 and 7). Indeed, since these are Shapley estimators, we do not expect the curves to vary dramatically; we verify that indeed they converge to a single curve. In this sense, sample efficiency positively impacts faithfulness.