# Multivariate Bernoulli Hoeffding Decomposition: From Theory to Sensitivity Analysis

Baptiste Ferrere, Nicolas Bousquet, Fabrice Gamboa, Jean-Michel Loubes, Joseph Muré

► To cite this version:

Baptiste Ferrere, Nicolas Bousquet, Fabrice Gamboa, Jean-Michel Loubes, Joseph Muré. Multivariate Bernoulli Hoeffding Decomposition: From Theory to Sensitivity Analysis. 2025. hal-05301780v3

HAL Id: hal-05301780

https://hal.science/hal-05301780v3

Preprint submitted on 4 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Multivariate Bernoulli Hoeffding Decomposition: From Theory to Sensitivity Analysis

Baptiste Ferrere*¹,², Nicolas Bousquet†¹, Fabrice Gamboa‡², Jean-Michel Loubes§²,³, and Joseph Muré¶¹

¹EDF R&D, Chatou, France

²IMT, Toulouse, France

³Regalia Team, INRIA, Toulouse, France

## Abstract

Understanding the behavior of predictive models with random inputs can be achieved through functional decompositions into sub-models that capture interpretable effects of input groups. Building on recent advances in uncertainty quantification, the existence and uniqueness of a generalized Hoeffding decomposition have been established for correlated input variables, using oblique projections onto suitable functional subspaces. This work focuses on the case of Bernoulli inputs and provides a complete analytical characterization of the decomposition. We show that, in this discrete setting, the associated subspaces are one-dimensional and that the decomposition admits a closed-form representation. One of the main contributions of this study is to generalize the classical Fourier-Walsh-Hadamard decomposition for pseudo-Boolean functions to the correlated case, yielding an oblique version when the underlying distribution is not a product measure, and recovering the standard orthogonal form when independence holds. This explicit structure offers a fully interpretable framework, clarifying the contribution of each input combination and theoretically enabling model reverse engineering. From this formulation, explicit sensitivity measures—such as Sobol' indices and Shapley effects—can be directly derived. Numerical experiments illustrate the practical interest of the approach for decision-support problems involving binary features. The paper concludes with perspectives on extending the methodology to high-dimensional settings and to models involving inputs with finite, non-binary support.

**Keywords** : Dependent inputs; Generalized Hoeffding decomposition; Multivariate Bernoulli distribution; Pseudo-Boolean functions; Fourier-Walsh-Hadamard decomposition; Sensitivity analysis.

## 1 Introduction

Let consider a square-integrable variable of interest $Y = G(X)$. Here, $X$ is a $d$-dimensional random vector that follows a multivariate Bernoulli distribution $P_X$, or Ising model [10]. Besides, $G$ may be a very complex function. Very broadly, a legitimate aim is to understand a prediction $Y$ by explaining the interplay between the components of $X$ and $G$. This basic formulation, sometimes described as an instance of the *multiparameter Bernoulli factory problem* [25], arises in numerous methodological and applied problems related to uncertainty quantification (UQ; [42]),

<footnote>*Corresponding author. baptiste.ferrere@edf.fr</footnote>

<footnote>†nicolas.bousquet@edf.fr</footnote>

<footnote>‡fabrice.gamboa@math.univ-toulouse.fr</footnote>

<footnote>§loubes@math.univ-toulouse.fr</footnote>

<footnote>¶joseph.mure@edf.fr</footnote>
including sensitivity analysis (SA; [9]). In statistical learning, for instance, X can encode the topological structure of a graph [11]. For instance, we refer to the representation of the molecular structure of a chemical compound [33]. It can also be a neural network trained to predict a specific property of interest (e.g., therapeutic efficiency), and Y can be the probability that X possesses this property. This formulation also applies to information retrieval [26], community detection [46, 40], vulnerability analysis of computer networks [34], risk analysis in actuarial contexts [2, 1] or the prediction of information or epidemic diffusion [45]. Other applications within this framework include probabilistic Boolean networks (e.g., used to describe genetic interactions [41] or cell behavior [6, 27]), binary decision diagrams [5] typically used for reliability analysis [3, 29, 14], or Bernoulli filters [37], which describe dynamic systems randomly switching on and off. Decision trees—and by extension random forests—admit a natural representation as functions of multivariate Bernoulli random vectors, with components corresponding to internal decision nodes and terminal leaf indicators.

This framework also emerges when attempting to explain predictions made by composite models through the stochastic activation of sub-models using Bernoulli distributions, as stochastic depth and DropPath [18]. These approaches are similar in spirit to the Monte Carlo Dropout [13], where neurons are randomly dropped at test time to estimate uncertainty. By drawing parallels to mixture-of-experts systems with probabilistic gating [39], such strategies aim to analyze how different activation patterns influence model output, thereby identifying influential model pathways.

Therefore, the ability to interpret predictions made by a (possibly black-box) model G with multivariate Bernoulli inputs can serve to multiple purposes in numerous applications, including improving user acceptability, supporting model design, and ensuring regulatory compliance.

Several approaches can be invoked to this end. Numerous tools producing explanations about a model behavior rely on the axiomatic of Cooperative Game Theory (CGT; [20]), as Shapley values [32], proportional marginal values [15] or Winter values [24]. Unfortunately, the application of the CGT to the study of numerical models still suffers from a lack of rigorous justification in terms of non-falsifiability [28]. Besides, as they typically focus on quantities like the variance of Y, such indicators bring somewhat crude, aggregated information on G. Functional decompositions, based on high-dimensional model representations (HDMR; [35]), theoretically allow for a deeper understanding of a model's behavior, without relying on the CGT axiomatic.

Rather, these decompositions often rely on representation theorems resulting on a decomposition of the $L^2$ space. For instance, classical $L^2$ decompositions are Karhunen-Loève's (leading to functional principal component analysis) [38], Cameron-Martin's theorem (generalized chaos expansion) [12], or Mercer's theorem (kernel representation) [44]. Without relying on prior choices; a seminal representative of such major tools is the Hoeffding Decomposition (HD), first introduced in [16], which gives the following decomposition for any square-integrable function

$$G(X) = \sum_{A \subseteq [d]} G_A(X_A), \quad (1)$$

where $[d] := \{1, \dots, d\}$ and $(G_A(X_A))_{A \subseteq [d]}$ is a unique orthogonal decomposition family.

Until recently, its use (as leading to Sobol's indices in SA, e.g. see [9]) suffered from strong restrictive assumptions: the components of X are assumed to be independent. This is often unrealistic in many applications [36]. Alternatively to Cooperative Game Theory (CGT) based approaches and kernel-based techniques that allow for dependency structures between the $X_i$'s to be taken into account [8], several authors proposed various HD relaxations. We refer to [21] for a review of such approaches. Relaxing the requirements of [7] on the structure of $P_X$, [21] finally provided the most refined characterization of the HD conditions, through oblique projections onto subspaces of $L^2$ that generalize the expressions of $G_A(X_A)$. Consequently, the subspaces defined by the set of square-integrable r.v. measurable with respect to $X_A$ in (1) are
hierarchically orthogonal. The Generalized Hoeffding Decomposition (GHD) promoted by these last authors increasingly appears as a practical and essential tool for explaining the behavior of a model with potentially dependent inputs.

While the aforementioned results are mostly theoretical, the present paper focuses on a specific and practically relevant setting where the input vector X follows a multivariate Bernoulli distribution. In this context, the projection subspaces associated with the GHD turn out to be one-dimensional, which yields a transparent and interpretable functional structure. Importantly, this property enables an explicit and tractable construction of the decomposition components.

We show that this construction can scale to high-dimensional settings, and demonstrate that the GHD framework can be effectively exploited for sensitivity and interpretability analyses when dealing with binary, possibly dependent, inputs. Beyond, this setting offers an ideal ground to investigate the computational and practical aspects of model decomposition on an oblique functional base.

The paper is organized as follows. Section 2 recalls the foundations of the GHD and introduces several new theoretical results, including some specific ones in the multivariate Bernoulli case. Section 3 derives sensitivity and interpretability indices from this decomposition. Section 4 presents numerical experiments on both synthetic and application-inspired examples, illustrating the practical benefits of our approach. We conclude with a discussion on the implications for model analysis in modern machine learning pipelines, and outline several directions for future research. All proofs are postponed to the appendix. A supplementary file provides result when the assumption of full support is violated.

# 2 Multivariate Bernoulli Hoeffding Decomposition

## 2.1 Generalized Hoeffding Decomposition (GHD)

The original form of the HD [17] applies to mutually independent real-valued input variables and provides a foundational framework for analyzing the contributions of individual and interacting features in complex predictive models. It forms the basis of many interpretability techniques, particularly those involving variance-based sensitivity analysis. A natural extension of the classical HD consists in relaxing the assumption of mutual independence among input variables. In this setting, the decomposition is no longer based on orthogonal projections onto $\mathbb{L}^2$ subspaces associated with independent features, but instead relies on a weaker notion of *hierarchical orthogonality*. As formalized in [7], this leads to the so-called GHD, which accommodates dependencies between input features under the following main assumption.

**Assumption 1.** Let $\mathcal{B}(\mathbb{R})$ denote the Borelian $\sigma$-algebra on $\mathbb{R}$. There exist $d$ measures $\nu_1, \dots, \nu_d$ on $(\mathbb{R}, \mathcal{B}(\mathbb{R}))$ such that $\mathbf{P}_X \ll \nu := \nu_1 \otimes \dots \otimes \nu_d$. Moreover, denoting by $f$ the density of $X$ with respect to $\nu$, and for any $A \subseteq [d]$, denoting by $f_A$ the density of $X_A$, there exists $M \in ]0, 1]$ such that

$$ \forall A \subseteq [d], \quad f \ge M f_A f_{A^c}. \qquad (2) $$

**Theorem 2.1** (GHD expression, defined in [7]). Let $G: \mathbb{R}^d \to \mathbb{R}$ be a square integrable function of a random input vector $X = (X_1, \dots, X_d)$, whose components may be arbitrarily dependent and satisfy Assumption (1). Then $G(X)$ can be expressed as a unique finite sum

$$ G(X) = \sum_{A \subseteq [d]} G_A(X_A) \qquad (3) $$

where each component $G_A(X_A)$ is a $\nu_A$-measurable function of $X_A$ satisfying the following hierarchical orthogonality condition

$$ \forall B, B' \subseteq [d], \quad B \subsetneq B' \implies \mathbb{E}[G_B(X_B) G_{B'}(X_{B'})] = 0. \qquad (4) $$
This decomposition is called the Generalized Hoeffding Decomposition (GHD) of G(X).

**Remark.** In the case where $X_1, \dots, X_d$ are mutually independent, the decomposition (3) is simply orthogonal and each component $G_A(X_A)$ captures the unique contribution of the feature subset $A \subseteq [d]$

$$ \forall A \subseteq [d], \quad G_A(X_A) := \sum_{B \subseteq A} (-1)^{|A|-|B|} \mathbb{E}[G(X) | X_B]. \qquad (5) $$

When dealing with an input vector supported on a finite space, the existence and uniqueness of the GHD can be ensured under a single assumption involving Assumption (1). Indeed, the following result, adapted from [7], establishes a sufficient condition for this property to hold.

**Corollary 2.1.** Let $X \in \mathbb{R}^d$ be a finite supported random vector, and define its support as

$$ \mathcal{X} := \{x \in \mathbb{R}^d \mid \mathbb{P}(X = x) > 0\}. $$

Assume that there exist non-empty subsets $\mathcal{X}_1, \dots, \mathcal{X}_d \subset \mathbb{R}$ such that $\mathcal{X} = \mathcal{X}_1 \times \dots \times \mathcal{X}_d$. Then, the assumption of Theorem 2.1 is satisfied and we can obtain the GHD of $G(X)$.

## 2.2 Main contribution

From now, we assume that $X = (X_1, \dots, X_d) \sim \mathbf{P}_X$ is a multivariate Bernoulli random vector of full support, meaning that all binary configurations occur with positive probability:

$$ \forall x \in \{0,1\}^d, \quad \mathbf{P}_X(x) > 0, \qquad (6) $$

with marginal probability parameters $\{q_i\}_{i=1,\dots,d} \in (0,1)^d$. Throughout, we assume that the marginal parameters of the Bernoulli random variables are nondegenerate, i.e., they never take the extreme values 0 or 1; in those cases, the input variables would be almost surely constant and yield no information to the model.

Owing to (6), the condition of Corollary 2.1 is satisfied, and the existence and uniqueness of the GHD are guaranteed. For any non-empty subset $A \subseteq [d]$, we write $X_A := (X_i)_{i \in A}$ for the corresponding subvector of $X$, with marginal distribution $\mathbf{P}_A$. We also define $X_\emptyset := 1$ almost surely. We denote by $\langle \cdot, \cdot \rangle := \mathbb{E}[\cdot \times \cdot]$ the scalar product induced by the expectation which is a finite sum of $2^d$ elements if it is computed between any real valued functions of $X$. The Euclidean norm which derives from this scalar product is denoted $\|\cdot\|$.

The main contribution of this paper consists of constructive and tractable representation Theorems 2.2 and 2.3 for the GHD, together with a set of useful derived results. They rely on the existence of the following family $e(X) := \{e_A(X_A)\}_{A \subseteq [d]}$ of random variables.

$$ \forall A \subseteq [d], \quad e_A(X_A) := \frac{(-1)^{\sum_{i \in A} X_i}}{\mathbf{P}_A(X_A)}. \qquad (7) $$

If (6) holds, this family is well defined and is an oblique transformation of the well-known Fourier-Walsh-Hadamard basis discussed in [30]. Remark also that $e_\emptyset(X_\emptyset) = 1$ a.s.

**Remark.** In the independent case, the canonical orthonormal basis which derives from (7) is given by

$$ e^{\perp} := \left\{ \prod_{i \in A} \left( \frac{(-1)^{X_i}}{q_i^{X_i-1/2} (1-q_i)^{1/2-X_i}} \right) \right\}_{A \subseteq [d]} \qquad (8) $$
**Theorem 2.2** (Multivariate Bernoulli Hoeffding Decomposition — MBHD).

Let $G: \{0, 1\}^d \to \mathbb{R}$. Let denote

$$g(X) := (e_A(X_A)G(X))_{A \subseteq [d]}.$$

Let $\mu$ be its mean vector and let $\Gamma = (\Gamma_{A,B})_{A,B \subseteq [d]}$ be the Gram matrix defined by

$$\Gamma_{A,B} := \mathbb{E}[e_A(X_A)e_B(X_B)], \quad A, B \subseteq [d].$$

Then, under the Assumption (6), the GHD of $G(X)$ is given by

$$G(X) = \sum_{A \subseteq [d]} \beta_A e_A(X_A) \qquad (9)$$

where the real $2^d$-vector $(\beta_A)_{A \subseteq [d]}$ is unique and satisfies the linear equation

$$\Gamma\beta = \mu. \qquad (10)$$

A direct consequence of (9) is that, for any $A \subseteq [d]$, $G_A(X_A) = \beta_A e_A(X_A)$.

This theorem is constructive in the sense that it gives a whole explicit functional decomposition of $G(x)$. Obviously $\Gamma$ is symmetric and positive definite, hence inverting the linear system (10) is a well-posed problem. Furthermore, a refined understanding of this theorem shows that the vectors $\{e_A(X_A)\}_{A \subseteq [d]}$ form a hierarchically orthogonal basis of the space of all real-valued functions of $X$, which is the finite dimensional Hilbert space of square integrable functions of $X$ in our specific study case (see Propositions B.2 and B.3 in Appendix), and that each subspace spanned by elements of this basis is one-dimensional : each $G_A(X_A)$ lives in a vectorial line for each $A \subseteq [d]$.

A natural corollary of this theorem, which enhances interpretability, is that if certain components of the $X_i$ do not contribute to the prediction $G(X)$ (no causal effect, although they can be correlated), then the components of $\beta$ associated with any coalition (i.e., subset of components) composed solely of these variables must necessarily be zero. This so-called exclusion property is important in that any classical sensitivity or interpretability measure, which can be derived from the MBHD (see Section 3), naturally inherits it. This helps address known shortcomings commonly observed in such indicators [22, 15]. The following corollary is a generalization of an unpublished result by Clément Benard (Thales cortAIx-Labs) obtained for the Shapley values.

**Corollary 2.2** (Exclusion property). Suppose there exists a strict subset $C \subseteq [d]$ with $|C| \le d$ and a function $H: \{0, 1\}^{|C|} \to \mathbb{R}$ such that $G(X) = H(X_C)$ a.s. Denote the MBHD of $G(X)$ and $H(X_C)$ by

$$G(X) = \sum_{A \subseteq [d]} \beta_A^{(G)} e_A(X_A) \quad \text{and} \quad H(X_C) = \sum_{A \subseteq C} \beta_A^{(H)} e_A(X_A).$$

Then, for any subset $A \subseteq [d]$,

$$\beta_A^{(G)} = \begin{cases} \beta_A^{(H)}, & \text{if } A \subseteq C, \\ 0, & \text{else.} \end{cases}$$

Adopting now a more geometric perspective, we can derive a second representation theorem. It enables us to interpret the coefficients $\beta_A$ in (9) as inner products, and to better understand the notion of *oblique projection*, which generalizes the classical orthogonal projection.

**Corollary 2.3** (Geometric MBHD). Let, for $A \subseteq [d]$, $e_A^*(X) = \sum_{B \subseteq [d]} (\Gamma^{-1})_{A,B} e_B(X_B)$ be the oblique dual vector of $e_A(X_A)$. Then, for any $A, B \subseteq [d]$, we have

$$\langle e_B(X_B), e_A^*(X) \rangle = \mathbf{1}_{\{A=B\}}, \qquad (11)$$

Furthermore, (9) rewrites as

$$G(X) = \sum_{A \subseteq [d]} \langle G(X), e_A^*(X) \rangle e_A(X_A). \qquad (12)$$
# 3 Derivation of feature importance indicators and interaction effects

The classical indicators used for model SA, as Sobol' indices or Shapley values [31] are tools offering a summarized, agglomerated view of the properties of {$G(X) : X \sim P_X$}. Besides, as highlighted by [19], it is worth noting that the use of such indices coming from the UQ framework has become popular to quantify uncertainty in machine learning models and to promote their explainability [36].

As a matter of fact, SA fundamentally investigates the relationship between input variability and output uncertainty [9]. Accordingly it is often focused on discriminating and ranking the input contributions to the variance of Y. By systematically decomposing output variance into functional components, one can identify which input factors—or combinations thereof—most significantly influence model predictions. This variance-based approach allows for a precise mathematical characterization of feature importance and interaction effects. Benefiting from the GHD, numerous indicators can be obviously derived for multivariate Bernoulli inputs. We give some examples in the remainder of this section.

## 3.1 Generalized Sobol' indices

Starting from the GHD (3), one can obtain the following canonical expression

$$ \mathrm{Var}[G(X)] = \mathrm{Cov}[G(X), G(X)] = \sum_{A \subseteq [d]} \mathrm{Cov}[G(X), G_A(X_A)] $$

This observation warrants further investigation of the quantities $\mathrm{Cov}[G(X), G_A(X_A)]$, which exhibit the desirable property of summing up to the total variance. This leads to the following definition.

**Definition 3.1** (Generalized Sobol' indices). For any subset $A \subseteq [d]$, we define the generalized Sobol' index as

$$ S_A := \frac{\mathrm{Cov}[G(X), G_A(X_A)]}{\mathrm{Var}[G(X)]}, \qquad (13) $$

The convention $S_\emptyset = 0$ still holds, and these indices yield a full variance decomposition, satisfying the normalization condition $\sum_{A \subseteq [d]} S_A = 1$. However, because this generalization is defined via a covariance term rather than a variance, negative Sobol' indices may occur: this happens when there is a strong negative correlation between $G_A(X_A)$ and $G(X)$.

Under the GHD assumptions, Sobol' indices $S_A$ can be decomposed as

$$ S_A = S_A^V + S_A^C \qquad (14) $$

where $S_A^V := \frac{\mathrm{Var}[G_A(X_A)]}{\mathrm{Var}[G(X)]}$ is the variance term and $S_A^C := \sum_{B \subseteq [d] : B \neq A} \frac{\mathrm{Cov}[G_A(X_A), G_B(X_B)]}{\mathrm{Var}[G(X)]}$ is the covariance term in the local contribution. Each term $S_A$ can be indeed interpreted as the fraction of the output variance to which $X_A$ contributes. Summing to 1, these new indices offer a clearer view of the sensitivity of $G(X)$ to any subset ("group of inputs") $X_A$ than usual tools in SA, as so-called grouped Sobol' indices [23, 4] which are limited to independent, distinct groups.

**Definition 3.2** (Sobol' matrix). For any pair of subsets $A, B \subseteq [d]$, we define the entry $(A, B)$ of the Sobol' matrix $S$ as

$$ S_{A,B} := \frac{\mathrm{Cov}[G_A(X_A), G_B(X_B)]}{\mathrm{Var}[G(X)]}. \qquad (15) $$

This matrix provides a refined decomposition that distinguishes between variance contributions (diagonal entries) and covariance contributions (off-diagonal entries). Under the GHD framework, the entries satisfy naturally the normalization property. This construction yields a matrix-based generalization of Sobol' indices for dependent input variables, as independently introduced by [7].
* Obviously, in the case where $X$ is a vector with mutually independent components, $S$ is diagonal and $S = \text{diag}((S_A^\perp)_{A \subseteq [d]})$ with $S_A^\perp$ the classical Sobol' index computed from the orthogonal HD

$$S_A^\perp = \sum_{B \subseteq A} (-1)^{|A|-|B|} \frac{\text{Var}[\mathbb{E}[G(X) | X_B]]}{\text{Var}[G(X)]}.$$

* Furthermore, considering a $d$-multivariate Bernoulli input $X \sim \mathbf{P}_X$ and under Assumption (6), then Equation (15) becomes straightforwardly

$$S_{A,B} = \frac{\beta_A \beta_B \Gamma_{A,B}}{\text{Var}[G(X)]}.$$

## 3.2 Generalized Shapley effects

Beyond the choice of the output variance to quantify variable importance, sensitivity indices based on CGT—particularly those employing Shapley values [31] —provide a unique and axiomatically justified attribution of *value* functions of interest by averaging marginal contributions across all variable coalitions. This approach offers a complementary perspective that has proven particularly relevant in interpretable machine learning. We now present an equivalent characterization of Shapley values that proves especially advantageous for our decomposition framework.

**Definition 3.3** (Shapley values via Harsanyi dividends). Consider a finite player set $[d]$ and let $h: 2^{|d|} \to \mathbb{R}$ be an arbitrary set function satisfying the normalization condition $h(\emptyset) = 0$. We define the value function $v$ via

$$v(A) := \sum_{B \subseteq A} h(B), \quad A \subseteq [d]. \qquad (16)$$

In CGT, the pair $(D, v)$ constitutes a cooperative game, where the family $(h(A))_{A \subseteq [d]}$ represents the Harsanyi dividends. Under this construction, the Shapley value of player $i \in [d]$ admits the following closed-form expression

$$\text{Sh}_i(v) = \sum_{\substack{A \subseteq [d] \\ i \in A}} \frac{h(A)}{|A|}. \qquad (17)$$

To bridge with sensitivity analysis (SA), we can focus on value functions $v$ (normalized such that $v(D) = 1$) where $v(A)$ represents the proportion of variance attributable to coalition $A$. The Shapley values associated with this cooperative game are coined Shapley effects [43].

**Definition 3.4** (Shapley effects via generalized Sobol' indices). In general case, one can construct the variance-based cooperative game $(D, v)$ by setting the Harsanyi dividends equal to the generalized Sobol' indices

$$h(A) := S_A, \quad A \subseteq [d].$$

Consequently, the value function becomes $v(A) = \sum_{B \subseteq A} S_B$. Since $S_\emptyset = 0$ and $\sum_{A \subseteq [d]} S_A = 1$ (as established previously), we obtain $v(\emptyset) = 0$ and $v(D) = 1$. The corresponding Shapley effect for variable $X_i$ is given by

$$\text{Sh}_i = \sum_{\substack{A \subseteq [d] \\ i \in A}} \frac{S_A}{|A|}. \qquad (18)$$

In the special case where the generalized Sobol' indices reduce to the classical Sobol' indices from the orthogonal Hoeffding decomposition, these Shapley effects coincide with the standard Shapley effects (see [43] for a review).

**Proposition 3.1.** In the multivariate Bernoulli case, under Assumption (6) with $X \sim \mathbf{P}_X$, then

$$\forall i \in [d], \quad \text{Sh}_i = \frac{1}{\text{Var}[G(X)]} \sum_{A \subseteq [d]: i \in A} \frac{\beta_A}{|A|} \sum_{B \subseteq [d]} \beta_B \Gamma_{A,B}, \qquad (19)$$
# 4 Numerical experiments

## 4.1 Weight importance in a standard Perceptron

To investigate the impact of feature dependence on functional decompositions and variance-based interpretability tools, we design a synthetic classifier inspired by a linear threshold function. More precisely, we define a binary classifier $G: \{0, 1\}^d \to \{-1, +1\}$ as

$$G(X) := \operatorname{sign}(W^\top X + b), \qquad (20)$$

where $X \in \{0, 1\}^d$ with $d = 10$ is a binary input vector, $W \in \mathbb{R}^d$ is a weight vector, and $b \in \mathbb{R}$ is a bias term. The parameters are set as

$$W = 0.1 \cdot (3 - 7 - 2 - 1 \ 5 - 1 \ 3 \ 8 \ 1 \ -9)^\top, \quad b = 0.12.$$

The weight vector is constructed such that the sum of its components is zero. Moreover, some coordinates have much larger absolute weights. For example, features $X_2, X_5, X_8$, and $X_{10}$ are expected to have a strong impact on the variance of $G(X)$, while $X_1, X_3$, and $X_7$ have moderate importance. The remaining variables are expected to have negligible or no effect. This setup enables a controlled evaluation of how feature dependence influences functional decomposition and variance attribution.

To introduce feature dependence, we generate $X$ by thresholding a multivariate Gaussian distribution. For a given correlation parameter $\rho \in [0, 1]$, we define

$$Z^\rho \sim \mathcal{N}_d(0, \Sigma^\rho), \quad \text{with} \quad \Sigma_{i,j}^\rho := \begin{cases} 1 & \text{if } i = j, \\ \rho & \text{if } i \neq j. \end{cases}$$

The binary vector $X^\rho \in \{0, 1\}^d$ is then obtained by thresholding

$$X_i^\rho := \mathbf{1}_{\{Z_i^\rho \le 0\}}, \quad \text{so that} \quad X_i^\rho \sim \text{Bernoulli}(0.5).$$

This construction ensures that all marginal distributions remain equiprobable Bernoulli, while $\rho$ directly controls the level of dependence between features. In the special case where $\rho = 0$, the components of $X$ are independent and identically distributed.

We consider three levels of dependence to evaluate the effect of $\rho$ on functional decomposition accuracy:

<table><caption>Table 1: Correlation values used in each synthetic study case to control feature dependence.</caption><thead><tr><td>Study case</td><td>Correlation ρ</td></tr></thead><tbody><tr><td>Case 1: High dependence</td><td>0.9</td></tr><tr><td>Case 2: Medium dependence</td><td>0.5</td></tr><tr><td>Case 3: Weak dependence</td><td>0.1</td></tr></tbody></table>

We apply the GHD of $G(X^\rho)$ under each of the three dependence values and compare the resulting quantities with those obtained under the (typically incorrect) assumption of independent features. This analysis allows us to quantify how dependence affects variance attribution and to evaluate the error introduced by relying on independence-based approximations commonly used in sensitivity analysis. We will conduct numerical comparisons of various quantities between the dependent models and the independent model to assess the approximation errors incurred by assuming input independence.

First, one can compute the true variance for the three cases and compare it with the variance of the independent case. Indeed, in the independent approximation, only a fraction of global variance is computed since all the covariances terms are supposed to be equal to zero.
<table><caption>Table 2: Absolute and relative variance errors between the dependent input model $X^\rho$ and the independent input model $X^0$ for three dependence scenarios.</caption><thead><tr><th>Metric</th><th>Case 1</th><th>Case 2</th><th>Case 3</th></tr></thead><tbody><tr><td>Var(G(X<sup>ρ</sup>)) - Var(G(X<sup>0</sup>))</td><td>-0.46</td><td>-0.09</td><td>-0.01</td></tr><tr><td>Var[G(X<sup>ρ</sup>)] - Var[G(X<sup>0</sup>)]</td><td>-0.87</td><td>-0.10</td><td>-0.01</td></tr><tr><td>Var[G(X<sup>ρ</sup>)]</td><td></td><td></td><td></td></tr></tbody></table>

Let denote $S^\rho$ the Sobol' matrix in our case and $S^\rho$ the vector of Sobol' indices. Generally, these quantities are not computed directly, but rather approximated under the independence assumption. In our Gaussian setting, this corresponds to $\rho = 0$, leading to the approximation

$$S^\rho \approx S_{\perp}^\rho \text{ with } S_{\perp}^\rho := \frac{\mathrm{Var}[G(X^0)]}{\mathrm{Var}[G(X^\rho)]} S^0$$

The errors arising from matrix approximations, as well as the resulting inaccuracies in the approximation of Sobol' indices, are summarized in the following two tables.

<table><caption>Table 3: Approximation errors under different dependence settings for (a) Sobol' matrices and (b) Sobol' indices.</caption><thead><tr><th>Metric</th><th>Case 1</th><th>Case 2</th><th>Case 3</th></tr></thead><tbody><tr><td>||S<sup>ρ</sup> - S<sub>⊥</sub><sup>ρ</sup>||<sub>1</sub></td><td>1.33</td><td>0.34</td><td>0.05</td></tr><tr><td>||S<sup>ρ</sup> - S<sub>⊥</sub><sup>ρ</sup>||<sub>1</sub></td><td>0.75</td><td>0.56</td><td>0.18</td></tr><tr><td>||S<sup>ρ</sup>||<sub>1</sub></td><td></td><td></td><td></td></tr><tr><td>||S<sup>ρ</sup> - S<sub>⊥</sub><sup>ρ</sup>||<sub>2</sub></td><td>0.96</td><td>0.24</td><td>0.04</td></tr><tr><td>||S<sup>ρ</sup> - S<sub>⊥</sub><sup>ρ</sup>||<sub>2</sub></td><td>0.78</td><td>0.58</td><td>0.15</td></tr><tr><td>||S<sup>ρ</sup>||<sub>2</sub></td><td></td><td></td><td></td></tr></tbody></table><table><thead><tr><th>Metric</th><th>Case 1</th><th>Case 2</th><th>Case 3</th></tr></thead><tbody><tr><td>||S<sup>ρ</sup> - S<sub>⊥</sub><sup>ρ</sup>||<sub>1</sub></td><td>1.60</td><td>0.55</td><td>0.12</td></tr><tr><td>||S<sup>ρ</sup> - S<sub>⊥</sub><sup>ρ</sup>||<sub>1</sub></td><td>1.45</td><td>0.51</td><td>0.12</td></tr><tr><td>||S<sup>ρ</sup>||<sub>1</sub></td><td></td><td></td><td></td></tr><tr><td>||S<sup>ρ</sup> - S<sub>⊥</sub><sup>ρ</sup>||<sub>2</sub></td><td>0.44</td><td>0.11</td><td>0.02</td></tr><tr><td>||S<sup>ρ</sup> - S<sub>⊥</sub><sup>ρ</sup>||<sub>2</sub></td><td>1.89</td><td>0.41</td><td>0.07</td></tr><tr><td>||S<sup>ρ</sup>||<sub>2</sub></td><td></td><td></td><td></td></tr></tbody></table>

(a) Sobol' matrices.

(b) Sobol' indices.

Here, $||.||_1$ denotes the matrix 1-norm, that is, the sum of the absolute values of the entries of the matrix, and $||.||_2$ denotes the matrix 2-norm, that is, the maximum absolute value of the spectrum of the matrix.

To further analyze the effect of feature dependence on interpretability, we compute the Shapley values associated with the model $G(X^\rho)$ for each of the three dependence regimes introduced above. These values reflect the contribution of each input feature to the model's output, taking into account the correlation structure induced by $\rho$. The plots in Figure 1 illustrate the

<figure><img src="image_1.png" /><figcaption>Figure 1: Shapley value bar plots for the three study cases.</figcaption></figure>

generalized Shapley values for each input feature $X_i$, under three different dependence scenarios. This decomposition mirrors the structure of generalized Sobol' indices and takes in account both direct and indirect effects. Ignoring dependence can lead to inflated importance scores, whereas the generalized Shapley framework accounts for amplification and cancellation effects due to input correlations, yielding a more accurate picture of feature relevance.
## 4.2 Two-dimensional parametrized study case

We conduct a controlled numerical experiment to investigate the influence of dependence structure on binary classification performance and subsequent sensitivity analysis. Our experimental framework employs a bivariate copula model with carefully calibrated dependence parameters.

Consider the bivariate random vector $(U, V)$ with uniform marginals on $[0, 1]$ and dependence structure characterized by the Farlie-Gumbel-Morgenstern (FGM) copula [2]:

$$C_{\theta}(u, v) = uv + \theta uv(1 - u)(1 - v), \quad \theta \in [-1, 1]. \qquad (21)$$

We parameterize the dependence through $\rho := \mathbb{P}(U \le 1/2, V \le 1/2)$, establishing the bijective relationship $\theta_{\rho} = 16\rho - 4$. Note that independence corresponds to $\rho = 1/4$ (hence $\theta_{\rho} = 0$). The binary response variable is defined as:

$$Y = \mathbf{1}_{\{U \le 1/2, V \le 1/2\}}, \qquad (22)$$

The Bayes-optimal decision boundary coincides with the rectangular region $[0, 1/2]^2$. We can now apply MBHD to $Y = G(X) := X_1X_2$ where $X_1 := \mathbf{1}_{\{U \le 1/2\}}$ and $X_2 := \mathbf{1}_{\{V \le 1/2\}}$ and we finally obtain:

$$G(X) = \rho + \frac{(-1)^{1-X_1}}{4} + \frac{(-1)^{1-X_2}}{4} + (-\rho)^{1_{\{X_1 \neq X_2\}}} \left(\frac{1}{2} - \rho\right)^{1_{\{X_1 = X_2\}}}, \qquad (23)$$

and corresponding generalized Sobol' indices quantify variance contributions under the true probability measure are given by

$$S_1 = S_2 = \frac{1}{4(1 - \rho)}, \quad S_{\{1,2\}} = \frac{\frac{1}{2} - \rho}{1 - \rho}. \qquad (24)$$

Figure 2 reports the evolution of the indices and of the variance parts. The decomposition supplies an explicit functional reconstruction of $G$, attributes variance consistently under dependence, and reveals the bias incurred when classical sensitivity analysis assumes independence. This

<figure><img src="image_2.png" /><figcaption>Figure 2: Behavior of Sobol' indices (left) and variance components (right) as functions of the dependence parameter $\rho$. The orange curve corresponds to the marginal contributions of $X_1$ and $X_2$, the blue curve to the joint contribution $X_{1,2}$, and the black dotted curve to the total variance.</figcaption></figure>

reconstruction provides interpretative value beyond the tree structure itself: it quantifies the relative importance of features and exposes systematic biases that arise when classical sensitivity analysis ignores dependence. The framework thus complements the classifier's predictive capability with rigorous attribution of variance contributions, demonstrating that accurate model interpretation requires accounting for the underlying probability structure.
## 4.3 Decision-tree rules as correlated Bernoulli inputs

We view a trained decision tree as a measurable map

$$Y = G(X_1, \dots, X_R),$$

where each internal node induces a binary decision rule $X_i \in \{0, 1\}$ (e.g., $X_i = 1_{\{U \le \tau\}}$ for some feature coordinate $U$ and threshold $\tau \in \mathbb{R}$). Because several rules may involve related features or nested thresholds, the vector $(X_1, \dots, X_R)$ is generally correlated rather than independent. This representation allows downstream probabilistic analysis directly on the rule space.

From an empirical standpoint, reliable estimation of the joint law of $(X_1, \dots, X_R)$ from finite samples requires the following controls:

1. *Dimensionality constraint.* Joint analysis should involve only a small subset of rules, since the number of configurations grows exponentially with $R$ relative to the sample size.

2. *Depth limitation.* Bounding tree depth mitigates combinatorial explosion and stabilizes empirical estimates of rule frequencies.

3. *Minimum leaf size.* Each terminal node must contain enough observations to yield statistically meaningful estimates and to avoid degenerate functional dependencies among rules (e.g., $X_i = 1 \Rightarrow X_j = 0$ almost surely), which would violate a nondegenerate Bernoulli structure.

In summary, a decision tree can be formalized as a function of correlated Bernoulli inputs, but practical inference demands tight control of the number of rules, the depth of the tree, and the statistical robustness of leaf-level estimates. The framework will be illustrated on a controlled simulation and on a standard benchmark dataset, enabling a direct application of our probabilistic analysis.

### Mushrooms Classification
We consider the Agaricus-Lepiota dataset from the UCI Machine Learning Repository, comprising $N = 8124$ observations of mushroom specimens with 22 categorical attributes. The binary response variable $Y \in \{0, 1\}$ encodes edibility (0: edible, 1: poisonous), with class near-perfect balanced proportions. Following standard preprocessing protocols, we apply one-hot encoding to the categorical features, yielding 95 binary indicators.

Guided by the theoretical considerations outlined previously, we impose 45 samples per leaf, which is approximately 0.7% of the training size, ensuring statistical validity of the decision rule. As expected, the trained tree achieves excellent performance on the held-out test set and this learned tree structure allows us to identify a very simple set of $R = 5$ binary rules that capture the essential decision logic, which are given by

$$ 
\begin{aligned}
X_1 &= 1_{\{\text{odor} \notin \{\text{none}\}\}}, \\
X_2 &= 1_{\{\text{stalk-root} \notin \{\text{club,rooted}\}\}}, \\
X_3 &= 1_{\{\text{gill-spacing} \neq \text{crowded}\}}, \\
X_4 &= 1_{\{\text{bruises} \neq \text{true}\}}, \\
X_5 &= 1_{\{\text{spore-print-color} \neq \text{green}\}},
\end{aligned}
$$

and the resulting classifier admits the compact representation:

$$Y_{\text{tree}} = G(X_1, \dots, X_5) := X_1 X_2 (X_3 + (1 - X_3) X_4) + (1 - X_1) (1 - X_5). \quad (25)$$

To empirically validate this simplified model, we evaluate its performance on the test set of 1,625 observations. The empirical joint distribution $\hat{P}_X$ exhibits pronounced departures from
<table><caption>Table 4: Performance metrics and correlation structure of the 5-rule Boolean model.</caption><thead><tr><th>Class</th><th>Prec.</th><th>Rec.</th><th>F1</th></tr></thead><tbody><tr><td>Edible</td><td>0.98</td><td>1.00</td><td>0.99</td></tr><tr><td>Poisonous</td><td>1.00</td><td>0.97</td><td>0.99</td></tr><tr><td>Accuracy</td><td colspan="3">0.987</td></tr><tr><td>Macro avg</td><td>0.99</td><td>0.99</td><td>0.99</td></tr></tbody></table>

(a) Classification performance on test set.

<table><thead><tr><th></th><th>X<sub>1</sub></th><th>X<sub>2</sub></th><th>X<sub>3</sub></th><th>X<sub>4</sub></th><th>X<sub>5</sub></th><th>q&#770;<sub>i</sub></th></tr></thead><tbody><tr><td>X<sub>1</sub></td><td>1.00</td><td></td><td></td><td></td><td></td><td>0.57</td></tr><tr><td>X<sub>2</sub></td><td>-0.27</td><td>1.00</td><td></td><td></td><td></td><td>0.91</td></tr><tr><td>X<sub>3</sub></td><td>0.37</td><td>-0.13</td><td>1.00</td><td></td><td></td><td>0.84</td></tr><tr><td>X<sub>4</sub></td><td>0.29</td><td>0.34</td><td>-0.30</td><td>1.00</td><td></td><td>0.58</td></tr><tr><td>X<sub>5</sub></td><td>0.11</td><td>-0.03</td><td>-0.04</td><td>0.11</td><td>1.00</td><td>0.99</td></tr></tbody></table>

(b) Rules correlations and marginals parameters.

independence, with moderate pairwise correlations between odor-related and morphological features.

The high precision and recall for both classes confirm that the function *G* successfully captures the essential decision boundary and the correlation matrix reveals the non-product nature of $\hat{\mathbf{P}}_X$.

We can now apply the MBHD to $G(X)$ under the empirical distribution $\hat{\mathbf{P}}_X$. We compute both generalized Sobol' indices under the empirical measure and Shapley values to quantify feature importance. The results, displayed in Figure 3, reveal a striking hierarchy of feature relevance. These indices demonstrate extreme concentration on the odor rule (none or not) accounting for 78.2% of the total variance. The stalk-root feature rule contributes for 17%, while the remaining features exhibit negligible main effects. All pairwise interaction terms are essentially zero, this is due to the marginal distribution, indeed, $X_2$ and $X_5$ are almost constant equal at 1. The negligible interaction terms suggest that the classifier's variance can be almost entirely attributed to additive main effects, simplifying the interpretation despite the multiplicative structure in the simplified representation.

<figure><img src="image_3.png" /><figcaption>(a) Generalized Sobol' indices</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Generalized Shapley effects</figcaption></figure>

<caption>Figure 3: Sensitivity analysis for mushroom toxicity classification. (a) Generalized Sobol' indices showing dominance of main effects and negligible interactions. (b) Shapley effects confirming the importance hierarchy $X_1 \gg X_2 > \{X_3, X_4, X_5\}$.</caption>

# 5 Discussion

This work examined the possibility of leveraging the generalized Hoeffding decomposition (GHD) to explain the behavior of models *G* whose inputs are multivariate Bernoulli variables. Our main result is the complete specification of the terms of this functional decomposition, which, beyond interpretability, enables transparency and full reverse-engineering of a model *G*, even when it is a black box.
From a theoretical perspective, the GHD considered here relies on the assumption from [7] of a full support for the input joint distribution of Bernoulli probabilities. In other words, all possible covariate configurations are allowed. In an exploratory setting, such an assumption may appear reasonable. However, as the dimension increases, the joint probabilities of variable coalitions decrease rapidly as coalition size grows, rendering the assumption too restrictive. Preliminary experiments in low dimensions—considering the behavior of the decomposition when one or more joint probabilities tend to zero—have been conducted and are summarized in the Supplementary Material. These results suggest the possibility of a practical relaxation of the assumption, consistent with the theoretical relaxation proposed in [21]. Work is ongoing to better understand the equivalences of that study’s results beyond the multivariate Bernoulli case.

From a practical standpoint, high dimensionality is reached in the vast majority of case studies, particularly in machine learning. It is therefore of interest to treat the GHD as a tool for deriving surrogate models whose error with respect to G can be formally controlled. Consequently, a key task is to develop numerically scalable approaches to such controlled approximations, namely approximations that are based on most informative $L^2$ subspaces. Other ideas to fight the curse of dimensionality could be finding clever strategies that exploit any available knowledge of the model G, for example through parallel exploration of dimension blocks or recursive exploration within a parameter hierarchy.

We believe that these two avenues of theoretical and practical improvement are fundamental to fully exploiting this explicit functional decomposition. Moreover, beyond the Bernoulli case, such a decomposition can likely be made explicit in settings where the input distribution is finitely supported. Validation of this intuition, as outlined in the Supplementary Material, would represent a significant step toward a decomposition that applies explicitly to the empirical measure, and thereby to the full range of machine learning problems.

# A General Notations

Recall that the positive integer $d$ denotes the fixed ambient dimension, and define $[d] := \{1, \dots, d\}$ as the set of feature indices. We denote by $|A|$ the cardinality of any subset $A \subseteq [d]$. Let $X = (X_1, \dots, X_d)$ be a vector under a multivariate Bernoulli distribution denoted $\mathbf{P}_X$. For any subset $A \subseteq [d]$, we define the subvector $X_A := (X_i)_{i \in A}$. Since $X$ takes a finite number of values, the Hilbert space of square integrable functions of $X$, denoted $L^2$, is just the space of real-valued measurable functions of $X$. More generally, one can define

$$ \forall A \subseteq [d], \quad L_A^2 := \{ f(X_A) \mid f: \{0,1\}^{|A|} \to \mathbb{R} \}, \qquad (26) $$

in the particular case where $A = [d]$, one has $L_A^2 = L^2$. In the multivariate Bernoulli setting, $L_A^2$ can also be identified with $\mathbb{R}^{2^{|A|}}$. Finally, given a closed subspace $F \subseteq L_A^2$, we denote by $F^{\perp_A}$ its orthogonal complement in $L_A^2$, which is exists since this vector space is a Hilbert space. For any family (or single) of random vectors $(Z_1, \dots, Z_n)$ of $L^2$, we denote by $\text{span}(Z_1, \dots, Z_n)$ the vector space generated by this family.

# B Proofs

*Proof of Corollary 2.1.* We assume only in this proof that $X$ is a general random vector which has finite support. The multivariate Bernoulli framework is so a particular case. We assume that its finite support can be written as a Cartesian product of $d$ non-empty sets, we will show that $X$ satisfies the Assumption (1). Recall that $\mathcal{X} \subsetneq \mathbb{R}^d$ denotes the support of $X$ and it is a
finite subset of $\mathbb{R}^d$. Consider now the following application

$$ \nu : \left\{ \begin{array}{ccc} \mathcal{B}(\mathbb{R}^d) & \to & \mathbb{R}_+ \\ B & \mapsto & |B \cap \mathcal{X}| \end{array} \right. $$

Since $\nu$ is constructed using the counting measure, it obvious that $\nu$ is a positive $\sigma$-finite measure on $(\mathbb{R}^d, \mathcal{B}(\mathbb{R}^d))$. Consider also the following measure on $(\mathbb{R}, \mathcal{B}(\mathbb{R}))$

$$ \nu_i(\cdot) := |\cdot \cap \mathcal{X}_i|. $$

Further, obviously using basic properties of the cardinality of intersections and of cartesian products we have that for $(B_1, \dots, B_d) \in \mathcal{B}(\mathbb{R})^d$,

$$ \nu(B_1 \times \dots \times B_d) = \nu_1(B_1) \times \dots \times \nu_d(B_d) = (\nu_1 \otimes \dots \otimes \nu_d)(B_1 \times \dots \times B_d). $$

So that, $\nu = (\nu_1 \otimes \dots \otimes \nu_d)$. Finally, let $B \in \mathcal{B}(\mathbb{R}^d)$ such that $\nu(B) = 0$. Then, $|B \cap \mathcal{X}| = 0$, which implies $B \subset \mathbb{R}^d \setminus \mathcal{X}$. By the definition of $\mathcal{X}$ and its complement, it follows that $\mathbb{P}(X \in B) = 0$. This allows us to conclude that $\mathbf{P}_X \ll \nu$. We so obtain that the density $f$ of $\mathbf{P}_X$ is given by

$$ \forall x \in \mathbb{R}^d, f(x) = \mathbb{P}(X = x) = \mathbf{P}_X(x) $$

More generally, we have

$$ \forall A \subseteq [d], f_A(\cdot) = \mathbb{P}(X_A = \cdot), $$

so to satisfy the assumption, we look for $M \in ]0, 1]$ such that

$$ \forall A \subseteq [d], f \geq M f_A f_{A^c}. $$

Let define the function that map any $A \subseteq [d]$ and $x \in \mathcal{X}$ to $\frac{f(x)}{f_A(x_A) f_{A^c}(x_{A^c})}$. This function is well-defined since, by the definition of $\mathcal{X}$, the denominator never vanishes. It is always positive and takes a finite number of values. Therefore, it has a positive minimum $M$ and this suffices to finally prove Assumption (1). $\square$

**Proof of Theorem 2.2.** The proof of this theorem requires establishing several intermediate results, which rely on the notion of Hoeffding decomposition spaces (Definition B.1). This concept allow us to characterize the Hoeffding decomposition in terms of $\mathbb{L}^2$ subspaces, denoted by $\{V_A\}_{A \subseteq [d]}$, generated by the effect of subsets $A$. Indeed, from [21], the following decomposition holds for the space $\mathbb{L}^2$, generated by the set of possible realizations of $G(X)$, at least under the conditions of Theorem 2.1

$$ \mathbb{L}^2 = \bigoplus_{A \subseteq [d]} V_A. \qquad (27) $$

Under Assumption (6), we subsequently show that:

(a) these subspaces $V_A$ are all one-dimensional (Proposition B.1);

(b) the family $e(X) = \{e_A(X_A)\}_{A \subseteq [d]}$ is hierarchically orthogonal (Proposition B.2), and that $V_A = \text{span}(e_A(X_A))$, i.e., $e(X)$ forms a basis of the whole space $\mathbb{L}^2$ (Proposition B.3).

This last result implies that $\mathbb{L}^2 = \text{span}((e_A(X_A))_{A \subseteq [d]})$. Hence, there exists a unique family of real numbers $\beta := (\beta_A)_{A \subseteq [d]}$ such that

$$ G(X) = \sum_{A \subseteq [d]} \beta_A e_A(X_A). $$
Then, for any $A \subseteq [d]$,

$$
\begin{align*}
\mu_A &= \mathbb{E}[G(X)e_A(X_A)], \\
&= \mathbb{E}\left[e_A(X_A) \times \sum_{B \subseteq [d]} \beta_B e_B(X_B)\right], \\
&= \sum_{B \subseteq [d]} \mathbb{E}[e_A(X_A)e_B(X_B)] \times \beta_B, \\
&= \sum_{B \subseteq [d]} \Gamma_{A,B}\beta_B, \\
&= (\Gamma\beta)_A
\end{align*}
$$

which leads finally to the equality $\Gamma\beta = \mu$. Moreover, since $\Gamma$ is the Gram matrix of the basis $e(X)$, it is well known that this matrix is symmetric positive definite. $\square$

**Definition B.1** (Hoeffding decomposition spaces (HDS)). Let $V_\emptyset$ denote the vector space of all almost surely constant real-valued random variables. We then recursively define the family of vector spaces $(V_A)_{A\subseteq[d]}$ as follows

$$
\forall A \subseteq [d], V_A := \left( \bigoplus_{B \subseteq A} V_B \right)^{\perp_A} \qquad (28)
$$

**Proposition B.1.** *Assume $X = (X_1, \dots, X_d) \sim \mathbf{P}_X$ is a multivariate Bernoulli random vector of full support (Assumption (6)). Then the associated HDS are all one-dimensional vector lines*

$$
\forall A \subseteq [d], \dim(V_A) = 1.
$$

*Proof of Proposition B.1.* From (27) and under Assumption (6), $\dim(\mathbb{L}^2) = 2^d$, which implies

$$
\sum_{A \subseteq [d]} \dim(V_A) = 2^d.
$$

Since the spaces $V_A$ form a direct sum, it follows that $\dim(V_A) \ge 1$ for every $A \subseteq [d]$. Therefore
the only possibility is having $\dim(V_A) = 1$ for every $A \subseteq [d]$. $\square$

**Proposition B.2.** The family of vectors $(e_A(X_A))_{A\subseteq[d]}$ is hierarchically orthogonal, i.e.

$$
\forall A, B \subseteq [d], B \subsetneq A \implies e_A(X_A) \perp e_B(X_B)
$$

*Proof of Proposition B.2.* It is obvious that $e_A(X_A) \in \mathbb{L}_A^2$. Moreover, since $B \subsetneq A$, we also have $e_B(X_B) \in \mathbb{L}_A^2$. Thus, the computation of the inner product in $\mathbb{L}^2$ becomes easier since it is now
the expectation of the product of two functions of $X_A$

$$ \langle e_A(X_A), e_B(X_B) \rangle = \mathbb{E}[e_A(X_A)e_B(X_B)] \quad (29) $$

$$ = \sum_{x \in \{0,1\}^{|A|}} e_A(x)e_B(x_B) \mathbb{P}(X_A = x), \quad (30) $$

$$ = \sum_{x \in \{0,1\}^{|A|}} \frac{(-1)^{\sum_{t \in A} x_t} \times (-1)^{\sum_{t \in B} x_t}}{\mathbb{P}(X_A = x) \times \mathbb{P}(X_B = x_B)} \times \mathbb{P}(X_A = x), \quad (31) $$

$$ = \sum_{x \in \{0,1\}^{|A|}} \frac{(-1)^{\sum_{t \in A \setminus B} x_t}}{\mathbb{P}(X_B = x_B)}, \quad (32) $$

$$ = \sum_{y \in \{0,1\}^{|B|}} \sum_{x \in \{0,1\}^{|A|} : x_B = y} \frac{(-1)^{\sum_{t \in A \setminus B} x_t}}{\mathbb{P}(X_B = x_B)}, \quad (33) $$

$$ = \sum_{y \in \{0,1\}^{|B|}} \frac{1}{\mathbb{P}(X_B = y)} \sum_{x \in \{0,1\}^{|A|} : x_B = y} (-1)^{\sum_{t \in A \setminus B} x_t}, \quad (34) $$

$$ = \sum_{y \in \{0,1\}^{|B|}} \frac{1}{\mathbb{P}(X_B = y)} \sum_{u \in \{0,1\}^{|A|-|B|}} (-1)^{\sum_{t=1}^{|A|-|B|} u_t}, \quad (35) $$

$$ = \sum_{y \in \{0,1\}^{|B|}} \frac{1}{\mathbb{P}(X_B = y)} \times 0 \quad \text{from Lemma B.2 below,} \quad (36) $$

Note that in the transition from Equation (34) to Equation (35), we use the following equality

$$ \sum_{x \in \{0,1\}^{|A|} : x_B = y} (-1)^{\sum_{t \in A \setminus B} x_t} = \sum_{u \in \{0,1\}^{|A|-|B|}} (-1)^{\sum_{t=1}^{|A|-|B|} u_t}. $$

It is observed that the summand on the left-hand side of the equality does not depend on $y$ but only on the cardinality of $A \setminus B$. Regardless of $y$, we iterate over all elements of $A$ to complete $y$ with the missing combinations from $A \setminus B$, in order to reconstruct the combinations of $A$. This requires fully traversing $A \setminus B$. $\square$

**Lemma B.2.** Let $k$ be a positive integer. Then $\sum_{u \in \{0,1\}^k} (-1)^{\sum_{t=1}^k u_t} = 0$.

*Proof of Lemma B.2.* Let $Z$ be a binomial random variable of parameters $(k, 0.5)$. It is enough to notice that

$$ \sum_{u \in \{0,1\}^k} (-1)^{\sum_{t=1}^k u_t} = 2^k \times \mathbb{E}[(-1)^Z] = 2^k \times \left(\frac{1-1}{2}\right)^k = 0. \quad \square $$

**Proposition B.3.** Under the assumptions of Theorem 2.2

$$ \forall A \subseteq [d], \quad V_A = \operatorname{span}(e_A(X_A)). $$

In other words, the family $(e_A(X_A))_{A \subseteq [d]}$ forms a basis with the hierarchical orthogonality property of $\mathbb{L}^2$.
Proof of Proposition B.3. We proceed by induction on the cardinality of $A \subseteq [d]$ to prove that $e_A(X_A) \in V_A$.

* Base case ($|A| = 0$): By definition, $e_\emptyset(X_\emptyset) = 1$ almost surely which obviously implies $e_\emptyset(X_\emptyset) \in V_\emptyset$.

* Inductive step: Let $c \in \mathbb{N}$ and assume the inductive hypothesis: for all $A \subseteq [d]$ such that $|A| \le c$, we have $e_A(X_A) \in V_A$. Now let $A \subseteq [d]$ with $|A| = c + 1$. We aim to show that $e_A(X_A) \in V_A$. By definition of $V_A$, and following [21], one has

$$V_A = \left( \bigoplus_{B \subseteq A} V_B \right)^{\perp_A} = \bigcap_{B \subseteq A} V_B^{\perp_A}.$$

Since $|B| \le c$ for all $B \subseteq A$, the inductive hypothesis gives $e_B(X_B) \in V_B$. Moreover, since each $V_B$ is a one-dimensional subspace (a line), $V_B = \text{span}(e_B(X_B))$. Therefore, we have

$$e_A(X_A) \in V_A \iff \forall B \subseteq A, e_A(X_A) \perp e_B(X_B).$$

This hierarchical orthogonality condition is satisfied by the family $e(X)$ thanks to the previous proposition, hence $e_A(X_A) \in V_A$. This concludes the induction, and thus the proof. $\square$

Proof of Corollary 2.2. Since we have $G(X) = H(X_C)$ almost surely, by applying MBHD to the both member, one can derive the following equality almost surely

$$\sum_{A \subseteq C} \beta_A^{(H)} e_A(X_A) = \sum_{A \subseteq [d]} \beta_A^{(G)} e_A(X_A),$$

by using basics of linear algebra, that necessarily leads to the result

$$\sum_{A \subseteq C} \underbrace{(\beta_A^{(G)} - \beta_A^{(H)})}_{\text{equal to 0}} e_A(X_A) + \sum_{A \subseteq [d] \setminus C} \underbrace{(\beta_A^{(G)} - 0)}_{\text{equal to 0}} e_A(X_A) = \underbrace{0}_{0 \text{ in } \mathbb{L}^2}.$$

Proof of Corollary 2.3. With $A, B \subseteq [d]$, it comes $e_A^*(X) = \sum_{C \subseteq [d]} \Gamma_{A,C}^{-1} e_C(X_C)$, which implies the following equalities

$$\begin{aligned}
\langle e_A^*(X), e_B(X_B) \rangle &= \mathbb{E}[e_A^*(X)e_B(X_B)] \\
&= \sum_{C \subseteq [d]} \Gamma_{A,C}^{-1} \mathbb{E}[e_C(X_C)e_B(X_B)] \\
&= \sum_{C \subseteq [d]} \Gamma_{A,C}^{-1} \Gamma_{C,B} \\
&= (\Gamma^{-1}\Gamma)_{A,B}.
\end{aligned}$$

Then (12) results from simple algebra. $\square$

## C Intuition for degenerate cases

We study pathological cases in which some configurations have vanishing probability, i.e., there exists $x_0 \in \{0, 1\}^d$ with $p_{x_0} := \mathbb{P}(X = x_0) = 0$. We present case studies and show how to handle the decomposition in such a case. In this section, we provide a detailed and pedagogical analysis of all possible pathologies arising in the cases $d = 2$ and $d = 3$. The aim is to illustrate to the reader the intrinsic practical difficulties associated with these low-dimensional settings. This discussion also serves as a basis for introducing a first heuristic approach to address the issue of vanishing probabilities.
### C.1 Fist case study: $d = 2$

Let $X = (X_1, X_2)$ with marginals $(q_1, q_2)$ and dependence parameter $\rho := \mathbb{E}[X_1X_2]$. Then, the probabilities $p_x$ for each configuration $x \in \{0, 1\}^2$ are given by

$$p_{00} = 1 - q_1 - q_2 + \rho, \quad p_{01} = q_2 - \rho, \quad p_{10} = q_1 - \rho, \quad p_{11} = \rho.$$

The random output $G(X)$ is viewed here as a vector $\mathbf{y} := (y_{00}, y_{01}, y_{10}, y_{11})^{\top}$ containing the possible configurations. Let $E$ be the configuration matrix of $e(X)$, with configurations ordered 00, 01, 10, 11 and the subsets of $[d]$ are given by $\{\emptyset, \{1\}, \{2\}, \{1, 2\}\}$. Since $\mathbf{y} = E\beta$ with $\beta$ a unique four dimensional vector, the configuration matrix for the Generalized Hoeffding Decomposition (GHD) components is given by

$$ (G_A(x_A))_{x \in \{0,1\}^2, A \subseteq [d]} = \begin{bmatrix} \emptyset & \{1\} & \{2\} & \{1,2\} \\ 00 & \beta_0 & \frac{\beta_1}{1-q_1} & \frac{\beta_2}{1-q_2} & \frac{\beta_{1,2}}{1-q_1-q_2+\rho} \\ 01 & \beta_0 & \frac{\beta_1}{1-q_1} & -\frac{\beta_2}{q_2} & -\frac{\beta_{1,2}}{q_2-\rho} \\ 10 & \beta_0 & -\frac{\beta_1}{q_1} & \frac{\beta_2}{1-q_2} & -\frac{\beta_{1,2}}{q_1-\rho} \\ 11 & \beta_0 & -\frac{\beta_1}{q_1} & -\frac{\beta_2}{q_2} & \frac{\beta_{1,2}}{\rho} \end{bmatrix}. $$

We treat the four single-zero cases, parameterized by $\rho \in \{0, q_1, q_2, q_1 + q_2 - 1\}$. In each matrix below, the zero-probability row is highlighted in red; undefined entries arise from division by zero, but only on that null-probability configuration, which is immaterial for evaluation.

<table><caption>Table 5: Degenerate 2 dimensional cases</caption><thead><tr><th colspan="4">Case 1: $\rho = 0 \Rightarrow p_{11} = 0$</th><th colspan="4">Case 2: $\rho = q_1 \Rightarrow p_{10} = 0$</th></tr></thead><tbody><tr><td>$\begin{bmatrix} q_1y_{10} + q_2y_{10} - y_{00}(q_1 + q_2 - 1) & q_1(y_{00} - y_{10}) & q_2(y_{00} - y_{10}) & 0 \\ q_1y_{10} + q_2y_{10} - y_{00}(q_1 + q_2 - 1) & q_1(y_{00} - y_{10}) & q_2y_{00} - q_2y_{10} - y_{00} + y_{10} & 0 \\ q_1y_{10} + q_2y_{10} - y_{00}(q_1 + q_2 - 1) & q_1y_{00} - q_1y_{10} - y_{00} + y_{10} & q_2y_{00} - q_2y_{10} - y_{00} + y_{10} & y_{00} - y_{10} - y_{10} + y_{11} \\ q_1y_{10} + q_2y_{10} - y_{00}(q_1 + q_2 - 1) & q_1y_{00} - q_1y_{10} - y_{00} + y_{10} & q_2y_{00} - q_2y_{10} - y_{00} + y_{10} & y_{00} - y_{10} - y_{10} + y_{11} \end{bmatrix}$</td><td></td><td></td><td></td><td>$\begin{bmatrix} q_1y_{11} - y_{00}(q_2 - 1) - y_{01}(q_1 - q_2) & q_1(y_{01} - y_{11}) & q_2(y_{01} - y_{11}) & 0 \\ q_1y_{11} - y_{00}(q_2 - 1) - y_{01}(q_1 - q_2) & q_1(y_{01} - y_{11}) & q_2y_{00} - q_2y_{11} - y_{00} + y_{11} & 0 \\ q_1y_{11} - y_{00}(q_2 - 1) - y_{01}(q_1 - q_2) & q_1y_{01} - q_1y_{11} - y_{01} + y_{11} & q_2(y_{01} - y_{11}) & undefined \\ q_1y_{11} - y_{00}(q_2 - 1) - y_{01}(q_1 - q_2) & q_1y_{01} - q_1y_{11} - y_{01} + y_{11} & q_2y_{00} - q_2y_{11} - y_{01} + y_{11} & undefined \end{bmatrix}$</td><td></td><td></td><td></td><td></td></tr><tr><th colspan="4">Case 3: $\rho = q_2 \Rightarrow p_{01} = 0$</th><th colspan="4">Case 4: $\rho = q_1 + q_2 - 1 \Rightarrow p_{00} = 0$</th></tr><tr><td>$\begin{bmatrix} q_2y_{11} - y_{00}(q_1 - 1) + y_{10}(q_1 - q_2) & q_1(y_{00} - y_{10}) & q_2(y_{10} - y_{11}) & 0 \\ q_2y_{11} - y_{00}(q_1 - 1) + y_{10}(q_1 - q_2) & q_1(y_{00} - y_{10}) & q_2y_{10} - q_2y_{11} - y_{10} + y_{11} & undefined \\ q_2y_{11} - y_{00}(q_1 - 1) + y_{10}(q_1 - q_2) & q_1y_{00} - q_1y_{10} - y_{00} + y_{10} & q_2(y_{10} - y_{11}) & 0 \\ q_2y_{11} - y_{00}(q_1 - 1) + y_{10}(q_1 - q_2) & q_1y_{00} - q_1y_{10} - y_{00} + y_{10} & q_2y_{10} - q_2y_{11} - y_{10} + y_{11} & 0 \end{bmatrix}$</td><td></td><td></td><td></td><td>$\begin{bmatrix} -y_{01}(q_1 - 1) - y_{10}(q_2 - 1) + y_{11}(q_1 + q_2 - 1) & q_1(y_{00} - y_{11}) & q_2(y_{10} - y_{11}) & undefined \\ -y_{01}(q_1 - 1) - y_{10}(q_2 - 1) + y_{11}(q_1 + q_2 - 1) & q_1(y_{00} - y_{11}) & q_2y_{10} - q_2y_{11} - y_{10} + y_{11} & 0 \\ -y_{01}(q_1 - 1) - y_{10}(q_2 - 1) + y_{11}(q_1 + q_2 - 1) & q_1y_{01} - q_1y_{11} - y_{01} + y_{11} & q_2(y_{10} - y_{11}) & 0 \\ -y_{01}(q_1 - 1) - y_{10}(q_2 - 1) + y_{11}(q_1 + q_2 - 1) & q_1y_{01} - q_1y_{11} - y_{01} + y_{11} & q_2y_{10} - q_2y_{11} - y_{10} + y_{11} & 0 \end{bmatrix}$</td><td></td><td></td><td></td><td></td></tr></tbody></table>

In all four cases, the only problematic column is that of $G_{1,2}(X_1, X_2)$. Divisions by zero occur only on the zero-probability configuration, which is never realized; the other configurations yield zeros for $G_{1,2}(X_1, X_2)$. Hence

$$G_{1,2}(X_1, X_2) = 0 \quad \text{a.s.}$$

Therefore, when exactly one cell probability is zero, the GHD of $G(X)$ reduces to

$$G(X) = \beta_0 + \beta_1 e_1(X_1) + \beta_2 e_2(X_2),$$

which is well defined.

For $d = 2$, more than one zero probability among the four joint probabilities forces degeneracy of one marginal, hence model collapse. For instance, if $p_{00} = p_{01} = 0$, then

$$\mathbb{P}(X_1 = 0) = p_{00} + p_{01} = 0,$$

so $X_1 = 1$ a.s. The model reduces to a function of $X_2$

$$G(X) = \tilde{G}(X_2) := \beta_0 + \beta_2 e_2(X_2).$$

In general, some joint probabilities may vanish, which is *not forbidden*. However, such cases must still satisfy the non-degeneracy conditions on the marginal parameters and preserve the absence
of functional dependence between variables, possibly after reconsidering the model structure or the effective dimension d.

For instance, in the bivariate case with two possibly correlated Bernoulli variables $X_1$ and $X_2$, at most one joint probability can be zero without constraining one of the two marginal degrees of freedom. Beyond this point, setting additional probabilities to zero directly affects the marginal distribution of one variable and induces a functional dependence between them.

## C.2 Second Case Study: $d = 3$

Let $X = (X_1, X_2, X_3) \in \{0, 1\}^3$ with joint cell probabilities $(p_{ijk})_{i,j,k \in \{0,1\}}$. The support of $X$ is defined as

$$ \mathcal{X} = \{(i, j, k) \in \{0, 1\}^3 : p_{ijk} > 0\}, \quad n_{\text{zero}} = 8 - |\mathcal{X}|. $$

Recall that $\mathbb{L}^2$ here is the space of all real functions of $X$ and can be identified as $\mathbb{R}^{|\mathcal{X}|}$. Hence the dimension of $\mathbb{L}^2$ is $8 - n_{\text{zero}}$. We order the subsets of $[d]$ as

$$ \{\emptyset, \{1\}, \{2\}, \{3\}, \{1, 2\}, \{1, 3\}, \{2, 3\}, \{1, 2, 3\}\}. $$

Let $E \in \mathbb{R}^{|\mathcal{X}| \times 8}$ denote the configuration matrix, where rows correspond to realized configurations in $\mathcal{X}$ and columns to the functions $\{e_A : A \subseteq [d]\}$ evaluated on those configurations.

When some cells have zero probability, the corresponding rows are removed. Terms $G_A$ that involve division by a zero-probability cell are undefined, but this is immaterial for evaluation since such configurations never occur. Identifiability depends only on the columns of $E$ restricted to $\mathcal{X}$. A valid generalized Hoeffding decomposition (GHD) on $\mathcal{X}$ exists if and only if one can select an index set $\mathcal{B}$ with $|\mathcal{B}| = |\mathcal{X}|$ such that the submatrix $E_{\mathcal{X},\mathcal{B}}$ has full rank. Under standard non-degeneracy (no perfect functional or stochastic dependence among the inputs), such a choice exists, though it is not unique.

Below we present canonical minimal selections $\mathcal{B}$ classified by the number of missing cells. “Counting” refers to the number of full-rank column sets $\mathcal{B}$ of size $|\mathcal{X}|$, referred to as representations.

We require the basis to preserve hierarchical orthogonality. Consequently, any vector that loses at least one modality must be discarded, as the absence of a modality violates the hierarchical orthogonality constraint. In addition, such a vector effectively resides in a lower-dimensional subspace and can thus be represented as a linear combination of the remaining basis vectors. Therefore, retaining it would introduce redundancy and compromise the uniqueness of the hierarchically orthogonal decomposition. Finally, the vector $e_{123}$ is always discarded, as it systematically loses $n_{\text{zero}}$ modalities.

### Subcase $n_{\text{zero}} = 1$

Without loss of generality, remove $p_{000}$. Then $|\mathcal{X}| = 7$. Generically, any $7 \times 7$ principal submatrix obtained by deleting a single column is nonsingular, yielding $\binom{8}{1} = 8$ valid representations. A natural choice discards the three-way interaction

$$ \mathcal{B} = \{\emptyset, \{1\}, \{2\}, \{3\}, \{1, 2\}, \{1, 3\}, \{2, 3\}\}. $$

### Subcase $n_{\text{zero}} = 2$

Three geometric patterns arise.

**(i) Opposite corners removed (same face)** Example: $p_{000} = p_{111} = 0$. Only the column for $e_{123}$ loses two modalities, while main and pairwise effects remain fully defined. A convenient choice is

$$ \mathcal{B} = \{\emptyset, \{1\}, \{2\}, \{3\}, \{1, 2\}, \{1, 3\}\}. $$
(ii) **Two cells sharing a face** Example: $p_{000} = p_{011} = 0$. Here, numerical checks show that $e_1, e_{12}, e_{13}$ cannot all be retained simultaneously, as $X_1$ appears twice. A convenient choice is

$$ \mathcal{B} = \{\emptyset, \{1\}, \{2\}, \{3\}, \{1, 2\}, \{2, 3\}\}. $$

(iii) **Two adjacent cells on a face** Example: $p_{000} = p_{001} = 0$. Here $e_{12}$ loses one modality. The only admissible representation is

$$ \mathcal{B} = \{\emptyset, \{1\}, \{2\}, \{3\}, \{1, 3\}, \{2, 3\}\}. $$

**Subcase $n_{\text{zero}} = 3$**

Three patterns are possible.

(i) **Three cells on a common face.** Example: $p_{000} = p_{001} = p_{101} = 0$. The term $e_{12}$ and $e_{23}$ each lose one modality. The only admissible representation is

$$ \mathcal{B} = \{\emptyset, \{1\}, \{2\}, \{3\}, \{1, 3\}\}. $$

(ii) **Right-triangle configuration.** Example: $p_{000} = p_{001} = p_{110} = 0$. The term $e_{12}$ loses one modality. The only admissible representation is

$$ \mathcal{B} = \{\emptyset, \{1\}, \{2\}, \{3\}, \{1, 3\}\}. $$

(iii) **Equilateral-triangle configuration.** Example: $p_{001} = p_{100} = p_{111} = 0$. All columns retain full modality coverage. A convenient choice is

$$ \mathcal{B} = \{\emptyset, \{1\}, \{2\}, \{3\}, \{1, 2\}\}. $$

**Subcase $n_{\text{zero}} = 4$**

Up to relabeling, the only nontrivial pattern that avoids functional dependence among $(X_1, X_2, X_3)$ is

$$ p_{011} = p_{101} = p_{110} = p_{111} = 0. $$

Then $|X| = 4$, and the space reduces to main marginal effects and the constant term

$$ \mathcal{B} = \{\emptyset, \{1\}, \{2\}, \{3\}\}. $$

### C.3 General discussion

We do not enumerate all possible zero-patterns. The examples above illustrate representative configurations and provide practical intuition for reconstructing the space $\mathbb{L}^2$ when zero-probability cells occur in dimension three.

Finally, we have restricted attention to cases without functional dependence among the components of $X$. Functional dependence can, however, be detected directly from the pattern of vanishing probabilities. Specifically, a variable $X_i$ is functionally independent of the others if and only if there exist two adjacent configurations—differing only in the coordinate $i$—that both have positive probability. Conversely, if $X_i$ is a deterministic function of the remaining variables, then for every configuration $(x_1, \dots, x_{i-1}, x_{i+1}, \dots, x_d)$, exactly one of the two cells

$$ (x_1, \dots, x_{i-1}, 0, x_{i+1}, \dots, x_d) \text{ and } (x_1, \dots, x_{i-1}, 1, x_{i+1}, \dots, x_d) $$

has nonzero probability, while the other must vanish. This provides an immediate combinatorial criterion for detecting functional dependence from the support pattern of the joint distribution.

In general, many zero-probability configurations are possible, and some are more pathological than others. A useful guiding intuition preserves two principles:
* Keep the model as simple as possible (favorize singleton terms over three-way interactions).

* Remove obviously redundant directions (if a direction loses an entire modality, it is natural to drop it).

We propose the following naive approach, which is guaranteed to terminate but is practical only in low dimension (typically $d < 15$).

<caption>**Algorithm 1: Rank-based MBHD construction**</caption>

**Require:** Joint law $P_X$ on $\{0,1\}^d$; support $\mathcal{X}$; arbitrary ordered family $(A_1, \dots, A_{2d})$ of all subsets $[d]$ with $A_1 = \emptyset$ and $A_{2d} = [d]$ by convention.

**Ensure:** $\mathcal{B}$ such that $\{e_A(X_A) : A \in \mathcal{B}\}$ is a basis of $\mathbb{L}^2$ identified as $\mathbb{R}^r$.

1: $r \leftarrow |\mathcal{X}|$; $\mathcal{B} \leftarrow \{\}$; $E \leftarrow []$; $m \leftarrow 1$

2: **while** $|\mathcal{B}| < r$ **do**

3: $c \leftarrow (e_{A_m}(x))_{x \in \mathcal{X}}$ {column vector of $\mathbb{R}^r$}

4: **if** $\text{rank}([E \ c]) = \text{rank}(E) + 1$ **then**

5: $E \leftarrow [E \ c]$

6: $\mathcal{B} \leftarrow \mathcal{B} \cup \{A_m\}$

7: **end if**

8: $m \leftarrow m + 1$

9: **end while**

10: **return** $\mathcal{B}$

We apply this algorithm to confirm the intuition proposed in the previous subsection for $d = 3$, and this intuition is validated in the table 6.

# D Geometric interpretation of MBHD in $\mathbb{R}^2$

In this section we restrict to the case $d = 2$, and consider the classical setting where the joint distribution of the vector $X := (X_1, X_2)$ is given by

$$ \mathbb{P}(X_1 = 0, X_2 = 0) = \mathbb{P}(X_1 = 1, X_2 = 1) = \rho, $$

$$ \mathbb{P}(X_1 = 0, X_2 = 1) = \mathbb{P}(X_1 = 1, X_2 = 0) = \frac{1}{2} - \rho. $$

In particular, both $X_1$ and $X_2$ are marginally distributed as Bernoulli(1/2), and the joint law of $X$ is completely determined by the correlation parameter $\rho \in [0, 1/2]$. In the degenerate case $\rho = 0$, one has $X_2 = 1 - X_1$ almost surely, while for $\rho = 1/2$, one has $X_2 = X_1$ almost surely. Independence occurs if and only if $\rho = 1/4$.

In this setting, we can express $e(X)$ as follows

$$ e_0(X_0) = 1, $$

$$ e_1(X_1) = 2(-1)^{X_1}, $$

$$ e_2(X_2) = 2(-1)^{X_2}, $$

$$ e_{12}(X_{12}) = \frac{(-1)^{X_1+X_2}}{\rho^{1_{\{X_1=X_2\}}} (1-\rho)^{1_{\{X_1 \neq X_2\}}}}. $$

Since the family $e(X)$ forms a basis of $\mathbb{L}^2$, the subfamily $(e_1(X_1), e_2(X_2))$ is two-dimensional. Consequently, any element in $\text{span}(e_1(X_1), e_2(X_2))$ enjoys the appealing property of being representable in the usual Euclidean plane.

We fix, without loss of generality, the direction and orientation of $e_1(X_1)$ by placing it canonically on the horizontal axis with the usual positive orientation; all subsequent objects are defined
<figure><img src="image_5.png" /><figcaption>Figure 4: Example of representation in $\mathbb{R}^2$ for 3 dependence levels.</figcaption></figure>

relative to this choice. Recall that we work in the Hilbert space $\mathbb{L}^2$ endowed with its standard inner product given by the expectation of the product. Consequently,

$$ \|e_1(X_1)\| = \|e_2(X_2)\| = 2. $$

Moreover, for two vectors $u$ and $v$ in a Hilbert space, the angle $\theta_{u,v} \in [0, \pi]$ is defined by

$$ \theta_{u,v} := \arccos\left(\frac{\langle u, v \rangle}{\|u\| \|v\|}\right). $$

In our setting, the (positively oriented, i.e., counterclockwise) angle from $e_1(X_1)$ to $e_2(X_2)$ is therefore

$$ \arccos(4\rho - 1). $$

Now, let consider an illustrative example. Let $G$ be the function such that

$$ G(X) := \frac{3}{4}e_1(X_1) + \frac{1}{4}e_2(X_2). $$

It is immediate that $\mathbb{E}[G(X)] = 0$ and $\mathrm{Var}[G(X)] = 1 + 6\rho$. Accordingly, we may depict $G(X)$ in the plane for the three dependence levels introduced above and we have the following graphical representation. It is apparent that $G(X)$ places substantially more weight on $e_1(X_1)$ than on

<figure><img src="image_6.png" /><figcaption>Figure 5: Representation of $G(X)$ in span $(e_1(X_1), e_2(X_2))$ for 3 dependence levels.</figcaption></figure>

$e_2(X_2)$, and a visual inspection confirms that its norm $\|G(X)\|$ increases with $\rho$. We now consider two carefully chosen cases that highlight two distinct phenomena, depicted below. These two configurations illustrate complementary behaviors. For $\rho_1 = 0.02$, the system approaches the degenerate regime where $X_2 = 1 - X_1$ almost surely. The vectors $e_1(X_1)$ and $e_2(X_2)$ become nearly anti-aligned, reflecting strong negative dependence. Yet the combination $G(X)$ remains well defined, with variance slightly above one.

At $\rho_2 = 1/6$, dependence persists but the geometry becomes orthogonal: $G(X)$ and $e_2(X_2)$ are uncorrelated in $\mathbb{L}^2$. This decorrelation results from an exact cancellation between the direct contribution of $e_2(X_2)$ and the indirect effect transmitted through $e_1(X_1)$.
<figure><img src="image_7.png" /><figcaption>Figure 6: Representation of 2 interesting cases : an extreme negative correlation between $e_1(X_1)$ and $e_2(X_2)$ (left) ; a special case where $G(X)$ and $e_2(X_2)$ are decorrelated.</figcaption></figure>

The first-order Sobol' indices are

$$S_1 = \frac{3}{2} \frac{1 + 2\rho}{1 + 6\rho}, \quad S_2 = \frac{1}{2} \frac{-1 + 6\rho}{1 + 6\rho}.$$

Their values mirror the geometry: for $\rho < 1/6$, $S_2 < 0$ since $e_2(X_2)$ acts against the dominant direction of $G(X)$; when $\rho = 1/6$, orthogonality implies $S_2 = 0$.

<figure><img src="image_8.png" /><figcaption>Figure 7: Behavior of generalized Sobol' indices (left) and variance components (right) as functions of the dependence parameter $\rho$. The blue curve corresponds to $S_1$, the orange curve to $S_2$ (note: $S_2 < 0$ for $\rho < 1/6$), and the black dotted curve to the total variance.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Orientation of the cube</figcaption></figure>

<table><thead><tr><th>Patterns</th><th>Base B</th></tr></thead><tbody><tr><td></td><td>B = {Ø, {1}, {2}, {3}, {1, 2}, {1, 3}, {2, 3}}</td></tr><tr><td></td><td>B = {Ø, {1}, {2}, {3}, {1, 2}, {1, 3}}</td></tr><tr><td></td><td>B = {Ø, {1}, {2}, {3}, {1, 2}, {2, 3}}</td></tr><tr><td></td><td>B = {Ø, {1}, {2}, {3}, {1, 3}, {2, 3}}</td></tr><tr><td></td><td>B = {Ø, {1}, {2}, {3}, {1, 3}}</td></tr><tr><td></td><td>B = {Ø, {1}, {2}, {3}, {1, 3}}</td></tr><tr><td></td><td>B = {Ø, {1}, {2}, {3}, {1, 2}}</td></tr><tr><td></td><td>B = {Ø, {1}, {2}, {3}}</td></tr></tbody></table>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

**Acknowledgements and Funding.** The authors warmly thank Clément Benard (Thales cortAIx-Labs) for kindly agreeing to the publication of Corollary 2.2, which was certainly inspired by his own earlier result on the exclusion property of the Shapley indices derived from the generalized Hoeffding decomposition. This work was supported by the French Association Nationale de la Recherche et de la Technologie (ANRT) through a CIFRE PhD project at Électricité de France (EDF). Fabrice Gamboa and Jean-Michel Loubes are supported by the ANR-3IA Artificial and Natural Intelligence Toulouse Institute (ANITI).
## References

[1] Christopher Blier-Wong, Hélène Cossette, Sebastien Legros, and Etienne Marceau. A new method to construct high-dimensional copulas with bernoulli and coxian-2 distributions. *Journal of Multivariate Analysis*, 201:105261, 2024.

[2] Christopher Blier-Wong, Hélène Cossette, and Etienne Marceau. Stochastic representation of fgm copulas using multivariate bernoulli random variables. *Computational Statistics & Data Analysis*, 173:107506, 2022.

[3] Marc Bouissou. An ordering heuristic for building binary decision diagrams from fault-trees. In *Proceedings of 1996 Annual Reliability and Maintainability Symposium*, pages 208–214. IEEE, 1996.

[4] Baptiste Broto, François Bachoc, Marine Depecker, and Jean-Marc Martinez. Sensitivity indices for independent groups of variables. *Mathematics and Computers in Simulation*, 163:19–31, 2019.

[5] Randal E Bryant. Binary decision diagrams. *Handbook of model checking*, pages 191–217, 2018.

[6] Laurence Calzone, Laurent Tournier, Simon Fourquet, Denis Thieffry, Boris Zhivotovsky, Emmanuel Barillot, and Andrei Zinovyev. Mathematical modelling of cell-fate decision in response to death receptor engagement. *PLoS computational biology*, 6(3):e1000702, 2010.

[7] Gaëlle Chastaing, Fabrice Gamboa, and Clémentine Prieur. Generalized Hoeffding-Sobol Decomposition for Dependent Variables – Application to Sensitivity Analysis. *Electronic Journal of Statistics*, 6:2420–2448, March 2012.

[8] Sébastien Da Veiga. Global sensitivity analysis with dependence measures. *Journal of Statistical Computation and Simulation*, 85(7):1283–1305, 2015.

[9] Sébastien Da Veiga, Fabrice Gamboa, Bertrand Iooss, and Clémentine Prieur. *Basics and Trends in Sensitivity Analysis: Theory and Practice in R*. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2021.

[10] Bin Dai, Shilin Ding, and Grace Wahba. Multivariate bernoulli distribution. *Bernoulli*, 19(4):1465–1483, 2013.

[11] Xiang Deng and Zhongfei Zhang. Graph-free knowledge distillation for graph neural networks. *Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence (IJCAI-21)*, abs/2105.07519, 2021.

[12] Oliver G Ernst, Antje Mugler, Hans-Jörg Starkloff, and Elisabeth Ullmann. On the convergence of generalized polynomial chaos expansions. *ESAIM: Mathematical Modelling and Numerical Analysis*, 46(2):317–339, 2012.

[13] Yarin Gal and Zoubin Ghahramani. Dropout as a bayesian approximation: Representing model uncertainty in deep learning. In *International Conference on Machine Learning*, pages 1050–1059. PMLR, 2016.

[14] Fausto Pedro García Márquez, Isaac Segovia Ramírez, Behnam Mohammadi-Ivatloo, and Alberto Pliego Marugán. Reliability dynamic analysis by fault trees and binary decision diagrams. *Information*, 11(6):324, 2020.
[15] Margot Herin, Marouane Il Idrissi, Vincent Chabridon, and Bertrand Iooss. Proportional marginal effects for global sensitivity analysis. *SIAM/ASA Journal on Uncertainty Quantification*, 12(2):667-692, 2024.

[16] Wassily Hoeffding. A Class of Statistics with Asymptotically Normal Distribution. *The Annals of Mathematical Statistics*, 19(3):293-325, 1948.

[17] Wassily Hoeffding. A class of statistics with asymptotically normal distribution. *Annals of Mathematical Statistics*, 19(3):293-325, 1948.

[18] Gao Huang, Yu Sun, Zhuang Liu, Daniel Sedra, and Kilian Q Weinberger. Deep networks with stochastic depth. In *European Conference on Computer Vision*, pages 646-661. Springer, 2016.

[19] Marouane Il Idrissi. *Development of interpretability methods for certifying machine learning models applied to critical systems*. Theses, Université de Toulouse, 2024.

[20] Marouane Il Idrissi, Nicolas Bousquet, Fabrice Gamboa, Bertrand Iooss, and Jean-Michel Loubes. On the coalitional decomposition of parameters of interest. *Comptes Rendus. Mathématique*, 361(G10):1653-1662, 2023.

[21] Marouane Il Idrissi, Nicolas Bousquet, Fabrice Gamboa, Bertrand Iooss, and Jean-Michel Loubes. Hoeffding decomposition of functions of random dependent variables. *Journal of Multivariate Analysis*, page 105444, 2025.

[22] Bertrand Iooss and Clémentine Prieur. Shapley effects for sensitivity analysis with cor-
related inputs: comparisons with Sobol' indices, numerical estimation and applications.
*International Journal for Uncertainty Quantification*, 9(5), 2019.

[23] Julien Jacques, Christian Lavergne, and Nicolas Devictor. Sensitivity analysis in presence of model uncertainty and correlated inputs. *Reliability Engineering and System Safety*, 91:1126-1134, 2006.

[24] Christophe Labreuche. Explanation with the Winter value: Efficient computation for hierarchical Choquet integrals. *International Journal of Approximate Reasoning*, 151:225-250, 2022.

[25] Paes Leme and Jon Schneider. Multiparameter bernoulli factories. *The Annals of Applied Probability*, 33(5):3987-4007, 2023.

[26] David E Losada and Leif Azzopardi. Assessing multivariate bernoulli models for information retrieval. *ACM Transactions on Information Systems (TOIS)*, 26(3):1-46, 2008.

[27] Gildas Mazo and Laurent Tournier. An inference method for global sensitivity analysis.
*hal-04199638*, 2023.

[28] Luke Merrick and Ankur Taly. The explanation game: Explaining machine learning models using shapley values. In *Machine Learning and Knowledge Extraction: 4th IFIP TC 5, TC 12, WG 8.4, WG 8.9, WG 12.9 International Cross-Domain Conference, CD-MAKE 2020, Dublin, Ireland, August 25-28, 2020, Proceedings*, page 17-38, Berlin, Heidelberg, 2020. Springer-Verlag.

[29] Albert Myers and Antoine Rauzy. Efficient reliability assessment of redundant systems subject to imperfect fault coverage using binary decision diagrams. *IEEE Transactions on Reliability*, 57(2):336-348, 2008.

[30] Ryan O'Donnell. *Analysis of boolean functions*. Cambridge University Press, 2014.
[31] Art B Owen. Sobol' indices and shapley value. *SIAM/ASA Journal on Uncertainty Quantification*, 2(1):245-251, 2014.

[32] Art B Owen and Clémentine Prieur. On shapley value for measuring importance of dependent inputs. *SIAM/ASA Journal on Uncertainty Quantification*, 5(1):986-1002, 2017.

[33] Gihan Panapitiya, Peiyuan Gao, C Mark Maupin, and Emily G Saldanha. Fragnet: A graph neural network for molecular property prediction with four layers of interpretability. *arXiv:2410.12156*, 2024.

[34] Andrea Ponti, Antonio Irpino, Antonio Candelieri, Anna Bosio, Ilaria Giordani, and Francesco Archetti. Network vulnerability analysis in wasserstein spaces. In Dimitris E. Simos, Varvara A. Rasskazova, Francesco Archetti, Ilias S. Kotsireas, and Panos M. Pardalos, editors, *Learning and Intelligent Optimization*, pages 263-277, Cham, 2022. Springer International Publishing.

[35] Herschel Rabitz and Ömer Faruk Aliş. General foundations of high-dimensional model representations. *Journal of Mathematical Chemistry*, 25(2):197-233, 1999.

[36] Saman Razavi, Anthony J Jakeman, Andrea Saltelli, Clémentine Prieur, Bertrand Iooss, Emanuele Borgonovo, Elmar Plischke, Samuele Lo Piano, Takuya Iwanaga, William Becker, Stefano Tarantola, Joseph Guillaume, John Jakeman, Hoshin Gupta, Nicola Melillo, Giovanni Rabitti, Vincent Chabridon, Qingyun Duan, Xifu Sun, Stefán Smith, Razi Sheikholeslami, Nasim Hosseini, Masoud Asadzadeh, Arnald Puy, Sergei Kucherenko, and Holger Maier. The future of sensitivity analysis: An essential discipline for systems modeling and policy support. *Environmental Modelling & Software*, 137:104954, 2021.

[37] Branko Ristic, Ba-Tuong Vo, Ba-Ngu Vo, and Alfonso Farina. A tutorial on bernoulli filters: theory, implementation and applications. *IEEE Transactions on Signal Processing*, 61(13):3406-3430, 2013.

[38] Han Lin Shang. A survey of functional principal component analysis. *AStA Advances in Statistical Analysis*, 98:121-142, 2014.

[39] Noam Shazeer, Azalia Mirhoseini, Piotr Maziarz, Andy Davis, Quoc Le, Geoffrey Hinton, and Jeff Dean. Outrageously large neural networks: The sparsely-gated mixture-of-experts layer. In *International Conference on Learning Representations (ICLR)*, 2017.

[40] Oleksandr Shchur and Stephan Günnemann. Overlapping community detection with graph neural networks. *Proceedings of ICLR 2019*, abs/1909.12201, 2019.

[41] Ilya Shmulevich, Edward R Dougherty, and Wei Zhang. From boolean to probabilistic boolean networks as models of genetic regulatory networks. *Proceedings of the IEEE*, 90(11):1778-1792, 2002.

[42] Ralph Smith. *Uncertainty quantification: theory, implementation, and applications*. SIAM, 2024.

[43] Eunhye Song, Barry L Nelson, and Jeremy Staum. Shapley effects for global sensitivity analysis: Theory and computation. *SIAM/ASA Journal on Uncertainty Quantification*, 4(1):1060-1083, 2016.

[44] Ingo Steinwart and Clint Scovel. Mercer's theorem on general domains: On the interaction between measures, kernels, and rkhss. *Constructive Approximation*, 35:363-417, 2012.
[45] Ding Wang, Wei Zhou, and Songil Hu. Information diffusion prediction with graph neural ordinary differential equation network. In *Proceedings of the 32nd ACM International Conference on Multimedia*, MM '24, page 9699–9708, New York, NY, USA, 2024. Association for Computing Machinery.

[46] Muhan Zhang and Yixin Chen. Link prediction based on graph neural networks. In *Proceedings of the 32nd International Conference on Neural Information Processing Systems*, NIPS'18, page 5171–5181, Red Hook, NY, USA, 2018. Curran Associates Inc.