Paper deep dive
Barycentric Fused Gromov-Wasserstein Balancing for Causal Inference under Multiple Treatments
Yuki Murakami, Takumi Hattori, Kohsuke Kubota
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%
Last extracted: 8/25/2026, 7:54:38 AM
Summary
The paper introduces CIHSI-Net, a deep learning framework for causal inference under multiple simultaneous treatments. It addresses limitations of existing pairwise balancing methods, such as quadratic computational complexity and inconsistent local proximity structures, by proposing Barycentric Fused Gromov-Wasserstein Balancing (BFG-WB). This objective aligns treatment representation distributions with a shared Wasserstein barycenter using Fused Gromov-Wasserstein discrepancy, achieving global alignment, preserving local geometry, and reducing complexity to linear. The method is validated through simulations and real-world marketing data applications.
Entities (8)
Relation Signals (7)
CIHSI-Net → uses → BFG-WB
confidence 98% · CIHSI-Net, a deep learning framework built on a novel Barycentric Fused Gromov-Wasserstein Balancing (BFG-WB) objective.
BFG-WB → alignsto → Wasserstein Barycenter
confidence 95% · BFG-WB aligns the representation distribution of each treatment pattern with a shared Wasserstein barycenter
BFG-WB → usesmetric → Fused Gromov-Wasserstein Discrepancy
confidence 94% · its Fused Gromov-Wasserstein discrepancy preserves the local proximity structures
BFG-WB → reducescomplexityto → linear complexity
confidence 93% · achieving global alignment while reducing the computational complexity from quadratic to linear
Pairwise Balancing → haslimitation → Quadratic Complexity
confidence 92% · such pairwise balancing scales quadratically with the number of treatment patterns
CIHSI-Net → estimates → CASE
confidence 90% · Simulation studies show that CIHSI-Net consistently outperforms state-of-the-art baselines... in estimating both heterogeneous single and interaction effects.
CIHSI-Net → estimates → CAIE
confidence 90% · Simulation studies show that CIHSI-Net consistently outperforms state-of-the-art baselines... in estimating both heterogeneous single and interaction effects.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Estimating heterogeneous single and interaction treatment effects from observational data under multiple simultaneous treatments is crucial for decision-making. To mitigate estimation variance, previous studies balance representation distributions between every pair of treatment patterns. However, such pairwise balancing scales quadratically with the number of treatment patterns and fails to preserve consistent local proximity structures across patterns, which degrades counterfactual estimation. To address these challenges, we propose the Causal Inference for Heterogeneous Single and Interaction Treatment Effects Network (CIHSI-Net), a deep learning framework built on a novel Barycentric Fused Gromov-Wasserstein Balancing (BFG-WB) objective. BFG-WB aligns the representation distribution of each treatment pattern with a shared Wasserstein barycenter, achieving global alignment while reducing the computational complexity from quadratic to linear, and its Fused Gromov-Wasserstein discrepancy preserves the local proximity structures essential for reliable heterogeneous effect estimation. Simulation studies show that CIHSI-Net consistently outperforms state-of-the-art baselines, and an application to real-world marketing data demonstrates its practical utility in complex multi-treatment scenarios.
Tags
Links
- Source: https://arxiv.org/abs/2608.22024v1
- Canonical: https://arxiv.org/abs/2608.22024v1
Trouble viewing inline? Open PDF directly →
Full Text
138,162 characters extracted from source content.
Expand or collapse full text
Barycentric Fused Gromov-Wasserstein Balancing for Causal Inference under Multiple Treatments Yuki Murakami Takumi Hattori Kohsuke Kubota Abstract Estimating heterogeneous single and interaction treatment effects from observational data under multiple simultaneous treatments is crucial for decision-making. To mitigate estimation variance, previous studies balance representation distributions between every pair of treatment patterns. However, such pairwise balancing scales quadratically with the number of treatment patterns and fails to preserve consistent local proximity structures across patterns, which degrades counterfactual estimation. To address these challenges, we propose the Causal Inference for Heterogeneous Single and Interaction Treatment Effects Network (CIHSI-Net), a deep learning framework built on a novel Barycentric Fused Gromov-Wasserstein Balancing (BFG-WB) objective. BFG-WB aligns the representation distribution of each treatment pattern with a shared Wasserstein barycenter, achieving global alignment while reducing the computational complexity from quadratic to linear, and its Fused Gromov-Wasserstein discrepancy preserves the local proximity structures essential for reliable heterogeneous effect estimation. Simulation studies show that CIHSI-Net consistently outperforms state-of-the-art baselines, and an application to real-world marketing data demonstrates its practical utility in complex multi-treatment scenarios. 1 Introduction Estimating heterogeneous single and interaction treatment effects under multiple simultaneous treatments is important in many domains, such as healthcare and marketing. For example, optimizing combination drug therapies requires understanding how these effects vary due to physiological differences (21; 46, e.g.). Similarly, in marketing, promotion effectiveness depends heavily on user attributes and their interactions (4; 29, e.g.). Given that randomized trials are often infeasible due to ethical concerns and the high costs associated with the combinatorial explosion of treatment options, analytical frameworks that can accurately estimate these effects from observational data are essential for decision-making in real-world applications. To enable such accurate estimation from observational data, representation balancing with deep learning has emerged as a promising approach to address the inherent selection bias. This approach mitigates distributional discrepancies between treatment and control groups by learning a balanced representation (47; 58, e.g.). Recently, Causal Inference for Single and Interaction Treatment Effects Network (CISI-Net) extended this approach to the multiple-treatment setting (36). To handle multiple treatments, this approach adopts a pairwise balancing strategy that aligns distributions between specific pairs of treatment patterns, while integrating a task embedding network to encode the unique contributions of single and interaction effects. However, existing pairwise balancing methods face three critical limitations. First, they struggle to achieve global alignment because reducing the discrepancy for one treatment-pattern pair may increase the discrepancies between other treatment-pattern pairs, leaving residual imbalance that can increase estimation variance. Second, pairwise constraints generally do not encourage consistent local proximity structures across all treatment patterns, and such inconsistent geometric distortions can degrade heterogeneous effect estimation (6; 54; 58). Third, the quadratic computational complexity from matching all pairs makes these methods impractical for large-scale applications with numerous treatment combinations. To address these challenges, we propose the Causal Inference for Heterogeneous Single and Interaction Treatment Effects Network (CIHSI-Net), which minimizes global distributional discrepancies while preserving local proximity structures. Central to CIHSI-Net is a novel Barycentric Fused Gromov-Wasserstein Balancing (BFG-WB) objective, which aligns the representation distributions of each treatment pattern to a shared Wasserstein barycenter (2) using the Fused Gromov-Wasserstein (FGW) discrepancy (50; 52). This approach achieves global alignment and preserves local proximity structures by accounting for both feature values and geometry, while reducing computational complexity from quadratic to linear. Simulation results demonstrate that CIHSI-Net consistently outperforms state-of-the-art baselines in estimating both heterogeneous single and interaction effects. Furthermore, ablation studies confirm the critical contribution of BFG-WB, and application to real-world marketing data validates the practical utility of the proposed framework. 2 Related Work We review existing research from two perspectives: representation balancing methods and deep learning-based methods for multiple treatments. First, deep learning-based representation balancing mitigates estimation variance arising from selection bias in observational data by reducing the discrepancy between the representations of the treatment and control groups. Optimal Transport (OT)-based regularization, which globally aligns the representation distributions across treatment groups, has proven effective for estimating the heterogeneous causal effect (47; 31). Motivated by the reliance of counterfactual inference on local smoothness, recent extensions encourage the preservation of local proximity structures within each distribution by keeping input-space neighbors adjacent after balancing (58; 59) or by guiding the cross-group alignment using local neighborhood geometries (54; 6). However, extending this pairwise balancing paradigm to the multiple-treatment setting faces fundamental limitations. For global alignment, balancing each pair of treatment patterns separately (36) allows the alignment of one pair to degrade that of another, leaving a residual imbalance (20; 32), with a cost quadratic in the number of patterns. For local structure preservation, pairwise objectives promote geometric consistency only within pairs, not across all treatment patterns (33; 35). Consequently, in a multiple-treatment setting, a unified framework achieving scalable global balancing with local geometry preservation remains an open challenge. Second, instead of relying on pseudo-samples via data augmentation (41, e.g.), we focus on architecture-driven approaches that intrinsically capture treatment interactions. These methods typically rely on either separate outcome networks (39, e.g.) or latent-variable generative models (44; 60, e.g.). However, separate networks often yield unstable estimates for rare treatment patterns due to isolated parameters (11), while latent-variable methods suffer reduced robustness when generative assumptions fail (42). Consequently, there remains a need for an architecture that achieves efficient parameter sharing without strong generative assumptions. To fulfill these needs, we propose a novel framework, CIHSI-Net. The following section defines the problem setting and the causal estimands necessary for our framework. 3 Preliminaries Heterogeneous Effects under Multiple Treatments We formulate the problem of estimating heterogeneous single and interaction treatment effects under multiple treatments with the potential outcomes framework (43). Let =(i,i,yi)i=1ND=\( x_i, t_i,y_i)\_i=1^N denote an observed dataset of N independent units. Here, ∈⊆ℝd X ^d denote the covariate vector and T∈=0,1KT =\0,1\^K denote the vector of K simultaneous binary treatments, with realizations x and t. For each unit i, a potential outcome Yi()∈ℝY_i( t) exists for every ∈ t , but only the factual outcome yi=Yi(i)y_i=Y_i( t_i) is observed. Let Y()Y( t) and Y represent the corresponding population-level random variables for the potential and factual outcomes, respectively. To identify causal effects of interest from the observed data, we adopt the following three assumptions commonly used in observational studies (23). Assumption 1 (Stable Unit Treatment Value Assumption). (1) no interference, meaning that the outcome of one unit is unaffected by the treatment assignments of other units; and (2) consistency of treatment, meaning that the potential outcomes correspond to well-defined and unique treatments (i.e., yi=Yi(i)y_i=Y_i( t_i)). Assumption 2 (Ignorability). For any treatment pattern, the potential outcome is independent of the assigned treatment T given the observed covariates X. Formally, for all t, Y()⟂|.Y( t) T X. Assumption 3 (Overlap). Every unit has a non-zero probability of receiving any treatment pattern given its observed covariates. Formally, for all t and x, 0<P(=∣=)<1.0<P( T= t X= x)<1. Based on these assumptions, we define our estimands of interest: the Conditional Average Single Effect (CASE) and the Conditional Average Interaction Effect (CAIE). Both estimands are defined using the conditional expected potential outcome given by μ(,)≔[Y()∣=].μ( x, t) [Y( t) X= x ]. (1) First, the CASE for treatment k quantifies the marginal effect of applying treatment k in isolation compared to the control group (no treatment). This estimand extends the standard conditional average treatment effect in the single-treatment setting (1) to the multi-treatment setting. The CASE is defined as τCASE(k,)≔μ(,+k)−μ(,), _CASE(k, x) μ( x, t_+k)-μ( x, 0), (2) where +k t_+k is the one-hot vector for treatment k. Second, the CAIE quantifies interactions among a subset of treatments S⊆1,…,KS \1,…,K\ (where |S|≥2|S|≥ 2). We adopt the additive scale, which aligns with decision-making based on absolute outcome differences (e.g., revenue), and extend the Average Interaction Effect (AIE) (17) to the conditional setting as follows. τCAIE(S,)≔∑Q⊆S(−1)|S|−|Q|μ(,(+Q)), _CAIE(S, x) _Q S(-1)^|S|-|Q|\,μ ( x, t_(+Q) ), (3) where (+Q) t_(+Q) denotes a vector with ones at indices in Q and zeros elsewhere. For example, in a two-treatment setting, τCAIE(1,2,) _CAIE(\1,2\, x) measures the deviation of the joint effect at =(1,1) t=(1,1) from the sum of the single effects of (1,0)(1,0) and (0,1)(0,1). We now establish that the estimands defined above are identifiable from the joint distribution of the observed data (,,y)( x, t,y) under the stated assumptions. Proposition 1. Under Assumptions 1-3, CASE and CAIE are identifiable from the observed data. The proof is provided in Appendix A.1. Motivation for Representation Balancing Figure 1: Illustration of neighborhood switching under Wasserstein-based representation balancing. The gray arrows show OT correspondences. Wasserstein minimization aligns the two distributions but can map neighboring points to distant locations, distorting within-distribution neighborhood relations; a structure-preserving match such as j↔j′j\! \!j would be desirable. Representation balancing mitigates estimation variance by minimizing distributional discrepancies across treatment groups. A widely adopted metric is the Wasserstein distance (8; 47, e.g.), which uses OT to measure the minimal transformation cost between distributions. However, because the Wasserstein distance minimizes only the aggregate point-to-point transport cost, it does not explicitly account for local proximity structures within each distribution (54). Here, the local proximity structures refer to the relative geometric relationships among units within a specific treatment pattern. Consequently, as illustrated in Figure 1, neighbors in one distribution may be mapped to scattered locations in another, potentially leading to structural collapse. This disruption of local proximity structures degrades the accuracy of causal effect estimation. Counterfactual prediction typically relies on the assumption of local smoothness, where units with similar representations yield similar outcomes (30; 55). Therefore, accurate estimation of counterfactual outcomes Y()Y( t) benefits from (i) reducing distributional discrepancies (global alignment) and (i) ensuring that neighborhood consistency is maintained across treatment patterns (local structure preservation). This issue is particularly critical in the multiple-treatment setting. When pairwise balancing is applied across numerous treatment patterns, local structure preservation becomes inconsistent, where the structural distortion may vary significantly depending on the pair being aligned. Such non-uniformity destabilizes the reference neighbors used for inference, thereby degrading the estimation of heterogeneous single and interaction effects. To address these challenges, we propose a novel framework that simultaneously achieves global alignment and local structure preservation. We introduce a new architecture incorporating the FGW discrepancy to satisfy these dual objectives in the following section. 4 CIHSI-Net: The Proposed Framework Barycentric Fused Gromov-Wasserstein Balancing To overcome the limitations of pairwise balancing, we propose BFG-WB. BFG-WB functions as a regularization term within our proposed CIHSI-Net, designed to achieve two critical objectives simultaneously: (i) reducing estimation variance induced by distributional imbalances globally by aligning all treatment representations to a common anchor, and (i) preserving local proximity structures to ensure stable heterogeneous effect estimation. Formally, BFG-WB minimizes the weighted sum of discrepancies between the representation distribution of each treatment pattern t and a shared Wasserstein barycenter b∗ R _b. Unlike standard approaches that use the Wasserstein distance, we employ the FGW discrepancy (50; 52). By integrating the feature-based Wasserstein distance with the structure-based Gromov-Wasserstein distance, the FGW discrepancy enables us to simultaneously achieve global distribution alignment (Objective (i)) and local structure preservation (Objective (i)). The BFG-WB regularization term ℒϕL_φ is defined as follows. ℒϕ=∑∈wF(,b∗),L_φ= _ t w_ tF ( R_ t, R _b ), (4) where b∗ R_b is the Wasserstein barycenter, and R_ t denotes the representation distribution on ℛR for treatment pattern t. F(⋅,⋅)F(·,·) measures the FGW discrepancy between two distributions, and w>0w_ t>0 is a weight satisfying ∑∈w=1 _ t w_ t=1. In the following, we detail the two core components: the Wasserstein barycenter b∗ R_b and the FGW discrepancy F(⋅,⋅)F(·,·). Wasserstein Barycenter. Following the formulation in (2), we define the Wasserstein barycenter b∗ R_b as the centroid distribution that minimizes the weighted sum of squared Wasserstein distances from the representation distributions R_ t as follows. b∗=argminb∑∈λW22(,b), R _b= R_bargmin _ t _ tW^2_2( R_ t, R_b), (5) where W2(⋅,⋅)W_2(·,·) is the 2-Wasserstein distance (see Appendix A.3.1 for its definition). b R_b denotes a barycenter candidate distribution, and λ≥0 _ t≥ 0 represents the weight for treatment pattern t satisfying ∑∈λ=1 _ t _ t=1. Crucially, the introduction of the Wasserstein barycenter b∗ R_b resolves the optimization conflicts in pairwise balancing. In standard pairwise approaches, minimizing the discrepancy for one pair of treatments may inadvertently increase it for another due to differing OT plans, which leads to persistent local imbalances associated with selection bias. In contrast, our approach adopts a “star-shaped” (See Figure 2) alignment strategy where all representation distributions are updated toward a single, fixed reference (barycenter) b∗ R_b . This alignment strategy promotes coherent alignment and effectively mitigates estimation variance across all treatment patterns (See Appendix A.3.2). Figure 2: Conceptual illustration of BFG-WB. The gray arrows show a subset of OT correspondences. By evaluating the FGW discrepancy, the alignment is encouraged to preserve within-distribution neighborhood relations. Fused Gromov-Wasserstein Discrepancy. With the barycenter b∗ R_b established as the global anchor, we employ the FGW discrepancy to measure the discrepancy between the representation distribution R_ t of each treatment pattern t and b∗ R_b . The FGW discrepancy integrates the Wasserstein distance (for global distribution alignment), which evaluates feature values, and the Gromov-Wasserstein distance (for local proximity structure preservation), which evaluates the geometric structure within distributions, under a single OT plan π()π^( t). Specifically, the FGW discrepancy F(,b∗)F( R_ t, R_b ) is defined as follows. F(,∗b)=infπ()∈Π(,b∗)[η∫c()(,)dπ()(,) F\! ( R_ t, R _b )= _π^( t)∈ ( R_ t, R _b) [η c^( t)\! ( r, z )\,dπ^( t)( r, z) (6) +(1−η)∬db()(,′,,′)dπ()(,)dπ()(′,′)], +\;(1-η) d_b^( t)( r, r , z, z )\,dπ^( t)( r, z)\,dπ^( t)( r , z ) ], where Π(,b∗) ( R_ t, R _b) denotes the set of all couplings on ℛ×ℛR×R with marginals R_ t and b∗ R _b. Here, ,′,,′∈ℛ r, r , z, z with ,′∼ r, r \! \! R_ t and ,′∼b∗ z, z \! \! R _b. The first term accounts for global distribution alignment, where c()(,)=‖−‖2c^( t)( r, z)=\| r- z\|_2 measures the Euclidean distance between features. The second term accounts for local proximity structure preservation. Specifically, we define db()(,′,,′):=(d()(,′)−d(b)(,′))2d_b^( t)( r, r , z, z ):= (d^( t)( r, r )-d^(b)( z, z ) )^2, where d()(,′)=‖−′‖2d^( t)( r, r )=\| r- r \|_2 and d(b)(,′)=‖−′‖2d^(b)( z, z )=\| z- z \|_2 represent the within-distribution distances (local geometry) within the treatment representation distribution and the barycenter, respectively. The hyperparameter η∈(0,1]η∈(0,1] controls the trade-off between the first term (global distribution alignment) and the second term (local proximity structure preservation). The second term in Equation (6) plays an important role in our framework. This term incurs a low cost only when pairs of points that are close in the representation space of treatment t are matched to pairs that are similarly close in the barycenter, thereby penalizing structure-blind matching based solely on feature similarity. This mechanism mitigates distortions in neighborhood geometry by aligning each distribution to the local proximity structure of the barycenter. As a result, it prevents counterfactual prediction from being performed using neighborhoods that mix units with different heterogeneity patterns. This consistency significantly improves the estimation stability of heterogeneous causal effects, because such estimation depends on covariate-conditioned counterfactual prediction. Beyond its methodological advantages in estimation variance reduction and structure preservation, BFG-WB offers a significant practical benefit: computational efficiency. Whereas pairwise balancing requires comparing all pairs of treatment patterns, our star-shaped strategy compares each pattern only with the barycenter. We formalize this advantage in the following proposition. Proposition 2 (Computational Efficiency). Let L=2KL=2^K denote the number of treatment patterns derived from K binary treatments, and let one evaluation denote the computation of a single OT-based discrepancy between two empirical distributions. Per training step, standard pairwise balancing requires (L2)=O(L2) L2=O(L^2) evaluations, whereas BFG-WB requires at most (Lb+1)L=O(L)(L_b+1)\,L=O(L) evaluations, where LbL_b denotes a fixed upper bound on the number of barycenter-update iterations and is independent of L. A detailed derivation is provided in Appendix A.2, and we empirically verify this scalability advantage through runtime comparisons in Appendix D.1 Theoretical Analysis In this section, we demonstrate that BFG-WB provides theoretically sound error control for estimating both CASE and CAIE from observational data. Although accurate causal effect estimation requires controlling for prediction errors under all possible treatment patterns, counterfactual prediction errors cannot be minimized directly because only factual outcomes are observed. To bridge this gap, we draw upon the theoretical framework of representation-based causal inference (54). A key insight from this literature is that the prediction error on the unobserved target domain (counterfactuals) is theoretically bounded by the error on the observed source domain (factuals) plus a discrepancy measure between their distributions. Motivated by this principle, we derive upper bounds for the estimation errors of CASE and CAIE that are expressible solely through observational quantities. Specifically, we demonstrate that these errors are bounded by the sum of the factual prediction error and the distributional discrepancies, which are explicitly minimized by the proposed BFG-WB. The accuracy of the estimated CASE and CAIE is assessed using the integrated squared error with respect to the marginal covariate distribution p()p( x). These metrics quantify the expected squared deviation between the estimated and true causal effects, defined analogously to the Precision in Estimation of Heterogeneous Effects (PEHE) (22) used in single-treatment studies. Formally, the estimation errors of CASE and CAIE are defined as follows. ϵCASE(k)=∫(τ^CASE(k,)−τCASE(k,))2p(), _CASE(k)= _X ( τ_CASE (k, x )- _CASE (k, x ) )^2p( x)d x, (7) ϵCAIE(S)=∫(τ^CAIE(S,)−τCAIE(S,))2p(), _CAIE(S)= _X ( τ_CAIE (S, x )- _CAIE (S, x ) )^2p( x)d x, (8) where τ^CASE(⋅) τ_CASE(·) and τ^CAIE(⋅) τ_CAIE(·) denote the estimated CASE and CAIE. For each treatment pattern t, the expected prediction error of the factual outcome is defined as follows. ϵF()=∫l(,)p(∣)p(), _F^( t)= _Xl ( x, t )\,p( x t)\,p( t)\,d x, (9) where l(,)=∫(Y()−μ^(,))2p(Y()∣)Y()l( x, t)= (Y( t)- μ( x, t) )^2\,p (Y( t) x )\,dY( t) denotes the conditional expected squared prediction loss, and μ^(,) μ( x, t) denotes the estimate of μ(,)μ( x, t). Formally, we establish the upper bounds for ϵCASE _CASE and ϵCAIE _CAIE as follows. Theorem 1 (Upper Bound for CASE). Suppose that Assumptions 1–3 and the auxiliary conditions in Appendix A.3.1 hold (in particular, the expected squared loss l(,)l( x, t) is BϕB_φ-Lipschitz in the representation space), and that the weights in Equation (4) are uniform, i.e., w=2−Kw_ t=2^-K. Then, for any k∈1,…,Kk∈\1,…,K\, ϵCASE(k)≤21p(+k)ϵF(+k)+1p()ϵF()+22KηBϕℒϕ. _CASE(k)≤ 2 \ 1p( t_+k) _F^( t_+k)+ 1p( 0) _F^( 0)+ 2^2KηB_φL_φ \. Theorem 2 (Upper Bound for CAIE). Under the same assumptions as Theorem 1, consider an interaction set S⊆1,…,KS \1,…,K\ with |S|≥2|S|≥ 2. Let (a)∈(a_ t)_ t be a constant vector reflecting the combinatorial structure of S (a∈−1,0,1a_ t∈\-1,0,1\). The estimation error ϵCAIE(S) _CAIE(S) is bounded by: ϵCAIE(S) _CAIE(S) ≤(∑∈a2)∑∈1p()ϵF() ≤ ( _ t a_ t^2 ) \ _ t 1p( t)\, _F^( t) +2K+1ηBϕ(2K−1)ℒϕ. + 2^K+1η\,B_φ\, (2^K-1 )\,L_φ \. The detailed proofs for Theorem 1 and Theorem 2 are provided in Appendix A.3. These theorems demonstrate that the estimation errors for CASE and CAIE are theoretically upper-bounded by the sum of the factual prediction error ϵF() _F^( t) and the distributional discrepancy term ℒϕL_φ. Crucially, these bounds explicitly highlight the inverse probability term 1/p()1/p( t), which formalizes the vulnerability of causal inference to rare treatment patterns. Furthermore, the dependence on 2K2^K in these bounds reflects the inherent combinatorial complexity of estimating causal effects over all treatment patterns. By aligning all treatment-pattern representation distributions to a single fixed anchor b∗ R_b , our star-shaped strategy efficiently minimizes all pairwise Wasserstein discrepancies via the triangle inequality. Simultaneously, the Gromov-Wasserstein component of the FGW discrepancy preserves local proximity structures within the representation space. This structural preservation makes the Lipschitz continuity assumption (BϕB_φ) more plausible in practice, supporting counterfactual inference from nearby factual samples. Together, these results show that BFG-WB provides both computational efficiency and a principled mechanism for controlling the estimation errors of CASE and CAIE. Network Architecture and Training Objective We implement our proposed method using CIHSI-Net, a neural network architecture designed to achieve both accurate outcome prediction and robust representation balancing. As illustrated in Figure 3, the network consists of three components: a representation learning network ϕ:ℝd→ℝdrφ:R^d ^d_r regularized by BFG-WB, a task embedding network tw:ℝK→ℝdtt_w:R^K ^d_t that captures treatment similarities, and the outcome prediction network h:ℝdr+dt→ℝh:R^d_r+d_t that estimates outcomes from the concatenated features. Figure 3: The architecture of CIHSI-Net. The model consists of three components: the representation learning network ϕφ, the task embedding network twt_w, and the outcome prediction network h. BFG-WB is a regularization term applied to the latent representation ϕ()φ( x). The entire network is optimized by minimizing a composite loss function ℒL that balances factual prediction accuracy with the distributional alignment constraints derived in our theoretical analysis as follows. ℒ=ℒy+αℒϕ+β‖w‖22,L=L_y+ _φ+β\|w\|_2^2, (10) where the first term ℒyL_y represents the factual outcome prediction error, the second term ℒϕL_φ is the BFG-WB regularization term, and the third term is an L2 regularization term applied to the network weights. The coefficients α and β are hyperparameters that control the strength of the corresponding terms. The outcome prediction loss ℒyL_y is designed to align with the error bounds established in Theorems 1 and 2. Guided by our theoretical analysis, which reveals that prediction errors on rare treatment patterns are amplified by the inverse probability 1/p()1/p( t), we employ an inverse-frequency weighted mean squared error to counteract this imbalance as follows. ℒy=1N∑i=1N1p^(i)(yi−y^i)2,L_y= 1N _i=1^N 1 p( t_i)(y_i- y_i)^2, (11) where p^(i)=1N∑j=1N[i=j] p( t_i)= 1N _j=1^NI[ t_i= t_j] represents the empirical frequency of the treatment pattern assigned to unit i. 5 Simulation Experiments We evaluate our CIHSI-Net using simulation datasets to (i) compare estimation accuracy with state-of-the-art baselines and (i) clarify the contributions of its components through ablation studies. Datasets. To assess estimation accuracy under selection bias and heterogeneous effects, we generate two types of simulation datasets following the protocol in 36. Simulation 1 (With Interactions) is generated with complex interaction effects among treatments and serves as the primary benchmark for evaluating CASE and CAIE estimation. Simulation 2 (Without Interactions) is generated without interaction effects and tests the model’s robustness by verifying that it does not produce spurious interaction effects when none exist. We fix the number of treatments at K=3K=3 and the sample size at N=50,000N=50,000. Detailed data generation processes are provided in Appendix B. Experimental Setup. We compare CIHSI-Net with three representative multi-treatment causal inference methods: Task Embedding–based Causal Effect Variational Autoencoder (TECE-VAE) (44), Neural Counterfactual Relation Estimation (NCoRE) (39), and CISI-Net (36). Implementation details for all baselines and the hyperparameters of CIHSI-Net are provided in Appendix C. The same network hyperparameter settings are applied across all experiments. For CIHSI-Net, we set α in Equation (10) to 1.0 and η in Equation (6) to 0.6 for comparative experiments, while fixing α=1.0α=1.0 and varying η for ablation studies. All models, including baselines, are trained using standard optimization protocols with 60% training, 10% validation, and 30% test data. A sensitivity analysis regarding the regularization coefficients is presented in Appendix D.5. Evaluation Metrics. Performance is evaluated using the integrated squared errors for CASE and CAIE, as defined in Equations (7) and (8). To ensure statistical reliability, we report average metrics over 100 independent runs. Accuracy Comparison against Baselines. Table 1 shows ϵCASE _CASE and ϵCAIE _CAIE of the proposed and baseline methods across two simulation datasets. Across all simulation settings, the proposed method consistently achieves the lowest estimation errors for both CASE and CAIE, outperforming all baseline methods. CIHSI-Net maintains high accuracy regardless of whether interaction effects are present, and its advantage over the baselines is most pronounced for ϵCAIE _CAIE, indicating that BFG-WB is particularly beneficial for interaction treatment effect estimation. Further analysis in Appendix D.4 demonstrates that our CIHSI-Net effectively reduces distributional discrepancies (measured by Wasserstein distance) while preserving local proximity structures (measured by Gromov-Wasserstein distance), thereby contributing to superior estimation performance. Moreover, CIHSI-Net maintains its CASE and CAIE estimation advantage as K increases to eight, demonstrating its scalability (see Appendix D.2), and achieves the best overall estimation performance in the semi-synthetic setting (see Appendix D.6). ϵCASE _CASE ϵCAIE _CAIE Sim. Method k=1k=1 k=2k=2 k=3k=3 S=1,2S=\1,2\ S=2,3S=\2,3\ S=1,3S=\1,3\ S=1,2,3S=\1,2,3\ 1 TECE-VAE 1.82 ± 0.09 1.91 ± 0.10 1.80 ± 0.14 3.58 ± 0.11 3.27 ± 0.10 3.30 ± 0.10 6.80 ± 0.18 NCoRE 0.28 ± 0.06 0.30 ± 0.06 0.25 ± 0.05 0.52 ± 0.10 0.33 ± 0.10 0.37 ± 0.10 0.56 ± 0.14 CISI-Net 0.22 ± 0.04 0.24 ± 0.03 0.19 ± 0.06 0.31 ± 0.07 0.11 ± 0.04 0.15 ± 0.08 0.31 ± 0.08 CIHSI-Net 0.19 ± 0.04 0.21 ± 0.04 0.18 ± 0.10 0.28 ± 0.05 0.08 ± 0.05 0.12 ± 0.04 0.24 ± 0.04 2 TECE-VAE 1.60 ± 0.05 1.65 ± 0.09 1.58 ± 0.16 3.01 ± 0.05 2.65 ± 0.08 2.69 ± 0.09 5.50 ± 0.15 NCoRE 0.26 ± 0.05 0.28 ± 0.09 0.24 ± 0.07 0.34 ± 0.12 0.30 ± 0.09 0.27 ± 0.10 0.35 ± 0.12 CISI-Net 0.16 ± 0.04 0.19 ± 0.07 0.16 ± 0.05 0.08 ± 0.08 0.07 ± 0.04 0.07 ± 0.05 0.11 ± 0.09 CIHSI-Net 0.14 ± 0.04 0.17 ± 0.03 0.16 ± 0.12 0.06 ± 0.04 0.06 ± 0.03 0.06 ± 0.03 0.09 ± 0.04 Table 1: Comparison of mean estimation errors and standard deviations for ϵCASE _CASE and ϵCAIE _CAIE in two simulation datasets. Here, k∈1,2,3k∈\1,2,3\ indexes individual treatments, and S⊆1,2,3S \1,2,3\ (with |S|≥2|S|≥ 2) denotes treatment combinations. Ablation Study. To verify the contribution of each component in our proposed method, we conducted an ablation study by varying the alignment strategy (with or without Wasserstein barycenter) and the discrepancy metric (Wasserstein distance (W), Gromov-Wasserstein discrepancy (GW), and FGW discrepancy (FGW)). Table 2 summarizes the results. Detailed results for ϵCASE _CASE and ϵCAIE _CAIE are reported in Appendix D.3. First, the combination of the Wasserstein barycenter and FGW discrepancy (No.7) consistently achieves the lowest average estimation errors for both CASE and CAIE. Second, methods using the Wasserstein barycenter generally outperform their pairwise counterparts (compare No.5-7 with No.2-4). This finding suggests that pairwise discrepancy minimization is insufficient for globally reducing representation imbalance across treatment patterns, whereas barycentric alignment more effectively mitigates imbalance-driven estimation variance. Third, regarding the discrepancy metric, while the FGW discrepancy fails to improve or even yields slightly inferior results compared to the Wasserstein distance in the pairwise setting (No.4 vs No.2), it improves performance when combined with the Wasserstein barycenter (No.7 vs No.5). These findings suggest that the structural preservation capability of FGW discrepancy is effective only when supported by the global alignment provided by the Wasserstein barycenter. No. WB Penalty type Average ϵCASE _CASE Average ϵCAIE _CAIE 1 ✗ None 0.27 ± 0.06 0.24 ± 0.02 2 ✗ W 0.22 ± 0.03 0.22 ± 0.03 3 ✗ GW 0.25 ± 0.06 0.24 ± 0.03 4 ✗ FGW 0.22 ± 0.02 0.23 ± 0.02 5 ✓ W 0.22 ± 0.04 0.21 ± 0.02 6 ✓ GW 0.22 ± 0.04 0.22 ± 0.03 7 ✓ FGW 0.19 ± 0.04 0.18 ± 0.02 Table 2: Ablation study results on simulation dataset 1. Average ϵCASE _CASE and Average ϵCAIE _CAIE denote the means of the reported ϵCASE _CASE and ϵCAIE _CAIE values, respectively. WB indicates whether the Wasserstein barycenter is used. Penalties: Wasserstein (W), Gromov-Wasserstein (GW), Fused GW (FGW), and None. 6 Application to Multiple Marketing Promotions We apply CIHSI-Net to a real-world dataset from a mobile payment platform to validate its effectiveness in capturing heterogeneous causal effects. The dataset includes three simultaneously conducted promotions: two offline promotions by the same merchant group (CP1, CP2) and one online promotion by a different merchant group (CP3). The outcome is defined as the total payment amount during the promotion period, standardized for confidentiality. User covariates include 71 features derived from service usage history and demographics. To analyze heterogeneity, we stratify users into 11 groups based on their payment amount in the pre-promotion month. Hyperparameters are set to the best-performing configuration from the simulation study (Section 5). Detailed dataset statistics and preprocessing procedures are provided in Appendix E.1. Figure 4 shows the estimated CASE and CAIE across user groups obtained from a single algorithm run. Figure 4: Estimated CASE (top) and CAIE (bottom) across user groups stratified by pre-promotion service usage. The top panel in Figure 4 shows that all estimated CASEs are positive. The online promotion CP3 exhibits the strongest effect for low-usage users. This result is consistent with findings that online incentives effectively encourage increased service spending, even among users with little service experience (4; 27). In contrast, the offline promotions (CP1, CP2) show modest increases for higher-usage groups. These results suggest that user groups for whom effects are strongly observed may differ across channels and promotion designs. The bottom panel in Figure 4 highlights the capability of CIHSI-Net to capture complex interaction structures. Whereas the estimated interaction between same-group promotions (τCAIE(1,2) _CAIE(\1,2\)) remains positive but declines with usage, cross-group interactions show mixed patterns. Notably, the estimated three-way interaction τCAIE(1,2,3) _CAIE(\1,2,3\) shifts from negative for low-usage users to positive for high-usage users. This result suggests that, for users who regularly use the service, conducting more promotions simultaneously may yield complementary effects, whereas for low-usage users, presenting multiple incentives may divide attention (like choice overload (9)) and result in a negative interaction effect. Additional sensitivity and uncertainty analyses in Appendix E.2 show that the qualitative CASE and CAIE patterns remain stable across the BFG-WB hyperparameters α and η and are broadly consistent under bootstrap confidence intervals. These results highlight the importance of capturing treatment effect heterogeneity and demonstrate the practical utility of CIHSI-Net in revealing complex behavioral mechanisms, such as coexisting synergistic and cannibalistic interactions that vary substantially across user groups. To further validate the generalizability of our CIHSI-Net, we extended our analysis to a second real-world scenario involving promotions from competing merchants. Due to space constraints, detailed results demonstrating the detection of cannibalization are provided in Appendix E.3. 7 Conclusion This study addresses the limitations of pairwise balancing in multi-treatment causal inference by proposing Causal Inference for Heterogeneous Single and Interaction Treatment Effects Network (CIHSI-Net). By introducing Barycentric Fused Gromov-Wasserstein Balancing (BFG-WB), our framework achieves scalable global alignment via a Wasserstein barycenter while preserving local proximity structures through the FGW discrepancy. Experiments on simulation datasets demonstrate that CIHSI-Net outperforms baselines in estimating heterogeneous single and interaction effects, and a real-world marketing application illustrates its practical utility in multi-treatment settings. Several directions remain for future work. First, a promising extension is to replace the Wasserstein barycenter with an FGW barycenter to unify global alignment and structural preservation, enabling more consistent alignment of feature distributions and within-distribution geometry. Second, it will be valuable to connect CIHSI-Net to downstream decision-making, such as uplift-based allocation for combinatorial treatment. References Abrevaya et al. (2015) J. Abrevaya, Y. Hsu, and R. P. Lieli Estimating conditional average treatment effects. Journal of Business & Economic Statistics 33 (4), p. 485–505. External Links: Document, Link, https://doi.org/10.1080/07350015.2014.975555 Cited by: §3. Agueh and Carlier (2011) M. Agueh and G. Carlier Barycenters in the wasserstein space. SIAM Journal on Mathematical Analysis 43 (2), p. 904–924. External Links: Document, Link, https://doi.org/10.1137/100805741 Cited by: §1, §4. Ascarza (2018) E. Ascarza Retention futility: targeting high-risk customers might be ineffective. Journal of Marketing Research 55 (1), p. 80–98. External Links: Document, Link, https://doi.org/10.1509/jmr.16.0163 Cited by: Appendix B. Blake et al. (2015) T. Blake, C. Nosko, and S. Tadelis Consumer heterogeneity and paid search effectiveness: a large-scale field experiment. Econometrica 83 (1), p. 155–174. Cited by: §1, §6. Breugelmans and Liu-Thompkins (2017) E. Breugelmans and Y. Liu-Thompkins The effect of loyalty program expiration policy on consumer behavior. Marketing Letters 28 (4), p. 537–550. Cited by: §E.3. Cao et al. (2026) F. Cao, J. Zhang, and X. Li PITE: multi-prototype alignment for individual treatment effect estimation. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 40, p. 19871–19879. Cited by: §1, §2. Chauhan et al. (2025) V. K. Chauhan, L. Clifton, G. Nigam, and D. A. Clifton Individualised treatment effects estimation with composite treatments and composite outcomes. In 2025 47th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), p. 1–7. External Links: Link, Document Cited by: 3rd item. Cheng et al. (2022) M. Cheng, X. Liao, Q. Liu, B. Ma, J. Xu, and B. Zheng Learning disentangled representations for counterfactual regression via mutual information minimization. In Proceedings of the 45th International ACM SIGIR Conference on Research and Development in Information Retrieval, SIGIR ’22, New York, NY, USA, p. 1802–1806. External Links: ISBN 9781450387323, Link, Document Cited by: §3. Chernev et al. (2015) A. Chernev, U. Böckenholt, and J. Goodman Choice overload: a conceptual review and meta-analysis. Journal of Consumer Psychology 25 (2), p. 333–358. External Links: ISSN 1057-7408, Document, Link Cited by: §6. Cho et al. (2022) H. J. Cho, J. Chae, S. Yoon, and D. Kim Aging and the prevalence of polypharmacy and hyper-polypharmacy among older adults in south korea: a national retrospective study during 2010–2019. Frontiers in Pharmacology Volume 13 - 2022. External Links: Link, Document, ISSN 1663-9812 Cited by: §D.2. Chu et al. (2022) Z. Chu, H. Ding, G. Zeng, Y. Huang, T. Yan, Y. Kang, and S. Li Hierarchical capsule prediction network for marketing campaigns effect. In Proceedings of the 31st ACM International Conference on Information & Knowledge Management, CIKM ’22, New York, NY, USA, p. 3043–3051. External Links: ISBN 9781450392365, Link, Document Cited by: §2. Clevert et al. (2016) D. Clevert, T. Unterthiner, and S. Hochreiter Fast and accurate deep network learning by exponential linear units (elus). External Links: 1511.07289, Link Cited by: 1st item. Cuturi and Doucet (2014) M. Cuturi and A. Doucet Fast computation of wasserstein barycenters. In Proceedings of the 31st International Conference on Machine Learning, E. P. Xing and T. Jebara (Eds.), Proceedings of Machine Learning Research, Vol. 32, Bejing, China, p. 685–693. External Links: Link Cited by: §A.2. Djulbegovic et al. (2014) B. Djulbegovic, S. Elqayam, T. Reljic, I. Hozo, B. Miladinovic, A. Tsalatsanis, A. Kumar, J. Beckstead, S. Taylor, and J. Cannon-Bowers How do physicians decide to treat: an empirical evaluation of the threshold model. BMC Medical Informatics and Decision Making 14 (1), p. 47. External Links: Document, Link, ISSN 1472-6947 Cited by: Appendix B. Dorotic et al. (2021) M. Dorotic, D. Fok, P. C. Verhoef, and T. H. A. Bijmolt Synergistic and cannibalization effects in a partnership loyalty program. Journal of the Academy of Marketing Science 49 (5), p. 1021–1042. External Links: ISSN 1552-7824, Document, Link Cited by: §E.3. Efron and Tibshirani (1994) B. Efron and R. J. Tibshirani An introduction to the bootstrap. Chapman and Hall/CRC. Cited by: §E.2. Egami and Imai (2019) N. Egami and K. Imai Causal interaction in factorial experiments: application to conjoint analysis. Journal of the American Statistical Association 114 (526), p. 529–540. External Links: Document, Link, https://doi.org/10.1080/01621459.2018.1476246 Cited by: §3. Freidlin et al. (2008) B. Freidlin, E. L. Korn, R. Gray, and A. Martin Multi-arm clinical trials of new agents: some design considerations. Clinical Cancer Research 14 (14), p. 4368–4371. External Links: ISSN 1078-0432, Document, Link, https://aacrjournals.org/clincancerres/article-pdf/14/14/4368/1924099/4368.pdf Cited by: §D.2, §D.2. Genevay et al. (2018) A. Genevay, G. Peyre, and M. Cuturi Learning generative models with sinkhorn divergences. In Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, A. Storkey and F. Perez-Cruz (Eds.), Proceedings of Machine Learning Research, Vol. 84, p. 1608–1617. External Links: Link Cited by: §D.1. Gong et al. (2022) F. Gong, Y. Nie, and H. Xu Gromov-wasserstein multi-modal alignment and clustering. In Proceedings of the 31st ACM International Conference on Information & Knowledge Management, CIKM ’22, New York, NY, USA, p. 603–613. External Links: ISBN 9781450392365, Link, Document Cited by: §2. Gradman et al. (2010) A. H. Gradman, J. N. Basile, B. L. Carter, and G. L. Bakris Combination therapy in hypertension. Journal of the American Society of Hypertension 4 (2), p. 90–98. External Links: ISSN 1933-1711, Document, Link Cited by: §1. Hill (2011) J. L. Hill Bayesian nonparametric modeling for causal inference. Journal of Computational and Graphical Statistics 20 (1), p. 217–240. External Links: Document, Link, https://doi.org/10.1198/jcgs.2010.08162 Cited by: §4. Imbens and Rubin (2015) G. W. Imbens and D. B. Rubin Causal inference in statistics, social, and biomedical sciences. Cambridge university press. Cited by: §3. Kingma and Ba (2017) D. P. Kingma and J. Ba Adam: a method for stochastic optimization. External Links: 1412.6980, Link Cited by: Appendix C. Kubota and Hoshino (2025) K. Kubota and T. Hoshino Impact of lottery promotion wins and losses: evidence from a promotion in a mobile payment service. Asia Pacific Journal of Marketing and Logistics, p. 1–19. External Links: ISSN 1355-5855, Document, Link, https://w.emerald.com/apjml/article-pdf/doi/10.1108/APJML-03-2025-0514/10974117/apjml-03-2025-0514en.pdf Cited by: §E.2. Kubota et al. (2025) K. Kubota, K. Ochiai, and T. Hoshino Causal effect of lottery promotions on post-win payments: evidence from a large field experiment. Innovative Marketing 21 (2), p. 225–236. External Links: Document Cited by: §E.3. Langen and Huber (2023) H. Langen and M. Huber How causal machine learning can leverage marketing strategies: assessing and improving the performance of a coupon campaign. Plos one 18 (1), p. e0278937. Cited by: §6. LeCun et al. (2015) Y. LeCun, Y. Bengio, and G. Hinton Deep learning. Nature 521 (7553), p. 436–444. External Links: Document, Link, ISSN 1476-4687 Cited by: Appendix C. Lesscher et al. (2021) L. Lesscher, L. Lobschat, and P. C. Verhoef Do offline and online go hand in hand? cross-channel and synergy effects of direct mailing and display advertising. International Journal of Research in Marketing 38 (3), p. 678–697. External Links: ISSN 0167-8116, Document, Link Cited by: §1. Li et al. (2019) J. Li, M. Jing, K. Lu, L. Zhu, and H. T. Shen Locality preserving joint transfer for domain adaptation. IEEE Transactions on Image Processing 28 (12), p. 6103–6115. External Links: Document Cited by: §3. Li et al. (2023) Q. Li, Z. Wang, S. Liu, G. Li, and G. Xu Causal optimal transport for treatment effect estimation. IEEE Transactions on Neural Networks and Learning Systems 34 (8), p. 4083–4095. External Links: Document Cited by: §2. Lian et al. (2021) X. Lian, K. Jain, J. Truszkowski, P. Poupart, and Y. Yu Unsupervised multilingual alignment using wasserstein barycenter. In Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence, IJCAI’20. External Links: ISBN 9780999241165 Cited by: §2. Liu et al. (2016) H. Liu, M. Shao, and Y. Fu Structure-preserved multi-source domain adaptation. In 2016 IEEE 16th International Conference on Data Mining (ICDM), Vol. , p. 1059–1064. External Links: Document Cited by: §2. Liu (2007) Y. Liu The long-term impact of loyalty programs on consumer purchase behavior and loyalty. Journal of Marketing 71 (4), p. 19–35. External Links: Document, Link, https://doi.org/10.1509/jmkg.71.4.019 Cited by: §E.3. Lu et al. (2025) M. Lu, L. Xing, and B. Chen Measuring generalized divergence for multiple distributions with application to deep clustering. Pattern Recognition 157, p. 110864. External Links: ISSN 0031-3203, Document, Link Cited by: §2. Murakami et al. (2025) Y. Murakami, T. Hattori, and K. Kubota Multiple treatments causal effects estimation with task embeddings and balanced representation learning. External Links: 2511.09814, Link Cited by: §A.1, 1st item, §E.1, §1, §2, §5, §5. Muralidharan et al. (2025) K. Muralidharan, M. Romero, and K. Wüthrich Factorial designs, model selection, and (incorrect) inference in randomized experiments. The Review of Economics and Statistics 107 (3), p. 589–604. External Links: ISSN 0034-6535, Document, Link, https://direct.mit.edu/rest/article-pdf/107/3/589/2075022/rest_a_01317.pdf Cited by: §D.2. Nair and Hinton (2010) V. Nair and G. E. Hinton Rectified linear units improve restricted boltzmann machines. In Proceedings of the 27th International Conference on International Conference on Machine Learning, ICML’10, Madison, WI, USA, p. 807–814. External Links: ISBN 9781605589077 Cited by: 2nd item. Parbhoo et al. (2021) S. Parbhoo, S. Bauer, and P. Schwab NCoRE: neural counterfactual representation learning for combinations of treatments. External Links: 2103.11175, Link Cited by: §2, §5. Parmar et al. (2008) M. K. B. Parmar, Friederike. M.-S. Barthel, M. Sydes, R. Langley, R. Kaplan, E. Eisenhauer, M. Brady, N. James, M. A. Bookman, A. Swart, W. Qian, and P. Royston Speeding up the evaluation of new agents in cancer. JNCI: Journal of the National Cancer Institute 100 (17), p. 1204–1214. External Links: ISSN 0027-8874, Document, Link, https://academic.oup.com/jnci/article-pdf/100/17/1204/7682916/djn267.pdf Cited by: §D.2. Qian et al. (2021) Z. Qian, A. Curth, and M. van der Schaar Estimating multi-cause treatment effects via single-cause perturbation. In Advances in Neural Information Processing Systems, M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. W. Vaughan (Eds.), Vol. 34, p. 23754–23767. External Links: Link Cited by: 2nd item, §2. Rissanen and Marttinen (2021) S. Rissanen and P. Marttinen A critical look at the consistency of causal estimation with deep latent variable models. In Proceedings of the 35th International Conference on Neural Information Processing Systems, NIPS ’21, Red Hook, NY, USA. External Links: ISBN 9781713845393 Cited by: §2. Rubin (2005) D. B. Rubin Causal inference using potential outcomes. Journal of the American Statistical Association 100 (469), p. 322–331. External Links: Document, Link, https://doi.org/10.1198/016214504000001880 Cited by: §3. Saini et al. (2019) S. K. Saini, S. Dhamnani, A. Srinivasan, A. A. Ibrahim, and P. Chavan Multiple treatment effect estimation using deep generative model with task embedding. In The World Wide Web Conference, W ’19, New York, NY, USA, p. 1601–1611. External Links: ISBN 9781450366748, Link, Document Cited by: 1st item, §2, §5. Schweisthal et al. (2023) J. Schweisthal, D. Frauen, V. Melnychuk, and S. Feuerriegel Reliable off-policy learning for dosage combinations. In Advances in Neural Information Processing Systems, A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine (Eds.), Vol. 36, p. 67900–67924. External Links: Link Cited by: 5th item. Sever et al. (2006) P. Sever, B. Dahlöf, N. Poulter, H. Wedel, G. Beevers, M. Caulfield, R. Collins, S. Kjeldsen, A. Kristinsson, G. McInnes, J. Mehlsen, M. Nieminem, E. O’Brien, and o. b. o. t. A. S. C. M. Östergren Potential synergy between lipid-lowering and blood-pressure-lowering in the anglo-scandinavian cardiac outcomes trial. European Heart Journal 27 (24), p. 2982–2988. External Links: ISSN 0195-668X, Document, Link, https://academic.oup.com/eurheartj/article-pdf/27/24/2982/9627980/ehl403.pdf Cited by: §1. Shalit et al. (2017) U. Shalit, F. D. Johansson, and D. Sontag Estimating individual treatment effect: generalization bounds and algorithms. In Proceedings of the 34th International Conference on Machine Learning, D. Precup and Y. W. Teh (Eds.), Proceedings of Machine Learning Research, Vol. 70, p. 3076–3085. External Links: Link Cited by: §A.3.3, 1st item, 1st item, Appendix C, §1, §2, §3. Sriperumbudur et al. (2010) B. K. Sriperumbudur, K. Fukumizu, A. Gretton, B. Schölkopf, and G. R. G. Lanckriet Non-parametric estimation of integral probability metrics. In 2010 IEEE International Symposium on Information Theory, Vol. , p. 1428–1432. External Links: Document Cited by: 3rd item. Thakral and Tô (2023) N. Thakral and L. T. Tô When are estimates independent of measurement units?. Boston University-Department of Economics. Cited by: §E.1. Titouan et al. (2019) V. Titouan, N. Courty, R. Tavenard, C. Laetitia, and R. Flamary Optimal transport for structured data with application on graphs. In Proceedings of the 36th International Conference on Machine Learning, K. Chaudhuri and R. Salakhutdinov (Eds.), Proceedings of Machine Learning Research, Vol. 97, p. 6275–6284. External Links: Link Cited by: §1, §4. Tsuboi et al. (2024) Y. Tsuboi, Y. Sakai, R. Shimizu, and M. Goto Multiple treatment effect estimation for business analytics using observational data. Cogent Engineering 11 (1), p. 2300557. External Links: Document, Link, https://doi.org/10.1080/23311916.2023.2300557 Cited by: §D.2. Vayer et al. (2020) T. Vayer, L. Chapel, R. Flamary, R. Tavenard, and N. Courty Fused gromov-wasserstein distance for structured objects. Algorithms 13 (9). External Links: Link, ISSN 1999-4893, Document Cited by: §1, §4. Villani et al. (2008) C. Villani et al. Optimal transport: old and new. Vol. 338, Springer. Cited by: §A.3.2. Wang et al. (2025) H. Wang, Z. Chen, Z. Liu, X. Chen, H. Li, and Z. Lin Proximity matters: local proximity enhanced balancing for treatment effect estimation. In Proceedings of the 31st ACM SIGKDD Conference on Knowledge Discovery and Data Mining V.2, KDD ’25, New York, NY, USA, p. 2927–2937. External Links: ISBN 9798400714542, Link, Document Cited by: §A.3.3, §1, §2, §3, §4. Xu and Mannor (2012) H. Xu and S. Mannor Robustness and generalization. Machine Learning 86 (3), p. 391–423. External Links: Document, Link Cited by: §3. Xu et al. (2020) J. Xu, Z. Li, B. Du, M. Zhang, and J. Liu Reluplex made more practical: Leaky ReLU. In 2020 IEEE Symposium on Computers and Communications (ISCC), Vol. , p. 1–7. External Links: Document Cited by: Appendix C. Yamada et al. (2022) M. Yamada, Y. Takezawa, R. Sato, H. Bao, Z. Kozareva, and S. Ravi Approximating 1-wasserstein distance with trees. External Links: 2206.12116, Link Cited by: §D.1. Yao et al. (2018) L. Yao, S. Li, Y. Li, M. Huai, J. Gao, and A. Zhang Representation learning for treatment effect estimation from observational data. In Advances in Neural Information Processing Systems, S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (Eds.), Vol. 31, p. . External Links: Link Cited by: §1, §1, §2. Yao et al. (2019) L. Yao, S. Li, Y. Li, M. Huai, J. Gao, and A. Zhang ACE: adaptively similarity-preserved representation learning for individual treatment effect estimation. In 2019 IEEE International Conference on Data Mining (ICDM), Vol. , p. 1432–1437. External Links: Document Cited by: §2. Zou et al. (2020) H. Zou, P. Cui, B. Li, Z. Shen, J. Ma, H. Yang, and Y. He Counterfactual prediction for bundle treatment. In Advances in Neural Information Processing Systems, H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (Eds.), Vol. 33, p. 19705–19715. External Links: Link Cited by: 4th item, §2. Barycentric Fused Gromov-Wasserstein Balancing for Causal Inference under Multiple Treatments (Supplementary Material) Appendix A Proofs of Theoretical Results A.1 Proof of Proposition 1 Lemma 1. Under Assumptions 1-3, the conditional expected potential outcome μ(,)μ( x, t) is identifiable from the observed data. Proof. The identifiability of the conditional expected potential outcome is a standard result in the potential outcome framework. The equality holds directly from Assumptions 1 and 2, while Assumption 3 ensures the conditional expected potential outcome is well-defined. For a detailed derivation, see, for example, 36. ∎ Using Lemma 1, we prove that CASE and CAIE are identifiable. Proof. Recall that the CASE for treatment k, denoted as τCASE(k,) _CASE(k, x), is defined as τCASE(k,)≔μ(,+k)−μ(,), _CASE(k, x) μ( x, t_+k)-μ( x, 0), where +k t_+k is the one-hot vector for treatment k. Similarly, recall that the CAIE for a set of treatments S, denoted as τCAIE(S,) _CAIE(S, x), is defined as τCAIE(S,)≔∑Q⊆S(−1)|S|−|Q|μ(,(+Q)). _CAIE(S, x) _Q S(-1)^|S|-|Q|\,μ ( x, t_(+Q) ). where (+Q) t_(+Q) corresponds to the treatment vector where treatments in the subset Q are active. According to Lemma 1, for any treatment vector ∈0,1K t∈\0,1\^K, the term μ(,)μ( x, t) is identifiable from the observed data as [Y∣=,=]E[Y X= x, T= t]. Both τCASE(k,) _CASE(k, x) and τCAIE(S,) _CAIE(S, x) are constructed as linear combinations of μ(,)μ( x, t) for specific values of t. Since every term μ(,)μ( x, t) in these equations is identifiable, any linear combination of them is also identifiable from the observed data. Therefore, CASE and CAIE are identifiable. ∎ A.2 Proof of Proposition 2 Proof. Let L=2KL=2^K denote the total number of treatment patterns derived from K binary treatments. Standard pairwise balancing evaluates an OT-based discrepancy for every unordered pair of distinct treatment-pattern representation distributions. Therefore, the total number of evaluations per training step is Npair=(L2)=L(L−1)2=O(L2).N_pair= L2= L(L-1)2=O(L^2). For BFG-WB, the total number of evaluations consists of two components. First, for Wasserstein barycenter estimation, each barycenter-update iteration solves one OT problem between the current barycenter and each of the L input treatment-pattern distributions (13). Therefore, under a fixed iteration budget of at most LbL_b barycenter-update iterations, the number of OT evaluations required for barycenter estimation satisfies Nbary≤LbLN_bary≤ L_bL where LbL_b is independent of L. Second, after the barycenter has been estimated, BFG-WB evaluates one FGW discrepancy between the estimated barycenter and each treatment-pattern representation distribution. Hence, this step requires NFGW=LN_FGW=L additional evaluations. Consequently, the total number of OT-based discrepancy evaluations required by BFG-WB per training step satisfies NBFG-WB=Nbary+NFGW≤LbL+L=(Lb+1)L.N_BFG -WB=N_bary+N_FGW≤ L_bL+L=(L_b+1)L. Because LbL_b is fixed independently of L, it follows that NBFG-WB=O(L).N_BFG -WB=O(L). Thus, standard pairwise balancing requires a quadratic number of evaluations in L, whereas BFG-WB requires at most a linear number of evaluations. ∎ A.3 Proof of Theorems 1 and 2 This appendix details the proofs of Theorems 1 and 2, organized into three parts. First, we summarize the notation, definitions, and key assumptions used throughout the proofs. Next, we present supporting lemmas that establish intermediate inequalities connecting counterfactual losses with distributional discrepancies. Finally, building on these results, we provide the main proofs of Theorems 1 and 2, which demonstrate how barycenter-based balancing yields the stated upper bounds. A.3.1 Setup and Definitions Notation and setup. Let ϕ:→ℛ⊂ℝdrφ:X ^d_r be the representation map, where ℛR is the representation space. For each treatment pattern ∈ t , we denote by p(∣)p( x t) the covariate distribution conditional on = T= t, and define R_ t as the distribution on ℛR of =ϕ() r=φ( x) when ∼p(∣) x p( x t). Throughout the theoretical analysis, R_ t denotes the population distribution on ℛR; in practice, we use its empirical distribution computed from the observed samples in the treatment pattern t. Given (,)( x, t), we form the predictor input by concatenating the representation and the task embedding tw()t_w( t), and define the point predictor of the potential outcome as μ^(,):=h([ϕ(),tw()]) μ( x, t)\;:=\;h\! ([φ( x),\,t_w( t)] ) where h is the outcome prediction network and [⋅,⋅][·,·] denotes concatenation. Definition 1. For covariates x and treatment pattern t, define the conditional expected squared loss l(,)=∫(Y()−μ^(,))2p(Y()∣)Y().l( x, t)= (Y( t)- μ( x, t) )^2\,p (Y( t) x )\,dY( t). The expected factual loss under t and the expected counterfactual loss for predicting outcomes under t using covariates drawn from other treatment patterns are: ϵF()=∫l(,)p(∣)p(), _F^( t)= _Xl( x, t)\,p( x t)\,p( t)\,d x, (A.1) ϵCF()=∑′∈,′≠∫l(,)p(∣′)p(′). _CF^( t)= _ t ,\, t ≠ t _Xl( x, t)\,p( x t )\,p( t )\,d x. (A.2) Definition 2. Let P and Q be probability measures on ℛ⊂ℝdrR ^d_r with finite p-th moments. The p-Wasserstein distance is defined as Wp(P,Q)≔(infπ∈Π(P,Q)∫ℛ×ℛ‖r−r′‖2pπ(r,r′))1/p,W_p(P,Q) ( _π∈ (P,Q) _R×R\|r-r \|_2^p\,dπ(r,r ) )^1/p, where Π(P,Q) (P,Q) denotes the set of all couplings on ℛ×ℛR×R with marginals P and Q. Definition 3. Let ℱF be a class of real-valued functions on ℛR. For two distributions P and Q on ℛR, the integral probability metric (IPM) induced by ℱF is defined as IPMℱ(P,Q)=supf∈ℱ|∫ℛf(r)(dP(r)−dQ(r))|.IPM_F(P,Q)= _f | _Rf(r)\,(dP(r)-dQ(r)) |. Assumption 1 (Stable Unit Treatment Value Assumption). (1) no interference, meaning that the outcome of one unit is unaffected by the treatment assignments of other units; and (2) consistency of treatment, meaning that the potential outcomes correspond to well-defined and unique treatments (i.e., yi=Yi(i)y_i=Y_i( t_i)). Assumption 2 (Ignorability). For any treatment pattern, the potential outcome is independent of the assigned treatment T given the observed covariates X. Formally, for all t, Y()⟂|.Y( t) T X. Assumption 3 (Overlap). Every unit has a non-zero probability of receiving any treatment pattern given its observed covariates. Formally, for all t and x, 0<P(=∣=)<1.0<P( T= t X= x)<1. Assumption 4 (Invertible representation map). The representation map ϕ:→ℛ⊂ℝdrφ:X ^d_r is one-to-one on X. Without loss of generality, we assume that ℛR is the image of X under ϕφ. Hence, there exists an inverse map Ψ:ℛ→ :R such that Ψ(ϕ())=(∀∈). (φ( x))= x (∀ x ). Assumption 5 (Lipschitz loss in representation space). There exists a constant Bϕ>0B_φ>0 such that, for any ∈ t , the function g():=1Bϕl(Ψ(),)g_ t( r):= 1B_φ\,l ( ( r), t ) belongs to the class of 11-Lipschitz functions on ℛR. A.3.2 Supporting Lemmas We present supporting lemmas and clarify their connections to the proofs of the main theorems. The lemmas can be grouped into three categories. First, Lemmas 2 and 3 relate counterfactual outcome losses to the corresponding factual outcome losses, plus Wasserstein distances between treatment-specific representation distributions. Lemma 2 serves as a technical preparation for Lemma 3: it shows that the increase in the expected outcome prediction loss for ′ t under covariates (and thus representations) from t is controlled by the Wasserstein distance between the representation distributions of t and ′ t . Lemma 2. Under Assumptions 1-5 , for any ,′∈ t, t , we have ∫l(,′)p(∣)≤∫l(,′)p(∣′)+BϕW1(,′). _Xl( x, t )p( x t)d x≤ _Xl( x, t )p( x t )d x+B_φW_1( R_ t, R_ t ). Proof. By the definition of R_ t and the Assumption 4, we obtain ∫l(,′)p(∣) _Xl( x, t )\,p( x t)\,d x =∫ℛl(Ψ(),′)dR(), = _Rl( ( r), t )\,dR_ t( r), ∫l(,′)p(∣′) _Xl( x, t )\,p( x t )\,d x =∫ℛl(Ψ(),′)dR′(). = _Rl( ( r), t )\,dR_ t ( r). Therefore, we can rewrite the difference of expectations as an integral over ℛR as follows. ∫l(,′)p(∣)−∫l(,′)p(∣′) _Xl( x, t )\,p( x t)\,d x- _Xl( x, t )\,p( x t )\,d x =∫ℛl(Ψ(),′)(dR()−dR′()) = _Rl ( ( r), t )\, (dR_ t( r)-dR_ t ( r) ) =Bϕ∫ℛg′()(dR()−dR′()) =B_φ _Rg_ t ( r)\, (dR_ t( r)-dR_ t ( r) ) ≤Bϕ|∫ℛg′()(dR()−dR′())| ≤ B_φ | _Rg_ t ( r)\, (dR_ t( r)-dR_ t ( r) ) | ≤BϕIPMℱ(,′)(by Assumption 5), ≤ B_φIPM_F( R_ t, R_ t ) (by Assumption~ assumption:lipschitz_loss), where ℱ=Lip1(ℛ)F=Lip_1(R) denotes the set of all 11-Lipschitz functions on ℛR. By the Kantorovich-Rubinstein duality (53), the IPM over 11-Lipschitz functions equals the 1-Wasserstein distance, IPMℱ(,′)=W1(,′).IPM_F( R_ t, R_ t )=W_1( R_ t, R_ t ). Accordingly, we have ∫l(,′)p(∣)≤∫l(,′)p(∣′)+BϕW1(,′). _Xl( x, t )\,p( x t)\,d x≤ _Xl( x, t )\,p( x t )\,d x+B_φW_1( R_ t, R_ t ). ∎ Building on this result, Lemma 3 upper-bounds the counterfactual outcome loss for treatment pattern t by its factual outcome loss plus Wasserstein discrepancies to other treatment patterns. Lemma 3. Under the Assumptions 1-5, for any treatment pattern ∈ t , the expected counterfactual loss satisfies ϵCF()≤(1p()−1)ϵF()+Bϕ∑′≠W1(,′). _CF^( t)≤ ( 1p( t)-1 ) _F^( t)\;+\;B_φ _ subarrayc t ≠ t subarray\,W_1\! ( R_ t, R_ t ). Proof. Applying Lemma 2 to each ′ t yields ∫l(,)p(∣′)≤∫l(,)p(∣)+BϕW1(,′). _Xl( x, t)\,p( x t )\,d x≤ _Xl( x, t)\,p( x t)\,d x+B_φ\,W_1\! ( R_ t, R_ t ). (A.3) Multiplying both sides of (A.3) by p(′)p( t ) and summing over ′≠ t ≠ t, we obtain ϵCF()≤∑′≠p(′)∫l(,)p(∣)d+Bϕ∑′≠p(′)W1(,′)=(∑′≠p(′))∫l(,)p(∣)d+Bϕ∑′≠p(′)W1(,′)=(1−p())∫l(,)p(∣)d+Bϕ∑′≠p(′)W1(,′)=(1p()−1)ϵF()+Bϕ∑′≠p(′)W1(,′)(by Definition in Equation (A.1))≤(1p()−1)ϵF()+Bϕ∑′≠W1(,′). split _CF^( t)&≤ _ t ≠ tp( t ) _Xl( x, t)\,p( x t)\,d x+B_φ _ t ≠ tp( t )W_1\! ( R_ t, R_ t )\\ &= ( _ t ≠ tp( t ) ) _Xl( x, t)\,p( x t)\,d x+B_φ _ t ≠ tp( t )W_1\! ( R_ t, R_ t )\\ &=(1-p( t)) _Xl( x, t)\,p( x t)\,d x+B_φ _ t ≠ tp( t )W_1\! ( R_ t, R_ t )\\ &= ( 1p( t)-1 ) _F^( t)+B_φ _ t ≠ tp( t )W_1\! ( R_ t, R_ t ) (by Definition in Equation~ app_eq:factual_loss)\\ &≤ ( 1p( t)-1 ) _F^( t)+B_φ _ t ≠ tW_1\! ( R_ t, R_ t ). split where the last inequality holds because p(′)≤1p( t )≤ 1 and W1(⋅,⋅)≥0W_1(·,·)≥ 0. ∎ Second, Lemma 4 provides an upper bound on aggregate discrepancies across multiple treatment patterns in terms of distances to a common Wasserstein barycenter. This bound demonstrates that the sum of all pairwise discrepancies can be controlled by reducing the discrepancies between each representation distribution and the Wasserstein barycenter. This is exactly what our barycenter-based regularization term ℒϕL_φ enforces. Lemma 4. Let T denote the set of all treatment patterns with ||=2K|T|=2^K. For any p≥1p≥ 1, the following inequality holds: ∑,′∈(2)Wp(,′)≤(2K−1)∑∈Wp(b∗,), _\ t, t \∈ T2W_p\! ( R_ t, R_ t )≤ (2^K-1 ) _ t W_p\! ( R_b , R_ t ), where b∗ R_b is a Wasserstein barycenter of ∈\ R_ t\_ t . Proof. By the triangle inequality for the p-Wasserstein distance defined in Definition 2, for any ,′∈ t, t , Wp(,′)≤Wp(,b∗)+Wp(b∗,′).W_p\! ( R_ t, R_ t )≤ W_p\! ( R_ t, R_b )+W_p\! ( R_b , R_ t ). Summing the above inequality over all ordered pairs (,′)( t, t ) with ≠′ t≠ t yields ∑≠′Wp(,′) _ t≠ t W_p\! ( R_ t, R_ t ) ≤∑≠′Wp(,b∗)+Wp(b∗,′) ≤ _ t≠ t \W_p\! ( R_ t, R_b )+W_p\! ( R_b , R_ t ) \ =(||−1)∑∈Wp(,b∗)+(||−1)∑∈Wp(b∗,), =(|T|-1) _ t W_p\! ( R_ t, R_b )+(|T|-1) _ t W_p\! ( R_b , R_ t ), since each ∈ t appears exactly ||−1|T|-1 times as the first element of an ordered pair. Using the symmetry of WpW_p, we have ∑∈Wp(,b∗)=∑∈Wp(b∗,) _ t W_p( R_ t, R_b )= _ t W_p( R_b , R_ t), and thus ∑≠′Wp(,′)≤2(||−1)∑∈Wp(b∗,). _ t≠ t W_p\! ( R_ t, R_ t )≤ 2(|T|-1) _ t W_p\! ( R_b , R_ t ). Finally, noting that ∑≠′Wp(,′)=2∑,′∈(2)Wp(,′) _ t≠ t W_p( R_ t, R_ t )=2 _\ t, t \∈ T2W_p( R_ t, R_ t ), we obtain ∑,′∈(2)Wp(,′)≤(||−1)∑∈Wp(b∗,)=(2K−1)∑∈Wp(b∗,), _\ t, t \∈ T2W_p\! ( R_ t, R_ t )≤(|T|-1) _ t W_p\! ( R_b , R_ t )= (2^K-1 ) _ t W_p\! ( R_b , R_ t ), which completes the proof. ∎ Although this lemma holds for any anchor distribution, using the Wasserstein barycenter can yield a smaller right-hand side, leading to a tighter anchor-based bound. Third, Lemma 5 provides an upper bound that relates the above Wasserstein discrepancies to the FGW discrepancy adopted in our objective, thereby allowing us to control the Wasserstein terms via FGW. Whereas replacing the Wasserstein terms with an FGW-based upper bound may loosen the theoretical bound, introducing FGW is crucial in practice because it encourages the preservation of local proximity structures in the representation space. Lemma 5. Let 0<η≤10<η≤ 1 and let F(⋅,⋅)F(·,·) be the FGW discrepancy defined in Equation (6). Then, for any treatment pattern ∈ t , W1(,b∗)≤1ηF(,b∗).W_1\! ( R_ t, R_b )≤ 1η\,F\! ( R_ t, R_b ). Proof. For any coupling π∈Π(,b∗)π∈ ( R_ t, R_b^*), write the FGW objective as J(π)=ηA(π)+(1−η)B(π)J(π)=η A(π)+(1-η)B(π), where A(π)A(π) is the feature-based transport cost and B(π)≥0B(π)≥ 0 is the structural term. By Definition 2, A(π)≥infπ~∈Π(,b∗)A(π~)=W1(R,Rb∗).A(π)≥ _ π∈ ( R_ t, R_b^*)A( π)=W_1(R_ t,R_b^*). Therefore, J(π)≥ηW1(,b∗)J(π)≥η W_1( R_ t, R_b^*). Taking the infimum over π∈Π(,b∗)π∈ ( R_ t, R_b^*) yields F(,b∗)≥ηW1(,b∗)F( R_ t, R_b^*)≥η W_1( R_ t, R_b^*) ,which proves the result. ∎ Combining Lemmas 2 - 5, we can upper-bound the counterfactual outcome losses terms in Theorems 1 and 2 by the corresponding factual outcome losses plus the barycenter-based regularizer (via Lemma 4) and the FGW discrepancy term in our objective (via Lemma 5). We then obtain the stated bounds by substituting these inequalities into the decompositions used in the proofs of Theorems 1 and 2. A.3.3 Main Proofs: Theorems 1 and 2 Using the above lemmas, we prove Theorems 1 and 2. Theorem 1 (Upper bound for CASE). Suppose that Assumptions 1–3 and the auxiliary conditions in Appendix A.3.1 hold (in particular, the expected squared loss l(,)l( x, t) is BϕB_φ-Lipschitz in the representation space), and that the weights in Equation (4) are uniform, i.e., w=2−Kw_ t=2^-K. Then, for any k∈1,…,Kk∈\1,…,K\, ϵCASE(k)≤2(1p(+k)ϵF(+k)+1p()ϵF()+22KηBϕℒϕ). _CASE(k)≤ 2 ( 1p( t_+k)\, _F^( t_+k)+ 1p( 0)\, _F^( 0)+ 2^2Kη\,B_φL_φ ). Proof. By definition in Equation (7), we rewrite ϵCASE(k) _CASE(k) =∫[μ^(,+k)−μ(,+k)−μ^(,)−μ(,)]2p(). = _X [ \ μ( x, t_+k)-μ( x, t_+k) \- \ μ( x, 0)-μ( x, 0) \ ]^2\,p( x)\,d x. Therefore, we have ϵCASE(k) _CASE(k) ≤2∫[μ^(,+k)−μ(,+k)2+μ^(,)−μ(,)2]p() ≤ 2 _X [ \ μ( x, t_+k)-μ( x, t_+k) \^2+ \ μ( x, 0)-μ( x, 0) \^2 ]\,p( x)\,d x =2∑′∈∫μ^(,+k)−μ(,+k)2p(∣′)p(′)+2∑′∈∫μ^(,)−μ(,)2p(∣′)p(′) =2 _ t _X \ μ( x, t_+k)-μ( x, t_+k) \^2\,p( x t )\,p( t )d x+2 _ t _X \ μ( x, 0)-μ( x, 0) \^2\,p( x t )\,p( t )d x ≤2(ϵF(+k)+ϵCF(+k))+2(ϵF()+ϵCF())(by the Definition of Equations (A.1) and (A.2)) ≤ 2 ( _F^( t_+k)+ _CF^( t_+k) )+2 ( _F^( 0)+ _CF^( 0) ) (by the Definition of Equations~ app_eq:factual_loss and app_eq:counterfactual_loss) ≤2(1p(+k)ϵF(+k)+1p()ϵF())+2Bϕ(∑′≠+kW1(+k,′)+∑′≠W1(,′))(by Lemma 3), ≤ 2 ( 1p( t_+k) _F^( t_+k)+ 1p( 0) _F^( 0) )+2B_φ ( _ t ≠ t_+kW_1( R_ t_+k, R_ t )+ _ t ≠ 0W_1( R_ 0, R_ t ) ) (by Lemma~ lem:cf_bound), where we use (μ^(,)−μ(,))2≤l(,)( μ( x, t)-μ( x, t))^2≤ l( x, t), which follows from the law of total variance. To connect the remaining Wasserstein terms to the BFG-WB regularization, we upper bound the pairwise sums by barycenter-based discrepancies. Let b∗ R_b denote the Wasserstein barycenter over ∈\ R_ t\_ t . For any fixed ∈ t , by the triangle inequality of W1W_1, ∑′≠W1(,′)≤∑′≠W1(,b∗)+W1(′,b∗)=(2K−1)W1(,b∗)+∑′≠W1(′,b∗)=(2K−2)W1(,b∗)+∑′∈W1(′,b∗). split _ t ≠ tW_1( R_ t, R_ t )&≤ _ t ≠ t \W_1( R_ t, R_b )+W_1( R_ t , R_b ) \\\ &=(2^K-1)\,W_1( R_ t, R_b )+ _ t ≠ tW_1( R_ t , R_b )\\ &=(2^K-2)\,W_1( R_ t, R_b )+ _ t W_1( R_ t , R_b ). split (A.4) Applying (A.4) to =+k t= t_+k and = t= 0 and summing them yields ∑′≠+kW1(+k,′)+∑′≠W1(,′) _ t ≠ t_+kW_1( R_ t_+k, R_ t )+ _ t ≠ 0W_1( R_ 0, R_ t ) ≤(2K−2)W1(+k,b∗)+W1(,b∗)+2∑′∈W1(′,b∗) ≤(2^K-2) \W_1( R_ t_+k, R_b )+W_1( R_ 0, R_b ) \+2 _ t W_1( R_ t , R_b ) ≤2K∑∈W1(,b∗) ≤ 2^K _ t W_1( R_ t, R_b ) ≤2Kη∑∈F(,b∗) ≤ 2^Kη _ t F( R_ t, R_b ) ≤22Kηℒϕ. ≤ 2^2KηL_φ. ∎ Theorem 2 (Upper Bound for CAIE). Under the same assumptions as Theorem 1, consider an interaction set S⊆1,…,KS \1,…,K\ with |S|≥2|S|≥ 2. Let (a)∈(a_ t)_ t be a constant vector reflecting the combinatorial structure of S (a∈−1,0,1a_ t∈\-1,0,1\). The estimation error ϵCAIE(S) _CAIE(S) is bounded by: ϵCAIE(S)≤(∑∈a2)∑∈1p()ϵF()+2K+1ηBϕ(2K−1)ℒϕ. _CAIE(S)\;≤\; ( _ t a_ t^2 ) \ _ t 1p( t)\, _F^( t)\;+\; 2^K+1η\,B_φ\,(2^K-1)\,L_φ \. Proof. ϵCAIE(S) _CAIE(S) =∫∑a(μ^(,)−μ(,))2p()d(by Definition in Equation (8)) = _X \ _ ta_ t ( μ( x, t)-μ( x, t) ) \^2p( x)\,d x (by Definition in Equation~ eq:def_caie_error) ≤(∑a2)∑∫(μ^(,)−μ(,))2p()d(by Cauchy–Schwarz) ≤ ( _ ta_ t^2 ) _ t _X ( μ( x, t)-μ( x, t) )^2p( x)\,d x (by Cauchy--Schwarz) =(∑a2)∑′∫(μ^(,)−μ(,))2p(∣′)p(′) = ( _ ta_ t^2 ) _ t _ t _X ( μ( x, t)-μ( x, t) )^2\,p( x t )\,p( t )\,d x =(∑a2)∑∫(μ^(,)−μ(,))2p(∣)p()+∑′≠∫(μ^(,)−μ(,))2p(∣′)p(′) = ( _ ta_ t^2 ) _ t \ _X ( μ( x, t)-μ( x, t) )^2\,p( x t)\,p( t)\,d x+ _ t ≠ t _X ( μ( x, t)-μ( x, t) )^2\,p( x t )\,p( t )\,d x \ ≤(∑a2)∑∫l(,)p(∣)p()+∑′≠∫l(,)p(∣′)p(′) ≤ ( _ ta_ t^2 ) _ t \ _Xl( x, t)\,p( x t)\,p( t)\,d x+ _ t ≠ t _Xl( x, t)\,p( x t )\,p( t )\,d x \ =(∑a2)∑(ϵF()+ϵCF())(by Definition in Equations (A.1) and (A.2)) = ( _ ta_ t^2 ) _ t ( _F^( t)+ _CF^( t) ) (by Definition in Equations~ app_eq:factual_loss and app_eq:counterfactual_loss) ≤(∑a2)∑ϵF()+(1p()−1)ϵF()+Bϕ∑′≠W1(,′) ≤ ( _ ta_ t^2 ) _ t \ _F^( t)+ ( 1p( t)-1 ) _F^( t)+B_φ\, _ t ≠ tW_1( R_ t, R_ t ) \ Note that ∑′≠W1(,′) _ t _ t ≠ tW_1( R_ t, R_ t ) counts each unordered pair ,′\ t, t \ twice. Hence, ∑′≠W1(,′)=2∑,′∈(2)W1(,′). _ t _ t ≠ tW_1( R_ t, R_ t )=2 _\ t, t \∈ T2W_1( R_ t, R_ t ). By substituting this identity into the previous inequality, we have ϵCAIE(S) _CAIE(S) =(∑a2)∑1p()ϵF()+2Bϕ∑,′∈(2)W1(,′) = ( _ ta_ t^2 ) \ _ t 1p( t)\, _F^( t)+2B_φ _\ t, t \∈ T2W_1( R_ t, R_ t ) \ ≤(∑a2)∑1p()ϵF()+2Bϕ(2K−1)∑W1(,b∗)(by Lemma 4) ≤ ( _ ta_ t^2 ) \ _ t 1p( t)\, _F^( t)+2B_φ(2^K-1) _ tW_1( R_ t, R_b ) \ (by Lemma~ lem:pairwise_to_barycenter_bound) ≤(∑a2)∑1p()ϵF()+2ηBϕ(2K−1)∑F(,b∗)(by Lemma 5) ≤ ( _ ta_ t^2 ) \ _ t 1p( t)\, _F^( t)+ 2ηB_φ(2^K-1) _ tF( R_ t, R_b ) \ (by Lemma~ lem:w1_upper_by_fgw) =(∑a2)∑1p()ϵF()+2K+1ηBϕ(2K−1)ℒϕ. = ( _ ta_ t^2 ) \ _ t 1p( t)\, _F^( t)+ 2^K+1η\,B_φ\,(2^K-1)\,L_φ \. ∎ Although the two theorems could be subsumed under a general linear-contrast bound, presenting them separately clarifies how BFG-WB controls the estimation errors of both single and interaction treatment effects. Together, these bounds provide a principled decomposition that motivates the training objective in Equation (10) by linking factual prediction errors and the BFG-WB discrepancy to the estimation errors of both estimands. It is important to note the connection between these theoretical bounds and our empirical design. Assumptions 4 and 5 are idealized conditions; they follow the standard setting adopted in representation-based causal inference (47; 54). Assumption 5 concerns the smoothness of the loss over the representation space, and is more plausible when nearby representations correspond to units with similar outcomes. By penalizing distortions of within-distribution neighborhoods, the proposed BFG-WB objective is designed to encourage representations of this kind, which we examine empirically in Appendix D.4 Appendix B Data Generating Process We describe the data-generating process (DGP) for the simulation datasets used in Section 5 to evaluate the ability of the models to estimate CASE and CAIE accurately and stably under selection bias. We generate two datasets under covariate-dependent treatment assignment and vary only whether the outcome function contains interaction effects. Simulation 1 includes interaction effects among treatments in the outcome generation, which serves as the primary benchmark for CASE and CAIE estimation. Simulation 2 excludes interaction effects to assess the model’s robustness and its ability to avoid inducing false interaction effects. In two scenarios, the simulation data follow the same functional form for covariate generation, treatment assignment, and outcome generation. The number of treatments K is fixed at 3, and the sample size N is set to 50,000. We introduce the indicator variables H to reflect realistic conditions where treatment assignment depends on covariates (3; 14). xi,n(j)∼N(cn(j),12),xi,u(j)∼U(−1,1),i=(i,n,i,u), x_i,n^(j) N(c_n^(j),1^2), x_i,u^(j) U(-1,1), x_i=( x_i,n, x_i,u), Hi(v)=(xi(v)+xi(v+1)>1),v=1,2,3, H_i^(v)=I (\,x_i^(v)+x_i^(v+1)>1\, ), v=1,2,3, i=(Hi(1),Hi(2),Hi(3)), H_i=(H_i^(1),H_i^(2),H_i^(3)), ti(k)∼Bern(σ(tk⊤i+H,k⊤i−1)),k=1,2,3, t_i^(k) (σ ( w_t_k x_i+ w_H,k H_i-1 ) ), k=1,2,3, =(ti(1),ti(2),ti(3)), t_i=(t_i^(1),t_i^(2),t_i^(3)), Yi∼(f(i,i,l), 12),l∈0,1, Y_i \! (f\! ( x_i, t_i,l ),\,1^2 ), l∈\0,1\, where j∈1,…,15j∈\1,…,15\ and cn(j)c_n^(j) are drawn from the uniform distribution U(−1,1)U(-1,1). (⋅)I(·) is the indicator function, σ(⋅)σ(·) is the sigmoid function defined as σ(x)=1/(1+exp(−x))σ(x)=1/(1+ (-x)). tk w_t_k and H,k w_H,k are the weight vectors whose elements are independently drawn from the uniform distribution U(−0.5,0.5)U(-0.5,0.5). The outcome-generating function f is defined as follows: f(i,i,l) f( x_i, t_i,l) =x⊤i+(xi(1)+1)ti(1)+1.2(xi(2)+1)ti(2) = w_x x_i+(x_i^(1)+1)t_i^(1)+1.2(x_i^(2)+1)t_i^(2) +0.8(xi(3)+1)ti(3)+l(xi(4)+0.5)ti(1)ti(2) +0.8(x_i^(3)+1)t_i^(3)+l \(x_i^(4)+0.5)t_i^(1)t_i^(2) −0.5(xi(5)−0.1)ti(1)ti(3)+0.3(xi(6)+0.1)ti(2)ti(3) -0.5(x_i^(5)-0.1)t_i^(1)t_i^(3)+0.3(x_i^(6)+0.1)t_i^(2)t_i^(3) +0.7(xi(7)+1)ti(1)ti(2)ti(3)+2, +0.7(x_i^(7)+1)t_i^(1)t_i^(2)t_i^(3) \+2, where vector x w_x is the weight vector whose elements are independently drawn from the uniform distribution U(−1,1)U(-1,1), and the parameter l∈0,1l∈\0,1\ determines whether interaction effect terms are included (l=1l=1) or excluded (l=0l=0). We define the case where l=1l=1 as the first scenario and the case where l=0l=0 as the second scenario. Appendix C Implementation Details This section provides the implementation details necessary to reproduce the simulation results reported in Section 5. We first specify the network architectures of the proposed CIHSI-Net, including representation learning, task embedding, and outcome prediction subnetworks. Subsequently, we outline the optimization protocols, hyperparameter settings, and the implementation configurations for the baseline methods. CIHSI-Net Architecture and Training. All neural networks within CIHSI-Net are constructed using fully connected (FC) layers (28) with Leaky ReLU activation functions (56). The specific architectures are as follows. • Representation Learning Network ϕφ: Consists of two hidden layers with 100 units each. The output dimension is set to 64. Crucially, we apply batch normalization with a fixed scale to the output of the representation learning network because the BFG-WB regularization is sensitive to feature scaling, similar to the balancing regularization in 47. • Task Embedding Network twt_w: Consists of two hidden layers with 100 units each, outputting a five-dimensional embedding vector. • Outcome Prediction Network h: Consists of two hidden layers with 100 units each. It takes the concatenated output of ϕφ and twt_w as input. During training, the Wasserstein barycenter is estimated for each mini-batch as a free-support W22W_2^2 barycenter with 16 uniformly weighted support points, using the treatment-pattern-specific empirical representation distributions. In all experiments, the weights λ _ t in Equation (5) are set uniformly across the treatment patterns included in the barycenter computation. The BFG-WB weights are also set uniformly as w=2−Kw_ t=2^-K for all treatment patterns. We compute the barycenter through iterative updates with a maximum of 100 iterations. The optimization procedure of CIHSI-Net is encapsulated in Algorithm 1. 1 Input: observed data =(i,i,yi)i=1ND=\( x_i, t_i,y_i)\_i=1^N; representation mapping ϕφ; task embedding mapping twt_w; outcome mapping h. 2 Parameters: α: strength of BFG-WB regularization; η: FGW trade-off parameter; β: strength of L2L_2 regularization; λ∈0,1K\ _ t\_ t∈\0,1\^K: Wasserstein barycenter weights; w∈0,1K\w_ t\_ t∈\0,1\^K: BFG-WB weights. 3 Output: ℒL: the learning objective of CIHSI-Net. 4 5 i←ϕ(i) R_i←φ( x_i); 6 7 ←i∣i= R_ t←\ R_i t_i= t\ for each ∈0,1K t∈\0,1\^K; 8 9 b∗←argmin∑∈bλW22(,b) R_b^*← _ R_b _ t _ tW^2_2( R_ t, R_b); 10 11 ℒϕ←∑∈wF(,b∗)L_φ← _ t w_ tF( R_ t, R_b^*); 12 13 p^()←1N∑i=1N[i=] p( t)← 1N _i=1^NI[ t_i= t] for each ∈0,1K t∈\0,1\^K; 14 15 Calculate ℒyL_y following (11) with p^() p( t) and outcome mapping h; 16 17 ℒ←ℒy+αℒϕ+β∥22L _y+ _φ+β w _2^2; 18 19 Algorithm 1 The computation workflow of CIHSI-Net with BFG-WB Hyperparameter Candidate values Selected value BFG-WB coefficient α 0.1,0.5,1.0,1.5,2.0,5.0\0.1,0.5,1.0,1.5,2.0,5.0\ 1.01.0 FGW trade-off η 10−5,0.1,0.2,…,0.9,1.0\10^-5,0.1,0.2,…,0.9,1.0\ 0.60.6 Number of units per FC hidden layer 20,50,100\20,50,100\ 100100 Table A1: Hyperparameter search ranges and selected values. All models are trained with Adam optimizer (24) for a maximum of 300 epochs, with an early stopping patience set to 15 epochs. Unless otherwise specified, we used a learning rate of 10−410^-4, a batch size of 256, and an L2 regularization of 10−510^-5. Data splitting is performed by random shuffling, allocating 60% for training, 10% for validation, and 30% for testing. Table A1 summarizes the candidate values considered for each hyperparameter. For CIHSI-Net, the final configuration are selected as the one that achieved the lowest validation loss on Simulation 1. Baseline Methods Selection and Implementation. To enable a comprehensive comparison, we implemented the following three baseline methods. These methods cover complementary modeling paradigms: TECE-VAE adopts a latent-variable generative approach, NCoRE explicitly models relations among treatment combinations, and CISI-Net employs pairwise Wasserstein balancing across treatment patterns. In particular, CISI-Net is the closest baseline to CIHSI-Net and enables us to directly evaluate the contribution of replacing pairwise balancing with the proposed barycenter-based BFG-WB objective. For each model, we set the hyperparameters as follows. • TECE-VAE: The latent dimension is set to 25, and the task embedding network has three hidden layers with 200 units and ELU activation (12). The task embedding network produces a five-dimensional embedding vector. • NCoRE: Each interaction subnetwork is implemented as two FC layers with 200 units per layer and ReLU activation (38). The representation learning network has two hidden FC layers with 200 units and outputs a 100-dimensional vector. • CISI-Net: CISI-Net has three neural networks, all of which are built using FC layers and leaky ReLU activation. Each network consists of two hidden layers with 100 units each. The task embedding network outputs a five-dimensional embedding vector. The balancing penalty coefficient α is set to 0.1, and the Integral Probability Metrics (IPM) used in the penalty is the 1-Wasserstein distance (48). All baseline models are tuned and trained under conditions comparable to CIHSI-Net to ensure a fair evaluation. We excluded other methods to ensure a focused and fair comparison for the following reasons: • Some methods have been empirically shown to underperform the baselines adopted in this study. In particular, methods originally developed for single-treatment settings, such as Treatment-Agnostic Representation Network (TARNet) (47) and Counterfactual Regression (CFR) (47), have been shown in previous experiments to be outperformed by methods specifically designed for multiple simultaneous treatments (36; 44). • Some methods rely on additional assumptions that are not adopted in our problem setting, such as the sequential ignorability assumption used in Single-cause Perturbation (41). • Some methods target different estimands or data structures, such as H-learner, which considers multiple outcomes (7). • Some methods share modeling paradigms and limitations similar to those already represented by the selected baselines, and therefore provide limited additional insight; Variational Sample Re-weighting (60), for example, represents a latent-variable approach similar to TECE-VAE. • Some methods employ architectures closely related to those of the selected baselines but are not applicable to binary treatments. For example, the Dosage Combination Network (45) is architecturally similar to CISI-Net but is specifically designed for combinations of continuous treatment dosages and therefore cannot be directly applied to the binary treatment setting considered here. Furthermore, conventional inverse probability weighting based and doubly robust estimators, are not included because previous empirical comparisons have shown them to be outperformed by TARNet or CFR (47). Hardware and Software Environments All experiments are conducted in a single-GPU computing environment running Ubuntu 22.04.5 LTS on an x86_64 architecture. The environment are equipped with one NVIDIA L4 GPU with 23 GB of GPU memory, 12 virtual CPUs corresponding to an Intel Xeon processor operating at 2.20 GHz, and 52 GB of system memory. CUDA compatibility version is 13.0. All methods are implemented in Python 3.12, and the neural network models are implemented primarily using TensorFlow 2.20.0. Appendix D Additional Simulation Results D.1 Comparison of Computing Speed Method K=2K=2 K=3K=3 K=4K=4 K=5K=5 K=6K=6 K=8K=8 CISI-Net 5.98 15.08 31.09 138.54 125.62 1529.83 CIHSI-Net 8.21 9.02 13.92 22.71 23.01 48.32 Table A2: Wall-clock training time per epoch (s/epoch) for CISI-Net (pairwise balancing) and CIHSI-Net (BFG-WB), evaluated with K∈2,3,4,5,6,8K∈\2,3,4,5,6,8\ treatments under the same hardware and training configuration. For each K, both methods are trained with the same batch size, while the batch size is adjusted across different values of K. Times are averaged across all observed epochs in 100 trials with early stopping (max 300 epochs, patience = 15). One advantage of implementing BFG-WB within CIHSI-Net is its lower computational cost relative to existing pairwise balancing methods. OT-based discrepancy evaluation, as used in representation balancing, requires solving optimization problems. This requirement often introduces a significant bottleneck during model training (19; 57). If the number of treatments is denoted as K, the total number of possible treatment patterns is L=2KL=2^K. The number of OT computations required for discrepancy evaluation in pairwise balancing increases quadratically with the number of treatment patterns, specifically following the combination formula (2K2)=O(L2)2^K 2=O(L^2). In contrast, BFG-WB reduces the required number of OT computations to O(L)O(L) and demonstrates a linear increase in computation count relative to the number of treatment patterns. To determine whether this theoretical reduction in computation translates into a measurable decrease in actual training time, we compare the required training time. We compare computational costs using the CISI-Net and CIHSI-Net with simulation dataset 1 (detailed in Appendix B). The number of treatments K is set to K∈2,3,4,5,6,8K∈\2,3,4,5,6,8\ to verify the scaling of training time as the number of treatment patterns increases to 2K2^K. For each setting, we conduct 100 independent trials using different random seeds for data generation. Each trial follows the hyperparameters and training protocol described in Appendix C. To ensure that each mini-batch contains at least four samples from each treatment pattern, we adjust the batch size according to K: 256 for K∈2,3,4,5K∈\2,3,4,5\, 1024 for K=6K=6, and 2048 for K=8K=8. We compute the training time per epoch (s/epoch) from the wall-clock measurements and compare the actual training time required for each method by averaging over all epochs across 100 trials. Table A2 summarizes the results. Table A2 demonstrates that the computational advantage of CIHSI-Net becomes clear as the number of treatments increases. For small treatment settings such as K=2K=2, this advantage is limited, and CIHSI-Net is slightly slower than CISI-Net due to the overhead of computing the barycenter-based discrepancy. However, from K=3K=3 onward, CIHSI-Net consistently achieves shorter training times than CISI-Net. The gap generally becomes more pronounced as K increases. For example, when K=5K=5, CISI-Net requires 138.54 [s/epoch], whereas CIHSI-Net requires only 22.71 [s/epoch]. Although the wall-clock time of CISI-Net decreases slightly from K=5K=5 to K=6K=6 because the batch size increases from 256 to 1024, thereby reducing the number of mini-batch updates per epoch, CIHSI-Net remains substantially faster in both settings. At K=8K=8, this computational advantage becomes particularly pronounced: CISI-Net requires 1529.83 [s/epoch], whereas CIHSI-Net requires only 48.32 [s/epoch]. These results suggest that CIHSI-Net offers practical advantages in training time, particularly in settings with an increased number of treatments. D.2 Scalability to a Larger Number of Treatments We investigate the scalability of the proposed method by evaluating whether its estimation accuracy remains superior to that of existing multiple-treatment methods as the treatment dimension increases. Multiple-treatment settings are practically relevant because many real-world studies involve the simultaneous evaluation of multiple treatments, which requires the estimation of numerous causal effects (40; 18). Therefore, we compare the estimation accuracy of each method under different values of K. In this scalability analysis, we set the number of treatments to K∈2,4,6,8K∈\2,4,6,8\. This range is chosen to reflect realistic multiple-treatment settings, since empirical applications rarely involve more than ten simultaneous treatments (10; 18; 37; 51). Therefore, setting the maximum number of treatments to K=8K=8 provides a practically meaningful evaluation range while still requiring the methods to handle a rapidly increasing number of treatment patterns. Across all settings, we use a common functional form for treatment assignment and outcome generation, and fix the sample size at N=50,000N=50,000. The covariates, treatment assignments, and outcomes are generated according to the following data-generating process: xi,n(j)∼N(cn(j),12),xi,u(j)∼U(−1,1),i=(i,n,i,u), x_i,n^(j) N(c_n^(j),1^2), x_i,u^(j) U(-1,1), x_i=( x_i,n, x_i,u), Hi(v)=(xi(v)+xi(v+1)>1),i=(Hi(1),Hi(2),Hi(3)), H_i^(v)=I (\,x_i^(v)+x_i^(v+1)>1\, ), H_i=(H_i^(1),H_i^(2),H_i^(3)), ti(k)∼Bern(σ(tk⊤i−H⊤i)),i=(ti(1),…,ti(K)), t_i^(k) (σ ( w_t_k x_i- w_H H_i ) ), t_i=(t_i^(1),…,t_i^(K)), Yi∼(g(i,i), 0.12), Y_i \! (g\! ( x_i, t_i ),\,0.1^2 ), where j∈1,…,15j∈\1,…,15\ denotes the covariate index, v∈1,2,3v∈\1,2,3\, and k∈1,…,Kk∈\1,…,K\. (⋅)I(·) is the indicator function, σ(⋅)σ(·) is the sigmoid function defined as σ(x)=1/(1+exp(−x))σ(x)=1/(1+ (-x)). The elements of tk w_t_k are independently drawn from U(−0.5,0.5)U(-0.5,0.5), and the elements of H w_H are independently drawn from U(−1,1)U(-1,1). The outcome generating function g is defined as follows: g(i,i)=x⊤i+∑∅≠S⊆1,…,K|S|≤KτS(i)∏k∈Sti(k)+2,τS(i)=wS(1+xi(jS)), g( x_i, t_i)= w_x x_i+ _ subarrayc ≠ S \1,…,K\\\ |S|≤ K subarray _S( x_i) _k∈ St_i^(k)+2,\ \ _S( x_i)=w_S (1+x_i^(j_S) ), where jS∈1,…,30j_S∈\1,…,30\ denotes a covariate index randomly assigned to each subset S. The vector x w_x is a weight vector whose elements are independently drawn from the uniform distribution U(−1,1)U(-1,1), and each wSw_S is independently drawn from U(0.2,1.0)U(0.2,1.0) with a randomly assigned sign. This setup induces treatment assignment bias through i H_i and heterogeneous treatment effects through τS(i) _S( x_i). All model hyperparameters are kept identical to those described in Appendix C. The only exception is the batch size, which is adjusted depending on K. Specifically, to ensure that each batch contains at least four samples from each treatment pattern, we set the batch size to 256 for K∈2,4K∈\2,4\, 1024 for K=6K=6, and 2048 for K=8K=8. This adjustment is necessary to obtain stable empirical estimates of the balancing regularization terms, because the number of treatment patterns increases as 2K2^K. We exclude TECE-VAE from this scalability analysis because it does not show competitive performance in the main simulation results reported in Section 5. Figure A1: Average ϵCASE _CASE (left) and ϵCAIE _CAIE (right) under varying numbers of treatments K∈2,4,6,8K∈\2,4,6,8\. Figure A1 shows the average ϵCASE _CASE and ϵCAIE _CAIE over 50 trials across varying numbers of treatments K∈2,4,6,8K∈\2,4,6,8\. CIHSI-Net consistently achieves the lowest estimation errors for both CASE and CAIE across all values of K, which demonstrates that the proposed method maintains high estimation accuracy even as K increases. Crucially, the performance gap between CIHSI-Net and CISI-Net becomes larger as K increases, and CIHSI-Net achieves this higher accuracy with substantially shorter training time than CISI-Net, as shown in Table A2. This tendency suggests that pairwise balancing becomes unstable as the number of treatment patterns grows, whereas the barycenter-based alignment in CIHSI-Net provides a more consistent global reference across treatment patterns. These results demonstrate that CIHSI-Net improves both estimation accuracy and computational efficiency in settings with many treatment patterns, supporting its scalability for multiple-treatment causal inference. D.3 Detailed Results of Ablation Study In this section, we provide the detailed values for ϵCASE _CASE and ϵCAIE _CAIE. Although the average results in Section 5 clarify the overall tendency, the detailed results further confirm whether the observed improvement is consistent across individual CASE and CAIE estimation errors. Table A3 shows the detailed ablation results on simulation dataset 1. The detailed results also show that the combination of the Wasserstein barycenter and the FGW discrepancy achieves the lowest estimation errors across all reported CASE and CAIE errors. These results are consistent with the discussion in Section 5 and provide additional support for the design of BFG-WB. ϵCASE _CASE ϵCAIE _CAIE No. WB Penalty type k=1k=1 k=2k=2 k=3k=3 S=1,2S=\1,2\ S=2,3S=\2,3\ S=1,3S=\1,3\ S=1,2,3S=\1,2,3\ 1 ✗ None 0.28 ± 0.08 0.29 ± 0.04 0.25 ± 0.09 0.36 ± 0.03 0.12 ± 0.04 0.16 ± 0.03 0.31 ± 0.04 2 ✗ W 0.22 ± 0.04 0.24 ± 0.03 0.19 ± 0.06 0.31 ± 0.07 0.11 ± 0.04 0.15 ± 0.08 0.31 ± 0.08 3 ✗ GW 0.25 ± 0.03 0.28 ± 0.04 0.23 ± 0.10 0.35 ± 0.03 0.12 ± 0.06 0.16 ± 0.04 0.32 ± 0.05 4 ✗ FGW 0.23 ± 0.04 0.25 ± 0.03 0.19 ± 0.05 0.33 ± 0.04 0.12 ± 0.06 0.16 ± 0.04 0.31 ± 0.04 5 ✓ W 0.22 ± 0.03 0.24 ± 0.03 0.19 ± 0.08 0.30 ± 0.04 0.10 ± 0.05 0.15 ± 0.09 0.28 ± 0.07 6 ✓ GW 0.22 ± 0.04 0.24 ± 0.03 0.19 ± 0.08 0.32 ± 0.06 0.10 ± 0.04 0.16 ± 0.05 0.29 ± 0.06 7 ✓ FGW 0.19 ± 0.04 0.21 ± 0.04 0.18 ± 0.10 0.28 ± 0.05 0.08 ± 0.05 0.12 ± 0.04 0.24 ± 0.04 Table A3: Ablation study result on simulation dataset 1. WB indicates whether the Wasserstein barycenter is used. Penalties: Wasserstein (W), Gromov-Wasserstein (GW), Fused GW (FGW), and None. D.4 Analysis of Representation Discrepancy We investigate the effectiveness of the proposed method in reducing distributional discrepancies across treatment patterns compared to existing approaches. Minimizing these discrepancies serves two critical objectives essential for accurate causal inference: (i) reducing selection bias arising from treatment assignment imbalances, which corresponds to minimizing distributional discrepancy (assessed via Wasserstein distance), and (i) preserving local proximity structures to stabilize the estimation of heterogeneous causal effects, which corresponds to maintaining geometric consistency (assessed via Gromov-Wasserstein distance). To empirically verify these properties, we compute these two metrics (the Wasserstein distance and the Gromov-Wasserstein distance) for every unique pair of treatment patterns on the dataset of Simulation 1 described in Section 5. We compare three methods: CISI-Net (without balancing), CISI-Net (with pairwise balancing), and our proposed CIHSI-Net (with BFG-WB). Figure A2 displays the boxplots of pairwise Wasserstein and Gromov-Wasserstein distances calculated across all treatment pairs. As shown in the left panel (Wasserstein distance), CIHSI-Net consistently achieves the lowest median and variance of discrepancies across representation distributions compared to the baselines. Whereas standard CISI-Net reduces discrepancies relative to the non-balanced model, it fails to match the compactness achieved by CIHSI-Net. This result suggests that simply extending pairwise balancing is insufficient for achieving comprehensive alignment in multiple-treatment settings. The right panel (Gromov-Wasserstein distance) further illustrates the advantage of our CIHSI-Net. CIHSI-Net exhibits lower structural discrepancies than CISI-Net. This result indicates that independent pairwise alignment can distort the local proximity structures within the representation space. In contrast, by aligning all distributions toward a shared Wasserstein barycenter, CIHSI-Net effectively mitigates both distributional divergence and local structural inconsistency. These results empirically demonstrate the superiority of CIHSI-Net in handling multiple treatments. Figure A2: Distribution of pairwise Wasserstein (left) and Gromov-Wasserstein distances (right) between representation distributions for all unique treatment pairs. We compare CISI-Net (without balancing), CISI-Net (with pairwise balancing), and CIHSI-Net (ours). Lower values indicate better alignment across treatment patterns. D.5 Sensitivity Analysis for α and η To examine the sensitivity of CIHSI-Net to the key hyperparameters in BFG-WB, we conduct a sensitivity analysis over the regularization strength α and the FGW trade-off parameter η. This analysis clarifies how the balance between global alignment and local structure preservation affects estimation accuracy, and assesses whether CIHSI-Net remains stable across a reasonable range of hyperparameter choices. Varying α∈0.1,0.5,1.0,1.5,2.0,5.0α∈\0.1,0.5,1.0,1.5,2.0,5.0\, we cover weak to strong regularization via BFG-WB, which examines the trade-off between representation balancing and outcome prediction accuracy. Additionally, to compare the relative contributions of local structure preservation and feature-wise distribution alignment, we vary η∈10−5,0.1,0.2,…,0.9,1.0η∈\10^-5,0.1,0.2,…,0.9,1.0\. Here, larger η places more weight on feature-wise Wasserstein alignment, whereas smaller η emphasizes the Gromov-Wasserstein-based structural term. We include η=10−5η=10^-5 as a near-boundary configuration that approximates GW-only alignment while remaining within the theoretical domain η∈(0,1]η∈(0,1]. Figure A3 summarizes the sensitivity results over α and η, averaged over 100 trials. Figure A3: Sensitivity analysis of hyperparameters α and η on simulation dataset 1. The top row shows ϵCASE(1) _CASE(1) - ϵCASE(3) _CASE(3), the middle row shows ϵCAIE(1,2) _CAIE(\1,2\) - ϵCAIE(1,3) _CAIE(\1,3\), and the bottom panel shows ϵCAIE(1,2,3) _CAIE(\1,2,3\). The results indicate that estimation performance depends on both α and η, which emphasizes the importance of appropriately selecting their combination. For 0.5<η<0.80.5<η<0.8, ϵCASE _CASE and ϵCAIE _CAIE remain small, suggesting that a moderate emphasis on feature-based alignment while retaining structural preservation is beneficial for causal effect estimation. At the near-boundary setting η=10−5η=10^-5, the estimation errors remain competitive for some estimands but are not consistently minimized across metrics and values of α. This pattern is consistent with the theoretical role of η: because the error bounds contain a factor of 1/η1/η, they become increasingly loose as η approaches zero and cease to provide finite control at the GW-only boundary η=0η=0. With respect to α, accurate estimation of both CASE and CAIE is achieved when 1.0≤α≤1.51.0≤α≤ 1.5. Weak balancing regularization, such as α=0.1α=0.1, tends to result in insufficient correction of selection bias and larger estimation errors across all metrics. Conversely, at α=5.0α=5.0, CASE estimation errors increase, while CAIE estimation errors are reduced. This may be interpreted as stronger selection bias affecting units receiving multiple treatments simultaneously, where strong balancing effectively estimated interaction treatment effects. D.6 Evaluation on Semi-Synthetic Dataset We conduct an additional semi-synthetic experiment to evaluate the proposed method under a more realistic setting while retaining access to ground-truth causal effects. Whereas the fully synthetic simulations allow us to evaluate the proposed method under controlled data-generating processes, they may not fully reflect the covariate distribution and treatment imbalance observed in real-world data. The purpose of this experiment is to verify whether CIHSI-Net can accurately estimate both CASE and CAIE when the covariates and treatment assignments follow an empirical distribution, while the outcome-generation mechanism is known. We construct a semi-synthetic dataset using the real covariates and observed treatment assignments from Dataset A used in Section 6. Specifically, we pair Dataset A’s real covariates and treatment assignments with a known synthetic outcome function that includes nonlinear heterogeneous effects for both CASE and CAIE. Let i∈ℝ71 x_i ^71 denote the observed covariates for unit i. The observed treatment vector i=(ti(1),ti(2),ti(3))∈0,13 t_i=(t_i^(1),t_i^(2),t_i^(3))∈\0,1\^3 is kept fixed from the original data. The observed outcome is then generated from the following semi-synthetic data-generating process: yi=g(i,i)+εi,εi∼(0,0.52).y_i=g( x_i, t_i)+ _i, _i (0,0.5^2). The outcome generating function g is defined as: g(i,)=tanh(i⊤g+bg)+∑∅≠S⊆1,2,3gS(i)∏k∈Stk,g( x_i, t)= ( x_i W_g+b_g)+ _ ≠ S \1,2,3\g_S( x_i) _k∈ St_k, where, for each non-empty subset S⊆1,2,3S \1,2,3\, gS(i)=wStanh(i⊤S+bS)g_S( x_i)=w_S ( x_i W_S+b_S) represents the heterogeneous treatment effect function. The elements of g W_g and S W_S are independently drawn from (0,12)N(0,1^2), and the bias terms bgb_g and bSb_S are drawn from (0,0.52)N(0,0.5^2). The coefficient wSw_S is independently drawn from U(−1,1)U(-1,1) for each non-empty subset S. This construction preserves the covariate distribution and treatment imbalance of the original data while providing known ground-truth CASE and CAIE values. We use the same model architectures, hyperparameters, and training protocol as those described in Appendix C. The data was randomly split into 60% training, 10% validation, and 30% test sets, with all reported evaluations performed on the test set. To ensure statistical reliability, we report average metrics over 100 independent runs. ϵCASE _CASE ϵCAIE _CAIE Method k=1k=1 k=2k=2 k=3k=3 S=1,2S=\1,2\ S=2,3S=\2,3\ S=1,3S=\1,3\ S=1,2,3S=\1,2,3\ TECE-VAE 1.23 ± 0.09 1.25 ± 0.06 1.28 ± 0.07 3.22 ± 0.46 1.82 ± 0.08 1.62 ± 0.17 4.00 ± 0.34 NCoRE 0.98 ± 0.16 0.98 ± 0.07 0.78 ± 0.05 2.94 ± 0.48 1.29 ± 0.28 1.37 ± 0.14 5.73 ± 1.20 CISI-Net 0.79 ± 0.01 1.13 ± 0.04 0.96 ± 0.12 1.49 ± 0.14 1.21 ± 0.21 0.96 ± 0.02 1.49 ± 0.11 CIHSI-Net 0.81 ± 0.02 1.02 ± 0.07 0.78 ± 0.07 1.45 ± 0.07 1.09 ± 0.08 0.96 ± 0.04 1.45 ± 0.06 Table A4: Semi-synthetic experiment results using Dataset A’s real covariates and observed treatment assignments. Bold values indicate the best performance, and underlined values indicate the second-best performance for each metric. Table A4 shows the results of the semi-synthetic experiment. Overall, CIHSI-Net achieves the strongest performance among the baseline methods under the empirical covariate distribution and observed treatment assignments of Dataset A. Although CISI-Net and NCoRE obtain the lowest error for some individual ϵCASE _CASE, CIHSI-Net remains competitive for CASE and achieves the best results for all ϵCAIE _CAIE. When the errors are averaged over the ϵCASE _CASE and ϵCAIE _CAIE, CIHSI-Net obtains the lowest overall errors, suggesting that the proposed barycenter-based balancing is effective for estimating both single and interaction effects in the semi-synthetic setting. This result suggests that CIHSI-Net remains robust even when the treatment assignment imbalance is inherited from real-world data. Appendix E Real-world Datasets and Additional Results In this section, we provide supplementary details and additional results for the real-world marketing application. We first describe the dataset specifications, preprocessing steps, and experimental protocols common to both Dataset A and Dataset B. Subsequently, we present the analysis results for Dataset B, which involves promotions from competing merchants. E.1 Details of Real-World Application Dataset Description and Treatment Assignment. The dataset is derived from a mobile payment platform and involves multiple concurrent marketing promotions. We utilize this proprietary dataset because existing publicly available causal inference benchmarks focus primarily on single-treatment scenarios and lack the complex, simultaneous multiple-treatment interactions observed in real-world marketing environments. • Dataset A (reported in Section 6): This dataset includes three marketing promotions that were simultaneously conducted: two offline promotions organized by the same merchant group (denoted as CP1 and CP2) and one online promotion conducted by another merchant group (denoted as CP3). Each promotion is represented by a binary treatment indicator, which results in 23=82^3=8 distinct treatment patterns. The observed sample proportions were approximately 2% for CP1 only, 15% for CP2 only, 42% for CP3 only, 0.2–1% for two-promotion combinations, 0.2% for all three promotions simultaneously, and the remainder as the control group. • Dataset B (reported in this appendix): This dataset includes two marketing promotions conducted by two different merchants in the same industry (denoted as CP4 and CP5), yielding 22=42^2=4 treatment patterns. The observed sample proportions were approximately 29% for CP4 only, 4% for CP5 only, 1% for both promotions, and the remainder as the control group. For both datasets, the control group consists of users who were not exposed to any of the promotion during the treatment period and who satisfied the positivity assumption, randomly drawn from users with at least one mobile payment transaction in the month preceding the promotions. Additionally, the preprocessing and experimental protocol described below are applied identically to both datasets. Outcome Definition and Preprocessing. The outcome variable Y is defined as the total payment amount during the one-month period following the promotion implementation. To ensure numerical stability and protect data confidentiality, Y was standardized to have zero mean and unit variance. This preprocessing preserves the sign of the estimated causal effect and ensures interpretational consistency between the standardized scale and the original scale because the standardization is an affine transformation (49; 36). The covariates consist of 71 variables, including service usage histories and user demographic attributes. User Stratification for Heterogeneity Analysis. To analyze heterogeneous treatment effects, we stratified users based on their total payment amount in the month before the promotions. Users were grouped into intervals of 5,000 JPY (e.g., 0–5,000, 5,000–10,000). Users with usage exceeding 50,000 JPY were aggregated into a single high-usage group, resulting in a total of 11 distinct strata. This grouping strategy allows us to examine how the sensitivity to marketing incentives varies with prior engagement levels. Experimental Configuration. For the real-world application, we employed the same hyperparameter configuration that achieved the best performance in the simulation experiments (Section 5). This approach was chosen to evaluate the method’s robust performance under a reproducible setting without dataset-specific tuning. The dataset was randomly split into 60% training, 10% validation, and 30% test sets, with all reported evaluations performed on the test set. E.2 Sensitivity and Uncertainty Analysis on Dataset A Figure A4: Sensitivity and uncertainty analysis of the estimated causal effects in the real-world application on Dataset A. The rows correspond to the BFG-WB regularization strength α∈0.5,1.0,1.5α∈\0.5,1.0,1.5\, and the columns correspond to the FGW trade-off parameter η∈0.5,0.6,0.7η∈\0.5,0.6,0.7\. In each panel, the top subpanel reports the estimated CASE and the bottom subpanel reports the estimated CAIE across user groups. The shaded regions represent bootstrap confidence intervals. This appendix reports additional sensitivity and uncertainty analyses for the estimated causal effects in the real-world application on Dataset A. In real-world observational data, the ground-truth causal effects are unavailable, and the reliability of the results cannot be assessed through direct estimation errors as in the simulation studies. Therefore, it is important to examine whether the estimated causal-effect patterns remain stable under different model configurations and sampling variability. We evaluate the sensitivity of the estimated CASE and CAIE to key BFG-WB hyperparameters and compute bootstrap confidence intervals for the estimated causal effects. These analyses provide complementary evidence on the robustness of the real-world application. For the sensitivity analysis, we focus on two key hyperparameters related to BFG-WB: the regularization strength α in Equation (10) and the FGW trade-off parameter η in Equation (6). We evaluate nine hyperparameter configurations by varying α∈0.5,1.0,1.5α∈\0.5,1.0,1.5\ and η∈0.5,0.6,0.7η∈\0.5,0.6,0.7\ while keeping the other model configurations fixed. For each configuration, the data are randomly split into 60% training, 10% validation, and 30% test sets, and all reported evaluations are performed on the test set. To quantify uncertainty, we compute bootstrap confidence intervals by repeating the entire estimation pipeline for each hyperparameter configuration (16; 25). We repeat this procedure 1,000 times for each hyperparameter configuration. We then aggregate the test-set estimates for each bootstrap replication and construct confidence intervals from the empirical distribution of the resulting aggregate statistics. Figure A4 shows the estimated CASE and CAIE with bootstrap confidence intervals under the nine hyperparameter configurations of α∈0.5,1.0,1.5α∈\0.5,1.0,1.5\ and η∈0.5,0.6,0.7η∈\0.5,0.6,0.7\. Although some individual confidence intervals include zero, the overall treatment-effect patterns remain stable across different hyperparameter settings. In particular, the estimated effects preserve the same qualitative trends observed in the main real-world analysis, such as the sign reversal of the three-way interaction effect from negative for low-usage users to positive for high-usage users. This suggests that the empirical findings are not driven by a specific choice of the BFG-WB regularization strength or the FGW trade-off parameter. Moreover, the consistency of the estimated CASE and CAIE patterns across the bootstrap replications suggests that CIHSI-Net captures the primary interaction structures in the real-world application. These results provide additional evidence that the real-world findings are robust to both hyperparameter variation and sampling variability. E.3 Additional Results on Dataset B Figure A5: In dataset B, estimated CASEs for CP4 and CP5 (left) and the estimated CAIE for CP4,CP5\CP 4,CP 5\ (right) across user groups stratified by pre-promotion service usage. The outcomes are standardized. Figure A5 shows the estimated causal effects of dataset B on a standardized outcome scale. The left panel in Figure A5 shows that both estimated single-treatment effects are positive across all user groups. The estimated effect of CP5 is consistently larger than that of CP4, and for both promotions, the effect sizes gradually decrease as the pre-promotion usage level increases. This pattern suggests that CP4 and CP5 are effective at inducing behavioral changes among low-usage users, whereas their incremental effects are limited for high-usage users whose service usage is already established. This heterogeneity is consistent with established findings in marketing research (5; 26; 34, e.g.), which suggests that light users typically possess a greater capacity to increase their transaction volume, whereas heavy users are often constrained by saturation in their activity levels. Collectively, the contrast with Dataset A (see Section 6) highlights that treatment effect heterogeneity is highly context-dependent, varying significantly based on promotion characteristics and competitive dynamics. The right panel in Figure A5 shows that the estimated interaction effect τCAIE(4,5) _CAIE(\4,5\) is consistently negative across all groups, which suggests that promotions offered by competing merchants in the same industry may induce a reallocation of spending across merchants (cannibalization). The results of this cannibalization are consistent with previous research (15). In contrast to the CASE results, the absolute magnitude of the negative interaction becomes smaller for higher-usage groups, indicating that the degree of cannibalization is mitigated as usage increases. This trend suggests that cannibalization may be weaker among high-usage users because their behavior of using multiple merchants is already established. Overall, these additional results suggest that CIHSI-Net serves as a practical analysis tool by flexibly uncovering heterogeneous single and interaction effects across user groups.