: Distributional Counterfactual Explanation With Optimal Transport

Lei Youdtu
Lele Caoeqt
Mattias Nilssonneko

dtuDepartment of Engineering Technology, Technical University of Denmark, Ballerup, Denmark.
eqtMotherbrain, EQT, Stockholm, Sweden
nekoNekoHealth, Stockholm, Sweden

Lei<EMAIL_ADDRESS>
Machine Learning, ICML

Counterfactual Explanations (CE) is the de facto method for providing insight and interpretability in black-box decision-making models by identifying alternative input instances that lead to different outcomes. This paper extends the concept of CEs to a distributional context, broadening the scope from individual data points to entire input and output distributions, named Distributional Counterfactual Explanation (DCE). In DCE, our focus shifts to analyzing the distributional properties of the factual and counterfactual, drawing parallels to the classical approach of assessing individual instances and their resulting decisions. We leverage Optimal Transport (OT) to frame a chance-constrained optimization problem, aiming to derive a counterfactual distribution that closely aligns with its factual counterpart, substantiated by statistical confidence. Our proposed optimization method, , strategically balances this confidence across both input and output distributions. This algorithm is accompanied by an analysis of its convergence rate. The efficacy of our proposed method is substantiated through a series of illustrative case studies, highlighting its potential in providing deep insights into decision-making models.

§ BACKGROUND

In the field of xai, ce has become the quintessential method for providing insight and interpretability in complex decision-making models. This prominence is particularly marked in domains where understanding the causal impact of variables is pivotal for informed decision-making and adherence to ethical standards [Verma et al., 2020]. The utility of ce is grounded in their ability to answer “what-if” scenarios, offering tangible insights into how slight alterations in input can lead to different model outputs. This not only aids in demystifying the decision process of black-box models but also helps in identifying potential biases and areas for model improvement.

Motivation. ce traditionally focuses on individual data points, examining how changes in specific variables can influence the model's output. While this method is beneficial, it neglects scenarios where predictions are generated for several comparable instances yielding a certain pattern of outcomes, a phenomenon highlighted in the work [Warren et al., 2023]. One example is the ride-hailing pricing models, aiming at making proper price to balance the demand and supply in the market meanwhile optimizing multiple metrics such as order completion rates and revenue etc [Yan et al., 2020]. For instance, consider a scenario where a stakeholder inquires, “What changes in the distribution of pick-up distances could potentially increase the order completion rate by 2% for a specific customer group?” Solving this query involves identifying a distributional counterpart for the group in question.

Addressing this gap, our work introduces the concept of dce, extending the conventional scope of ce by considering the entire spectrum of input and output distributions, rather than limiting the analysis to isolated instances. This shift enables a more comprehensive understanding of the model's behavior, which is particularly relevant in scenarios where the model's decisions are influenced by the distributional properties of the input rather than just individual data points.

Related Work. The pioneering study by Wachter et al. [Wachter et al., 2017] is a seminal work in the field of ce and has gained notable recognition. This study approaches ce through a minimization problem, laying the foundation for subsequent explorations in the field.
A comprehensive overview of the CE literature is presented in the works of Guidotti et al. [Guidotti, 2022] and Bodria et al. [Bodria et al., 2023]. Our focus aligns with the model-agnostic optimization approach for generating counterfactuals as proposed by Wachter et al. [Wachter et al., 2017]. Several studies have explored the generation of multiple counterfactuals in one-shot [Ustun et al., 2019, Karimi et al., 2020, Mothilal et al., 2020, Pawelczyk et al., 2020, Cheng et al., 2020, Brughmans et al., 2023], providing insights into the explanation of individual factual instances.

Recent developments in the field have seen a growing interest in group-based ce. Rawal et al. [Rawal & Lakkaraju, 2020] introduced a framework for generating ce for sub-populations within datasets. Plumb et al. [Plumb et al., 2020] proposed a counterfactual-summary method for explaining subgroup data points. Furthermore, Warren et al. [Warren et al., 2023] developed a novel approach for creating gce, focusing on common feature differences among similar individuals. [Carrizosa et al., 2024] takes a stakeholder's point of view to explain the necessity of considering GCE.

Despite these advancements, a notable gap remains in bridging these concepts to the realm of dce, discussed below.

Challenges and the Gap. The conceptual leap from the classic ce to dce represents a significant paradigm shift. Unlike traditional approaches that generate multiple counterfactual instances, dce treats these data points as part of an empirical distribution, capturing their inherent statistical properties. This distinction is crucial; in classic methods, each data point is optimized independently (except for gce), resulting in a set of isolated instances, while dce considers these points as an entire group, on which the optimization is performed while maintaining certain statistical properties of the group. Therefore, the optimization is point-wisely dependent. It is important to note a fundamental distinction between dce and GCE. GCE primarily focuses on ensuring a direct match between each generated counterfactual and its corresponding factual instance. In contrast, DCE operates at the distributional level, where a specific data point in the counterfactual distribution may not necessarily correspond to any individual factual data point at the instance level. This approach is more aligned with practical real-world scenarios, where resources may not suffice to facilitate individual-level shifts for each person.

Although metrics exist to quantify the extent of difference between distributions, interpreting what these differences mean in practical terms is not always clear. In the realm of ce or its generalized gce, counterfactuals are identified by adjusting the values in each dimension of the vector representation of a data point. dce, however, extends this concept from the modification of “quantity” to the adaptation of “quantiles”, fitting well within the data science and analysis paradigm where quantile plots and histogram comparisons are commonly used to assess similarities or divergences between distributions.

In classical ce and gce, it is straightforward to determine the difference between a counterfactual and an observed data point, given a predefined threshold. However, dce faces the challenge of accounting for the sampling error, as the evaluation of the similarity between two distributions typically relies on sample data rather than the entire population. This reliance on sampling introduces additional complexity. The need for sampling in dce raises the issue of sample complexity and the curse of dimensionality, which needs to addressed to ensure accurate and computationally feasible counterfactual generation.

trig format plots=rad,
hide axis,

3 [surf, fill=white, domain=0:6, domain y=0:6, samples=30,
samples y=50]

3 [
only marks,
ball color=red,
mark size=0.6mm,
nodes near coords=$x'$
] coordinates (3, 3, 0.3);

3 [
only marks,
ball color=blue,
mark size=0.6mm,
nodes near coords=$x$
] coordinates (4, 3.5, 0.1);

3 [
point meta = 1,
z buffer=sort,

domain y=2.5:5

at (axis cs:6.3, 1, 1) [pin=[pin distance=1cm, pin edge=black, dashed, thick] below: Confidence] ;

at (axis cs:8.1, 1.2, 0) [pin=[pin distance=3cm, pin edge=black, dashed, thick] right: ] ;

axis lines=middle,
hide axis,
xmin=-1, xmax=3.6,
ymin=0, ymax=0.4,

legend style=draw=none, empty legend, at=(1.1, 1.1),
legend cell align=left

[blue, thick] 1/(1.1*sqrt(2*pi)) * exp(-((x-1)^2)/(2*0.5^2));

[BrickRed, thick] 1/(1.2*sqrt(2*pi)) * exp(-((x-1.2)^2)/(2*0.6^2));

[yellow!50!black, fill=LimeGreen, fill opacity=0.5] min(1/(1.1*sqrt(2*pi)) * exp(-((x-1)^2)/(2*0.5^2)), 1/(1.2*sqrt(2*pi)) * exp(-((x-1.2)^2)/(2*0.6^2))) ;

In this depiction, $x'$ represents the factual distribution. The objective of dce is to identify a counterfactual distribution $x$ on a manifold (that is embedded in a high-dimensional Euclidean space) that stays close to the factual within a pre-specified confidence region. The density plot below delineates the confidence that the model output $y$ ($y=b(x)$) stays close to the target $y^*$.

Approach Sketch.
We propose to utilize ot (and its 1d linear projection form known as the sliced version of ot) to evaluate the divergence between distributions via their quantile functions. This methodology is particularly effective in high-dimensional spaces because the sliced ot maintains the same sample complexity and offering a dimension-free convergence rate.

We conceptualize dce as a chance-constrained optimization problem. In this framework, the objective is to optimize for a distributional counterfactual that is statistically close to the factual distribution, while ensuring that the resulting model output distribution aligns closely with a predefined target distribution. We employ confidence intervals of both the Wasserstein distance and its sliced version, thereby providing statistical guarantee for the validity and proximity of the generated counterfactuals. Our optimization balances two critical factors: the similarity of the input distribution to the factual, and the necessity for the input to serve as a counterfactual. This trade-off is managed by optimizing these two components within a bcd framework on a manifold. See a depiction in  <ref>.

To cater to diverse requirements, we explore both discrete and continuous optimization strategies for performing the above mentioned balance. Furthermore, we have rigorously proven the convergence rates of our  approach in solving the bi-level optimization problem, for each of these strategies. Through this approach, our goal is not only  to provide a statistically sound method to generate distributional counterfactuals, but also to lay the foundation linking ot with ce specifically tailored for distributional contexts, serving as an initial foray into this research area. It is our aspiration that this preliminary exploration will catalyze further, more extensive investigations, leading to richer and more profound advancements in the field.

§ PRELIMENARIES AND PROBLEM SETUP

We first give a brief introduction of ot and its sliced version. Then we formulate the dce problem as a chance constrained optimization problem based on ot.

ot and the Wasserstein Distance. Denote $\P$ the set of probability measures with finite second moments. Let $\gamma_1, \gamma_2\in\P$ and let $\Pi_{\gamma_1, \gamma_2}$ denote the set of probability measures in $\P$ with marginal distributions equal to $\gamma_1$ and $\gamma_2$. We call $\Pi(\gamma_1, \gamma_2)$ the set of transport plans between $\gamma_1$ and $\gamma_2$. The 2-Wasserstein distance, i.e., ot, is defined as

\begin{equation}
\ot(\gamma_1, \gamma_2) \triangleq \inf_{\pi\in\Pi(a_1, a_2)}  \int_{\RR^{d\times d}} \norm{a_1-a_2}^2 \d{\pi(a_1,a_2)}
\end{equation}

Sliced Wasserstein Distance. Different from its original version, the sliced Wasserstein distance first performs 1d linear projections of two high-dimensional measures $\gamma_1$ and $\gamma_2$, then it computes the Wasserstein distance between the projected measures. Formally, let $\SS^{d-1}\subset\RR^d$ stand for the $d$-dimensional unit sphere:

\begin{equation}
\sot(\gamma_1, \gamma_2) \triangleq \int_{\SS^{d-1}} \ot(\theta \sharp \gamma_1, \theta \sharp \gamma_2)\d{\sigma(\theta)} \\
\end{equation}

where $\sigma$ is the uniform distribution on $\SS^{d-1}$ and $\sharp$ is the push-forward operator, which projects (high-dimensional) measrues $\gamma_1$ and $\gamma_2$ onto a 1d space.

dce Problem Formulation. Denote $b$ a black-box model. Given any data distribution $x'$ as a factual (i.e. an observed distribution) input, denote by $y'=b(x')$ the factual model output distribution. Denote by $y^*$ a target distribution, and $y^*\neq y'$.  dce tries to identify a counterfactual input distribution $x$ that is as close as possible to $x'$, so that the model output $b(x)$ is close enough to the target $y^*$, formulated as a chance-constrained optimization problem shown in (<ref>). The inequality (<ref>) imposes that the sliced Wasserstein distance between the high-dimensional $x'$ and its counterfactual distribution $x$ is subject to a chance constraint such that the distance is smaller than a threshold $U_{\mu}$ with a probability no lower than $\varphi$. The inequality (<ref>) is for the same goal but for the model output[This formulation generalizes straightforwardly to multi-dimensional output by using sliced wasserstein distance in (<ref>). The method and conclusions proposed in this paper still hold.]. We set a confidence level $\alpha/2$, such that a feasible solution needs to make both chance constraints hold beyond the level probabilistically.

\begin{align}
& \max_{x, \varphi} \varphi \\
\text{s.t.} \quad & \varphi \leq  \PP\left[\sot(x,x')\leq U_{\mu}\right] \label{eq:main_problem_C1} \\
& \varphi \leq  \PP\left[\ot(b(x), y^*)\leq U_{\nu}\right] \label{eq:main_problem_C2} \\
& \varphi \geq 1-\alpha/2
\end{align}

The 2-Wasserstein distance that is used in this paper is special in that its geometry admits a formal Riemannian structure. It is shown later in Section <ref> that such a property eases optimization.

§ THEORETICAL FOUNDATIONS OF DCE

We emphasize the formulation (<ref>) possess several important properties. We highlight that dce extends ce from quantity comparison to quantile comparison, followed by a discussion on the sample complexity of the slicing approximation.

Quantile-Based Expression. Let $y = b(x)$. The Wasserstein distance in (<ref>) can alternatively be written as follows.

\begin{align}
\ot(y, y^*) & = \int_{0}^{1} \norm{F^{-1}_{y}(q) - F^{-1}_{y^*}(q)}^2 \d{q} \label{eq:y-quantile}
\end{align}

In the equation above, $F^{-1}_y(q)$ and $F^{-1}_{y^*}(q)$ are the inverse cumulative distribution functions (i.e. quantile functions) of $y$ and $y^{*}$.
The equation holds due to both $y$ and $y^*$ being 1d distributions, providing an alternative explanation of the Wasserstein distance: It is the integral of the square of the difference between the corresponding quantiles in the interval $[0,1]$. Essentially, this aligns the model output distribution $y$ ($y=b(x)$) and the target distribution $y^*$ in a way that we are comparing the location of equivalent “percentiles” of mass in each distribution. In fact, it accounts for the full range of variability in the outputs, not just their central tendency (e.g. in contrast to comparing $\EE(y)$ and $\EE(y^*)$ only), which is crucial for many practical applications where distributional characteristics are important.

Sample Complexity of Slicing Approximation.
In the objective of formulation (<ref>), the distance between $x$ and $x'$ is evaluated using the sliced Wasserstein distance, an approximation of their true Wasserstein distance. In real-world applications, we typically have finite samples drawn from the underlying distributions of $x$ and $x'$. The effectiveness of this approximation hinges on its capacity to provide accurate distance estimates using a limited number of samples, a characteristic known as “good sample complexity”.  A natural question is whether the slicing approximation possesses the same sample complexity as its origin. By <cit.> directly, we have the following result.

Consider any measure $z\in\P(\RR)$ with an empirical measure $\vec{z}=\{z_i\}_{i=1}^n$ and measure $z'\in\P(\RR)$ with an empirical measure $\vec{z}'=\{z'_j\}_{j=1}^n$. Let $\gamma$ be a function of $n$. Suppose the Wasserstein distance $\ot$ admits the following sample complexity

\[
\EE[\ot(z, z') - \ot(\vec{z},\vec{z}')] \leq \gamma(n)
\]

Then for any measures $x, x'\in\P(\RR^d)$ with corresponding $\vec{x}=\{\vec{x}_i\}_{i=1}^n$ and $\vec{x}'=\{\vec{x}'_i\}_{i=1}^n$, the sample complexity is

\[
\EE[\sot(x, x') - \sot(\vec{x}, \vec{x}')] \leq \gamma(n).
\]

Proposition <ref> asserts that the sample complexity of the slicing approximation is directly proportional to that of its original form. The convergence rate is independent of dimensionality for $\vec{x}, \vec{x}'\in\RR^d$ with respect to $d$. Hence formulation (<ref>) scales to high-dimensional input well.

§ OPTIMIZATION FOR PROBLEM SOLVING

Preliminaries. Otto [Otto, 2001] showed that the space metrized by the 2-Wasserstein distance admits the structure of a formal Riemannian metric.  Let $\M$ be a smooth Riemannian manifold of $\RR^{d}$, closed and non-empty, where $d$ is the dimension of $\vec{x}$.

\begin{align}
Q_{\mu}(\vec{x}, \bm{\mu}) & \triangleq \sum_{\bm{\theta}\in\Theta}\sum_{i=1}^{n}\sum_{j=1}^{n}\abs{\bm{\theta}^\top\vec{x}_i - \bm{\theta}^{\top}\vec{x}'_j}^2\mu^{(\theta)}_{ij} \label{eq:Q_mu} \\
Q_{\nu}(\vec{x}, \bm{\nu}) & \triangleq   \sum_{i=1}^{n}\sum_{j=1}^{n}\abs{b(\vec{x}_i) - y^{*}_j}^2 \nu_{ij} \label{eq:Q_nu}%
\end{align}

Both $Q_{\mu}$ and $Q_{\nu}$ are smooth in $\vec{x}$, $\bm{\mu}$ and $\bm{\nu}$, respectively.
The empirical version of (<ref>) reads:

\begin{align}
& \max_{\vec{x}\in\M, \varphi\geq 0} \varphi \\
\text{s.t.} \quad & \varphi \leq   \PP_n \left[ \inf_{\bm{\mu}\in \Pi} Q_{\mu}(\vec{x},\bm{\mu})\leq U_{\mu} \right] \label{eq:main_problem_empirical-U_1}\\
& \varphi \leq   \PP_n \left[\inf_{\bm{\nu}\in \Pi} Q_{\nu}(\vec{x}, \bm{\nu})\leq U_{\nu}\right] \label{eq:main_problem_empirical-U_2}\\
& \varphi \geq 1-\alpha / 2 \label{eq:main_problem_empirical-alpha_2}
\end{align}

Initially, we establish a theorem to manage the chance constraints (<ref>) and (<ref>).
Subsequently, we demonstrate a partial optimality condition for resolving (<ref>), uncovering that achieving the optimum (or local optima) is contingent upon identifying an appropriate balance between $Q_{\mu}$ and $Q_{\nu}$. Lastly, we introduce our algorithm , which supports both discrete and continuous strategies for this balancing, and come with a guarantee of convergence. The discrete approach, termed Set Shrinking, optimizes the weights by selecting from a finite set of candidate values. In contrast, the continuous method, known as Interval Narrowing, utilizes a bisection strategy applied to a progressively narrowing domain of variables.

The convergence rate of integrating both strategies is elaborated in Section <ref>.

Theorem <ref> below, which follows <cit.>, yields ucl for $ \inf_{\bm{\mu}\in \Pi} Q_{\mu}(\vec{x},\bm{\mu})$ and $\inf_{\bm{\nu}\in \Pi} Q_{\nu}(\vec{x}, \bm{\nu})$ in (<ref>) and (<ref>). We let $F^{-1}$ be the quantile function, with its empirical version defined as $ F^{-1}_n(x)=\frac{1}{n}\sum_{i=1}^n I(X_i\leq x)$. For example, $F^{-1}_{x, n}$ is the empirical quantile function of the distribution $x$ with $n$ data points sampled from $x$. We also remark that $N$ stands for the size of the set $\Theta$, i.e. $N=|\Theta|$. Sketch proof is in Appendix <ref>.

Denote by $\downq_{\alpha, n, N}$ and  $\upq_{\alpha, n, N}$ some sequences of functions $(0,1)\rightarrow \RR$ that satisfies

\begin{align}
\inf \PP_n\big(F_n^{-1}(\downq_{\alpha, n, N}(u)) \leq & F^{-1}(u)\leq  F_n^{-1}(\upq_{\alpha, n, N}(u)), \nonumber \\
& \forall u\in(0,1) \big) \geq 1-\frac{\alpha}{2N}.
\label{eq:alpha/2N}
\end{align}

Let $\delta\in(0,1/2)$ be a trimming constant. The following inequality holds.

\begin{align}
\PP_n\left[\inf_{\bm{\nu}\in\Pi}Q_{\nu}(\vec{x},\bm{\nu}) \leq \frac{1}{1-2\delta} \int_{\delta}^{1-\delta} D_1(u)\d u \right]  \nonumber \\
\geq 1-\frac{\alpha}{2} \label{eq:int_D_1}
\end{align}

\begin{align}
D_1(u) \triangleq  & \max \big\{F_{y,n}^{-1}(\upq_{\alpha, n,1}(u)) - F_{y^{*},n}^{-1}(\downq_{\alpha, n,1}(u)),  \nonumber \\
& F_{y^{*},n}^{-1}(\upq_{\alpha, n,1}(u)) - F_{y,n}^{-1}(\downq_{\alpha, n,1}(u))\big\}. \label{eq:D_1}
\end{align}

Let $\theta_1\,\ldots,\theta_N$ be an i.i.d. sample from a distribution $\omega$. Let $\omega_{N}$ denote the corresponding empirical measure, then

\begin{align}
\PP_n\bigg[ \inf_{\bm{\mu}\in\Pi} Q_{\mu}(\vec{x}, \bm{\nu}) \leq \frac{1}{1-2\delta}\int_{\SS^{d-1}}\!\!\int_{\delta}^{1-\delta}\!\!\!\! & D_{\theta, N}(u)\d u \d\omega_N(\theta) \bigg] \nonumber \\
& \geq 1-\frac{\alpha}{2} \label{eq:int_D_N}
\end{align}

\begin{align}
D_{\theta,N}(u) \triangleq & \max\big\{ F_{\theta,x,n}^{-1}(\upq_{\alpha, n,N}(u)) - F_{\theta, x',n}^{-1}(\downq_{\alpha, n,N}(u)),  \nonumber \\
& F_{\theta, x',n}^{-1}\!(\upq_{\alpha, n,N}(u))\! - \!F_{\theta, x,n}^{-1}\!(\downq_{\alpha, n,N}(u))\big\}
\label{eq:D_N}
\end{align}

This theorem facilitates the optimization process while simultaneously ensure the chance constraints. The inequalities (<ref>) and (<ref>) provide methods to validate and substantiate the chance constraints (<ref>) and (<ref>), respectively. The right-hand-side of (<ref>) is no less than $1-\alpha/2$ if and only if we have the corresponding ucl no larger than $U_{\nu}$, i.e.

\begin{equation}
\overline{Q}_{\nu}\triangleq  \frac{1}{1-2\delta} \int_{\delta}^{1-\delta} D_1(u)\d u  \leq U_{\nu}.
\label{eq:overline_Q_nu}
\end{equation}

The computation of (<ref>) lies on $D_1(u)$, defined in (<ref>), which quantifies the disparity between the quantile functions of $x$ and $x'$ by utilizing pre-specified function sequences $\downq_{\alpha, n,N}:(0,1)\rightarrow\RR$ and $\upq_{\alpha, n,N}$. The critical task is to formulate suitable $\downq_{\alpha, n,N}$ and $\upq_{\alpha, n,N}$ functions that are grounded in statistical theory. The variable $u$ in $D_1(u)$ corresponds to the quantile level, while the sequences $\downq_{\alpha, n,N}$ and $\upq_{\alpha, n,N}$ serve to transform this quantile level $u$ into a real number. This transformation adheres to the principles of uniform quantile bounds, a topic extensively discussed in the work of Shorack & Wellner, 2009. One classic example is $\downq_{\alpha, n,N}=u-\beta_n$,  $\upq_{\alpha, n,N}=u+\beta_n$, and $\beta_n=\sqrt{\log(4/\alpha)/2n}$. Then (<ref>) holds by dkw inequalities [Dvoretzky et al., 1956]. One could refer to <cit.> for more examples. Similarly, the right-hand-side of (<ref>) is no less than $1-\alpha/2$ if and only if

\begin{equation}
\overline{Q}_{\mu}\triangleq \frac{1}{1-2\delta}\int_{\SS^{d-1}}\!\!\int_{\delta}^{1-\delta}\!\!D_{\theta, N}(u)\d u\d\omega_N(\theta)  \leq U_{\mu}
\label{eq:overline_Q_mu}
\end{equation}

with $D_{\theta,N}$ defined in (<ref>). Compared to (<ref>), the difference in (<ref>) is that we average $u_{N}(\theta)$ on $\theta$ that distributed on the $d$-dimensional unit sphere[Remark $Q_{\nu}$ is finite sample inference whereas $Q_{\mu}$ is not. One needs to enlarge $Q_{\mu}$ with a factor that is inversely proportional to $N=|\Theta|$ to obtain the finite version [Manole et al., 2022].].

Alternatively, the clt is applicable to both ot and its sliced counterpart, allowing $\overline{Q}_{\mu}$ and $\overline{Q}_{\nu}$ to be estimated via the bootstrap method. For further information, refer to the theorems presented in <cit.> and <cit.>. Regardless of the approach, both chance constraints in the dce formulation (<ref>) are rigorously respected.

Algorithm Design. Next, we investigate how to optimize $\vec{x}$ (which determines $u$ essentially in $\overline{Q}_{\mu}$ and $\overline{Q}_{\nu}$). Define

\begin{equation}
Q(\vec{x},\bm{\mu}, \bm{\nu} | \eta) \triangleq  (1-\eta)\cdot Q_{\mu}(\vec{x},\bm{\mu}) +  \eta \cdot Q_{\nu}(\vec{x}, \bm{\nu})
\label{eq:Q-eta}
\end{equation}

Suppose (<ref>) is feasible. There exists $\eta^*$ ($0\leq \eta^* \leq 1$) such that the optimum of $Q(\vec{x},\bm{\mu}, \bm{\nu} | \eta)$, i.e.,

\begin{equation}
\vec{x}^{*} \triangleq \argmin_{\vec{x},\bm{\mu},\bm{\nu}}  Q(\vec{x},\bm{\mu}, \bm{\nu} | \eta^*),
\label{eq:main_problem_empirical_alternative}
\end{equation}

is the optimum of (<ref>), namely,

\[
\vec{x}^* = \argmin_{\vec{x},\bm{\mu},\bm{\nu}} \varphi \text{ s.t. } \eqref{eq:main_problem_empirical-U_1}-\eqref{eq:main_problem_empirical-alpha_2}.
\]

The proof of Theorem <ref> is in Appendix <ref>. The theorem demonstrates  addressing (<ref>) necessitates finding a proper $\eta$. We perform alternating optimization, where one side searches for an optimal $\eta$, and the other focuses on optimizing $\vec{x}$, $\bm{\mu}$, and $\bm{\nu}$, with $\eta$ held constant.

Algorithm <ref> is designed following the framework of bcd [Peng & Vidal, 2023, Huang et al., 2021, Gutman & Ho-Nguyen, 2023], shown as follows.

tributional erfactual ()
[1]
$\{\vec{x}'_i\}_{i=1}^n$, $\{y^*_j\}_{j=1}^n$ black-box model $b$, projections $\Theta$, bounds $U_{\mu},U_{\nu}$ and significance level $\alpha$.
Counterfactual distribution $\{\vec{x}_i\}_{i=1}^n$ or $\varnothing$.
$\vec{x}^{0}\leftarrow \vec{x}' + \sigma$; $t\leftarrow 0$
$\bm{\mu}^{t}\leftarrow \argmin_{\bm{\mu}:\bm{\mu}^{(\theta)}\in\Pi}Q_{\mu}(\bm{\mu}|\Theta, \vec{x}^{t},\vec{x}')$
$\bm{\nu}^{t}\leftarrow \argmin_{\bm{\nu}\in\Pi}Q_{\nu}(\bm{\nu} | \vec{x}^{t},\vec{y}^*)$

$\overline{Q}_{\nu}\leftarrow$ Eq. (<ref>)
$\overline{Q}_{\mu}\leftarrow$ Eq. (<ref>)
$\eta^{t}\leftarrow$ Algorithm <ref> or <ref> using $U_{\mu}-\overline{Q}_{\mu}$, $U_{\nu}-\overline{Q}_{\nu}$
$\renabla Q\leftarrow\renabla_{\vec{x}}Q\left(\vec{x}|\eta^{t}, \bm{\mu}^{t}, \bm{\nu}^{t}\right)$

$\vec{x}^{t+1}\leftarrow \retr({-\tau \renabla Q})$
$t\leftarrow t + 1$
$\norm{\vec{x}^{t+1}-\vec{x}^{t}}\leq \varepsilon_{\discount}$
$\overline{Q}_{\mu}\leq U_{\mu}$ and $\overline{Q}_{\nu}\leq U_{\nu}$
return $\{\vec{x}^{t+1}_i \}_{i=1}^n$
return $\varnothing$

Lines <ref> and <ref> compute the ot distances defined in (<ref>) and (<ref>) and obtain the optimal transportation plans $\bm{\mu}$ and $\bm{\nu}$, respectively. Here we use the symbol “$|$” to seperate between variables (e.g. $\bm{\mu}$ and $\bm{\nu}$) that are subject to optimization in the current line and parameters that stay fixed during this optimization (e.g. $\Theta$, $\vec{x}'$ etc.). In lines <ref> and <ref>, we compute the two corresponding ucl $\overline{Q}_{\nu}$ and $\overline{Q}_{\mu}$ that are defined in (<ref>) and (<ref>). Their values are determined by the quantile $u$ of current distributions $\{b(\vec{x}_i^{t})\}_{i=1}^n$ and  $\{\vec{x}_i^{t}\}_{i=1}^n$.

Because $\vec{x}^t=[\vec{x}^t_1\ldots,\vec{x}^t_n]$ is updated in every iteration, both ucl are re-computed per iteration. Then we use the two ucl to compute an $\eta$ as shown in line <ref>. The computed $\eta$ influences the optimization direction of the counterfactual distribution $\{\vec{x}_i\}_{i=1}^n$ (shown by line <ref>) to balance the satisfaction of the two chance constraints (<ref>) and (<ref>), ensuring that both are adequately addressed in the solution. Line <ref> computes the Riemannian gradient with respect to $\vec{x}$, and line <ref> makes a descent step using the retraction mapping [Absil et al., 2008] (denoted by $\retr$) from the Riemannian gradient.

Algorithm <ref> provides a way to select $\eta$ from a pre-defined set $\K$ consisting of multiple candidate values of $\eta$. If one value is expected to be selected more than once, then it is duplicated with corresponding copies in $\K$.

Set Shrinking, $\K\!=\!\{\eta_{(k)}\}_{k=1}^K$ ($0\leq\eta_{(k)}\!\!\leq 1$)
[1]
$\overline{Q}_{\mu}$, $\overline{Q}_{\nu}$, $U_{\mu}$,$U_{\nu}$, $\K$
$\eta \leftarrow$  Balance the gaps $U_{\mu}\!-\!\overline{Q}_{\mu}$ and $U_{\nu}\!-\!\overline{Q}_{\nu}$, $\eta\in \K$
$\K \leftarrow \K \backslash \{\eta\}$
Save $\K$ as the input of the next run
return $\eta$

Line <ref> optimizes $\eta$ towards balancing the satisfactory of ucl (<ref>) and (<ref>) (hence (<ref>) and (<ref>)). See Appendix <ref> for more details. The set $\K$ shrinks in every run as the selected $\eta$ gets removed from it, until $K=1$ and the last $\eta$  is used until convergence.

Algorithm <ref> optimizes $\eta$ within a pre-defined interval $[l,r]$, where the optimization criteria stays the same with Algorithm <ref>, shown by line <ref> (See Appendix <ref>).

Interval Narrowing, $\eta \in [l,r]$ ($0\leq l<r\leq 1$)
[1]
$\overline{Q}_{\mu}$, $\overline{Q}_{\nu}$, $U_{\mu}$,$U_{\nu}$, $[l,r]$, and $\kappa$ ($0<\kappa<1$)
$\eta \leftarrow$  Balance the gaps $U_{\mu}\!-\!\overline{Q}_{\mu}$ and $U_{\nu}\!-\!\overline{Q}_{\nu}$, $\eta\in[l,r]$
$\eta > (l+r)/2$
$l \leftarrow l + \kappa (r-l) $
$r \leftarrow r - \kappa (r-l) $
Save $[l,r]$ and $\kappa$ as the input for the next run
return $\eta$

In each run, the interval gets narrowed by a small proportion of $\kappa$, following a bi-section strategy, shown by lines <ref>–<ref>.

§ CONVERGENCE RATE ANALYSIS

To ease the presentation, we let $\bm{\rho} = [\vec{x}, \bm{\mu}, \bm{\nu}]$, and perform analysis for the optimization of $\bm{\rho}$ and $\eta$. The optimization problem in (<ref>) reads $\min_{\bm{\rho}} Q(\bm{\rho})$. Let $\bm{\rho}\in\M$ where $\M$ is a Riemannian manifold. Let $\vec{x}\in\M_1, \bm{\mu}\in\M_2, \bm{\nu}\in\M_3$ ($\M=\M_1\cup\M_2\cup\M_3$), where $\M_i$ ($i=1,2,3$) are compact smooth submanifolds and $\M_1$ is bounded. Assume $b$ is continuously differentiable and Lipschitz smooth with respect to $\vec{x}$. Let $L$ be the Lipschitz smoothness constant of $Q$ with respect to $\vec{x}$. It is shown by <cit.> that there exists a constant $\bL$ (dependent on $L$, $M_1$, and ) such that $Q(\vec{\retr(x), \bm{\mu}, \bm{\nu}}|\eta)$ is bounded by the Riemannian gradient of $\renabla_{\vec{x}}Q$ and $\bL$. Define

\begin{align*}
C & \triangleq  \sqrt{2\bL} + \gamma_1 L\cdot\sqrt{\frac{2}{\bL}}, \\
B & \triangleq \sup_{\bm{\rho}\in\M}\{Q_{\mu}(\bm{\rho}),Q_{\nu}(\bm{\rho})\}.
\end{align*}

where $B$ is ascertained finite (Appendix <ref>).
Theorem <ref> and Theorem <ref> below show the convergence of .

The iterates $\{\bm{\rho}^t\}_t$ of Algorithm <ref> (combined with Algorithm <ref>) with stepsize $\tau = 1/\bL$ satisfy

\begin{align*}
& \min_{t=0,1,\ldots T} \norm{\tnabla Q(\bm{\rho}^t)}  \leq  C\sqrt{\frac{1}{T+1}\bigg[Q_{\mu}(\bm{\rho}^0)+Q_{\nu}(\bm{\rho}^0)\bigg]}
\end{align*}

The iterates $\{\bm{\rho}^t\}_t$ of Algorithm <ref> (combined with Algorithm <ref>) with stepsize $\tau = 1/\bL$ satisfy

\begin{align*}
& \min_{t=0,1,\ldots T} \norm{\tnabla Q(\bm{\rho}^t)} \\
& \leq  C\sqrt{\frac{1}{T+1}\bigg[ Q_{\mu}(\bm{\rho}^{0})  + Q_{\nu}(\bm{\rho}^{0})+ \frac{r-l}{\kappa}B \bigg] }
\end{align*}

The proofs of the two theorems are provided in Appendix <ref> and <ref> respectively, based on the convergence proof for bcd on smooth manifolds in <cit.>. Briefly, both theorems state that the magnitude of the Riemannian gradient goes towards zero with the number of iterations. For Set Shrinking, the gradient is bounded by the values of $Q_{\mu}$ and $Q_{\nu}$ at the starting point. For Interval Narrowing, the upper bound of $Q_{\mu}$ and $Q_{\nu}$ matters on top of the size and narrowing rate of the interval.

§ NUMERICAL RESULTS

[German-Credit, MLPNet] The first plot “Risk” shows the model's output distribution. The second plot “All Features” shows the quantiles of the 1D projected (by $\Theta$) factual and counterfactual distributions.  The other plots show marginal distributions for each feature, where numerical ones are shown by quantile and categorical by histogram. Factual risk (average) is 31.3% and counterfactual 17.5%.

[German-Credit, MLPNet] Risk with respect to Credit amount, age, and duration (indicated by the size of each point).

[Cardiovascular Disease]  The models are trained on all features whereas the dce optimization is performed only on age, weight, and height.

[Cardiovascular Disease] The setup follows  <ref>. The difference is computed as counterfactual minus factual. The BMI is computed as “kg/m$^2$”.

We demonstrate how dce can be used to account for a model's decision making. The experiments are performed on two datasets. Three models are considered: three layers mlpnet and rbfnet, and svm. Throughout our experiments, $\alpha$ is set to $0.1$. This setting implies that the significance level for testing the distinctness of two distributions is $5\%$ as $\alpha/2$ corresponds to the two-tailed significance threshold. The code is accessible for experiment reproducibility[The code is available on
<https://anonymous.4open.science/r/distributional-counterfactual-explanation-F00F> ].

German-Credit. This dataset contains information about individuals who took a loan from a particular bank, describing people by a set of attributes as being with good or bad credit risks [Hofmann, 1994].

Cardiovascular Disease. This dataset consists of a significant number of patient records with a variety of features that are considered relevant for predicting cardiovascular diseases [Halder, 2020].

[Cardiovascular Disease, RBFNet] The model is trained on all features, whereas the dce optimization is performed only on the two features Diastolic/Systolic blood pressure. The four counterfactuals are obtained by optimizing from the same factual (show left). The darker color indicates the higher risk (see the color bar in  <ref>).

[Cardiovascular Disease, MLPNet] Convergence of  with Interval Narrowing. At the beginning  $\eta$ stays at $r=1$ due to the violation of the chance constraint of $y$, such that the optimization leans entirely towards bringing $y$ to $y^{*}$ closer gradually. When $\overline{Q}_{\nu}$ is below $U_2$ at iteration $57$, a feasible solution is found. Then $\eta$ is optimized within $[l,r]$ to balance the gaps $U_{\mu}-\overline{Q}_{\mu}$ and $U_{\nu}-\overline{Q}_{\nu}$, until the algorithm converges.

Validity and Proximity. Validity measures if a generated counterfactual indeed corresponds to a different outcome from the factual. Proximity refers to how close the generated counterfactual are to the original factual in terms of features. It is important to note that the choice of $U_{\mu}$ and $U_{\nu}$
depends on the specific features being explained. For instance, when conducting research on two distinct groups of individuals characterized by differing risk levels, one could calculate the Wasserstein distance between these groups for calibration purposes. This approach allows for a tailored assessment of the distributional differences in risk factors pertinent to the groups under study.
We train an mlpnet on the German-Credit dataset with all features. We sample 100 data points from the testing dataset as the factual distribution, then apply the proposed  on all training features to obtain the corresponding counterfactual distribution for explanation. Our target $y^*$ is generated by a Beta distribution ($0.1,0.9$) that represent a group of customer with very low risk.  Both the features and the model output (risk) of the factual and the counterfactual are shown in  <ref>. The counterfactual features stay close to its factual counterpart, yet resulting in a significant lower risk. Example of counterfactual points are shown in Appendix <ref>.

Counterfactual Reasoning. Our investigation into the German Credit dataset, depicted in  <ref> and  <ref>, explores the model's distinct responses to groups with similar feature distributions. Key observations include the model's sensitivity to age and credit amount, where minor variations can result in significant risk alterations, particularly in the 40–50 age group with around 8000 in credit amount. Notably, we identified a counterintuitive behavior where increasing credit amounts paradoxically lowers risk for younger individuals. This suggests a misinterpretation of risk factors by the model, mistaking age-related risk indicators as credit amount related. This is because there are quite many young (age 20–30) people who took a small credit amount for loan in the dataset, and their risk is marked as high because of their age (rather than the fact of taking small credit). Yet the model learns its relation to credit amount, failing to realize this being not causation.
Employing dce, we analyze quantile shifts across features, allowing us to delve into the model's decision-making process and discern substantial risk profile discrepancies between groups, even with identical distributions in certain features. This approach yields deeper statistical insights into the model’s distributional behavior and its implications.

In our study with the Cardiovascular Disease dataset, we employed three models: mlpnet, svm, and rbfnet, training them on all available features. We then selected a sample of 100 data points from the test set to represent our empirical factual distribution. Specifically, we curated a subset of 100 individuals with low cardiovascular risk to form the target output distribution $y^*$. The quantile-based comparison of the factual $y$ and target $y^*$, as shown in  <ref>, underscores the effectiveness of dce. Notably, the risk profiles of these counterfactual distributions, when converted to binary form, closely align with the target distribution.  <ref> sheds light on how the models respond to variations in age and bmi. The scatter plot and regression curve together reveal the nuanced relationships these models learn among risk, age, weight, and height.

We remark such insights for both dataset above are uniquely accessible through a dce framework, which goes beyond the capabilities of classical or group-based ce by focusing on distributional shifts and relationships. Additionally, we discussed how the transportation plan $\bm{\mu}$ and $\bm{\nu}$ can be used for reasoning in Appendix <ref>. We also showcase counterfactual samples as data point instances in Appendix <ref>.

Diversity.  <ref> shows the capability of  for generating diverse counterfactual distributions with different setups of $U_{\mu}$ given fixed $U_{\nu}$. These  counterfactuals differ in their features distributions, yet all yield very close prediction distributions. A small $U_{\mu}$ constraints the flexibility of counterfactual optimization but ascertains the generated distribution is close to observations. One could use $U_{\mu}$ to control the trade-off between validity and proximity.

Convergence.  <ref> demonstrates the convergence behavior of , where $\eta$ is computed by Interval Narrowing in Algorithm <ref>. We argue that the convergence behavior of Set Shrinking in Algorithm <ref> is similar (guaranteed by Theorem <ref>). We remark the set $\K$ is however more difficult to visualize than the interval $[l,r]$, hence its convergence behavior is not (repeatedly) shown here.

§ CONCLUSIONS

We have broadened the scope of ce by introducing it in a distributional context, termed as dce, under a statistically rigorous framework. Some explorations on geometric interpretations as well as unresolved challenges are discussed in Appendix <ref> for further reading.

[Absil et al., 2008]
Absil, P.-A., Mahony, R., and Sepulchre, R.
Optimization algorithms on matrix manifolds.
Princeton University Press, 2008.

[Amari et al., 2018]
Amari, S.-i., Karakida, R., and Oizumi, M.
Information geometry connecting wasserstein distance and kullback–leibler divergence via the entropy-relaxed transportation problem.
Information Geometry, 1:0 13–37, 2018.

[Bodria et al., 2023]
Bodria, F., Giannotti, F., Guidotti, R., Naretto, F., Pedreschi, D., and Rinzivillo, S.
Benchmarking and survey of explanation methods for black box models.
Data Mining and Knowledge Discovery, pp.  1–60, 2023.

[Brughmans et al., 2023]
Brughmans, D., Leyman, P., and Martens, D.
Nice: an algorithm for nearest instance counterfactual explanations.
Data Mining and Knowledge Discovery, pp.  1–39, 2023.

[Carrizosa et al., 2024]
Carrizosa, E., Ramírez-Ayerbe, J., and Morales, D. R.
Mathematical optimization modelling for group counterfactual explanations.
European Journal of Operational Research, 2024.

[Cédric, 2008]
Cédric, V.
Optimal transport: Old and new (grundlehren der mathematischen wissenschaften, 338).
English. Springer, pp.  998, 2008.

[Cheng et al., 2020]
Cheng, F., Ming, Y., and Qu, H.
Dece: Decision explorer with counterfactual explanations for machine learning models.
IEEE Transactions on Visualization and Computer Graphics, 270 (2):0 1438–1447, 2020.

[Cuturi, 2013]
Cuturi, M.
Sinkhorn distances: Lightspeed computation of optimal transport.
In Advances in Neural Information Processing Systems, volume 26. Curran Associates, Inc., 2013.

[Del Barrio & Loubes, 2019]
Del Barrio, E. and Loubes, J.-M.
Central limit theorems for empirical transportation cost in general dimension.

[del Barrio et al., 2023]
del Barrio, E., Sanz, A. G., Loubes, J.-M., and Niles-Weed, J.
An improved central limit theorem and fast convergence rates for entropic transportation costs.
SIAM Journal on Mathematics of Data Science, 50 (3):0 639–669, 2023.

[Dvoretzky et al., 1956]
Dvoretzky, A., Kiefer, J., and Wolfowitz, J.
Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator.
The Annals of Mathematical Statistics, pp.  642–669, 1956.

[Goldfeld et al., 2022]
Goldfeld, Z., Kato, K., Rioux, G., and Sadhu, R.
Statistical inference with regularized optimal transport.
arXiv preprint arXiv:2205.04283, 2022.

[Guidotti, 2022]
Guidotti, R.
Counterfactual explanations and how to find them: literature review and benchmarking.
Data Mining and Knowledge Discovery, pp.  1–55, 2022.

[Gutman & Ho-Nguyen, 2023]
Gutman, D. H. and Ho-Nguyen, N.
Coordinate descent without coordinates: Tangent subspace descent on riemannian manifolds.
Mathematics of Operations Research, 480 (1):0 127–159, 2023.

[Halder, 2020]
Halder, R. K.
Cardiovascular disease dataset, 2020.
URL <https://dx.doi.org/10.21227/7qm5-dz13>.

[Hofmann, 1994]
Hofmann, H.
Statlog (German Credit Data).
UCI Machine Learning Repository, 1994.
DOI: https://doi.org/10.24432/C5NC77.

[Huang et al., 2021]
Huang, M., Ma, S., and Lai, L.
A riemannian block coordinate descent method for computing the projection robust wasserstein distance.
In International Conference on Machine Learning, pp.  4446–4455. PMLR, 2021.

[Karakida & Amari, 2017]
Karakida, R. and Amari, S.-i.
Information geometry of wasserstein divergence.
In International Conference on Geometric Science of Information, pp.  119–126. Springer, 2017.

[Karimi et al., 2020]
Karimi, A.-H., Barthe, G., Balle, B., and Valera, I.
Model-agnostic counterfactual explanations for consequential decisions.
In International Conference on Artificial Intelligence and Statistics, pp.  895–905. PMLR, 2020.

[Khan & Zhang, 2022]
Khan, G. and Zhang, J.
When optimal transport meets information geometry.
Information Geometry, 50 (1):0 47–78, 2022.

[Manole et al., 2022]
Manole, T., Balakrishnan, S., and Wasserman, L.
Minimax confidence intervals for the sliced wasserstein distance.
Electronic Journal of Statistics, 160 (1):0 2252–2345, 2022.

[Mothilal et al., 2020]
Mothilal, R. K., Sharma, A., and Tan, C.
Explaining machine learning classifiers through diverse counterfactual explanations.
In Proceedings of the 2020 conference on fairness, accountability, and transparency, pp.  607–617, 2020.

[Nadjahi et al., 2020]
Nadjahi, K., Durmus, A., Chizat, L., Kolouri, S., Shahrampour, S., and Simsekli, U.
Statistical and topological properties of sliced probability divergences.
Advances in Neural Information Processing Systems, 33:0 20802–20812, 2020.

[Oizumi et al., 2014]
Oizumi, M., Albantakis, L., and Tononi, G.
From the phenomenology to the mechanisms of consciousness: integrated information theory 3.0.
PLoS computational biology, 100 (5):0 e1003588, 2014.

[Otto, 2001]
Otto, F.
The geometry of dissipative evolution equations: the porous medium equation.

[Pawelczyk et al., 2020]
Pawelczyk, M., Broelemann, K., and Kasneci, G.
Learning model-agnostic counterfactual explanations for tabular data.
In Proceedings of the web conference 2020, pp.  3126–3132, 2020.

[Peng & Vidal, 2023]
Peng, L. and Vidal, R.
Block coordinate descent on smooth manifolds.
arXiv preprint arXiv:2305.14744, 2023.

[Plumb et al., 2020]
Plumb, G., Terhorst, J., Sankararaman, S., and Talwalkar, A.
Explaining groups of points in low-dimensional representations.
In International Conference on Machine Learning, pp.  7762–7771. PMLR, 2020.

[Rawal & Lakkaraju, 2020]
Rawal, K. and Lakkaraju, H.
Beyond individualized recourse: Interpretable and interactive summaries of actionable recourses.
Advances in Neural Information Processing Systems, 33:0 12187–12198, 2020.

[Shorack & Wellner, 2009]
Shorack, G. R. and Wellner, J. A.
Empirical processes with applications to statistics.
SIAM, 2009.

[Ustun et al., 2019]
Ustun, B., Spangher, A., and Liu, Y.
Actionable recourse in linear classification.
In Proceedings of the conference on fairness, accountability, and transparency, pp.  10–19, 2019.

[Verma et al., 2020]
Verma, S., Boonsanong, V., Hoang, M., Hines, K. E., Dickerson, J. P., and Shah, C.
Counterfactual explanations and algorithmic recourses for machine learning: A review.
arXiv preprint arXiv:2010.10596, 2020.

[Wachter et al., 2017]
Wachter, S., Mittelstadt, B., and Russell, C.
Counterfactual explanations without opening the black box: Automated decisions and the gdpr.
Harv. JL & Tech., 31:0 841, 2017.

[Warren et al., 2023]
Warren, G., Keane, M. T., Gueret, C., and Delaney, E.
Explaining groups of instances counterfactually for xai: A use case, algorithm and user study for group-counterfactuals.
arXiv preprint arXiv:2303.09297, 2023.

[Yan et al., 2020]
Yan, C., Zhu, H., Korolko, N., and Woodard, D.
Dynamic pricing and matching in ride-hailing platforms.
Naval Research Logistics (NRL), 670 (8):0 705–724, 2020.

§ APPENDIX

§.§ Proof of Theorem <ref>

The proof of the theorem follows [Manole et al., 2022], which proves the confidence intervals for Wasserstein distance and its sliced version respectively for $\delta$-trimmed ($\delta\in (0,1/2)$) empirical distributions, without assumptions on the knowledge of the distributions. The trimming is necessary, such that certain amount of data points are left out for capturing the tail behavior of the distributions.

The ucl in (<ref>) comes from <cit.> with $r=2$, noticing that $\inf_{\bm{\nu}\in\Pi}Q_{\nu}(\vec{x}, \bm{\nu})$ defined in (<ref>) is the Wasserstein distance between the empirical distributions $\vec{x}$ and $\vec{x}'$. Consider $\downq_{\alpha, n, N}$ and $\upq_{\alpha, n, N}$ with $N=1$, Eq. (<ref>) reads

\[
\inf \PP_n\big(F_n^{-1}(\downq_{\alpha, n, 1}(u)) \leq  F^{-1}(u)\leq  F_n^{-1}(\upq_{\alpha, n, 1}(u)), \forall u\in(0,1) \big) \geq 1-\frac{\alpha}{2}.
\]

This indeed establishes a confidence band for the two sequences of functions, $\downq_{\alpha, n, 1}$ and $\upq_{\alpha, n, 1}$. The proof then pivots on the assertion that the 1d Wasserstein distance can be characterized as the $L^2$-norm of the quantile functions of $y$ ($y=b(\vec{x})$) and $\vec{x}'$ ($y'=b(\vec{x}')$). The two sequences of functions $\downq_{\alpha, n, 1}$ and $\upq_{\alpha, n, 1}$ that satisfy the above inequality then guarantee the following inequality:

\[
|F^{-1}_{y}(u) - F^{-1}_{y^*}(u)|\leq D_1(u)
\]

where $D_1(u)$ is defined in (<ref>). We recall that the left-hand side is the Wasserstein distance betweeen $y$ and $y^*$ (see (<ref>). Hence

\[
\inf_{\bm{\nu}\in\Pi}Q_{\nu}(\vec{x}, \bm{\nu}) \leq \frac{1}{1-2\delta} \int_{\delta}^{1-\delta} D_1(u)\d u
\]

where the correction factor $1/(1-2\delta)$ on the right-hand side makes it a trimmed estimator, such that a certain percentage (determined by $\delta$) of the tail values of the distributions can be disregarded in the computation. This inequality holds probabilistically with confidence level $\alpha$.

The ucl in (<ref>) follows <cit.>, which constructs $\downq_{\alpha, n, 1}$ and $\upq_{\alpha, n, 1}$ independently for each $\bm{\theta}$ ($\bm{\theta}\in\Theta$) projected 1d distributions $\bm{\theta}^\top\vec{x}$ and $\bm{\theta}^\top\vec{x}'$, based on the confidence band below

\[
\inf \PP_n\big(F_{\theta, n}^{-1}(\downq_{\alpha, n, N}(u)) \leq  F_{\theta}^{-1}(u)\leq  F_{\theta, n}^{-1}(\upq_{\alpha, n, N}(u)), \forall u\in(0,1) \big) \geq 1-\frac{\alpha}{2N}.
\]

where $F^{-1}_{\theta, n}$ is the (empirical) quantile function for the $\theta$-projected distribution and $N=|\Theta|$. Therefore,

\[
\sum_{i=1}^{n}\sum_{j=1}^{n}\abs{\bm{\theta}^\top\vec{x}_i - \bm{\theta}^{\top}\vec{x}'_j}^2\mu^{(\theta)}_{ij} \leq \frac{1}{1-2\delta} \int_{\delta}^{1-\delta} D_{\theta,N}(u)\d u
\]

where $D_N(u)$ is defined in (<ref>). Then taking summation over $\bm{\theta}\in\Theta$ on the left-hand side and  Lebesgue integral over $\theta\in\SS^{d-1}$ on the right-hand side yields

\[
\inf_{\bm{\mu}\in\Pi} Q_{\mu}(\vec{x}, \bm{\nu}) \leq \frac{1}{1-2\delta}\int_{\SS^{d-1}}\int_{\delta}^{1-\delta} D_{\theta, N}(u)\d u \d\omega_N(\theta)
\]

This inequality holds probabilisitically with confidence level $\alpha$.

Hence the conclusion.

§.§ Proof of Theorem <ref>

Remark that the minimization of $Q(\vec{x},\bm{\mu}, \bm{\nu} | \eta)$ is equivalent to the optimization below with respect to $\vec{x}$:

\[
\min_{\vec{x}} Q(\vec{x}|\eta) =  (1-\eta)\cdot \inf_{\bm{\mu}\in \Pi} Q_{\mu}(\vec{x},\bm{\mu}) + \eta\cdot \inf_{\bm{\nu}\in \Pi} Q_{\nu}(\vec{x}, \bm{\nu}).
\]

Since the problem (<ref>) is feasible, there exists at least one pair $(\vec{x}, \varphi) $ that satisfies the constraints (<ref>)-(<ref>) for any $ \eta $ in the interval [0, 1].
Given that $ Q_{\mu} $ and $ Q_{\nu} $ are continuous in their respective variables, the above definition ensures that $ Q(\vec{x}|\eta) $ is continuous in $\vec{x}$ for any fixed $ \eta $.

Let us consider the set $S$ of all $ \eta $ such that the corresponding $ \vec{x}_{\eta} = \arg\min_x Q(\vec{x}|\eta) $ is feasible for (<ref>). Since the constraints (<ref>)-(<ref>) are continuous in $ x $, and $ x_{\eta} $ is continuous in $\eta$, the set $S$ is nonempty and compact. Thus, there exists a minimum and maximum $ \eta $ in $ S $, denoted as $\eta_{\min}$ and $\eta_{\max}$ respectively.

We now show that at least for one $ \eta^* \in [\eta_{\min}, \eta_{\max}] $, the $ \vec{x}_{\eta^*} $ is optimal for (<ref>). Note that $\eta$ is a parameter that interpolates between $Q_{\mu}(\vec{x},\bm{\mu})$ and $  Q_{\nu}(\vec{x}, \bm{\nu})$. Since $ Q(\vec{x}|\eta) $ is continuous in $\vec{x}$ and $ \eta $, and the feasible region is closed and bounded, by the Extreme Value Theorem, $ Q(\vec{x}|\eta) $ achieves its minimum for some $ \vec{x}^{*} $ in the feasible region of (<ref>).

Hence, there exists an $ \eta^* \in [\eta_{\min}, \eta_{\max}] $ such that $ \vec{x}^* = \arg\min_x Q(\vec{x}|\eta^*) $ satisfies all the constraints (<ref>)-(<ref>) and thus is optimal for (<ref>). This concludes the proof.

§.§ Proof of Theorem <ref>

Consider the iteration process of  where $\eta$ sequentially gets to be $\eta_{(1)},\eta_{(2)},\ldots\eta_{(K)}$, with $\K=\{\eta_{(1)},\eta_{(2)},\ldots\eta_{(K)}\}$. Correspondingly, the function $Q$ defined in (<ref>) given each of these $\eta$ values are denoted by $Q_{(1)}, Q_{(2)}, \ldots Q_{(K)}$. Given that each value of $\eta$ may remain constant for several iterations before transitioning to the next, we denote $T_1, T_2, \ldots, T_K$ as the respective number of iterations for which it persists.

Consider $Q_{(1)}$ under $\eta_{(1)}$. It is shown in <cit.> that the sum of the square of the Riemannian gradient of $Q_{(1)}$ on $\bm{\rho}$ is bounded by $(b-1)\cdot C_{b-1}^2\cdot [Q_{(1)}(\bm{\rho}^{0})-Q_{(1)}(\bm{\rho}^{T_1})]$, where $b$ is the number of blocks of which the variables are subject to optimization. We have two blocks in , $\vec{x}$ and $[\bm{\mu}, \bm{\nu}]$, where the former is subject to Riemannian Gradient Descent (i.e. lines <ref> and <ref>) and the latter Exact Minimization (i.e. solving the ot problem to the exact minimum in lines <ref> and <ref>). Hence $b=2$, and we can use $C$ without index $b-1$ to represent the bound. This gives the first inequality below.

\begin{align}
\setcounter{ineq}{0}
\sum_{t=0}^{T_1-1} \norm{\renabla Q_{(1)}(\bm{\rho}^t)}^2 & \leq C^2\bigg[Q_{(1)}(\bm{\rho}^{0})-Q_{(1)}(\bm{\rho}^{T_1})\bigg] \nonumber \\
& \ineq{=} C^2 \bigg\{\bigg[ \bigg(1-\eta_{(1)}\bigg) Q_{\mu}(\bm{\rho}^0)+ \eta_{(1)} Q_{\nu}(\bm{\rho}^0) \bigg] - \bigg[\bigg(1-\eta_{(1)}\bigg)Q_{\mu}(\bm{\rho}^{T_1})+ \eta_{(1)} Q_{\nu}(\bm{\rho}^{T_1})\bigg]\bigg\} \nonumber \\
& = C^2 \bigg\{ \bigg(1-\eta_{(1)}\bigg)\bigg[ Q_{\mu}(\bm{\rho}^0) - Q_{\mu}(\bm{\rho}^{T_1}) \bigg] + \eta_{(1)} \bigg[ Q_{\nu}(\bm{\rho}^0) - Q_{\nu}(\bm{\rho}^{T_1}) \bigg] \bigg\} \nonumber  \\
& \ineq{\leq} C^2 \bigg\{ \bigg[ Q_{\mu}(\bm{\rho}^0) - Q_{\mu}(\bm{\rho}^{T_1}) \bigg] + \bigg[ Q_{\nu}(\bm{\rho}^0) - Q_{\nu}(\bm{\rho}^{T_1}) \bigg] \bigg\}
\label{eq:Q_(1)}
\end{align}

Note that (i) holds by the definition of $Q$ in (<ref>), and (ii) by $0\leq \eta \leq 1$, as well as the non-negativity of $Q_{\mu}(\bm{\rho}^0) - Q_{\mu}(\bm{\rho}^{T_1})$ and  $ Q_{\nu}(\bm{\rho}^0) - Q_{\nu}(\bm{\rho}^{T_1})$. This non-negativity is due to the fact that the Riemannian Gradient Descent yields descreasing objective.

Similarly, we can apply the same conclusions to $Q_{(2)},Q_{(3)}\ldots Q_{(K)}$ to obtain the following inequalities.

\begin{align}
\sum_{t=T_1}^{T_2-1} \norm{\renabla Q_{(2)}(\bm{\rho}^t)}^2 &\leq C^2\bigg[Q_{(2)}(\bm{\rho}^{T_1})-Q_{(2)}(\bm{\rho}^{T_2})\bigg] \nonumber \\
& \leq C^2 \bigg\{ \bigg[ Q_{\mu}(\bm{\rho}^{T_1}) - Q_{\mu}(\bm{\rho}^{T_2}) \bigg] + \bigg[ Q_{\nu}(\bm{\rho}^{T_1}) - Q_{\nu}(\bm{\rho}^{T_2}) \bigg] \bigg\} \nonumber\\
&  \quad\quad\quad \vdots \nonumber\\
\sum_{t=T_{K-1}}^{T_K-1} \norm{\renabla Q_{(K)}(\bm{\rho}^t)}^2 &\leq C^2\bigg[Q_{(K)}(\bm{\rho}^{T_{K-1}})-Q_{(K)}(\bm{\rho}^{T_K})\bigg]  \nonumber\\
& \leq C^2 \bigg\{ \bigg[ Q_{\mu}(\bm{\rho}^{T_{K-1}}) - Q_{\mu}(\bm{\rho}^{T_K}) \bigg]  + \bigg[ Q_{\nu}(\bm{\rho}^{T_{K-1}}) - Q_{\nu}(\bm{\rho}^{T_K}) \bigg] \bigg\}
\label{eq:Q_(K)}
\end{align}

Denote $T=T_K-1$. For any $t=0,1,\ldots T$, let $\varrho_t$ be the corresponding index of $\eta$ at the iteration $t$. Summing up the inequalities above, we obtain

\begin{align}
\setcounter{ineq}{0}
\sum_{t=0}^{T} \norm{\renabla Q_{(\varrho_t)}}^2
& \leq C^2 \bigg\{
\bigg[Q_{\mu}(\bm{\rho}^0) - Q_{\mu}(\bm{\rho}^{T_1}) +  Q_{\mu}(\bm{\rho}^{T_1}) - Q_{\mu}(\bm{\rho}^{T_2}) + \cdots +  Q_{\mu}(\bm{\rho}^{T_{K-1}}) - Q_{\mu}(\bm{\rho}^{T_K}) \bigg] \nonumber \\
&\quad \quad \quad  + \bigg[Q_{\nu}(\bm{\rho}^0) - Q_{\nu}(\bm{\rho}^{T_1}) +  Q_{\nu}(\bm{\rho}^{T_1}) - Q_{\nu}(\bm{\rho}^{T_2}) + \cdots +  Q_{\nu}(\bm{\rho}^{T_{K-1}}) - Q_{\nu}(\bm{\rho}^{T_K}) \bigg]
\bigg\} \nonumber \\
& = C^2 \bigg\{ \bigg[ Q_{\mu}(\bm{\rho}^0) - Q_{\mu}(\bm{\rho}^{T_K}) \bigg] + \bigg[ Q_{\nu}(\bm{\rho}^0) - Q_{\nu}(\bm{\rho}^{T_K}) \bigg]\bigg\} \nonumber \\
& \ineq{\leq} C^2 \bigg[ Q_{\mu}(\bm{\rho}^0) + Q_{\nu}(\bm{\rho}^0) \bigg]
\label{eq:sum_Q_discrete}
\end{align}

The inequality (i) holds due to the non-negativity of $Q_{\mu}$ and $Q_{\nu}$, as defined in (<ref>) and (<ref>) respectively.
Notice that the minimum of $\renabla Q_{\varrho_t}$ is bounded by this summation. Therefore,

\[
\min_{t=0,1,\ldots T} \norm{\renabla Q_{(\varrho_t)}}^2 \leq  \frac{1}{T+1} \sum_{t=1}^{T} \norm{\renabla Q_{(\varrho_t)}}^2 \leq \frac{C^2}{T+1}\bigg[ Q_{\mu}(\bm{\rho}^0) + Q_{\nu}(\bm{\rho}^0) \bigg]
\]

Taking root square of both side yields the theorem.

§.§ Proof of Theorem <ref>

We use the same notations in Appendix <ref>. Suppose $\eta$ sequentially gets to be $\eta_{(1)},\eta_{(2)},\ldots$, corresponding to functions $Q_{(1)}, Q_{(2)}, \ldots$
Let the initial interval in Interval Narrowing be $[l, r]$. At each iteration of the algorithm, the interval is reduced by a fixed proportion \( \kappa \) where \( 0 < \kappa < 1 \).
After \( t \) iterations, the length of the interval will be \( (1 - \kappa)^t (r - l) \).

Consider two arbitrary iteration steps corresponding to $\eta_{(k)}$ and $\eta_{(h)}$, where $h>k$ in the sequence of $\eta$. By the definition of $Q$ in (<ref>), we have

\[
Q_{(k)}(\bm{\rho}) =  (1-\eta_{(k)})\cdot Q_{\mu}(\bm{\rho}) +  \eta_{(k)} \cdot Q_{\nu}(\bm{\rho})
\]

\[
Q_{(h)}(\bm{\rho}) =  (1-\eta_{(h)})\cdot Q_{\mu}(\bm{\rho}) +  \eta_{(h)} \cdot Q_{\nu}(\bm{\rho})
\]

Next, we show that the gap between $Q_{(k)}(\bm{\rho})$ and $Q_{(h)}(\bm{\rho})$ is bounded by the interval length and $B$. Let $T_k$ be the number of elapsed iterations by the moment of $\eta_{(k)}$.

\begin{align}
\setcounter{ineq}{0}
\bigg|Q_{(k)}(\bm{\rho}) - Q_{(h)}(\bm{\rho})\bigg|
& = \bigg|  (1-\eta_{(k)})\cdot Q_{\mu}(\bm{\rho}) +  \eta_{(k)} \cdot Q_{\nu}(\bm{\rho}) - (1-\eta_{(h)})\cdot Q_{\mu}(\bm{\rho}) - \eta_{(h)} \cdot Q_{\nu}(\bm{\rho}) \bigg| \nonumber \\
& = \bigg|(\eta_{(h)}-\eta_{(k)}) Q_{\mu}(\bm{\rho}) + (\eta_{(k)}-\eta_{(h)}) Q_{\nu}(\bm{\rho})\bigg| \nonumber\\
& \ineq{\leq} |\eta_{(h)}-\eta_{(k)}|\cdot \bigg|Q_{\mu}(\bm{\rho}) - Q_{\nu}(\bm{\rho})\bigg| \nonumber\\
& \ineq{\leq} (1 - \kappa)^{T_k}(r - l) \bigg|Q_{\mu}(\bm{\rho}) - Q_{\nu}(\bm{\rho})\bigg| \nonumber\\
& \ineq{\leq} (1 - \kappa)^{T_k}(r - l)B
\label{eq:bound_by_B}
\end{align}

The step (i) comes from the triangle inequality. The step (ii) holds due to the interval narrowing mechanism. Namely, the interval narrows down to $(1-\kappa)^{T_k}(r-l)$ when we encounter $\eta_{(k)}$. Note that $\eta_{(h)}$ is encountered later than $\eta_{(k)}$, hence the interval at the moment of $\eta_{(h)}$ is a subset of that of $\eta_{(k)}$. Therefore, the gap between $\eta_{(k)}$ and $\eta_{(h)}$ is no larger than $(1-\kappa)^{T_k}(r-l)$. The step (iii) holds because of the definition of $B$ ($B\triangleq\sup_{\bm{\rho}\in\M}\{Q_{\mu}(\bm{\rho}),Q_{\nu}(\bm{\rho})\}$), as well as the fact that $Q_{\mu}$ and $Q_{\nu}$ are non-negative.

Consider the iteration process in . For any $t=0,1,\ldots T$, let $\varrho_t$ be the corresponding index of $\eta$ at the iteration $t$. Similarly to the proof in Appendix <ref>, we could also obtain the same inequalities as (<ref>)–(<ref>). The difference is that there are infinite number of such inequalities in Interval Narrowing, rather than the fixed number $K$ as for Set Shrinking. But one could consider an arbitrary $K$ in the $\eta$ sequence $\eta_{(1)},\eta_{(2)}\ldots\eta_{(K)}\ldots$ then let $T=T_K-1$ and sum up all these inequalities, which yields:

\begin{align*}
\setcounter{ineq}{0}
\sum_{t=0}^{T} \norm{\renabla Q_{(\varrho_t)}(\bm{\rho}^t)}
& \leq C^2 \bigg\{ \bigg[ Q_{(1)}(\bm{\rho}^{0}) -  Q_{(1)}(\bm{\rho}^{T_{1}})\bigg] + \bigg[ Q_{(2)}(\bm{\rho}^{T_{1}}) -  Q_{(2)}(\bm{\rho}^{T_{2}})\bigg]\bigg\} \\
& \quad\quad\quad + \cdots + \bigg[ Q_{(K)}(\bm{\rho}^{T_{K-1}}) -  Q_{(K)}(\bm{\rho}^{T_{K}})\bigg] \\
& \leq C^2\bigg\{ Q_{(1)}(\bm{\rho}^{0}) + \bigg[ Q_{(2)}(\bm{\rho}^{T_{1}}) -  Q_{(1)}(\bm{\rho}^{T_{1}})\bigg] + \bigg[ Q_{(3)}(\bm{\rho}^{T_{2}}) -  Q_{(2)}(\bm{\rho}^{T_{2}})\bigg] \\
& \quad\quad\quad + \cdots + \bigg[ Q_{(K)}(\bm{\rho}^{T_{K-1}}) -Q_{(K-1)}(\bm{\rho}^{T_{K-1}})  \bigg] -  Q_{(K)}(\bm{\rho}^{T_{K}})
\bigg\} \\
& \leq C^2\bigg\{ Q_{(1)}(\bm{\rho}^{0}) + \bigg|  Q_{(1)}(\bm{\rho}^{T_{1}}) - Q_{(2)}(\bm{\rho}^{T_{1}}) \bigg| + \bigg|   Q_{(2)}(\bm{\rho}^{T_{2}}) - Q_{(3)}(\bm{\rho}^{T_{2}}) \bigg| \\
& \quad\quad\quad + \cdots + \bigg| Q_{(K-1)}(\bm{\rho}^{T_{K-1}}) -  Q_{(K)}(\bm{\rho}^{T_{K-1}}) \bigg| -  Q_{(K)}(\bm{\rho}^{T_{K}})
\bigg\} \\
& \ineq{\leq} C^2\bigg\{ Q_{(1)}(\bm{\rho}^{0}) + (1 - \kappa)^{T_{1}}(r - l)B + (1 - \kappa)^{T_{2}}(r - l)B \\
& \quad\quad\quad + \cdots + (1 - \kappa)^{T_{K-1}}(r - l)B-  Q_{(K)}(\bm{\rho}^{T_{K}})
\bigg\} \\
& \ineq{\leq}  C^2\bigg\{Q_{(1)}(\bm{\rho}^{0}) + (r - l)B\bigg[ 1 + (1-\kappa) + (1-\kappa)^2 + \cdots  \bigg]  -  Q_{(K)}(\bm{\rho}^{T_{K}}) \bigg\} \\
& \ineq{\leq} C^2\bigg[ Q_{(1)}(\bm{\rho}^{0}) - Q_{(K)}(\bm{\rho}^{T_{K}}) + \frac{1}{\kappa}(r - l)B
\bigg] \\
& \ineq{\leq} C^2\bigg[  Q_{(1)}(\bm{\rho}^{0}) + \frac{1}{\kappa}(r - l)B
\bigg] \\
& \ineq{=} C^2 \bigg\{\bigg[ \big(1-\eta_{(1)}\big) Q_{\mu}(\bm{\rho}^{0})+ \eta_{(1)} Q_{\nu}(\bm{\rho}^{0}) \bigg]  + \frac{1}{\kappa}(r - l)B  \bigg\} \\
& \ineq{\leq} C^2 \bigg[ Q_{\mu}(\bm{\rho}^{0})  + Q_{\nu}(\bm{\rho}^{0})+ \frac{1}{\kappa}(r - l)B \bigg]
\end{align*}

The step (i) holds because of the inequality (<ref>) for arbitrary two values $\eta_{(k)}$ and $\eta_{(h)}$ ($h>k$), following which we obtain the summation of the geometric series in (ii) ($0<\kappa<1$). In step (iii) we take the summation of the series. Step (iv) is because the non-negativity of $Q_{(K)}$. Step (v) is by the definition of $Q$ in (<ref>). Finally, step (vi) is because of $0<\eta<1$ as well as the non-negativity of $Q_{\mu}$ and $Q_{\nu}$.

Notice that the minimum of $\renabla Q_{\varrho_t}$ is bounded by this summation $\sum_{t=0}^{T} \norm{\renabla Q_{(\varrho_t)}(\bm{\rho}^t)}$. Therefore,  we have

\[
\min_{t=0,1,\ldots T} \norm{\renabla Q_{(\varrho_t)}}^2\leq \frac{1}{T+1}\sum_{t=0}^{T}\norm{\renabla Q_{(\varrho_t)}(\bm{\rho}^t)} \leq \frac{C^2}{T+1}\bigg[ Q_{\mu}(\bm{\rho}^{0})  + Q_{\nu}(\bm{\rho}^{0})+ \frac{1}{\kappa}(r - l)B \bigg]
\]

We remark that $B$ is finite because of its definition $B\triangleq\sup_{\bm{\rho}\in\M}\{Q_{\mu}(\bm{\rho}),Q_{\nu}(\bm{\rho})\}$: Note that $Q_{\mu}$ and $Q_{\nu}$ are continuous on the bounded and compact manifold $\M_1$ with respect to $\vec{x}$. Since $\bm{\mu}, \bm{\nu}\in\Pi$ (namely, each being non-negative and with elements summed up to 1) during the optimization, the two functions $Q_{\mu}(\bm{\rho})$ and $Q_{\nu}(\bm{\rho})$ are therefore bounded.

Hence the conclusion.

§.§ Optimization of $\eta$

We explore the optimization strategy for $\eta$, as outlined in line 1 of both Algorithm <ref> and Algorithm <ref>. The parameter $\eta$, as defined in (<ref>), plays a crucial role in determining the descent direction during the optimization process, as seen in line <ref> of Algorithm <ref>. The optimization aims to adjust this direction to fulfill the dual chance constraints (<ref>) and (<ref>), effectively maximizing the probability of meeting both constraints.

We first consider the Interval Narrowing in Algorithm <ref>. The general strategy involves the following steps:

* When either $U_{\mu} - \overline{Q}_{\mu}$ or $U_{\nu} - \overline{Q}_{\nu}$ is negative (indicating constraint violation), while the other is non-negative, priority is given to addressing the constraint that is violated (i.e., the one with the negative value). In this case, optimizing the non-violated constraint is less urgent.

* If both $U_{\mu} - \overline{Q}_{\mu}$ and $U_{\nu} - \overline{Q}_{\nu}$ are negative or positive, $\eta$ is assigned a value based on these gaps. Let $a = U_{\mu} - \overline{Q}_{\mu}$ and $b = U_{\nu} - \overline{Q}_{\nu}$. The value of $\eta$ is then determined as follows:

\begin{equation}
\eta = \left\{\begin{array}{ll}
\frac{b}{a+b} & \text{if $a$ and $b$ are both negative} \\
\frac{a}{a+b} & \text{if $a$ and $b$ are both non-negative (but not both zero)} \\
0.5 & \text{if both $a$ and $b$ are exactly zero (unlikely to happen numerically in computation)}
\end{array}\right.
\label{eq:eta_strategy}
\end{equation}

In the case where both $a$ and $b$ are negative, the term with the larger absolute value (hence, more significantly violating its constraint) receives a greater weight in the optimization of (<ref>). Conversely, if both are non-negative (and not both zero), it indicates that both constraints are satisfied, but the one with the smaller absolute value is closer to its limit, and thus receives more attention in the optimization.

Then, we consider the Set Shrinking in Algorithm <ref>. The main different is that $\eta$ must be selected from $\K$. One could perform sorting on $\K$ first to obtain a $\K_{\text{sorted}}$, and then use bisection search to find the element in $\K_{\text{sorted}}$ that is close to the $\eta$ in (<ref>).

§.§ Geometric Explanations and Open Problems Discussions

The 2-Wasserstein distance effectively determines how “far apart” two distributions are in a geometric sense. The confidence interval defined by the Wasserstein ball then represents a geometric region on this manifold, encompassing counterfactual distributions that are within a certain 'distance' (as measured by the 2-Wasserstein metric) from the central factual distribution.

trig format plots=rad,
hide axis,

3 [surf, opacity=0.3, domain=0:6, domain y=0:6, samples=30,
samples y=50]

3 [
only marks,
ball color=red,
mark size=0.6mm,
nodes near coords=$x'$,
name path=node1
] coordinates (3, 3, 0.3);

3 [
only marks,
ball color=blue,
mark size=0.6mm,
nodes near coords=$x$,
name path=node2
] coordinates (4, 3.5, 0.1);

3 [
no marks,
] coordinates
(3, 3, 0.3)
(3.5, 3.25,  (1+sin(deg(3.5)/10))*sin(deg(3.5)/2)*cos(deg(3.25)) )
(3.75, 3.375,  (1+sin(deg(3.75)/10))*sin(deg(3.75)/2)*cos(deg(3.375)) )
(4, 3.5,  (1+sin(deg(4)/10))*sin(deg(4)/2)*cos(deg(3.5)) )

axis lines=middle,
zmin=-0.5, zmax=0.5,
ymin=-0.5, ymax=0.5,
xmin=-0.5, xmax=0.5,
axis equal,
axis line style=draw=none,
scale uniformly strategy=units only,

3 [
only marks,
ball color=red,
mark size=1.25,
nodes near coords=$x'$,
nodes near coords style=above,
] coordinates (0, 0, 0);

[dashed, thick, blue] (axis cs:0,0,0) – (axis cs:0.5,0,0) node[midway, below right] $U_{\mu}$;

3 [
y domain=0:180,
variable y=,̌
point meta=u,
] (
0.5 * sin(v) * cos(u),
0.5 * sin(v) * sin(u),
0.5 * cos(v)

Geometric illustration of geodesic distance between two distributions $x$ and $x'$ and the Wasserstein ball centered at $x'$.

In the context of a Riemannian manifold, geodesics are the shortest paths between two points (in this case, between two probability distributions). The 2-Wasserstein distance is a geodesic distance on the manifold <cit.>. This framework allows for a geometric interpretation of how distributions relate to each other within this space. For any given factual distribution $x'$, the region of counterfactual distribution $x$ defined by the probability constraint can geometrically be visualized as a kind of radius $U_{\mu}$ sphere (i.e. Wasserstein Ball) centered at $x'$ where a significant majority (more than the probability $\varphi$) strictly adheres to the boundary constraint of the sphere. See  <ref> for an illustration.

Transportation Plan. The transportation plan can be conceptualized as a geometrical mapping that delineates the optimal transfer of probability mass across the manifold from the factual input distribution $x'$ to the counterfactual input distribution $x$. Embedded within the manifold's structure, the values of $\bm{\mu}$ and $\bm{\nu}$ represent the relative extent of mass reallocation necessary between specific points or regions in the probability space. Large values in $\vec{\Pi}$ indicate areas requiring substantial mass movement, reflective of significant distributional changes, while smaller values suggest minor adjustments. Zero values denote regions where the existing distribution aligns with the target, requiring no modification. This transportation plan, therefore, becomes an integral tool in understanding and visualizing the nature of distributional shifts required in counterfactual scenarios.

[Cardiovascular Disease, MLPNet] Visualization of the transportation plan for $\bm{\mu}$ and $\bm{\nu}$, as computed in lines <ref> and <ref> of Algorithm <ref>, respectively. This heatmap depicts the averaged values of $\bm{\mu}$ over $\bm{\theta}\in\Theta$, normalized to form a distribution with elements summing to one. Factual and counterfactual points are aligned (in their original order) along the x-axis and y-axis, respectively.

<ref> demonstrates the utilization of the transportation plans $\bm{\mu}$ and $\bm{\nu}$ in dce. A key observation from the heatmap of $\bm{\mu}$ is the pronounced diagonal, suggesting that the predominant mass transfer occurs between each factual point and its direct counterpart in the counterfactual dataset. This implies that individuals in the factual group are most likely to correspond to their original counterparts in the counterfactual group, rather than to others. Such a pattern is typical when the distribution distance between $\vec{x}$ and $\vec{x}'$ is constrained to be very close. In the extreme case of identical distributions, ot would occur strictly along the diagonal (and hence the heatmap $\bm{\nu}$ also looks strictly diagonal).

Conversely, the heatmap of $\bm{\nu}$ displays a less distinct diagonal pattern, indicating a notable shift in the risk ranking of individuals from the factual to the counterfactual group. This suggests that the individual with the highest cardiovascular risk in the factual group might not maintain a comparably high risk in the counterfactual group. While this is not inherently problematic, it highlights that the model may be highly sensitive to certain features, with minor changes in $\vec{x}$ potentially leading to significant alterations in $\vec{y}$. However, if the diagonal is prominently visible in $\bm{\mu}$ but not in $\bm{\nu}$, further investigation is warranted to determine whether the model is behaving as intended. This discrepancy may suggest that while the model maintains individual correspondences between the factual and counterfactual groups (as indicated by $\bm{\mu}$), the implications or outputs (reflected in $\bm{\nu}$) do not preserve these correspondences. Such a scenario could reveal potential issues or sensitivities in the model's learning process that merit closer examination.

eot and Information Geometry There is also a link between ot and information geometry, established by introducing an entropy-regularization term for $\ot$ [Cuturi, 2013], i.e.

\begin{equation}
\ot_{\varepsilon}(\gamma_1, \gamma_2) =  \inf_{\pi\in\Pi(\gamma_1, \gamma_2)} \bigg\{ \int_{\RR^{d\times d}} \norm{a_1-a_2}^2 \d{\pi(a_1,a_2)}  + \varepsilon H(\pi|\gamma_1 \otimes \gamma_2)\bigg\}
\end{equation}

where the regularization term $H(\pi|\gamma_1 \otimes \gamma_2)$ denotes the relative entropy with respect to the product measure $\gamma_1\otimes \gamma_2$,  defined as

\begin{equation}
H(\pi|\gamma_1 \otimes \gamma_2) = \int_{\RR^d} \log\left(\frac{\d{\pi}}{\d{\gamma_1}\d{\gamma_2}}\right).
\end{equation}

Using $\ot_{\varepsilon}$ as replacement for its counterparts in either / both the objective and the constraint of the formulation (<ref>) reduces the computational effort to obtain the transportation plan $\pi$. Furthermore, this regularization term is equivalent to using kl divergence as a regularizer and naturally defines certain geometrical structures from the information geometry point of view [Karakida & Amari, 2017, Khan & Zhang, 2022]. The role of the product of marginal distributions $\gamma_1\otimes \gamma_2$ in the kl divergence term is to establish a baseline of independence between the distributions being compared. By measuring the divergence of the transport plan $\pi$ from this baseline, the entropy term encourages solutions that respect the underlying structure of the individual distributions, rather than overly concentrating the mass transfer in a few specific regions. By increasing the value of $\varepsilon$, one tends to seek for a counterfactual $x$ that is drawn "more independently" from the distribution of the orignial input $x'$, leading to solutions that consider a broader range of possibilities, resulting in more diverse and potentially more insightful counterfactuals that better capture the overall structure of the data.

The integration of the entropic term makes $\sot$ a quasi-distance (due to loss of symmetry). The entropic Wasserstein quasi-distance may be regarded as such a metric in scenarios where the expense associated with the transfer of information across distinct nodes is contingent upon the spatial coordinates of these nodes [Oizumi et al., 2014], and the entropic term controls the amount of information integration [Amari et al., 2018], which is measured by the amount of interactions of information among different nodes.

Open Problems of dce with eot. The solution of the chance-constrained optimization problem (<ref>) is essentially based on a statistically trustworthy estimation of the upperbounds of $\ot$ and $\sot$. Remark that Theorem <ref> hinges on the fact that the
1d Wasserstein distance may be expressed as quantile functions, so it does not generalize to entropic regularized $\ot$ or $\sot$ (because they are not rigorously tied to the quantile expressions anymore). A clt is proved for the eot cost that is centered at the population cost, and this yields asymptotically valid confidence interval for eot <cit.>. Additionally, bootstrap is shown to be valid for eot <cit.>, making asymptotic inference for the interval efficient and straightforward. To the best of our knowledge, whether the sliced version of eot accepts an (even asymptotic) inference for the confidence interval remains open. Hence if one uses eot in (<ref>), its theoretical foundation calls for further investigations.

§.§ Examples of Counterfactual Data

[German-Credit] Example data pairs for factual (odd row, white) and counterfactual (even row, grey). Model training and counterfactual optimization are on all features.

Age     Sex     Job     Housing     Saving     Checking     Credit     Duration     Purpose     Risk

39     male     2     own     little     Unknown     1473     18     radio/TV     0.0722

39     male     2     own     little     Unknown     1442     17     radio/TV     0.0241

54     male     2     own     little     little     3051     48     domestic appliances     0.9197

54     male     2     own     little     little     3029     48     domestic appliances     0.9178

27     male     2     own     quite rich     moderate     2326     15     business     0.0625

27     male     2     own     quite rich     moderate     2586     14     business     0.0368

40     male     1     own     little     moderate     1155     12     radio/TV     0.0864

40     male     1     own     little     moderate     1385     10     radio/TV     0.0477

36     male     2     own     Unknown     Unknown     3079     36     car     0.1539

36     male     2     own     Unknown     Unknown     3169     36     car     0.0218

35     male     1     own     moderate     little     4380     18     car     0.0635

36     male     1     own     moderate     little     4472     17     car     0.0239

25     female     1     own     little     moderate     1355     24     car     0.8184

25     female     1     own     little     moderate     1332     24     car     0.8179

24     female     2     rent     little     little     652     12     furniture/equipment     0.7827

24     female     2     rent     little     little     662     12     furniture/equipment     0.7811

24     male     1     own     Unknown     little     3763     21     car     0.3846

25     male     1     own     Unknown     little     4104     17     car     0.1197

32     male     2     rent     moderate     moderate     11760     39     education     0.9343

32     male     2     rent     moderate     moderate     11766     39     education     0.9354

43     male     2     rent     little     little     2625     16     car     0.5006

45     male     2     rent     moderate     little     2489     10     car     0.0259

74     male     3     own     little     Unknown     4526     24     business     0.2316

74     male     3     own     little     Unknown     3975     22     business     0.0235

61     male     1     rent     moderate     moderate     2767     21     business     0.1984

62     male     1     rent     moderate     moderate     3127     18     business     0.0219

30     male     2     own     Unknown     moderate     2028     12     car     0.0877

31     male     2     own     Unknown     moderate     1901     11     car     0.0766

39     female     1     own     Unknown     moderate     932     6     education     0.0949

42     female     1     own     Unknown     moderate     1578     4     education     0.0312

22     male     2     own     little     moderate     2039     24     radio/TV     0.4255

23     male     2     own     little     moderate     3092     22     radio/TV     0.1678

[Cardiovascular Disease] Example data pairs for factual (odd row, white) and counterfactual (even row, grey). Model training and counterfactual optimization are on all features.

Age     Gender     Height     Weight     Systolic      Diastolic     Chol     Gluc     Smoke     Alco     Active     Risk

17427     2     159     58.0     120     80     1     1     0     0     1     0.2225

17438     2     159     57.7     114     75     1     1     0     0     1     0.1399

20406     1     162     89.0     150     70     1     3     0     0     1     0.7296

20364     1     162     88.7     166     80     1     3     0     0     1     0.9224

17733     2     166     70.0     150     90     1     1     0     0     0     0.9133

17754     2     166     69.9     151     87     1     1     0     0     0     0.9222

18065     2     165     67.0     120     80     1     1     0     0     1     0.2618

18059     2     165     67.0     107     72     1     1     0     0     1     0.1403

18169     2     167     70.0     100     70     1     1     0     0     1     0.1548

18151     2     167     70.2     101     72     1     1     0     0     1     0.1419

19776     1     160     55.0     120     80     3     3     0     0     1     0.5167

19709     1     160     54.7     102     64     3     3     0     0     1     0.3027

18371     2     182     121.0     130     90     3     1     1     1     1     0.8096

18383     2     182     120.8     138     97     3     1     1     1     1     0.9222

20224     2     174     102.0     120     80     2     2     1     0     1     0.4643

20193     2     174     102.0     98     64     2     2     1     0     1     0.2004

16789     1     158     67.0     170     100     1     1     0     0     1     0.9838

16825     1     158     67.1     171     97     1     1     0     0     1     0.9842

21894     2     155     55.0     130     90     1     1     0     0     1     0.6086

21898     2     155     55.1     149     100     1     1     0     0     1     0.9226

14496     1     150     68.0     120     80     2     1     0     0     1     0.3394

14490     1     150     67.9     109     71     2     1     0     0     1     0.1403

23264     1     170     53.0     140     80     3     3     0     0     1     0.8115

23280     1     170     52.9     149     80     3     3     0     0     1     0.8844

16676     1     163     65.0     120     80     1     1     0     0     1     0.2272

16659     1     163     64.9     114     78     1     1     0     0     1     0.1402

20476     1     167     59.0     130     80     1     1     0     0     0     0.5574

20397     1     167     58.8     106     67     1     1     0     0     0     0.2061

15511     2     176     68.0     120     80     1     1     0     0     0     0.2036

15513     2     176     67.8     113     78     1     1     0     0     0     0.1400

21454     1     150     53.0     140     80     1     1     0     0     1     0.7558

21485     1     150     53.1     151     88     1     1     0     0     1     0.9224

18812     1     160     63.0     125     80     1     1     0     0     1     0.3501

18784     1     160     63.1     109     72     1     1     0     0     1     0.1491

14626     1     166     67.0     120     80     1     1     0     0     1     0.1644

14630     1     166     67.1     119     81     1     1     0     0     1     0.1403