On Least Squares Approximations for Shapley Values and Applications to Interpretable Machine Learning

1. Introduction

The original KernelSHAP paper by Lundberg and Lee [1] from 2017 was a landmark success for interpretable machine learning. Its primary achievement was making the Shapley value [2] — a theoretically optimal but computationally intractable [3,4,5] concept from cooperative game theory [6] — practical for explaining machine learning models. Lundberg and Lee [1] discovered that the exact Shapley values for a prediction could be recovered as the solution to a specially weighted least squares regression problem. By sampling a manageable number of coalitional inputs and solving this approximate regression, KernelSHAP provided a single, consistent metric for feature importance. The accompanying open-source shap library SHAP [1] fueled its widespread use, making it a predominant standard for model interpretation, particularly in the model-agnostic setting.

The follow-up work on Unbiased KernelSHAP by Covert and Lee [7] represents a success in methodological rigor. It addresses the problem that the original KernelSHAP algorithm has not been proven to be unbiased and does not provide uncertainty estimates. While not as publicly visible as the initial breakthrough, this refinement solidified SHAP’s foundation as a robust, statistically sound tool.

According to Lundberg and Lee [1], KernelSHAP was developed without knowledge of prior articles on cooperative game theory like Charnes et al. [8] and Ruiz et al. [9], who decades earlier had already derived the precise weighting scheme that makes a least-squares regression yield the Shapley value when complete data is available. This outlines a fascinating case of parallel innovation across disciplines. Cooperative game theory, with the Shapley value [2] as its most important solution concept, had long been a mature and applied field providing rigorous solutions to real-world problems of fair division and coalitional analysis, from allocating airport landing fees among airlines based on runway use [10] to measuring the voting power of political blocs in a legislature [11,12,13]. Yet, the field of machine learning, increasingly reliant on inscrutable black-box models, had not seen its own predictions as precisely such a system — where input features act as cooperating agents whose joint effort produces a model’s output. Lundberg and Lee’s pivotal contribution [1] was to bridge this conceptual gap. They recognized that the abstract "game" defined by a model’s prediction function was a natural domain for Shapley’s theory. By doing so, they inadvertently reinvented a specific computational tool — the weighted least-squares approximation — from the game theory literature, but with a transformative new purpose: not to analyze economic coalitions, but to explain the inner logic of artificial intelligence.

Whereas KernelSHAP [1] and its unbiased variant [7] are very widely used, the methods from Benati et al. [14] — which are equally based on the least squares formulation of the Shapley value — have received comparatively little attention. This article represents, to our knowledge, their first application in explainable artificial intelligence. We adapt the ideas from Benati et al. [14] as general TU game approximation algorithms, terming their weighted sampling strategy LSS (Least Squares Sampling). Regarding their proposal for stratification, we differentiate S-LSS, an algorithm without sample reuse across strata, from SRS-LSS, an algorithm with sample reuse across strata suggested in Benati et al. [14].

The objective of this paper is neither to introduce a novel Shapley value estimator nor to perform a broad comparison of approximation algorithms. Instead, the contribution of this work is a detailed structural, theoretical, and algorithmic analysis of the methods LSS, S-LSS and SRS-LSS from Benati et al. [14]. We compare the variances of the LSS and S-LSS estimators in detail. We point out how sample reuse across strata may introduce non-zero covariance terms between strata for SRS-LSS. We prove that LSS and UKS approximate the same underlying problem, thereby only differing in their respective sampling strategies. Finally, we test how the unbiased Shapley estimators LSS, S-LSS and SRS-LSS perform in comparison to KernelSHAP and UKS for both classical cooperative games and real-world applications from interpretable machine learning.

The remainder of this article is structured as follows. In Section 2, we summarize basic ideas from cooperative game theory, including linear solution concepts and the Shapley value, and introduce the BShap (Baseline Shapley) model for interpretable machine learning. Section 3 briefly reviews Monte Carlo methods, along with stratified sampling and importance sampling for variance reduction. A summary of results on approximating linear solution concepts by means of importance sampling on the coalition space from our previous work [15] is presented in Section 4. The central developments of this paper are detailed in Section 5. We formally introduce the LSS, S-LSS and SRS-LSS estimators, compare variances of LSS and S-LSS, investigate the emergence of non-zero covariance terms between strata for SRS-LSS and clarify the similarities and differences between LSS and UKS. The empirical performance of the algorithms introduced in Section 5 is analyzed in Section 6 using two types of cooperative games and three real-world explainability scenarios, thereby numerically substantiating our previous analytical claims. The paper concludes in Section 7 with a summary and recommendations.

2. Preliminaries on Cooperative Game Theory, the Shapley Value and the BShap Model for Interpretable Machine Learning

Following the expositions in Chakravarty et al. [6], Peters [16], and Benati et al. [14], this section offers a concise introduction to transferable utility (TU) cooperative games and point-valued solution concepts such as the Shapley value. Additionally, we briefly describe two representative TU games — the airport game and the weighted voting game — and the model from interpretable machine learning which we utilize in our analysis and our numerical experiments in Section 6.

2.1. Transferable Utility Games and Their Characteristic Functions

We investigate transferable utility games (TU games), cooperative games where a coalition’s earnings are expressed as a scalar [6,16]. This implies that the utility can be distributed without restriction among the players.

Following [6,16], a TU game is formally defined as a pair $(N,v)$. Here, $N=\{1,\dots ,n\}$ is the player set and $v:2^N\to \mathbb{R}$ is the characteristic function, which assigns a real value to each coalition $S\subseteq N$ representing its total utility from cooperation. We adopt the standard normalization $v(\varnothing )=0$.

It is frequently convenient to represent a coalition S by an indicator vector $\mathbf{z}\in {\{0,1\}}^n$. We define $\mathbf{z}\left(S\right)={\left[{\mathbb{1}}_{i\in S}\right]}_{i=1}^n$ and its inverse $S\left(\mathbf{z}\right)=\{i\in N\mid z_i=1\}$. This mapping allows functions and operators defined on $2^N$ or ${\{0,1\}}^n$ to be used interchangeably; for instance, $v\left(S\right)$ and $v\left(\mathbf{z}\right)$ denote the same value, and $\left|\mathbf{z}\right|=\left|S\right(\mathbf{z}\left)\right|$.

2.2. Linear Solution Concepts for TU Games

Let $G\left(N\right)$ denote the set of all TU games on player set N. A point-valued solution concept is a function $\alpha :G\left(N\right)\to {\mathbb{R}}^n$ that returns a vector $\alpha (N,v)\in {\mathbb{R}}^n$ for each game $(N,v)\in G\left(N\right)$. Each entry ${\alpha }_i$ quantifies the value assigned to player i. These concepts are employed to evaluate individual influence or to distribute a coalition’s payoff $v\left(S\right)$ among its members $S\subseteq N$.

A key property of a point-valued solution concept $\alpha$ is linearity [14]. This requires that $\alpha$ can be written as

$$\alpha (N,v)=\sum _{S\subseteq N}\mathbf{w}\left(S\right)\odot \nu (S,v),$$

(1)

with ⊙ denoting the Hadamard product of vectors, $\mathbf{w}\left(S\right)={[w_1\left(S\right),\dots ,w_n\left(S\right)]}^{\top }$ standing for weights depending only on S and the player, and $\nu (S,v)={[{\nu }_1(S,v),\dots ,{\nu }_n(S,v)]}^{\top }$ denoting values depending upon S, v, and the player.

2.3. The Shapley Value

The Shapley value [2] represents the most prominent solution concept in cooperative game theory. In recent years, it has been extensively adopted in machine learning and explainable artificial intelligence [1,17,18,19]. Formally, for any player $i\in N$, the Shapley value is the expected marginal contribution ${\mathsf{\Delta }}_i(S,v)=v(S\cup i)-v\left(S\right)$, taken over all sets S of players that precede i in a uniformly random permutation $\pi$ of N, i.e.,

$$\begin{array}{cc}\hfill {\varphi }_i(N,v)& ={\mathbb{E}}_{\pi \sim \mathcal{U}\left(\mathsf{\Pi }\right(N\left)\right)}\left[{\mathsf{\Delta }}_i\left({Pre}_i\left(\pi \right),v\right)\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{n!}}\sum _{\pi \in \mathsf{\Pi }\left(N\right)}{\mathsf{\Delta }}_i\left({Pre}_i\left(\pi \right),v\right)\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & =\sum _{\begin{array}{c}S\subseteq N\\ i\notin S\end{array}}{\displaystyle \frac{\left|S\right|!(n-|S|-1)!}{n!}}{\mathsf{\Delta }}_i(S,v),\hfill \end{array}$$

(3)

where $\mathcal{U}$ stands for a uniform distribution, $\mathsf{\Pi }\left(N\right)$ denotes the set of all permutations of N, and ${Pre}_i\left(\pi \right)$ gives the set of players preceding i in a permutation $\pi \in \mathsf{\Pi }\left(N\right)$.

Let $\varphi ={[{\varphi }_1,\dots ,{\varphi }_n]}^{\top }$ be the vector of Shapley values, and let $\widehat{\varphi }$ represent any approximation, with ${\widehat{\varphi }}_i$ its i-th entry. Where no ambiguity arises, we use ${\varphi }_i$ to mean ${\varphi }_i(N,v)$.

It follows directly from () that the Shapley value belongs to the class of linear solution concepts introduced in SubSection 2.2. The coefficients $\mathbf{w}\left(S\right)$ are specified by ${\nu }_i(S,v)={\mathsf{\Delta }}_i(S,v)$ and

$$w_i\left(S\right)=\left\{\begin{array}{cc}{\displaystyle \frac{\left|S\right|!(n-|S|-1)!}{n!}}\hfill & \mathrm{if}i\notin S\hfill \\ 0\hfill & \mathrm{if}i\in S.\hfill \end{array}\right.$$

(4)

Note that the Shapley value exhibits symmetry, i.e., two players that contribute equally to each coalition have the same Shapley value meaning there holds

$$\mathrm{if}\forall S\subseteq N\setminus \{i,j\},v(S\cup \left\{i\right\})=v(S\cup \left\{j\right\}),\mathrm{then}{\varphi }_i(N,v)={\varphi }_j(N,v),$$

for any game $(N,v)$ and any players $i,j\in N$ with $i\ne j$. It also exhibits efficiency, i.e., for all games $(N,v)$ there holds ${\sum }_{i\in N}{\varphi }_i(N,v)=v\left(N\right)$.

2.4. Airport Games

Airport games, proposed by Littlechild and Thompson [10], address the problem of distributing runway construction costs among players whose aircrafts require different runway lengths. Costs are encoded in a vector $\mathbf{c}={[c_1,\dots ,c_n]}^{\top }$, where $c_i$ corresponds to player i. The characteristic function is

$$v\left(S\right)=\underset{i\in S}{max}{c}_i.$$

(5)

Classified as line-graph maintenance problems, these games have closed-form Shapley values [20], making them ideal benchmark instances for the numerical studies in Section 6.

2.5. Weighted Voting Games

Weighted voting games [6,16,21] represent voting bodies in which players possess different voting weights (e.g., parliamentary seats), collected in $\mathbf{c}={[c_1,\dots ,c_n]}^{\top }$. A coalition S succeeds when ${\sum }_{i\in S}{c}_i\ge C$ for a predetermined quota C. This yields a characteristic function defined by

$$v\left(S\right)=\left\{\begin{array}{cc}0\hfill & \mathrm{if}{\sum }_{i\in S}{c}_i<C\hfill \\ 1\hfill & \mathrm{if}{\sum }_{i\in S}{c}_i\ge C.\hfill \end{array}\right.$$

(6)

In our theoretical and numerical investigations in Section 5 and Section 6, we frequently use a family of weighted voting games with $n=7$ players and

$$\mathbf{c}={[1,2,3,1,1,1,1]}^{\top },C\in \{1,\dots ,10\}.$$

(7)

There exist fast methods for computing point-valued solutions of weighted voting games for large player numbers, for instance dynamic programming algorithms [22,23]. Hence these games serve as excellent test games for our experiments in Section 6.

2.6. Baseline Shapley (BShap) for Interpretable Machine Learning

Our numerical experiments in Section 6 apply the Shapley value framework to attribute a model’s prediction to its input features. Here, the prediction task is interpreted as a cooperative game with features as players. Given a model $\mathsf{o}$ and an input $\mathsf{x}$, the characteristic function $v\left(S\right)$ specifies the value of a feature subset S.

We adopt Baseline Shapley (BShap), a method originally suggested in [1] and formalized by Sundararajan and Najmi [24]. The value $v\left(S\right)$ is defined deterministically relative to a single baseline vector ${\mathsf{z}}^{\prime }$, where features not in S are set to their baseline counterparts:

$$v_{\mathrm{BShap}}\left(S\right)=\mathsf{o}\left([{\mathsf{z}}_S,{\mathsf{z}}_{\overline{S}}^{\prime }]\right)-\mathsf{o}\left({\mathsf{z}}^{\prime }\right)$$

(8)

with $\overline{S}$ denoting the set of features not in S. The subtrahend $\mathsf{o}\left({\mathsf{z}}^{\prime }\right)$ in (8) ensures the normalization $v_{\mathrm{BShap}}(\varnothing )=0$ as introduced in SubSection 2.1. Note that this normalization would not be needed for the algorithms based on marginal contributions investigated in our article [15], but becomes crucial for the algorithms based on least square formulations of the Shapley value discussed in Section 5.

Our specific baseline: In the applications in Section 6, the baseline ${\mathsf{z}}^{\prime }$ is taken as the expected feature vector over the training distribution:

$${\mathsf{z}}^{\prime }=\overline{\mathsf{z}}={\mathbb{E}}_{\mathsf{z}\sim \mathcal{D}}\left[\mathsf{z}\right].$$

with $\mathcal{D}$ representing the training data distribution. This anchors Shapley values to the prediction $\mathsf{o}\left(\overline{\mathsf{z}}\right)$ of an average sample, providing a natural reference point. The resulting attributions thus explain the deviation of $\mathsf{o}\left(\mathsf{x}\right)$ from this baseline.

3. Preliminaries on Monte Carlo Methods for Estimation

In this section, we first introduce the use of Monte Carlo methods for estimating expectations together with two techniques for variance reduction. We draw mainly from the short overviews provided by Rubinstein and Kroese [25] and Botev and Ridder [26] as well as from the textbook by Rubinstein and Kroese [27].

Let $\mathbf{X}$ be a d-dimensional discrete (or continuous) random variable with a sample space $\mathcal{X}\subseteq {\mathbb{R}}^d$ and a probability mass function $q_{\mathbf{X}}$ (or probability density function, respectively), which may be known or unknown. Denote a specific realization of $\mathbf{X}$ by $\mathbf{x}\in \mathcal{X}$. For a function $H:\mathcal{X}\to \mathbb{R}$, we want to estimate the expectation

$$\mu =\mathbb{E}\left[H\right(\mathbf{X}\left)\right].$$

(9)

This expectation is frequently intractable in closed form, either because $q_{\mathbf{X}}$ is unknown (and only i.i.d. samples are given) or because the defining sum or integral is computationally infeasible. Monte Carlo methods circumvent this analytical hurdle through simulation.

We begin with Crude Monte Carlo as the foundational method, then present importance sampling and stratified sampling for variance reduction which we examine in this paper. A finite variance of $H\left(\mathbf{X}\right)$ is assumed throughout this manuscript.

3.1. Crude Monte Carlo Method

Crude Monte Carlo proceeds by drawing an i.i.d. sample ${\mathbf{X}}_1,\dots ,{\mathbf{X}}_{\tau }{\sim }_{\mathrm{iid}}{q}_{\mathbf{X}}$ and averaging the $H\left({\mathbf{X}}_k\right)$ values. This yields the standard estimator [26]:

$$\widehat{\mu }={\displaystyle \frac{1}{\tau }}\sum _{k=1}^{\tau }H\left({\mathbf{X}}_k\right).$$

(10)

According to [26], the estimator is unbiased, i.e. $\mathbb{E}\left[\widehat{\mu }\right]=\mu$, and its variance is

$$Var\left[\widehat{\mu }\right]={\displaystyle \frac{1}{\tau }}Var\left[H\left(\mathbf{X}\right)\right].$$

(11)

3.2. Stratified Sampling for Variance Reduction

Stratified sampling is a well-established Monte Carlo variance reduction technique [26]. It separates the sample space $\mathcal{X}$ into ℓ disjoint strata $\{{\mathcal{X}}_1,\dots ,{\mathcal{X}}_{\ell }\}$ such that ${\bigcup }_{l=1}^{\ell }{\mathcal{X}}_l=\mathcal{X}$ with ${\mathcal{X}}_l\cap {\mathcal{X}}_{l^{\prime }}=\varnothing$ for $l\ne l^{\prime }$. Suppose L is a discrete random variable taking values in $\{1,\dots ,\ell \}$ with known probabilities $p_L\left(l\right)=\mathbb{P}(L=l)=\mathbb{P}(\mathbf{X}\in {\mathcal{X}}_l)$. We can rewrite $\mu$ as

$$\mu ={\mathbb{E}}_{p_L}\left[{\mathbb{E}}_{q_{\mathbf{X}}}[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_L]\right]=\sum _{l=1}^{\ell }{p}_L\left(l\right)\mathbb{E}[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_l]=\sum _{l=1}^{\ell }{p}_L\left(l\right){\mu }_l,$$

with ${\mu }_l=\mathbb{E}[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_l]$ standing for the expectation over the conditional probability distribution of $\mathbf{X}$ given that $\mathbf{X}\in {\mathcal{X}}_l$. Our stratified estimator of $\mu$ becomes

$${\widehat{\mu }}^{st}=\sum _{l=1}^{\ell }{p}_L\left(l\right){\widehat{\mu }}_l,$$

where ${\widehat{\mu }}_l$ denotes the estimated value of $H\left(\mathbf{X}\right)$ in stratum ${\mathcal{X}}_l$, i.e.,

$${\widehat{\mu }}_l={\displaystyle \frac{1}{{\tau }_l}}\sum _{k=1}^{{\tau }_l}H\left({\mathbf{X}}_{l,k}\right),$$

(12)

with the sample ${\mathbf{X}}_{l,1},\dots ,{\mathbf{X}}_{l,{\tau }_l}$ of size ${\tau }_l$ being drawn i.i.d. from the conditional probability distribution of $\mathbf{X}$ given that $\mathbf{X}\in {\mathcal{X}}_l$. To ensure fair comparisons to other Monte Carlo techniques, we always assume ${\sum }_{l=1}^{\ell }{\tau }_l\approx \tau$ (up to rounding errors), where $\tau$ stands for the overall sample budget. The estimator is unbiased [26], i.e.,

$$\mathbb{E}\left[{\widehat{\mu }}^{st}\right]=\mu .$$

Assuming that the sample sizes per stratum are proportionally assigned, i.e., ${\tau }_l=p_L\left(l\right)\tau$, the variance of the crude Monte Carlo estimator (11) has been established as an upper bound [26], i.e.

$$Var\left[{\widehat{\mu }}^{st}\right]\le Var\left[\widehat{\mu }\right],$$

meaning that the stratified estimator ${\widehat{\mu }}^{st}$ never performs worse than its crude counterpart $\widehat{\mu }$.

3.3. Importance Sampling for Variance Reduction

Importance sampling introduces another probability mass (or density) function $p_{\mathbf{X}}$, such that $p_{\mathbf{X}}\left(\mathbf{x}\right)=0\Rightarrow H\left(\mathbf{x}\right)q_{\mathbf{X}}\left(\mathbf{x}\right)=0$. This allows (9) to be rewritten as

$$\mu ={\mathbb{E}}_{p_{\mathbf{X}}}\left[{\displaystyle \frac{q_{\mathbf{X}}\left(\mathbf{X}\right)}{p_{\mathbf{X}}\left(\mathbf{X}\right)}}H\left(\mathbf{X}\right)\right].$$

(13)

A Monte Carlo approximation based on (13) is

$$\widehat{\mu }={\displaystyle \frac{1}{\tau }}\sum _{k=1}^{\tau }{\displaystyle \frac{q_{\mathbf{X}}\left({\mathbf{X}}_k\right)}{p_{\mathbf{X}}\left({\mathbf{X}}_k\right)}}H\left({\mathbf{X}}_k\right),$$

(14)

where ${\mathbf{X}}_1,\dots ,{\mathbf{X}}_{\tau }{\sim }_{\mathrm{iid}}{p}_{\mathbf{X}}$. As shown in [25], this estimator is unbiased, i.e. ${\mathbb{E}}_{p_{\mathbf{X}}}\left[\widehat{\mu }\right]=\mu$, and its variance is analogous to the crude Monte Carlo case

$${Var}_{p_{\mathbf{X}}}\left[\widehat{\mu }\right]={\displaystyle \frac{1}{\tau }}{Var}_{p_{\mathbf{X}}}\left[{\displaystyle \frac{q_{\mathbf{X}}\left(\mathbf{X}\right)}{p_{\mathbf{X}}\left(\mathbf{X}\right)}}H\left(\mathbf{X}\right)\right].$$

(15)

With a properly selected $p_{\mathbf{X}}$, the variance of the importance sampling estimator is less than or equal to that of the crude Monte Carlo estimator. This follows from comparing (11) and (15), yielding

$${\displaystyle \frac{1}{\tau }}{Var}_{q_{\mathbf{X}}}\left[H\left(\mathbf{X}\right)\right]\le \underset{p_{\mathbf{X}}}{min}{\displaystyle \frac{1}{\tau }}{Var}_{p_{\mathbf{X}}}\left[{\displaystyle \frac{q_{\mathbf{X}}\left(\mathbf{X}\right)}{p_{\mathbf{X}}\left(\mathbf{X}\right)}}H\left(\mathbf{X}\right)\right],$$

which holds because $p_{\mathbf{X}}=q_{\mathbf{X}}$ can always be chosen as the trivial case. We apply importance sampling in the context of TU games in Section 4.

4. Approximating Linear Solution Concepts via Importance Sampling

The exact calculation of a linear solution concept, i.e., (1), for a TU game normally requires summing up terms whose number grows exponentially in the number of players n. Therefore, approximation algorithms are needed to estimate these values in real-world situations, especially in the context of applications in interpretable machine learning [1,17,18,19] where one can typically not exploit any special structure of the underlying game for computing Shapley values exactly. In this section, we briefly summarize some ideas from [14] and the general framework for importance sampling on the coalition space for linear solution concepts which we recently proposed in [15] and will apply in Section 5.

Benati et al. [14] consider a uniform sampling strategy on the coalition space for approximating linear solution concepts, i.e. they simply adopt the uniform distribution $q=\mathcal{U}\left(2^N\right)$. Let us regard it as the crude Monte Carlo method on the coalition space. Although sampling subsets from the uniform distribution is both straightforward and unbiased, it is obviously not an optimal — or even recommendable — sampling strategy for all linear solution concepts or problem settings. For example, when approximating the Shapley value by sampling coalitions using (), small or large coalitions obtain larger weights resulting in a strong influence on the estimator. However, when sampling uniformly from the coalition space, they may extract only a few samples resulting in a less accurate estimator.

To enable more efficient sampling, we apply the importance sampling technique from SubSection 3.3. By appropriately reweighting the estimator, importance sampling allows us to draw samples from a non-uniform, user-defined distribution while maintaining an unbiased estimate of the underlying linear solution concept.

Theorem 1.

[15] For all $i\in N$, let

$$p_i:2^N\to [0,1]\mathrm{with}{p}_i\left(S\right)=0\Rightarrow w_i\left(S\right){\nu }_i\left(S\right)=0$$

be a probability distribution and ${\mathbb{S}}_i{\sim }_{\mathrm{iid}}{p}_i$ with ${\mathbb{S}}_i={[{\mathcal{S}}_{i,1},\dots ,{\mathcal{S}}_{i,\tau }]}^{\top }$ be a sample of size τ generated by sampling with replacement according to $p_i$. Then,

$$\widehat{\alpha }(N,v)={\left[{\widehat{\alpha }}_i(N,v)\right]}_{i\in N}={\left[{\displaystyle \frac{1}{\tau }}\sum _{\mathcal{S}\in {\mathbb{S}}_i}{\displaystyle \frac{w_i\left(\mathcal{S}\right)}{p_i\left(\mathcal{S}\right)}}{\nu }_i(\mathcal{S},v)\right]}_{i\in N}$$

(16)

is an importance sampling estimator of the linear solution concept $\alpha (N,v)$.

The following proposition connects our importance sampling estimator from Theorem 1 established in Pollmann and Staudacher [15] to a finding from Benati et al. [14], p. 95.

Proposition 1.

The importance sampling estimator from Theorem 1 has the following properties:

(a)

It is unbiased, i.e.,

$$\mathbb{E}\left[{\widehat{\alpha }}_i(N,v)\right]={\alpha }_i(N,v)(\forall i\in N).$$

(17)

(b)

Its variance is given by

$$Var\left[{\widehat{\alpha }}_i(N,v)\right]={\displaystyle \frac{1}{\tau }}\left(\sum _{S\subseteq N}{\displaystyle \frac{w_i{\left(S\right)}^2}{p_i\left(S\right)}}{\nu }_i{(S,v)}^2-{\alpha }_i{(N,v)}^2\right)(\forall i\in N).$$

(18)

(c)

It is consistent in probability, i.e.,

$${\widehat{\alpha }}_i(N,v)\stackrel{\tau \to \infty }{\to }{\alpha }_i(N,v)(\forall i\in N).$$

(19)

We note that Theorem 1 and Proposition 1 subsume both the crude Monte Carlo method (setting $p_i\left(S\right)=q\left(S\right)=2^{-n}$) and stratified sampling. The latter is covered because any stratum estimator for a linear solution concept can be written in the form of (16), allowing direct application of these results.

5. Least Squares Approaches for Approximating Shapley Values

This section covers the idea of approaching the approximation of Shapley values as a least squares problem. First, we give an introduction to the family of least square values in Section 5.1, followed by an overview about Shapley value approximation algorithms in this setting in Section 5.2. In Section 5.3, Section 5.4, and Section 5.5, we compare these algorithms and put them into the overall context.

5.1. The Family of least Square Values

The family of least square values was proposed by Ruiz et al. [9] and comprises a generalization of the formulation of the Shapley value as a least squares problem, which was already proposed earlier by Charnes et al. [8].

Essentially, this family consists of all solution concepts that can be obtained as minimizers of a weighted least squares problem, where the weights determine the relative importance of different coalition sizes. These weights are given by the function $m:\{0,\dots ,n\}\to {\mathbb{R}}_{\ge 0}$ that assigns a weight to each coalition size $s\in \{0,\dots ,n\}$ with $m\left(s\right)>0$ for at least one s. Clearly, per definition, m is symmetric in a sense of assigning the same weight to coalitions of the same cardinality in order to retrieve symmetric solutions like the Shapley value. We denote by $\mathbf{y}\in {\mathbb{R}}^n$ a payoff vector for a given game with $y_i$ being the payoff for player i. In the following, we use the shorthand notation $\mathbf{y}\left(S\right)={\sum }_{i\in S}{y}_i$. Recall that a payoff vector is called efficient if $\mathbf{y}\left(N\right)=v\left(N\right)$.

With these definitions established, any solution concept belonging to the family of least square values can be expressed as the optimizer of the following weighted least squares problem:

$$\begin{array}{cc}\hfill \underset{\mathbf{y}\in {\mathbb{R}}^n}{min}& \sum _{S\subseteq N}m\left(\right|S\left|\right){\left(v\left(S\right)-\mathbf{y}\left(S\right)\right)}^2\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{y}\left(N\right)=v\left(N\right),\hfill \end{array}$$

(20)

with m being the weight function that is specific to the particular solution concept.

The closed-form solution — again, for a fixed m — to problem (20) is

$$y_i={\displaystyle \frac{v\left(N\right)}{n}}+{\displaystyle \frac{1}{\gamma n}}\left(n{a}_i-\sum _{j\in N}{a}_j\right)(\forall i\in N),$$

(21)

with

$$a_j=\sum _{\begin{array}{c}S\subseteq N\\ j\in S\end{array}}m\left(\right|S\left|\right)v\left(S\right)\mathrm{and}\gamma =\sum _{s=1}^{n-1}m\left(s\right)\left({\displaystyle \genfrac{}{}{0pt}{}{n-2}{s-1}}\right).$$

(22)

These expressions illustrate how the resulting allocation depends on the weight function m.

Least square values as linear solution concepts

Let us briefly point out why all members of the family of least square values are linear solution concepts in the sense of the definition from SubSection 2.2. Benati et al. [14] showed that (21) can be rewritten as

$$\mathbf{y}=\mathbf{1}{\displaystyle \frac{v\left(N\right)}{n}}+\sum _{\begin{array}{c}S\subseteq N\\ 0<\left|S\right|<n\end{array}}\mathbf{w}\left(S\right)v\left(S\right),$$

(23)

where $\mathbf{1}\in {\mathbb{R}}^n$ is the all-ones vector and the individual elements of $\mathbf{w}\left(S\right)$ are given by

$$w_i\left(S\right)=\left\{\begin{array}{cc}{\displaystyle \frac{(n-|S\left|\right)m\left(\right|S\left|\right)}{n\gamma }}\hfill & \mathrm{if}i\in S\hfill \\ -{\displaystyle \frac{\left|S\right|m\left(\right|S\left|\right)}{n\gamma }}\hfill & \mathrm{if}i\notin S.\hfill \end{array}\right.$$

(24)

From (23), it follows directly that this representation satisfies the definition of linear solution concepts provided in SubSection 2.2.

The Shapley value in the context of least square values

The family of least square values was introduced by Ruiz et al. [9] after Charnes et al. [8] had demonstrated that the Shapley value can be expressed as a least squares optimization problem. In detail, the vector of Shapley values $\varphi$ is the solution to problem (20) when m is defined as

$$m\left(s\right)={\displaystyle \frac{1}{n-1}}{\left({\displaystyle \genfrac{}{}{0pt}{}{n-2}{s-1}}\right)}^{-1},$$

(25)

and therefore, via (22), one obtains

$$\gamma =\sum _{s=1}^{n-1}m\left(s\right)\left({\displaystyle \genfrac{}{}{0pt}{}{n-2}{s-1}}\right)=\sum _{s=1}^{n-1}{\displaystyle \frac{1}{n-1}}{\left({\displaystyle \genfrac{}{}{0pt}{}{n-2}{s-1}}\right)}^{-1}\left({\displaystyle \genfrac{}{}{0pt}{}{n-2}{s-1}}\right)=1.$$

(26)

Thus, by inserting (26) into (21) and reformulating the result, one obtains the closed-form solution for the Shapley value as

$${\varphi }_i={\displaystyle \frac{v\left(N\right)}{n}}+{\displaystyle \frac{n-1}{n}}{a}_i-{\displaystyle \frac{1}{n}}\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}{a}_j(\forall i\in N),$$

(27)

with $a_j$ given by (22) for all $j\in N$. We illustrate the computation of the Shapley value as a least squares optimization problem in Figure 1.

5.2. Shapley Value Approximation Algorithms

In this subsection, we provide an introduction to Shapley value approximation algorithms that are based on the formulation of the Shapley value as a least squares problem. The algorithms discussed below are not meant to be an exhaustive list, but rather represent those algorithms we will analyse and test in the rest of the article. Specifically, the first three algorithms have been proposed by Benati et al. [14] in the context of stochastic approximations of cooperative games and are based on the family of least square values defined in SubSection 5.1, while the latter two have their origins in the machine learning world.

Additionally, we would like to draw attention to the absence of the projection method proposed by Benati et al. [14], which ensures efficiency and symmetry of the obtained estimators. Although this approach potentially improves the accuracy of the estimators (or at least never makes their accuracy worse), it is not agnostic to the underlying characteristic function, as it requires prior knowledge about symmetric players, which is an information that we assume to be unavailable — in particular, in the context of interpretable machine learning. Nevertheless, one could still achieve estimator efficiency through a single evaluation of $v\left(N\right)$, similar to Castro et al. [28], who proposed filling the “efficiency gap” post hoc. However, since such additional considerations lie outside the scope of this work, we exclude them equally from all algorithms.

Least Squares Sampling (LSS)

This paragraph deals with an algorithm proposed by Benati et al. [14] that approximates (23) by using importance sampling as defined in Theorem 1. Specifically, we refer to the estimator in Equation (22) on p. 97 in the paper [14] combined with their weighted sampling strategy. In the following, we name this algorithm LSS.

The idea of LSS is to sample coalitions with replacement according to

$$p\left(S\right)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{(n-1)\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\hfill & \mathrm{if}0<\left|S\right|<n\hfill \\ 0\hfill & \mathrm{else},\hfill \end{array}\right.$$

(28)

such that, by combining (24), (25), (26), and (28), one obtains

$${\displaystyle \frac{w_i\left(S\right)}{p\left(S\right)}}=\left\{\begin{array}{cc}{\displaystyle \frac{n-1}{\left|S\right|}}\hfill & \mathrm{if}i\in S\hfill \\ -{\displaystyle \frac{n-1}{n-\left|S\right|}}\hfill & \mathrm{if}i\notin S.\hfill \end{array}\right.$$

(29)

Then, based on Theorem 1 as well as (23) and (29), an approximation of the Shapley value of player i is given by

$$\begin{array}{cc}\hfill {\widehat{\varphi }}_i& ={\displaystyle \frac{v\left(N\right)}{n}}+{\displaystyle \frac{1}{\tau }}\sum _{\mathcal{S}\in \mathbb{S}}{\displaystyle \frac{w_i\left(\mathcal{S}\right)}{p\left(\mathcal{S}\right)}}v\left(\mathcal{S}\right)\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{v\left(N\right)}{n}}+{\displaystyle \frac{1}{\tau }}\left(\sum _{\begin{array}{c}\mathcal{S}\in \mathbb{S}\\ i\in \mathcal{S}\end{array}}{\displaystyle \frac{n-1}{\left|\mathcal{S}\right|}}v\left(\mathcal{S}\right)-\sum _{\begin{array}{c}\mathcal{S}\in \mathbb{S}\\ i\notin \mathcal{S}\end{array}}{\displaystyle \frac{n-1}{n-\left|\mathcal{S}\right|}}v\left(\mathcal{S}\right)\right),\hfill \end{array}$$

(31)

where $\mathbb{S}$ is obtained by sampling $\tau$ times with replacement according to p.

Equation (31) serves as the basis for Algorithm 1. Note that one evaluation of v can be used to update all players, reducing the number of evaluations of v and increasing convergence speed.

Algorithm 1 LSS

1: $\widehat{\varphi }\leftarrow \mathbf{1}\frac{v\left(N\right)}{n}$ 2: for$\tau$times do 3:     Draw $\mathcal{S}\sim p$ 4:     $v_{\mathcal{S}}\leftarrow v\left(\mathcal{S}\right)$ 5:     $w_{\mathrm{in}}\leftarrow {\displaystyle \frac{n-1}{\left|\mathcal{S}\right|}}$ 6:     $w_{\mathrm{out}}\leftarrow {\displaystyle \frac{n-1}{n-\left|\mathcal{S}\right|}}$ 7:     for all $i\in N$ do 8:         if $i\in \mathcal{S}$ then 9:            ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i+{\displaystyle \frac{1}{\tau }}{w}_{\mathrm{in}}{v}_{\mathcal{S}}$ 10:         else 11:            ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i-{\displaystyle \frac{1}{\tau }}{w}_{\mathrm{out}}{v}_{\mathcal{S}}$ 12:         end if 13:     end for 14: end for

Since (31) is an importance sampling estimator of a linear solution concept in the sense of Theorem 1, which can easily be obtained via (23) and (30), LSS inherits the properties stated in Proposition 1. Thus, the estimator is unbiased, and we can use (18) to calculate its variance:

Proposition 2.

The variance of the LSS estimator is given by

$$Var\left[{\widehat{\varphi }}_i\right]={\displaystyle \frac{1}{\tau }}\left(\sum _{\begin{array}{c}S\subseteq N\\ 0<\left|S\right|<n\\ i\in S\end{array}}{\displaystyle \frac{n-\left|S\right|}{n\left|S\right|\left(\genfrac{}{}{0pt}{}{n-2}{\left|S\right|-1}\right)}}v{\left(S\right)}^2+\sum _{\begin{array}{c}S\subseteq N\\ 0<\left|S\right|<n\\ i\notin S\end{array}}{\displaystyle \frac{\left|S\right|}{n(n-|S\left|\right)\left(\genfrac{}{}{0pt}{}{n-2}{\left|S\right|-1}\right)}}v{\left(S\right)}^2-{\alpha }_i^2\right)$$

(32)

with

$${\alpha }_i={\varphi }_i-{\displaystyle \frac{v\left(N\right)}{n}}$$

for all $i\in N$.

Proof. 

Using (18) and setting ${\nu }_i=v$ results in

$$Var\left[{\widehat{\varphi }}_i\right]={\displaystyle \frac{1}{\tau }}\left(\sum _{\begin{array}{c}S\subseteq N\\ 0<\left|S\right|<n\end{array}}{\displaystyle \frac{w_i{\left(S\right)}^2}{p\left(S\right)}}v{\left(S\right)}^2-{\alpha }_i^2\right),$$

(33)

where ${\alpha }_i={\varphi }_i-{\displaystyle \frac{v\left(N\right)}{n}}$. Here, the term ${\displaystyle \frac{v\left(N\right)}{n}}$ is excluded from the variance calculation since it is added once at the beginning of the LSS algorithm, but it is not part of the sampling procedure and thus, there is no randomness to it, and it does not influence the variance of ${\widehat{\varphi }}_i$. As a result, we construct a temporary solution concept of the game, ${\alpha }_i$, that does not contain the constant term ${\displaystyle \frac{v\left(N\right)}{n}}$.

With that given, one obtains (32) by calculating

$${\displaystyle \frac{w_i{\left(S\right)}^2}{p\left(S\right)}}=\left\{\begin{array}{cc}{\displaystyle \frac{n-\left|S\right|}{n\left|S\right|}}{\left({\displaystyle \genfrac{}{}{0pt}{}{n-2}{\left|S\right|-1}}\right)}^{-1}\hfill & \mathrm{if}i\in S\hfill \\ {\displaystyle \frac{\left|S\right|}{n(n-|S\left|\right)}}{\left({\displaystyle \genfrac{}{}{0pt}{}{n-2}{\left|S\right|-1}}\right)}^{-1}\hfill & \mathrm{if}i\notin S,\hfill \end{array}\right.$$

and inserting the result into (33).    □

Stratified Least Squares Sampling (S-LSS)

This paragraph discusses another algorithm proposed by Benati et al. [14], which approximates the Shapley value through what the authors introduce as their stratification strategy. Specifically, we refer to equations (23) and (24) on p. 97 in the paper [14] without any reuse of samples across strata. In the following, we refer to this method as S-LSS.

The stratification according to [14] comes into play when approximating the individual addends $a_j$ in (21), where a reformulation of (22) results in

$$a_j=\sum _{s=1}^{n-1}\underset{a_{j,s}}{\underbrace{{\sum }_{\begin{array}{c}S\subseteq N\\ s=\left|S\right|\\ j\in S\end{array}}m\left(s\right)v\left(S\right)}}.$$

(34)

Based on that formulation, approximations of $a_{j,s}$, denoted by ${\widehat{a}}_{j,s}$, are obtained as

$${\widehat{a}}_{j,s}={\displaystyle \frac{1}{{\tau }_{j,s}}}\sum _{\mathcal{S}\in {\mathbb{S}}_{j,s}}{\displaystyle \frac{m\left(s\right)}{p_{j,s}\left(\mathcal{S}\right)}}v\left(\mathcal{S}\right)={\displaystyle \frac{1}{{\tau }_{j,s}(n-s)}}\sum _{\mathcal{S}\in {\mathbb{S}}_{j,s}}v\left(\mathcal{S}\right),$$

(35)

with ${\mathbb{S}}_{j,s}$ being an i.i.d. sample of size ${\tau }_{j,s}$ from $\{S\subseteq N\mid (s=\left|S\right|)\wedge (j\in S\left)\right\}$ according to the uniform probability distribution $p_{j,s}\left(S\right)={\left({\displaystyle \genfrac{}{}{0pt}{}{n-1}{s-1}}\right)}^{-1}$.

Then, based on (27), one obtains an estimator of ${\varphi }_i$ as

$${\widehat{\varphi }}_i={\displaystyle \frac{v\left(N\right)}{n}}+{\displaystyle \frac{n-1}{n}}\sum _{s=1}^{n-1}{\widehat{a}}_{i,s}-{\displaystyle \frac{1}{n}}\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}\sum _{s=1}^{n-1}{\widehat{a}}_{j,s},$$

forming the basis for Algorithm 2.

Algorithm 2 S-LSS

1: for$i\in N$do 2:     for $s=1$to$n-1$ do 3:         ${\widehat{a}}_{i,s}\leftarrow 0$ 4:         for ${\tau }_{i,s}$ times do 5:            Draw $\mathcal{S}\sim p_{i,s}$ 6:            ${\widehat{a}}_{i,s}\leftarrow {\widehat{a}}_{i,s}+{\displaystyle \frac{1}{{\tau }_{i,s}(n-s)}}v\left(\mathcal{S}\right)$ 7:         end for 8:     end for 9: end for 10: ${\widehat{\varphi }}_i={\displaystyle \frac{v\left(N\right)}{n}}+{\displaystyle \frac{n-1}{n}}{\sum }_{s=1}^{n-1}{\widehat{a}}_{i,s}-{\displaystyle \frac{1}{n}}{\sum }_{j\in N:j\ne i}{\sum }_{s=1}^{n-1}{\widehat{a}}_{j,s}(\forall i\in N)$

Clearly, by using the formulation given in (35), the estimators ${\widehat{a}}_{j,s}$ can be seen as estimators of individual solution concepts $a_{j,s}$ in the sense of Theorem 1 by setting

$$w_j\left(S\right)=\left\{\begin{array}{cc}m\left(s\right)\hfill & \mathrm{if}(s=|S\left|\right)\wedge (j\in S)\hfill \\ 0\hfill & \mathrm{else},\hfill \end{array}\right.$$

as well as

$$p_j\left(S\right)=\left\{\begin{array}{cc}{p}_{j,s}\left(S\right)\hfill & \mathrm{if}(s=|S\left|\right)\wedge (j\in S)\hfill \\ 0\hfill & \mathrm{else},\hfill \end{array}\right.$$

and ${\nu }_j=v$.

Therefore, the individual stratum estimators are unbiased, and (18) can be used to calculate their variances:

Proposition 3.

In the context of the S-LSS algorithm, the variance of ${\widehat{a}}_{j,s}$ is given by

$$Var\left[{\widehat{a}}_{j,s}\right]={\displaystyle \frac{1}{{\tau }_{j,s}}}\left({\displaystyle \frac{1}{(n-s)(n-1)\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}}\sum _{\begin{array}{c}S\subseteq N\\ s=\left|S\right|\\ j\in S\end{array}}v{\left(S\right)}^2-a_{j,s}^2\right)$$

(36)

for all $j\in N$ and $s\in \{1,\dots ,n-1\}$.

Proof. 

By using the fact that

$${\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{n-1}{s-1}\right)}{\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}}={\displaystyle \frac{(n-1)!}{(s-1)!(n-s)!}}{\displaystyle \frac{(s-1)!(n-s-1)!}{(n-2)!}}={\displaystyle \frac{n-1}{n-s}},$$

(37)

we can use (18) to obtain

$$\begin{array}{cc}\hfill \sum _{\begin{array}{c}S\subseteq N\\ s=\left|S\right|\\ j\in S\end{array}}{\displaystyle \frac{m{\left(s\right)}^2}{p_{j,s}\left(S\right)}}v{\left(S\right)}^2& ={\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{n-1}{s-1}\right)}{{(n-1)}^2{\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}^2}}\sum _{\begin{array}{c}S\subseteq N\\ s=\left|S\right|\\ j\in S\end{array}}v{\left(S\right)}^2\hfill \\ \hfill & ={\displaystyle \frac{n-1}{(n-s){(n-1)}^2\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}}\sum _{\begin{array}{c}S\subseteq N\\ s=\left|S\right|\\ j\in S\end{array}}v{\left(S\right)}^2\hfill \\ \hfill & ={\displaystyle \frac{1}{(n-s)(n-1)\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}}\sum _{\begin{array}{c}S\subseteq N\\ s=\left|S\right|\\ j\in S\end{array}}v{\left(S\right)}^2.\hfill \end{array}$$

Inserting this result back into (18) results in (36).    □

Based on the individual stratum estimator variances given by Proposition 3, Benati et al. [14] stated the overall variance of the estimated Shapley value of player i as

$$\begin{array}{cc}\hfill Var\left[{\widehat{\varphi }}_i\right]& =Var\left[{\displaystyle \frac{v\left(N\right)}{n}}+{\displaystyle \frac{n-1}{n}}\sum _{s=1}^{n-1}{\widehat{a}}_{i,s}-{\displaystyle \frac{1}{n}}\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}\sum _{s=1}^{n-1}{\widehat{a}}_{j,s}\right]\hfill \\ \hfill & ={\left({\displaystyle \frac{n-1}{n}}\right)}^2\sum _{s=1}^{n-1}Var\left[{\widehat{a}}_{i,s}\right]+{\displaystyle \frac{1}{n^2}}\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}\sum _{s=1}^{n-1}Var\left[{\widehat{a}}_{j,s}\right].\hfill \end{array}$$

(38)

Perhaps surprisingly and contrary to what might be expected from the definition of stratified sampling in SubSection 3.2, S-LSS does not guarantee a smaller variance compared to LSS, as we discuss in SubSection 5.3 and show empirically in Section 6.

Sample Reuse Stratified Least Squares Sampling (SRS-LSS)

Benati et al. [14] proposed a third algorithm, which is based on the S-LSS algorithm and reuses samples across strata. Specifically, we refer to the final three sentences in subsection 4.2 on p. 97 in the paper [14]. In the following, we refer to this method as SRS-LSS.

The central idea is to reuse one evaluation $v\left(\mathcal{S}\right)$ in (35) to update all stratum estimators ${\widehat{a}}_{j,s}$ where $s=\left|\mathcal{S}\right|$ and $j\in \mathcal{S}$, hence the name. To achieve that, samples ${\mathbb{S}}_s$ of size ${\tau }_s$ are drawn from $\{S\subseteq N\mid s=|S\left|\right\}$ with probabilities $p_s\left(S\right)={\left({\displaystyle \genfrac{}{}{0pt}{}{n}{s}}\right)}^{-1}$. Based on the elements of $\mathcal{S}\in {\mathbb{S}}_s$, all stratum estimators ${\widehat{a}}_{j,s}$ with $j\in \mathcal{S}$ are updated with a single evaluation $v\left(\mathcal{S}\right)$. As a side effect, the number of evaluations of v that are used to update player j’s stratum, i.e., ${\tau }_{j,s}$, are dependent on the generated samples, which is why, in Algorithm 3, we need to keep track of their values as well. This essentially means that all ${\tau }_{j,s}$ become estimators themselves, which is why we denote them by ${\widehat{\tau }}_{j,s}$. Additionally, we use ${\tilde{a}}_{j,s}$ to denote unnormalized estimators of $a_{j,s}$ before dividing by ${\widehat{\tau }}_{j,s}(n-s)$.

Algorithm 3 SRS-LSS

1: for$s=1$to$n-1$do 2:     ${\widehat{\tau }}_{i,s}\leftarrow 0(\forall i\in N)$ 3:     ${\tilde{a}}_{i,s}\leftarrow 0(\forall i\in N)$ 4:     for ${\tau }_s$ times do 5:         Draw $\mathcal{S}\sim p_s$ 6:         $v_{\mathcal{S}}\leftarrow v\left(\mathcal{S}\right)$ 7:         for all $i\in \mathcal{S}$ do 8:            ${\tilde{a}}_{i,s}\leftarrow {\tilde{a}}_{i,s}+v_{\mathcal{S}}$ 9:            ${\widehat{\tau }}_{i,s}\leftarrow {\widehat{\tau }}_{i,s}+1$ 10:         end for 11:     end for 12:     ${\widehat{a}}_{i,s}\leftarrow {\displaystyle \frac{1}{{\widehat{\tau }}_{i,s}(n-s)}}{\tilde{a}}_{i,s}(\forall i\in N)$ 13: end for 14: ${\widehat{\varphi }}_i={\displaystyle \frac{v\left(N\right)}{n}}+{\displaystyle \frac{n-1}{n}}{\sum }_{s=1}^{n-1}{\widehat{a}}_{i,s}-{\displaystyle \frac{1}{n}}{\sum }_{j\in N:j\ne i}{\sum }_{s=1}^{n-1}{\widehat{a}}_{j,s}(\forall i\in N)$

One needs to guarantee that ${\widehat{\tau }}_{j,s}>0$ for all $j\in N$ and $s\in \{1,\dots ,n-1\}$ in order for SRS-LSS to function properly and avoid divisions by zero. One straightforward approach is to repeatedly generate samples until this condition is satisfied, or alternatively, to employ a conditional probability distribution that accounts for which elements have not yet been sampled. In the remainder of this work, we rely on the results from the following proposition, and therefore, we do not incorporate additional safeguards into Algorithm 3.

Proposition 4.

Let ${\tau }_s\to \infty$ as $\tau \to \infty$ for all $s\in \{1,\dots ,n-1\}$. Then, with probability approaching one asymptotically, each stratum estimator receives at least one sample, i.e.,

$$\underset{\tau \to \infty }{lim}\mathbb{P}(\exists j\in N,s\in \{1,\dots ,n-1\}:{\widehat{\tau }}_{j,s}=0)=0.$$

Proof. 

The probability of a player $j\in N$ belonging to a random coalition of size $s\in \{1,\dots ,n-1\}$ is given by ${\displaystyle \frac{s}{n}}>0$. Thus, we have

$$\mathbb{P}({\widehat{\tau }}_{j,s}=0)={\left(1-{\displaystyle \frac{s}{n}}\right)}^{{\tau }_s}\stackrel{{\tau }_s\to \infty }{\to }0.$$

(39)

Since

$$\mathbb{P}(\exists j\in N,s\in \{1,\dots ,n-1\}:{\widehat{\tau }}_{j,s}=0)=\mathbb{P}\left(\bigcup _{j\in N}\bigcup _{s=1}^{n-1}\{{\widehat{\tau }}_{j,s}=0\}\right),$$

we can use the union bound [29] and (39) to derive

$$\mathbb{P}\left(\bigcup _{j\in N}\bigcup _{s=1}^{n-1}\{{\widehat{\tau }}_{j,s}=0\}\right)\le \sum _{j\in N}\sum _{s=1}^{n-1}\mathbb{P}({\widehat{\tau }}_{j,s}=0)\stackrel{\tau \to \infty }{\to }0,$$

which holds as long as ${\tau }_s\to \infty$ with $\tau \to \infty$. □

Since Benati et al. [14] did not include a variance analysis of the SRS-LSS estimator in their work, we refer to SubSection 5.4 for a detailed discussion of the variance of the obtained estimator.

KernelSHAP (KS)

Let us now introduce the KernelSHAP algorithm, one of the predominant Shapley approximation algorithms in the machine learning world. It was developed by Lundberg and Lee [1], without knowledge of prior works like [8] and [9]. The optimization objective is

$$\begin{array}{cc}\hfill \varphi =\underset{\mathbf{y}\in {\mathbb{R}}^n}{arg\; min}& \sum _{\begin{array}{c}\mathbf{z}\in {\{0,1\}}^n\\ 0<\left|\mathbf{z}\right|<n\end{array}}p\left(\mathbf{z}\right){\left(\mathbf{y}\left(\mathbf{z}\right)-v\left(\mathbf{z}\right)\right)}^2\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{y}\left(N\right)=v\left(N\right)\hfill \end{array}$$

(40)

to calculate the Shapley values for the cooperative game $(N,v)$ [7]. Here, $p\propto \omega$ denotes a probability mass function satisfying $p\left(\mathbf{0}\right)=0$ and $p\left(\mathbf{1}\right)=0$ such that

$$p\left(\mathbf{z}\right)=\left\{\begin{array}{cc}{Q}^{-1}\omega \left(\mathbf{z}\right)\hfill & \mathrm{if}0<\left|\mathbf{z}\right|<n\hfill \\ 0\hfill & \mathrm{else},\hfill \end{array}\right.$$

(41)

where

$$\omega \left(\mathbf{z}\right)={\displaystyle \frac{n-1}{\left(\genfrac{}{}{0pt}{}{n}{\left|\mathbf{z}\right|}\right)\left|\mathbf{z}\right|(n-\left|\mathbf{z}\right|)}}$$

(42)

is the Shapley kernel and

$$Q=(n-1)\sum _{s=1}^{n-1}{\displaystyle \frac{1}{s(n-s)}}$$

(43)

is a normalization constant [7]. As noted in SubSection 2.1, all functions defined on ${\{0,1\}}^n$ are implicitly also defined on N.

Since solving (40) requires optimizing over the elements in ${\{0,1\}}^n$, a set that grows exponentially in n, KernelSHAP approximates Shapley values by generating a sample $\mathbb{S}{\sim }_{\mathrm{iid}}p$ of size $\tau$ and solving the approximate optimization problem defined over the sampled subsets only, i.e.,

$$\begin{array}{cc}\hfill \widehat{\varphi }=\underset{\mathbf{y}\in {\mathbb{R}}^n}{arg\; min}& {\displaystyle \frac{1}{\tau }}\sum _{\mathcal{S}\in \mathbb{S}}{\left(\mathbf{y}\left(\mathcal{S}\right)-v\left(\mathcal{S}\right)\right)}^2\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \mathbf{y}\left(N\right)=v\left(N\right).\hfill \end{array}$$

(44)

Following [7], the Lagrangian of the approximate problem given by (44) with multiplier $\lambda \in \mathbb{R}$ is

$$\begin{array}{cc}\hfill \widehat{\mathcal{L}}(\mathbf{y},\lambda )& ={\mathbf{y}}^{\top }\left(\underset{\widehat{A}}{\underbrace{{\displaystyle \frac{1}{\tau }}{\sum }_{\mathcal{S}\in \mathbb{S}}\mathbf{z}\left(\mathcal{S}\right)\mathbf{z}{\left(\mathcal{S}\right)}^{\top }}}\right)\mathbf{y}-2{\mathbf{y}}^{\top }\left(\underset{\widehat{\mathbf{b}}}{\underbrace{{\displaystyle \frac{1}{\tau }}{\sum }_{\mathcal{S}\in \mathbb{S}}\mathbf{z}\left(\mathcal{S}\right)v\left(\mathcal{S}\right)}}\right)\hfill \\ \hfill & +{\displaystyle \frac{1}{\tau }}\sum _{\mathcal{S}\in \mathbb{S}}v{\left(\mathcal{S}\right)}^2+2\lambda \left(\mathbf{y}\left(N\right)-v\left(N\right)\right).\hfill \end{array}$$

(45)

Based on the Lagrangian, Covert and Lee [7] derived the closed-form solution

$$\widehat{\varphi }={\widehat{A}}^{-1}\left(\widehat{\mathbf{b}}-\mathbf{1}{\displaystyle \frac{{\mathbf{1}}^{\top }{\widehat{A}}^{-1}\widehat{\mathbf{b}}-v\left(N\right)}{{\mathbf{1}}^{\top }{\widehat{A}}^{-1}\mathbf{1}}}\right).$$

The resulting estimator is consistent in probability, but assessing its unbiasedness or variance is challenging due to terms like ${\widehat{A}}^{-1}\widehat{\mathbf{b}}$. While the unbiasedness or variance of $\widehat{A}$ or $\widehat{\mathbf{b}}$ can be analyzed individually [7], their interaction complicates the analysis of the final Shapley value estimator.

The KernelSHAP algorithm will not be further analyzed theoretically in this work. Due to its setup, it is unclear how the obtained estimator fits the framework defined by Theorem 1 or Monte Carlo estimators introduced in Section 3 in general. To the best of our knowledge, no analytic solution for its variance or a proof of its unbiasedness has been derived up today, which is why the unbiased KernelSHAP algorithm was introduced, guaranteeing unbiasedness and facilitating the variance analysis of the obtained Shapley value estimator.

Unbiased KernelSHAP (UKS)

The UKS algorithm [7] is a continued development of the KernelSHAP algorithm from the previous paragraph, where we pointed out that it is unclear if the KernelSHAP estimator is unbiased due to the interactions between $\widehat{A}$ and $\widehat{\mathbf{b}}$.

Considering the original optimization problem defined in (40), its Lagrangian with multiplier $\lambda \in \mathbb{R}$ is given by

$$\mathcal{L}(\mathbf{y},\lambda )={\mathbf{y}}^{\top }\underset{A}{\underbrace{\mathbb{E}\left[Z{Z}^{\top }\right]}}\mathbf{y}-2{\mathbf{y}}^{\top }\underset{\mathbf{b}}{\underbrace{\mathbb{E}\left[Zv\right(Z\left)\right]}}+\mathbb{E}\left[v{\left(Z\right)}^2\right]+2\lambda \left({\mathbf{1}}^{\top }\mathbf{y}-v\left(N\right)\right),$$

whereby Z is a random variable taking values from ${\{0,1\}}^n$ and being distributed according to p defined by (41). With that formulation given, the closed-form solution of the original problem can be derived as

$$\varphi =A^{-1}\left(\mathbf{b}-\mathbf{1}{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{b}-v\left(N\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right).$$

Clearly, $\mathbf{b}$ cannot be computed exactly for large n, since it depends on v, and thus, needs $2^n-2$ evaluations of v in order to be exact. Nevertheless, since A does not depend on v and the probability distribution p is known, A can be pre-computed and is exact, which is the main improvement introduced by UKS in contrast to the original KernelSHAP algorithm.

Thus, UKS estimates all players’ Shapley values via

$$\widehat{\varphi }=A^{-1}\left(\widehat{\mathbf{b}}-\mathbf{1}{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\widehat{\mathbf{b}}-v\left(N\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right),$$

(46)

where $\widehat{\mathbf{b}}$ is defined in (45) as for the KernelSHAP algorithm, i.e., $\widehat{\mathbf{b}}={\displaystyle \frac{1}{\tau }}{\sum }_{\mathcal{S}}\mathbf{z}\left(\mathcal{S}\right)v\left(\mathcal{S}\right)$, and $A\in {\mathbb{R}}^{n\times n}$ is given by

$$A=\left[\begin{array}{cccc}{\displaystyle \frac{1}{2}}& \rho & \dots & \rho \\ \rho & {\displaystyle \frac{1}{2}}& \dots & \rho \\ ⋮& ⋮& \ddots & ⋮\\ \rho & \rho & \dots & {\displaystyle \frac{1}{2}}\end{array}\right]$$

(47)

with

$$\rho ={\displaystyle \frac{1}{n(n-1)}}{\displaystyle \frac{{\sum }_{s=2}^{n-1}\frac{s-1}{n-s}}{{\sum }_{s=1}^{n-1}\frac{1}{s(n-s)}}}.$$

(48)

Following the results provided in [7], the resulting Shapley value estimator is unbiased and consistent in probability. Furthermore, its covariance is given by

$${\Sigma }_{\widehat{\varphi }}=E{\Sigma }_{\widehat{\mathbf{b}}}{E}^{\top }$$

(49)

with

$${\Sigma }_{\widehat{\mathbf{b}}}=Cov\left[Zv\left(Z\right)\right]$$

and

$$E=A^{-1}-{\displaystyle \frac{A^{-1}\mathbf{1}{\mathbf{1}}^{\top }{A}^{-1}}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}.$$

5.3. Comparison of LSS and S-LSS

In this subsection, we first clarify the relationship between the two algorithms LSS and S-LSS which we introduced at the beginning of SubSection 5.2, followed by a detailed comparison of variances.

Algorithm Comparison

For comparing the algorithms LSS and S-LSS, we start with the following proposition:

Proposition 5.

In the sense of the definition of stratified sampling provided in Section 3.2, S-LSS is not a valid stratification scheme of LSS.

Proof. 

As stated in Section 3.2, stratification divides the sample space $\mathcal{X}$ into strata $\{{\mathcal{X}}_1,\dots ,{\mathcal{X}}_{\ell }\}$, such that ${\mathcal{X}}_l\cap {\mathcal{X}}_{l^{\prime }}=\varnothing$ for $l,l^{\prime }\in \{1,\dots ,\ell \}$ and $l\ne l^{\prime }$.

By observing the strata definitions given in (34), one obtains that their sample spaces have the intersection

$$\{S\subseteq N\mid (s=\left|S\right|)\wedge (j\in S)\wedge (j^{\prime }\in S)\}(\forall j\ne j^{\prime }\in N).$$

Clearly, for any coalition size $s\in \{2,\dots ,n-1\}$, this set is nonempty, since it contains at least one element of the form $\{j,j^{\prime },\dots \}$, whereby the dots indicate the presence of $s-2$ distinct players other than j and $j^{\prime }$. Therefore, S-LSS does not form a valid partition of the sample space as required by our definition in SubSection 3.2. □

We note in passing that the structure of the overlap between strata is such that it does not introduce a bias in the Shapley estimator defined by S-LSS.

In order to generate a deeper understanding regarding the relationship between LSS and S-LSS, we explain how Benati et al. [14] derived the LSS algorithm from the initial formulation given by (21). Via (21) and (22), one directly obtains

$$\begin{array}{cc}\hfill y_i& ={\displaystyle \frac{v\left(N\right)}{n}}+{\displaystyle \frac{1}{\gamma n}}\left(n\underset{a_i}{\underbrace{{\sum }_{\begin{array}{c}S\subseteq N\\ i\in S\end{array}}m\left(\right|S\left|\right)v\left(S\right)}}-\sum _{j\in N}\underset{a_j}{\underbrace{{\sum }_{\begin{array}{c}S\subseteq N\\ j\in S\end{array}}m\left(\right|S\left|\right)v\left(S\right)}}\right)\hfill \\ \hfill & ={\displaystyle \frac{v\left(N\right)}{n}}+\sum _{S\subseteq N}{\displaystyle \frac{n{\mathbb{1}}_{i\in S}}{\gamma n}}m\left(\right|S\left|\right)v\left(S\right)-\underset{{\sum }_{S\subseteq N}{\displaystyle \frac{\left|S\right|}{\gamma n}}m\left(\right|S\left|\right)v\left(S\right)}{\underbrace{{\sum }_{j\in N}{\sum }_{\begin{array}{c}S\subseteq N\\ j\in S\end{array}}{\displaystyle \frac{1}{\gamma n}}m\left(\right|S\left|\right)v\left(S\right)}}\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{v\left(N\right)}{n}}+\sum _{S\subseteq N}\underset{w_i\left(S\right)}{\underbrace{{\displaystyle \frac{n{\mathbb{1}}_{i\in S}-\left|S\right|}{\gamma n}}m\left(\right|S\left|\right)}}v\left(S\right),\hfill \end{array}$$

(51)

which matches the formulation proposed in (23) with weights given by (24). Therefore, one evaluation of v for a given S in (51) corresponds to updating all $a_j$ where $j\in S$ in (50). This implicit update of all $a_j$ where $j\in S$ is encapsulated by the weight function $w_i$.

As a result, stratifying by s and j, as done in the case of S-LSS, cannot be derived from LSS in the same way as defined in Section 3.2, since LSS is an algorithm that simultaneously updates multiple strata of S-LSS with a single sampled element, i.e., with one evaluation of v.

We conclude that the statement from SubSection 3.2, which points out that stratified sampling always yields a variance less than or equal to that of the crude Monte Carlo method, does not apply here, since S-LSS is not a valid stratification scheme of LSS. Therefore, in order to understand how their variances compare, we calculate both algorithms’ variances on two selected examples.

Variance Comparison in the Absence of Variance Within Strata

First, let us compare the variances of the LSS and S-LSS estimator in the context of an example game characterized by large variance between strata and no variance within strata. In particular, we consider a cardinality game with $n\ge 2$ players and the characteristic function

$$v\left(S\right)=\left|S\right|.$$

(52)

Lemma 1.

On the game defined by (52) with $n\ge 2$, the variance of the LSS estimator is greater zero for all players, i.e.,

$$Var\left[{\widehat{\varphi }}_i^{\mathrm{LSS}}\right]>0(\forall i\in N).$$

Proof. 

First, we obtain that the Shapley values of the game are given by ${\displaystyle \frac{v\left(N\right)}{n}}$ for all $i\in N$, since all i are symmetric. Without loss of generality, we fix one player $i\in N$ in order to simplify notations. Then, by using Proposition 2, we derive

$${\alpha }_i^2={\left({\varphi }_i-{\displaystyle \frac{v\left(N\right)}{n}}\right)}^2=0.$$

Inserting ${\alpha }_i^2$ into (32), multiplying by $\tau$, and using the identities defined by (37) as well as

$${\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{n-1}{s}\right)}{\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}}={\displaystyle \frac{(n-1)!}{s!(n-s-1)!}}{\displaystyle \frac{(s-1)!(n-s-1)!}{(n-2)!}}={\displaystyle \frac{n-1}{s}}$$

results in

$$\begin{array}{cc}\hfill \tau Var\left[{\widehat{\varphi }}_i^{\mathrm{LSS}}\right]& =\sum _{\begin{array}{c}S\subseteq N\\ 0<\left|S\right|<n\\ i\in S\end{array}}{\displaystyle \frac{n-\left|S\right|}{n\left|S\right|\left(\genfrac{}{}{0pt}{}{n-2}{\left|S\right|-1}\right)}}v{\left(S\right)}^2+\sum _{\begin{array}{c}S\subseteq N\\ 0<\left|S\right|<n\\ i\notin S\end{array}}{\displaystyle \frac{\left|S\right|}{n(n-|S\left|\right)\left(\genfrac{}{}{0pt}{}{n-2}{\left|S\right|-1}\right)}}v{\left(S\right)}^2\hfill \\ \hfill & =\sum _{s=1}^{n-1}\left(\sum _{\begin{array}{c}S\subseteq N\\ s=\left|S\right|\\ i\in S\end{array}}{\displaystyle \frac{n-s}{sn\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}}+\sum _{\begin{array}{c}S\subseteq N\\ s=\left|S\right|\\ i\notin S\end{array}}{\displaystyle \frac{s}{n(n-s)\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}}\right)s^2\hfill \\ \hfill & =\sum _{s=1}^{n-1}\left(\left({\displaystyle \genfrac{}{}{0pt}{}{n-1}{s-1}}\right){\displaystyle \frac{n-s}{sn\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}}+\left({\displaystyle \genfrac{}{}{0pt}{}{n-1}{s}}\right){\displaystyle \frac{s}{n(n-s)\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}}\right)s^2\hfill \\ \hfill & =\sum _{s=1}^{n-1}\left(\underset{>0}{\underbrace{{\displaystyle \frac{n-1}{n-s}}\phantom{{\displaystyle \frac{s}{(n-s)n}}}}}\underset{>0}{\underbrace{{\displaystyle \frac{n-s}{sn}}\phantom{{\displaystyle \frac{s}{(n-s)n}}}}}+\underset{>0}{\underbrace{{\displaystyle \frac{n-1}{s}}\phantom{{\displaystyle \frac{s}{(n-s)n}}}}}\underset{>0}{\underbrace{{\displaystyle \frac{s}{n(n-s)}}}}\right)\underset{>0}{\underbrace{s^2\phantom{{\displaystyle \frac{s}{n(n-s)}}}}}>0.\hfill \end{array}$$

Dividing by $\tau$, we obtain

$$Var\left[{\widehat{\varphi }}_i^{\mathrm{LSS}}\right]>0(\forall i\in N).$$

□

Lemma 2.

On the game defined by (52) with $n\ge 2$, the variance of the S-LSS estimator is zero for all players, i.e.,

$$Var\left[{\widehat{\varphi }}_i^{\mathrm{S}-\mathrm{LSS}}\right]=0(\forall i\in N),$$

as long as ${\tau }_{j,s}>0$ for all $j\in N$ and $s\in \{1,\dots ,n-1\}$.

Proof. 

Without loss of generality, we fix $j\in N$ and $s\in \{1,\dots ,n-1\}$ in order to simplify notations.

To calculate the variances of the individual stratum estimators via (36), the exact values of the stratum estimators $a_{j,s}$ are required. By using the identity (37), these are given by

$$\begin{array}{cc}\hfill a_{j,s}& =\sum _{\begin{array}{c}S\subseteq N\\ s=\left|S\right|\\ j\in S\end{array}}m\left(s\right)v\left(S\right)={\displaystyle \frac{1}{(n-1)\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}}\sum _{\begin{array}{c}S\subseteq N\\ s=\left|S\right|\\ j\in S\end{array}}v\left(S\right)={\displaystyle \frac{1}{(n-1)\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}}\left({\displaystyle \genfrac{}{}{0pt}{}{n-1}{s-1}}\right)s\hfill \\ \hfill & ={\displaystyle \frac{n-1}{(n-1)(n-s)}}s={\displaystyle \frac{s}{n-s}}.\hfill \end{array}$$

(53)

Furthermore, by using (37), the sum of squared evaluations of v in (36) is given by

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{(n-s)(n-1)\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}}\sum _{\begin{array}{c}S\subseteq N\\ s=\left|S\right|\\ j\in S\end{array}}v{\left(S\right)}^2& ={\displaystyle \frac{1}{(n-s)(n-1)\left(\genfrac{}{}{0pt}{}{n-2}{s-1}\right)}}\left({\displaystyle \genfrac{}{}{0pt}{}{n-1}{s-1}}\right)s^2\hfill \\ \hfill & ={\displaystyle \frac{n-1}{(n-1){(n-s)}^2}}{s}^2={\displaystyle \frac{s^2}{{(n-s)}^2}}.\hfill \end{array}$$

(54)

Inserting (53) and (54) into (36) results in

$$Var\left[{\widehat{a}}_{j,s}\right]={\displaystyle \frac{1}{{\tau }_{j,s}}}\left({\displaystyle \frac{s^2}{{(n-s)}^2}}-a_{j,s}^2\right)={\displaystyle \frac{1}{{\tau }_{j,s}}}\left({\displaystyle \frac{s^2}{{(n-s)}^2}}-{\displaystyle \frac{s^2}{{(n-s)}^2}}\right)=0,$$

as long as ${\tau }_{j,s}>0$. Then, one can easily verify via (38) that

$$Var\left[{\widehat{\varphi }}_i^{\mathrm{S}-\mathrm{LSS}}\right]=0(\forall i\in N).$$

□

Corollary 1.

There exist games for which the variances of all players’ S-LSS estimators are smaller than those the corresponding LSS estimators, i.e.,

$$Var\left[{\widehat{\varphi }}_i^{\mathrm{S}-\mathrm{LSS}}\right]<Var\left[{\widehat{\varphi }}_i^{\mathrm{LSS}}\right](\forall i\in N),$$

as long as the sample allocation scheme of S-LSS satisfies ${\tau }_{j,s}>0$ for all $j\in N$ and $s\in \{1,\dots ,n-1\}$.

Proof. 

The proof is straightforward and follows directly from Lemmas 1 and 2. □

Variance comparison in the presence of variance within strata

Let us now consider a weighted voting game as defined in (6) with

$$\mathbf{c}={[1,2,3]}^{\top }\mathrm{and}C=5$$

(55)

as a second example. Here, we use concrete numbers in order to simplify calculations. In detail, we are interested in the variances of the Shapley value estimators of player 2 with Shapley value ${\varphi }_2={\displaystyle \frac{1}{2}}$.

To simplify our later reasoning, we first calculate the population and estimator variances for all strata of the S-LSS estimator in Table 1.

Lemma 3.

On the game defined by (55), the variance of player 2’s LSS estimator is given by

$$Var\left[{\widehat{\varphi }}_2^{\mathrm{LSS}}\right]={\displaystyle \frac{5}{36\tau }}.$$

Proof. 

The proof is straightforward by inserting (55) into (32). □

Lemma 4.

On the game defined by (55), the variance of player 2’s S-LSS estimator with equal sample allocation is given by

$$Var\left[{\widehat{\varphi }}_2^{\mathrm{S}-\mathrm{LSS}}\right]={\displaystyle \frac{5}{6\tau }}.$$

Proof. 

The proof is straightforward and is done by calculating (36), first, to obtain values for all $Var\left[{\widehat{a}}_{j,s}\right]$. Note that these $Var\left[{\widehat{a}}_{j,s}\right]$ are given in the “estimator variances” column from Table 1. Then, one obtains the final estimator’s variance by inserting all $Var\left[{\widehat{a}}_{j,s}\right]$ into (38). □

Corollary 2.

There exist games for which the variance of a player’s S-LSS estimator with equal sample allocation is larger than that of the LSS estimator, i.e., $Var\left[{\widehat{\varphi }}_i^{\mathrm{S}-\mathrm{LSS}}\right]>Var\left[{\widehat{\varphi }}_i^{\mathrm{LSS}}\right]$ for some $i\in N$.

Proof. 

The proof is straightforward and directly follows from Lemmas 3 and 4. □

Although an equal sample allocation scheme might be a reasonable choice, especially when calculating all players’ Shapley values, one might still ask whether alternative allocation procedures could further reduce the variance of the S-LSS estimator. While our goal is to treat characteristic functions — and, in the context of explainable machine learning, the machine learning models behind them — as black boxes without relying on internal knowledge, let us examine a sample allocation strategy that makes use of a priori information about the variances within strata. Interestingly, we will see that even this variance-informed allocation performs worse than the LSS method on the weighted voting game defined by (55).

Lemma 5.

On the game defined by (55), the variance of player 2’s S-LSS estimator with an optimum sample allocation, i.e., a sample allocation that is tailored to player 2 in order to minimize the variance of their estimator, is given by

$$Var\left[{\widehat{\varphi }}_2^{\mathrm{S}-\mathrm{LSS}}\right]={\displaystyle \frac{1}{4\tau }}.$$

Proof. 

By inserting the estimator variances from Table 1 into (38), we obtain

$$\begin{array}{cc}\hfill Var\left[{\widehat{\varphi }}_2^{\mathrm{S}-\mathrm{LSS}}\right]& ={\displaystyle \frac{4}{9}}{\displaystyle \frac{1}{4{\tau }_{2,2}}}+{\displaystyle \frac{1}{9}}{\displaystyle \frac{1}{4{\tau }_{3,2}}}={\displaystyle \frac{1}{9}}\left({\displaystyle \frac{1}{{\tau }_{2,2}}}+{\displaystyle \frac{1}{4{\tau }_{3,2}}}\right)\hfill \\ \hfill & ={\displaystyle \frac{1}{9}}\left({\displaystyle \frac{1}{x\tau }}+{\displaystyle \frac{1}{4(1-x)\tau }}\right)={\displaystyle \frac{1}{9\tau }}\left({\displaystyle \frac{1}{x}}+{\displaystyle \frac{1}{4(1-x)}}\right),\hfill \end{array}$$

(56)

whereby $x\in (0,1)$ is the fraction of $\tau$ that is used for estimating the stratum with $j=2,s=2$ and $(1-x)$ is the fraction of $\tau$ that is used for estimating the stratum with $j=3,s=2$.

Therefore, we get the proxy minimization problem

$$\underset{x\in (0,1)}{min}{\displaystyle \frac{1}{x}}+{\displaystyle \frac{1}{4(1-x)}}.$$

Let us now set

$$f\left(x\right)={\displaystyle \frac{1}{x}}+{\displaystyle \frac{1}{4(1-x)}}$$

with its derivative given by

$$f^{\prime }\left(x\right)=-{\displaystyle \frac{1}{x^2}}+{\displaystyle \frac{1}{4{(1-x)}^2}}.$$

Setting $f^{\prime }\left(x\right)=0$, we get

$$x={\displaystyle \frac{4\pm 2}{3}},$$

with $x={\displaystyle \frac{2}{3}}$ being the only solution in $(0,1)$. Furthermore, since

$$f^{\prime \prime }\left(x\right)={\displaystyle \frac{2}{x^3}}+{\displaystyle \frac{1}{2{(1-x)}^3}}>0\mathrm{on}(0,1),$$

we obtain that $x={\displaystyle \frac{2}{3}}$ indeed is a minimum with value $f\left({\displaystyle \frac{2}{3}}\right)={\displaystyle \frac{9}{4}}.$ Thus, via (56), the minimum variance of the S-LSS estimator is given by

$$Var\left[{\widehat{\varphi }}_2^{\mathrm{S}-\mathrm{LSS}}\right]={\displaystyle \frac{1}{9\tau }}{\displaystyle \frac{9}{4}}={\displaystyle \frac{1}{4\tau }}.$$

□

Corollary 3.

There exist games for which the variance of a player’s S-LSS estimator with optimum sample allocation, i.e., a sample allocation that is tailored to that specific player in order to minimize the variance of their estimator, is larger than that of the LSS estimator, i.e., $Var\left[{\widehat{\varphi }}_i^{\mathrm{S}-\mathrm{LSS}}\right]>Var\left[{\widehat{\varphi }}_i^{\mathrm{LSS}}\right]$ for some $i\in N$.

Proof. 

The proof is straightforward and directly follows from Lemmas 3 and 5. □

Concluding Variance Comparison of LSS and S-LSS

Recapitulating the results from the previous two paragraphs, we obtain:

Theorem 2.

Depending on the specific cooperative game, the variance of a player’s Shapley value estimator may be smaller under either LSS or S-LSS, even when the S-LSS algorithm uses an optimal, player-specific sample allocation.

Proof. 

The proof is straightforward and directly follows from Corollaries 1 and 3. □

To summarize the results in this subsection so far, Proposition 5 demonstrates that S-LSS is not the stratified variant of LSS in a sense of the definition provided in Section 3.2, and Theorem 2 states that neither of both algorithms is preferable over the other one on all cooperative games. Since the proof of Theorem 2 builds on the calculation of both algorithms’ theoretical variances on two distinct games, it does not provide an intuitive explanation why S-LSS may perform worse than LSS in some settings. Therefore, we provide an explanation in the following.

Since Benati et al. [14] did not provide a sample allocation scheme for the S-LSS algorithm, we first derive a proportional sample allocation similar to SubSection 3.2, which is motivated by the fact that all coalition sizes are equally likely in the context of LSS, i.e.,

$$\mathbb{P}(s=|\mathcal{S}\left|\right)={\displaystyle \frac{1}{n-1}}(\forall s\in \{1,\dots ,n-1\}),$$

with $\mathcal{S}\sim p^{\mathrm{LSS}}$, whereby $p^{\mathrm{LSS}}$ is the probability mass function from (28).

Furthermore, once a random coalition $\mathcal{S}\sim p^{\mathrm{LSS}}$ of some size $s\in \{1,\dots ,n-1\}$ is sampled, the probability of any player being in it is the same for all players, i.e.,

$$\mathbb{P}(j\in \mathcal{S}\mid s=|\mathcal{S}\left|\right)={\displaystyle \frac{s}{n}}(\forall j\in N).$$

Then, the overall probability that a random coalition $\mathcal{S}\sim p^{\mathrm{LSS}}$ belongs to a specific stratum of the S-LSS algorithm — which is the case when $\mathcal{S}$ meets a specific size s and, additionally, some specific player j is in it — is given by

$$\mathbb{P}(s=|\mathcal{S}\left|\right)\mathbb{P}(j\in \mathcal{S}\mid s=|\mathcal{S}\left|\right)={\displaystyle \frac{s}{n(n-1)}}(\forall j\in N,\forall s\in \{1,\dots ,n-1\}).$$

(57)

As we established at the beginning of this subsection, the statement on variance reduction for a stratified estimator with stratum sample sizes proportional to the probabilities of sampled elements belonging to a stratum from SubSection 3.2 does not apply in this context since S-LSS is not a valid stratified variant of LSS. However, the proportional sample allocation will come out as a useful starting point for our subsequent explanations. Moreover, as established in Theorem 2, the exact choice of the applied sample allocation scheme does not alter the fact that S-LSS can perform worse than LSS, since the increased variance of S-LSS was also observed under an optimal allocation scheme. Additionally, as we will show later, sample allocation strategies other than the proportional sample allocation do not change the outcome and conclusions of our following reasoning.

Hence we set all ${\tau }_{j,s}^{\mathrm{S}-\mathrm{LSS}}\propto \mathbb{P}(s=|\mathcal{S}\left|\right)\mathbb{P}(j\in \mathcal{S}\mid s=|\mathcal{S}\left|\right)$. Then, we obtain

$$\begin{array}{cc}\hfill {\tau }_{j,s}^{\mathrm{S}-\mathrm{LSS}}& =\tau {\displaystyle \frac{\mathbb{P}(s=|\mathcal{S}\left|\right)\mathbb{P}(j\in \mathcal{S}\mid s=|\mathcal{S}\left|\right)}{{\sum }_{s^{\prime }=1}^{n-1}{\sum }_{j^{\prime }\in N}\mathbb{P}(s^{\prime }=\left|\mathcal{S}\right|\left)\mathbb{P}\right(j^{\prime }\in \mathcal{S}\mid s^{\prime }=\left|\mathcal{S}\right|)}}\hfill \\ \hfill & ={\displaystyle \frac{s\tau }{n(n-1){\sum }_{s^{\prime }=1}^{n-1}{\sum }_{j^{\prime }\in N}\frac{s^{\prime }}{n(n-1)}}}={\displaystyle \frac{s\tau }{{\sum }_{s^{\prime }=1}^{n-1}{\sum }_{j^{\prime }\in N}{s}^{\prime }}}={\displaystyle \frac{s\tau }{n{\sum }_{s^{\prime }=1}^{n-1}{s}^{\prime }}}\hfill \\ \hfill & ={\displaystyle \frac{2s\tau }{n^2(n-1)}}(\forall j\in N,\forall s\in \{1,\dots ,n-1\}).\hfill \end{array}$$

(58)

On the other hand, let us now find out how often a sample is expected to belong to a stratum $a_{j,s}$ of S-LSS when actually running LSS, which is denoted by $\mathbb{E}\left[{\tau }_{j,s}^{\mathrm{LSS}}\right]$. Via (57), we obtain

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\tau }_{j,s}^{\mathrm{LSS}}\right]& =\tau \mathbb{P}(s=|\mathcal{S}\left|\right)\mathbb{P}(j\in \mathcal{S}\mid s=|\mathcal{S}\left|\right)\hfill \\ \hfill & ={\displaystyle \frac{s\tau }{n(n-1)}}(\forall j\in N,\forall s\in \{1,\dots ,n-1\}).\hfill \end{array}$$

(59)

Finally, by comparing (58) and (59), we observe

$${\tau }_{j,s}^{\mathrm{S}-\mathrm{LSS}}={\displaystyle \frac{2s\tau }{n^2(n-1)}}\stackrel{n>2}{<}{\displaystyle \frac{s\tau }{n(n-1)}}=\mathbb{E}\left[{\tau }_{j,s}^{\mathrm{LSS}}\right]$$

(60)

for all $j\in N$ and $s\in \{1,\dots ,n-1\}$. Clearly, (60) shows that the expected sample size for stratum $a_{j,s}$ is larger by a factor of ${\displaystyle \frac{n}{2}}$ when executing LSS in comparison to the corresponding sample size of S-LSS with a proportional sample allocation. This is due to the fact that LSS actually updates multiple strata with one evaluation of v, see the beginning of this subsection. We conclude that S-LSS might eliminate the variance between the estimators of different strata (similar to the definition provided in SubSection 3.2), but taking into account the results from Proposition 3, one obtains the dependency on ${\tau }_{j,s}^{\mathrm{S}-\mathrm{LSS}}$ when quantifying the stratum estimators’ variances, which results in larger variances of the individual stratum estimators due to smaller ${\tau }_{j,s}^{\mathrm{S}-\mathrm{LSS}}$.

This underscores our observations from the previous paragraphs. Whenever the variances within strata are small, e.g., in the context of the game defined by (52), the small sample sizes for each stratum are enough for S-LSS to provide accurate stratum estimators, and the variance between the strata is eliminated by algorithm design. Contrary, if the variances within strata are larger, e.g., in the context of the game defined by (55), the reduced amount of samples available for individual strata of the S-LSS algorithm might increase the stratum estimators’ variances more than the stratified algorithm design reduces the overall variance.

This result does not only apply to our chosen sample allocation scheme, i.e., the proportional sample allocation. Clearly, from (60), one obtains that there are too little samples available for the individual stratum estimators in comparison to LSS, and any other sample allocation scheme would also have the drawback that there exists at least one stratum for which there are fewer samples available in comparison to LSS. Additionally, we highlight that even an optimized sample allocation cannot account for the reduced amount of samples available, see Lemma 5 and Corollary 3.

In summary, we recommend using LSS whenever one expects large variances within the strata of S-LSS and a small variance between the stratum estimators. On the other hand, whenever one expects a large variance between the stratum estimators and small variances within strata, we recommend choosing S-LSS. In the latter case, if one has information about the exact variances of individual strata or a good approximation of them, one could distribute the sample sizes across strata similar to the procedure from Lemma 5. Additionally, we highlight that several improvements can be made in order to reduce the mutual shortcomings of both algorithms. One approach is the reuse of samples across strata in order to eliminate the sample size inequality from (60), which was proposed by Benati et al. [14] and introduced as SRS-LSS in SubSection 5.2. However, as we will point out in SubSection 5.4, this approach adds covariance terms to the overall variances of the Shapley value estimators. Another improvement could be some kind of two-stage algorithm, where the first stage is used to decide whether the LSS or S-LSS algorithm should be executed in the second stage.

Finally, we visualize the behavior of LSS and S-LSS when applied to different cooperative games. Figure 2 illustrates the previously obtained results by plotting the variances of LSS and S-LSS against weighted voting games with different quotas. As expected, the variances show a clear dependence to the underlying game, determining whether LSS or S-LSS performs better.

5.4. Analysis of the Variance of SRS-LSS

Since Benati et al. [14] did not address the variance of the SRS-LSS algorithm, we derive:

Proposition 6.

For all $i\in N$, the variance of the SRS-LSS estimator is given by

$$\begin{array}{cc}\hfill Var\left[{\widehat{\varphi }}_i\right]& ={\left({\displaystyle \frac{n-1}{n}}\right)}^2\sum _{s=1}^{n-1}Var\left[{\widehat{a}}_{i,s}\right]\hfill \\ \hfill & \phantom{=}+{\displaystyle \frac{1}{n^2}}\left(\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}\sum _{s=1}^{n-1}Var\left[{\widehat{a}}_{j,s}\right]+2\sum _{\begin{array}{c}j,k\in N\\ j<k\\ j,k\ne i\end{array}}\sum _{s=1}^{n-1}Cov[{\widehat{a}}_{j,s},{\widehat{a}}_{k,s}]\right)\hfill \\ \hfill & \phantom{=}-2{\displaystyle \frac{n-1}{n^2}}\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}\sum _{s=1}^{n-1}Cov[{\widehat{a}}_{i,s},{\widehat{a}}_{j,s}],\hfill \end{array}$$

(61)

with $Var\left[{\widehat{a}}_{j,s}\right]$ given in (36),

$$Cov[{\widehat{a}}_{j,s},{\widehat{a}}_{k,s}]={\displaystyle \frac{1}{{(n-s)}^2}}Cov\left[{\displaystyle \frac{1}{{\widehat{\tau }}_{j,s}}}\sum _{\mathcal{S}\in {\mathbb{S}}_s}{\mathbb{1}}_{j\in \mathcal{S}}v\left(\mathcal{S}\right),{\displaystyle \frac{1}{{\widehat{\tau }}_{k,s}}}\sum _{\mathcal{S}\in {\mathbb{S}}_s}{\mathbb{1}}_{k\in \mathcal{S}}v\left(\mathcal{S}\right)\right],$$

(62)

$${\widehat{\tau }}_{j,s}=\sum _{\mathcal{S}\in {\mathbb{S}}_s}{\mathbb{1}}_{j\in \mathcal{S}}\mathrm{and}{\widehat{\tau }}_{k,s}=\sum _{\mathcal{S}\in {\mathbb{S}}_s}{\mathbb{1}}_{k\in \mathcal{S}},$$

and ${\mathbb{S}}_s{\sim }_{\mathrm{iid}}\mathcal{U}\left(\{S\subseteq N\mid s=|S\left|\right\}\right)$.

Proof. 

Recall the closed-form formulation of the Shapley value as a member of the family of least square values in (27). By using the stratification scheme introduced in (34), one obtains the estimator

$${\widehat{\varphi }}_i={\displaystyle \frac{v\left(N\right)}{n}}+{\displaystyle \frac{n-1}{n}}{\widehat{a}}_i-{\displaystyle \frac{1}{n}}\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}{\widehat{a}}_j$$

with

$${\widehat{a}}_j=\sum _{s=1}^{n-1}{\widehat{a}}_{j,s}.$$

(63)

With these definitions in place, the variance of the Shapley value estimator is given by

$$\begin{array}{cc}\hfill Var\left[{\widehat{\varphi }}_i\right]& =Var\left[{\displaystyle \frac{v\left(N\right)}{n}}+{\displaystyle \frac{n-1}{n}}{\widehat{a}}_i-{\displaystyle \frac{1}{n}}\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}{\widehat{a}}_j\right]\hfill \\ \hfill & ={\left({\displaystyle \frac{n-1}{n}}\right)}^{2}Var\left[{\widehat{a}}_i\right]+{\displaystyle \frac{1}{n^2}}Var\left[\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}{\widehat{a}}_j\right]-2{\displaystyle \frac{n-1}{n^2}}Cov\left[{\widehat{a}}_i,\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}{\widehat{a}}_j\right],\hfill \end{array}$$

(64)

where

$$Var\left[\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}{\widehat{a}}_j\right]=\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}Var\left[{\widehat{a}}_j\right]+2\sum _{\begin{array}{c}j,k\in N\\ j<k\\ j,k\ne i\end{array}}Cov[{\widehat{a}}_j,{\widehat{a}}_k],$$

(65)

and

$$Cov\left[{\widehat{a}}_i,\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}{\widehat{a}}_j\right]=\sum _{\begin{array}{c}j\in N\\ j\ne i\end{array}}Cov[{\widehat{a}}_i,{\widehat{a}}_j].$$

(66)

Then, via (63), we rewrite

$$Cov[{\widehat{a}}_j,{\widehat{a}}_k]=Cov\left[\sum _{s=1}^{n-1}{\widehat{a}}_{j,s},\sum _{s^{\prime }=1}^{n-1}{\widehat{a}}_{k,s^{\prime }}\right]=\sum _{s=1}^{n-1}\sum _{s^{\prime }=1}^{n-1}Cov[{\widehat{a}}_{j,s},{\widehat{a}}_{k,s^{\prime }}],$$

for all $j\ne k\in N$, but since samples are independent for $s\ne s^{\prime }$, the equation can be simplified to

$$Cov[{\widehat{a}}_j,{\widehat{a}}_k]=\sum _{s=1}^{n-1}Cov[{\widehat{a}}_{j,s},{\widehat{a}}_{k,s}].$$

Next, from line 12 in Algorithm 3, we derive

$$\begin{array}{cc}\hfill Cov[{\widehat{a}}_{j,s},{\widehat{a}}_{k,s}]& =Cov\left[{\displaystyle \frac{1}{{\widehat{\tau }}_{j,s}(n-s)}}\sum _{\mathcal{S}\in {\mathbb{S}}_s}{\mathbb{1}}_{j\in \mathcal{S}}v\left(\mathcal{S}\right),{\displaystyle \frac{1}{{\widehat{\tau }}_{k,s}(n-s)}}\sum _{\mathcal{S}\in {\mathbb{S}}_s}{\mathbb{1}}_{k\in \mathcal{S}}v\left(\mathcal{S}\right)\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{{(n-s)}^2}}Cov\left[{\displaystyle \frac{1}{{\widehat{\tau }}_{j,s}}}\sum _{\mathcal{S}\in {\mathbb{S}}_s}{\mathbb{1}}_{j\in \mathcal{S}}v\left(\mathcal{S}\right),{\displaystyle \frac{1}{{\widehat{\tau }}_{k,s}}}\sum _{\mathcal{S}\in {\mathbb{S}}_s}{\mathbb{1}}_{k\in \mathcal{S}}v\left(\mathcal{S}\right)\right]\hfill \end{array}$$

with ${\widehat{\tau }}_{j,s}={\sum }_{\mathcal{S}\in {\mathbb{S}}_s}{\mathbb{1}}_{j\in \mathcal{S}}$, ${\widehat{\tau }}_{k,s}={\sum }_{\mathcal{S}\in {\mathbb{S}}_s}{\mathbb{1}}_{k\in \mathcal{S}}$, and ${\mathbb{S}}_s{\sim }_{\mathrm{iid}}\mathcal{U}\left(\{S\subseteq N\mid s=|S\left|\right\}\right)$.

Finally, setting

$$Var\left[{\widehat{a}}_j\right]=\sum _{s=1}^{n-1}Var\left[{\widehat{a}}_{j,s}\right],$$

since $Cov[{\widehat{a}}_{j,s},{\widehat{a}}_{j,s^{\prime }}]=0$ for $s\ne s^{\prime }$ and inserting everything back into (64) results in (61). □

Corollary 4.

For S-LSS, (61) is equal to (38).

Proof. 

In the case of S-LSS, sampling within each stratum is independent of sampling in all other strata. Therefore, the covariance terms defined by (62) are zero. Consequently, setting these covariance terms to zero in (61) results in (38). □

Theorem 3.

For SRS-LSS, the covariance terms in (61) do not vanish in general.

Proof. 

To simplify notations, we first obtain from (62) that the constant, non-zero term ${(n-s)}^{-2}>0$ can be ignored when demonstrating that $Cov[{\widehat{a}}_{j,s},{\widehat{a}}_{k,s}]\ne 0$ might occur. Therefore, we define

$$Cov[{\overline{a}}_{j,s},{\overline{a}}_{k,s}]=Cov\left[{\displaystyle \frac{1}{{\widehat{\tau }}_{j,s}}}\sum _{\mathcal{S}\in {\mathbb{S}}_s}{\mathbb{1}}_{j\in \mathcal{S}}v\left(\mathcal{S}\right),{\displaystyle \frac{1}{{\widehat{\tau }}_{k,s}}}\sum _{\mathcal{S}\in {\mathbb{S}}_s}{\mathbb{1}}_{k\in \mathcal{S}}v\left(\mathcal{S}\right)\right].$$

Now, we demonstrate by example that $Cov[{\overline{a}}_{j,s},{\overline{a}}_{k,s}]$ may be non-zero. For this purpose, we define a weighted voting game with $\mathbf{c}=[1,0,1]$ and quota $C=2$, such that the only winning coalition of size 2 is $\{1,3\}$. Let us now fix $s=2$, such that the corresponding sample is obtained as ${\mathbb{S}}_2{\sim }_{\mathrm{iid}}\mathcal{U}\left(\left\{\right\{1,2\},\{1,3\},\{2,3\left\}\right\}\right)$. Without loss of generality, we set the sample size to ${\tau }_2=3$.

In Proposition 4, we established that the probability of at least one value ${\widehat{\tau }}_{j,s}$ being zero converges to zero as ${\tau }_s\to \infty$. For ${\tau }_2=3$, this probability is given by

$$\mathbb{P}({\mathcal{S}}_1={\mathcal{S}}_2={\mathcal{S}}_3)=3{\left({\displaystyle \frac{1}{3}}\right)}^3={\displaystyle \frac{1}{9}},$$

(67)

where ${\mathbb{S}}_2={[{\mathcal{S}}_1,{\mathcal{S}}_2,{\mathcal{S}}_3]}^{\top }$. We denote this event by B and its complement by $\overline{B}$, such that $\mathbb{P}\left(B\right)=\mathbb{P}({\mathcal{S}}_1={\mathcal{S}}_2={\mathcal{S}}_3)$.

Although in this example the probability of failing to sample at least one element from each stratum is relatively large, we proceed under the standard assumption that this probability is negligible, which is justified for large ${\tau }_2$. Therefore, we calculate the conditional expectations in Table 2 based on the assumption that B does not occur.

Based on Table 2, the conditional covariance between ${\overline{a}}_{1,2}$ and ${\overline{a}}_{3,2}$ can be calculated as

$$Cov[{\overline{a}}_{1,2},{\overline{a}}_{3,2}\mid \overline{B}]=\mathbb{E}[{\overline{a}}_{1,2}{\overline{a}}_{3,2}\mid \overline{B}]-\mathbb{E}[{\overline{a}}_{1,2}\mid \overline{B}]\mathbb{E}[{\overline{a}}_{3,2}\mid \overline{B}]={\displaystyle \frac{5}{16}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2}}={\displaystyle \frac{1}{16}}>0.$$

Although conditioning on $\overline{B}$ is necessary for computing the covariance, it may slightly affect its numerical value. The true (unconditional) covariance could differ marginally. Accounting for the undefined cases would require specifying a particular strategy to handle such edge cases, which was not proposed for the SRS-LSS algorithm and lies beyond the scope of this investigation. Nevertheless, the key point remains that the covariance does not vanish. As can be observed from the ${\overline{a}}_{1,2}$ and ${\overline{a}}_{3,2}$ columns in Table 2, these quantities have a clear positive correlation. Since players 1 and 3 together form the only winning coalition of size 2, large realizations of ${\overline{a}}_{1,2}$ necessarily coincide with large realizations of ${\overline{a}}_{3,2}$. This dependence confirms that $Cov[{\overline{a}}_{1,2},{\overline{a}}_{3,2}]>0$, and as a result, we have $Cov[{\widehat{a}}_{1,2},{\widehat{a}}_{3,2}]\ne 0$. □

From Theorem 3, it follows that reusing evaluations of the characteristic function introduces non-zero covariance terms. These terms are hard to characterize, which complicates the exact computation of the SRS-LSS estimator’s theoretical variance and its comparison to LSS and S-LSS. A detailed analysis of (62) together with the resulting influence on (61) might be a question for further research. In the meantime, we rely on empirical observations of the SRS-LSS estimator’s variance in the numerical experiments in Section 6, which indicate a significantly reduced variance of SRS-LSS in comparison to LSS and S-LSS.

5.5. Comparison of LSS and UKS

In the following, we compare LSS (Algorithm 1) with UKS (see the end of SubSection 5.2). First, we show how UKS fits into the general importance sampling framework defined by Theorem 1 and discuss the relationship between both algorithms within this context. Subsequently, we compare their variances.

Before going into detail, we note that the formulation of UKS derived in the following was already proposed by Fumagalli et al. [30]. The authors identify UKS as a special case of their SHAPley Interaction Quantification framework, which enables a reformulation of the algorithm similar to ours. Nevertheless, we present our derivations for several reasons. First, our derivation of UKS is, from our perspective, more intuitive. Second, while Fumagalli et al. [30] suggest that their framework enables a particular formulation of UKS that arises specifically from their approach, this form of UKS can also be directly obtained by expanding the matrices and vectors in (46) without using the indirection via their SHAPley Interaction Quantification framework. Third, UKS has not yet been formally integrated into the importance sampling framework defined by Theorem 1, which allows us to directly apply our results from Proposition 1. Finally, and most importantly, to the best of our knowledge, we are the first to demonstrate that UKS approximates the same underlying problem as LSS, establishing a novel link between the two algorithms.

We refer to Appendix A, which shows that our following results are consistent with those of Fumagalli et al. [30].

Algorithm Comparison

First, recall the definition of LSS given in (30) with $w_i$ being defined in (24), m in (25), $\gamma =1$ in (26), and p in (28). Since we have only derived $w_i\left(S\right)/p\left(S\right)$ so far and not $w_i\left(S\right)$ directly, inserting m into $w_i$ results in

$$\begin{array}{cc}\hfill {\displaystyle \frac{n-\left|S\right|}{n(n-1)\left(\genfrac{}{}{0pt}{}{n-2}{\left|S\right|-1}\right)}}& ={\displaystyle \frac{(n-|S\left|\right)\left(\right|S|-1)!(n-|S|-1)!}{n(n-1)(n-2)!}}={\displaystyle \frac{\left(\right|S|-1)!(n-|S\left|\right)!}{n!}}\hfill \\ \hfill & ={\displaystyle \frac{1}{\left|S\right|}}{\displaystyle \frac{\left|S\right|!(n-|S\left|\right)!}{n!}}={\displaystyle \frac{1}{\left|S\right|\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\hfill \end{array}$$

for the case $i\in S$, and

$$\begin{array}{cc}\hfill -{\displaystyle \frac{\left|S\right|}{n(n-1)\left(\genfrac{}{}{0pt}{}{n-2}{\left|S\right|-1}\right)}}& =-{\displaystyle \frac{\left|S\right|\left(\right|S|-1)!(n-|S|-1)!}{n(n-1)(n-2)!}}=-{\displaystyle \frac{\left|S\right|!(n-|S|-1)!}{n!}}\hfill \\ \hfill & =-{\displaystyle \frac{1}{n-\left|S\right|}}{\displaystyle \frac{\left|S\right|!(n-|S\left|\right)!}{n!}}=-{\displaystyle \frac{1}{(n-|S\left|\right)\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\hfill \end{array}$$

for the case $i\notin S$. Thus, $w_i$ can be rewritten as

$$w_i^{\mathrm{LSS}}\left(S\right)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{\left|S\right|\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\hfill & \mathrm{if}i\in S\hfill \\ -{\displaystyle \frac{1}{(n-|S\left|\right)\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\hfill & \mathrm{if}i\notin S.\hfill \end{array}\right.$$

(68)

Now, let us turn our attention to the UKS algorithm:

Proposition 7.

UKS can be brought into the form

$${\widehat{\varphi }}_i=t+{\displaystyle \frac{1}{\tau }}\sum _{\mathcal{S}\in {\mathbb{S}}^{\mathrm{UKS}}}{\displaystyle \frac{w_i^{\mathrm{UKS}}\left(\mathcal{S}\right)}{p^{\mathrm{UKS}}\left(\mathcal{S}\right)}}v\left(\mathcal{S}\right)$$

with ${\mathbb{S}}^{\mathrm{UKS}}{\sim }_{\mathrm{iid}}{p}^{\mathrm{UKS}}$ being a sample according to (41),

$$t={\left[A^{-1}\mathbf{1}{\displaystyle \frac{v\left(N\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right]}_i$$

being a constant, and

$${\displaystyle \frac{w_i^{\mathrm{UKS}}\left(S\right)}{p^{\mathrm{UKS}}\left(S\right)}}={\left[A^{-1}\mathbf{z}\left(S\right)\right]}_i-{\left[A^{-1}\mathbf{1}\right]}_i{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{z}\left(S\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}.$$

(69)

Proof. 

Starting from (46) and substituting $\widehat{\mathbf{b}}$ as defined in (45) results in

$$\begin{array}{cc}\hfill {\widehat{\varphi }}_i& ={\left[A^{-1}\left(\widehat{\mathbf{b}}-\mathbf{1}{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\widehat{\mathbf{b}}-v\left(N\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right)\right]}_i\hfill \\ \hfill & ={\left[A^{-1}\widehat{\mathbf{b}}-A^{-1}\mathbf{1}{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\widehat{\mathbf{b}}}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}+A^{-1}\mathbf{1}{\displaystyle \frac{v\left(N\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right]}_i\hfill \\ \hfill & =\underset{t}{\underbrace{{\left[A^{-1}\mathbf{1}{\displaystyle \frac{v\left(N\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right]}_i}}+{\left[A^{-1}\widehat{\mathbf{b}}-A^{-1}\mathbf{1}{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\widehat{\mathbf{b}}\right]}_i\hfill \\ \hfill & =t+{\left[\left(A^{-1}-A^{-1}\mathbf{1}{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right)\widehat{\mathbf{b}}\right]}_i\hfill \\ \hfill & =t+{\left[\left(A^{-1}-A^{-1}\mathbf{1}{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right){\displaystyle \frac{1}{\tau }}\sum _{\mathcal{S}\in {\mathbb{S}}^{\mathrm{UKS}}}\mathbf{z}\left(\mathcal{S}\right)v\left(\mathcal{S}\right)\right]}_i\hfill \\ \hfill & =t+{\left[{\displaystyle \frac{1}{\tau }}\sum _{\mathcal{S}\in {\mathbb{S}}^{\mathrm{UKS}}}\left(A^{-1}-A^{-1}\mathbf{1}{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right)\mathbf{z}\left(\mathcal{S}\right)v\left(\mathcal{S}\right)\right]}_i\hfill \\ \hfill & =t+{\left[{\displaystyle \frac{1}{\tau }}\sum _{\mathcal{S}\in {\mathbb{S}}^{\mathrm{UKS}}}\left(A^{-1}\mathbf{z}\left(\mathcal{S}\right)-A^{-1}\mathbf{1}{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{z}\left(\mathcal{S}\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right)v\left(\mathcal{S}\right)\right]}_i\hfill \\ \hfill & =t+{\displaystyle \frac{1}{\tau }}\sum _{\mathcal{S}\in {\mathbb{S}}^{\mathrm{UKS}}}\underset{w_i^{\mathrm{UKS}}\left(\mathcal{S}\right)/p^{\mathrm{UKS}}\left(\mathcal{S}\right)}{\underbrace{\left({\left[A^{-1}\mathbf{z}\left(\mathcal{S}\right)\right]}_i-{\left[A^{-1}\mathbf{1}\right]}_i{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{z}\left(\mathcal{S}\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right)}}v\left(\mathcal{S}\right).\hfill \end{array}$$

(70)

□

We now aim to show that both algorithms share the same weight function w and the same initial term, i.e., $w_i^{\mathrm{LSS}}\left(S\right)=w_i^{\mathrm{UKS}}\left(S\right)$ and $t={\displaystyle \frac{v\left(N\right)}{n}}$. To prepare for this, we first carry out some reformulations in order to simplify the upcoming calculations. Readers only interested in the final result may skip ahead to Theorem 4.

Lemma 6.

The expression for Q given in (43) can be rewritten as

$$Q={\displaystyle \frac{2(n-1)H_{n-1}}{n}},$$

(71)

where $H_{n-1}$ denotes the $(n-1)$-th harmonic number.

Proof. 

The proof is straightforward by inserting the identity ${\sum }_{s=1}^{n-1}{\displaystyle \frac{1}{s(n-s)}}={\displaystyle \frac{2H_{n-1}}{n}}$ with $H_{n-1}$ being the $(n-1)$-th harmonic number [30] into (43), which results in (71). □

Lemma 7.

The inverse of A as defined in (47) is given by

$$A^{-1}=\left[\begin{array}{cccc}2H_{n-1}+\overline{\rho }& \overline{\rho }& \dots & \overline{\rho }\\ \overline{\rho }& 2H_{n-1}+\overline{\rho }& \dots & \overline{\rho }\\ ⋮& ⋮& \ddots & ⋮\\ \overline{\rho }& \overline{\rho }& \dots & 2H_{n-1}+\overline{\rho }\end{array}\right]$$

with

$$\overline{\rho }=-{\displaystyle \frac{2H_{n-1}(H_{n-1}-1)}{1+n(H_{n-1}-1)}}.$$

Proof. 

The article by Fumagalli et al. [30] shows that the off-diagonal entries of the inverse of A as defined in (47) are given by

$$\overline{\rho }=-2H_{n-1}{\displaystyle \frac{\frac{H_{n-1}-1}{2H_{n-1}}}{\frac{1}{2}+(n-1)\frac{H_{n-1}-1}{2H_{n-1}}}}.$$

Its denominator can be reformulated as

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{2}}+(n-1){\displaystyle \frac{H_{n-1}-1}{2H_{n-1}}}& ={\displaystyle \frac{H_{n-1}+(n-1)(H_{n-1}-1)}{2H_{n-1}}}\hfill \\ \hfill & ={\displaystyle \frac{1+(H_{n-1}-1)+(n-1)(H_{n-1}-1)}{2H_{n-1}}}\hfill \\ \hfill & ={\displaystyle \frac{1+n(H_{n-1}-1)}{2H_{n-1}}},\hfill \end{array}$$

(72)

and thus,

$$\overline{\rho }=-2H_{n-1}{\displaystyle \frac{\frac{H_{n-1}-1}{2H_{n-1}}}{\frac{1+n(H_{n-1}-1)}{2H_{n-1}}}}=-{\displaystyle \frac{2H_{n-1}(H_{n-1}-1)}{1+n(H_{n-1}-1)}}.$$

For the diagonal entries, Fumagalli et al. [30] provides

$$2H_{n-1}{\displaystyle \frac{\frac{1}{2}+(n-2)\frac{H_{n-1}-1}{2H_{n-1}}}{\frac{1}{2}+(n-1)\frac{H_{n-1}-1}{2H_{n-1}}}},$$

and similarly to (72), its numerator can be derived as

$${\displaystyle \frac{1}{2}}+(n-2){\displaystyle \frac{H_{n-1}-1}{2H_{n-1}}}={\displaystyle \frac{1+(n-1)(H_{n-1}-1)}{2H_{n-1}}}.$$

(73)

Putting (72) and (73) together, we retrieve

$$\begin{array}{cc}\hfill 2H_{n-1}{\displaystyle \frac{\frac{1}{2}+(n-2)\frac{H_{n-1}-1}{2H_{n-1}}}{\frac{1}{2}+(n-1)\frac{H_{n-1}-1}{2H_{n-1}}}}& =2H_{n-1}{\displaystyle \frac{\frac{1+(n-1)(H_{n-1}-1)}{2H_{n-1}}}{\frac{1+n(H_{n-1}-1)}{2H_{n-1}}}}\hfill \\ \hfill & =2H_{n-1}{\displaystyle \frac{1+(n-1)(H_{n-1}-1)}{1+n(H_{n-1}-1)}}\hfill \\ \hfill & =2H_{n-1}-2H_{n-1}\left(1-{\displaystyle \frac{1+(n-1)(H_{n-1}-1)}{1+n(H_{n-1}-1)}}\right)\hfill \\ \hfill & =2H_{n-1}-2H_{n-1}{\displaystyle \frac{2H_{n-1}(1+n(H_{n-1}-1)-1-(n-1)(H_{n-1}-1))}{1+n(H_{n-1}-1)}}\hfill \\ \hfill & =2H_{n-1}-{\displaystyle \frac{2H_{n-1}(H_{n-1}-1)}{1+n(H_{n-1}-1)}}=2H_{n-1}+\overline{\rho }.\hfill \end{array}$$

□

Proposition 8.

The constant terms of both algorithms are equal, i.e.,

$$t={\displaystyle \frac{v\left(N\right)}{n}}.$$

Proof. 

From (47) one directly obtains that $\mathbf{1}$ is an eigenvector of A with eigenvalue $\lambda =\frac{1}{2}+(n-1)\rho$. Thus, we have $A^{-1}\mathbf{1}={\lambda }^{-1}\mathbf{1}$. With that given, rewriting t from Proposition 7 results in

$$\begin{array}{cc}\hfill t& ={\left[A^{-1}\mathbf{1}{\displaystyle \frac{v\left(N\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right]}_i\hfill \\ \hfill & ={\left[A^{-1}\mathbf{1}{\displaystyle \frac{1}{{\mathbf{1}}^{\top }\left(A^{-1}\mathbf{1}\right)}}\right]}_{i}v\left(N\right)={\left[{\lambda }^{-1}\mathbf{1}{\displaystyle \frac{1}{{\mathbf{1}}^{\top }\left({\lambda }^{-1}\mathbf{1}\right)}}\right]}_{i}v\left(N\right)\hfill \\ \hfill & ={\left[{\lambda }^{-1}\mathbf{1}{\displaystyle \frac{1}{n{\lambda }^{-1}}}\right]}_{i}v\left(N\right)={\displaystyle \frac{{\lambda }^{-1}}{n{\lambda }^{-1}}}v\left(N\right)={\displaystyle \frac{v\left(N\right)}{n}}.\hfill \end{array}$$

□

Lemma 8.

For any $S\subseteq N$ and $i\in N$, there holds

$${\left[A^{-1}\mathbf{1}\right]}_i{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{z}\left(S\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}={\displaystyle \frac{2\left|S\right|H_{n-1}}{n}}+\left|S\right|\overline{\rho }.$$

Proof. 

Using Lemma 7, we derive

$${\mathbf{1}}^{\top }{A}^{-1}=\left[2H_{n-1}+n\overline{\rho },\dots ,2H_{n-1}+n\overline{\rho }\right]\in {\mathbb{R}}^{1\times n},$$

and thus,

$${\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{z}\left(S\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}={\displaystyle \frac{\left|S\right|(2H_{n-1}+n\overline{\rho })}{n(2H_{n-1}+n\overline{\rho })}}={\displaystyle \frac{\left|S\right|}{n}}.$$

Furthermore, also based on Lemma 7,

$${\left[A^{-1}\mathbf{1}\right]}_i=2H_{n-1}+n\overline{\rho },$$

and therefore,

$${\left[A^{-1}\mathbf{1}\right]}_i{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{z}\left(S\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}={\displaystyle \frac{2\left|S\right|H_{n-1}}{n}}+{\displaystyle \frac{\left|S\right|n\overline{\rho }}{n}}={\displaystyle \frac{2\left|S\right|H_{n-1}}{n}}+\left|S\right|\overline{\rho }.$$

□

Proposition 9.

For all $S\subseteq N$ with $0<\left|S\right|<n$ and $i\in S$,

$$w_i^{\mathrm{UKS}}\left(S\right)={\displaystyle \frac{1}{\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)\left|S\right|}}.$$

Proof. 

From Proposition 7, we obtain

$${\displaystyle \frac{w_i^{\mathrm{UKS}}\left(S\right)}{p^{\mathrm{UKS}}\left(S\right)}}={\left[A^{-1}\mathbf{z}\left(S\right)\right]}_i-{\left[A^{-1}\mathbf{1}\right]}_i{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{z}\left(S\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}.$$

(74)

Using Lemma 7, we derive

$${\left[A^{-1}\mathbf{z}\left(S\right)\right]}_i=2H_{n-1}+\left|S\right|\overline{\rho }.$$

(75)

Rewriting the Shapley kernel from (42) to follow our general notation for cooperative games results in

$$\omega \left(S\right)={\displaystyle \frac{n-1}{\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)\left|S\right|(n-\left|S\right|)}}.$$

(76)

Combining the results from Lemmas 6 and 8 as well as (41), (74), (75), and (76) results in

$$\begin{array}{cc}\hfill w_i^{\mathrm{UKS}}\left(S\right)& =p^{\mathrm{UKS}}\left(S\right){\displaystyle \frac{w_i^{\mathrm{UKS}}\left(S\right)}{p^{\mathrm{UKS}}\left(S\right)}}\hfill \\ \hfill & =Q^{-1}\omega \left(S\right)\left({\left[A^{-1}\mathbf{z}\left(S\right)\right]}_i-{\left[A^{-1}\mathbf{1}\right]}_i{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{z}\left(S\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right)\hfill \\ \hfill & ={\displaystyle \frac{n}{2(n-1)H_{n-1}}}{\displaystyle \frac{n-1}{\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)\left|S\right|(n-\left|S\right|)}}\left(2H_{n-1}+\left|S\right|\overline{\rho }-{\displaystyle \frac{2\left|S\right|H_{n-1}}{n}}-\left|S\right|\overline{\rho }\right)\hfill \\ \hfill & ={\displaystyle \frac{n}{2H_{n-1}}}{\displaystyle \frac{1}{\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)\left|S\right|(n-\left|S\right|)}}\left(2H_{n-1}-{\displaystyle \frac{2\left|S\right|H_{n-1}}{n}}\right)\hfill \\ \hfill & ={\displaystyle \frac{n}{2H_{n-1}}}{\displaystyle \frac{1}{\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)\left|S\right|(n-\left|S\right|)}}{\displaystyle \frac{2(n-|S\left|\right)H_{n-1}}{n}}={\displaystyle \frac{1}{\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)\left|S\right|}}.\hfill \end{array}$$

□

Proposition 10.

For all $S\subseteq N$ with $0<\left|S\right|<n$ and $i\notin S$,

$$w_i^{\mathrm{UKS}}\left(S\right)=-{\displaystyle \frac{1}{\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)(n-|S\left|\right)}}.$$

Proof. 

From Proposition 7, we have

$${\displaystyle \frac{w_i^{\mathrm{UKS}}\left(S\right)}{p^{\mathrm{UKS}}\left(S\right)}}={\left[A^{-1}\mathbf{z}\left(S\right)\right]}_i-{\left[A^{-1}\mathbf{1}\right]}_i{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{z}\left(S\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}.$$

(77)

Using Lemma 7, we derive

$${\left[A^{-1}\mathbf{z}\left(S\right)\right]}_i=\left|S\right|\overline{\rho }.$$

(78)

Combining the results from Lemmas 6 and 8 as well as (41), (76), (77), and (78) results in

$$\begin{array}{cc}\hfill w_i^{\mathrm{UKS}}\left(S\right)& =p^{\mathrm{UKS}}\left(S\right){\displaystyle \frac{w_i^{\mathrm{UKS}}\left(S\right)}{p^{\mathrm{UKS}}\left(S\right)}}\hfill \\ \hfill & =Q^{-1}\omega \left(S\right)\left({\left[A^{-1}\mathbf{z}\left(S\right)\right]}_i-{\left[A^{-1}\mathbf{1}\right]}_i{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}z\left(S\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right)\hfill \\ \hfill & ={\displaystyle \frac{n}{2(n-1)H_{n-1}}}{\displaystyle \frac{n-1}{\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)\left|S\right|(n-\left|S\right|)}}\left(\left|S\right|\overline{\rho }-{\displaystyle \frac{2\left|S\right|H_{n-1}}{n}}-\left|S\right|\overline{\rho }\right)\hfill \\ \hfill & =-{\displaystyle \frac{n}{2H_{n-1}}}{\displaystyle \frac{1}{\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)\left|S\right|(n-\left|S\right|)}}{\displaystyle \frac{2\left|S\right|H_{n-1}}{n}}=-{\displaystyle \frac{1}{\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)(n-|S\left|\right)}}.\hfill \end{array}$$

□

Theorem 4.

The estimators obtained via LSS and UKS are both importance sampling estimators for the same linear solution concept in a sense of Theorem 1, differing only in their respective sampling strategies.

Proof. 

Per definition given in (30), LSS has the form

$${\widehat{\varphi }}_i={\displaystyle \frac{v\left(N\right)}{n}}+{\displaystyle \frac{1}{\tau }}\sum _{\mathcal{S}\in \mathbb{S}}{\displaystyle \frac{w_i\left(\mathcal{S}\right)}{p\left(\mathcal{S}\right)}}v\left(\mathcal{S}\right),$$

with weights given by (68).

Proposition 7 shows that UKS can also be written in this form, sharing the same constant term $t={\displaystyle \frac{v\left(N\right)}{n}}$, see Proposition 8. Moreover, by combining propositions 9 and 10, one obtains

$$w_i^{\mathrm{UKS}}\left(S\right)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{\left|S\right|\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\hfill & \mathrm{if}i\in S\hfill \\ -{\displaystyle \frac{1}{(n-|S\left|\right)\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\hfill & \mathrm{if}i\notin S\hfill \end{array}\right.$$

(79)

for $0<\left|S\right|<n$, which is identical to the weight function of LSS defined in (68).

Therefore, we conclude that both algorithms are importance sampling estimators of the same problem

$${\varphi }_i={\displaystyle \frac{v\left(N\right)}{n}}+\sum _{\begin{array}{c}S\subseteq N\\ 0<\left|S\right|<n\end{array}}{w}_i\left(S\right)v\left(S\right)$$

with $w_i\left(S\right)=w_i^{\mathrm{LSS}}\left(S\right)=w_i^{\mathrm{UKS}}\left(S\right)$.

Finally, we highlight that both algorithms differ only in their respective sampling distributions, since (28) and the distribution that is proportional to (76) are clearly different sampling distributions. □

We illustrate the result from Theorem 4 in Figure 3 which visualizes the different sampling strategies of LSS and UKS for $n=9$.

Variance Comparison

In the previous section, we showed that UKS fits the importance sampling framework introduced by Theorem 1, see Proposition 7 and Theorem 4. Therefore, although Covert and Lee [7] already provided (49) as the variance for the UKS estimator, we can now use a different formulation:

Proposition 11.

The variance of the UKS estimator is given by

$$Var\left[{\widehat{\varphi }}_i\right]={\displaystyle \frac{1}{\tau }}\left(\sum _{\begin{array}{c}S\subseteq N\\ 0<\left|S\right|<n\end{array}}{\displaystyle \frac{w_i{\left(S\right)}^2}{p\left(S\right)}}v{\left(S\right)}^2-{\left({\varphi }_i-{\displaystyle \frac{v\left(N\right)}{n}}\right)}^2\right)(\forall i\in N)$$

with

$${\displaystyle \frac{w_i{\left(S\right)}^2}{p\left(S\right)}}=\left\{\begin{array}{cc}{\displaystyle \frac{2(n-|S\left|\right)H_{n-1}}{n\left|S\right|\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\hfill & \mathrm{if}i\in S\hfill \\ {\displaystyle \frac{2\left|S\right|H_{n-1}}{n(n-|S\left|\right)\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\hfill & \mathrm{if}i\notin S.\hfill \end{array}\right.$$

(80)

Proof. 

First, let us analyze the case $i\in S$. Using Lemma 8 as well as (69), (75), and (79), we derive

$$\begin{array}{cc}\hfill {\displaystyle \frac{w{\left(S\right)}^2}{p\left(S\right)}}& =w\left(S\right){\displaystyle \frac{w\left(S\right)}{p\left(S\right)}}\hfill \\ \hfill & ={\displaystyle \frac{1}{\left|S\right|\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\left({\left[A^{-1}\mathbf{z}\left(S\right)\right]}_i-{\left[A^{-1}\mathbf{1}\right]}_i{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{z}\left(S\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right)\hfill \\ \hfill & ={\displaystyle \frac{1}{\left|S\right|\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\left(2H_{n-1}+\left|S\right|\overline{\rho }-{\displaystyle \frac{2\left|S\right|H_{n-1}}{n}}-\left|S\right|\overline{\rho }\right)\hfill \\ \hfill & ={\displaystyle \frac{1}{\left|S\right|\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\left(2H_{n-1}-{\displaystyle \frac{2\left|S\right|H_{n-1}}{n}}\right)={\displaystyle \frac{2(n-|S\left|\right)H_{n-1}}{n\left|S\right|\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}.\hfill \end{array}$$

Similarly, for the case $i\notin S$, by using Lemma 8 as well as (69), (78), and (79), we obtain

$$\begin{array}{cc}\hfill {\displaystyle \frac{w{\left(S\right)}^2}{p\left(S\right)}}& =w\left(S\right){\displaystyle \frac{w\left(S\right)}{p\left(S\right)}}\hfill \\ \hfill & =-{\displaystyle \frac{1}{(n-|S\left|\right)\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\left({\left[A^{-1}\mathbf{z}\left(S\right)\right]}_i-{\left[A^{-1}\mathbf{1}\right]}_i{\displaystyle \frac{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{z}\left(S\right)}{{\mathbf{1}}^{\top }{A}^{-1}\mathbf{1}}}\right)\hfill \\ \hfill & =-{\displaystyle \frac{1}{(n-|S\left|\right)\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}\left(\left|S\right|\overline{\rho }-{\displaystyle \frac{2\left|S\right|H_{n-1}}{n}}-\left|S\right|\overline{\rho }\right)={\displaystyle \frac{2\left|S\right|H_{n-1}}{n(n-|S\left|\right)\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}.\hfill \end{array}$$

Combining both cases results in (80), and since UKS fits the framework introduced by Theorem 1, we can use (18) to compute its variance. Note that we construct a temporary solution concept ${\varphi }_i-{\displaystyle \frac{v\left(N\right)}{n}}$ that excludes ${\displaystyle \frac{v\left(N\right)}{n}}$, similar to the variance of LSS given by Proposition 2. □

In Figure 4, we use Proposition 11 to compute the variances of UKS and Proposition 2 for those of LSS. As expected, neither LSS nor UKS outperforms the other, since they implement the same algorithm while differing only in their respective sampling strategies. The effectiveness of any given sampling strategy depends on the characteristics of the underlying problem, such as the parameters of the weighted voting game. Obviously, there exist problems for which sampling strategies other than the one employed in LSS or UKS achieve superior performance. More broadly, this illustrates that no single sampling strategy is universally optimal, and selecting the most effective one requires tailoring it to the specific problem.

6. Empirical Results

We validate our results by applying the algorithms introduced in SubSection 5.2 to approximate Shapley values for airport games, weighted voting games, and three real-world interpretable machine learning problems.

The algorithms introduced in SubSection 5.2 were implemented in Python. Both the implementations and test problems can be freely accessed via the first author’s GitHub repository

https://github.com/tim-pollmann/shapley-least-squares (accessed on 27 January 2026).

We consistently assume that all players’ Shapley values must be estimated. Although not universal, this is common in explainable machine learning, where one seeks to explain a prediction using all features.

Before proceeding, we note that each algorithm uses different parameters to determine the number of sampled coalitions and thus, evaluations of the characteristic function v. For fair comparisons, we introduce a unified sample budget T, ensuring each algorithm performs T evaluations of v, up to negligible rounding errors.

The overall sample budget T, introduced above, sets the total number of v evaluations per algorithm. We now express the parameters of each algorithm from SubSection 5.2 in terms of T.

For LSS (Algorithm 1) and UKS (see the end of SubSection 5.2), this is a one-to-one mapping, i.e., $\tau =T$. Note that we ignore the single evaluation of v needed for calculating the initial term ${\displaystyle \frac{v\left(N\right)}{n}}$, since it is negligible in comparison to the rounding errors of other algorithms, especially for large T.

S-LSS (Algorithm 2) divides the sample space based on $n-1$ distinct values of s, for each $i\in N$. We use the heuristically motivated proportional sample allocation from (58) to obtain

$${\tau }_{i,s}=⌈{\displaystyle \frac{2sT}{n^2(n-1)}}⌉(\forall i\in N,\forall s\in \{1,\dots ,n-1\}).$$

Finally, SRS-LSS (Algorithm 3) divides the sample space into $n-1$ distinct sets. Thus, for simplicity, noting that alternative sample allocation schemes across strata may provide better results, we set

$${\tau }_s=⌈{\displaystyle \frac{T}{n-1}}⌉(\forall s\in \{1,\dots ,n-1\}).$$

We conclude that the deviation from the true total sample budget T is 1 for LSS and UKS, while the upper bounds of their deviations are given by $n^2-n-1$ for S-LSS and $n-2$ for SRS-LSS. We consider these deviations to be negligible in our subsequent analysis, in particular for large T.

Note that for SRS-LSS (Algorithm 3), it is not guaranteed that the algorithm runs successfully in the sense that every stratum receives at least one sample, compare Proposition 4. In our mean squared error comparisons, we mandate for any $\tau$ that at least half of all runs must be successful for the results to be displayed in the final figure.

With these definitions established, we proceeded with our experiments.

Figure 5 confirms the theoretical variances of LSS and S-LSS presented in Propositions 2 and 3, respectively. Moreover, it underscores our results from Theorem 2, where we demonstrated that one or the other algorithm might achieve a smaller variance, depending on the underlying cooperative game. Furthermore, we validate that LSS and UKS perform as expected with respect to the theoretical variances presented in Propositions 2 and 11 as Figure 5 clearly illustrates that the derived theoretical variances are accurate and that, depending on the specific problem, either LSS or UKS may slightly outperform the other, as one sampling strategy might be better suited to the given task than the other method, see Theorem 4.

As for mean squared error comparisons, we first look at an airport game with 100 players specified in Castro et al. [31] and compare the approximation methods from SubSection 5.2 in Figure 6. Note that the graph for SRS-LSS begins only at a sample size of 50,000, as Proposition 4 fails too frequently with smaller sample budgets. With $n=100$ players, guaranteeing at least one sample per stratum becomes more challenging.

Figure 7 compares mean squared errors for our Monte Carlo estimators from SubSection 5.2 for a weighted voting game with 50 players and normally distributed weights, which was previously employed in the software EPIC [22,23]. It can be found on the GitHub page of the second author via https://github.com/jhstaudacher/EPIC/blob/master/test_cases/normal_sqrd/normal_sqrd.n50.q249646.csv.

Figure 8 compares mean squared errors for our Monte Carlo estimators from SubSection 5.2 for a weighted voting game with 150 players and uniformly distributed weights, which was previously employed in the software EPIC [22,23]. It can be found on the GitHub page of the second author via https://github.com/jhstaudacher/EPIC/blob/master/test_cases/uniform/uniform.n150.q35951.csv. Note that the graph for SRS-LSS begins only at a sample size of 150,000, as Proposition 4 fails too frequently with smaller sample budgets. With $n=150$ players, guaranteeing at least one sample per stratum becomes even more challenging than in the airport game from Figure 6.

Figure 9 uses the standard diabetes dataset (442 patients, 10 baseline features) where the target is disease progression after one year. We train a Gradient Boosting Regressor and evaluate the approximation methods from SubSection 5.2 by their mean squared error against exact reference Shapley values.

In Figure 10, the California housing dataset (20,640 entries, 8 features) predicts median house value in hundreds of thousands of dollars. Using an MLP Regressor, we compare approximation methods by their mean squared error on Shapley value estimation against exact references.

Figure 11 uses the classic wine dataset (178 instances, 13 features, 3 wine classes). We train a Random Forest Classifier and evaluate the approximation methods from SubSection 5.2 via mean squared error in estimating Shapley values for predicting the probability of class 0 only.

Let us summarize the comparisons of mean squared errors of all algorithms for approximating the Shapley value from SubSection 5.2 in Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11 succinctly. As expected, LSS and UKS perform more or less equally for all six test problems. SRS-LSS outperforms the other three provably unbiased algorithms LSS, S-LSS and UKS for all test problems. This observation is consistent with the results shown in Figure 1 in Benati et al. [14], where only LSS and SRS-LSS were compared. Therefore, we conclude that the covariance terms established in Theorem 3 do not significantly affect the overall variances of individual players’ Shapley value estimators negatively. S-LSS performs worst for five test problems with the notable exception of the airport game with 100 players in Figure 6 where S-LSS is even faster than KernelSHAP (KS). Not surprisingly, KS is consistently faster than UKS and LSS. For the three large cooperative games in Figure 6, Figure 7 and Figure 8 SRS-LSS converges faster than KS, whereas that comparison reverses for the machine learning tasks in Figure 9, Figure 10 and Figure 11.

7. Summary, Conclusions and Outlook

The celebrated KernelSHAP approach by Lundberg and Lee [1] as well as its later refinement Unbiased KernelSHAP (UKS) proposed by Covert and Lee [7] compute the Shapley value as a least squares optimization problem. While these two algorithms are extremely well established in the machine learning community, the methods from the paper by Benati et al. [14] — which are also based on the least square formula for the Shapley value — have received fairly little attention. To our knowledge, we are the first to apply them in the context of explainable artificial intelligence. We formulate the ideas from Benati et al. [14] as approximation algorithms for general TU games. We refer to their weighted sampling strategy as LSS (Least Square Sampling). As for their approach towards stratification, we distinguish S-LSS with no reuse of samples across strata and SRS-LSS with sample reuse across strata as proposed in Benati et al. [14].

As the result of our thorough and detailed analysis of LSS, S-LSS and SRS-LSS, we presented three key findings.

First, in Proposition 5, we showed that S-LSS as proposed in [14] is not a valid stratified variant of LSS in the sense of the definition provided SubSection 3.2,i.e., for S-LSS the strata overlap. Therefore, we demonstrated that S-LSS might reduce the variance of the obtained estimator, but it could also lead to an increase in comparison to LSS, see Theorem 2.

Second, in Theorem 3, we showed that the SRS-LSS approach proposed by Benati et al. [14] introduces covariance terms between stratum estimators, making its theoretical variance difficult to analyze. Although empirical results suggest that the variance is significantly reduced in comparison to LSS and S-LSS, a theoretical analysis remains an open research question.

Third, in Theorem 4, we established that LSS and UKS are importance sampling estimators in the sense of Theorem 1, addressing the same underlying problem but differing in their respective sampling strategies. Therefore, neither LSS nor UKS is superior over the other one. Which algorithm’s variance is smaller depends on the underlying problem and whether the respective sampling strategy is suitable for this problem or not.

As noted in the introduction, this paper’s aim is neither to propose a new Shapley value estimator nor to compare numerous approximation algorithms. Rather, the focus is on providing structural, theoretical, and algorithmic insight into the methods from Benati et al. [14]. While antithetic sampling — as discussed, for instance, in [32] — could be integrated into the SRS-LSS algorithm, this would have vastly exceeded the scope of this study. In the future, it could be worthwhile to investigate whether there are algorithmic approaches comparable to sophisticated stratification strategies [28,33,34] which could be successfully incorporated to enhance the performance of SRS-LSS. Definitely, our theoretical and numerical results suggest that the SRS-LSS estimator — which is unbiased — justifies more attention. Another, perhaps even more pressing question, is why SRS-LSS outperforms KernelSHAP for weighted voting games and airport games while KernelSHAP exhibits superior performance for our test problems from interpretable machine learning. We have yet to identify the properties of the characteristic function responsible for this phenomenon.