# Collective Counterfactual Explanations via Optimal Transport

Ahmad-Reza Ehyaei
Max Planck Institute for Intelligent Systems, Tübingen AI Center, Tübingen,
Germany
<EMAIL_ADDRESS>
&Ali Shirali
University of California, Berkeley
<EMAIL_ADDRESS>&Samira Samadi
Max Planck Institute for Intelligent Systems, Tübingen AI Center, Tübingen,
Germany
<EMAIL_ADDRESS>
###### Abstract

Counterfactual explanations provide individuals with cost-optimal actions that
can alter their labels to desired classes. However, if substantial instances
seek state modification, such individual-centric methods can lead to new
competitions and unanticipated costs. Furthermore, these recommendations,
disregarding the underlying data distribution, may suggest actions that users
perceive as outliers. To address these issues, our work proposes a collective
approach for formulating counterfactual explanations, with an emphasis on
utilizing the current density of the individuals to inform the recommended
actions. Our problem naturally casts as an optimal transport problem.
Leveraging the extensive literature on optimal transport, we illustrate how
this collective method improves upon the desiderata of classical
counterfactual explanations. We support our proposal with numerical
simulations, illustrating the effectiveness of the proposed approach and its
relation to classic methods.

_K_ eywords Counterfactual Explanations $\cdot$ Optimal Transport $\cdot$
Collective Counterfactual Explanations

##  1 Introduction

For individuals adversely affected by algorithmic decisions, such as a
rejection from a loan application, counterfactual explanations (CE) Wachter et
al. (2017) provide the least costly action (e.g., increasing savings or salary
to a specific threshold) to attain a favorable outcome, such as securing loan
approval. To provide a mathematical formulation, consider a population
characterized by a feature space $\mathcal{X}\subseteq\mathbb{R}^{d}$, with
the decision process modeled by a classifier
$h:\mathcal{X}\rightarrow\mathcal{Y}$, where $h(x)$ maps instances to a
discrete set, such as $\mathcal{Y}=\\{\pm 1\\}$. In this context, an
individual $x$ dissatisfied with her label ($h(x)=-1$) is advised by
$\textbf{CE}(x)$ to take a minimum-cost action to achieve the desired label:

$\textbf{CE}(x)=\underset{x^{\prime}\in\mathcal{A}}{\text{arg\,min}}\,\\{c(x,x^{\prime}):\quad
h(x^{\prime})=+1\\},$ (1)

here, $c:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ is the
modification cost function, with $c(x,x^{\prime})$ representing the cost of
modifying $x$ to $x^{\prime}$, and $\mathcal{A}\subset\mathcal{X}$ is the
subset of feasible actions Ustun et al. (2019). Let’s denote the decision
boundary of $h$ as $L$. In a finely-grained feature space with a continuous
cost function, we expect that $\textbf{CE}(x)$ will lie on $L$.

Figure 1: CE often transitions a data point from a negative decision space to
a position near the decision boundary $L$. In classifiers with a margin, this
can result in the point being located in a low-probability region, an outlier
concerning the underlying distribution. Such mappings may also accumulate data
points in certain areas, creating new competition.

As depicted in Fig. 1, while the $\textbf{CE}(x)$ achieves cost optimization,
it fails to consider the underlying distribution of the feature space, denoted
by $P$. As a result, the proposed explanations, concentrated in low-density
areas proximal to the decision boundary, may represent an outlier concerning
$P$. This oversight is particularly significant in classifiers such as SVM,
which emphasizes the margin between the decision boundary $L$ and the nearest
data point.

Neglecting the underlying distribution also questions the robustness of the
solution to misspecified costs or feasibility constraints. In fact, social
structure or unobserved costs may have already pushed away individuals from
low-density areas. Recommending to such areas might no longer be beneficial to
the users. To address the above issues, it is often recommended to include
data manifold closeness among the CE’s desiderata Karimi et al. (2022);
Guidotti (2022); Verma et al. (2020).

Another issue with the current approach to CE is its individual-centric
formulation. Since CE chooses the actions for each individual independently,
if a substantial portion of the population seeks to change their states, an
accumulation of modified states may be observed near the boundary. This
externality can lead to increased competition, thereby escalating the cost of
modification. Such phenomenon exemplifies the tragedy of the commons Gross and
De Dreu (2019), similar to a navigation algorithm directing users to
individually optimal routes, which could be avoided or minimized with a
collective approach.

###  1.1 Our Contribution

Towards addressing the aforementioned issues, we propose collective
counterfactual explanations (CCE). In this approach, we guide individuals to
the denser areas of the feature space while ensuring that following these
guidelines does not introduce new competition. This involves a collective
approach to recommendation and cost minimization. We choose target regions
using the underlying distribution, classifier level curves, and the cost
function. To prevent CE from concentrating in a specific area during
optimization, we add a constraint to spread CE across the target region.
Optimal transport (OT) theory, which efficiently moves a source distribution
to a target, is a natural fit for this approach and is central to our study.
Our main contributions are:

  * •

Introduce CCE for more robust CE and externality minimization, aligning with
data manifold closeness ($\S$ 3.1).

  * •

Propose a new interpretation for the cost of modification in an equilibrium
market ($\S$ 3.1).

  * •

Demonstrate the improvement of standard CE through CCE across various
desiderata. ($\S$3.2).

  * •

Extend CE with a temporal dimension through path-guided counterfactuals
($\S$4).

  * •

Present a fast computation approach for estimating CCE with amortized
inference properties ($\S$5).

#### Organization.

In $\S$ 3, we introduce CCE under the OT framework and demonstrate the
existence of CCE under minimal conditions on the cost function and the
classifier. Leveraging dual OT, we interpret CE as a competition equilibrium,
providing a valuable perspective for estimating the cost function in
practical, real-world applications. In $\S$ 4, we propose that beyond the
state of the final recommendation in CE, the trajectory from $x$ to a
$\textbf{CE}(x)$ should also be optimal. We illustrate how a dynamic extension
of OT can provide path-guided CE. Additionally, we demonstrate that OT
provides a solution for ordered classifiers and utility-based explanations,
two natural extensions of the CE problem. Finally, in $\S$ 5, we explore the
Sinkhorn and back-and-forth methods for solving OT problems. Through
simulations, we demonstrate the effectiveness of CCE and compare it with the
standard CE solution. The appendix provides more detailed information about
the simulations. To have a self-contained document, all definitions are
included in the appendix and are linked to its first appearance.

Figure 2: The OT map offers a collective approach for CE, moving the negative
part of the distribution to the $\delta$-confidence region alongside decision
boundary.

###  1.2 Further Related Works

Poyiadzi et al. (2020); Kanamori et al. (2020) encode the population
probability into the cost function, effectively resulting in finding CE in
denser areas near the boundary $L$. While these approaches are beneficial in
mitigating outlier recommendations, they remain individual-centric and are
prone to overcrowding. Additionally, these approaches lack a principled
foundation, as the integration of density into the cost function appears
somewhat arbitrary.

Recently, Carrizosa et al. (2024) introduced a notion of collective CE,
focusing on optimizing the modification cost for a group of instances rather
than for individuals. This approach aims to harmonize the behavior of CE
within a group, thereby mitigating the occurrence of cost-outlier CE. However,
by overlooking the underlying density, this method might still result in CE
that are outliers with respect to probability measures. Additionally, coupling
CEs within a group could amplify the externalities.

##  2 Optimal Transport

This section provides an overview of OT and its extensions utilized in this
study. OT initially introduced by Monge (1781), focused on cost-efficient
transportation of soil for fortifications, but in general, OT aims to transfer
probability measures from space $\mathcal{X}$ to $\mathcal{X}^{\prime}$. While
$\mathcal{X}$ and $\mathcal{X}^{\prime}$ are generally Polish spaces, in the
scope of this work, they are considered open or closed bounded subsets of
$\mathbb{R}^{d}$.

#### Monge Problem.

Let $P\in\mathcal{P}(\mathcal{X}),Q\in\mathcal{P}(\mathcal{X}^{\prime})$ be
two probability measures, and let
$c:\mathcal{X}\times\mathcal{X}^{\prime}\rightarrow\mathbb{R}$ be a ground-
cost function representing the cost of transporting a unit mass from $x$ to
$y$. The _Monge’s OT problem_ is to find a map
$T:\mathcal{X}\rightarrow\mathcal{X}^{\prime}$ that push-forwards $P$ to $Q$
(i.e., $T_{*}P=Q$) and minimizes the effort of transportation:

$\underset{T:T_{*}P=Q}{\text{min}}\,\quad\underset{\mathbf{X}\sim
P}{\mathbb{E}}[c(\mathbf{X},T(\mathbf{X}))].$ (2)

An optimal map is defined as $T$ that minimizes the given equation. The
function $T$ satisfying the property $T_{*}P=Q$ is called push-forward map
(Ambrosio et al., 2021, $\S$ 1.2).

Regardless of the cost function $c$, Monge’s problem may be ill-posed due to
the nonexistence of a push-forward map and the lack of weak sequential closure
issues w.r.t. the weak topology Ambrosio et al. (2013). After 150 years,
Kantorovich (1942) addressed these weaknesses by relaxing the problem.

#### Kantorovich Problem.

Instead of finding an optimal map, Kantorovich proposed to minimize the
transportation cost for coupling $\pi\in\Gamma(P,Q)$:

$\underset{\pi\in\Gamma(P,Q)}{\text{min}}\,\quad\underset{(\mathbf{X},\mathbf{X}^{\prime})\sim\pi}{\mathbb{E}}[c(\mathbf{X},\mathbf{X}^{\prime})].$
(3)

The solution of the Kantorovich problem, when it exists, is called the optimal
plan (Ambrosio et al., 2021, $\S$ 2.1). The set of push-forward maps, denoted
by $\Gamma_{0}(P,Q)$ and known as Monge couplings, are specific instances of
couplings characterized by $\pi\sim(X,T(X))$.

#### Kantorovich–Rubinstein Duality.

The OT problem, as introduced by Kantorovich, is formulated as a convex
problem on $\mathcal{P}(\mathcal{X}\times\mathcal{X}^{\prime})$, and its dual
is presented as a constrained concave maximization problem (Ambrosio et al.,
2021, $\S$ 3.1). The Kantorovich duality states that the minimum of the
Kantorovich problem is equal to the maximum of the dual problem in terms of
bounded and continuous Kantorovich potentials
$\varphi:\mathcal{X}\rightarrow\mathbb{R}$ and
$\psi:\mathcal{X}^{\prime}\rightarrow\mathbb{R}$:

$\underset{(\varphi,\psi)\in\Phi_{c}}{\text{sup}}\,\quad\underset{\mathbf{X}\sim
P}{\mathbb{E}}[\varphi(\mathbf{X})]+\underset{\mathbf{X}^{\prime}\sim
Q}{\mathbb{E}}[\psi(\mathbf{X}^{\prime})],$ (4)

over the constraint
$\Phi_{c}=\\{(\varphi,\psi):\varphi(x)+\psi(x^{\prime})\leq
c(x,x^{\prime})\\}$. To understand the duality concept, consider a scenario
where a logistics company aims to transport goods from a resource point $x$ to
a manufacturing location $x^{\prime}$. The company sets a loading fee
$\varphi(x)$ at $x$ and an unloading fee $\psi(x^{\prime})$ at $x^{\prime}$.
To be competitive, their profit margin, represented by
$\varphi(x)-\psi(x^{\prime})$, should not exceed the cost $c(x,x^{\prime})$.
To maximize their overall profits, the firm strategically modifies the pricing
functions $\varphi$ and $\psi$ in Eq. 4. (For more examples, see Galichon,
2018, $\S$ 2).

####  OT with Linear Constraints.

In practical applications, solutions to the OT problem often need to meet
certain constraints. In Zaev (2015) work, these conditions are incorporated by
introducing linear constraints. The OT problem with linear constraints seeks
optimal couplings with additional constraints over the family
$\mathcal{W}\subset\mathcal{C}_{L}(P,Q)$ of continuous and $L^{1}$-integrable
function on product space $\mathcal{X}\times\mathcal{X}^{\prime}$ such that:

$\displaystyle\underset{\pi\in\Gamma(P,Q)}{\text{min}}\,\quad\underset{(\mathbf{X},\mathbf{X}^{\prime})\sim\pi}{\mathbb{E}}[c(\mathbf{X},\mathbf{X}^{\prime})]\
\text{s.t.}\
\underset{(\mathbf{X},\mathbf{X}^{\prime})\sim\pi}{\mathbb{E}}[w(\mathbf{X},\mathbf{X}^{\prime})]=0,\quad\forall
w\in\mathcal{W}.$ (5)

Invariant and martingale OT are examples of this constrained problem.

#### Unbalanced OT.

In standard OT, the total supply equals the total demand. Unbalanced OT
extends the classical OT problem to scenarios where the total mass of the
source and target distributions are different or not necessarily preserved
during transport. Unbalanced optimal transport (UOT) incorporates additional
terms to account for the creation and annihilation of mass, providing a more
flexible framework for real-world applications Séjourné et al. (2022). Let
$\mu\in\mathcal{M}_{+}(\mathcal{X}),\nu\in\mathcal{M}_{+}(\mathcal{X}^{\prime})$
be two positive measures and
$\pi\in\mathcal{M}_{+}(\mathcal{X}\times\mathcal{X}^{\prime})$ be positive
measure on product space. The UOT problem can be formulated as:

$\displaystyle\underset{\pi\in\mathcal{M}_{+}(\mathcal{X},\mathcal{X}^{\prime})}{\text{min}}\,\int_{\mathcal{X}\times\mathcal{X}^{\prime}}c(x,x^{\prime})\,d\pi(x,x^{\prime})+\lambda_{1}D_{\varphi_{1}}(\mu,\pi_{1})+\lambda_{2}D_{\varphi_{2}}(\nu,\pi_{2}),$

where $D_{\varphi_{1}}$ and $D_{\varphi_{2}}$ are $\varphi-$divergence terms
for mass creation and annihilation, $(\lambda_{1},\lambda_{2})$ their hyper-
parameters, and $(\pi_{1},\pi_{2})$ are $\pi$’s marginal distributions on
spaces $\mathcal{X}$ and $\mathcal{X}^{\prime}$.

Figure 3: OT enhances security in counterfactual explanations by relying on
distribution and $\delta$-confidence, reducing the risk of decision boundary
manipulation. (a) Boundary manipulation using classical CE through circles.
(b) The manipulation strategy yields weaker results in identifying the
boundary in CCE.

#### Dynamic OT.

Dynamic optimal transport (DOT) is an extension of the classical OT problem,
incorporating a temporal dimension to consider the evolution of mass
distributions over time Benamou and Brenier (2000). It’s a framework for
studying how a distribution of mass or probability evolves optimally from one
state to another over a continuous time interval. Let
$\mu_{t}(\mathbf{x})=\mu(t,\mathbf{x})$ be path of probability measures such
that $\mu_{0}=P$ and $\mu_{1}=Q$ and $v_{t}(x)=v(t,\mathbf{x})$ be a velocity
field, then DOT can be formulated as follows:

$\underset{(\mu_{t},v_{t})}{\text{arg\,min}}\,\int_{0}^{T}\int_{\mathbb{R}^{n}}c(x,v_{t}(x,t),t)\,d\mu_{t}(x)\,dt,$

subject to the continuity equation $\partial\mu_{t}/\partial
t+\nabla\cdot(\mu_{t}v_{t})=0$, where, $c(x,x^{\prime},t)$ represents the cost
of transporting a unit mass from point $x$ to point $x^{\prime}$ at time $t$,
and the integral calculates the total cost used in moving the mass along the
path $\mu_{t}$. The continuity equation ensures the conservation of mass
during the transport.

##  3 Collective Counterfactual Explanations

This section introduces CCE, which addresses the side effects of
individualized CE by considering both the collective transition cost and
population probability measures.

###  3.1 Counterfactual Explanations Via OT

In practical applications, we can assume the feature space $\mathcal{X}$ is a
closed and bounded subset, with feature variables $\mathbf{X}$ and a known or
empirically estimated population probability measure $P$. The classifier $h$
partitions $\mathcal{X}$ into two subsets:
$\mathcal{X}^{-}=\\{x\in\mathcal{X}:h(x)=-1\\}$ and
$\mathcal{X}^{\textbf{+}}=\\{x\in\mathcal{X}:h(x)=+1\\}$. For simplicity, we
assume $h$ is constructed by $h(x)=\text{sign}(f(x))$, where
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a non-constant, Lipschitz
continuous score function that avoids assigning a single label to all points.
So, the decision boundary is obtained by the inverse image of $f$, i.e.,
$L=f^{-1}(0)$. In establishing OT, with no assumption about who wants to
change its state, naturally $\mathcal{X}^{-}$ can be considered as the source
space. Traditional CE aims to map points in $\mathcal{X}^{-}$ to a minimal
cost point on $L$. However, treating $L$ as the target space leads to a
concentration of points in low-density areas, altering the cost function. To
circumvent this issue, we expand the boundary $L$ to ensure individuals moving
there face no competition for limited resources. This approach leads to the
following definition.

######  Definition 1 ($\delta$-confidence region).

For classifier $h$ with above assumptions, the $\delta$-confidence region
consists of smallest super-level set with probability more than $\delta$,
i.e., $L^{\textbf{+}}_{\delta}:=\\{x\in\mathcal{X}^{\textbf{+}}:0<f(x)<r\\}$,
where $r=\inf\\{\lambda:P(0<f(\mathbf{X})<\lambda))>\delta\\}$. If
$\delta<P(\mathcal{X}^{\textbf{+}})$, then $L^{\textbf{+}}_{\delta}$ exists.

The choice of $\delta$ is guided by several intuitions: moving away from low-
density areas along $L$, ensuring societal stability and minimal impact on the
cost function when individuals move to the target area, and considerations of
robustness. By determining the source and target, we are now prepared to set
up a OT framework.

######  Definition 2 (Collective Counterfactual Explanations).

By the assumptions of Def. 1, let $P_{\textbf{-}}$ and $P_{\delta}$ be
truncated probability measures of $P$ over subsets of $\mathcal{X}^{-}$ and
$L^{\textbf{+}}_{\delta}$. A collective counterfactual explanations (CCE) plan
is a coupling $\pi\in\Gamma(P_{\textbf{-}},P_{\delta})$ that minimizes the
total cost of modification,

$\pi=\underset{\pi\in\Gamma(P_{\textbf{-}},P_{\delta})}{\text{arg\,min}}\,\quad\underset{(\mathbf{X},\mathbf{X}^{\prime})\sim\pi}{\mathbb{E}}[c(\mathbf{X},\mathbf{X}^{\prime})].$
(6)

When there is a plan CCE, CE for a given instance $x$ is produced by sampling
the conditional distribution $\pi$ given $x$, i.e.,
$x^{\prime}\sim\pi(\cdot|x)$.

Similarly to Monge’s problem, the CCE map is the measurable function
$T:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ that pushes forward
$P_{\textbf{-}}$ to $P_{\delta}$ and minimizes the total cost of modification:

$T=\underset{T:T_{*}P_{\textbf{-}}=P_{\delta}}{\text{arg\,min}}\,\quad\underset{\mathbf{X}\sim
P_{\textbf{-}}}{\mathbb{E}}[c(\mathbf{X},T(\mathbf{X}))].$ (7)

If $T$ exists, CE are obtained via $T(x)$.

The following theorem investigates the conditions for the existence of a CCE
plan and its relationship to the dual problem.

######  Theorem 1 (Villani, 2009, $\S$ 5).

If modification cost function
$c:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}\cup\\{\infty\\}$ is lower
semi-continuous and bounded from below, then the $CCE$ plan $\pi$ exists.
Moreover, the collective cost of the $CCE$ plan is equal to the supremum of
the dual problem,

$\underset{(\varphi,\psi)\in\Phi_{c}}{\text{sup}}\,\underset{\mathbf{X}\sim
P_{\textbf{-}}}{\mathbb{E}}[\varphi(\mathbf{X})]+\underset{\mathbf{X}^{\prime}\sim
P_{\delta}}{\mathbb{E}}[\psi(\mathbf{X}^{\prime})].$ (8)

Thm. 1’s assumption regarding the cost function is more general and
encompasses commonly used costs in the CE literature. Moreover, the dual
problem offers a competitive interpretation of CCE as a market equilibrium.

Figure 4: Using displacement interpolation, OT provides path-guided
counterfactual explanations.

####  CCE as Market Equilibrium.

Remember the loan example in the introduction: an individual in state $x$ is
rejected for a loan and must change to state $x^{\prime}=\textbf{CE}(x)$ to be
eligible. Let’s assume there are intermediate financial institutions that
purchase loans from banks at a price $\psi(x^{\prime})$ where
$\psi(x^{\prime})$ denotes the monetary risk value of issuing a loan to an
individual in state $x^{\prime}$ and selling these loans to applicants in
state $x$ at a price $\varphi(x)$. Applicants have the option to either
purchase loans directly or alter their state to become eligible for a loan.
The incentive to purchase rather than modify their state hinges on ensuring
the profit $\varphi(x)-\psi(x^{\prime})$ is lower than the cost of change
$c(x,x^{\prime})$. Otherwise, they would prefer altering their state over
buying a loan.

In a competitive market, financial institutions aim to maximize collective
profit, so reaching an equilibrium where profit is at its peak, expressed as:

$\sup_{(\varphi,\psi)\in\Phi_{c}}\underset{\mathbf{X}\sim
P_{\textbf{-}}}{\mathbb{E}}[\varphi(\mathbf{X})]+\underset{\mathbf{X}^{\prime}\sim
P_{\delta}}{\mathbb{E}}[-\psi(\mathbf{X}^{\prime})].$

According to Thm. 1, this equilibrium reflects the collective efficient effort
required for state modification. This fact is instrumental in evaluating the
subjective cost function concept. It suggests that in market equilibrium, the
cost can be inferred from the seller profit, which means
$c(x,x^{\prime})=\varphi(x)-\psi(x^{\prime})$.

While a CCE plan exists under minimal assumptions, CCE maps are more suited
for practical applications. The upcoming theorem delineates conditions under
which an optimal plan can be approximated by an optimal map, even when a
global optimal map solution is not available.

######  Proposition 1 (Gangbo, 1999).

By the assumptions of Def. 2, if $P_{\textbf{-}}$ is continuous measure then
optimal maps are dense in the space of couplings:

$\Gamma(P_{\textbf{-}},P_{\delta})=\overline{\Gamma_{0}(P_{\textbf{-}},P_{\delta})},$

where $\Gamma_{0}$ is set of Monge couplings and $\overline{\Gamma_{0}}$ is
its set closure. If $c$ is continuous and bounded then:

$\underset{\pi\in\Gamma}{\text{min}}\,\
\underset{(\mathbf{X},\mathbf{X}^{\prime})\sim\pi}{\mathbb{E}}[c(\mathbf{X},\mathbf{X}^{\prime})]=\underset{T\in\Gamma_{0}}{\text{min}}\,\
\underset{\mathbf{X}\sim
P_{\textbf{-}}}{\mathbb{E}}[c(\mathbf{X},T(\mathbf{X}))].$

An immediate consequence of Prop. 1 is that CCE plans and estimated CCE maps
can converge arbitrarily closely in terms of total cost. Finally, we outline
the sufficient conditions for the existence of a CCE map.

######  Theorem 2 (Villani, 2009, $\S$ 10).

By Def. 2 assumptions, let $P$ be continuous measure, and assume that $c$ is
twice continuously differentiable and for every compact set
$K\subseteq\mathbb{R}^{d}$, there is $c_{K}>0$ such that for every
$x_{1},x_{2}\in\mathcal{X}$,

$\sup_{x^{\prime}\in\mathcal{X}^{\prime}}|c(x_{1},x^{\prime})\textendash
c(x_{2},x^{\prime})|\leq c_{K}|x_{1}\textendash x_{2}|,$

and also has the twist condition holds, i.e., for every $x\in\mathcal{X}$ and
$x^{\prime}_{1},x^{\prime}_{2}\in\mathcal{X}^{\prime}$,

$\nabla_{x}c(x,x^{\prime}_{1})=\nabla_{x}c(x,x^{\prime}_{2})\Rightarrow
x^{\prime}_{1}=x^{\prime}_{2}.$

Then by the assumptions of Def. 2 we have:

  1. 1.

Both CCE plan and map exist and the collective cost of them are equal.
Moreover, the plan is unique and is generated by map $T$;

  2. 2.

The solution $(\varphi,\psi)$ of the dual problem is unique (up to an additive
constant), and $\varphi$ is almost-everywhere differentiable. Moreover
$\psi(x^{\prime})=\underset{x\in\mathcal{X}}{\inf}(c(x,x^{\prime})\textendash\varphi(x))$;

  3. 3.

$T$ and $\varphi$ are related by $T(x)=e_{x}(\nabla\varphi(x)),$ where
$e_{x}(\cdot):=\nabla_{x}c(x,\cdot)^{\textendash 1}$.

###  3.2 Collective Counterfactual Explanations Desiderata

In this section, we explore how CCE augments traditional CE by integrating a
list of additional desiderata as mentioned in surveys Verma et al. (2020);
Karimi et al. (2022).

#### Actionability.

Actionability refers to the feasibility of implementing proposed changes in
counterfactual scenarios Ustun et al. (2019). Variables such as birthplace,
immutable, or age, characterized by limited modifiability and a naturally
increasing trajectory, should be regarded as constraints when designing CE.
One method to incorporate these constraints is to embed them into the cost
function. For instance, the cost of changing gender can be set sufficiently
high. Alternatively, constraints can be formulated as functions and integrated
to CCE by using extensions of OT with linear constraints. This approach is
feasible because conditions, encoded as Boolean functions, which are the
limits of bounded continuous functions, therefore can be employed as linear
constraints. The following proposition suggests that the existence of a CCE
plan is assured, as long as there is at least one coupling that meets the
specified linear constraints.

######  Proposition 2 (Zaev, 2015).

With conditions of Prop. 1, the OT problem with linear constraints has a
solution if and only if the set $\Gamma_{\mathcal{W}}:=\\{\pi\in\Gamma:\int
wd\pi=0,w\in\overline{\mathcal{W}}\\}$ is not empty, where
$\overline{\mathcal{W}}$, the closure of $\mathcal{W}$ in the $C_{L}$
topology.

#### Sparsity.

To achieve sparsity in CE—which focuses on minimizing feature alterations in
recommendations state $x^{\prime}$—a sparsity objective function
$g(x^{\prime}-x)$, such as $L_{0}/L_{1}$ norms as suggested by Miller (2019),
can be employed. This approach is also applicable to CCE by integrating the
sparsity objective into the cost function, a concept further detailed in
Cuturi et al. (2023).

Figure 5: For the Gaussian mixture with means $\mu_{1}=[0.3,0.3]$,
$\mu_{2}=[0.7,0.7]$, and a covariance of $\text{diag}(0.2,0.2)$ and classifier
with $f(x)=x_{1}+x_{2}+\frac{1}{5}sin(2\pi x_{1})-1$. (a) displays the
probability measure $P$ alongside the $\delta$-confidence region, (b)
illustrates the negatively induced probability measure $P_{-}$, and (c) shows
the probability measure induced by the $\delta$-confidence region
$L^{\textbf{+}}_{\delta}$.

#### Amortized Inference.

Traditional CE generation typically involves an optimization process for each
instance. In contrast, amortized inference enhances CE generation by utilizing
patterns learned from previous instances, eliminating the need for individual
optimization for each new data point. When a CCE map $T$ exists, numerical
methods can be used to approximate it. This approximation function is not
limited to instances present in the training data; it can be extended to other
cases as well. $\S$ 5 presents an efficient algorithm for estimating the CCE
map.

#### Data Manifold Closeness.

To ensure credibility and avoid outlier positions, a CE should resemble the
training dataset and preserve feature correlations. In CCE, recommendations
are spread across the data manifold by targeting a significant part of the
positive region, thereby avoiding minimum-cost outliers. Additionally,
adjusting UOT hyperparameters can steer recommendations toward denser areas.

#### Security & Privacy.

The classical CE API poses a risk of facilitating surrogate model construction
through low-budget queries, as illustrated in Fig. 3(a). Each CE suggests that
instances within a radius $c(x,\textbf{CE}(x))$ centered at $x$ belong to the
negative region, enabling swift manipulation of the decision boundary.
Conversely, in CCE, the recommendations are influenced by the choice of
$L^{\textbf{+}}_{\delta}$, the probability measure $P$, and the collective
modification cost considered in the optimization process. Combining these
diverse elements in the search for CCE complicates the boundary manipulation,
as shown in Fig. 3(b).

#### Robustness & Individual Fairness.

In real-world decision-making applications, robustness and individual fairness
are crucial aspects of the algorithm. They ensure that individuals with
similar features receive similar recommendations. In recent work, Otto et al.,
2021 (Cor. 1.4) shows that under mild conditions on the probability density
function $P$ and for a general family of cost functions, the CCE map
constitutes a diffeomorphism. Given that the space $\mathcal{X}$ is compact,
it implies that minor modifications in state result in only slight variations
in the explanations. The classical CE does not consider these two properties.
Regarding the boundary shape, it may result in jumps in $\textbf{CE}(x)$
recommendations.

Figure 6: Demonstrates controlling additional costs in CCE using UOT by
adjusting the parameter $\lambda_{2}$. When $\lambda_{2}=0$ (d), the solution
resembles standard CE (b), while setting $\lambda_{2}=1$ (f) aligns the
estimation of CCE with the Sinkhorn method (c).

##  4 Counterfactual Explanations Extensions

In this section, we discuss some extensions of original CE and provide
solutions through the OT framework.

###  4.1 Path-Guided Counterfactual Explanations

Applying the standard CE method faces challenges due to the uncertain
trajectory towards the final modification state, especially in prolonged
processes like obtaining a loan. Path-guided CE overcomes this by
demonstrating the progression towards the goal through temporal steps, similar
to displacement interpolation in dynamic OT. This interpolation, a continuous-
time version of OT, connects initial and final probability distributions. The
subsequent theorem emphasizes that in a metric space $(\mathcal{X},d)$ with
geodesics, if a CCE map exists, displacement interpolation generates geodesics
in the probability measure space, relying on the geodesics of the underlying
space $\mathcal{X}$ between $x$ and the CCE map $T(x)$ Khan and Zhang (2022).

######  Theorem 3 (Villani, 2009, $\S$ 7).

Let $(\mathcal{X},d)$ be a Polish, locally compact length space. Then, given
any two measures $P,Q\in\mathcal{P}(\mathcal{X})$, and a continuous curve
$(\mu_{t})_{0\leq t\leq 1}\in$ $\mathcal{P}(\mathcal{X})_{p}$, the following
properties are equivalent:

i. $\mu_{t}$ is the push-forward by the flow of minimizing, constant-speed
geodesics from $x$ to $T(x)$, where $T$ is the solution to the OT map.

ii. $(\mu_{t})_{0\leq t\leq 1}$ is a geodesic curve in the space
$\mathcal{P}(\mathcal{X})_{p}$.

In the case of the squared-distance cost function, the displacement
interpolation $\mu_{t}$ for $t\in[0,1]$ is defined by:

$\mu_{t}=((1-t)\text{id}+tT)_{*}P,$

where id is the identity map on $\mathcal{X}$. Let
$T_{t}(x):[0,1]\to\mathcal{X}$ denote the constant speed geodesic connecting
$x$ to $T(x)$. Then, according to Thm. 3, the optimal geodesic in the
probability space is achieved by $\mu_{t}=(T_{t})_{*}P$. This means that the
optimal flow maps $x$ along the geodesics $T_{t}(x)$ towards its target. Such
a parametrization is instrumental in constructing path-guided CE based on
discrete temporal steps that can be seen in Fig. 4 and Fig. 7.

###  4.2 CE for Ordered Classifiers

In traditional CE, it is commonly assumed that eligibility is determined by a
single classifier. However, in real-world scenarios, like loan applications,
eligibility may vary across different loan amounts. For example, a loan
request for the amount $l$ is rejected, but for a lower amount $l^{\prime}$,
it will be accepted. As a result, distinct classifiers are employed at each
loan value to make individualized decisions. It means for each loan amount
$l$, there exists a classifier
$h_{l}:\mathcal{X}\subset\mathbb{R}^{d}\rightarrow\mathcal{Y}=\\{\pm 1\\}$.
This family of classifiers follows a natural order, i.e., $l_{1}<l_{2}\implies
h_{l_{1}}(x)\geq h_{l_{2}}(x)$.

Figure 7: Utilizing the back-and-forth method, the CCE map is estimated, and
through constructing the displacement interpolation, the optimal flow that
pushes forward the negative area into the $\delta$-confidence region is
visualized.

In this scenario, the CE problem aims to identify the optimal state
modification $x^{\prime}$ for each $x$ in $\mathcal{X}$ and $l$ in space of
loan values $\mathcal{T}=[0,B]$:

$\textbf{CE}_{l}(x):=\underset{x^{\prime}\in\mathcal{X}}{\text{inf}}\
c(x,x^{\prime})\ \text{s.t.}\ h_{l}(x^{\prime})=1.$

To formulate the ordered classifier version of CE within the OT framework, in
addition to the feature space probability $P$, it is essential to know $T$,
the loan amounts probability distribution. Let $S$ represent the joint
probability of features and loans over the space
$\mathcal{S}=(\mathcal{X},\mathcal{T})$.

For simplicity assume $h_{l}$ is defined as the sign of a real-valued,
continuously differentiable function
$f:\mathcal{X}\times\mathcal{T}\rightarrow\mathbb{R}$ at level $l$, i.e.,
$h_{l}(x)=\text{sign}(f(x,l))$. To maintain the properties of an ordered
classifier, $h$ should be a decreasing function with respect to $l$. This is
equivalent to assuming that $\frac{\partial f}{\partial l}<0$. Under this
formulation, the decision boundary for the family of classifiers is determined
as follows:

$L=\\{(x,t)\in(\mathcal{X},\mathcal{T}):f(x,t)=0\\}.$

Given that $\nabla f$ is non-zero everywhere, the implicit function theorem
(Lang, 2012, $\S$ 5) states that $L$ is a $C^{1}$ $n$-dimensional manifold in
an $(n+1)$-dimensional space $\mathcal{S}$. It is presumed that the
probability measure $S$ is absolutely continuous, hence it possesses a
probability density function $g$. We select $L$ as the target space. Given the
presence of a continuous family of ordered classifiers, it is reasonable to
assume that the margins are sufficiently enriched to map points to these
classifiers. Furthermore, the concept of buffering the decision boundary of
$L^{\textbf{+}}_{\delta}$ by a margin $\delta$ is applicable in this context.
To derive the target probability measure $S_{L}$ on surface $L$, we apply the
coarea formula (Evans, 2018, $\S$ 3) with the density function $g(s)$ leading
to:

$S_{L}(A)=\frac{1}{Z}\int_{A}g(s)\frac{1}{|\nabla
f(s)|}\,d\mathcal{H}^{n}(s)\quad\text{for }A\subset L,$

here $\mathcal{H}^{n}$ is the (n)-dimensional Hausdorff measure, $|\nabla
f(s)|$ is the norm of the gradient of $f$ at $s$, and $Z$ is a normalization
factor ensuring that $S_{L}$ is a probability measure on $L$, with
$S_{L}(L)=1$. The term $\frac{1}{|\nabla f(s)|}$ ensures proper scaling along
the surface $L$.

To construct the source measure, let $S^{-}$ be the restriction of the
probability measure $S$ to the space $\\{(x,t)\in\mathcal{S}:f(x,t)<0\\}$. OT
can now be applied to the probability measures $S^{-}$ and $S_{L}$ with the
cost function $c^{*}((x,t),(x^{\prime},t^{\prime}))=c(x,x^{\prime})$. To
ensure that in the explanations $(x^{\prime},t^{\prime})$, $t^{\prime}\geq t$,
a constraint $1_{\\{T^{\prime}>T\\}}$ is added to OT problem where
$T,T^{\prime}$ are the values of a loan in source and target densities.
Utilizing OT with linear constraints enables the identification of the most
efficient plan and the creation of CCE for ordered classifiers.

###  4.3 Utility-based Counterfactuals Explanations

To better align CE with human decision-making contexts like loan applications,
it is advantageous to restructure it using scenarios based on human utility
functions. In classical CE, the utility function of an individual at state $x$
who wishes to modify their state to $x^{\prime}$ for loan approval can be
expressed as follows:

$u(x,x^{\prime})=\frac{1}{2}(1+h(x^{\prime}))l-c(x,x^{\prime}),$ (9)

here, $l$ is the additional utility from a positive decision. It is a natural
assumption that for all $x^{\prime}$ satisfying $h(x^{\prime})=1$, the
expression $l-c(x,x^{\prime})$ is not negative. By this formulation it is easy
to see that the optimal solution of Eq. 9 is found by minimizing the cost
function, $\text{argmin}\ c(x,x^{\prime});\ h(x^{\prime})=1$, aligning with
the counterfactual optimization problem as described in Eq. 1. The utility
function in Eq. 9 represents a simplified real-world scenario by assuming that
transitioning to state $x^{\prime}$ does not alter individual utility and
ignores the fact that certain actions, such as increasing saving, enhance
utility independently of the loan’s approval. Therefore, we can assign utility
to each state: $u:\mathcal{X}\rightarrow\mathbb{R}$. Then moving to a better
creditworthy condition $x^{\prime}$, not only entails a cost
$-c(x,x^{\prime})$ but also yields an additional utility
$\eta(x,x^{\prime})\coloneqq u(x^{\prime})-u(x)$. Consequently, a more nuanced
model can be expressed as follows:

$u(x,x^{\prime})=\frac{1}{2}(1+h(x^{\prime}))l-c(x,x^{\prime})+\eta(x,x^{\prime}).$
(10)

In Eq. 10, it is assumed that the individual requires a loan to address their
financial challenges, i.e., for all $x^{\prime}$ where $h(x^{\prime})=-1$, it
holds that $-c(x,x^{\prime})+\eta(x,x^{\prime})<0$. Additionally, within a
finite time frame (e.g., one year), the net benefit of changing states, given
by $\eta(x,x^{\prime})-c(x,x^{\prime})$, is not always positive. For example,
the cost of acquiring savings of 1 million dollars for an average worker goes
to infinity, while the utility $\eta(x,x^{\prime})$ remains finite. Hence, the
optimal solution can be formulated as:

$\underset{x^{\prime}:h(x^{\prime})=1}{\text{arg\,max}}\,l-c(x,x^{\prime})+\eta(x,x^{\prime}),$
(11)

This advanced formulation accounts for both the costs and benefits associated
with changing states, providing a more nuanced understanding of the decision-
making process in financial contexts. In Eq. 11, reporting a single point as a
CE is not obligatory. However, if a customer aims for a minimum utility of
$u_{o}$, then any points of set
$\mathcal{X}^{\prime}_{x}=\\{x^{\prime}:l-c(x,x^{\prime})+\eta(x,x^{\prime})>u_{0}\\}$
can be selected as a target. If the
$\mathcal{X}^{\prime}=\underset{x\in\mathcal{X}^{-}}{\bigcap}\mathcal{X}^{\prime}_{x}$
is non-zero measurable set (assume a case that the value of $u^{+}$ is more
than cost $c(x,x^{\prime})$ alongside decision boundary $L$, so
$L\subset\mathcal{X}^{\prime}$), then we can define OT problem for
$\mathcal{X}^{-}$ and $\mathcal{X}^{\prime}$ spaces.

We extend this model with a temporal dimension in $\S$ A.3. In particular, we
consider a setting where utilities are not observable and demonstrate that
taking the underlying density of individuals into account, as we do in our
optimal transport formulation, gives the optimal action.

Figure 8: Blue curve: Represents the percentage increase in cost relative to
standard CE as $\lambda_{2}$ varies from 0 to 1. Red curve: Depicts the
logarithmic value of the KL-divergence between the UOT solution and the target
distribution.

##  5 Computational Methods

This section explores the estimation of the CCE map and plan for the OT and
UOT problems, incorporating a comparative analysis of the incremental costs
associated with CCE solutions relative to the classical CE approach. To see
the results more concretely, we use a $2$-D Gaussian mixture model with means
$\mu_{1}=[0.3,0.3]$ and $\mu_{2}=[0.7,0.7]$, and a covariance matrix of
$\text{diag}(0.2,0.2)$ for synthetic data. Our study examines three distinct
scenarios, each characterized by a unique classifier function:
$f_{1}(x)=x_{1}+x_{2}+\frac{1}{5}\sin(2\pi x_{1})-1$,
$f_{2}(x)=x_{1}+x_{2}-1$, and
$f_{3}(x)=(x_{1}-0.5)^{2}+(x_{2}-0.5)^{2}-\frac{1}{9}$. For each classifier,
we analyze both the negative region and the $\delta$-confidence region for
$\delta=0.2$, along with their respective probability densities, as shown in
Fig. 5.

In each experiment, 100 samples are generated across source $P_{-}$ and target
$P_{+\delta}$ distributions (see Fig. 6.a). The CE for each negative area
point is determined by identifying the nearest point using the $L_{2}$
Euclidean cost function, as depicted in Fig. 6.b. The CCE plan is efficiently
estimated using the Sinkhorn algorithm, recognized as the one of the fastest
method for resolving the OT problem, with a complexity of $O(n^{2})$ per
iteration (detailed in $\S$ B.1). Additionally, Fig. 6.c illustrates the
solution of the CCE problem.

The OT solution’s total cost increases due to the necessity to encompass all
points within the $\delta$-confidence region. Employing UOT with the tunable
$\lambda_{2}$ parameter enables a balance between total cost and coverage of
the target source. We use the unbalanced Sinkhorn algorithm for solving UOT
problem (see details in $\S$ B.2). Fig. 6.d shows that setting $\lambda_{2}=0$
yields a result similar to classical CE. Increasing $\lambda_{2}$ makes the
solution more balanced, as shown in parts e and f of Fig. 6.

In Fig. 8 on the left, we examine the UOT solution’s behavior when
$\lambda_{2}$ varies from $0$ to $1$. The blue curve represents the extra cost
percentage compared to the classical CE solution. The red curve indicates the
logarithmic KL-divergence between the UOT solution’s marginal and the target
distribution. As we expect, the Fig. 8 shows that increasing $\lambda_{2}$
makes the solution more closely match the target distribution, but at a higher
total cost.

If the cost function is $c(x,y)=h(y-x)$, where $h$ is a strictly convex and
even function, the CCE map can be computed using the back-and-forth algorithm
(BFM). The BFM , is a very fast method that requires $O(n)$ storage and
$O(n\log(n))$ computation per iteration. For more details, refer to $\S$ B.3.
The solution obtained from BFM allows for displacement interpolation
parameterized by $t\in[0,1]$, as shown in the Fig. 4 alongside the CCE map
solution. The details of other two scenarios can be found in $\S$ B.4.

##  6 Discussion and Future works

CCE is specifically designed to address the limitations of traditional CE. It
takes into account the probability distribution of the population, effectively
preventing the significant performative impacts that may arise from
substantial shifts in individual behavior on the cost function. Our research
demonstrates that OT theory, along with its computational tools, is highly
effective in addressing various concerns in the design of counterfactuals.

To compute the CE solution, the estimation of $P$ is not crucial. The Sinkhorn
methods facilitate working with sample distributions rather than entire
distributions. Moreover, the assumption regarding $h$ is not obligatory; we
can employ black-box models for classifiers. A further challenge involves
implementing CCE with real-world datasets. As we can only report additional
costs relative to classical CE, we confine our analysis to synthetic examples.
These examples intuitively demonstrate the impact of CCE on recommendations.

Selecting the appropriate value for $\delta$ requires a mathematical framework
for modeling the externality’s impact on the cost function, which requires a
separate research endeavor. In our exploration of the various desiderata of
classical CE, we omitted a discussion on causality. Given the evolving theory
of optimal transport within a causal framework, this subject will be the focus
of our subsequent research.

## References

  * Ambrosio et al. (2013) Luigi Ambrosio, Alberto Bressan, Dirk Helbing, Axel Klar, Enrique Zuazua, Luigi Ambrosio, and Nicola Gigli.  A user’s guide to optimal transport.  _Modelling and Optimisation of Flows on Networks: Cetraro, Italy 2009, Editors: Benedetto Piccoli, Michel Rascle_ , pages 1–155, 2013.
  * Ambrosio et al. (2021) Luigi Ambrosio, Elia Brué, Daniele Semola, et al.  _Lectures on optimal transport_.  Springer, 2021.
  * Benamou and Brenier (2000) Jean-David Benamou and Yann Brenier.  A computational fluid mechanics solution to the monge-kantorovich mass transfer problem.  _Numerische Mathematik_ , 84(3):375–393, 2000.
  * Carrizosa et al. (2024) Emilio Carrizosa, Jasone Ramírez-Ayerbe, and Dolores Romero Morales.  Generating collective counterfactual explanations in score-based classification via mathematical optimization.  _Expert Systems with Applications_ , 238:121954, 2024.
  * Cuturi et al. (2023) Marco Cuturi, Michal Klein, and Pierre Ablin.  Monge, bregman and occam: Interpretable optimal transport in high-dimensions with feature-sparse maps.  _arXiv preprint arXiv:2302.04065_ , 2023.
  * Evans (2018) LawrenceCraig Evans.  _Measure theory and fine properties of functions_.  Routledge, 2018.
  * Galichon (2018) Alfred Galichon.  _Optimal transport methods in economics_.  Princeton University Press, 2018.
  * Gangbo (1999) Wilfrid Gangbo.  The monge mass transfer problem and its applications.  _Contemporary Mathematics_ , 226:79–104, 1999.
  * Gross and De Dreu (2019) Jörg Gross and Carsten KW De Dreu.  Individual solutions to shared problems create a modern tragedy of the commons.  _Science advances_ , 5(4):eaau7296, 2019.
  * Guidotti (2022) Riccardo Guidotti.  Counterfactual explanations and how to find them: literature review and benchmarking.  _Data Mining and Knowledge Discovery_ , pages 1–55, 2022.
  * Jacobs and Léger (2020) Matt Jacobs and Flavien Léger.  A fast approach to optimal transport: The back-and-forth method.  _Numerische Mathematik_ , 146(3):513–544, 2020.
  * Kanamori et al. (2020) Kentaro Kanamori, Takuya Takagi, Ken Kobayashi, and Hiroki Arimura.  Dace: Distribution-aware counterfactual explanation by mixed-integer linear optimization.  In _IJCAI_ , pages 2855–2862, 2020.
  * Kantorovich (1942) L Kantorovich.  On the transfer of masses (in russian).  In _Doklady Akademii Nauk_ , volume 37, page 227, 1942.
  * Karimi et al. (2022) Amir-Hossein Karimi, Gilles Barthe, Bernhard Schölkopf, and Isabel Valera.  A survey of algorithmic recourse: contrastive explanations and consequential recommendations.  _ACM Computing Surveys_ , 55(5):1–29, 2022.
  * Khan and Zhang (2022) Gabriel Khan and Jun Zhang.  When optimal transport meets information geometry.  _Information Geometry_ , 5(1):47–78, 2022.
  * Lang (2012) Serge Lang.  _Fundamentals of differential geometry_ , volume 191.  Springer Science & Business Media, 2012.
  * Miller (2019) Tim Miller.  Explanation in artificial intelligence: Insights from the social sciences.  _Artificial intelligence_ , 267:1–38, 2019.
  * Monge (1781) G Monge.  Mémoire sur la théorie des déblais et des remblais, histoire de l’académie royale des sciences de paris, avec les mémoires de mathématique et de physique pour la même année, 666-704, 1781.  _Cited on_ , page 8, 1781.
  * Otto et al. (2021) Felix Otto, Maxime Prod’Homme, and Tobias Ried.  Variational approach to regularity of optimal transport maps: general cost functions.  _Annals of PDE_ , 7(2):17, 2021.
  * Peyré et al. (2017) Gabriel Peyré, Marco Cuturi, et al.  Computational optimal transport.  _Center for Research in Economics and Statistics Working Papers_ , 2017.
  * Poyiadzi et al. (2020) Rafael Poyiadzi, Kacper Sokol, Raul Santos-Rodriguez, Tijl De Bie, and Peter Flach.  Face: feasible and actionable counterfactual explanations.  In _Proceedings of the AAAI/ACM Conference on AI, Ethics, and Society_ , pages 344–350, 2020.
  * Séjourné et al. (2022) Thibault Séjourné, Gabriel Peyré, and François-Xavier Vialard.  Unbalanced optimal transport, from theory to numerics.  _arXiv preprint arXiv:2211.08775_ , 2022.
  * Ustun et al. (2019) Berk Ustun, Alexander Spangher, and Yang Liu.  Actionable recourse in linear classification.  In _Proceedings of the conference on fairness, accountability, and transparency_ , pages 10–19, 2019.
  * Verma et al. (2020) Sahil Verma, Varich Boonsanong, Minh Hoang, Keegan E Hines, John P Dickerson, and Chirag Shah.  Counterfactual explanations and algorithmic recourses for machine learning: A review.  _arXiv preprint arXiv:2010.10596_ , 2020.
  * Villani (2009) Cédric Villani.  _Optimal transport: old and new_ , volume 338.  Springer, 2009.
  * Wachter et al. (2017) Sandra Wachter, Brent Mittelstadt, and Chris Russell.  Counterfactual explanations without opening the black box: Automated decisions and the gdpr.  _Harv. JL & Tech._, 31:841, 2017.
  * Zaev (2015) Danila A Zaev.  On the monge–kantorovich problem with additional linear constraints.  _Mathematical Notes_ , 98:725–741, 2015.

#### Notational Conventions.

Random variables are represented by bold letters, for instance, $\mathbf{X}$,
their spaces by calligraphic letters like $\mathcal{X}$, and probability
measures by uppercase letters, such as $P$. The space of probability measures
on $\mathcal{X}$ is denoted as $\mathcal{P(X)}$, and positive measures on
$\mathcal{X}$ as $\mathcal{M}_{+}(\mathcal{X})$. For enhanced clarity in
presenting Optimal Transport definitions and theorems, probability notations
are employed rather than measure theory notations.

Symbol | Notion
---|---
CE |  Counterfactual Explanations
OT |  Optimal Transport
$\mathcal{X}$ |  Feature space s.t. $\mathcal{X}\subseteq\mathbb{R}^{d}$
$\mathcal{Y}$ |  Label space
$h$ |  Decision classifier $h:\mathcal{X}\rightarrow\mathcal{L}$
$c:\mathcal{X}\times\mathcal{X}^{\prime}\rightarrow\mathbb{R}$ |  Modification cost function
$\textbf{CE}(\mathbf{x})$ |  Counterfactual explanations
$\mathcal{A}$ |  Space of feasible actions
$L$ |  Decision boundary of $h$
$P$ |  Population probability measure
$\mathcal{P}(\mathcal{X})$ |  Space of probability measures
$\mathcal{M}_{+}(\mathcal{X})$ |  Space of positive measures
$c:\mathcal{X}\times\mathcal{X}^{\prime}\rightarrow\mathbb{R}$ |  ground-cost function
$T:\mathcal{X}\rightarrow\mathcal{X}^{\prime}$ |  Measurable map from $\mathcal{X}$ to $\mathcal{X}^{\prime}$
$T_{*}P$ |  Push-forward probability of $Q$
$\Gamma(P,Q)$ |  Couplings between $P$ and $Q$
$\Gamma_{0}(P,Q)$ |  Monge couplings
$\varphi:\mathcal{X}\rightarrow\mathbb{R}$ |  Kantorovich potential function
$\psi:\mathcal{X}^{\prime}\rightarrow\mathbb{R}$ |  Kantorovich potential function
$\Phi_{c}$ |  Space of constraint potentials
$\mathcal{C}_{L}(P)$ |  Space of integrable continuous functions
$\mathcal{C}_{L}(P,Q)$ |  Space of linear constraints on $\mathcal{X}\times\mathcal{X}^{\prime}$
UOT |  Unbalanced optimal transport
$D_{\varphi_{1}}$ |  $\varphi-$divergence measure
$\pi_{1},\pi_{2}$ |  Marginals distributions of $\pi$ over $\mathcal{X}$ and $\mathcal{X}^{\prime}$
$\lambda_{1},\lambda_{2}$ |  $\varphi-$divergence hyper parameters
DOT |  Dynamic optimal transport
$\mu_{t}$ |  Temporal path on $\mathcal{P}(\mathcal{X})$
$v_{t}$ |  Velocity field on $\mathcal{P}(\mathcal{X})$
$(\mathcal{X},d)$ |  Metric space
$\mathcal{P}_{p}(\mathcal{X})$ |  Probability measures with finite $p$-th moment
$W_{p}(P,Q)$ |  $p$-Wasserstein distance
$\mathcal{X}^{-}$ |  Negative part of $\mathcal{X}$
$\mathcal{X}^{\textbf{+}}$ |  Positive part of $\mathcal{X}$
$f(x)$ |  Score function s.t. $h=\text{sign}(f)$
$L^{\textbf{+}}_{\delta}$ |  $\delta$-confidence region
$P_{\textbf{-}}$ |  Truncated probability on $\mathcal{X}^{-}$
$P_{\delta}$ |  Truncated probability on $L^{\textbf{+}}_{\delta}$
$\pi^{*}$ |  Optimal $CCE$ plan
$T^{*}$ |  Optimal $CCE$ map
$\nabla$ |  Gradient operator
$\mathcal{W}$ |  Linear constraint functions
$T_{t}(x)$ |  Constant speed geodesic from $x$ to $T(x)$
$\mu_{t}$ |  Displacement interpolation
$l$ |  Levels of classifier function
$h_{l}$ |  Classifier at level $l$
$\mathcal{T}$ |  The space of classifier levels
$T$ |  Probability measure over $\mathcal{T}$
$\mathcal{S}$ |  product space of $\mathcal{X}$ and $\mathcal{T}$
$S$ |  Probability measure over $\mathcal{S}$
$\textbf{CE}_{l}(x)$ |  Ordered CE problem
$g$ |  $S$ probability density function
$\mathcal{H}^{n}$ |  (n)-dimensional Hausdorff measure
$S_{L}$ |  Probability measure over $L$
$u(x,x^{\prime})$ |  Individual utility function of going from state $x$ to $x^{\prime}$
BFM |  Back-and-Forth method
Fig. |  Figure
Thm. |  Theorem
Eq. |  Equation
Prop. |  Proposition
Def. |  Definition

##  Appendix A Supplementary Materials

###  A.1 Definitions

######  Definition 3 (Push-forward Measure).

Let $P$, $Q$ be two probability measures in $\mathcal{P}(\mathbb{R}^{d})$ and
$T:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is map, the measure $Q$ is called
the push-forward of $P$ through $T$ is denoted by $T_{*}P$ if:

$Q(B)=P(T^{-1}(B)),\quad\forall B\subset\mathbb{R}^{d}$

######  Definition 4 (Weak Topology).

Weak topology on a space of probability measures on $\mathbb{R}^{d}$, denoted
by $\mathcal{P}(\mathbb{R}^{d})$, is defined by convergence in distribution. A
sequence of probability measures $(P_{n})_{n\in\mathbb{N}}$ in
$\mathcal{P}(\mathbb{R}^{d})$ converges weakly to a probability measure $P$ if
for every bounded continuous function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$,

$\displaystyle\lim_{n\to\infty}\underset{\mathbf{X}_{n}\sim
P_{n}}{\mathbb{E}}[f(\mathbf{X}_{n})]=\underset{\mathbf{X}\sim
P}{\mathbb{E}}[f(\mathbf{X})].$

######  Definition 5 (Set of Couplings).

The set $\Gamma(P,Q)$ represents the couplings of probability distributions
$P,Q\in\mathcal{P}(\mathbb{R}^{d})$, comprising distributions over
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ with margins $P$ and $Q$. A measure $\pi$
belongs to $\Gamma(P,Q)$ if and only if

$\pi(A\times\mathbb{R}^{d})=P(A)\quad\text{and}\quad\pi(\mathbb{R}^{d}\times
B)=Q(B)\quad\forall A,B\subset\mathbb{R}^{d}$

By extension, a random pair $(X,Y)\sim\pi$, where $\pi\in\Gamma(P,Q)$, will
also be called a coupling of $P$ and $Q$.

######  Definition 6 (Linear Constraints).

Let $\mathcal{C}_{L}(P)$ and $\mathcal{C}_{L}(Q)$ be the space of continuous
functions on $\mathcal{X}$ and $\mathcal{X}^{\prime}$ and $L^{1}$-integrable
respect to $P$ and $Q$ measures. Let $\mathcal{C}_{L}(P,Q)$ be the family of
continuous functions on $\mathcal{X}\times\mathcal{X}^{\prime}$ such that:

$\mathcal{C}_{L}(P,Q)=\\{h\in\mathcal{C}(\mathcal{X}\times\mathcal{X}^{\prime}):\exists
f=f_{1}+f_{2}\ \ \text{s.t.}\ \ |h|\leq f\\},$

where $f_{1}\in\mathcal{C}_{L}(P)$, $f_{2}\in\mathcal{C}_{L}(Q)$.

######  Definition 7 (c-cyclical monotonicity (Ch. 4 Ambrosio et al. (2021))).

A probability measure $\pi$ on $X\times Y$ satisfies $c$-cyclical monotonicity
if, for any finite set $(x_{i},y_{i})_{i=1}^{n}$ in its support and any
permutation $\sigma$, the inequality
$\sum_{i=1}^{n}c(x_{i},y_{i})\leq\sum_{i=1}^{n}c(x_{i},y_{\sigma(i)})$ holds.

######  Definition 8 (c-transforms; c-superdifferential ( Ambrosio et al.,
2021, $\S$.4)).

Let $\varphi:\mathbb{R}^{d}\rightarrow\mathbb{R}\cup\\{\pm\infty\\}$ be a
function and $c$ a cost function. The $c_{+}$-transform of $\varphi$, denoted
$\varphi^{c_{+}}:\mathbb{R}^{d}\rightarrow\mathbb{R}\cup\\{-\infty\\}$, is
defined by:

$\varphi^{c_{+}}(x^{\prime}):=\inf_{x\in\mathbb{R}^{d}}\\{c(x,x^{\prime})-\varphi(x)\\}$

The $c$-superdifferential of $\varphi$,
$\partial^{c+}\varphi\subset\mathbb{R}^{d}\times\mathbb{R}^{d}$, is given by:

$\partial^{c+}\varphi:=\\{(x,x^{\prime})\in\mathbb{R}^{d}\times\mathbb{R}^{d}:\varphi(x)+\varphi^{c_{+}}(x^{\prime})=c(x,x^{\prime})\\}$

######  Definition 9 (The Wasserstein Metric).

In the metric space $(\mathcal{X},d)$, OT naturally defines a metric known as
the Wasserstein or Earth Mover’s distance. To define the Wasserstein metric,
we consider the space of probability measures with finite $p$-th moment
(Wasserstein space):

$\mathcal{P}_{p}(\mathcal{X})=\\{P\in\mathcal{P}(\mathcal{X}):\
\underset{\mathbf{X}\sim
P}{\mathbb{E}}[d^{p}(\mathbf{X},x^{\prime})]<\infty,\forall
x^{\prime}\in\mathcal{X}\\}$

For two probability measures $P,Q\in\mathcal{P}_{p}(\mathcal{X})$, the
$p$-Wasserstein distance is defined as:

$\displaystyle
W_{p}(P,Q)=\bigg{(}\underset{\pi\in\Gamma(P,Q)}{\text{inf}}\,\quad\underset{(\mathbf{X},\mathbf{X}^{\prime})\sim\pi}{\mathbb{E}}[d^{p}(\mathbf{X},\mathbf{X}^{\prime})]\bigg{)}^{\frac{1}{p}}.$

This metric is positive-definite, finite, symmetric, and adheres to the
triangle inequality (Villani, 2009, $\S$ 6).

######  Definition 10 (Truncated Probability Measure).

Let $P\in\mathcal{P}(\mathbb{R}^{d})$ be a probability measure and
$A\subseteq\mathbb{R}^{d}$ be Borel subset of $\mathbb{R}^{d}$ such that
$P(A)>0$, the truncated probability measure
$P_{A}\in\mathcal{P}(\mathbb{R}^{d})$ is defined as:

$P_{A}(B)=\frac{P(A\cap B)}{P(A)}$

for all $B\in\mathcal{B}(\mathbb{R}^{d})$.

######  Definition 11 (lower semi-continuous).

A function $f:\mathbb{R}^{d}\to\mathbb{R}\cup\\{+\infty\\}$ is said to be
lower semi-continuous (l.s.c) at a point $x_{0}\in rd$ if for every
$\epsilon>0$ there exists a neighborhood $U$ of $x_{0}$ such that for all
$x\in U$, $f(x)>f(x_{0})-\epsilon$. In economic terms, this concept implies
that small perturbations in the input do not lead to a substantial decrease in
the function value, signifying a form of stability or predictability in
economic models, such as cost functions in production processes.

######  Definition 12 (Closure).

Given a set $A$ in a topological space $X$, the _closure_ of $A$, denoted by
$\overline{A}$, is the smallest closed set in $X$ that contains $A$.
Equivalently, it includes all the points of $A$ along with all its limit
points (i.e., points that can be approached arbitrarily closely by points in
$A$).

######  Definition 13 (Continuous Measure).

A measure $\mu$ on a measurable space $(X,\mathcal{F})$ is said to be
_continuous_ if for every $A\in\mathcal{F}$, $\mu(A)=0$ whenever $A$ is a set
of a single point. In other words, $\mu(\\{x\\})=0$ for every $x\in X$.

######  Definition 14 (c-Conjugate).

Given a cost function $c:\mathcal{X}\times\mathcal{X}^{\prime}\to\mathbb{R}$,
where $\mathcal{X}$ and $\mathcal{X}^{\prime}$ are two metric spaces, the
$c$-conjugate of a function $\varphi:\mathcal{X}\to\mathbb{R}$ is a function
$\varphi^{c}:\mathcal{X}^{\prime}\to\mathbb{R}$ defined as follows:

$\varphi^{c}(x^{\prime})=\inf_{x\in\mathcal{X}}\\{c(x,x^{\prime})-\varphi(x)\\}$

This transformation is crucial in optimal transport theory, as it relates the
potential functions on the spaces $\mathcal{X}$ and $\mathcal{X}^{\prime}$
concerning the cost function $c$. The $c$-conjugate is used to establish
duality theorems and to find solutions to the Monge-Kantorovich problem.

######  Definition 15 ($\varphi$-divergence).

The $\varphi$-divergence between two probability distributions $P$ and $Q$
over the same probability space, for a convex function $\varphi$, is defined
as

$D_{\varphi}(P||Q)=\int\varphi\left(\frac{dP}{dQ}\right)dQ$

where $\frac{dP}{dQ}$ is the Radon-Nikodym derivative of $P$ with respect to
$Q$. The divergence measures the difference between the two distributions,
with different choices of $\varphi$ leading to different divergence measures.
Common examples include the Kullback-Leibler divergence for $\varphi(x)=x\log
x$, the Total Variation distance for $\varphi(x)=\frac{1}{2}|x-1|$, and the
squared Hellinger distance for $\varphi(x)=(\sqrt{x}-1)^{2}$.

######  Definition 16 (Length Space).

A length space is a metric space $(\mathcal{X},d)$ such that for every pair of
points $x,x^{\prime}\in\mathcal{X}$, the distance $d(x,x^{\prime})$ between
them is equal to the infimum of the lengths of continuous curves in $X$
connecting $x$ and $x^{\prime}$. Formally, this can be written as:

$d(x,x^{\prime})=\inf\\{L(\gamma)\mid\gamma\text{ is a curve connecting
}x\text{ and }x^{\prime}\\}$

where $L(\gamma)$ denotes the length of the curve $\gamma$, and a curve
$\gamma:[a,b]\to\mathcal{X}$ is a continuous map from a closed interval
$[a,b]$ in the real line $\mathbb{R}$ to the space $\mathcal{X}$.

######  Definition 17 (Locally Compact Space).

A topological space $\mathcal{X}$ is called locally compact if, for every
point $x\in\mathcal{X}$, there exists a neighborhood $U$ of $x$ such that the
closure of $U$, denoted as $\overline{U}$, is compact in $\mathcal{X}$. In
other words, for each $x\in\mathcal{X}$, there is an open set $U$ containing
$x$ such that the smallest closed set containing $U$ is compact.

######  Definition 18 (constant speed geodesic).

A constant speed geodesic in a metric space $(\mathcal{X},d)$ is a curve
$\gamma:[0,1]\to\mathcal{X}$ that satisfies the following properties:

  1. 1.

Geodesic Property: For any $t_{1},t_{2}\in[0,1]$, the distance
$d(\gamma(t_{1}),\gamma(t_{2}))$ is proportional to $|t_{2}-t_{1}|$. This
means $\gamma$ is the shortest path between its endpoints.

  2. 2.

Constant Speed: There exists a constant $c\geq 0$ such that for all
$t_{1},t_{2}\in[0,1]$, we have
$d(\gamma(t_{1}),\gamma(t_{2}))=c\cdot|t_{2}-t_{1}|$. This implies that the
curve $\gamma$ traverses distance at a uniform rate.

In simple terms, a constant speed geodesic is the shortest path in a metric
space that is traveled at a uniform speed.

######  Definition 19 (n-dimensional Hausdorff Measure).

Let $(X,d)$ be a metric space. The n-dimensional Hausdorff measure
$\mathcal{H}^{n}$ of a subset $A\subseteq X$ is defined as follows:

$\mathcal{H}^{n}(A)=\lim_{\delta\to
0}\inf\left\\{\sum_{i=1}^{\infty}\left(\operatorname{diam}(U_{i})\right)^{n}:\\{U_{i}\\}\text{
is a }\delta\text{-cover of }A\right\\}$

where a $\delta\text{-cover}$ of $A$ is a countable collection of sets
$\\{U_{i}\\}$ with $\operatorname{diam}(U_{i})<\delta$ such that
$A\subseteq\bigcup_{i}U_{i}$, and $\operatorname{diam}(U_{i})$ is the diameter
of the set $U_{i}$. For $n\in\mathbb{N}$, $\mathcal{H}^{n}$ generalizes the
notion of n-dimensional volume, with $\mathcal{H}^{1}$ representing length,
$\mathcal{H}^{2}$ area, and $\mathcal{H}^{3}$ volume.

######  Definition 20 ($\delta$-Confidence Positive Region).

$\delta$-confidence positive region is denoted by $L^{\textbf{+}}_{\delta}$
consists of smallest buffer of $L$ in $\mathcal{X}^{\textbf{+}}$ i.e.
$B=\\{x\in\mathcal{X}^{\textbf{+}}:\text{dist}(x,L)\leq r\\}$ such that
$P(B)\geq\epsilon$. where
$\text{dist}(x,L):=\inf\\{c(x,x^{\prime}):x^{\prime}\in L\\}$.

###  A.2 General Facts about Optimal Transport

######  Theorem 4 (Fundamental theorem of optimal transport).

Assume that $c:\mathcal{X}\times\mathcal{X}^{\prime}\rightarrow\mathbb{R}$ is
continuous and bounded from below and let
$\mu\in\mathcal{P}(\mathcal{X}),\nu\in\mathcal{P}(\mathcal{X}^{\prime})$ be
such that $c(x,x^{\prime})\leq a(x)+b(x^{\prime})$, for some $a\in
L^{1}(\mu),b\in L^{1}(\nu)$. Also, let $\pi\in\Gamma(\mu,\nu)$. Then the
following three are equivalent:

  * •

the plan $\pi$ is optimal,

  * •

the set $supp(\pi)$ is c-cyclically monotone,

  * •

there exists a c-concave function $\varphi$ such that $max\\{\varphi,0\\}\in
L_{1}(\mu)$ and $supp(\pi)\subset\partial^{c+}\varphi$

  * •

Then the minimum of the Kantorovich problem is equal to the supremum of the
dual problem 4. Furthermore, the supremum of the dual problem is attained, and
the maximizing couple $(\varphi,\psi)$ is of the form $(\varphi,\varphi^{c+})$
for some c-concave function $\varphi$.

###  A.3 Motivating Example

Assume individuals live in a finite set of states $\mathcal{X}$. A binary
classifier $h:\mathcal{X}\rightarrow\\{-1,1\\}$ divides these states into
positive states $\mathcal{X}^{+}$ and negative states $\mathcal{X}^{-}$. By
paying a cost of $c(x,x^{\prime})$, one can move from state $x$ to
$x^{\prime}$. In the original formulation of CE, the cost function is assumed
observable, and individuals are willing to pay any cost to get into
$\mathcal{X}^{+}$ from $\mathcal{X}^{-}$, so the question is about which
positive state costs less. This formulation simply treats all states in
$\mathcal{X}^{+}$ as equally beneficial, neglects the temporal dimension of
improvement, and is unaware of the incentives underlying individual progress
even in the absence of the decision made by $h$.

In a more realistic formulation, we assume state $x$ provides a utility of
$u(x)$ to an individual at $x$ at each time step. Unlike the cost function,
these utilities may not be directly observable to us. The individuals live for
an infinite horizon however they discount future additions to the utility by a
factor of $\gamma<1$ per time step (or effectively live for $1/(1-\gamma)$
steps). We may assume individuals are not capable of incurring a cost more
than$B$ at every step. So, to move from $x$ to $x^{\prime}$, since there is no
intermediate state in our finite state model, they have to wait for roughly
$c(x,x^{\prime})/B$ time steps at $x$. Without this movement, the effective
utility of an individual staying at $x$ forever is $u(x)/(1-\gamma)$. The
movement from $x$ to $x^{\prime}$ comes with an additional utility of

$u(x^{\prime})\frac{\gamma^{c(x,x^{\prime})/B}}{1-\gamma}+u(x)\frac{1-\gamma^{c(x,x^{\prime})/B}}{1-\gamma}-\frac{u(x)}{1-\gamma}.$

Define effective cost as
$\tilde{c}(x,x^{\prime})=1-\gamma^{c(x,x^{\prime})/B}$. Then an individual
maximizing her utility from movement should solve

$\underset{x^{\prime}\in\mathcal{X}}{\text{arg\,max}}\,u_{0}(x,x^{\prime})\coloneqq(1-\tilde{c}(x,x^{\prime}))(u(x^{\prime})-u(x)).$
(12)

Now in the presence of a decision made by $h$ (e.g., a loan), we assume all
the positive states come with an additional utility of $l$. We assume
individuals are not aware of $h$ and therefore in the absence of CE the best
action would come from maximizing $u_{0}$. This solution, although aligned
with the general interest of the individual, might miss the opportunity that
$h$ has generated. Optimally, CE should solve

$\underset{x^{\prime}\in\mathcal{X}^{+}}{\text{arg\,max}}\,u_{l}(x,x^{\prime})\coloneqq(1-\tilde{c}(x,x^{\prime}))(u(x^{\prime})+l-u(x))=u_{0}(x,x^{\prime})+l(1-\tilde{c}(x,x^{\prime})).$
(13)

However, the standard method of CE only partially solves this problem by
maximizing $(1-\tilde{c}(x,x^{\prime}))$, or equivalently neglecting
$u(x^{\prime})-u(x)$. We next discuss how the current actions of individuals
may guide CE to solve Eq. 12 instead of this simplification.

#### CE to Partially Rational Individuals.

In fact, although utilities may not be directly observable, assuming a level
of rationality for individuals, their actions partially reveal their perceived
utilities. Formally, let $\pi_{\beta}(x^{\prime}|x)$ be the distribution of
the current actions taken by individuals at state $x$. Recall that in our
simplified model we assumed the knowledge of $f$ is not publicly available and
individuals aim at maximizing $u_{0}$. We assume a $\lambda$-rational
individual chooses the action from

$\pi_{\beta}(x^{\prime}|x)=\exp(\lambda\,u_{0}(x,x^{\prime}))/z,$

where $z$ is a normalizing factor. A fully rational individual corresponds to
$\lambda\rightarrow\infty$. Using this model, we can write $u_{0}$ as
$\log(\pi_{\beta})/\lambda$ up to an additive constant. This allows us to
write the optimization of Eq. 13 as

$\underset{x^{\prime}\in\mathcal{X}^{+}}{\text{arg\,min}}\,\tilde{c}(x,x^{\prime})-\frac{1}{\lambda
l}\log(\pi_{\beta}(x^{\prime}|x)).$ (14)

Unfortunately, $\pi_{\beta}$ may not be directly observable and requires
longitudinal studies. Alternatively, under the stationary assumption, we can
approximate $\pi_{\beta}$ with a coupling that moves the mass from
$\mathcal{X}^{-}$ to $\mathcal{X}^{+}$. The Eq. 14 turns into an unbalanced
Sinkhorn’s algorithm ($\S$ B.2) with appropriately chosen entropy
regularization.

##  Appendix B Computational Experiments Supplementary Materials

###  B.1 Sinkhorn’s algorithm

The Sinkhorn algorithm (Peyré et al., 2017, $\S$ 4) approximates optimal
transport plans between distributions by using entropy regularization, which
speeds up computation and stabilizes solutions in large datasets. This
iterative method adjusts scaling vectors to minimize regularized transport
costs, beneficial for machine learning and data analysis in optimal transport
contexts. The computational complexity per Sinkhorn iteration is $O(n^{2})$,
with total complexity dependent on the number of iterations for convergence.

For two discrete distributions $\mathbf{P}$ and $\mathbf{Q}$, with $N$ and $M$
points respectively, and a cost matrix $\mathbf{C}$, the algorithm solves the
entropic regularized optimal transport problem. This problem is formulated as
minimizing:

$\min_{\mathbf{T}\geq 0}\langle\mathbf{T},\mathbf{C}\rangle_{F}-\epsilon\cdot
H(\mathbf{T})\quad\text{subject
to}\quad\mathbf{T}\mathbf{1}_{M}=\mathbf{P},\quad\mathbf{T}^{\top}\mathbf{1}_{N}=\mathbf{Q}$

where $\langle\cdot,\cdot\rangle_{F}$ denotes the Frobenius dot product,
$\epsilon$ is a small regularization parameter, $H(\mathbf{T})$ is the entropy
of the transport plan $\mathbf{T}$, and $\mathbf{1}_{M}$, $\mathbf{1}_{N}$ are
vectors of ones of appropriate sizes.

Algorithm 1 Sinkhorn Method for Optimal Plan

Input: probability measures $P=(p_{i})_{i}\in\mathbb{R}^{n}$ and
$Q=(q_{j})_{j}\in\mathbb{R}^{m}$, cost matrix
$C=(c_{ij})_{ij}\in\mathbb{R}^{n\times m}$, regularization parameter
$\epsilon$ and number of iterations $N$;

Output: approximated optimal transport matrix $\pi$;

Initialize: $u_{0}=\mathbf{1}_{n}$, $v_{0}=\mathbf{1}_{m}$;

Compute: $K=e^{-\epsilon C}$

for $n=0$ to $N-1$ do

Update $u_{n+1}=\frac{P}{Kv_{n}}$

Update $v_{n+1}=\frac{Q}{K^{\top}u_{n}}$

end for

Return: $\mathbf{T}=\text{diag}(u_{N})K\text{diag}(v_{N})$;

###  B.2 Unbalanced Sinkhorn’s algorithm

The unbalanced Sinkhorn algorithm (Peyré et al., 2017, $\S$ 10) is an
extension of the classic Sinkhorn algorithm, adapted for solving optimal
transport problems where the mass of the distributions does not necessarily
match. The Unbalanced Sinkhorn algorithm modifies this problem to allow for
differences in mass between $\mathbf{P}$ and $\mathbf{Q}$. The constraints are
relaxed using so-called Kullback-Leibler (KL) divergence terms, leading to the
unbalanced optimal transport problem:

$\min_{\mathbf{T}\geq 0}\langle\mathbf{T},\mathbf{C}\rangle_{F}-\epsilon\cdot
H(\mathbf{T})+\lambda_{1}\cdot\text{KL}(\mathbf{T}\mathbf{1}_{M}\|\mathbf{P})+\lambda_{2}\cdot\text{KL}(\mathbf{T}^{\top}\mathbf{1}_{N}\|\mathbf{Q})$

Here, $\lambda_{1}$ and $\lambda_{2}$ are regularization parameters for the
marginal constraints, and the KL divergence terms $\text{KL}(\cdot\|\cdot)$
measure the discrepancy between the marginals of the transport plan
$\mathbf{T}$ and the given distributions

In empirical settings, the algorithm deals with discrete distributions often
derived from data samples. This involves computing a transport plan between
empirical distributions, which are represented as sums of Dirac masses. The
empirical part of the algorithm refers to its application to empirical
distributions, i.e., distributions represented by samples (data points), which
is a common scenario in practical applications.

Algorithm 2 Unbalanced Sinkhorn Optimal Transport

Input: probability measures $P=(p_{i})_{i}\in\mathbb{R}^{n}$ and
$Q=(q_{j})_{j}\in\mathbb{R}^{m}$, cost matrix
$C=(c_{ij})_{ij}\in\mathbb{R}^{n\times m}$, regularization parameter
$\epsilon$, $\lambda_{1}$ and $\lambda_{2}$ regularization parameters and
number of iterations $N$;

Output: approximated optimal transport matrix $\pi$;

Initialize: $u_{0}=\mathbf{1}_{n}$, $v_{0}=\mathbf{1}_{m}$;

Compute: $K=e^{-\epsilon C}$

for $n=0$ to $N-1$ do

Update
$u_{n+1}=\big{(}\frac{P}{Kv_{n}}\big{)}^{\frac{\lambda_{1}}{\epsilon+\lambda_{1}}}$

Update
$v_{n+1}=\big{(}\frac{Q}{K^{\top}u_{n}})^{\frac{\lambda_{2}}{\epsilon+\lambda_{2}}}$

end for

Return: $\mathbf{T}=\text{diag}(u_{N})K\text{diag}(v_{N})$;

###  B.3 The back-and-forth method

The Back-and-Forth Jacobs and Léger (2020) method offers a robust solution for
computing optimal transport maps with strictly convex costs, including p-power
costs, for probability densities $P$ and $Q$ on an $n$-point grid. This
method, characterized by its computational efficiency, requires $O(n)$ storage
and $O(n\log(n))$ computation per iteration. The iteration count needed to
achieve $\epsilon$ accuracy is proportional to
$O(\max(\|P\|_{\infty},\|Q\|_{\infty})\log(\frac{1}{\epsilon}))$, showcasing
the method’s effectiveness in both storage and computational resource
optimization.

In the back-and-forth method $\Omega$, is consider as a convex and compact
subset of $\mathbb{R}^{d}$, and focuses on a specific cost function
$c\colon\Omega\times\Omega\to\mathbb{R}$ defined as $c(x,y)=h(y-x)$. Here,
$h\colon\mathbb{R}^{d}\to\mathbb{R}$ is a strictly convex and even function.
The dual Kantorovich problem is considered in the two following equivalent
forms:

$I(\psi)=\int\psi\,d\mu+\int\psi^{c}\,d\nu\quad\text{or}\quad
J(\phi)=\int\varphi\,d\nu+\int\varphi^{c}\,d\mu,$

where $\varphi^{c},\psi^{c}$ are $c-$transformation of $\varphi,\psi$. The
gradient of $J(\phi)$ in the space of functions from $\Omega$ to $\mathbb{R}$
can be written as ( Jacobs and Léger, 2020, Lemma 3):

$\nabla J(\phi)=(-\Delta)^{-1}\big{(}\nu-T_{\varphi\,*}P\big{)},$

where the $\Delta$ is Laplacian operator and

$T_{\varphi}(x)=x-(\nabla h)^{-1}(\nabla\varphi^{c}(x))$

is the explicit solution of c-transformation (Jacobs and Léger, 2020, Lemma
1). We now introduce the gradient descent back-and-forth method.

Algorithm 3 Back-and-Forth Method for Optimal Map

Input: probability measures $P$ and $Q$, cost function $c$ and number of
iterations $N$;

Output: approximated Kantorovich potential functions;

Initialize: set $\varphi_{0}=0,\psi_{0}=0$

for $n=0$ to $N-1$ do

$\varphi_{n+\frac{1}{2}}=\varphi_{n}+\sigma\nabla J(\varphi_{n})$,

$\psi_{n+\frac{1}{2}}=(\varphi_{n+\frac{1}{2}})^{c}$,

$\psi_{n+1}=\psi_{n+\frac{1}{2}}+\sigma\nabla I(\psi_{n+\frac{1}{2}})$,

$\varphi_{n+1}=(\psi_{n+1})^{c}$.

end for

Return: $\varphi_{N},\psi_{N}$;

Theorem 2 establishes that if $\varphi$ is the solution to the dual
Kantorovich problem, then $T_{\varphi}$ constitutes the optimal map.

###  B.4 Supplementary Simulation

Figure 9:  (left) For the Gaussian mixture with means $\mu_{1}=[0.3,0.3]$ and
$\mu_{2}=[0.7,0.7]$, and a covariance matrix of $\text{diag}(0.2,0.2)$ and
classifier with $f(x)=x_{1}+x_{2}-1$: (a) displays the probability measure $P$
alongside the $\delta$-confidence region, (b) illustrates the negatively
induced probability measure $P_{-}$, and (c) shows the probability measure
induced by the $\delta$-confidence region $L^{+}_{\delta}$. (Right) Blue
curve: represents the percentage increase in cost relative to standard CE as
$\lambda_{2}$ varies from 0 to 1. Red curve: Depicts the logarithmic value of
the KL-divergence between the UOP solution and the target distribution.

Figure 10: (top) This figure demonstrates controlling additional costs in CCE
using unbalanced OT by adjusting the parameter $\lambda_{2}$. When
$\lambda_{2}=0$ (d), the solution resembles standard CE (b) while setting
$\lambda_{2}=1$ (f) aligns the estimation of CCE with the Sinkhorn method (c).
(down) Utilizing the back-and-forth method, the CCE map is estimated, and
through constructing the displacement interpolation, the optimal flow that
pushes forward the negative area into the $\delta$-confidence region is
visualized.

Figure 11:  (left) For the Gaussian mixture with means $\mu_{1}=[0.3,0.3]$ and
$\mu_{2}=[0.7,0.7]$, and a covariance matrix of $\text{diag}(0.2,0.2)$ and
classifier with implicit classifier formula
$f(x)=(x_{1}-0.5)^{2}+(x_{2}-0.5)^{2}-\frac{1}{9}$. (a) displays the
probability measure $P$ alongside the $\delta$-confidence region, (b)
illustrates the negatively induced probability measure $P_{-}$, and (c) shows
the probability measure induced by the $\delta$-confidence region
$L^{+}_{\delta}$." (Right) Blue curve: represents the percentage increase in
cost relative to standard CE as $\lambda_{2}$ varies from 0 to 1. Red curve:
Depicts the logarithmic value of the KL-divergence between the UOP solution
and the target distribution.

Figure 12: (top) This figure demonstrates controlling additional costs in CCE
using unbalanced OT by adjusting the parameter $\lambda_{2}$. When
$\lambda_{2}=0$ (d), the solution resembles standard CE (b) while setting
$\lambda_{2}=1$ (f) aligns the estimation of CCE with the Sinkhorn method (c).
(down) Utilizing the back-and-forth method, the CCE map is estimated, and
through constructing the displacement interpolation, the optimal flow that
pushes forward the negative area into the $\delta$-confidence region is
visualized.