On Importance Sampling and Multilinear Extensions for Approximating Shapley Values with Applications to Explainable Artificial Intelligence

1. Introduction

Cooperative game theory [1] models situations where entities (players) collaborate to achieve a common objective. Example scenarios (games) include voting situations in parliaments [2,3,4], sharing of infrastructure costs [5,6,7] or the analysis of genetic networks [8,9]. More recently, cooperative game theory has garnered plenty of interest in the field of explainable artificial intelligence (XAI) and interpretable machine learning [10,11,12,13], where the players are the features of a machine learning model contributing to a certain prediction. In each case, cooperative game theory provides tools — primarily in the form of point-valued solution concepts — to analyze how the total benefit or cost resulting from cooperation should be fairly distributed between the participating players. Among these, the Shapley value [14] is both the most widely used point-valued solution concept in classical cooperative game theory applications and the predominant game theory based tool in XAI.

For n players, there are $2^n$ coalitions in a cooperative game. Calculating the Shapley value of a general cooperative game exactly was shown to be NP-hard [15,16,17]. Various concepts for approximating Shapley values for large n have been proposed — far too numerous to even attempt listing them all in this introduction — with plenty of recent approaches sparked by the desire to employ Shapley values for determining the significance of features with respect to the outcome of neural networks [10,11,12,13]. In our view, the large variety of approximation algorithms is at least partly due to the fact that the methods are tailored to specific formulations of the Shapley value, which can be expressed in several mathematically equivalent ways. For example, the Shapley value can be regarded as an average over permutations [18], as the optimal solution to a least squares problem [19], as a weighted sum over coalitions [14] or as an integral of the game’s multilinear extension [20].

To overcome this fragmentation, Benati et al. [21] proposed a unified framework as a “stochastic approach to approximate values in cooperative games”. In detail, their work defines a general approach for estimating linear solution concepts, a subclass of point-valued solution concepts including the Shapley value, via Monte Carlo methods. We extend their approach by introducing a novel connection to importance sampling, i.e., a variance reduction technique for Monte Carlo methods, thereby providing a stronger theoretical foundation for their framework as well as promoting it — since, as for our for our knowledge this article is the first work applying the ideas from [21] in XAI. Building on this connection, we show that several well-known approximation algorithms for the Shapley value, in particular some of those based on multilinear extensions, fit into this framework. This unifying perspective enables a more systematic comparison of algorithms, leading to new insights and recommendations for choosing among different approximation methods.

Importantly, the goal of this paper is neither to propose any completely novel approach for estimating Shapley values nor to perform comparisons over a range of approximation algorithms. Rather, beyond the importance sampling framework on the coalition space mentioned in the previous paragraph, we also provide structural and algorithmic insight into very widely used methods for computing Shapley values based on multilinear extensions. Specifically, we take a closer look at an algorithm originally called Owen sampling (OS) published by Okhrati and Lipani in 2021 [22], which, according to Google Scholar on Dec. 31, 2025, was already cited 76 times. The seminal paper Chen et al. [12] reports that OS is biased, but does not give a formal proof or perform a detailed numerical investigation into the bias. Further, Chen et al. [12] sketch an unbiased version for multilinear extension sampling while conceding that the biased approach may be more common as it improves convergence and later briefly mentioning in one sentence that their unbiased version was equivalent to sampling appropriately from the coalition space. This paper fills the theoretical and algorithmic gaps by

• establishing the bias in the original algorithm from [22] formally and confirming it in numerical experiments,
• formally specifying and analyzing the unbiased version for multilinear extension sampling sketched in [12] and pointing out its equivalence to a weighted sampling strategy on the coalition space rigorously using the aforementioned importance sampling framework and
• establishing and analyzing an unbiased and consistent estimator based on stratification, which outperforms the multilinear extension sampling algorithms from [12,22] on practical problems from cooperative game theory and XAI.

The paper is organized as follows. In Section 2, we introduce the basic ideas from cooperative game theory, including linear solution concepts and the Shapley value, whereas Section 3 provides a brief introduction to Monte Carlo methods as well as importance sampling and stratified sampling for variance reduction. We establish a novel framework for importance sampling on the coalition space for linear solution concepts based upon Benati et al. [21] in Section 4. The core of the paper follows in Section 5. We introduce and develop unbiased sampling algorithms based on multilinear extensions. Using these algorithms and the importance sampling framework from Section 4, we then derive new algorithms, which are closely related to the multilinear-extension-based algorithms, employing similar distributions and incorporating analogical architectures including sample reuse across players. We offer coalition-based sampling algorithms with more favourable properties than the multilinear-extension-based concepts from [12,22]. In Section 6, we analyze the performance of the algorithms from Section 5 empirically for two different types of cooperative games as well as three real-world XAI applications and confirm the claims from our previous analysis numerically. We end with a summary, our conclusions and some recommendations in Section 7.

2. Preliminaries on Cooperative Game Theory

This section provides a brief introduction to cooperative games with transferable utility and point-valued solution concepts such as the Shapley value along the lines of the two textbooks by Chakravarty et al. [1] and by Peters [23] as well as the research paper by Benati et al. [21].

2.1. Transferable Utility Games and Their Characteristic Functions

In this work, we study transferable utility games (TU games), which are cooperative games satisfying the restriction that the earnings of a coalition can be expressed by a scalar [1,23]. This comes with the implicit assumption that this value — the amount of utility — can be freely distributed and transferred among players.

Following the definitions in [1,23], a TU game is formally given by a pair $(N,v)$ with $N=\{1,\dots ,n\}$ being the set of players and $v:2^N\to \mathbb{R}$ being the characteristic function defined on the subsets of N. The characteristic function v assigns a real value to each coalition $S\subseteq N$ that represents its total utility when its members cooperate. We always follow the normalization $v(\varnothing )=0$.

It sometimes proves to be convenient to identify a coalition S with an indicator vector $\mathbf{z}\in {\{0,1\}}^n$. Thus, we define $\mathbf{z}\left(S\right)={\left[{\mathbb{1}}_{i\in S}\right]}_{i=1}^n$, and, conversely, $S\left(\mathbf{z}\right)=\{i\in N\mid z_i=1\}$. Under this mapping, all functions and operators defined on either $2^N$ or ${\{0,1\}}^n$ are implicitly also defined on its counterpart, e.g., $v\left(S\right)$ and $v\left(\mathbf{z}\right)$ are used interchangeably and $\left|\mathbf{z}\right|$ is equivalent to $\left|S\right(\mathbf{z}\left)\right|$.

2.2. Linear Solution Concepts for TU Games

Let $G\left(N\right)$ be the set of TU games with player set N. A point-valued solution concept is a map $\alpha :G\left(N\right)\to {\mathbb{R}}^n$ that assigns a vector $\alpha (N,v)\in {\mathbb{R}}^n$ to each game $(N,v)\in G\left(N\right)$. Each element ${\alpha }_i$ of that vector $\alpha$ denotes the worth of player i according to the underlying solution concept. In general, point-valued solution concepts are used when the influences of players need to be measured, or when the payoff $v\left(S\right)$ of a coalition of cooperating players $S\subseteq N$ needs to be allocated between the members of S.

A desirable property of a solution concept $\alpha$ is linearity [21], which means that $\alpha$ can be written as

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

(1)

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

2.3. The Shapley Value

The Shapley value [14] is the most prominent solution concept for TU games. Recently, it received plenty of attention in the context of machine learning and explainable artificial intelligence [10,11,12]. For each player $i\in N$, its Shapley value is defined as the expected marginal contribution ${\Delta }_i(S,v)=v(S\cup \left\{i\right\})-v\left(S\right)$ to the set S of all players before player i in a random permutation $\pi$ of N, i.e.,

$$\begin{array}{cc}\hfill {\varphi }_i(N,v)& ={\mathbb{E}}_{\pi \sim \mathcal{U}\left(\Pi \right(N\left)\right)}\left[{\Delta }_i\left({Pre}_i\left(\pi \right),v\right)\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{n!}}\sum _{\pi \in \Pi \left(N\right)}{\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!}}{\Delta }_i(S,v),\hfill \end{array}$$

(3)

where $\mathcal{U}$ denotes a uniform distribution, $\Pi \left(N\right)$ is the set of all permutations of N, and ${Pre}_i\left(\pi \right)$ returns the set of all players before player i in a given permutation $\pi \in \Pi \left(N\right)$. Note that we will write ${\Delta }_i\left(S\right)={\Delta }_i(S,v)$ as long as the characteristic function v is clear from context.

The vector of all players’ Shapley values is given by $\varphi ={[{\varphi }_1,\dots ,{\varphi }_n]}^{\top }$, and any approximation thereof is denoted by $\widehat{\varphi }$, with its i-th component given by ${\widehat{\varphi }}_i$. As long as there is no ambiguity, we will write ${\varphi }_i={\varphi }_i(N,v)$.

Clearly, via (3), one obtains that the Shapley value fits into the framework of linear solution concepts defined in Section 2.2, with the individual elements of $\mathit{w}\left(S\right)$ being defined as

$$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)

and ${\nu }_i(S,v)={\Delta }_i(S,v)$.

2.4. The Banzhaf Value

The Banzhaf value [24,25] is another prominent point-valued solution concept for TU games. In the fields of machine learning and artificial intelligence, it can be regarded the second most important solution concept from cooperative game theory [26,27]. For a player $i\in N$, the Banzhaf value is defined as the expected marginal contribution of i to a random coalition, i.e.,

$${\beta }_i(N,v)={\mathbb{E}}_{S\sim \mathcal{U}\left(2^{N\setminus \left\{i\right\}}\right)}\left[{\Delta }_i(S,v)\right]={\displaystyle \frac{1}{2^{n-1}}}\sum _{\begin{array}{c}S\subseteq N\\ i\notin S\end{array}}{\Delta }_i(S,v),$$

(5)

with $\mathcal{U}$ denoting the uniform distribution and ${\Delta }_i(S,v)$ being the marginal contribution of i to S, as in the previous subsection. The vector of all players’ Banzhaf values is given by $\beta ={[{\beta }_1,\dots ,{\beta }_n]}^{\top }$. Similar to the Shapley value, the Banzhaf value also fits into the framework of linear solution concepts from Section 2.2, with the individual components of $\mathit{w}\left(S\right)$ being defined as

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

and ${\nu }_i(S,v)={\Delta }_i(S,v)$.

3. Preliminaries on Monte Carlo Methods for Estimation

In this section, we give a brief introduction to the application of Monte Carlo methods for estimating expectations. We mostly follow the two short overviews provided by Rubinstein and Kroese [28] and by Botev and Ridder [29].

Let $\mathbf{X}$ be a discrete (or continuous) d-dimensional random variable taking values from $\mathcal{X}\subseteq {\mathbb{R}}^d$ and following a known or unknown probability mass function $q_{\mathbf{X}}$ (or probability density function, respectively). Furthermore, let $\mathbf{x}\in \mathcal{X}$ represent a concrete realization of $\mathbf{X}$. For a given function $H:\mathcal{X}\to \mathbb{R}$, our goal is to estimate

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

(6)

For many practical problems, this expectation cannot be computed in closed form, because often $q_{\mathbf{X}}$ itself is unknown (we only have i.i.d. samples, e.g., in the form of observations or measurements), or because the sum (or integral, respectively) defining the expectation is intractable to compute. In the latter case, Monte Carlo methods replace the analytical challenge by simulation.

We start by introducing the crude Monte Carlo method, which is the simplest approach, followed by two variance reduction techniques studied in this paper. We emphasize that we always assume finite variances of $H\left(\mathbf{X}\right)$ throughout this article.

3.1. Crude Monte Carlo Method

When using the Crude Monte Carlo method, a sample ${\mathbf{X}}_1,\dots ,{\mathbf{X}}_{\tau }{\sim }_{\mathrm{iid}}{q}_{\mathbf{X}}$ of size $\tau$ is drawn, and $\mu$ is approximated by averaging all $H\left({\mathbf{X}}_k\right)$. Thus, Botev and Ridder [29] define an estimator of $\mu$ as

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

(7)

This estimator is unbiased [29], i.e.,

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

and its variance, as stated in [29], is given by

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

(8)

Furthermore, one can easily obtain via Chebyshev’s inequality [31] that $\widehat{\mu }$ is consistent in probability, i.e., for all $\epsilon >0$, there holds

$$\mathbb{P}\left(|\widehat{\mu }-\mu |\ge \epsilon \right)\le {\displaystyle \frac{Var\left[\widehat{\mu }\right]}{{\epsilon }^2}}={\displaystyle \frac{Var\left[H\right(\mathbf{X}\left)\right]}{\tau {\epsilon }^2}}\stackrel{\tau \to \infty }{\to }0.$$

3.2. Importance Sampling for Variance Reduction

Following the introductions of importance sampling in [28,30], we provide a brief unified treatment that covers both discrete and continuous probabilities.

In the context of importance sampling, another probability mass function (or probability density function, respectively) $p_{\mathbf{X}}$ with $p_{\mathbf{X}}\left(\mathbf{x}\right)=0\Rightarrow H\left(\mathbf{x}\right)q_{\mathbf{X}}\left(\mathbf{x}\right)=0$ is introduced, such that (6) can 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].$$

(9)

Based on (9), a Monte Carlo approximation is given by

$$\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),$$

(10)

with ${\mathbf{X}}_1,\dots ,{\mathbf{X}}_{\tau }{\sim }_{\mathrm{iid}}{p}_{\mathbf{X}}$. The resulting estimator is unbiased [28], i.e.,

$${\mathbb{E}}_{p_{\mathbf{X}}}\left[\widehat{\mu }\right]=\mu ,$$

(11)

and clearly, as for the crude Monte Carlo method, its variance is given by

$${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].$$

(12)

As long as $p_{\mathbf{X}}$ is properly selected, the variance of the estimator is always less or equal compared to the variance of the crude Monte Carlo estimator from Section 3.1. This can easily be seen by comparing (8) and (12), resulting in

$${\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 always holds since one could always choose $p_{\mathbf{X}}=q_{\mathbf{X}}$ in the worst case.

The importance sampling estimator is also consistent in probability, which can be confirmed via Chebyshev’s inequality [31]. For all $\epsilon >0$, there holds

$$\mathbb{P}\left(|\widehat{\mu }-\mu |\ge \epsilon \right)\le {\displaystyle \frac{Var\left[\widehat{\mu }\right]}{{\epsilon }^2}}={\displaystyle \frac{{Var}_{p_{\mathbf{X}}}\left[\frac{q_{\mathbf{X}}\left(\mathbf{X}\right)}{p_{\mathbf{X}}\left(\mathbf{X}\right)}H\left(\mathbf{X}\right)\right]}{\tau {\epsilon }^2}}\stackrel{\tau \to \infty }{\to }0.$$

(13)

3.3. Stratified Sampling for Variance Reduction

Stratified sampling [29] is a well-known variance reduction technique for Monte Carlo methods. It partitions 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 }$. Let L be a discrete random variable taking values from $\{1,\dots ,\ell \}$ with known probabilities $p_L\left(l\right)=\mathbb{P}(L=l)=\mathbb{P}(\mathbf{X}\in {\mathcal{X}}_l)$. Then, $\mu$ can be rewritten 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]$ being the expectation over the conditional probability distribution of $\mathbf{X}$ given that $\mathbf{X}\in {\mathcal{X}}_l$. Based on this formulation of $\mu$, one obtains an estimator of $\mu$ as

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

where ${\widehat{\mu }}_l$ is 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),$$

(14)

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

$$\mathbb{E}\left[\widehat{\mu }\right]=\mu ,$$

and its variance is discussed in the following. Unlike in the previous subsections, our goal here is not to simply define the variance of the estimator, but instead to derive it in order to understand its origins. This deeper understanding will be crucial for our subsequent analysis in Section 5.

For the two random variables $H\left(\mathbf{X}\right)$ and L defined earlier, Eve’s law [31], also referred to as the law of total variance or the variance decomposition formula in the literature, states

$$Var\left[H\left(\mathbf{X}\right)\right]={\mathbb{E}}_{p_L}\left[{Var}_{q_{\mathbf{X}}}[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_L]\right]+{Var}_{p_L}\left[{\mathbb{E}}_{q_{\mathbf{X}}}[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_L]\right],$$

(15)

and therefore, by using (8), the variance of the crude Monte Carlo estimator from Section 3.1 is given by

$${\displaystyle \frac{1}{\tau }}{\mathbb{E}}_{p_L}\left[{Var}_{q_{\mathbf{X}}}[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_L]\right]+{\displaystyle \frac{1}{\tau }}{Var}_{p_L}\left[{\mathbb{E}}_{q_{\mathbf{X}}}[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_L]\right].$$

(16)

Since the individual estimators ${\widehat{\mu }}_l$ are independent, the variance of the estimator based on stratified sampling is given by

$$Var\left[\widehat{\mu }\right]=Var\left[\sum _{l=1}^{\ell }{p}_L\left(l\right){\widehat{\mu }}_l\right]=\sum _{l=1}^{\ell }{p}_L{\left(l\right)}^{2}Var\left[{\widehat{\mu }}_l\right]=\sum _{l=1}^{\ell }{p}_L{\left(l\right)}^2{\displaystyle \frac{Var[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_l]}{{\tau }_l}}.$$

(17)

Once we assume that the sample sizes per stratum are proportionally assigned, i.e., ${\tau }_l=p_L\left(l\right)\tau$, Equation (17) simplifies to

$$\begin{array}{cc}\hfill Var\left[\widehat{\mu }\right]& =\sum _{l=1}^{\ell }{p}_L{\left(l\right)}^2{\displaystyle \frac{Var[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_l]}{p_L\left(l\right)\tau }}={\displaystyle \frac{1}{\tau }}\sum _{l=1}^{\ell }{p}_L\left(l\right)Var[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_l]\hfill \\ \hfill & ={\displaystyle \frac{1}{\tau }}{\mathbb{E}}_{p_L}\left[{Var}_{q_{\mathbf{X}}}[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_L]\right].\hfill \end{array}$$

(18)

Therefore, by comparing the variance from (18) to the variance of the crude Monte Carlo estimator given in (16), it is easy to obtain that the variance of the stratified estimator is always less or equal compared to that of the crude Monte Carlo estimator, as long as the stratum sample sizes are proportionally assigned with respect to $p_L$, see also [29]. In detail, we highlight that stratification removes the variance of the expectations between strata, i.e., the non-negative term ${\displaystyle \frac{1}{\tau }}{Var}_{p_L}\left[{\mathbb{E}}_{q_{\mathbf{X}}}[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_L]\right]$ in (16) vanishes.

For completeness, we note that the strata sample sizes ${\tau }_l$ which minimize the variance of the overall estimator [29] are given by

$${\tau }_l={\displaystyle \frac{p_L\left(l\right){\sigma }_l}{{\sum }_{l^{\prime }=1}^{\ell }{p}_L\left(l^{\prime }\right){\sigma }_{l^{\prime }}}}\tau \mathrm{with}{\sigma }_l=\sqrt{Var[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_l]}.$$

(19)

Since in most settings the standard deviations ${\sigma }_l$ are not known, one might use two-stage algorithms that use pilot runs in the first stage to approximate all ${\sigma }_l$ and calculate the final sample sizes via (19) based on those ${\sigma }_l$.

Concluding this subsection, we note that the estimator obtained via stratified sampling with a sample allocation proportional to $p_L$ is consistent in probability, which can again be shown via Chebyshev’s inequality [31]. For all $\epsilon >0$, there holds

$$\mathbb{P}\left(|\widehat{\mu }-\mu |\ge \epsilon \right)\le {\displaystyle \frac{Var\left[\widehat{\mu }\right]}{{\epsilon }^2}}={\displaystyle \frac{{\mathbb{E}}_{p_L}\left[{Var}_{q_{\mathbf{X}}}[H\left(\mathbf{X}\right)\mid \mathbf{X}\in {\mathcal{X}}_L]\right]}{\tau {\epsilon }^2}}\stackrel{\tau \to \infty }{\to }0.$$

4. An Importance Sampling Framework for Linear Solution Concepts and Marginal Contribution Importance Sampling for the Shapley Value

Calculating the value of a linear solution concept, i.e., (1), exactly usually requires summing up $2^n$ or $2^{n-1}$ terms. Thus, approximation approaches are needed to estimate these values in real-world scenarios. The latter statement is particularly relevant in the context of applications in explainable artificial intelligence [10,11,12,13] where one can frequently not exploit any special structure of the underlying game in order to compute Shapley values exactly. In the following two subsections, we revisit some ideas from Benati et al. [21] and present them in a general framework for importance sampling on the coalition space for linear solution concepts. We first introduce the crude Monte Carlo method on the coalition space based on the uniform sampling strategy from [21]. Afterwards, we develop our general importance sampling framework for approximating linear solution concepts and point out that the weighted sampling strategy for estimating the Shapley value from (3) introduced in [21] can be regarded as marginal contribution importance sampling.

4.1. Crude Monte Carlo Method on the Coalition Space

Benati et al. [21] discuss a uniform sampling strategy on the coalition space. We will in this work consider it as the crude Monte Carlo method on the coalition space.

First, by reformulating (1) for a specific player $i\in N$, we obtain

$${\alpha }_i(N,v)=\sum _{S\subseteq N}{w}_i\left(S\right){\nu }_i(S,v)=\sum _{S\subseteq N}{\displaystyle \frac{1}{2^n}}{2}^n{w}_i\left(S\right){\nu }_i(S,v)={\mathbb{E}}_{S\sim q}\left[2^n{w}_i\left(S\right){\nu }_i(S,v)\right]$$

(20)

with $q=\mathcal{U}\left(2^N\right)$ being a uniform distribution over all subsets, i.e., $q\left(S\right)=2^{-n}$.

Thus, a simple approximation of (1) via the crude Monte Carlo method in a sense of (7) is given by generating a sample $\mathbb{S}{\sim }_{\mathrm{iid}}q$ of size $\tau$ with $\mathbb{S}={[{\mathcal{S}}_1,\dots ,{\mathcal{S}}_{\tau }]}^{\top }$, and, based on that, calculating

$$\widehat{\alpha }(N,v)={\left[{\widehat{\alpha }}_i(N,v)\right]}_{i\in N}={\left[{\displaystyle \frac{2^n}{\tau }}\sum _{\mathcal{S}\in \mathbb{S}}{w}_i\left(\mathcal{S}\right){\nu }_i(\mathcal{S},v)\right]}_{i\in N}.$$

Note that throughout this work, we use S to denote subsets in general, while sampled subsets are denoted by $\mathcal{S}$.

We emphasize that using the uniform distribution over all subsets, q, is the most natural and simplest choice for the crude Monte Carlo approximation of (1), since more sophisticated probability distributions are not justified in this setting. Now, one might ask why the weights $w_i$ cannot be interpreted as probabilities of different coalitions forming, and therefore, serve as the initial underlying distribution that defines the distribution for the crude Monte Carlo method. For instance, for the well-known formulation of the Shapley value given by (3), which follows directly from the assumption that all player permutations are equally likely as defined in (2), these weights indeed can be interpreted as the probability of a subset $S\subseteq N\setminus \left\{i\right\}$ forming. However, this is not the case for all formulations of the Shapley value, and certainly not for all linear solution concepts [21]. Thus, in general, we cannot use $w_i$ as q in (20).

With that in mind, one might argue that any valid probability distribution (other than the uniform distribution over subsets) could be used in (20), with the inner expectation adapted accordingly, which would also be mathematically correct. However, TU games as introduced in Section 2.1 are specified only by the tuple $(N,v)$ and do not include any probabilities for the formation of any specific coalitions. (Note in passing that different solution concepts may account for the latter. For example, the Banzhaf value introduced in Section 2.4 assumes that each coalition is equally likely, while the Shapley value introduced in Section 2.3 assumes that each permutation is equally likely, which implies a non-uniform distribution over coalitions.) Therefore, since there is no basis for assigning specific probabilities to different coalitions, we adopt the uniform distribution $q=\mathcal{U}\left(2^N\right)$ as a natural and straightforward choice for our crude Monte Carlo method, providing a reasonable starting point for all linear solution concepts.

Readers familiar with the ApproShapley algorithm [18] might think of ApproShapley as the crude Monte Carlo method for approximating the Shapley value contradicting our results above. However, the ApproShapley algorithm addresses a different formulation, i.e., the permutation-based definition of the Shapley value given by (2), where the sampling space is given by all permutations of N. In contrast, when the Shapley value is expressed within the coalition-based framework of linear solution concepts as in (1), and the only requirement is that $\mathit{w}\in {\mathbb{R}}^n$, the appropriate crude Monte Carlo approximation — based on our previous arguments — is given by sampling according to the uniform distribution over all coalitions and approximating $2^n{w}_i\left(S\right){\nu }_i(S,v)$.

4.2. Importance Sampling on the Coalition Space

In the previous subsection, we focused exclusively on sampling subsets from the uniform distribution. While this approach is simple and unbiased, it is not necessarily optimal for all linear solution concepts or problem settings. For example, when computing the Shapley value using (3), small or large coalitions tend to have larger weights. Thus, despite having a strong influence on the estimator, they may retrieve only a few samples when sampling uniformly from the coalition space, potentially reducing the accuracy of the estimator.

To address this issue and enable more efficient sampling schemes, we now consider the use of importance sampling introduced in Section 3.2 as a variance reduction technique. By appropriately reweighting the estimator, importance sampling allows drawing subsets from a non-uniform, user-defined distribution, while still providing an unbiased estimator of the underlying linear solution concept.

Theorem 1. 

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}$$

(21)

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

Proof. 

In the following, without loss of generality, we fix one player $i\in N$ to simplify notations.

First, recall the formulation of linear solution concepts via the uniform distribution over coalitions denoted by q, as defined in (20). Now, we define

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

which simplifies to $p_i\left(S\right)=0\Rightarrow w_i\left(S\right){\nu }_i\left(S\right)=0$. This serves as our new sampling distribution, and it clearly satisfies the requirements for a sampling distribution in the context of importance sampling as specified in Section 3.2. Then, based on (9), we derive

$${\alpha }_i(N,v)={\mathbb{E}}_{S\sim q}\left[2^n{w}_i\left(S\right){\nu }_i(S,v)\right]={\mathbb{E}}_{S\sim p_i}\left[{\displaystyle \frac{q\left(S\right)}{p_i\left(S\right)}}{2}^n{w}_i\left(S\right){\nu }_i(S,v)\right]={\mathbb{E}}_{S\sim p_i}\left[{\displaystyle \frac{w_i\left(S\right)}{p_i\left(S\right)}}{\nu }_i(S,v)\right].$$

Now, via (10), one obtains

$${\widehat{\alpha }}_i(N,v)={\displaystyle \frac{1}{\tau }}\sum _{\mathcal{S}\in {\mathbb{S}}_i}{\displaystyle \frac{q\left(\mathcal{S}\right)}{p_i\left(\mathcal{S}\right)}}{2}^n{w}_i\left(\mathcal{S}\right){\nu }_i(\mathcal{S},v)={\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)$$

as the importance sampling estimator of ${\alpha }_i(N,v)$ as long as ${\mathbb{S}}_i{\sim }_{\mathrm{iid}}{p}_i$.    □

The following proposition relates our importance sampling estimator from Theorem 1 to a result from Benati et al. [21], 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).$$

(22)

(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).$$

(23)

(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).$$

(a)

Proof. 

The proof is straightforward. (22) directly follows from (11), while (24) can be obtained from (13).

By using (12), the variance of the estimator of player i is given by

$$\begin{array}{cc}\hfill Var\left[{\widehat{\alpha }}_i(N,v)\right]& ={\displaystyle \frac{1}{\tau }}{Var}_{S\sim p_i}\left[{\displaystyle \frac{q\left(S\right)}{p_i\left(S\right)}}{2}^n{w}_i\left(S\right){\nu }_i(S,v)\right]={\displaystyle \frac{1}{\tau }}{Var}_{S\sim p_i}\left[{\displaystyle \frac{w_i\left(S\right)}{p_i\left(S\right)}}{\nu }_i(S,v)\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{\tau }}\left({\mathbb{E}}_{S\sim p_i}\left[{\displaystyle \frac{w_i{\left(S\right)}^2}{p_i{\left(S\right)}^2}}{\nu }_i{(S,v)}^2\right]-{\mathbb{E}}_{S\sim p_i}{\left[{\displaystyle \frac{w_i\left(S\right)}{p_i\left(S\right)}}{\nu }_i(S,v)\right]}^2\right)\hfill \\ \hfill & ={\displaystyle \frac{1}{\tau }}\left(\sum _{S\subseteq N}{p}_i\left(S\right){\displaystyle \frac{w_i{\left(S\right)}^2}{p_i{\left(S\right)}^2}}{\nu }_i{(S,v)}^2-{\alpha }_i{(N,v)}^2\right)\hfill \\ \hfill & ={\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).\hfill \end{array}$$

   □

Remark 1. 

It is worth noting that the results from Theorem 1 and Proposition 1 implicitly cover the crude Monte Carlo method from Section 3.1 (and Section 4.1, respectively) by setting $p_i\left(S\right)=q\left(S\right)=2^{-n}$. Additionally, it also covers stratified sampling from Section 3.3. Clearly, each individual stratum estimator of a linear solution concept can be brought into the form of (21) and treated as an estimator of an individual linear solution concept, which allows for applying the results from Theorem 1 and Proposition 1.

We note that the goal of our results in this subsection is not to perform better than algorithms tailored to specific solution concepts. Instead, our results provide a general framework that is theoretically grounded in importance sampling on the coalition space, which allows for employing any suitable sampling distribution to approximate any linear solution concept from cooperative game theory. Furthermore, as stated in Remark 1, our approach implicitly covers the crude Monte Carlo method as well as stratified sampling, making it even more general.

Let us finally define a special class of importance sampling algorithms, which we call marginal contribution importance sampling algorithms. As the name suggests, these algorithms approximate linear solution concepts by sampling marginal contributions via importance sampling. Thus, all algorithms that admit the form proposed in Theorem 1 with ${\nu }_i(S,v)={\Delta }_i(S,v)=v(S\cup \left\{i\right\})-v\left(S\right)$ belong to this class.

A simple estimator named MCIS belonging to the class of marginal contribution importance sampling algorithms for approximating the Shapley value of player i is given by

$${\widehat{\varphi }}_i={\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)}}{\Delta }_i\left(S\right),$$

(25)

with $w_i$ given by (4) and ${\mathbb{S}}_i{\sim }_{\mathrm{iid}}{p}_i$, whereby

$$p_i\left(S\right)=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.$$

(26)

is the sampling distribution. Note that we can use $w_i$ as the sampling distribution since $w_i:2^N\to [0,1]$ and ${\sum }_{S\subseteq N}{w}_i\left(S\right)=1$. Thus,

$${\widehat{\varphi }}_i={\displaystyle \frac{1}{\tau }}\sum _{\mathcal{S}\in {\mathbb{S}}_i}{\Delta }_i\left(\mathcal{S}\right).$$

(27)

We will revisit and further analyze the MCIS estimator for the Shapley value in Section 5.6. We finally note that MCIS estimator for the Shapley value derived via (25), (26) and (27) was discussed in Benati et al. [21], p.96, under the name weighted sampling strategy. We end this section by stressing that the least squares strategy for Shapley values which is promoted and analyzed in detail in the article by Benati et al. [21], is covered by Theorem 1 and Proposition 1. While any approximations based on the least squares formula for the Shapley value are beyond the scope of this manuscript, our key takeaways are both the generality of our importance sampling framework and the relevance of the specific representation of the Shapley value for its Monte Carlo estimation.

5. Multilinear Extensions for Approximating Shapley Values and Connections to Importance Sampling on the Coalition Space

This section covers the idea of multilinear extensions of cooperative games and how this formulation of games can be used to approximate Shapley values of cooperative games. First, we provide a definition of multilinear extensions in Section 5.1 and a representation of the Shapley value in their context in Section 5.2. Afterwards, in Section 5.3, we discuss an existing multilinear-extension-based approximation algorithm for the Shapley value which goes by the name Owen Sampling, was introduced in Okhrati and Lipani [22] and is frequently applied in XAI applictions. Although it was previously noticed that Owen Sampling is biased, e.g. in [12], we are the first to provide a detailed analysis of the bias. Then, we propose unbiased multilinear-extension-based algorithms in Section 5.4 and  Section 5.5, with the latter subsection offering a stratified variant. These two unbiased multilinear-extension-based approaches serve as intermediary results for our investigations on the relations to our importance sampling framework on the coalition space from Section 4.2. We discuss our results, conclusions and recommendations in the subsequent Section 5.6 and  Section 5.7.

5.1. Definition of Multilinear Extensions

Multilinear extensions have been proposed by Owen [20]. The idea is to extend the characteristic function of a cooperative game, v, which is defined on the subsets of N, to the unit n-hypercube by the function f. As mentioned in Section 2.1, v is implicitly defined on ${\{0,1\}}^n$. Based on that formulation, Owen [20] extended the domain to ${[0,1]}^n$ and defined the extension of v as

$$f\left(\mathbf{x}\right)=\sum _{S\subseteq N}\prod _{i\in S}{x}_i\prod _{i\notin S}(1-x_i)v\left(S\right).$$

(28)

It is easy to see that this function coincidences with v at the corners of the hypercube, i.e., $f\left(\mathbf{x}\right)=v\left(\mathbf{x}\right)$ for all $\mathbf{x}\in {\{0,1\}}^n\subset {[0,1]}^n$, which means that f is indeed an extension of v. Furthermore, f is linear in each variable, i.e., multilinear, hence the name multilinear extension. We note that Owen [20] showed that f as defined in (28) is the unique multilinear extension of v.

Furthermore, Owen [20] highlighted that f can be thought of as the expected value of v when each player i joins the coalition with their respective probability $x_i$, i.e.,

$$f\left(\mathbf{x}\right)={\mathbb{E}}_{\mathbf{z}\sim \mathcal{B}\left(\mathbf{x}\right)}\left[v\left(\mathbf{z}\right)\right],$$

whereby $\mathcal{B}\left(\mathbf{x}\right)$ is the Bernoulli distribution with probabilities $\mathbf{x}\in {\mathbb{R}}^n$. Note that all elements of $\mathbf{z}$ are mutually independent events and are obtained via $z_i\sim \mathcal{B}\left(x_i\right)$ for all $i\in N$.

5.2. The Shapley Value in the Context of Multilinear Extensions

In the setting of multilinear extensions, Owen [20] provided a definition of the Shapley value as

$${\varphi }_i={\int }_0^1{\displaystyle \frac{\partial f}{\partial x_i}}\left(t\mathbf{1}\right)dt,$$

(29)

with $\mathbf{1}\in {\mathbb{R}}^n$ being the all-ones vector and

$${\displaystyle \frac{\partial f}{\partial x_i}}\left(\mathbf{x}\right)=\sum _{\begin{array}{c}S\subseteq N\\ i\notin S\end{array}}\prod _{j\in S}{x}_j\prod _{\begin{array}{c}j\notin S\\ j\ne i\end{array}}(1-x_j)\left(v(S\cup \{i\left\}\right)-v\left(S\right)\right).$$

(30)

This formulation can also be interpreted as the expected marginal contribution of player i when all other players $j\in N\setminus \left\{i\right\}$ join the coalition with probability $x_j$, i.e.,

$${\displaystyle \frac{\partial f}{\partial x_i}}\left(\mathbf{x}\right)={\mathbb{E}}_{\mathbf{z}\sim \mathcal{B}(\mathbf{x}\odot {\overline{\mathbf{e}}}_i)}[v(\mathbf{z}+{\mathbf{e}}_i)-v\left(\mathbf{z}\right)],$$

(31)

where ⊙ denotes the Hadamard product of vectors, and ${\overline{\mathbf{e}}}_i=\mathbf{1}-{\mathbf{e}}_i$ with ${\mathbf{e}}_i\in {\mathbb{R}}^n$ being the i-th standard basis vector. Note that $\mathbf{z}\sim \mathcal{B}(\mathbf{x}\odot {\overline{\mathbf{e}}}_i)$ means sampling each $z_j\sim \mathcal{B}\left(x_j\right)$ except the i-th component, which is forced to be 0, such that player i is explicitly excluded.

The process of calculating the Shapley value via (29) and (30) (or (31), respectively) is visualized in Figure 1.

5.3. Owen Sampling (OS) for Computing Shapley Values According to Okhrati and Lipani

Since calculating (30) or (31) exactly requires $2^{n-1}$ evaluations of v, multilinear extensions cannot be used as an exact computation method for large n when no a priori knowledge about the characteristic function is given, but instead serve as a starting point for the development of new approximation algorithms.

In this section, we introduce the OS algorithm proposed by Okhrati and Lipani [22]. To avoid confusion, we highlight that the name of this algorithm is based on Guillermo Owen, who published the idea of multilinear extensions of games [20,25] and defined the Shapley value in this setting [20]. In particular, we emphasize that the algorithm is not related to the Owen value [32], a solution concept for cooperative games with a priori unions.

Okhrati and Lipani [22] approximate (29) by

$${\displaystyle \frac{1}{Q}}\sum _{q=0}^{Q-1}{\displaystyle \frac{\partial f}{\partial x_i}}\left({\displaystyle \frac{q}{Q-1}}\mathbf{1}\right)$$

(32)

with $Q\in \mathbb{N}$ being the discretization level. Note in passing that, contrary to the original algorithm in [22], we adjusted the upper bound of the sum in (32) from Q to $Q-1$ to ensure that the sum includes exactly Q terms, rather than $Q+1$. Our adjustment is consistent with the source code provided with the original paper [22] which is publicly available via

https://github.com/aldolipani/OwenShap.

Since calculating $\partial f/\partial x_i$ exactly in (32) requires summing over all subsets $S\subseteq N\setminus \left\{i\right\}$, the authors approximate (31) by

$${\displaystyle \frac{1}{{\tau }_q}}\sum _{k=1}^{{\tau }_q}\left(v({\mathbf{z}}_{q,i,k}+{\mathbf{e}}_i)-v\left({\mathbf{z}}_{q,i,k}\right)\right),$$

(33)

whereby ${\mathbf{z}}_{q,i,k}$ denotes the k-th sampled vector that is drawn according to the distribution from (31) with all probabilities $x_j$ for $j\in N\setminus \left\{i\right\}$ being $q/(Q-1)$, i.e., ${\mathbf{z}}_{q,i,k}\sim \mathcal{B}\left({\overline{\mathbf{e}}}_{i}q/(Q-1)\right)$. ${\tau }_q$ hereby is an additional parameter of the algorithm that specifies the number of evaluations per q. Okhrati and Lipani [22] fixed ${\tau }_q=2$ for all $q\in \{0,\dots ,Q-1\}$.

Thus, by combining (32) and (33), the estimator of the Shapley value for player i is given by

$${\widehat{\varphi }}_i={\displaystyle \frac{1}{Q{\tau }_q}}\sum _{q=0}^{Q-1}\sum _{k=1}^{{\tau }_q}\left(v({\mathbf{z}}_{q,i,k}+{\mathbf{e}}_i)-v\left({\mathbf{z}}_{q,i,k}\right)\right).$$

(34)

It has been shown in the original paper Okhrati and Lipani [22] that this estimator is consistent, i.e., for any $i\in N$ and $\epsilon >0:{lim}_{Q,M\to \infty }\mathbb{P}\left(\right|{\widehat{\varphi }}_i-{\varphi }_i|>\epsilon )=0$. Unfortunately, Okhrati and Lipani [22] did not include a theoretical investigation of bias or variance. The bias of OS is discussed in [12], where the authors concluded that OS is biased without giving a formal proof. We provide the latter in

Proposition 2. 

The OS estimator is biased, i.e., $\mathbb{E}\left[{\widehat{\varphi }}_i\right]={\varphi }_i$ does not hold in the general case.

Proof. 

Let us fix an arbitrary player $i\in N$ whose Shapley value we want to estimate. We state that (32) does not provide an unbiased approximation of (29), resulting in a biased estimator ${\widehat{\varphi }}_i$. To see this, consider the case $Q=3$, where samples are drawn from $\mathcal{B}\left(\mathbf{0}\right)$, $\mathcal{B}\left({\displaystyle \frac{1}{2}}{\overline{\mathbf{e}}}_i\right)$, and $\mathcal{B}\left({\overline{\mathbf{e}}}_i\right)$, such that we have

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\widehat{\varphi }}_i\right]& =\mathbb{E}\left[{\displaystyle \frac{1}{Q{\tau }_q}}\sum _{q=0}^{Q-1}\sum _{k=1}^{{\tau }_q}\left(v({\mathbf{z}}_{q,i,k}+{\mathbf{e}}_i)-v\left({\mathbf{z}}_{q,i,k}\right)\right)\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{Q}}\sum _{q=0}^{Q-1}{\displaystyle \frac{1}{{\tau }_q}}\sum _{k=1}^{{\tau }_q}{\mathbb{E}}_{\mathbf{z}\sim \mathcal{B}({\overline{\mathbf{e}}}_{i}q/(Q-1))}\left[v(\mathbf{z}+{\mathbf{e}}_i)-v\left(\mathbf{z}\right)\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{Q}}\sum _{q=0}^{Q-1}{\mathbb{E}}_{\mathbf{z}\sim \mathcal{B}({\overline{\mathbf{e}}}_{i}q/(Q-1))}\left[v(\mathbf{z}+{\mathbf{e}}_i)-v\left(\mathbf{z}\right)\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{3}}\left(v\left({\mathbf{e}}_i\right)-v\left(\mathbf{0}\right)+v\left(\mathbf{1}\right)-v\left({\overline{\mathbf{e}}}_i\right)+{\mathbb{E}}_{\mathbf{z}\sim \mathcal{B}\left({\displaystyle \frac{1}{2}}{\overline{\mathbf{e}}}_i\right)}\left[v(\mathbf{z}+{\mathbf{e}}_i)-v\left(\mathbf{z}\right)\right]\right)\hfill \\ \hfill & ={\displaystyle \frac{1}{3}}\left(v\left(\left\{i\right\}\right)+v\left(N\right)-v(N\setminus \left\{i\right\})+{\beta }_i\right).\hfill \end{array}$$

Clearly, the last equality does not represent a valid relationship between the Shapley value and the Banzhaf value. E.g., consider a game with $n=3$ players and the characteristic function

$$v(\varnothing )=0,v\left(\left\{j\right\}\right)=0(\forall j\in N),v\left(\{j,j^{\prime }\}\right)=1(\forall j\ne j^{\prime }\in N),v\left(N\right)=1,$$

such that the Shapley values are given by

$${\varphi }_j={\displaystyle \frac{1}{3}}(\forall j\in N).$$

Then, we obtain

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\widehat{\varphi }}_i\right]& ={\displaystyle \frac{1}{3}}\left(v\left(\right\{i\left\}\right)+v\left(N\right)-v(N\setminus \{i\left\}\right)\right)+{\displaystyle \frac{1}{3}}{\displaystyle \frac{1}{2^2}}\sum _{S\subseteq N\setminus \left\{i\right\}}\left(v(S\cup \{i\left\}\right)-v\left(S\right)\right)\hfill \\ \hfill & ={\displaystyle \frac{1}{12}}\sum _{S\subseteq N\setminus \left\{i\right\}}\left(v(S\cup \{i\left\}\right)-v\left(S\right)\right)={\displaystyle \frac{1}{6}}\ne {\displaystyle \frac{1}{3}}={\varphi }_i.\hfill \end{array}$$

   □

Proposition 2 establishes that the OS estimator is biased. (Formally, one might argue that Proposition 2 establishes the existence of a bias for the case $Q=3$ by counterexample. Although our proof is specific to this setting, the form of the OS estimator indicates that the bias persists for any finite Q, see also the numerical results in [12]). The origin of that bias can be observed by visualizing the expected number of sampled coalitions per coalition size across different values of Q, as illustrated in Figure 2. Clearly, as Q increases, the expected number of samples per coalition size approaches an equal distribution, as required by the Shapley definition. However, for the two finite values of Q shown in Figure 2, this equal distribution is never fully reached.

Additionally, we refer to Figure 3 in Section 6, which empirically validates the existence of a bias when executing OS.

Due to its biasedness, OS does not fit into the importance sampling framework defined by Theorem 1 since estimators of that framework are unbiased, see Proposition 1. Thus, we cannot use (23) to calculate the variance of OS. As a solution, we propose an unbiased variant of OS in Section 5.4, for which we will calculate the variance.

Algorithm 1 implements (34) as a ready-to-use algorithm with some enhancements. When approximating the Shapley value of all players, these additional improvements reduce the number of evaluations of v, and thus, increase convergence. In detail, by looking at Algorithm 1, one obtains that evaluations of $v\left(\mathbf{z}\right)$ can be reused to update the Shapley value of each individual player, either as the minuend, if player i is already included in the sampled coalition, or as the subtrahend, if excluded. Note that we force the exclusion of player i in the subtrahend to avoid adding marginal contributions of zero. This approach reduces the number of evaluations of the characteristic function per iteration in line 3 from $2n$ to $n+1$. We use this implementation when comparing variances and convergence speeds in Section 6.2.

Algorithm 1:OS

1: ${\widehat{\varphi }}_i\leftarrow 0(\forall i\in N)$ 2: for$q\mathbf{in}0,\dots ,Q-1$do 3:     for ${\tau }_q$ times do 4:         Draw $\mathbf{z}\sim \mathcal{B}\left(\mathbf{1}q/(Q-1)\right)$ 5:         $v_{\mathbf{z}}\leftarrow v\left(\mathbf{z}\right)$ 6:         for $i\mathbf{in}N$ do 7:            if $z_i=1$ then 8:                ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i+v_{\mathbf{z}}-v(\mathbf{z}-{\mathbf{e}}_i)$ 9:            else 10:                ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i+v(\mathbf{z}+{\mathbf{e}}_i)-v_{\mathbf{z}}$ 11:            end if 12:         end for 13:     end for 14: end for 15: ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i/\left(Q{\tau }_q\right)(\forall i\in N)$

Concluding this subsection, we highlight two more details. First, we note that approximations of the integral in (29) other than (32) have also been proposed. Mitchell et al. [33] mentioned that they used the trapezoidal rule in their implementation of OS. On the other hand, Chen et al. [12] mentioned that sampling $q\sim \mathcal{U}\left(\right[0,1\left]\right)$ at random to approximate (29) provides an unbiased estimator. While we assume the former variant to be biased as well, we will further study the latter approach in Section 5.4.

Second, another algorithm proposed by Okhrati and Lipani [22] called halved Owen sampling uses antithetic sampling as a variance reduction technique [12,28,34]. However, this method lies outside the scope of our paper. Similar to the OS estimator, the halved Owen sampling estimator remains biased. Rather than extending all following developments like the construction of an unbiased variant of OS to the antithetic case as well, we focus exclusively on the OS estimator. From the results of this paper, one could easily derive unbiased or stratified antithetic versions of OS. We refer to [34] for more details on how unbiased antithetic variants of existing algorithms can be constructed.

5.4. Multilinear Extension Sampling for the Shapley Value (MES)

As discussed in the previous subsection, the OS estimator is biased, which motivates the introduction of two new unbiased variants of the estimator in this and the following subsection. Chen et al. [12] proposed an algorithm called random q as an alternative to OS. The idea is to sample points across the interval $[0,1]$ at random, resulting in an unbiased estimator of the integral in (29) even for small sample sizes. Although the authors refer to this algorithm at several points in their work, Chen et al. [12], do not provide any theoretical analysis of their algorithm. This lack of theoretical grounding motivates a closer examination of random q along with its advantages and limitations in comparison to other approaches.

From now on, we refer to this algorithm as MES (for Multilinear Extension Sampling). In our implementation, we fix ${\tau }_q=1$ due to the fact that we fail to see any justification in [22] to set ${\tau }_q=2$. In particular, the authors stated in [22] that the optimal choice of ${\tau }_q$ is yet to be explored. In our opinion, there is no clear reason why the algorithm benefits from sampling more than one coalition for a given, sampled q. Therefore, our proposed algorithm is controlled by one parameter $\tau$ only, which specifies the number of random draws of q and thus, implicitly, the number of sampled coalitions. The overall number of evaluations of v is then given by $\tau (n+1)$.

Proposition 3. 

The MES estimator is unbiased, i.e., $\mathbb{E}\left[{\widehat{\varphi }}_i\right]={\varphi }_i,\forall i\in N$.

Proof. 

Let us fix an arbitrary player $i\in N$ whose Shapley value we want to estimate. The sampling procedure of Algorithm 2 is equivalent to sampling $\mathbf{z}\sim \mathcal{B}\left({\overline{\mathbf{e}}}_{i}q\right)$. Therefore, we obtain the estimator

$${\widehat{\varphi }}_i={\displaystyle \frac{1}{\tau }}\sum _{k=1}^{\tau }\left(v({\mathbf{z}}_{i,k}+{\mathbf{e}}_i)-v\left({\mathbf{z}}_{i,k}\right)\right),$$

(35)

where ${\mathbf{z}}_{i,k}$ denotes the k-th sample obtained by first drawing $q\sim \mathcal{U}\left(\right[0,1\left]\right)$ at random, and then sampling ${\mathbf{z}}_{i,k}\sim \mathcal{B}\left({\overline{\mathbf{e}}}_{i}q\right)$.

Now, let $X\sim \mathcal{U}\left(\right[0,1\left]\right)$ denote a random variable that is uniformly distributed across $[0,1]$ with density $g_X\left(q\right)=1$ and concrete realizations $q\in [0,1]$. Then,

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\widehat{\varphi }}_i\right]& =\mathbb{E}\left[{\displaystyle \frac{1}{\tau }}\sum _{k=1}^{\tau }\left(v({\mathbf{z}}_{i,k}+{\mathbf{e}}_i)-v\left({\mathbf{z}}_{i,k}\right)\right)\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{\tau }}\sum _{k=1}^{\tau }\mathbb{E}\left[\left(v({\mathbf{z}}_{i,k}+{\mathbf{e}}_i)-v\left({\mathbf{z}}_{i,k}\right)\right)\right]\hfill \\ \hfill & ={\mathbb{E}}_{X\sim \mathcal{U}\left(\right[0,1\left]\right)}\left[{\mathbb{E}}_{\mathbf{z}\sim \mathcal{B}\left({\overline{\mathbf{e}}}_{i}X\right)}\left[v(\mathbf{z}+{\mathbf{e}}_i)-v\left(\mathbf{z}\right)\right]\right]\hfill \\ \hfill & ={\int }_0^1{\mathbb{E}}_{\mathbf{z}\sim \mathcal{B}\left({\overline{\mathbf{e}}}_{i}q\right)}\left[v(\mathbf{z}+{\mathbf{e}}_i)-v\left(\mathbf{z}\right)\right]\underset{=1}{\underbrace{g_X\left(q\right)}}dq\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill & ={\int }_0^1\sum _{\begin{array}{c}S\subseteq N\\ i\notin S\end{array}}\mathbb{P}(S\mid X=q)\left(v(S\cup \{i\left\}\right)-v\left(S\right)\right)dq\hfill \\ \hfill & =\sum _{\begin{array}{c}S\subseteq N\\ i\notin S\end{array}}{\int }_0^1\mathbb{P}(S\mid X=q)dq\left(v(S\cup \{i\left\}\right)-v\left(S\right)\right)\hfill \\ \hfill & =\sum _{\begin{array}{c}S\subseteq N\\ i\notin S\end{array}}{\int }_0^1{q}^{\left|S\right|}{(1-q)}^{n-\left|S\right|-1}dq\left(v(S\cup \{i\left\}\right)-v\left(S\right)\right)\hfill \\ \hfill & =\sum _{\begin{array}{c}S\subseteq N\\ i\notin S\end{array}}\mathrm{Beta}\left(\right|S|+1,n-\left|S\right|)\left(v(S\cup \{i\left\}\right)-v\left(S\right)\right)\hfill \\ \hfill & =\sum _{\begin{array}{c}S\subseteq N\\ i\notin S\end{array}}{\displaystyle \frac{\left|S\right|!(n-|S|-1)!}{n!}}\left(v(S\cup \{i\left\}\right)-v\left(S\right)\right)\hfill \\ \hfill & ={\varphi }_i,\hfill \end{array}$$

(37)

whereby the closed-form solution of the integral directly comes from the beta function [35,36]. Note that we implicitly switched from vector notation to set notation in the step between (36) and (37), as described in Section 2.1.    □

Algorithm 2:MES

1: ${\widehat{\varphi }}_i\leftarrow 0(\forall i\in N)$ 2: for$\tau$times do 3:     Draw $q\sim \mathcal{U}\left(\right[0,1\left]\right)$ 4:     Draw $\mathbf{z}\sim \mathcal{B}\left(\mathbf{1}q\right)$ 5:     $v_{\mathbf{z}}\leftarrow v\left(\mathbf{z}\right)$ 6:     for $i\mathbf{in}N$ do 7:         if $z_i=1$ then 8:            ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i+v_{\mathbf{z}}-v(\mathbf{z}-{\mathbf{e}}_i)$ 9:         else 10:            ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i+v(\mathbf{z}+{\mathbf{e}}_i)-v_{\mathbf{z}}$ 11:         end if 12:     end for 13: end for 14: ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i/\tau (\forall i\in N)$

Empirical data given in Figure 3 in Section 6 supports the claim that MES is unbiased. Additionally, it is worth noting that MES may perform worse than OS, as shown in Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9. We attribute this observation to the fact that OS generates its q-values systematically, which probably helps reduce its variance. Although we will later derive a theoretical variance for the MES estimator (see Corollary 1), we cannot compare it to OS since a theoretical variance analysis for OS has not been performed yet.

5.5. Stratified Multilinear Extension Sampling for the Shapley Value (S-MES)

In order to reduce the variance of MES without reintroducing a bias, we propose a new algorithm which splits the interval $[0,1]$ into $\tau$ strata, such that the strata are defined as

$$I_k=\left[{\displaystyle \frac{k-1}{\tau }},{\displaystyle \frac{k}{\tau }}\right)(\forall k\in \{1,\dots ,\tau \}),$$

(38)

hence its name S-MES. To avoid confusion, we highlight that the strata are defined on the interval space, not on the subset space. As we demonstrate in Section 5.7, these are different stratification schemes.

In each stratum $I_k$ from (38), a random value $q_k\in I_k$ is drawn uniformly, which serves as the basis to generate a new random coalition based on the Bernoulli distribution $\mathcal{B}\left(\mathbf{1}{q}_k\right)$. Again, similar to MES from the previous section, S-MES is controlled by one parameter $\tau$ only, which defines the number of sampled coalitions. The total amount of evaluations of v is then given by $\tau (n+1)$.

Proposition 4. 

The S-MES estimator is unbiased, i.e., $\mathbb{E}\left[{\widehat{\varphi }}_i\right]={\varphi }_i,\forall i\in N$.

Proof. 

Let us fix an arbitrary player $i\in N$ whose Shapley value we want to estimate. The sampling procedure defined in line 4 of Algorithm 3 is equivalent to sampling $\mathbf{z}\sim \mathcal{B}\left({\overline{\mathbf{e}}}_{i}q\right)$. Therefore, we obtain the estimator

$${\widehat{\varphi }}_i={\displaystyle \frac{1}{\tau }}\sum _{k=1}^{\tau }\left(v({\mathbf{z}}_{i,k}+{\mathbf{e}}_i)-v\left({\mathbf{z}}_{i,k}\right)\right),$$

where ${\mathbf{z}}_{i,k}$ denotes the k-th sample obtained by first drawing $q_k\sim \mathcal{U}\left(I_k\right)$, and then sampling ${\mathbf{z}}_{i,k}\sim \mathcal{B}\left({\overline{\mathbf{e}}}_i{q}_k\right)$.

Now, let us denote by $X_k\sim \mathcal{U}\left(I_k\right)$ random variables that are uniformly distributed across $I_k$ with concrete realizations $q_k\in I_k$ and densities $g_{X_k}\left(q_k\right)=\tau$ whenever $q_k\in I_k$. Then,

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\widehat{\varphi }}_i\right]& =\mathbb{E}\left[{\displaystyle \frac{1}{\tau }}\sum _{k=1}^{\tau }\left(v({\mathbf{z}}_{i,k}+{\mathbf{e}}_i)-v\left({\mathbf{z}}_{i,k}\right)\right)\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{\tau }}\sum _{k=1}^{\tau }{\mathbb{E}}_{X_k\sim \mathcal{U}\left(I_k\right)}\left[{\mathbb{E}}_{\mathbf{z}\sim \mathcal{B}\left({\overline{\mathbf{e}}}_i{X}_k\right)}\left[v(\mathbf{z}+{\mathbf{e}}_i)-v\left(\mathbf{z}\right)\right]\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{\tau }}\sum _{k=1}^{\tau }{\int }_{I_k}{\mathbb{E}}_{\mathbf{z}\sim \mathcal{B}\left({\overline{\mathbf{e}}}_i{q}_k\right)}\left[v(\mathbf{z}+{\mathbf{e}}_i)-v\left(\mathbf{z}\right)\right]\underset{=\tau }{\underbrace{g_{X_k}\left(q_k\right)}}d{q}_k\hfill \\ \hfill & =\sum _{k=1}^{\tau }{\int }_{I_k}{\mathbb{E}}_{\mathbf{z}\sim \mathcal{B}\left({\overline{\mathbf{e}}}_i{q}_k\right)}\left[v(\mathbf{z}+{\mathbf{e}}_i)-v\left(\mathbf{z}\right)\right]d{q}_k\hfill \\ \hfill & ={\int }_0^1{\mathbb{E}}_{\mathbf{z}\sim \mathcal{B}\left({\overline{\mathbf{e}}}_{i}q\right)}\left[v(\mathbf{z}+{\mathbf{e}}_i)-v\left(\mathbf{z}\right)\right]dq.\hfill \end{array}$$

Since the last formula is identical to (36), the rest of the proof proceeds exactly as in Proposition 3 and is therefore omitted here. From there, we obtain

$$\mathbb{E}\left[{\widehat{\varphi }}_i\right]={\int }_0^1{\mathbb{E}}_{\mathbf{z}\sim \mathcal{B}\left({\overline{\mathbf{e}}}_{i}q\right)}\left[v(\mathbf{z}+{\mathbf{e}}_i)-v\left(\mathbf{z}\right)\right]dq={\varphi }_i.$$

   □

Algorithm 3:S-MES

1: ${\widehat{\varphi }}_i\leftarrow 0(\forall i\in N)$ 2: for$k\in \{1,\dots ,\tau \}$do 3:     Draw $q\sim \mathcal{U}\left(I_k\right)$ 4:     Draw $\mathbf{z}\sim \mathcal{B}\left(\mathbf{1}q\right)$ 5:     $v_{\mathbf{z}}\leftarrow v\left(\mathbf{z}\right)$ 6:     for $i\mathbf{in}N$ do 7:         if $z_i=1$ then 8:            ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i+v_{\mathbf{z}}-v(\mathbf{z}-{\mathbf{e}}_i)$ 9:         else 10:            ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i+v(\mathbf{z}+{\mathbf{e}}_i)-v_{\mathbf{z}}$ 11:         end if 12:     end for 13: end for 14: ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i/\tau (\forall i\in N)$

Again, we refer to Figure 3 in Section 6, which validates the unbiasedness of S-MES empirically. Additionally, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 in Section 6 indicate that the variance was reduced compared to MES and appears to be very similar to that of OS, at least in the context of the exemplary games.

Note that we do not derive a theoretical variance for S-MES, but instead propose an improved variant in Section 5.7, for which we will provide the theoretical variance.

5.6. Comparison of MES and Marginal Contribution Importance Sampling (MCIS)

Let us revisit the marginal contribution importance sampling (MCIS) approach for the Shapley value defined by (27), which we introduced at the end of Section 4.2 in the context of our importance sampling framework on the coalition space. We observe

Theorem 2. 

MES is the same importance sampling estimator as MCIS.

Proof. 

Per definition, MCIS is an importance sampling estimator.

Now, let us fit the MES algorithm into the importance sampling framework from Theorem 1. From (35), one obtains that it is based on marginal contributions as well, such that we have ${\nu }_i={\Delta }_i$.

Furthermore, from (35), we retrieve $p_i^{\mathrm{MES}}\left(S\right)=0$, if $i\in S$. For the case $i\notin S$, we derive

$$\begin{array}{cc}\hfill p_i^{\mathrm{MES}}\left(S\right)& ={\mathbb{E}}_{X\sim \mathcal{U}\left(\right[0,1\left]\right)}\left[\mathbb{P}(S\mid X)\right]\hfill \\ \hfill & ={\int }_0^1\mathbb{P}(S\mid X=q)\underset{=1}{\underbrace{g_X\left(q\right)}}dq\hfill \\ \hfill & ={\displaystyle \frac{\left|S\right|!(n-|S|-1)!}{n!}},\hfill \end{array}$$

whereby the last equality directly follows from the proof of Proposition 3, which is why we omit any intermediate steps here.

Since MES updates the estimator by adding ${\Delta }_i\left(S\right)$ with factor 1 (see Algorithm 2), we have $p_i^{\mathrm{MES}}/w_i^{\mathrm{MES}}=1$, and thus, $p_i^{\mathrm{MES}}=w_i^{\mathrm{MES}}$.

As a result, via (26), one obtains $p_i^{\mathrm{MES}}=w_i^{\mathrm{MES}}=p_i^{\mathrm{MCIS}}=w_i^{\mathrm{MCIS}}$ and that completes the proof.    □

Proposition 5. 

The variance of the MCIS estimator is given by

$$Var\left[{\widehat{\varphi }}_i(N,v)\right]={\displaystyle \frac{1}{\tau }}\left(\sum _{\begin{array}{c}S\subseteq N\\ i\notin S\end{array}}{w}_i\left(S\right){\Delta }_i{(S,v)}^2-{\varphi }_i{(N,v)}^2\right)(\forall i\in N).$$

Proof. 

The proof is straightforward by inserting $w_i$, $p_i$, and ${\Delta }_i$ into (23).    □

Corollary 1. 

The variance of the MES estimator is given by

$$Var\left[{\widehat{\varphi }}_i(N,v)\right]={\displaystyle \frac{1}{\tau }}\left(\sum _{\begin{array}{c}S\subseteq N\\ i\notin S\end{array}}{w}_i\left(S\right){\Delta }_i{(S,v)}^2-{\varphi }_i{(N,v)}^2\right)(\forall i\in N).$$

Proof. 

The proof is straightforward taking into account that MES is the same estimator as MCIS (see Theorem 2) and the variance of MCIS is given by Proposition 5.    □

We conclude: MCIS is unbiased, which can be obtained directly from Proposition 1 since MCIS belongs to the importance sampling framework defined by Theorem 1, and MES is unbiased, see Proposition 3. Additionally, via Proposition 5 and Corollary 1, we see that both algorithms share the same variance. As a result, we argue that the indirection introduced by Algorithm 2, i.e., sampling from the multilinear extension of v, does not have any advantages over sampling directly from the coalition space with probabilities $p_i=w_i$.

We refer to Section 6, where Figure 4 validates the proposed variances from Proposition 5 and Corollary 1.

To allow fair comparisons to other algorithms, we reuse evaluations of v across players, as for all other algorithms. We refer to Algorithm 4 for a ready-to-use implementation of this improved variant, where the number of evaluations of v is given by $\tau (n+1)$ and

$$p\left(S\right)={\displaystyle \frac{1}{(n+1)\left(\genfrac{}{}{0pt}{}{n}{\left|S\right|}\right)}}$$

(39)

is the adapted sampling distribution to achieve this improvement.

Proposition 6. 

In the context of Algorithm 4, (39) recovers the original sampling probabilities as defined in (26).

Proof. 

When sampling according to (39), i.e., $\mathcal{S}\sim p$, the probability of sampling a subset of size $s\in \{0,\dots ,n-1\}$ that can be used to update (27) is given by

$$\begin{array}{cc}\hfill \phantom{=}& \mathbb{P}\left((s=|\mathcal{S}\left|\right)\wedge (i\notin \mathcal{S})\right)+\mathbb{P}\left((s=|\mathcal{S}\setminus \left\{i\right\}\left|\right)\wedge (i\in \mathcal{S})\right)\hfill \\ \hfill =& \mathbb{P}(s=|\mathcal{S}\left|\right)\mathbb{P}(i\notin \mathcal{S}\mid s=|\mathcal{S}\left|\right)+\mathbb{P}(s+1=|\mathcal{S}\left|\right)\mathbb{P}(i\in \mathcal{S}\mid s+1=|\mathcal{S}\left|\right)\hfill \\ \hfill =& {\displaystyle \frac{1}{n+1}}{\displaystyle \frac{n-s}{n}}+{\displaystyle \frac{1}{n+1}}{\displaystyle \frac{s+1}{n}}={\displaystyle \frac{1}{n}}.\hfill \end{array}$$

(40)

From (40), we obtain that each coalition size without player i is equally likely, as in (26). Furthermore, from (39), we see that each coalition has the same probability within equal-sized coalitions, such that sampling according to (39) in Algorithm 4 is indeed equivalent to (26). That completes the proof.    □

Algorithm 4:MCIS

1: ${\widehat{\varphi }}_i\leftarrow 0(\forall i\in N)$ 2: for$\tau$times do 3:     Draw $\mathcal{S}\sim p$ 4:     $v_{\mathcal{S}}\leftarrow v\left(\mathcal{S}\right)$ 5:     for $i\mathbf{in}N$ do 6:         if $i\in \mathcal{S}$ then 7:            ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i+v_{\mathcal{S}}-v(\mathcal{S}\setminus \left\{i\right\})$ 8:         else 9:            ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i+v(\mathcal{S}\cup \left\{i\right\})-v_{\mathcal{S}}$ 10:         end if 11:     end for 12: end for 13: ${\widehat{\varphi }}_i\leftarrow {\widehat{\varphi }}_i/\tau (\forall i\in N)$

We acknowledge that Chen et al. [12] already observed that MES is equivalent to sampling coalitions according to (26), which aligns with our Theorem 2. This, however, does not diminish the novelty of our contribution. First, we establish MES within the general importance sampling framework proposed in Theorem 1. Second, the results in this subsection serve only as intermediate steps. While Chen et al. [12] concluded that sampling q at fixed intervals, as in the OS algorithm, outperforms MES and, therefore, MCIS, our contributions go beyond this: We develop an unbiased, stratified version of MES (named S-MES, see Algorithm 3) and compare it against other stratified sampling approaches (see Section 5.7), resulting in new recommendations regarding the usage of multilinear extensions for approximating Shapley values.

5.7. Comparison of S-MES and Stratified Marginal Contribution Importance Sampling (S-MCIS)

In the previous subsection, we argued that there is no benefit of choosing the MES algorithm over the MCIS algorithm. In detail, we stated that sampling from the multilinear extension just mimics the true distribution, and thus, does not have any advantages.

Now, we consider the stratified case, i.e., S-MES from Section 5.5, where the stratification is applied to the interval $[0,1]$. This stratification on the continuous interval $[0,1]$ does not translate directly to the discrete space of coalitions. In S-MES, we first draw $q_k\sim \mathcal{U}\left(I_k\right)$ for each stratum indexed by $k\in \{1,\dots ,\tau \}$. Then, from the perspective of a fixed player $i\in N$, subsets are obtained by sampling i.i.d. according to $\mathcal{B}\left({\overline{\mathbf{e}}}_i{q}_k\right)$.

The discrete stratification scheme most consistent with the continuous scheme outlined above is stratifying by the size $s^{\prime }\in \{0,\dots ,n-1\}$ of coalitions excluding player i. Clearly, sampling from $\mathcal{B}\left({\overline{\mathbf{e}}}_i{q}_k\right)$ is always symmetric in a sense of coalitions of the same size having the same probability of being sampled, and therefore, we argue that a further, more fine-grained stratification in the discrete case would not be supported by S-MES. Conversely, it is also unclear how a less granular stratification could be justified by S-MES, since, clearly, each stratum $I_k$ implies distinct values $q_k$, which imply unequal probabilities for different coalition sizes.

Thus, our proposed approximation named S-MCIS is given by

$${\widehat{\varphi }}_i={\displaystyle \frac{1}{n}}\sum _{s^{\prime }=0}^{n-1}\underset{{\widehat{\varphi }}_{i,s}}{\underbrace{{\displaystyle \frac{1}{{\widehat{\tau }}_{i,s^{\prime }}}}\left({\sum }_{\begin{array}{c}\mathcal{S}\in {\mathbb{S}}_{s^{\prime }}\\ i\notin \mathcal{S}\end{array}}{\Delta }_i\left(\mathcal{S}\right)+{\sum }_{\begin{array}{c}\mathcal{S}\in {\mathbb{S}}_{s^{\prime }+1}\\ i\in \mathcal{S}\end{array}}{\Delta }_i(\mathcal{S}\setminus \left\{i\right\})\right)}},$$

(41)

where

$${\widehat{\tau }}_{i,s^{\prime }}=\sum _{\mathcal{S}\in {\mathbb{S}}_{s^{\prime }}}{\mathbb{1}}_{i\notin \mathcal{S}}+\sum _{\mathcal{S}\in {\mathbb{S}}_{s^{\prime }+1}}{\mathbb{1}}_{i\in \mathcal{S}}$$

is the number of marginal contributions where i joins a coalition of size $s^{\prime }$. Hereby, for all $s\in \{0,\dots ,n\}$, the sample ${\mathbb{S}}_s$ is of size ${\tau }_s$ and taken with uniform probability $p_s=\mathcal{U}\left(\{S\subseteq N\mid s=|S\left|\right\}\right)$.

From the definition above, one obtains that each stratum estimator ${\widehat{\varphi }}_{i,s^{\prime }}$ is the average marginal contribution of player i to subsets of size $s^{\prime }$. Since this algorithm is meant to be the counterpart to S-MES and the stratified version of MCIS, we reuse evaluations of v across players, similar to Algorithm 3 and 4. For simplicity, we aim for a proportional sample allocation (compare Section 3.3) with respect to (26), which is the case when all strata receive the same amount of samples ${\widehat{\tau }}_{i,s^{\prime }}$. Unfortunately, the values ${\widehat{\tau }}_{i,s^{\prime }}$ from (41) are estimators themselves and cannot be specified directly. Instead, we can only control the values ${\tau }_s$. Thus, to retrieve an approximate equal sample allocation across strata, we allocate the same amount of samples to each ${\tau }_s$, which is backed by the following proposition:

Proposition 7. 

When ${\tau }_s=c$ for all $s\in \{0,\dots ,n\}$ and some constant $c\in \mathbb{N}$, the expected sample sizes of all strata are equal, i.e.,

$$\mathbb{E}\left[{\widehat{\tau }}_{i,s^{\prime }}\right]=\mathbb{E}\left[{\widehat{\tau }}_{j,s^{\prime \prime }}\right](\forall i,j\in N,\forall s^{\prime },s^{\prime \prime }\in \{0,\dots ,n-1\}).$$

Proof. 

Let ${\mathcal{S}}_s\subseteq N$ be a random subset of size s obtained via $p_s$. Then, for any $i\in N$ and $s^{\prime }\in \{0,\dots ,n-1\}$, we observe

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\widehat{\tau }}_{i,s^{\prime }}\right]& =\mathbb{E}\left[\sum _{{\mathcal{S}}_{s^{\prime }}\in {\mathbb{S}}_{s^{\prime }}}{\mathbb{1}}_{i\notin {\mathcal{S}}_{s^{\prime }}}\right]+\mathbb{E}\left[\sum _{{\mathcal{S}}_{s^{\prime }+1}\in {\mathbb{S}}_{s^{\prime }+1}}{\mathbb{1}}_{i\in {\mathcal{S}}_{s^{\prime }+1}}\right]\hfill \\ \hfill & ={\tau }_{s^{\prime }}\mathbb{P}(i\notin {\mathcal{S}}_{s^{\prime }})+{\tau }_{s^{\prime }+1}\mathbb{P}(i\in {\mathcal{S}}_{s^{\prime }+1})\hfill \\ \hfill & ={\tau }_{s^{\prime }}{\displaystyle \frac{n-s^{\prime }}{n}}+{\tau }_{s^{\prime }+1}{\displaystyle \frac{s^{\prime }+1}{n}}.\hfill \end{array}$$

Thus, whenever ${\tau }_s=c$ for all $s\in \{0,\dots ,n\}$ with $c\in \mathbb{N}$, we obtain

$$\mathbb{E}\left[{\widehat{\tau }}_{i,s^{\prime }}\right]=c{\displaystyle \frac{n-s^{\prime }}{n}}+c{\displaystyle \frac{s^{\prime }+1}{n}}=c{\displaystyle \frac{n+1}{n}},$$

(42)

such that all $\mathbb{E}\left[{\widehat{\tau }}_{i,s^{\prime }}\right]$ are equal.    □

We refer the reader to Algorithm 5, which implements a ready-to-use implementation including the reuse of samples across players.

Algorithm 5:S-MCIS

1: ${\tilde{\varphi }}_{i,s^{\prime }}\leftarrow 0(\forall i\in N,\forall s^{\prime }\in \{0,\dots ,n-1\})$ 2: ${\widehat{\tau }}_{i,s^{\prime }}\leftarrow 0(\forall i\in N,\forall s^{\prime }\in \{0,\dots ,n-1\})$ 3: for$s\in \{0,\dots ,n\}$do 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 $i\mathbf{in}N$ do 8:            if $i\in \mathcal{S}$ then 9:                ${\tilde{\varphi }}_{i,s-1}\leftarrow {\tilde{\varphi }}_{i,s-1}+v_{\mathcal{S}}-v(\mathcal{S}\setminus \left\{i\right\})$ 10:                ${\widehat{\tau }}_{i,s-1}\leftarrow {\widehat{\tau }}_{i,s-1}+1$ 11:            else 12:                ${\tilde{\varphi }}_{i,s}\leftarrow {\tilde{\varphi }}_{i,s}+v(\mathcal{S}\cup \left\{i\right\})-v_{\mathcal{S}}$ 13:                ${\widehat{\tau }}_{i,s}\leftarrow {\widehat{\tau }}_{i,s}+1$ 14:            end if 15:         end for 16:     end for 17: end for 18: ${\widehat{\varphi }}_{i,s^{\prime }}\leftarrow {\tilde{\varphi }}_{i,s^{\prime }}/{\widehat{\tau }}_{i,s^{\prime }}(\forall i\in N,\forall s^{\prime }\in \{0,\dots ,n-1\})$ 19: ${\widehat{\varphi }}_i\leftarrow {\displaystyle \frac{1}{n}}{\sum }_{s^{\prime }=0}^{n-1}{\widehat{\varphi }}_{i,s^{\prime }}(\forall i\in N)$

We note that S-MCIS is not the first stratified estimator for approximating Shapley values. In particular, we mention the St-ApproShapley algorithm proposed by Castro et al. [37], where stratification was applied based on the position of players in a random permutation. The idea is very similar to the stratification based on coalition sizes, since being in position s in a permutation just means joining a coalition of size $s-1$. Additionally, we highlight the work of Maleki et al. [38], where the authors proposed to stratify by coalition sizes similar to us. However, this manuscript does not aim to compare our S-MCIS algorithm to those and other stratified algorithms. Instead, we developed the S-MCIS algorithm solely for our analysis of sampling algorithms in the context of multilinear extensions.

It is worth noting that it is not guaranteed that ${\widehat{\tau }}_{i,s^{\prime }}>0$ holds for all strata, such that additional checks are required when implementing S-MCIS in the context of real-world applications. Our implementations do not include such checks. Instead, we rely on the following result:

Proposition 8. 

Let τ denote the overall sample budget across all strata. Thus, we have ${\tau }_s\to \infty ,\forall s\in \{0,\dots ,n\}$ as $\tau \to \infty$. Then, with probability approaching one asymptotically, each stratum estimator receives at least one sample, i.e.,

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

Proof. 

For $s^{\prime }\in \{0,\dots ,n-1\}$, the probability of a player $i\in N$ belonging to a random coalition of size $s^{\prime }+1$ is given by $0<{\displaystyle \frac{s^{\prime }+1}{n}}\le 1$, while the probability of not belonging to a random coalition of size $s^{\prime }$ is given by $0<{\displaystyle \frac{n-s^{\prime }}{n}}\le 1$. Thus, we have

$$\begin{array}{cc}\hfill \mathbb{P}({\widehat{\tau }}_{i,s^{\prime }}=0)& =\mathbb{P}\left(\sum _{\mathcal{S}\in {\mathbb{S}}_{s^{\prime }}}{\mathbb{1}}_{i\notin \mathcal{S}}=0\right)+\mathbb{P}\left(\sum _{\mathcal{S}\in {\mathbb{S}}_{s^{\prime }+1}}{\mathbb{1}}_{i\in \mathcal{S}}=0\right)\hfill \\ \hfill & ={\left(1-{\displaystyle \frac{n-s^{\prime }}{n}}\right)}^{{\tau }_{s^{\prime }}}{\left(1-{\displaystyle \frac{s^{\prime }+1}{n}}\right)}^{{\tau }_{s^{\prime }+1}}\stackrel{\begin{array}{c}{\tau }_{s^{\prime }}\to \infty \\ {\tau }_{s^{\prime }+1}\to \infty \end{array}}{\to }0.\hfill \end{array}$$

(43)

Since

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

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

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

□

Now, let us analyze the properties of the S-MCIS estimator:

Proposition 9. 

The S-MCIS estimator is unbiased, i.e., $\mathbb{E}\left[{\widehat{\varphi }}_i\right]={\varphi }_i,\forall i\in N$.

Proof. 

Without loss of generality, we fix one player $i\in N$. We denote by

$${\mathbb{S}}_{i,s^{\prime }}^{\prime }=\left[\mathcal{S}:(\mathcal{S}\in {\mathbb{S}}_{s^{\prime }})\wedge (i\notin \mathcal{S})\right]$$

the vector of all sampled coalitions of size $s^{\prime }$ excluding player i. Similarly, we define

$${\mathbb{S}}_{i,s^{\prime }}^{\prime \prime }=\left[\mathcal{S}\setminus \left\{i\right\}:(\mathcal{S}\in {\mathbb{S}}_{s^{\prime }+1})\wedge (i\in \mathcal{S})\right]$$

as the vector of all sampled coalitions of size $s^{\prime }+1$ that include player i, with player i subsequently removed.

Furthermore, the concatenation of those vectors is denoted by

$${\mathbb{S}}_{i,s^{\prime }}=[{\mathbb{S}}_{i,s^{\prime }}^{\prime },{\mathbb{S}}_{i,s^{\prime }}^{\prime \prime }],$$

with ${\widehat{\tau }}_{i,s^{\prime }}$ being its length.

Thus, with ${\varphi }_{i,s^{\prime }}$ being the true average marginal contribution of player i to all coalitions of size $s^{\prime }$, we get

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\widehat{\varphi }}_i\right]& =\mathbb{E}\left[{\displaystyle \frac{1}{n}}\sum _{s^{\prime }=0}^{n-1}{\widehat{\varphi }}_{i,s^{\prime }}\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{n}}\sum _{s^{\prime }=0}^{n-1}\mathbb{E}\left[{\widehat{\varphi }}_{i,s^{\prime }}\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{n}}\sum _{s^{\prime }=0}^{n-1}\mathbb{E}\left[{\displaystyle \frac{1}{{\widehat{\tau }}_{i,s^{\prime }}}}\left(\sum _{\begin{array}{c}\mathcal{S}\in {\mathbb{S}}_{s^{\prime }}\\ i\notin \mathcal{S}\end{array}}{\Delta }_i\left(\mathcal{S}\right)+\sum _{\begin{array}{c}\mathcal{S}\in {\mathbb{S}}_{s^{\prime }+1}\\ i\in \mathcal{S}\end{array}}{\Delta }_i(\mathcal{S}\setminus \left\{i\right\})\right)\right]\hfill \\ \hfill & ={\displaystyle \frac{1}{n}}\sum _{s^{\prime }=0}^{n-1}\mathbb{E}\left[{\displaystyle \frac{1}{{\widehat{\tau }}_{i,s^{\prime }}}}\left(\sum _{\mathcal{S}\in {\mathbb{S}}_{i,s^{\prime }}^{\prime }}{\Delta }_i\left(\mathcal{S}\right)+\sum _{\mathcal{S}\in {\mathbb{S}}_{i,s^{\prime }}^{\prime \prime }}{\Delta }_i\left(\mathcal{S}\right)\right)\right]\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{1}{n}}\sum _{s^{\prime }=0}^{n-1}\mathbb{E}\left[{\displaystyle \frac{1}{{\widehat{\tau }}_{i,s^{\prime }}}}\sum _{\mathcal{S}\in {\mathbb{S}}_{i,s^{\prime }}}{\Delta }_i\left(\mathcal{S}\right)\right]\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{1}{n}}\sum _{s^{\prime }=0}^{n-1}{\varphi }_{i,s^{\prime }}\hfill \\ \hfill & ={\varphi }_i.\hfill \end{array}$$

Note that the equality between (44) and (45) holds since the expressions inside the expectations in (44) are sample means and the values ${\varphi }_{i,s^{\prime }}$ in (45) are population means. □

Proposition 10. 

The variance of the S-MCIS estimator is approximately given by

$$Var\left[{\widehat{\varphi }}_i\right]\approx {\displaystyle \frac{1}{n^2}}\sum _{s^{\prime }=0}^{n-1}\widehat{Var}\left[{\widehat{\varphi }}_{i,s^{\prime }}\right](\forall i\in N)$$

(46)

with

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

(47)

as long as all ${\tau }_s$ are equal.

Proof. 

Let each ${\widehat{\varphi }}_{i,s^{\prime }}$ be an estimator of a linear solution concept ${\varphi }_{i,s^{\prime }}$ as defined by Theorem 1, i.e.,

$${\varphi }_{i,s^{\prime }}={\displaystyle \frac{1}{\left(\genfrac{}{}{0pt}{}{n-1}{s^{\prime }}\right)}}\sum _{\begin{array}{c}S\subseteq N\\ s^{\prime }=\left|S\right|\\ i\notin S\end{array}}{\Delta }_i\left(S\right),$$

such that, by using (23), the variance of ${\widehat{\varphi }}_{i,s^{\prime }}$ is given by

$$Var\left[{\widehat{\varphi }}_{i,s^{\prime }}\right]={\displaystyle \frac{1}{{\widehat{\tau }}_{i,s^{\prime }}}}\left({\displaystyle \frac{1}{\left(\genfrac{}{}{0pt}{}{n-1}{s^{\prime }}\right)}}\sum _{\begin{array}{c}S\subseteq N\\ s^{\prime }=\left|S\right|\\ i\notin S\end{array}}{\Delta }_i{\left(S\right)}^2-{\varphi }_{i,s^{\prime }}^2\right).$$

(48)

Unfortunately, ${\widehat{\tau }}_{i,s^{\prime }}$ itself is an estimator and not a fixed value. Therefore, we use its expectation given by (42) with $c={\tau }_s$ in order to approximate (48), which results in (47).

Finally, we have

$$Var\left[{\widehat{\varphi }}_i\right]=Var\left[{\displaystyle \frac{1}{n}}\sum _{s^{\prime }=0}^{n-1}{\widehat{\varphi }}_{i,s^{\prime }}\right]={\displaystyle \frac{1}{n^2}}\sum _{s^{\prime }=0}^{n-1}Var\left[{\widehat{\varphi }}_{i,s^{\prime }}\right]\approx {\displaystyle \frac{1}{n^2}}\sum _{s^{\prime }=0}^{n-1}\widehat{Var}\left[{\widehat{\varphi }}_{i,s^{\prime }}\right].$$

□

In the following, we compare the variances of MCIS and S-MCIS.

Proposition 11. 

By assuming that all ${\widehat{\tau }}_{i,s^{\prime }}$ are close enough to $\mathbb{E}\left[{\widehat{\tau }}_{i,s^{\prime }}\right]$ and all ${\tau }_s$ are equal, the variance of the S-MCIS estimator is always less or equal compared to the variance of the MCIS estimator, i.e.,

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

Proof. 

The proof is straightforward taking into account that S-MCIS constitutes a valid stratification of MCIS for each individual player in a sense of the definition provided in Section 3.3. Additionally, by assuming that all ${\widehat{\tau }}_{i,s^{\prime }}$ are close enough to $\mathbb{E}\left[{\widehat{\tau }}_{i,s^{\prime }}\right]$, we derive that S-MCIS uses a proportional sample allocation, since each subset size is equally likely in MCIS (compare (26) and Proposition 6) and the expected sample size per stratum is equal over all strata in S-MCIS (see Proposition 7). Consequently, the variance reduction proof from Section 3.3 holds, which completes the proof. □

Now, we compare the S-MES and S-MCIS estimator. From Propositions 4 and 9, we obtain that both algorithms are unbiased. Proposition 10 provides an approximated theoretical variance for the S-MCIS estimator, but we have not derived a theoretical variance for the S-MES algorithm. Therefore, in order to compare these algorithms and determine how to select among them, we state

Theorem 3. 

By assuming that all ${\widehat{\tau }}_{i,s^{\prime }}$ are close enough to $\mathbb{E}\left[{\widehat{\tau }}_{i,s^{\prime }}\right]$ and all ${\tau }_s$ are equal, the variance of the S-MCIS estimator is always less or equal compared to the variance of the S-MES estimator, i.e.,

$$Var\left[{\widehat{\varphi }}_i^{\mathrm{S}-\mathrm{MCIS}}\right]\le Var\left[{\widehat{\varphi }}_i^{\mathrm{S}-\mathrm{MES}}\right](\forall i\in N).$$

Proof. 

First, we note that the variance of MES and MCIS can be expressed similar to (16), i.e.,

$$\begin{array}{cc}\hfill {Var}_{p_i^{\mathrm{MCIS}}}\left[{\widehat{\varphi }}_i^{\mathrm{MES}}\right]& ={Var}_{p_i^{\mathrm{MCIS}}}\left[{\widehat{\varphi }}_i^{\mathrm{MCIS}}\right]\hfill \end{array}$$

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{1}{\tau }}{\mathbb{E}}_{s^{\prime }\sim \mathcal{U}\left(\{0,\dots ,n-1\}\right)}\left[{Var}_{S\sim p_i^{\mathrm{MCIS}}}[{\Delta }_i\left(S\right)\mid s^{\prime }=\left|S\right|]\right]\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill & \phantom{=}+{\displaystyle \frac{1}{\tau }}{Var}_{s^{\prime }\sim \mathcal{U}\left(\{0,\dots ,n-1\}\right)}\left[{\mathbb{E}}_{S\sim p_i^{\mathrm{MCIS}}}[{\Delta }_i\left(S\right)\mid s^{\prime }=\left|S\right|]\right],\hfill \end{array}$$

(50)

with $p_i^{\mathrm{MCIS}}$ being defined in (26).

Recall that stratification with proportional sample allocation removes the variance between stratum estimators, i.e., (50), such that the overall variance is just the expectation over all stratum estimators’ variances, i.e., (49). We refer to Section 3.3 for more details.

As a result, the variance of the S-MCIS estimator is given by (49) only, since S-MCIS uses a proportional sample allocation. We refer to Propositions 7 and 11, which demonstrate that S-MCIS uses a proportional sample allocation scheme with respect to MCIS when assuming that all ${\widehat{\tau }}_{i,s^{\prime }}$ are close enough to $\mathbb{E}\left[{\widehat{\tau }}_{i,s^{\prime }}\right]$.

In contrast, we argue that the same does not hold for S-MES. Since the stratification is performed over the continuous interval $[0,1]$ and coalitions are sampled in a second stage according to an i.i.d. Bernoulli distribution with parameter $q_k\in I_k$, any coalition size $s^{\prime }$ can theoretically occur in any stratum. Although this approach tends to reduce the variance — because coalition sizes follow a binomial distribution, assigning higher probabilities to subset sizes around $q_k(n-1)$, which could be interpreted as some sort of “weak stratification” — it does not completely eliminate the variance between stratum estimators, i.e., (50). Consequently, the non-negative term given by (50) may indeed be smaller than for its non-stratified variant MES, but unlike in the context of S-MCIS, this term does not vanish in general.

Thus, since both estimators’ variances additionally share the same term given by (49), one obtains that the overall variance of the S-MCIS estimator is always less or equal compared to the variance of the S-MES estimator. □

We refer to Section 6, where we show in Figure 4 that our approximated theoretical variance of S-MCIS from Proposition 10 is very close to the observed empirical variance. Additionally, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 validate Proposition 11 and Theorem 3.

Our main finding in this and the previous subsection is that the indirections via multilinear extensions do not provide any benefit when seeing the characteristic function as a black box model and using Monte Carlo approximations. In detail, we found that one can easily derive the same estimator without the indirections caused by multilinear extensions (see Theorem 2) or derive an estimator with a variance that is less or equal in comparison to that of the respective multilinear-extension-based method, see Theorem 3. Furthermore, dropping the indirections via multilinear extensions allowed us to use the importance sampling framework from Theorem 1 to derive properties like the unbiasedness or variance of an obtained estimator, which in our view also speaks in favor of sampling directly on the coalition space.

6. Empirical Results

To validate our results, we apply the algorithms from Section 5 to approximate the Shapley values of the exemplary games defined in Section 6.1. In detail, these are airport games and weighted voting games from cooperative game theory. Furthermore, we provide a very brief introduction to Baseline Shapley [40] as we later want look at three real-world explainable machine learning problems. In Section 6.2, we report our numerical experiments and generate insights regarding the behavior of different algorithms.

6.1. Example Games from Classical Cooperative Game Theory and our Model for Explainable Artificial Intelligence

In this section, we provide a brief introduction to two selected TU games, i.e., the airport game and the weighted voting game from cooperative game theory, and then provide a very brief introduction to the XAI model we employ.

Airport Game

Littlechild and Thompson [5] proposed airport games in order to distribute the costs of building a runway at an airport between players, whereby each player owns an airplane with a different required runway length. Clearly, these runway lengths imply differing costs, which are represented by the cost vector $\mathbf{c}={[c_1,\dots ,c_n]}^{\top }$, with the i-th component of this vector representing the building (or maintenance) costs of the runway required by the airplane of player i. Formally, the characteristic function of an airport game is given by

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

(51)

Airport games are maintenance problems where the underlying graph is a line. Their closed form solutions for the Shapley value make them ideal test games [41].

For our empirical bias comparison in Figure 3, we define an airport game with $n=23$ players and

$$\mathbf{c}={[1,2,3,7,2,1,2,3,4,4,4,2,1,2,3,1,1,1,8,2,2,3,4]}^{\top }.$$

(52)

6.1.0.2. Weighted Voting Game

Weighted voting games [1,23,36] model scenarios where players with different numbers of seats in a parliament aim to reach a successful vote. These resources are defined by the vector $\mathbf{c}={[c_1,\dots ,c_n]}^{\top }$, with the i-th component of this vector representing the number of seats of player i. If the total number of seats of a coalition is greater or equal than the predefined quota C, the coalition wins, i.e., the vote is successful, and unsuccessful otherwise. Formally, the characteristic function of a weighted voting game is given 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.$$

(53)

There are fast algorithms for computing point-valued solutions of weighted games. We refer to [42,43] for the dynamic programming algorithms and software we employ for the Shapley values.

For our empirical variance validation of MES and MCIS and S-MCIS in Figure 4, we define several weighted voting games with $n=7$ players and

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

(54)

6.1.0.3. Baseline Shapley (BShap) for XAI computations

In the numerical experiments in the following subsection, we adopt the Shapley value framework to attribute a model’s prediction to its input features. The prediction task is treated as a cooperative game where features are treated as players. For a model o and an input $\mathsf{x}$, the characteristic function $v\left(S\right)$ defines the payoff for a coalition of features S.

We employ Baseline Shapley (BShap) (which was originally suggested in [10]) as formalized by Sundararajan and Najmi [40]. In this approach, the value $v\left(S\right)$ is defined deterministically using a single, fixed baseline vector ${\mathsf{z}}^{\prime }$. Missing features are replaced by their corresponding values from this baseline:

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

(55)

with $\overline{S}$ denoting the set of features not in S.

Our specific baseline: In our concrete applications in the following subsection, we choose the baseline ${\mathsf{z}}^{\prime }$ to be the expected value over the training data distribution, i.e., the vector of per-feature averages

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

with $\mathcal{D}$ standing for the training data distribution. This choice centers the Shapley values on the model’s prediction for an “average” data point $o\left(\overline{\mathsf{z}}\right)$ providing a natural reference for interpreting feature contributions. The resulting Shapley values explain the deviation of the prediction $o\left(\mathsf{x}\right)$ from the baseline prediction $o\left(\overline{\mathsf{z}}\right)$.

It is important to distinguish our approach BShap from Random Baseline Shapley (RBShap), where $v\left(S\right)={\mathbb{E}}_{\mathsf{z}\sim \mathcal{D}}\left[o\left([{\mathsf{x}}_S,{\mathsf{z}}_{\overline{S}}]\right)\right]$, see [40]. While RBShap accounts for the full joint distribution when marginalizing features, our use of BShap with the mean baseline offers a computationally efficient and interpretable alternative, defining a clear, data-conditional reference point for explanation. A detailed discussion on how to choose between RBShap and BShap or alternative models is beyond the scope of this article. We refer the reader to [40,44].

6.2. Numerical Experiments

We implemented our algorithms introduced in Section 5 in Python. The implementations and test problems are freely available via the GitHub page of the first author

https://github.com/tim-pollmann/shapley-mcis-mes (accessed on 31 December 2025).

We always assume that we want to estimate all players’ Shapley values. While this may not be the setting in all real-world scenarios, it is a common setting in the context of explainable machine learning, where one wants to explain the current prediction in terms of all features.

Before going into detail, we note that each algorithm has different parameters controlling the number of sampled coalitions and, consequently, the number of evaluations of v. To ensure fair comparisons, we introduce a unified variable T representing the overall sample budget, such that every algorithm performs T evaluations of v, up to negligible rounding errors.

As stated in the previous paragraph, the overall sample budget T defines the total number of evaluations of v that each algorithm is allowed to perform. We now express the algorithm-specific parameters in terms of T for all algorithms from Section 5.

For OS (Algorithm 1), Q specifies the discretization level when approximating (29). For each value $q\in \{0,\dots ,Q-1\}$, ${\tau }_q=2$ coalitions depending on q are sampled, as proposed by Okhrati and Lipani [22]. Since each iteration of those ${\tau }_q$ iterations further requires $n+1$ evaluations of v in order to update all players’ Shapley values, we derive

$$Q=⌈{\displaystyle \frac{T}{{\tau }_q(n+1)}}⌉=⌈{\displaystyle \frac{T}{2(n+1)}}⌉.$$

For MES (Algorithm 2), S-MES (Algorithm 3) and MCIS (Algorithm 4), the parameter $\tau$ controls the number of evaluations of v. Per iteration of $\tau$, $n+1$ evaluations of v are needed to update all players’ Shapley or values, resulting in

$$\tau =⌈{\displaystyle \frac{T}{n+1}}⌉.$$

Finally, in case of S-MCIS (Algorithm 5), there are $n+1$ iterations for different values of $s\in \{0,\dots ,n\}$, and for each s, ${\tau }_s$ iterations are executed, whereby each of those inner iterations requires $n+1$ evaluations of v to update all players’ Shapley values, resulting in

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

We conclude that the maximum deviation from the true total sample budget T is $2n+1$ for OS, n for MES, S-MES and MCIS, as well as ${(n+1)}^2-1$ for S-MCIS. We consider these deviations to be negligible in our subsequent analysis, in particular for large T.

Note that for S-MCIS (Algorithm 5), it is not guaranteed that the algorithm runs successfully in the sense that every stratum receives at least one sample, compare Proposition 8. In the mean squared error comparisons, we require for any $\tau$ that at least half of all runs must be successful for the results to be displayed in the final figure. On the other hand, we do not account for that behavior when executing the variance comparisons in Figure 4. Instead, we rely on large $\tau$ and small n to assume that the probability of S-MCIS failing is close to 0, see Proposition 8.

With the conventions stated above, we conducted our experiments. From Figure 3, one obtains that OS is not unbiased, as stated in Proposition 2. As expected, the bias decreases as T increases, since OS converges to the true Shapley values as $Q\to \infty$ and ${\tau }_q\to \infty$. Additionally, we see that MES as well as S-MES are unbiased, matching the results of Propositions 3 and 4, respectively.

Examining Figure 4, we obtain that the empirical variances of MCIS and MES closely match the theoretical variances, supporting our results from Proposition 5 and Corollary 1, respectively.

A similar result is observed for the S-MCIS algorithm in Figure 4, where the empirically obtained variances closely match the approximated theoretical values stated in Proposition 10. Additionally, we highlight that the theoretical as well as the empirical variances of S-MCIS are always less than those of MCIS and MES in Figure 4, regardless of the concrete parameterization of the underlying cooperative game, which matches the result from Proposition 11.

As for mean squared error comparisons, we first look at an airport game with 100 players specified in Castro et al. [18] and compare the approximation methods from Section 5 in Figure 5. Note that the graph for S-MCIS only starts with a sample size of 70,000 or higher since Proposition 8 fails to succeed frequently enough beforehand. In other words: Since the airport game has $n=100$ players, it becomes harder to guarantee that every stratum receives at least one sample.

Figure 6 compares mean squared errors for our Monte Carlo estimators from Section 5 for a weighted voting game with 50 players, which was previously employed in the software EPIC [42,43]. 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.

In Figure 7 we test the diabetes dataset which is a standard regression benchmark comprising physiological measurements from 442 patients. Each instance includes 10 baseline features. The target variable is a quantitative measure of disease progression one year after baseline. In this XAI experiment, we train a Gradient Boosting Regressor on this dataset and evaluate the approximation methods from Section 5 by comparing their mean squared errors when estimating the Shapley values relative to exact or high-precision reference attributions. The dataset’s moderate size, continuous features, and clinical interpretability make it suitable for benchmarking attribution fidelity in regression tasks.

The California housing dataset employed in the experiments for Figure 8 contains 20,640 census block group entries from the 1990 California census. Each instance includes 8 predictive features. The target variable is the median house value for each block group in hundreds of thousands of dollars. In this XAI evaluation, we train a Multi-Layer Perceptron Regressor (MLP Regressor) on this dataset and benchmark the approximation methods from Section 5 by comparing their mean squared errors in estimating the Shapley values relative to exact or high-precision reference attributions. The dataset’s geographical nature, mixed feature types, and real-world socioeconomic relevance make it appropriate for testing attribution methods on neural networks in regression settings.

Figure 9 deals with the wine dataset, a classic classification benchmark containing 178 instances of chemical analyses from three grape varieties (cultivars) of wine grown in the same region of Italy. Each instance includes 13 continuous features. The target variable is the cultivar type (class 0, 1, or 2). In this XAI experiment, we train a Random Forest Classifier on this dataset and evaluate our Monte Carlo approximation methods from Section 5 by comparing their mean squared errors when estimating feature importance values for predicting the probability of class 0 only. The dataset’s multivariate chemical profiles, clear class structure, and moderate dimensionality make it well-suited for benchmarking attribution fidelity in probability-based classification explanations.

We can summarize the comparisons of mean squared errors of all algorithms for approximating the Shapley value from Section 5 in Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 succinctly. As expected, MCIS has an equal performance compared to MES since they are the same importance sampling estimator, see Theorem 2. Although being biased, OS clearly outperforms those two algorithms, probably due to its systematic sampling strategy, which should help reduce its variance. As one might expect, our proposed unbiased, stratified estimator S-MES performs equally well compared to OS without reintroducing a bias. The algorithm achieving the best performance in all five examples is S-MCIS, which, following Proposition 11, has a less or equal variance (and therefore, mean squared error) in comparison to MCIS, and, according to Theorem 3, is expected to always outperform or at least match S-MES.

7. Summary, Conclusions and Outlook

This manuscript makes several key contributions that advance the theoretical understanding and practical analysis of Monte-Carlo-based Shapley value approximation methods. We summarize our the most important findings and contextualize them critically in the following two paragraphs.

We establish a novel connection between the unified approach proposed by Benati et al. [21] and Monte Carlo importance sampling on the coalition space, as formalized in Theorem 1. While we deem this importance sampling framework to be precious beyond the algorithms discussed in this paper — which are all based on sampling marginal contributions — we most certainly do not mean to claim it was the one and only way to perform importance sampling for computing Shapley values.

We demonstrate that the multilinear-extensions-based sampling algorithms from [12,22] offer no theoretical advantages over algorithms that are based on sampling coalitions, as shown in Theorems 2 and 3. Nevertheless, we wish to emphasize that our results on multilinear extensions do not extend to tensor-based approaches such as [45] which we identify as a promising research direction for fast algorithms for Shapley values. Likewise, our findings on multilinear extensions do not diminish the potential of non-atomic cooperative games [46] for XAI and other machine learning applications, see e.g. [40].

Beyond these two core contributions, this article presents several supporting and individual results that deepen the analysis and understanding of the algorithms covered in this paper. These results include:

• a formal proof that OS [22] is biased (Proposition 2)
• an analysis of MES [12] establishing a theoretical variance for the algorithm (Proposition 3 and Corollary 1)
• the introduction of a new stratified multilinear-extension-based algorithm (S-MES, see Algorithm 3) along with a proof of its unbiasedness (Proposition 4)
• a multilinear-extension-inspired algorithm (MCIS, see Algorithm 4) along with its analysis (Propositions 5 and 6)
• a stratified version of MCIS (S-MCIS, see Algorithm 5) along with a detailed analysis (Propositions 7, 8, 9, 10, and 11).

As we emphasized in the introduction, the goal of this paper is neither to establish a new approach for estimating Shapley values nor to perform comparisons over an array of approximation algorithms. Instead, this work is about structural, theoretical and algorithmic insight. We are aware that it is possible to incorporate the idea of antithetic sampling for variance reduction along the lines of [34] into all the algorithms discussed in this work. However, that would exceed the scope this study. Likewise, we assumed that sample sizes per stratum are proportionally allocated throughout this work as incorporating more sophisticated stratification strategies [37,47,48] would clearly go beyond the scope of our article. The least squares strategy for Shapley values which is actually the primary subject matter of the work by Benati et al. [21], is covered by our importance sampling framework via Theorem 1 and Proposition 1 and warrants further study for explainable artificial intelligence and related machine learning applications.