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

HAL
open science

# Optimal Transport on Categorical Data for Counterfactuals using Compositional Data and Dirichlet Transport

Agathe Fernandes Machado, Arthur Charpentier, Ewen Gallic

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

## To cite this version:

Agathe Fernandes Machado, Arthur Charpentier, Ewen Gallic. Optimal Transport on Categorical Data for Counterfactuals using Compositional Data and Dirichlet Transport. 2025. hal-05085147

HAL Id: hal-05085147

https://hal.science/hal-05085147v1

Preprint submitted on 26 May 2025

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Optimal Transport on Categorical Data
for Counterfactuals using Compositional Data
and Dirichlet Transport

Agathe Fernandes Machado<sup>1</sup>, Arthur Charpentier<sup>1</sup>, and Ewen Gallic<sup>2,3</sup>

<sup>1</sup>Département de Mathématiques, Université du Québec à Montréal, Montréal, Québec, Canada

<sup>2</sup>CNRS - Université de Montréal CRM - CNRS

<sup>3</sup>Aix Marseille Univ, CNRS, AMSE, Marseille, France

May 22, 2025

Abstract

Recently, optimal transport-based approaches have gained attention for deriving coun-
terfactuals, e.g., to quantify algorithmic discrimination. However, in the general multi-
variate setting, these methods are often opaque and difficult to interpret. To address this,
alternative methodologies have been proposed, using causal graphs combined with itera-
tive quantile regressions (Plečko and Meinshausen, 2020) or sequential transport (Fernan-
des Machado et al., 2025) to examine fairness at the individual level, often referred to as
"counterfactual fairness." Despite these advancements, transporting categorical variables
remains a significant challenge in practical applications with real datasets. In this paper,
we propose a novel approach to address this issue. Our method involves (1) converting cat-
egorical variables into compositional data and (2) transporting these compositions within
the probabilistic simplex of ℝ<sup>d</sup>. We demonstrate the applicability and effectiveness of this
approach through an illustration on real-world data, and discuss limitations.

# 1 Introduction

## 1.1 Counterfactuals

Counterfactual analysis, the third level in Pearl (2009)'s causal hierarchy, is widely used in ma-
chine learning, policy evaluation, economics and causal inference. It involves reasoning about
"what could have happened" under alternative scenarios, providing insights into causality and
decision-making effectiveness. An example could be the concept of counterfactual fairness, as
introduced by Kusner et al. (2017), that ensures fairness by evaluating how decisions would

Agathe Fernandes Machado acknowledges that the project leading to this publication has received funding from OBVIA. Arthur Charpentier acknowledges funding from the SCOR Foundation for Science and the National Sciences and Engineering Research Council (NSERC) for funding (RGPIN-2019-07077). Ewen Gallic acknowledges funding from the French government under the "France 2030" investment plan managed by the French National Research Agency (reference: ANR-17-EURE-0020) and from Excellence Initiative of Aix-Marseille University - A*MIDEX.

Replication codes and companion e-book: https://github.com/fer-agathe/transport-simplex

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

Corresponding author: fernandes_machado.agathe@courrier.uqam.ca
change under alternative, counterfactual conditions. Counterfactual fairness focuses on mitigating bias by ensuring that sensitive attributes, such as race, gender, or socioeconomic status, do not unfairly influence outcomes.

In the counterfactual problem, we consider data $\{(s_i, x_i), i = 1, \dots, n\}$, where $s$ denotes a binary “treatment” (taking values in $\{0, 1\}$). With generic notations, the counterfactual version of $(0, x)$ can be constructed as $(1, T^*(x))$, where $T^*$ is the optimal transport (OT) mapping from $X|S=0$ to $X|S=1$, as discussed in Black et al. (2020), Charpentier et al. (2023) and De Lara et al. (2024). Unfortunately, this multivariate mapping is usually both complicated to estimate, and hard to interpret. If $x$ is univariate, it is simply a quantile interpretation: if $x$ is associated to rank probability $u$ within group $s=0$, then its counterfactual version should be associated with the same rank probability in group $s=1$ (mathematically, $T^* = F_1^{-1} \circ F_0$, where $F_j : \mathbb{R} \to [0, 1], j = \{0, 1\}$ denotes the cumulative distribution in group $j$, and $F_j^{-1}$ is the generalized inverse, i.e., the quantile function). In higher dimensions, one could consider multivariate quantiles, as in Hallin et al. (2021) or Hallin and Konen (2024), but the heuristics is still hard to interpret. While OT-based counterfactual methods have been proposed to assess counterfactual fairness Black et al. (2020); De Lara et al. (2024), an alternative approach introduced by Plečko and Meinshausen (2020) is grounded in causal graphs (DAGs). In this framework, the outcome $y$ depends on variables $(s, x)$, where the sensitive attribute $s$ “is a source” (a vertex without parents) and $y$ is a “sink” (a vertex without outgoing edges). Recently, Fernandes Machado et al. (2025) unified these approaches by introducing sequential transport aligned with the “topological ordering” of a DAG.

For example, to test whether a predictor $\hat{m}(x)$ is gender-neutral; let the sensitive attribute $s$ be gender (binary genders for simplicity); compare its output on a woman's features $x$ with that on her *mutatis mutandis* male counterpart. Unlike a *ceteris paribus* change, which flips $s$ while holding all other features fixed, a *mutatis mutandis* intervention also adjusts any $x_j$ causally influenced by $s$. Thus, if $x_1$ is height, the counterfactual of a 5'4" woman would not be a 5'4" man but, say, a 5'10" man, via an OT map. While OT handles continuous attributes naturally, categorical features (e.g. occupation or neighbourhood) lack a canonical distance. As a result, generating counterfactuals (e.g. the male counterpart of a female *nurse*, or where a Black *resident of X* would live if they were White) becomes particularly challenging.

## 1.2 The Case of Categorical Variables

For absolutely continuous variables, the approaches of Plečko and Meinshausen (2020); Plečko et al. (2024) on the one hand (based on quantile regressions) and Black et al. (2020); Charpentier et al. (2023); De Lara et al. (2024); Fernandes Machado et al. (2025) (based on OT) are quite similar.

If Plečko and Meinshausen (2020) considered quantile regressions for absolutely continuous variables, the case of ordered categorical variables is considered (at least with some sort of meaningful ordering) in the section related to “Practical aspects and extensions.” Discrete optimal transport between two marginal multinomial distributions is considered, but as discussed, it suffers multiple limitations. Here, we will consider an alternative approach, based on the idea of transforming categorical variables into continuous ones, coined “compositional variables” in Chayes (1971), and then, using “Dirichlet optimal transport,” on those compositions.

While motivated by counterfactual fairness, the primary aim of this study is to present the core of a method for deriving counterfactuals for categorical data, applicable to any context requiring counterfactual analysis. Here, for simplicity, we have set aside considerations related to the assumption of a known Structural Causal Model (SCM).<fnref n="1" />

<footnote n="1">¹Details on how the method can be integrated within an SCM are discussed in Appendix C.</footnote>
## 1.3 Agenda

After recalling notations on OT in Section 2, we discuss how to transform categorical variables with $d$ categories into variables taking values in the simplex $S_d$ in $\mathbb{R}^d$, i.e., compositional variables, in Section 3. In Section 4, we review the topological and geometrical properties of the probability simplex $S_d \subset \mathbb{R}^d$. In Section 5, we introduce the first methodology, which transports distributions within $S_d$ via Gaussian OT. This approach relies on an alternative representation of probability vectors in the Euclidean space $\mathbb{R}^{d-1}$ and assumes approximate normality in the transformed space. In Section 6, we present a second methodology, which operates directly on $S_d$ using a tailored cost function instead of the standard quadratic cost. Theoretical aspects of this “Dirichlet transport” framework are discussed in Section 6.1, while empirical strategies for matching categorical observations are developed in Section 6.2. Section 7 provides two empirical illustrations using the German Credit and Adult datasets.

Our main contributions can be summarized as follows:

* We propose a novel method to handle categorical variables in counterfactual modeling by using optimal transport directly on the simplex. This approach transforms categorical variables into compositional data, enabling the use of probabilistic representations that preserve the geometric structure of the simplex.

* By integrating optimal transport techniques on this domain, the method ensures consistency with the properties of compositional data and offers a robust framework for counterfactual analysis in real-world scenarios.

* Our approach does not require imposing an arbitrary order on the labels of categorical variables.

## 2 Optimal Transport

Given two metric spaces $\mathcal{X}_0$ and $\mathcal{X}_1$, consider a measurable map $T: \mathcal{X}_0 \to \mathcal{X}_1$ and a measure $\mu_0$ on $\mathcal{X}_0$. The push-forward of $\mu_0$ by $T$ is the measure $\mu_1 = T_{\#}\mu_0$ on $\mathcal{X}_1$ defined by $T_{\#}\mu_0(B) = \mu_0(T^{-1}(B))$, $\forall B \subset \mathcal{X}_1$. For all measurable and bounded $\varphi: \mathcal{X}_1 \to \mathbb{R}$,

$$ \int_{\mathcal{X}_1} \varphi(\mathbf{x}_1) T_{\#} \mu_0(\mathrm{d}\mathbf{x}_1) = \int_{\mathcal{X}_0} \varphi(T(\mathbf{x}_0)) \mu_0(\mathrm{d}\mathbf{x}_0). $$

For our applications, if we consider measures $\mathcal{X}_0 = \mathcal{X}_1$ as a compact subset of $\mathbb{R}^d$, then there exists $T$ such that $\mu_1 = T_{\#}\mu_0$, when $\mu_0$ and $\mu_1$ are two measures, and $\mu_0$ is atomless, as shown in Villani (2003) and Santambrogio (2015). Out of those mappings from $\mu_0$ to $\mu_1$, we can be interested in “optimal” mappings, satisfying Monge problem, from Monge (1781), i.e., solutions of

$$ \inf_{T_{\#}\mu_0=\mu_1} \int_{\mathcal{X}_0} c(\mathbf{x}_0, T(\mathbf{x}_0)) \mu_0(\mathrm{d}\mathbf{x}_0), \qquad (1) $$

for some positive ground cost function $c: \mathcal{X}_0 \times \mathcal{X}_1 \to \mathbb{R}_+$. In general settings, however, such a deterministic mapping $T$ between probability distributions may not exist (in particular if $\mu_0$ and $\mu_1$ are not absolutely continuous, with respect to Lebesgue measure). This limitation motivates the Kantorovich relaxation of Monge's problem, as considered in Kantorovich (1942),

$$ \inf_{\pi \in \Pi(\mu_0, \mu_1)} \int_{\mathcal{X}_0 \times \mathcal{X}_1} c(\mathbf{x}_0, \mathbf{x}_1) \pi(\mathrm{d}\mathbf{x}_0, \mathrm{d}\mathbf{x}_1), \qquad (2) $$

with our cost function $c$, where $\Pi(\mu_0, \mu_1)$ is the set of all couplings of $\mu_0$ and $\mu_1$. This problem focuses on couplings rather than deterministic mappings. It always admits solutions referred to
**Algorithm 1** From categorical variables into compositions.

Input: training dataset $\mathcal{D} = \{(s_i, \mathbf{x}_i)\}$

Input: new observation $(s, \mathbf{x})$, with $\mathbf{x}_j$'s either in $\mathbb{R}$ or $[[d_j]]$

Output: $(s, \tilde{\mathbf{x}})$, with $\tilde{\mathbf{x}}_j$'s either in $\mathbb{R}$ or $S_{d_j}$

for $j \in \{1, \dots, k\}$ do

if $\mathbf{x}_j \in [[d_j]]$ then

estimate a MLR to predict categorical $\mathbf{x}_j$ using $\mathcal{D}$

get estimates $\hat{\beta}_2, \dots, \hat{\beta}_{d_j}$

$\tilde{\mathbf{x}}_j \leftarrow \mathcal{C}(1, e^{\mathbf{x}_j^\top \hat{\beta}_2}, \dots, e^{\mathbf{x}_j^\top \hat{\beta}_d})$

else

$\tilde{\mathbf{x}}_j \leftarrow \mathbf{x}_j$

end if

end for

as OT plans. Observe that $T^*$ is an “increasing mapping,” in the sense of being the gradient of a convex function, from Brenier (1991)). Finally, one should have in mind the cost function $c$ is related to the geometry of sets $\mathcal{X}$.

## 3 From Categorical to Compositional Data

Using the notations of the introduction, consider a dataset $\{s, \mathbf{x}\}$ where features $\mathbf{x}$ are either numerical (assumed to be “continuous”), or categorical. In the latter case, suppose that $\mathbf{x}_j$ takes values in $\{x_{j,1}, \dots, x_{j,d_j}\}$, or more conveniently, $[[d_j]] = \{1, \dots, d_j\}$, corresponding to the $d_j$ categories (as in the standard “One Hot” encoding).

The aim is to transform a categorical variable $x$, which takes values in $[[d]]$, into a numerical one in the simplex $S_d$. To achieve this, we suggest using a probabilistic classifier. This classifier is based on the other features in $\mathbf{x}$, denoted by $\mathcal{X}_{-x}$. Mathematically, we consider a mapping from $\mathcal{X}_{-x}$ to $S_d$ (and not to $[[d]]$ as in a standard multiclass classifier). The most natural model for this transformation is the Multinomial Logistic Regression (MLR), which is based on the “softmax” loss function. To normalize the output of the classifier into the simplex, we define the closure operator $\mathcal{C}: \mathbb{R}_+^d \to S_d$ as

$$ \mathcal{C}[x_1, x_2, \dots, x_d] = \left[ \frac{x_1}{\sum_{i=1}^d x_i}, \frac{x_2}{\sum_{i=1}^d x_i}, \dots, \frac{x_d}{\sum_{i=1}^d x_i} \right], $$

or shortly

$$ \mathcal{C}(\mathbf{x}) = \frac{\mathbf{x}}{\mathbf{x}^\top 1}, $$

where $1$ is a vector of ones in $\mathbb{R}^d$. Then, in the MLR model, the transformation $\hat{T}: \mathcal{X}_{-x} \to S_d$ is given by

$$ \hat{T}(\mathbf{x}) = \mathcal{C}(1, e^{\mathbf{x}^\top \hat{\beta}_2}, \dots, e^{\mathbf{x}^\top \hat{\beta}_d}) \in S_d, $$

where $\hat{\beta}_2, \dots, \hat{\beta}_d$ are the estimated coefficients for each category, and the first category is taken as the reference. This procedure is described in Algorithm 1.

As an illustration, consider the purpose variable from the German dataset. For simplicity, this variable has been reduced to three categories: C, E, O (representing cars, equipment, and other, respectively). More details on the dataset are provided in Section 7.1. The purpose variable is converted into a continuous variable using four models: (i) a GAM-MLR with splines
<table><caption>Table 1: Mappings from the purpose categorical variable *x* to the compositional one $\tilde{x}$, (in the german credit dataset), for the first four individuals of the dataset. The first two models are GAM-MLR (multinomial model with splines for continuous variables), then, a random forest, and a boosting algorithm.</caption><thead><tr><th colspan="4">GAM-MLR (1)</th><th colspan="4">random forest</th></tr><tr><th>x</th><th>x&#772;<sub>C</sub></th><th>x&#772;<sub>E</sub></th><th>x&#772;<sub>O</sub></th><th>x&#772;<sub>C</sub></th><th>x&#772;<sub>E</sub></th><th>x&#772;<sub>O</sub></th></tr></thead><tbody><tr><td>E</td><td>18.38%</td><td><strong>61.56%</strong></td><td>20.06%</td><td>23.68%</td><td><strong>46.32%</strong></td><td>30.00%</td></tr><tr><td>C</td><td>40.86%</td><td><strong>42.38%</strong></td><td>16.76%</td><td>34.68%</td><td><strong>36.42%</strong></td><td>28.90%</td></tr><tr><td>E</td><td>19.41%</td><td><strong>70.82%</strong></td><td>9.77%</td><td>16.87%</td><td><strong>76.51%</strong></td><td>6.63%</td></tr><tr><td>C</td><td><strong>47.04%</strong></td><td>26.83%</td><td>26.13%</td><td><strong>53.16%</strong></td><td>26.84%</td><td>20.00%</td></tr></tbody></table><table><thead><tr><th colspan="4">GAM-MLR (2)</th><th colspan="4">gradient boosting model</th></tr><tr><th>x</th><th>x&#772;<sub>C</sub></th><th>x&#772;<sub>E</sub></th><th>x&#772;<sub>O</sub></th><th>x&#772;<sub>C</sub></th><th>x&#772;<sub>E</sub></th><th>x&#772;<sub>O</sub></th></tr></thead><tbody><tr><td>E</td><td>9.22%</td><td><strong>75.92%</strong></td><td>14.86%</td><td>11.25%</td><td><strong>68.51%</strong></td><td>20.24%</td></tr><tr><td>C</td><td><strong>46.80%</strong></td><td>24.06%</td><td>29.14%</td><td><strong>61.14%</strong></td><td>13.10%</td><td>25.76%</td></tr><tr><td>E</td><td>11.23%</td><td><strong>79.07%</strong></td><td>9.71%</td><td>12.48%</td><td><strong>75.58%</strong></td><td>11.94%</td></tr><tr><td>C</td><td><strong>50.74%</strong></td><td>26.98%</td><td>22.28%</td><td><strong>51.12%</strong></td><td>25.17%</td><td>23.71%</td></tr></tbody></table>

for three continuous variables, (ii) a GAM-MLR incorporating these variables and seven categorical ones, (iii) a random forest, and (iv) a gradient boosting model. Table 1 presents the observed values in the first column for each model, along with the estimated scores for each category in the three remaining columns, corresponding to the transformed values $T^*(x)$.

Note that if we want to go back from compositions to categories, the standard approach is based on the majority (or argmax) rule.

In the rest of the paper, given a dataset {$s, x$}, all categorical variables are transformed into compositions, so that $\mathcal{X}$ is a product space of sets that are either $\mathbb{R}$ for numerical variables or $S_d$ (type) for compositions ($d$ will change according to the number of categories).

In fact, for privacy issues, a classical strategy is to consider aggregated data on small groups (usually on a geographic level, per block, or per zip code), even if there is an ecological fallacy issue (that occurs when conclusions about individual behaviour or characteristics are incorrectly drawn based on aggregate data for a group, see King et al. (2004)). Hence, using “compositional data” is quite natural in many cases, as unobserved categorical variables can often be represented as compositions predicted from observed variables serving as proxies. For example, in U.S. datasets, racial information about individuals may not always be available. However, the proportions of groups such as “White and European,” “Asian,” “Hispanic and Latino,” “Black or African American,” etc., within a neighbourhood may be observed instead (see, e.g., Cheng et al. (2010), Naeini et al. (2015) and Zadrozny and Elkan (2001) for more general discussions, or Imai et al. (2022) about the use of predicted probabilities when categories are not observed).

# 4 Topology and Geometry of the Simplex

The standard simplex of $\mathbb{R}^d$ is the regular polytope $S_d = \{x \in \mathbb{R}_+^d \mid x^\top 1 = 1\}$, but for convenience, consider the open version of that set,

$$S_d = \{x \in (0, 1)^d \mid x^\top 1 = 1\}.$$

Following Aitchison (1982), define the inner product

$$\langle \mathbf{x}, \mathbf{y} \rangle = \frac{1}{d} \sum_{i<j} \log \frac{x_i}{x_j} \log \frac{y_i}{y_j} \quad \forall \mathbf{x}, \mathbf{y} \in S_d, \qquad (3)$$
and the simplex becomes a metric vector space if we consider the associated “Aitchison distance,” as coined in Pawlowsky-Glahn and Egozcue (2001). Figure 1 shows $n = 61$ points in $S_3$. Each point $x$ can be seen as a probability vector over $\{A, B, C\}$, drawn either from a distribution $P_0$ for red points or $P_1$ for blue points.

<figure><img src="image_5.png" /><figcaption>Figure 1: $n = 61$ points in $S_3$, with a toy dataset.</figcaption></figure>

If we define the binary operator $\diamond$ on $S_d$,

$$ \mathbf{x} \diamond \mathbf{y} = \left[ \frac{x_1 y_1}{\sum_{i=1}^{d} x_i y_i}, \dots, \frac{x_d y_d}{\sum_{i=1}^{d} x_i y_i} \right], $$

then $(S_d, \diamond)$ is a commutative group, with identity element $d^{-1}1$, and the inverse of $\mathbf{x}$ is

$$ \mathbf{x}^{-1} = \left[ \frac{1/x_1}{\sum_{i=1}^{d} 1/x_i}, \dots, \frac{1/x_d}{\sum_{i=1}^{d} 1/x_i} \right] = C(1/\mathbf{x}). $$

# 5 Using an Alternative Representation of Simplex Data

A first strategy to define a transport mapping could be to use some isomorphism, $h: S_d \to \mathcal{E}$ and then define the inverse mapping $h^{-1}: \mathcal{E} \to S_d$, where $\mathcal{E}$ is some Euclidean space, classically $\mathbb{R}^{d-1}$, where the standard quadratic cost can be considered. This idea corresponds to the dual transport problem in Pal and Wong (2018).

## 5.1 Classical Transformations

The additive log ratio (alr) transform is an isomorphism where $\text{alr}: S_d \to \mathbb{R}^{d-1}$, given by

$$ \text{alr}(\mathbf{x}) = \left[ \log \frac{x_1}{x_d}, \dots, \log \frac{x_{d-1}}{x_d} \right]. $$
Its inverse is, for any $\mathbf{z} \in \mathbb{R}^{d-1}$,

$$ \operatorname{clr}^{-1}(\mathbf{z}) = \mathcal{C}(\exp(z_1), \dots, \exp(z_{d-1}), 1) = \mathcal{C}(\exp([\mathbf{z}, 0])) $$

Such a map, from $S_d$ to $\mathbb{R}^{d-1}$ is related to the so-called “exponential coordinate system” of the unit simplex, in Pal (2024). The center log ratio (clr) transform is both an isomorphism and an isometry where $\operatorname{clr}: S^d \to \mathbb{R}^d$,

$$ \operatorname{clr}(\mathbf{x}) = \left[ \log \frac{x_1}{\bar{\mathbf{x}}_g}, \dots, \log \frac{x_D}{\bar{\mathbf{x}}_g} \right], $$

where $\bar{\mathbf{x}}_g$ denotes the geometric mean of $\mathbf{x}$. Observe that the inverse of this function is the softmax function, i.e.,

$$ \operatorname{clr}^{-1}(\mathbf{z}) = \mathcal{C}(\exp(z_1), \dots, \exp(z_d)) = \mathcal{C}(\exp(\mathbf{z})), \mathbf{z} \in \mathbb{R}^d. $$

Finally, the isometric log ratio (ilr) transform, defined in Egozcue et al. (2003), is both an isomorphism and an isometry where $\operatorname{ilr}: S_d \to \mathbb{R}^{d-1}$,

$$ \operatorname{ilr}(\mathbf{x}) = [\langle \mathbf{x}, \vec{e}_1 \rangle, \dots, \langle \mathbf{x}, \vec{e}_{d-1} \rangle] $$

for some orthonormal base $\{\vec{e}_1, \dots, \vec{e}_{d-1}, \vec{e}_d\}$ of $\mathbb{R}^d$. One can consider some matrix $\mathbf{M}$, $d \times (d-1)$ such that $\mathbf{M}\mathbf{M}^\top = \mathbb{I}_{d-1}$ and $\mathbf{M}^\top\mathbf{M} = \mathbb{I}_d + \mathbf{1}_{d \times d}$. Then

$$ \operatorname{ilr}(\mathbf{x}) = \operatorname{clr}(\mathbf{x})\mathbf{M} = \log(\mathbf{x})\mathbf{M}, $$

and

$$ \operatorname{ilr}^{-1}(\mathbf{z}) = \mathcal{C}((\exp(\mathbf{z}\mathbf{M}^\top)), \mathbf{z} \in \mathbb{R}^{d-1}). $$

## 5.2 Gaussian Mapping in the Euclidean Representation

Given a random vector $\mathbf{X}$ in $S_d$, we say that $\mathbf{x}$ follows a “normal distribution on the simplex” if, for some isomorphism $h$, the vector of orthonormal coordinates, $\mathbf{Z} = h(\mathbf{X})$ follows a multivariate normal distribution on $\mathbb{R}^{d-1}$. If we suppose that both $\mathbf{X}_0$ and $\mathbf{X}_1$, taking values in $S_d$, follow “normal distributions on the simplex,” then we can use standard Gaussian optimal transport, between $\mathbf{Z}_0$ and $\mathbf{Z}_1$. For convenience, suppose that the same isomorphism is used for both distributions (but that assumption can easily be relaxed). Hence, if $\mathbf{Z}_0 \sim \mathcal{N}(\boldsymbol{\mu}_0, \boldsymbol{\Sigma}_0)$ and $\mathbf{Z}_1 \sim \mathcal{N}(\boldsymbol{\mu}_1, \boldsymbol{\Sigma}_1)$, the optimal mapping is linear,

$$ \mathbf{z}_1 = T^*(\mathbf{z}_0) = \boldsymbol{\mu}_1 + \mathbf{A}(\mathbf{z}_0 - \boldsymbol{\mu}_0), \quad (4) $$

where $\mathbf{A}$ is a symmetric positive matrix that satisfies $\mathbf{A}\boldsymbol{\Sigma}_0\mathbf{A} = \boldsymbol{\Sigma}_1$, which has a unique solution given by $\mathbf{A} = \boldsymbol{\Sigma}_0^{-1/2}(\boldsymbol{\Sigma}_0^{1/2}\boldsymbol{\Sigma}_1\boldsymbol{\Sigma}_0^{1/2})^{1/2}\boldsymbol{\Sigma}_0^{-1/2}$, where $\mathbf{M}^{1/2}$ is the square root of the square (symmetric) positive matrix $\mathbf{M}$ based on the Schur decomposition ($\mathbf{M}^{1/2}$ is a positive symmetric matrix), as described in Higham (2008). Interestingly, it is possible to derive McCann’s displacement interpolation, from McCann (1997), to have some sort of continuous mapping $T_t^*$ such that $T_t^* = T^*$ and $T_0 = \mathbb{I}_d$, and so that $\mathbf{Z}_t = T_t^*(\mathbf{Z}_0)$ has distribution $\mathcal{N}(\boldsymbol{\mu}_t, \boldsymbol{\Sigma}_t)$ where $\boldsymbol{\mu}_t = (1-t)\boldsymbol{\mu}_0 + t\boldsymbol{\mu}_1$ and

$$ \boldsymbol{\Sigma}_t = \boldsymbol{\Sigma}_0^{-1/2} \left( (1-t)\boldsymbol{\Sigma}_0 + t \left( \boldsymbol{\Sigma}_0^{1/2} \boldsymbol{\Sigma}_1 \boldsymbol{\Sigma}_0^{1/2} \right)^{1/2} \right)^2 \boldsymbol{\Sigma}_0^{-1/2}. $$
<figure><img src="image_6.png" /><figcaption>Figure 2: Counterfactuals using the ilr transformation, and Gaussian optimal transports, $\mu_0 \mapsto \mu_1$ on the left, and $\mu_1 \mapsto \mu_0$ on the right. Below are the averages of $x_{0,i}$'s and $x_{1,i}$'s, and of the transported points. The lines are geodesics in the dual spaces, mapped in the simplex. Optimal transport in $\mathbb{R}^2$, on $z_{0,i}$'s and $z_{1,i}$'s, can be visualized at the bottom (with linear mapping since Gaussian assumptions are made).</figcaption></figure>

Empirically, this can be performed using Algorithm 2, and a simulation can be visualized in Figure 2, where $h = \text{clr}$. On the left, we can visualize the mapping of red points to the blue distribution, and on the right, the “inverse mapping” of blue points to the red distribution. Transformed points $z = h(x)$, that are plotted at the bottom, are supposed to be normally distributed, and a multivariate Gaussian Optimal Transport mapping is used. Hence, $T_t^*$ is linear in $\mathbb{R}^{d-1}$, as given by expression 4, as well as displacement interpolation, corresponding to red and blue segments. But, as we can see on top of Figure 2, in the original space, $t \mapsto x_t := h^{-1}(z_t)$ will be nonlinear. Tables are average values of the three components of $x$'s and
**Algorithm 2 Gaussian Based Transport of $\mathbf{x}_0$ on $S_d$**

**Input:** $\mathbf{x}_0 \in S_d$

**Parameter:** $\{\mathbf{x}_{0,1}, \dots, \mathbf{x}_{0,n_0}\}$ and $\{\mathbf{x}_{1,1}, \dots, \mathbf{x}_{1,n_1}\}$ in $S_d$; isomorphic transformation $h: S_d \to \mathbb{R}^{d-1}$

**Output:** $\mathbf{x}_1$

for $i \in \{1, \dots, n_0\}$ do
$\mathbf{z}_{0,i} \leftarrow h(\mathbf{x}_{0,i})$

end for

for $i \in \{1, \dots, n_1\}$ do
$\mathbf{z}_{1,i} \leftarrow h(\mathbf{x}_{1,i})$

end for

$\mathbf{m}_0 \leftarrow \text{average of } \{\mathbf{z}_{0,1}, \dots, \mathbf{z}_{0,n_0}\}$

$\mathbf{m}_1 \leftarrow \text{average of } \{\mathbf{z}_{1,1}, \dots, \mathbf{z}_{1,n_1}\}$

$\mathbf{S}_0 \leftarrow \text{empirical variance matrix of } \{\mathbf{z}_{0,1}, \dots, \mathbf{z}_{0,n_0}\}$

$\mathbf{S}_1 \leftarrow \text{empirical variance matrix of } \{\mathbf{z}_{1,1}, \dots, \mathbf{z}_{1,n_1}\}$

$\mathbf{A} \leftarrow \mathbf{S}_0^{-1/2}(\mathbf{S}_0^{1/2}\mathbf{S}_1\mathbf{S}_0^{1/2})^{1/2}\mathbf{S}_0^{-1/2}$

$\mathbf{x}_1 \leftarrow h^{-1}(\mathbf{m}_1 + \mathbf{A}(h(\mathbf{x}_0) - \mathbf{m}_0))$

$T^*(\mathbf{x})$'s.

# 6 Optimal Transport for Measures on $S_d$

## 6.1 Theoretical Properties

A function $\psi: S_d \to \mathbb{R}$ is exponentially concave if $\exp[\psi]: S_d \to \mathbb{R}_+$ is concave. As a consequence, such a function $\psi$ is differentiable almost everywhere. Let $\nabla\psi$ and $\nabla_{\vec{u}}\psi$ denote, respectively, its gradient, and its directional derivative. Following Pal and Wong (2016, 2018, 2020), define an allocation map generated by $\psi$, $\pi_\psi: S_d \to S_d$ defined as

$$ \pi_{\psi}(\mathbf{x}) = [x_1(1 + \nabla_{\vec{e}_1-\mathbf{x}}\psi(\mathbf{x})), \dots, x_d(1 + \nabla_{\vec{e}_d-\mathbf{x}}\psi(\mathbf{x}))], $$

where $\{\vec{e}_1, \dots, \vec{e}_d\}$ is the standard orthonormal basis of $\mathbb{R}^d$. Consider the optimal transport problem with the following cost function, on $S_d \times S_d$, i.e., the L-divergence corresponding to the cross-entropy,

$$ c(\mathbf{x}, \mathbf{y}) = \log \left( \frac{1}{d} \sum_{i=1}^{d} \frac{y_i}{x_i} \right) - \frac{1}{d} \sum_{i=1}^{d} \log \left( \frac{y_i}{x_i} \right), \qquad (5) $$

called “Dirichlet transport” in Baxendale and Wong (2022). See Pistone and Shoaib (2024) for a discussion about the connections with the distance induced by Aitchsion’s inner product of Equation (3). From Theorem 1 in Pal and Wong (2020), for this cost function, there exists an exponentially concave function $\psi^*: S_d \to \mathbb{R}$ such that

$$ T^*(\mathbf{x}) = \mathbf{x} \diamond \pi_{\psi^*}(\mathbf{x}^{-1}) $$

defines a push-forward from $\mathbb{P}_0$ to $\mathbb{P}_1$, and the coupling $(\mathbf{x}, T^*(\mathbf{x}))$ is optimal for problem (1), and is unique if $\mathbb{P}_0$ is absolutely continuous. Observe that if $\mathbf{y} = T^*(\mathbf{x})$,

$$ \mathbf{y} = C(\pi_{\psi^*}(\mathbf{z})_1/z_1, \dots, \pi_{\psi^*}(\mathbf{z})_d/z_d), $$

where $\mathbf{z} = \mathbf{x}^{-1}$.
One can also consider an interpolation,

$$T_t^*(\mathbf{x}) = \mathbf{x} \diamond \pi_t(\mathbf{x}^{-1})$$

where $\pi_t = (1-t)d^{-1}\mathbf{1} + t\pi_{\psi^*}$ (even if this approach differs from McCann's displacement interpolation).

Note that a classical distribution on $S_d$ is Dirichlet distribution, with density

$$f(x_1, \dots, x_d; \alpha) = \frac{1}{B(\alpha)} \prod_{i=1}^{d} x_i^{\alpha_i - 1}$$

for some $\alpha = (\alpha_1, \dots, \alpha_d) \in \mathbb{R}_+^d$, and a normalizing constant denoted $B(\alpha)$. Level curves of the density of Dirichlet distributions fitted on our toy dataset can be visualized in Figure 3. Unfortunately, unlike the multivariate Gaussian distribution, there is no explicit expression for the optimal mapping between Dirichlet distribution (regardless of the cost). Therefore, to remain within $S_d$ and avoid the $\mathbb{R}^{d-1}$ representation, numerical techniques should be considered.

<figure><img src="image_7.png" /><figcaption>Figure 3: Densities of Dirichlet distributions in $S_3$ fitted on observations of the toy dataset of Figure 1.</figcaption></figure>

## 6.2 Matching

Consider two samples in the $S_d$ simplex, $\{\mathbf{x}_{0,1}, \dots, \mathbf{x}_{0,n_0}\}$ and $\{\mathbf{x}_{1,1}, \dots, \mathbf{x}_{1,n_1}\}$. The discrete version of the Kantorovich problem (corresponding to Equation 2) is

$$\min_{P \in U(n_0, n_1)} \left\{ \sum_{i=1}^{n_0} \sum_{j=1}^{n_1} P_{i,j} C_{i,j} \right\} \qquad (6)$$

where, as in Brualdi (2006), $U(n_0, n_1)$ is the set of $n_0 \times n_1$ matrices corresponding to the convex transportation polytope

$$U(n_0, n_1) = \left\{ P : P \mathbf{1}_{n_1} = \mathbf{1}_{n_0} \text{ and } P^{\top} \mathbf{1}_{n_0} = \frac{n_0}{n_1} \mathbf{1}_{n_1} \right\},$$
**Algorithm 3 Coupling samples on S<sub>d</sub>**

Input: {x<sub>0,1</sub>, ..., x<sub>0,n<sub>0</sub></sub>} and {x<sub>1,1</sub>, ..., x<sub>1,n<sub>1</sub></sub>} in S<sub>d</sub>;

Output: weight matching matrix n<sub>0</sub> × n<sub>1</sub> P<sup>*</sup>

C ← matrix n<sub>0</sub> × n<sub>1</sub>, C<sub>i,j</sub> = c(x<sub>i</sub>, x<sub>j</sub>) using (5)

P<sup>*</sup> ← solution of Equation (6), using LP libraries

and where C denotes the n<sub>0</sub> × n<sub>1</sub> cost matrix, C<sub>i,j</sub> = c(x<sub>i</sub>, x<sub>j</sub>), associated with cost from Equation (5).

In Algorithm 3, we recall how this procedure works, which is the one explained in Peyré et al. (2019), with a specific cost function (from Equation (5)). In the toy dataset, this can be visualized for two specific observations x<sub>0,i</sub>. If n<sub>0</sub> ≠ n<sub>1</sub>, it is not a one-to-one coupling, and “the counterfactual” is actually a weighted average of x<sub>1,j</sub>’s, where weights are given in row P<sup>*</sup><sub>i</sub> = [P<sup>*</sup><sub>i,1</sub>, ..., P<sup>*</sup><sub>i,n<sub>1</sub></sub>] ∈ S<sub>n<sub>1</sub></sub>.

<figure><img src="image_8.png" /><figcaption>Figure 4: Getting empirical counterfactuals using matching techniques, with x<sub>0,i</sub> in red (on the top-left hand-side), and counterfactuals x<sub>1,j</sub>’s in blue (bottom-right hand-side), with size proportional to P<sup>*</sup><sub>i</sub> = [P<sup>*</sup><sub>i,1</sub>, ..., P<sup>*</sup><sub>i,n<sub>1</sub></sub>] ∈ S<sub>n<sub>1</sub></sub>.</figcaption></figure>

# 7 Application on Sequential Transport for Counterfactuals

Variables x<sub>j</sub> in tabular data are either continuous or categorical. If x<sub>j</sub> is continuous, since x<sub>j</sub> ∈ ℝ, transporting from observed x<sub>j</sub>|s = 0 to counterfactual x<sub>j</sub>|s = 1 is performed using standard (conditional) monotonic mapping, as discussed in Fernandes Machado et al. (2025), using classical F<sub>1</sub><sup>-1</sup> ◦ F<sub>0</sub>. If x<sub>j</sub> is categorical, with d categories, consider some fitted model ̂m(x<sub>j</sub>|x<sub>-j</sub>), using some multinomial loss, and let ̂x<sub>j</sub> = ̂m(x<sub>j</sub>|x<sub>-j</sub>) denote the predicted scores, so that ̂x<sub>j</sub> ∈ S<sub>d</sub>. Then use Algorithm 2, with a Gaussian mapping in an Euclidean representation space, to transport from observed ̂x<sub>j</sub>|s = 0 to counterfactual ̂x<sub>j</sub>|s = 1, in S<sub>d</sub>.

## 7.1 German Credit: Purpose

In the popular German Credit dataset, from Hofmann (1994), the variable Purpose described the reason an individual took out a loan. This variable is an important predictor for
explaining potential defaults. The original variable is based on ten categories, that are merged here into three main classes, cars, equipment and other, in order to visualize the transport in a ternary plot — or Gibbs triangle. The sensitive variable *s* is here *Sex*.

We aim to construct a counterfactual value for the loan purpose, assuming the individuals were of a different sex. To achieve this, we apply our suggested procedure from to represent the purpose categorical variable as a compositional variable, using the same four models outlined in Section 3 and then apply Gaussian mapping from Section 5.2. The results provided by all of the models, shown in Figure 5, suggest that, had the individuals been of a different sex, the purpose of the loan would have changed. Specifically, if the average scores in each group (cars, equipment, and other) were approximately [35%, 45%, 20%] in the female population, after transporting to obtain the counterfactuals, the average scores become [31%, 52%, 18%], which closely resemble the actual frequencies of each category in the original male population.

One can also consider our second approach, using matching in $S_3$. Consider individual *i* among women, e.g., the left of Figure 6, $x_{0,i}$ = “equipment.” Using a MLR model, we obtain composition $\mathbf{x}_{0,i}$, here [11.38%, 79.30%, 9.32]. Using Algorithm 3, three points $\mathbf{x}_{1,j}$’s are matched, respectively with weights [0.453, 0.094, 0.453]. The first and the third individuals are such that $\mathbf{x}_{1,j}$ = “equipment” too, the second one “other.” So it would make sense to suppose that the counterfactual version of woman *i* with an “equipment” credit is a man with the same purpose. Actually, using Gaussian transport, $T^*(\mathbf{x}_{0,i}) = [15.78\%, 69.54\%, 14.68\%]$.

## 7.2 Adult: Marital Status

Following the numerical applications in Plečko et al. (2024) and Fernandes Machado et al. (2025), we consider here the Adult dataset, from Becker and Kohavi (1996). We regrouped categories of the Marital Status variable to create three generic ones (that can be visualized in a ternary plot, as in Figure 7), namely Married, Never-married and Separated. This example is interesting because if we compare status with respect to the Sex variable, proportions are quite different. In the dataset, proportions for married, never married, and separated are (roughly) [62%, 27%, 12%] for men, [14%, 44%, 41%] for women (more precise values are at the top of the table in Figure 7). Thus, the counterfactual of a “separated” woman is more likely to be a “married” man than a “separated” man. Four models are used to convert the categorical variable Marital Status into a composition, as previously. The first MLR is based on three variables: a categorical variable, occupation, and two continuous ones, age, and hours_per_week, modeled nonlinearly using b-splines (hence, it is referred to as a logistic GAM). This model is clearly underfitted. Therefore, observations $\mathbf{x}_{0,i}$’s for women and $\mathbf{x}_{1,i}$’s for men clearly are in the interior of $S_d$. In contrast, the more complex MLR (which uses additional features), as well as the random forest and boosting models, can produce predictions near the simplex boundary, $\partial S_d$.

For the underfitted model (top left), transported scores have a distribution very close to the ones in the population of men. For the more accurate MLR model (top right), proportions are very close to the actual proportions (which is not surprising since GLMs are usually well calibrated), but the transported scores are slightly different than the proportions of categories (proportions were [62%, 27%, 12%] while average transported scores are [67%, 25%, 8%]). At least, we are different from the original ones, but the mapping is not as accurate as it should be. This might come from the fact that when the points $\mathbf{x}_i$ are close to the border $\partial S_d$, it is quite unlikely that the sample $\mathbf{z}_i$ is Gaussian.
<figure><img src="image_9.png" /><figcaption>Figure 5: Optimal transport using the *clr* transformation, and Gaussian optimal transports, on the *purpose* scores in the *German Credit* database, with two logistic *GAM* models to predict scores, on top, and below a random forest (left) and a boosting model (right). Points in red are compositions for women, while points in blue are for men. Lines indicate the displacement interpolation when generating counterfactuals.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 6: Empirical matching of two women, in red, from the German Credit dataset, with 2 or 3 men, in blue. Size of blue dots are proportional to the weights $\mathbf{P}_i^*$.</figcaption></figure>

# 8 Conclusion

In this article, we introduce a novel approach for constructing counterfactuals for categorical data by transforming them into compositional data using a probabilistic classifier. Our approach avoids imposing arbitrary assumptions about label ordering. However, our methodology is not without limitations. OT computations, particularly on the simplex, can be computationally intensive for large-scale datasets, posing challenges in high-dimensional settings. Additionally, the reliance on a probabilistic classifier in the initial step introduces potential vulnerabilities. Biases may arise from a poorly calibrated or inaccurate classifier, impacting the quality of the subsequent analysis—especially with scarce categories that may need grouping to apply the proposed method.

# A Complexity of the Main Algorithms

The complexity of Algorithm 2 is $O(d^3+(n_0+n_1)d^2)$ primarily arising from the computation of the class-wise covariance matrices $S_0, S_1$, and the transformation matrix $A$. When the number of classes $d$ becomes large, this cost becomes prohibitive; in such scenarios, alternative OT formulations (beyond the Gaussian OT mapping) are advisable.

Alternatively, Algorithm 3 provides a more scalable approach in high-dimensional settings. Here, $d$ only affects the cost matrix computation with complexity $O(n_0n_1d)$, while the dominant computational burden lies in solving the OT problem, which scales as $O((n_0n_1)^{3/2})$ when using the Operator Splitting Quadratic Program (OSQP) solver.

# B Choice of the Cost Function in Algorithm 3

Section 6 introduces an approach to transporting probability measures directly on the simplex $S_d$. More precisely, Section 6.1 defines the existence of an OT map $T^*$ from a source measure $P_0$ on $S_d$ to a target measure $P_1$ on $S_d$, using the cross-entropy as the cost function, which is referred to as “Dirichlet transport” (Baxendale and Wong, 2022). While this cost function
<figure><img src="image_11.png" /><figcaption>Figure 7: Optimal transport using the *clr* transformation, and Gaussian optimal transports, on the *Marital Status* scores in the *Adult* database, with two logistic *GAM-MLR* models to predict scores, on top, and below a random forest (left) and a boosting model (right). Points in red are compositions for women, while points in blue are for men. Lines indicate the displacement interpolation when generating counterfactuals.</figcaption></figure>
guarantees the theoretical existence of an optimal map, no closed-form expression for $T^*$ is available in practice. Instead, numerical optimization is performed via the Kantorovich formulation of the OT problem to match individuals from $P_0$ to $P_1$ on $S_d$, as described in Section 6.2. Although Euclidean costs, such as the Wasserstein-1 or Wasserstein-2 distances, could be considered, they violate the geometry of the simplex by relying on absolute differences, unlike the cross-entropy, which accounts for relative proportions.

## C Integration of the methodology into a SCM for counterfactual fairness assessment

In Section 3, to transform a categorical variable $x$ into a numerical one, we suggest using a probabilistic classifier based on the features $X_{-x}$, i.e., all features except the categorical variable. In this appendix, we consider a more global approach, within an SCM, which requires positing a structural causal model in advance. Consider the DAG shown in Figure 8, whose topological order is

$$S \to X_1 \to X_2 \to X_3 \to Y,$$

where $S \in \{0, 1\}$ is the sensitive treatment (for example, binary gender), $X_1 \in \mathbb{R}$ is a numeric feature, $X_2 \in [d_2]$ and $X_3 \in [d_3]$ are categorical features with $d_2$ and $d_3$ categories, respectively, and $Y$ is the outcome sink. To generate a counterfactual for an individual with $S = 0$, we first flip the treatment by setting $S = 1$. We then obtain the numeric counterfactual $X_1^*$ using a chosen mapping mechanism (e.g., using optimal transport).

Subsequently, we construct each categorical counterfactual in sequence. For $X_2$, we fit a classifier on its parents $(S, X_1)$ to estimate the conditional distribution $\hat{P}(X_2 | S, X_1)$. This yields probability vectors in the simplex $S_{d_2}$ for both the factual group ($S = 0$) and the counterfactual group ($S = 1$), which needs to be converted to a single label $X_2^*$ in $\{1, \dots, d_2\}$. We suggest either sampling according to the probabilities or selecting the highest-probability category (top-label). We repeat the same process for $X_3$, estimating $\hat{P}(X_3 | S, X_2)$ and mapping its output into $\{1, \dots, d_3\}$. Finally, we predict the counterfactual outcome $Y^*$ with a model conditioned on $(S, X_1^*, X_2^*, X_3^*)$. This sequential, topologically ordered procedure embeds each probabilistic result into its required discrete space, ensuring consistency with the a priori SCM.

<figure><img src="image_12.png" /><figcaption>Figure 8: Example of a Structural Causal Model with a sensitive attribute $S$, a numeric variable $X_1$, two categorical variables $X_2$ and $X_3$, and an output variable $Y$.</figcaption></figure>

Matching individuals using their probability vectors on $S_d$ allows to uniquely determine them as it involves continuous distributions on $[0, 1]$ for each category. In contrast, if counterfactuals were computed directly from categorical data, one would need to rely on the non-deterministic counterfactual framework described in (Pearl et al., 2016, Chap. 4).
References

Aitchison, J. (1982). The statistical analysis of compositional data. *Journal of the Royal Statistical Society: Series B (Methodological)* 44: 139–160.

Baxendale, P. and Wong, T.-K. L. (2022). Random concave functions. *The Annals of Applied Probability* 32: 812–852.

Becker, B. and Kohavi, R. (1996). Adult. UCI Machine Learning Repository, doi:10.24432/C5XW20.

Black, E., Yeom, S. and Fredrikson, M. (2020). Fliptest: fairness testing via optimal transport. In *Proceedings of the 2020 Conference on Fairness, Accountability, and Transparency*, 111–121.

Brenier, Y. (1991). Polar factorization and monotone rearrangement of vector-valued functions. *Communications on pure and applied mathematics* 44: 375–417, doi:10.1002/cpa.3160440402.

Brualdi, R. A. (2006). *Combinatorial matrix classes*, 13. Cambridge University Press.

Charpentier, A., Flachaire, E. and Gallic, E. (2023). Optimal transport for counterfactual estimation: A method for causal inference. In *Optimal Transport Statistics for Economics and Related Topics*. Springer, 45–89.

Chayes, F. (1971). *Ratio correlation*. University of Chicago Press.

Cheng, W., Hüllermeier, E. and Dembczynski, K. J. (2010). Bayes optimal multilabel classification via probabilistic classifier chains. In *Proceedings of the 27th international conference on machine learning (ICML-10)*, 279–286.

De Lara, L., González-Sanz, A., Asher, N., Risser, L. and Loubes, J.-M. (2024). Transport-based counterfactual models. *Journal of Machine Learning Research* 25: 1–59.

Egozcue, J. J., Pawlowsky-Glahn, V., Mateu-Figueras, G. and Barcelo-Vidal, C. (2003). Isometric logratio transformations for compositional data analysis. *Mathematical geology* 35: 279–300.

Fernandes Machado, A., Charpentier, A. and Gallic, E. (2025). Sequential conditional transport on probabilistic graphs for interpretable counterfactual fairness. *Proceedings of the AAAI Conference on Artificial Intelligence* 39: 19358–19366, doi:10.1609/aaai.v39i18.34131.

Hallin, M., Del Barrio, E., Cuesta-Albertos, J. and Matrán, C. (2021). Distribution and quantile functions, ranks and signs in dimension d: A measure transportation approach. *The Annals of Statistics* 49: 1139–1165.

Hallin, M. and Konen, D. (2024). Multivariate quantiles: Geometric and measure-transportation-based contours. In *Applications of Optimal Transport to Economics and Related Topics*. Springer, 61–78.

Higham, N. J. (2008). *Functions of matrices: theory and computation*. SIAM.

Hofmann, H. (1994). German Credit Data. UCI Machine Learning Repository, doi:10.24432/C5NC77.

Imai, K., Olivella, S. and Rosenman, E. T. (2022). Addressing census data problems in race imputation via fully bayesian improved surname geocoding and name supplements. *Science Advances* 8: eadc9824, doi:10.1126/sciadv.adc9824.

Kantorovich, L. V. (1942). On the translocation of masses. In *Doklady Akademii Nauk USSR, 37*, 199–201.

King, G., Tanner, M. A. and Rosen, O. (2004). *Ecological inference: New methodological strategies*. Cambridge University Press.

Kusner, M. J., Loftus, J., Russell, C. and Silva, R. (2017). Counterfactual fairness. In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S. and Garnett, R. (eds), *Advances in Neural Information Processing Systems 30*. NIPS, 4066–4076.

McCann, R. J. (1997). A convexity principle for interacting gases. *Advances in mathematics* 128: 153–179.

Monge, G. (1781). *Mémoire sur la théorie des déblais et des remblais. Histoire de l'Académie*
Royale des Sciences de Paris.

Naeini, M. P., Cooper, G. and Hauskrecht, M. (2015). Obtaining well calibrated probabilities using bayesian binning. In *Proceedings of the AAAI conference on artificial intelligence*, 29.

Pal, S. (2024). On the difference between entropic cost and the optimal transport cost. *The Annals of Applied Probability* 34: 1003–1028.

Pal, S. and Wong, T.-K. L. (2016). The geometry of relative arbitrage. *Mathematics and Financial Economics* 10: 263–293.

Pal, S. and Wong, T.-K. L. (2018). Exponentially concave functions and a new information geometry. *The Annals of probability* 46: 1070–1113.

Pal, S. and Wong, T.-K. L. (2020). Multiplicative schrödinger problem and the dirichlet transport. *Probability Theory and Related Fields* 178: 613–654.

Pawlowsky-Glahn, V. and Egozcue, J. J. (2001). Geometric approach to statistical analysis on the simplex. *Stochastic Environmental Research and Risk Assessment* 15: 384–398.

Pearl, J. (2009). *Causality*. Cambridge university press.

Pearl, J., Glymour, M. and Jewell, N. P. (2016). *Causal inference in statistics: A primer*. John Wiley Sons.

Peyré, G., Cuturi, M. et al. (2019). Computational optimal transport: With applications to data science. *Foundations and Trends in Machine Learning* 11: 355–607.

Pistone, G. and Shoaib, M. (2024). A unified approach to aitchison's, dually affine, and transport geometries of the probability simplex. *Axioms* : 823.

Plečko, D. and Meinshausen, N. (2020). Fair data adaptation with quantile preservation. *Journal of Machine Learning Research* 21: 1–44.

Plečko, D., Bennett, N. and Meinshausen, N. (2024). Fairadapt: Causal reasoning for fair data preprocessing. *Journal of Statistical Software* 110: 1–35.

Santambrogio, F. (2015). *Optimal transport for applied mathematicians*. Springer.

Villani, C. (2003). *Topics in optimal transportation*, 58. American Mathematical Society.

Zadrozny, B. and Elkan, C. (2001). Obtaining calibrated probability estimates from decision trees and naive Bayesian classifiers. In *Proceedings of the Eighteenth International Conference on Machine Learning*, ICML ’01. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc., 609–616.