# On optimal transport with f-divergence regularization

Maxime Nicaise, Yaiza Bermudez, Samir M. Perlaza

► To cite this version:

Maxime Nicaise, Yaiza Bermudez, Samir M. Perlaza. On optimal transport with f-divergence regularization. RR-9607, Inria. 2026. hal-05454647

HAL Id: hal-05454647

https://inria.hal.science/hal-05454647v1

Submitted on 12 Jan 2026

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
informatics mathematics
Inria

# On optimal transport with f-divergence regularization

Maxime Nicaise, Yaiza Bermudez, and Samir M. Perlaza

## RESEARCH REPORT

### N° 9607

January 2026

Project-Team NEO

ISSN 0249-6399 ISRN INRIA/RR-9607-FR+ENG
<figure><img src="image_2.png" /></figure>
# On optimal transport with f-divergence regularization

Maxime Nicaise, Yaiza Bermudez, and Samir M. Perlaza

Project-Team NEO

Research Report n° 9607 — January 2026 — 36 pages

**Abstract:** This report establishes that an optimal transport (OT) problem regularized by a given f-divergence admits the same solution as another OT problem regularized by a different g-divergence, under an appropriate transformation of the cost function. This structural equivalence between OT problems regularized by distinct divergences, in the sense of sharing the same unique minimizer, is demonstrated within the framework of Polish spaces with bounded cost functions. Under these general assumptions, the existence of optimal potentials for the dual problem and the uniqueness of the optimal coupling for the primal problem are rigorously established.

**Key-words:** Optimization, Optimal Transport, Entropic Optimal Transport, Regularized Optimal Transport.

RESEARCH CENTRE
SOPHIA ANTIPOLIS – MÉDITERRANÉE

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex

<footnote>Maxime Nicaise, Yaiza Bermudez and Samir M. Perlaza are with INRIA, Centre Inria d'Université Côte d'Azur, Sophia Antipolis 06902, France. Samir M. Perlaza is also with the GAATI Laboratory at the Université de la Polynésie Française, Faaa 98702, French Polynesia; and also with the Electrical and Computer Engineering Department at Princeton University, Princeton N.J. 08544, USA.</footnote>
Sur le transport optimal avec régularisation par
$f$-divergence

**Résumé** : Ce rapport établit qu’un problème de transport optimal (TO) régularisé par une *f*-divergence donnée admet la même solution qu’un autre problème de TO régularisé par une *g*-divergence différente, moyennant une transformation appropriée de la fonction de coût. Cette équivalence structurelle entre problèmes de TO régularisés par des divergences distinctes, au sens où ils partagent le même unique minimiseur, est démontrée dans le cadre des espaces polonais avec des fonctions de coût bornées. Sous ces hypothèses générales, l’existence de potentiels optimaux pour le problème dual et l’unicité du couplage optimal pour le problème primal sont rigoureusement établies.

**Mots-clés** : Optimisation, Transport Optimal, Transport Optimal Entropique, Transport Optimal Régularisé.
# Contents

<table><tr><td>1</td><td>Introduction</td><td>4</td></tr><tr><td>1.1</td><td>Notations and preliminaries</td><td>5</td></tr><tr><td>1.2</td><td>φ-divergences</td><td>6</td></tr><tr><td>2</td><td>Optimal transport with φ-divergence regularization</td><td>8</td></tr><tr><td>2.1</td><td>Primal problem</td><td>9</td></tr><tr><td>2.1.1</td><td>Classical formulation</td><td>9</td></tr><tr><td>2.1.2</td><td>Equivalent formulation</td><td>9</td></tr><tr><td>2.2</td><td>Derivation of the dual problem</td><td>12</td></tr><tr><td>2.3</td><td>Weak duality</td><td>14</td></tr><tr><td>2.4</td><td>Attainment in the dual problem</td><td>16</td></tr><tr><td>2.4.1</td><td>(c, λ, φ)-transform</td><td>16</td></tr><tr><td>2.4.2</td><td>Existence of maximizing potentials</td><td>22</td></tr><tr><td>2.5</td><td>Attainment in the primal problem</td><td>26</td></tr><tr><td>2.5.1</td><td>Closed-form expression of the optimal coupling</td><td>26</td></tr><tr><td>2.5.2</td><td>Sufficient conditions</td><td>28</td></tr><tr><td>3</td><td>Equivalence of the problem to regularization on the family of φ-divergences</td><td>29</td></tr><tr><td>3.1</td><td>Main results</td><td>29</td></tr><tr><td>3.2</td><td>Examples</td><td>31</td></tr><tr><td>3.2.1</td><td>From KL to reverse KL regularization</td><td>32</td></tr><tr><td>3.2.2</td><td>From KL to Jensen-Shannon regularization</td><td>32</td></tr><tr><td>4</td><td>Conclusions</td><td>33</td></tr></table>
# 1 Introduction

Optimal transport (OT) is established as a fundamental tool for comparing probability distributions, with applications ranging from machine learning to computer graphics [1–6]. The choice of a suitable discrepancy measure between distributions is essential for tasks such as regression, classification, and generative modeling. In this respect, OT-based distances offer a significant advantage over classical divergences, such as the Kullback–Leibler divergence or other $f$-divergences, as demonstrated in [2, 4]. Indeed, the topology induced through Wasserstein distances is weaker and better adapted to comparing measures with non-overlapping supports or those supported on lower-dimensional sets [2, 7]. This property makes OT particularly attractive for modern data analysis.

Despite these appealing properties, from a computational perspective, solving the Kantorovich formulation of optimal transport involves tackling a large-scale linear program. This task often becomes intractable for high-dimensional or large-sample problems, as shown in [8]. To mitigate the computational complexity of the original problem and improve the sample complexity of the estimators, regularization is commonly employed. This ranges from the classic entropic regularization leading to the celebrated Sinkhorn algorithm [9], to quadratic [10] or structured regularizers [3]. However, more recent developments highlight the benefits of broadening this framework to general $f$-divergences [11, 12]. These schemes yield desirable properties such as sparsity or robustness to outliers. Parallel to these advances in optimal transport, similar variational problems are extensively studied within the field of statistical learning theory, particularly the statistical empirical risk minimization (ERM) framework [13–16]. Specifically, theoretical equivalence between variational problems regularized with two different $f$-divergences is established in [16, 17]. Therein, it is proved that an ERM problem regularized by a specific $f$-divergence is reformulable using a different divergence and a suitably transformed cost function, yielding the exact same minimizer.

This work establishes an equivalence between optimal transport problems with $f$-divergence regularizations, which is analogous to [16, Theorem 2] in the context of empirical risk minimization with $f$-divergence regularizations. More specifically, it is proved that two optimal transport problems regularized by different divergences admit the same unique optimal coupling, up to a suitable transformation of the cost function. This structural equivalence is demonstrated within the framework of Polish spaces with bounded cost functions. Under these general assumptions, the existence of optimal potentials for the dual problem and the uniqueness of the optimal coupling for the primal problem are rigorously established. Furthermore, a closed-form expression for the optimal coupling is obtained.
## 1.1 Notations and preliminaries

Given a function $f: \mathcal{X} \to \mathcal{Y}$, its effective domain is denoted by $\text{dom}(f) \triangleq \{x \in \mathcal{X}: f(x) < +\infty\}$; and its range is denoted by $\text{Ran}(f) \triangleq \{f(x): x \in \mathcal{X}\}$. Given a set $\mathcal{A}$, its interior is denoted by $\text{int}(\mathcal{A})$. Given two measures $\mu$ and $\nu$, the notation $\mu \ll \nu$ signifies that $\mu$ is absolutely continuous with respect to $\nu$. Across this report, $\mathcal{X}$ and $\mathcal{Y}$ are Polish spaces. The set of probability measures over $\mathcal{X}$ equipped with the $\sigma$-field of Borel sets is denoted by $\Delta(\mathcal{X})$. The set of $p$-integrable functions with respect to $P_X$ is denoted by $L^p(\mathcal{X}, P_X)$. Given two probability measures $P_X \in \Delta(\mathcal{X})$ and $P_Y \in \Delta(\mathcal{Y})$, the notation $\Pi(P_X, P_Y)$ denotes the set of probability measures on $\mathcal{X} \times \mathcal{Y}$ with marginals $P_X$ and $P_Y$ in $\Delta(\mathcal{X})$ and $\Delta(\mathcal{Y})$, respectively. The probability measures in $\Pi(P_X, P_Y)$ are called couplings. The product measure of $P_X \in \Delta(\mathcal{X})$ and $P_Y \in \Delta(\mathcal{Y})$ in $\Delta(\mathcal{X} \times \mathcal{Y})$ is denoted by $P_X P_Y$. Given two functions $f: \mathcal{X} \to \mathbb{R}$ and $g: \mathcal{X} \to \mathbb{R}$, the notation $f \stackrel{\text{a.s.}}{=} g$ means that $f$ and $g$ are almost surely equal with respect to $P_X$.

Similarly, the notation $f \stackrel{\text{a.s.}}{\ge} g$ means that $f$ is greater than $g$ almost surely with respect to $P_X$. Given $\mathcal{A}$ and $\mathcal{B}$, two subsets of $\mathcal{X}$, the notation $\mathcal{A} \stackrel{\text{a.s.}}{\subset} \mathcal{B}$ means that $P_X(\mathcal{A} \setminus \mathcal{B}) = 0$. Let $\mathcal{E}$ be a Banach space and $\mathcal{E}'$ be its topological dual space. Given $x \in \mathcal{E}$ and $y \in \mathcal{E}'$, $\langle y, x \rangle$ denotes the dual pairing between $y$ and $x$. If $f: \mathcal{X} \to \mathbb{R}$ and $g: \mathcal{Y} \to \mathbb{R}$ are two real functions then $f \oplus g: \mathcal{X} \times \mathcal{Y} \to \mathbb{R}$ is the function defined through $f \oplus g(x, y) = f(x) + g(y)$. The following definition formalizes the notion of admissibility for an optimization problem, ensuring that candidates yield a finite value within the functional domain.

**Definition 1.1.** Given a set $\mathcal{M}$ and a functional $R: \mathcal{M} \to \mathbb{R} \cup \{+\infty\}$, an instance $\theta \in \mathcal{M}$ is said to be admissible for the following problem

$$ \inf_{\nu \in \mathcal{M}} R(\nu), \qquad (1) $$

if $R(\theta) < +\infty$.

Analysis of optimality conditions and the derivation of dual problems rely heavily on the transformation of functions into a dual space representation. This operation is performed through the convex conjugate, providing a bridge between primal and dual variables in Banach spaces. The convex conjugate is defined below.

**Definition 1.2.** Let $\mathcal{E}$ be a Banach space and $\mathcal{E}'$ be its topological dual space. Let $\phi: \mathcal{E} \to \mathbb{R} \cup \{-\infty, +\infty\}$; the convex conjugate (or Fenchel conjugate) of $\phi$ is the function $\phi^*: \mathcal{E}' \to \mathbb{R} \cup \{-\infty, +\infty\}$ defined for $y \in \mathcal{E}'$ through (see [18, Section 2.3])

$$ \phi^*(y) = \sup_{x \in \mathcal{E}} (\langle y, x \rangle - \phi(x)). \qquad (2) $$
Complementary to the convex conjugate, the indicator function serves as a fundamental tool through which constraints are encoded directly into the objective functional. This formulation is particularly suited for optimal transport problems, as it ensures the strict enforcement of marginal constraints. The following definition formalises this concept, providing a mechanism to assign an infinite cost to any instance outside the feasible set.

**Definition 1.3.** Let $C$ be a subset of a set $\mathcal{E}$. The indicator function of $C$, denoted through $\iota_C$, is defined as:

$$ \iota_C(x) = \begin{cases} 0 & \text{if } x \in C; \\ +\infty & \text{if } x \notin C. \end{cases} \qquad (3) $$

## 1.2 $\phi$-divergences

The purpose of this section is to define the class of $\phi$-divergences endowed with sufficient regularity to guarantee satisfactory theoretical properties when used to regularize optimal transport problems. Let $\phi$ be a convex, proper, lower semi-continuous function from $\mathbb{R}$ to $(-\infty, +\infty]$. Throughout this report, the real numbers $a_\phi$ and $b_\phi$ denote the following quantities

$$ a_{\phi} \triangleq \inf \mathrm{dom}(\phi); \quad \text{and} \qquad (4) $$

$$ b_{\phi} \triangleq \sup \mathrm{dom}(\phi). \qquad (5) $$

As stated in [19, Lemma 2.1], the following limits exist

$$ \alpha_{\phi} \triangleq \lim_{x \to -\infty} \frac{\phi(x)}{x}; \quad \text{and} \qquad (6) $$

$$ \beta_{\phi} \triangleq \lim_{x \to +\infty} \frac{\phi(x)}{x}. \qquad (7) $$

Note that if $\mathrm{dom}(\phi) = (0, +\infty)$, then by convention $\alpha_\phi = -\infty$. It is shown in [19, Lemma 2.1] that the convex conjugate of $\phi$, which is denoted by $\phi^*$, verifies

$$ \mathrm{int}(\mathrm{dom}(\phi^*)) = (\alpha_{\phi}, \beta_{\phi}). \qquad (8) $$

In the analysis of regularized variational problems, the analytical properties of the convex conjugate $\phi^*$ play a pivotal role. To ensure that the dual problem admits well-defined solutions and to facilitate the reconstruction of the unique primal optimizer via the gradient of the conjugate, it is standard to restrict attention to a class of functions that ensures a smooth one-to-one correspondence between primal and dual variables. This is achieved by the class of Legendre functions, which combines strict convexity with essential smoothness. The steep behavior of their gradients at the boundary of the domain naturally enforces constraints, while the invertibility of the gradient map is instrumental in deriving closed-form expressions for the optimal coupling. The following definition characterizes functions of Legendre type.
**Definition 1.4** ([19], Definition 2.5). A function $\phi$ is said to be of Legendre type if it is strictly convex on $\text{dom}(\phi)$ and differentiable on $\text{int}(\text{dom}(\phi))$, and its derivative $\phi^{*'}$: $\text{int}(\text{dom}(\phi)) \to \mathbb{R}$ satisfies

$$ \lim_{x \to a_{\phi}} \phi'(x) = -\infty, \text{ if } a_{\phi} > -\infty; \text{ and} \quad (9) $$

$$ \lim_{x \to b_{\phi}} \phi'(x) = +\infty, \text{ if } b_{\phi} < +\infty. \quad (10) $$

The specific class of Legendre functions possesses certain properties regarding conjugation (Definition 1.2), as detailed in the following lemma.

**Lemma 1.1** ([19], Lemma 2.6). If $\phi: \mathbb{R} \to (-\infty, +\infty]$ is of Legendre type, then so is its convex conjugate $\phi^*$. Moreover $\phi': (a_{\phi}, b_{\phi}) \to (\alpha_{\phi}, \beta_{\phi})$ and $\phi^{*'}: (\alpha_{\phi}, \beta_{\phi}) \to (a_{\phi}, b_{\phi})$ are continuous, strictly increasing and mutually inverse maps.

In the context of this work, the focus is on a specific subset of Legendre functions that are suitable for defining divergences between probability measures, which are defined as follows.

**Definition 1.5** (Legendre-type generator). A Legendre type function $\phi: \mathbb{R} \to (-\infty, +\infty]$ is said to be a Legendre-type generator when $\phi(1) = 0$; $a_{\phi} = 0$ and $b_{\phi} = +\infty$, with $a_{\phi}$ in (4) and $b_{\phi}$ in (5); and its derivative satisfies $\phi'(1) = 0$.

Given a Legendre-type generator $\phi$, it follows from (8) that the interior of the domain of $\phi^*$ is $(-\infty, \beta_{\phi})$, and $\phi^{*'}: (-\infty, \beta_{\phi}) \to (0, +\infty)$ is a one to one strictly increasing map, the inverse of which is $\phi': (0, +\infty) \to (-\infty, \beta_{\phi})$. Note that a Legendre-type generator $\phi$ is always positive, as it is convex and verifies $\phi'(1) = 0$, thus its minimum is achieved at 1, and the function satisfies $\phi(1) = 0$. Moreover, note that the condition $\phi(1) = \phi'(1) = 0$ is not restrictive, as $\phi$-divergences are invariant up to an affine transformation [20]. The following definition specifies the notion of $\phi$-divergence adopted throughout this work.

**Definition 1.6.** Given a Legendre-type generator $\phi$, a measurable space $(\mathcal{E}, \mathcal{F})$, and two probability measures $P$ and $Q$ in $\Delta(\mathcal{E})$, the $\phi$-divergence of $P$ with respect to $Q$ is

$$ D_{\phi}(P||Q) = \begin{cases} \int \phi(\frac{dP}{dQ}(\theta))dQ(\theta), & \text{if } P \ll Q; \\ +\infty, & \text{otherwise.} \end{cases} \quad (11) $$

Although Definition 1.6 is not standard when the function $\phi$ is not assumed to be superlinear (i.e., $\lim_{x \to +\infty} \phi(x)/x = +\infty$). In the typical definition, certain probability measures may yield a finite $\phi$-divergence value even when they are not absolutely continuous with respect to the reference measure; see [21–23]. However, when a $\phi$-divergence is used as a regularization term in an optimization problem, allowing such non-absolutely continuous measures may lead to situations in which the optimizer (the measure minimizing the objective function) no longer admits a closed-form expression, as there may exist a non-zero singular part, as illustrated in [12]
**Example 1.1.** The function $\phi(x)$ is the **generator function** (or kernel) that defines the properties of the resulting $\phi$-divergence. Table 1 lists the divergences that satisfy these criteria, while Table 2 covers those that do not, typically due to a lack of strict convexity or the absence of the steepness property.

<table><caption>Table 1: Divergences satisfying Definition 1.5.</caption><thead><tr><td>Divergence</td><td>Generator φ(x)</td></tr></thead><tbody><tr><td>KL</td><td>x log(x) - x + 1</td></tr><tr><td>Reverse KL</td><td>-log(x) + x - 1</td></tr><tr><td>Power/Csiszár (0 ≤ α ≤ 1)</td><td><sup>x<sup>α</sup> - α(x-1) - 1</sup>⁄<sub>α(α-1)</sub></td></tr><tr><td>Neyman χ<sup>2</sup></td><td><sup>(x-1)<sup>2</sup></sup>⁄<sub>x</sub></td></tr><tr><td>Hellinger squared</td><td>(√x - 1)<sup>2</sup></td></tr><tr><td>Jeffreys</td><td>(x - 1) log(x)</td></tr><tr><td>Jensen-Shannon</td><td>x log(x) - (x + 1) log (<sup>x+1</sup>⁄<sub>2</sub>)</td></tr></tbody></table>

<table><caption>Table 2: Divergences failing to satisfy Definition 1.5.</caption><thead><tr><td>Divergence</td><td>Generator φ(x)</td></tr></thead><tbody><tr><td>Pearson χ<sup>2</sup></td><td>(x - 1)<sup>2</sup></td></tr><tr><td>Power/Csiszár (α > 1)</td><td><sup>x<sup>α</sup> - α(x-1) - 1</sup>⁄<sub>α(α-1)</sub></td></tr><tr><td>Triangular discrimination</td><td><sup>(x-1)<sup>2</sup></sup>⁄<sub>x+1</sub></td></tr><tr><td>Total variation</td><td>|x - 1|</td></tr></tbody></table>

## 2 Optimal transport with $\phi$-divergence regularization

This section establishes the mathematical formulation of optimal transport problem regularized by a $\phi$-divergence. The primal problem is introduced, as well as its dual counterpart, and the necessary functional analytic tools. Throughout this section, the function $\phi$ is a Legendre-type generator (Definition 1.5). The measure $P_{XY}$ in $\Delta(\mathcal{X} \times \mathcal{Y})$ is assumed to exhibit given marginal measures denoted by $P_X \in \Delta(\mathcal{X})$ and $P_Y \in \Delta(\mathcal{Y})$. The constant $\lambda$ is a strictly
positive real number, and the function

$$c: \mathcal{X} \times \mathcal{Y} \to \mathbb{R} \qquad (12)$$

is a bounded Borel measurable cost function. These elements constitute the parameters of the optimal transport problem with φ-divergence regularization, presented in [11] and [12].

## 2.1 Primal problem

This section is dedicated to establishing a rigorous mathematical setting for the primal problem. First, the classical formulation involving measures is presented, before introducing an equivalent formulation on Lebesgue function spaces that facilitates the derivation of the dual problem.

### 2.1.1 Classical formulation

The optimal transport problem regularized by a φ-divergence is defined in [11] and [12] as

$$\inf_{P_{XY} \in \Pi(P_X, P_Y)} \int_{\mathcal{X} \times \mathcal{Y}} c(x, y) dP_{XY}(x, y) + \lambda D_{\phi}(P_{XY} \| P_X P_Y), \quad (13)$$

where the function c is defined in (12). From Definition 1.6, any measure P_XY that is not absolutely continuous with respect to P_X P_Y yields an infinite divergence. Consequently, the present analysis is restricted to measures P_XY that are absolutely continuous with respect to the product measure P_X P_Y. To derive the associated dual problem directly, it is advantageous to express this formulation within the setting of [18, Section 2.8], which is the purpose of the following section.

### 2.1.2 Equivalent formulation

While the classical formulation in (13) is intuitive, it is mathematically more convenient to work with an equivalent formulation on Lebesgue function spaces. This change of space allows for a simplified use of standard results from convex analysis to derive the dual problem. The problem will be formulated in the general form given in [18, Section 2.8], with two functionals F and G and a linear operator A. The subsequent analysis explicitly constructs these objects, demonstrating how the operator A and the functionals F and G naturally emerge from the marginal constraints and the regularized objective function, respectively. A non-negative function r ∈ L¹(\mathcal{X} × \mathcal{Y}, P_X P_Y) represents the Radon-Nikodym derivative of a measure P_XY ∈ Π(P_X, P_Y) with respect to the product measure P_X P_Y if and only if it satisfies the following specific marginal constraints

$$\int_{\mathcal{Y}} r(x, y)dP_{\mathcal{Y}}(y) \stackrel{\text{a.s.}}{=} 1; \text{ and} \qquad (14)$$

$$\int_{\mathcal{X}} r(x, y)dP_{\mathcal{X}}(x) \stackrel{\text{a.s.}}{=} 1. \qquad (15)$$
See for instance [24, Theorem 10]. To formalize these constraints, the marginal operator A is defined as follows

$$A: \begin{cases} \mathrm{L}^1(\mathcal{X} \times \mathcal{Y}, P_X P_Y) \to \mathrm{L}^1(\mathcal{X}, P_X) \times \mathrm{L}^1(\mathcal{Y}, P_Y) \\ r \mapsto \left( \int_{\mathcal{Y}} r(\cdot, y) \mathrm{d}P_Y(y), \int_{\mathcal{X}} r(x, \cdot) \mathrm{d}P_X(x) \right), \end{cases} \quad (16)$$

where the notation $r(\cdot, y)$ (resp. $r(x, \cdot)$) denotes the function $x \mapsto r(x, y)$ (resp. $y \mapsto r(x, y)$). The product space $\mathrm{L}^1(\mathcal{X}, P_X) \times \mathrm{L}^1(\mathcal{Y}, P_Y)$ is equipped with the norm

$$\|(f,g)\|_{\mathrm{L}^1(\mathcal{X},P_X)\times\mathrm{L}^1(\mathcal{Y},P_Y)} = \|f\|_{\mathrm{L}^1(\mathcal{X},P_X)} + \|g\|_{\mathrm{L}^1(\mathcal{Y},P_Y)}. \quad (17)$$

**Lemma 2.1.** The linear operator A in (16) is continuous with respect to the norm in (17).

*Proof:* The continuity of the operator A follows from Fubini's Theorem [25, Theorem 8.8], which guarantees that the $\mathrm{L}^1$-norms of the resulting marginals are bounded by the $\mathrm{L}^1$-norm of the input function $r$. ■

From [26, Theorem 1.10.13], the dual of the product space endowed with the norm in (17) is the space

$$\mathrm{L}^\infty(\mathcal{X}, P_X) \times \mathrm{L}^\infty(\mathcal{Y}, P_Y), \quad (18)$$

equipped with the operator norm

$$\|(f,g)\|_{\infty} = \max (\|f\|_{\infty}, \|g\|_{\infty}). \quad (19)$$

Let the functional G be defined as

$$G: \begin{cases} \mathrm{L}^1(\mathcal{X}, P_X) \times \mathrm{L}^1(\mathcal{Y}, P_Y) \to \mathbb{R} \cup \{+\infty\} \\ (f,g) \mapsto \iota_{(1,1)}(f,g), \end{cases} \quad (20)$$

where $\iota_{(1,1)}$ denotes the indicator function of the singleton $\{(1,1)\}$ (Definition 1.3). The functional G serves as a rigorous mechanism to enforce the marginal constraints. Since it acts as an indicator function, it takes the value 0 if its argument is exactly the pair of constant functions $(1,1)$ and $+\infty$ otherwise. Consequently, when composed with the operator A, the term $G(A(r))$ is finite if and only if $r$ possesses the correct marginals. Additionally, the functional $F: \mathrm{L}^1(\mathcal{X} \times \mathcal{Y}, P_X P_Y) \to \mathbb{R} \cup \{+\infty\}$ is defined as follows

$$F(r) = \int_{\mathcal{X} \times \mathcal{Y}} c(x, y) r(x, y) \mathrm{d}P_X P_Y(x, y) \\ + \lambda \int_{\mathcal{X} \times \mathcal{Y}} \phi(r(x, y)) \mathrm{d}P_X P_Y(x, y). \quad (21)$$

The functional F encapsulates the objective function of the problem in (13). The first integral represents the linear transportation cost associated with the coupling density $r$, while the second integral corresponds to the regularization penalty defined by the divergence $\phi$, as defined in Definition 1.6. Using the
definition of A in (16), the definition of G in (20) and the definition of F in (21), the primal problem is formulated on L¹($\mathcal{X}$ × $\mathcal{Y}$, $P_X P_Y$) as

$$ \mathcal{P} \triangleq \inf_{r \in L^1(\mathcal{X} \times \mathcal{Y}, P_X P_Y), r \ge 0} \mathbf{P}(r), \quad (22) $$

$$ \text{with } \mathbf{P}(r) \triangleq F(r) + G(A(r)). \quad (23) $$

The structure of $\mathbf{P}(r)$ as the sum $F(r) + G(A(r))$ is canonical in convex analysis for handling constrained optimization. Minimizing this sum effectively seeks a function $r$ that minimizes the objective $F(r)$ while satisfying the constraint $A(r) = (1, 1)$ imposed by the indicator $G$, as any violation of the constraint results in an infinite objective value. This section mainly focuses on the problem in (22). The following lemma rigorously justifies that the two problems in (13) and (22) are equivalent by establishing a bijection between the feasible sets of the two formulations.

**Lemma 2.2.** There exists a one-to-one correspondence between the admissible functions of the problem in (22), in the sense of Definition 1.1, and couplings in $\Pi(P_X, P_Y)$ that are absolutely continuous with respect to $P_X P_Y$.

*Proof:* Let $r \in L^1(\mathcal{X} \times \mathcal{Y}, P_X P_Y)$ be admissible for the functional $\mathbf{P}$ in (23), in the sense of Definition 1.1. It will be shown that $r$ represents the Radon-Nikodym derivative of some probability measure in $\Pi(P_X, P_Y)$. The function $r$ being admissible for $\mathbf{P}$ is equivalent to requiring that $r$ satisfies $r \stackrel{\text{a.s.}}{\ge} 0$ and

$$ \int r(x, y) \, dP_Y(y) \stackrel{\text{a.s.}}{=} 1; \text{ and} \quad (24) $$

$$ \int r(x, y) \, dP_X(x) \stackrel{\text{a.s.}}{=} 1. \quad (25) $$

Consider the probability measure $P_{XY}$ defined by

$$ P_{XY}(\mathcal{A} \times \mathcal{B}) = \int_{\mathcal{A} \times \mathcal{B}} r(x, y) \, dP_X P_Y(x, y). \quad (26) $$

Then $P_{XY} \ll P_X P_Y$, and for any measurable set $\mathcal{A} \subset \mathcal{X}$,

$$ P_{XY}(\mathcal{A} \times \mathcal{Y}) = \int_{\mathcal{A}} \left( \int_{\mathcal{Y}} r(x, y) \, dP_Y(y) \right) dP_X(x) = P_X(\mathcal{A}), \quad (27) $$

which follows from the equality in (25). Similarly, the constraint over $\mathcal{Y}$ is verified: $P_{XY}(\mathcal{X} \times \mathcal{B}) = P_Y(\mathcal{B})$. Hence $P_{XY}$ is in $\Pi(P_X, P_Y)$. Conversely, let $P_{XY} \in \Pi(P_X, P_Y)$ with $P_{XY} \ll P_X P_Y$. The Radon-Nikodym theorem [25, Theorem 6.10] implies the existence of an almost surely unique nonnegative function

$$ r(x, y) = \frac{dP_{XY}}{dP_X P_Y}(x, y), \quad (28) $$
which belongs to L¹($\mathcal{X} \times \mathcal{Y}$, $P_X P_Y$). It will be shown that such an $r$ is admissible for the primal functional $P$ in 22 in the sense of Definition 1.1. Using (27) and (28), for any measurable $A \subset \mathcal{X}$, it holds that

$$P_X(A) = P_{XY}(A \times \mathcal{Y}) \qquad (29)$$

$$= \int_{A \times \mathcal{Y}} r(x, y) dP_X P_Y(x, y) = \int_A \left( \int_{\mathcal{Y}} r(x, y) dP_Y(y) \right) dP_X(x), (30)$$

which yields $\int_{\mathcal{Y}} r(x, y) dP_Y(y) \stackrel{\text{a.s.}}{=} 1$. A symmetric argument leads to the same result over $\mathcal{Y}$, namely $\int_{\mathcal{X}} r(x, y) dP_X(x) \stackrel{\text{a.s.}}{=} 1$. Thus $r$ is admissible for the functional $P$ in (23), in the sense of Definition 1.1. ■

Note that the correspondence in Lemma 2.2 is between the measure $P_{XY}$ and the equivalence class of functions $r$ in L¹ (uniqueness almost everywhere).

## 2.2 Derivation of the dual problem

Having established the primal problem on L¹ spaces, this section introduces the corresponding dual problem in the sense of [18, Section 2.8]. The derivation involves computing the convex conjugates of the functionals constituting the primal objective $P$ in (23). From [18, Section 2.8], the dual problem of the primal in (22) is formulated in its canonical Fenchel-Rockafellar form as follows

$$\mathcal{D} \triangleq \sup_{(f,g) \in L^\infty(\mathcal{X}, P_X) \times L^\infty(\mathcal{Y}, P_Y)} \mathrm{D}(f,g), \qquad (31)$$

$$\text{with } \mathrm{D}(f,g) \triangleq -G^*(f,g) - F^*(-A^*(f,g)), \qquad (32)$$

where $A^*: L^\infty(\mathcal{X}, P_X) \times L^\infty(\mathcal{Y}, P_Y) \to L^\infty(\mathcal{X} \times \mathcal{Y}, P_X P_Y)$ is the adjoint operator of $A$, defined in (16) and $G^*$ and $F^*$ are respectively the convex conjugates (definition 1.2) of $G$ and $F$, defined in (20) and (21), respectively. Any pair of functions solving the problem in (31) is referred to as a pair of optimal potentials. The explicit form of the dual problem is obtained by computing the conjugates and the adjoint operator, as detailed in the following lemma.

**Lemma 2.3.** For all $(f,g) \in L^\infty(\mathcal{X}, P_X) \times L^\infty(\mathcal{Y}, P_Y)$ and for all $s \in L^\infty(\mathcal{X} \times \mathcal{Y}, P_X P_Y)$, the operator $A$ in (16), the functional $G$ in (20) and the functional $F$ in (21), satisfy

$$A^*(f,g) = f \oplus g; \qquad (33)$$

$$G^*(f,g) = \int_{\mathcal{X}} f dP_X + \int_{\mathcal{Y}} g dP_Y; \qquad (34)$$

$$F^*(s) = \lambda \int_{\mathcal{X} \times \mathcal{Y}} \phi^* \left( \frac{s(x, y) - c(x, y)}{\lambda} \right) dP_X P_Y(x, y), \qquad (35)$$

where the function $c$ is defined in (12).
*Proof:* The proof is divided in three parts, corresponding to the proofs of each of the three equalities in (33), (34) and (35). To compute the adjoint of the linear operator A in (16), let $r$ be a function in $L^1(\mathcal{X} \times \mathcal{Y}, P_X P_Y)$ and $(f, g)$ be two functions in $L^\infty(\mathcal{X}, P_X) \times L^\infty(\mathcal{Y}, P_Y)$. The dual pairing of $(f, g)$ with A($r$) is given by

$$ \langle (f, g), A(r) \rangle = \int_{\mathcal{X}} f(x) \left( \int_{\mathcal{Y}} r(x, y) dP_Y(y) \right) dP_X(x) + \int_{\mathcal{Y}} g(y) \left( \int_{\mathcal{X}} r(x, y) dP_X(x) \right) dP_Y(y) \quad (36) $$

$$ = \int_{\mathcal{X} \times \mathcal{Y}} (f(x) + g(y)) r(x, y) dP_X P_Y(x, y) \quad (37) $$

$$ = \langle f \oplus g, r \rangle, \quad (38) $$

thus $A^*(f, g) = f \oplus g$, which yields (33).

The functional $G$ in (20) is the indicator of the singleton $\{(1,1)\}$. Classical results in convex analysis (e.g., [27, Example 13.2]) imply that for a pair $(f, g)$ of functions in $L^\infty(\mathcal{X}, P_X) \times L^\infty(\mathcal{Y}, P_Y)$,

$$ G^*(f, g) = \langle (f, g), (1, 1) \rangle = \int_{\mathcal{X}} f dP_X + \int_{\mathcal{Y}} g dP_Y, \quad (39) $$

which yields (34).

To compute the conjugate of the functional $F$, let $s$ be a function in $L^\infty(\mathcal{X} \times \mathcal{Y}, P_X P_Y)$. From Definition 1.2, the conjugate of the functional $F$ evaluated at $s$ is

$$ F^*(s) = \sup_{r \in L^1(\mathcal{X} \times \mathcal{Y}, P_X P_Y)} \left[ \begin{aligned} & \int_{\mathcal{X} \times \mathcal{Y}} s(x, y) r(x, y) dP_X P_Y(x, y) \\ & - \int_{\mathcal{X} \times \mathcal{Y}} c(x, y) r(x, y) dP_X P_Y(x, y) \\ & - \lambda \int_{\mathcal{X} \times \mathcal{Y}} \phi(r(x, y)) dP_X P_Y(x, y) \end{aligned} \right]. \quad (40) $$

Consider an element $(x, y)$ in $\mathcal{X} \times \mathcal{Y}$; the integrand $s(x, y)r(x, y) - c(x, y)r(x, y) - \lambda\phi(r(x, y))$ in (40) can be maximized pointwise, which yields

$$ \sup_{t \ge 0} [s(x, y)t - c(x, y)t - \lambda\phi(t)] = \lambda \sup_{t \ge 0} \left[ t \frac{s(x, y) - c(x, y)}{\lambda} - \phi(t) \right] \quad (41) $$

$$ = \lambda \phi^* \left( \frac{s(x, y) - c(x, y)}{\lambda} \right), \quad (42) $$

where $\phi^*$ is the convex conjugate of $\phi$ (Definition 1.2). The function obtained in (42) is measurable as $\phi^*$ is a sup of measurable functions. Thus, it holds that

$$ F^*(s) = \lambda \int_{\mathcal{X} \times \mathcal{Y}} \phi^* \left( \frac{s(x, y) - c(x, y)}{\lambda} \right) dP_X P_Y(x, y), \quad (43) $$
which yields (35). This concludes the proof.

The following corollary allows for a more explicit expression of the dual objective in (32).

**Corollary 2.1.** The dual functional D in (32) can be written more explicitly as

$$D(f,g) = \int_X f(x)dP_X(x) + \int_Y g(y)dP_Y(y) \\ -\lambda \int_{X \times Y} \phi^* \left( \frac{f(x) + g(y) - c(x,y)}{\lambda} \right) dP_X P_Y(x,y), \quad (44)$$

where $\phi^*$ is the convex conjugate of the function $\phi$ (Definition 1.2).

*Proof:* Plugging the expressions in (33), (34) and (35) in (32) yields (44).

## 2.3 Weak duality

Establishing the validity of weak duality requires the functionals constituting the primal objective in (22) to be proper and convex. These properties ensure that the optimization landscape is regular enough to apply standard tools from convex analysis and that the problem is not degenerate (i.e., the objective function is not identically equal to either $+\infty$ or $-\infty$). This is proven in the following lemma.

**Lemma 2.4.** The functionals $G$ in (20) and $F$ in (21) are proper and convex.

*Proof:* The functional $G$ is the indicator of a nonempty convex set, that is the singleton $\{(1,1)\}$, which implies it is proper and convex. Regarding the functional $F$, note first that it is proper. Indeed, since the function $c$ in (12) is bounded and the function $\phi$ is nonnegative (Definition 1.5), monotonicity of the integral guarantees that $F$ never attains the value $-\infty$. Moreover, evaluating $F$ at the constant function 1 yields a finite value, so that $F$ is not identically equal to $+\infty$. Convexity of $F$ then follows directly from the linearity of the integral and the convexity of $\phi$.

The following lemma establishes the relationship between the optimal values of the primal and dual problems, respectively in (22) and (31). It acts as a consistency check for the optimization setting, confirming that weak duality holds (the dual value lower-bounds the primal value) and, crucially, that both values are finite numbers, ensuring the problem is not ill-posed.

**Lemma 2.5.** Both $\mathcal{P}$ in (23) and $\mathcal{D}$ in (32) are finite numbers and weak duality holds between the two, that is $\mathcal{P} \ge \mathcal{D}$.

*Proof:* The inequality between $\mathcal{P}$ and $\mathcal{D}$ follows from Proposition 1.1 in [28], which requires $A$ in (16) to be continuous and the functionals $G$ in (20) and $F$ in (21) to be convex and proper. The first requirement is established in Lemma 2.1, while the second follows from Lemma 2.4. Therefore, [28, Proposition 1.1] implies that $\mathcal{P} \ge \mathcal{D}$.
The second part of the proof demonstrates that both values are finite. Regarding the primal value, every admissible $r \in L^1(\mathcal{X} \times \mathcal{Y}, P_X P_Y)$ must be non-negative and satisfy the marginal constraints (24) and (25) imposed by $G$ (Lemma 2.2); in particular, $r$ must be a probability density function with respect to $P_X P_Y$. Given such a function $r$, denote by $P^{(r)}$ the associated probability measure; the following identity then holds

$$ \int_{\mathcal{X} \times \mathcal{Y}} \phi(r(x, y)) dP_X P_Y(x, y) = D_{\phi}(P^{(r)} \| P_X P_Y) \geq 0. \quad (45) $$

This follows from Definition 1.6 of $D_{\phi}$, and the non-negativity follows from [29, Theorem 4.1]. Given that the function $c$ in (12) is bounded, there exists a constant $m$ such that

$$ \int_{\mathcal{X} \times \mathcal{Y}} c(x, y) r(x, y) dP_X P_Y(x, y) \geq m \int_{\mathcal{X} \times \mathcal{Y}} r(x, y) dP_X P_Y(x, y) \quad (46) $$

$$ \geq m, \quad (47) $$

where (47) holds because $r$ is a probability density function, as it satisfies the marginal constraints (24) and (25). Consequently, $P(r) \geq m$ for all admissible $r$, with $P$ in (23), which shows that $\mathcal{P} > -\infty$. Moreover, the admissible set is nonempty, since it contains the constant function 1. Consequently, the primal value $\mathcal{P}$ is finite, and the primal problem is feasible and bounded.

Regarding the dual value, the bound $\mathcal{D} \leq \mathcal{P}$ implies $\mathcal{D} < +\infty$. To ensure the dual problem does not degenerate to $-\infty$, a specific dual candidate with a finite value is exhibited. Given that $\phi$ is of Legendre type and $c$ is bounded, it is possible to choose constant functions $f$ and $g$ such that the range of $(x, y) \mapsto (f(x) + g(y) - c(x, y))/\lambda$ is contained in a bounded subset of $\text{dom}(\phi^*)$ of the form $(\tau, \gamma)$ with $\gamma < \beta_{\phi}$, where $\beta_{\phi}$ is in (7). From Lemma 1.1 and Definition 1.5, the function $\phi^*$ is continuous on $(-\infty, \beta_{\phi})$. Thus, with such choices of $f$ and $g$, the map $(x, y) \mapsto \phi^*((f(x)+g(y)-c(x, y))/\lambda)$ is bounded, and therefore $\text{D}(f, g)$ is a finite real number, where $\text{D}$ is in (44). This proves that $\mathcal{D} > -\infty$.

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

However, the required conditions in [18, Theorem 2.8.3] to prove that strong duality holds (i.e., $\mathcal{P} = \mathcal{D}$) are not fully met in this context. Indeed, it has been demonstrated in [12] that without stronger growth assumptions on $\phi$ (such as superlinearity), ensuring strong duality typically requires extending the domain of the divergence to assign finite values even to some measures that are not absolutely continuous with respect to the reference. However, in the latter, as mentioned in the introduction of this work, the optimal coupling no longer has a closed form expression in general. Since the present framework restricts the primal problem to densities in $L^1$, an alternative analytical approach is required to establish the existence of optimal potentials and the equality of optimal values.
## 2.4 Attainment in the dual problem

This section demonstrates that the supremum in the dual problem is always attained. The proof relies on a regularization-specific transform, analogous to the c-transform in classical optimal transport, which allows for the construction of maximizing sequences with favorable compactness properties.

### 2.4.1 (c, λ, φ)-transform

The concept of the (c, λ, φ)-transform relies on maximizing the dual functional with respect to one variable while holding the other fixed. This transform was originally introduced in [11] for superlinear φ and subsequently generalized in [12] to functions φ of Legendre type, assuming compact spaces X, Y and a Lipschitz cost c. The following analysis extends the results of [12] to Polish spaces and bounded cost functions. The introduction of the (c, λ, φ)-transform is critical for the study of the dual problem (31). Fundamentally, given f ∈ L∞(X, P_X) and g ∈ L∞(Y, P_Y), the corresponding (c, λ, φ)-transforms, denoted f(c,λ,φ) and g(c,λ,φ), are designed to satisfy

$$f^{(c,\lambda,\phi)} \in \arg \max_{g \in L^{\infty}(Y, P_Y)} \mathcal{D}(f,g); \text{ and} \quad (48)$$

$$g^{(c,\lambda,\phi)} \in \arg \max_{f \in L^{\infty}(X, P_X)} \mathcal{D}(f,g), \quad (49)$$

where the dual functional D is in (44). In the subsequent lemmas, properties are stated for the transform over L∞(X, P_X), however analogous properties hold for the transform over L∞(Y, P_Y). The following lemma formally defines this transform.

**Lemma 2.6.** Given a function $f \in L^\infty(X, P_X)$, its (c, λ, φ)-transform, denoted by $f^{(c,\lambda,\phi)}$, is implicitly defined at $y \in Y$ by

$$\int_X \phi^{*'} \left( \frac{f(x) + f^{(c,\lambda,\phi)}(y) - c(x,y)}{\lambda} \right) dP_X(x) = 1, \quad (50a)$$

where $\phi^*$ is the convex conjugate of $\phi$ (Definition 1.2), provided such a number $f^{(c,\lambda,\phi)}(y)$ exists. Otherwise, the value $f^{(c,\lambda,\phi)}(y)$ is explicitly defined by

$$f^{(c,\lambda,\phi)}(y) = \inf_{x \in X} \left[ \lambda \beta_{\phi} + c(x,y) - f(x) \right], \quad (50b)$$

where $\beta_\phi$ is defined in (7). Moreover, the function $f^{(c,\lambda,\phi)}$ satisfies (48).

*Proof:* Let $g \in L^\infty(Y, P_Y)$ be an integrable function. The dual functional in (44) can be rewritten as

$$\mathcal{D}(f,g) = \int_X f(x) dP_X(x) + \int_Y \Phi_y(g(y)) dP_Y(y), \quad (51)$$

where $\Phi_y$ is defined by

$$\Phi_y(\tau) = \tau - \lambda \int_X \phi^* \left( \frac{f(x) + \tau - c(x,y)}{\lambda} \right) dP_X(x), \quad (52)$$
where $\phi^*$ is the convex conjugate (Definition 1.2) of the function $\phi$. The property required for the $(c, \lambda, \phi)$-transform is stated in (48). From (51), a function $g$ maximizes the functional $D(f, \cdot)$ if and only if it maximizes the functional $h \mapsto \int_Y \Phi_y(h(y))dP_Y(y)$. According to [30, Theorem 14.60], $g$ maximizes the functional $h \mapsto \int_Y \Phi_y(h(y))dP_Y(y)$ if and only if, for all $y$, $g(y)$ maximizes $\Phi_y$ in (52). Consider an element $y$ in $\mathcal{Y}$, and let the function $\bar{g}$ be defined at $y$ by

$$ \bar{g}(y) \triangleq \inf_{x \in \mathcal{X}} [\lambda\beta_{\phi} + c(x, y) - f(x)], \quad (53) $$

where $\beta_{\phi}$ is defined in (7). Note that for all $\tau < \bar{g}(y)$, and for all $x \in \mathcal{X}$, the inequality $\frac{f(x)+\tau-c(x,y)}{\lambda} < \beta_{\phi}$ holds. Consequently, the following quantity is well defined

$$ S_y(\tau) \triangleq \int_{\mathcal{X}} \phi^{*'} \left( \frac{f(x) + \tau - c(x, y)}{\lambda} \right) dP_X(x), \quad (54) $$

where $\phi^{*'}$ is the derivative of the convex conjugate of $\phi$. As $\phi$ is a Legendre-type generator (Definition 1.5), Lemma 1.1 ensures that $\phi^*$ is differentiable on $(-\infty, \beta_{\phi})$; therefore, $S_y$ is defined on the open interval $(-\infty, \bar{g}(y))$. Moreover, strict convexity of $\phi^*$ implies strict monotonicity of $\phi^{*'}$, which in turn renders $S_y$ strictly increasing. Continuity of $S_y$ follows from the dominated convergence theorem [25, Theorem 1.34], as the monotonicity of $\phi^{*'}$ ensures that the integrand is locally dominated by an integrable function. Pointwise convergence at $-\infty$ is obtained from Lemma 1.1 as follows

$$ \phi^{*'} \left( \frac{f(x) + \tau - c(x, y)}{\lambda} \right) \underset{\tau \to -\infty}{\longrightarrow} 0. \quad (55) $$

As the function $\phi^{*'}$ is strictly increasing, a number $\tau_0$ satisfying $\frac{\tau_0 + \|f-c\|_{\infty}}{\lambda} < \beta_{\phi}$ can be chosen such that, for all $x$ in $\mathcal{X}$, the following holds

$$ 0 \le \phi^{*'} \left( \frac{f(x) + \tau - c(x, y)}{\lambda} \right) \le \phi^{*'} \left( \frac{\tau_0 + \|f - c\|_{\infty}}{\lambda} \right). \quad (56) $$

The upper bound obtained is thus a finite constant, acting as an integrable dominating function. Dominated convergence [25, Theorem 1.34] applied to (55) implies that $S_y(\tau) \underset{\tau \to -\infty}{\longrightarrow} 0$. Two cases arise:

A. $\lim_{\tau \uparrow \bar{g}(y)} S_y(\tau) > 1$. In this scenario, since $S_y$ is continuous and strictly increasing, the intermediate value theorem guarantees the existence of a unique $\tau_y < \bar{g}(y)$ such that $S_y(\tau_y) = 1$. This is equivalent to $\Phi'_y(\tau_y) = 0$. Concavity of $\Phi_y$ implies that the maximum is attained at $\tau_y$; hence $f^{(c, \lambda, \phi)}(y) = \tau_y$ satisfies the implicit equation (50a).

B. $\lim_{\tau \uparrow \bar{g}(y)} S_y(\tau) \le 1$. In this case, for all $\tau < \bar{g}(y)$, $\Phi'_y(\tau) \ge 0$. Maximizing $\Phi_y$ therefore requires choosing the largest feasible value, namely $\bar{g}(y)$; it remains to verify that this value is feasible. By definition of $\bar{g}(y)$ in (53), for all $x \in \mathcal{X}$, $\frac{f(x)+\bar{g}(y)-c(x,y)}{\lambda} \le \beta_{\phi}$. Assume, for the sake of contradiction,
that there exists a subset $\Omega \subset \mathcal{X}$ of strictly positive $P_X$-measure such that $\frac{f(x)+\bar{g}(y)-c(x,y)}{\lambda} = \beta_{\phi}$ on $\Omega$. Let $(\tau_n)$ be a sequence in $(-\infty, \bar{g}(y))$ such that $\tau_n \uparrow \bar{g}(y)$. As the function $\phi^{*'}$ is positive and strictly increasing, it holds that

$$S_y(\tau_n) \geq \int_{\Omega} \phi^{*'} \left( \frac{f(x) + \tau_n - c(x, y)}{\lambda} \right) dP_X(x). \quad (57)$$

For all $x \in \Omega$, since $\frac{f(x)+\tau_n-c(x,y)}{\lambda} \xrightarrow{n \to +\infty} \beta_{\phi}$, and the function $\phi$ is a Legendre-type generator (Definition 1.5), Lemma 1.1 implies that

$$\phi^{*'} \left( \frac{f(x) + \tau_n - c(x, y)}{\lambda} \right) \xrightarrow{n \to +\infty} +\infty. \quad (58)$$

Consequently, the monotone convergence theorem [25, Theorem 1.26] implies

$$\int_{\Omega} \phi^{*'} \left( \frac{f(x) + \tau_n - c(x, y)}{\lambda} \right) dP_X(x) \xrightarrow{n \to +\infty} +\infty. \quad (59)$$

Combined with (57), (59) leads to $S_y(\tau_n) \to +\infty$, which contradicts the assumption $\lim_{\tau \uparrow \bar{g}(y)} S_y(\tau) \leq 1$. Therefore:

$$\frac{f(x) + \bar{g}(y) - c(x, y)}{\lambda} <_{P_X}^{\text{a.s.}} \beta_{\phi}. \quad (60)$$

Thus $f^{(c,\lambda,\phi)}(y) \triangleq \bar{g}(y)$ is a feasible value, and (50b) is satisfied.

Note that in the case of a superlinear $\phi$ (i.e. $\beta_{\phi} = +\infty$), the $(c, \lambda, \phi)$-transform is defined by (50a) for all $y \in \mathcal{Y}$, since $\bar{g}(y) = +\infty$, and monotone convergence [25, Theorem 1.26] yields $\lim_{\tau \uparrow +\infty} S_y(\tau) = +\infty$.

The following lemma establishes a critical bound on the range of the $(c, \lambda, \phi)$-transform. This bound ensures that the transform remains within the domain where the conjugate function $\phi^*$ behaves regularly, a property essential for subsequent differentiability arguments.

**Lemma 2.7.** Consider a bounded function $f \in L^{\infty}(\mathcal{X}, P_X)$. Then

$$\frac{f \oplus f^{(c, \lambda, \phi)} - c}{\lambda} <_{P_X P_Y}^{\text{a.s.}} \beta_{\phi}, \quad (61)$$

where $f^{(c,\lambda,\phi)}$ is the $(c, \lambda, \phi)$-transform of $f$, defined in (50).

*Proof:* The proof of Lemma 2.6 demonstrates that for a fixed $y \in \mathcal{Y}$, $f^{(c,\lambda,\phi)}(y)$ satisfies either (50a) or (50b). The proof proceeds by analyzing these two cases. In the first case, there exists $\tau_y < \bar{g}(y)$ solving (50a), where $\bar{g}(y)$ is defined in (53). Consequently,

$$\frac{f(x) + \tau_y - c(x, y)}{\lambda} < \beta_{\phi}, \quad (62)$$
implying that the strict inequality holds for all $x \in \mathcal{X}$.

In the second case, the $(c, \lambda, \phi)$-transform of $f$ at $y$ satisfies $f^{(c, \lambda, \phi)}(y) = \bar{g}(y)$. The argument by contradiction presented in the proof of item B of Lemma 2.6 establishes (60), which is equivalent to

$$P_X(\{x \in \mathcal{X} : \frac{f(x)+\bar{g}(y)-c(x,y)}{\lambda} = \beta_{\phi}\}) = 0. \qquad (63)$$

Hence (61) is satisfied, and for all $y \in \mathcal{Y}$, in both cases A and B, the equality in (63) is satisfied. Define the set $\mathcal{E}$ by

$$\mathcal{E} \triangleq \{(x, y) \in \mathcal{X} \times \mathcal{Y} : \frac{f(x) + f^{(c, \lambda, \phi)}(y) - c(x, y)}{\lambda} = \beta_{\phi}\}. \qquad (64)$$

According to the Fubini-Tonelli Theorem [25, Theorem 8.8] and (63), the following property is satisfied

$$P_X P_Y(\mathcal{E}) = \int_{\mathcal{Y}} P_X(\{x : (x, y) \in \mathcal{E}\}) dP_Y(y) = \int_{\mathcal{Y}} 0 dP_Y(y) = 0. \qquad (65)$$

Therefore, the set $\mathcal{E}$ has zero $P_X P_Y$-measure, or equivalently

$$\frac{f(x) + f^{(c, \lambda, \phi)}(y) - c(x, y)}{\lambda} \underset{P_X P_Y}{\overset{\text{a.s.}}{<}} \beta_{\phi}. \qquad (66)$$

This completes the proof.

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

The following lemma plays a central role in establishing the existence of maximizing potentials. Statement (i) generalizes [11, Lemma 3.1(i)], and its proof follows a similar strategy. Statement (ii) is adapted from [12, Proposition 1.15(iv)]; while the assumptions on $\mathcal{X}$, $\mathcal{Y}$, and c differ slightly, the core argument remains applicable.

**Lemma 2.8.** Given a function $f \in L^\infty(\mathcal{X}, P_X)$, its $(c, \lambda, \phi)$-transform defined in (50) satisfies the following inequalities:

$$\mathrm{osc}(f^{(c, \lambda, \phi)}) \le 2\|c\|_{\infty}; \quad \text{and} \qquad (67)$$

$$\|f^{(c, \lambda, \phi)}\|_{\infty} \le \begin{cases} \|f\|_{\infty} + \|c\|_{\infty} & \text{if } \beta_{\phi} = +\infty, \\ \|f\|_{\infty} + \|c\|_{\infty} + \beta_{\phi} & \text{if } \beta_{\phi} < +\infty. \end{cases} \qquad (68)$$

*Proof:* The proof of (67) is done first in three parts, and then the proof of (68) is presented. Let $y$ and $z$ be two elements in the set $\mathcal{Y}$. The proof of the inequality in (67) is divided into three parts; the first part assumes that both $f^{(c, \lambda, \phi)}(y)$ and $f^{(c, \lambda, \phi)}(z)$ are defined by the implicit equation (50a). The second part assumes that both $f^{(c, \lambda, \phi)}(y)$ and $f^{(c, \lambda, \phi)}(z)$ are defined by the explicit equation (50b), and the third part assumes that $f^{(c, \lambda, \phi)}(y)$ is defined by (50a) and $f^{(c, \lambda, \phi)}(z)$ is defined by (50b). Since c is bounded, for all $x \in \mathcal{X}$, the inequality $c(x, z) \le c(x, y) + 2\|c\|_{\infty}$ holds, which implies

$$f(x) + f^{(c, \lambda, \phi)}(z) - c(x, z) \ge f(x) + f^{(c, \lambda, \phi)}(z) - c(x, y) - 2\|c\|_{\infty}. \qquad (69)$$
As the function $\phi^{*'}$ is strictly increasing (Lemma 1.1), the inequality in (69) yields

$$1 = S_y \left( f^{(c,\lambda,\phi)}(y) \right) \qquad (70)$$

$$= S_z \left( f^{(c,\lambda,\phi)}(z) \right) \qquad (71)$$

$$= \int_X \phi^{*'} \left( \frac{f(x) + f^{(c,\lambda,\phi)}(z) - c(x, z)}{\lambda} \right) dP_X(x) \qquad (72)$$

$$\geq \int_X \phi^{*'} \left( \frac{f(x) + f^{(c,\lambda,\phi)}(z) - c(x, y) - 2\|c\|_{\infty}}{\lambda} \right) dP_X(x) \qquad (73)$$

$$\geq S_y \left( f^{(c,\lambda,\phi)}(z) - 2\|c\|_{\infty} \right), \qquad (74)$$

where the equality in (72) follows from the definition of $S_z$ in (54), and the inequality in (73) follows from the monotonicity of $\phi^{*'}$ combined with (69). Since $\phi^{*'}$ is monotone, $S_y$ is a non-decreasing function; the previous inequality therefore implies that $f^{(c,\lambda,\phi)}(y) \geq f^{(c,\lambda,\phi)}(z) - 2\|c\|_{\infty}$. Exchanging the role of $y$ and $z$ leads to

$$|f^{(c,\lambda,\phi)}(y) - f^{(c,\lambda,\phi)}(z)| \leq 2\|c\|_{\infty}, \qquad (75)$$

which yields (67). The second part of the proof of the inequality in (67) is as follows. Consider the case in which the two values are defined by (50b). Let $\epsilon$ be a strictly positive real; from (75), the definition of the infimum implies the existence of some $x_{\epsilon} \in \mathcal{X}$ such that

$$f^{(c,\lambda,\phi)}(z) \geq \lambda\beta_{\phi} + c(x_{\epsilon}, z) - f(x_{\epsilon}) - \epsilon. \qquad (76)$$

Furthermore, from Lemma 2.7, $f^{(c,\lambda,\phi)}(y)$ verifies

$$f^{(c,\lambda,\phi)}(y) \leq \lambda\beta_{\phi} + c(x_{\epsilon}, y) - f(x_{\epsilon}). \qquad (77)$$

Combining both inequalities in (76) and in (77) yields

$$f^{(c,\lambda,\phi)}(y) - f^{(c,\lambda,\phi)}(z) \leq c(x_{\epsilon}, y) - c(x_{\epsilon}, z) + \epsilon \qquad (78)$$

$$\leq 2\|c\|_{\infty} + \epsilon. \qquad (79)$$

Exchanging the roles of $y$ and $z$ yields $|f^{(c,\lambda,\phi)}(y) - f^{(c,\lambda,\phi)}(z)| \leq 2\|c\|_{\infty} + \epsilon$. This holds for all $\epsilon$, thus the equality in (75) is satisfied.

The third part of the proof of the inequality in (67) is as follows. Consider the case in which $f^{(c,\lambda,\phi)}(y)$ is defined by (50a) and $f^{(c,\lambda,\phi)}(z)$ is defined by (50b). Let $\epsilon$ be strictly positive; from Lemma 2.7, the quantity $f^{(c,\lambda,\phi)}(y)$ satisfies $f^{(c,\lambda,\phi)}(y) \leq \lambda\beta_{\phi} + c(x_{\epsilon}, y) - f(x_{\epsilon})$. Consequently, the same reasoning as above leads to the inequality in (79), and since this holds for all $\epsilon > 0$,

$$f^{(c,\lambda,\phi)}(y) - f^{(c,\lambda,\phi)}(z) \leq 2\|c\|_{\infty}. \qquad (80)$$
In the proof of Lemma 2.6, for the case in which $f^{(c,\lambda,\phi)}(z)$ is defined by (50b), the corresponding integral functional does not saturate the constraint, meaning that $S_z(f^{(c,\lambda,\phi)}(z)) < 1$. Combining this strict inequality with the definition of the implicit case for $y$, the following relations hold

$$1 = S_y(f^{(c,\lambda,\phi)}(y)) \qquad (81)$$

$$> S_z(f^{(c,\lambda,\phi)}(z)) \qquad (82)$$

$$= \int_X \phi^{*'} \left( \frac{f(x) + f^{(c,\lambda,\phi)}(z) - c(x, z)}{\lambda} \right) dP_X(x) \qquad (83)$$

$$\geq \int_X \phi^{*'} \left( \frac{f(x) + f^{(c,\lambda,\phi)}(z) - c(x, y) - 2\|c\|_{\infty}}{\lambda} \right) dP_X(x) \qquad (84)$$

$$= S_y(f^{(c,\lambda,\phi)}(z) - 2\|c\|_{\infty}). \qquad (85)$$

where the equality in (83) is the definition of $S_z$ in (54) and the inequality in (84) is justified by the boundedness of the cost function, specifically $c(x, z) \le c(x, y) + 2\|c\|_{\infty}$, which implies $-c(x, z) \ge -c(x, y) - 2\|c\|_{\infty}$, coupled with the monotonicity of $\phi^{*'}$. Since $S_y$ is strictly increasing, the inequality between the first term and the last term in (85) allows to conclude that $f^{(c,\lambda,\phi)}(y) \ge f^{(c,\lambda,\phi)}(z) - 2\|c\|_{\infty}$. Combining this with (80) yields the equality in (75), which concludes the proof of (67). The proof of (68) is identical to the proof of (iv) in [12, Proposition 15]. ■

The following Lemma serves as a crucial step in establishing the existence of optimal potentials for the dual problem in (31). In variational analysis, a standard proof strategy involves considering a maximizing sequence of potentials. However, ensuring the weak* ([31, Section 3.4]) convergence of a subsequence of such a sequence requires it to be bounded. This lemma provides the necessary tools to overcome this difficulty. It demonstrates that the search for the supremum can be restricted, without loss of optimality, to a bounded subset of the functional space. Specifically, any arbitrary pair of potentials can be improved, or at least matched, by a pair of potentials whose $L^\infty$ norms are uniformly bounded by a constant depending solely on the problem parameters. This result parallels [11, Lemma 3.3].

**Lemma 2.9.** Given two bounded functions $f \in L^\infty(\mathcal{X}, P_X)$ and $g \in L^\infty(\mathcal{Y}, P_Y)$, there exist two functions $\hat{f} \in L^\infty(\mathcal{X}, P_X)$ and $\hat{g} \in L^\infty(\mathcal{Y}, P_Y)$ that satisfy

$$D(\hat{f}, \hat{g}) \geq D(f, g), \qquad (86)$$

where $D$ is in (44). Moreover such functions $\hat{f}$ and $\hat{g}$ satisfy

$$\|\hat{f}\|_{\infty} \leq M; \text{ and} \qquad (87)$$

$$\|\hat{g}\|_{\infty} \leq M, \qquad (88)$$

where $M$ exclusively depends on the functions $c$ and $\phi$, and the real number $\lambda$.
*Proof:* From the definition of the $(c, \lambda, \phi)$-transform (48), it follows that

$$D(f, g) \leq D(f, f^{(c, \lambda, \phi)}), \qquad (89)$$

and Lemma 2.8 (i) yields $\text{osc}(f^{(c,\lambda,\phi)}) \leq 2\|c\|_{\infty}$. Hence, there exists a constant $a \in \mathbb{R}$ such that $\hat{g} \triangleq f^{(c,\lambda,\phi)} - a$ verifies

$$\|\hat{g}\|_{\infty} \leq 2\|c\|_{\infty}. \qquad (90)$$

Then, $D(f, f^{(c,\lambda,\phi)}) = D(f + a, \hat{g})$. Let $\hat{f}$ be the $(c, \lambda, \phi)$-transform of $\hat{g}$, that is $\hat{f} \triangleq \hat{g}^{(c,\lambda,\phi)}$. The definition of the $(c, \lambda, \phi)$-transform implies

$$D(f, g) \leq D(f, f^{(c, \lambda, \phi)}) = D(f + a, \hat{g}) \leq D(\hat{f}, \hat{g}), \qquad (91)$$

which completes the proof of (86). The two inequalities in (68) and (90) yield

$$\|\hat{f}\|_{\infty} \leq \|\hat{g}\|_{\infty} + \|c\|_{\infty} + K_{\phi}; \qquad (92)$$

$$\leq 3\|c\|_{\infty} + K_{\phi}, \qquad (93)$$

where $K_{\phi} = 0$ if $\beta_{\phi} = +\infty$ and $K_{\phi} = \beta_{\phi}$ otherwise. Finally, setting $M = 3\|c\|_{\infty} + K_{\phi}$ in (93) yields (87) and (88), which concludes the proof. ■

### 2.4.2 Existence of maximizing potentials

While the primal problem in (22) may not always admit a minimizer without further assumptions, specifically regarding the behavior of measures at infinity, optimal potentials can always be found for the dual problem in (31). This asymmetry highlights the robustness of the dual formulation. Relying on the bounds established in Lemma 2.9, the search for optimal potentials is effectively restricted to a bounded set, allowing for the application of compactness arguments. Furthermore, this subsection presents the first-order sufficient conditions for optimality, providing a practical criterion to identify these solutions. The following theorem asserts the existence of maximizers for the dual problem in (31).

**Theorem 2.1.** The supremum in problem (31) is attained for a unique pair of potentials $(f_0, g_0)$ up to the transformation $(f, g) \mapsto (f + a, g - a)$, with $a \in \mathbb{R}$.

*Proof:* Using Lemma 2.5, and Lemma 2.9, the proof of this result follows the same reasoning as the proof of [11, Theorem 3.4]. ■

The uniqueness is stated "up to a constant" because the dual functional depends on the potentials primarily through their sum $f(x)+g(y)$. Consequently, adding a scalar $a$ to $f$ and subtracting it from $g$ leaves the objective value unchanged. A fundamental property of optimal potentials arises from the structure of the dual maximization problem. Since the $(c, \lambda, \phi)$-transform acts as a coordinate-wise maximization operator, updating one potential to its optimal value given the other, it follows that at any global optimum, neither potential can be improved individually. Therefore, optimal potentials are necessarily fixed points of these
alternating updates. This interdependency provides a strong characterization of the solutions to the dual problem, reducing the search for a pair of functions to the search for a single coupled state.

**Lemma 2.10.** If two functions $f \in L^\infty(X, P_X)$ and $g \in L^\infty(Y, P_Y)$ are maximizing potentials for the dual problem in (31), then

$$f = g^{(c,\lambda,\phi)}; \text{ and} \qquad (94)$$

$$g = f^{(c,\lambda,\phi)}, \qquad (95)$$

where $g^{(c,\lambda,\phi)}$ and $f^{(c,\lambda,\phi)}$ are the corresponding $(c, \lambda, \phi)$-transform of $g$ and $f$.

*Proof:* This is a direct consequence of the definition of the $(c, \lambda, \phi)$-transform in (48) and (49). ■

Lemma 2.10 provides a necessary condition for optimality. However, sufficiency is established when the potentials belong to a specific regularity set defined below. The construction of this set is not arbitrary; it is specifically tailored to ensure that the dual functional in (32) is Gâteaux differentiable. To establish the first-order optimality conditions for the dual problem in (31), the integral term involving $\phi^*$ in the dual functional in (32) has to be differentiated. The interchange of differentiation and integration via the dominated convergence theorem requires the derivative of the integrand to be bounded, in particular the function $\phi^{*'}$ needs to be bounded. However, since $|\phi^{*'}(u)|$ diverges as $u$ approaches the boundary $\beta_\phi$ (Lemma 1.1), differentiation is only guaranteed if the argument of $\phi^*$ stays strictly bounded away from this singularity. The set $A_{\phi,c}$ enforces precisely this safety margin, prohibiting the potentials from reaching the boundary where the derivative cannot be bounded. Moreover, within this set, the optimality condition is consistently expressed through the implicit equation of the $(c, \lambda, \phi)$-transform. The set is defined as follows

$$A_{\phi,c} \triangleq \{(f,g) \in L^\infty(X, P_X) \times L^\infty(Y, P_Y): \\ \exists \gamma < \beta_\phi, \text{Ran} \left( \frac{f \oplus g - c}{\lambda} \right) \underset{P_X P_Y}{\overset{a.s.}{\subset}} (-\infty, \gamma)\}, \quad (96)$$

where $\text{Ran}(f \oplus g - c)$ denotes the range of the function $f \oplus g - c$. The following lemma establishes that for potentials within $A_{\phi,c}$, the maximization condition is equivalent to satisfying specific marginal constraints involving the derivative of $\phi^*$.

**Lemma 2.11.** A pair of functions $(\hat{f}, \hat{g}) \in A_{\phi,c}$, where $A_{\phi,c}$ is defined in (96), are maximizing potentials for the dual problem in (31) if and only if:

$$\int_Y \phi'^{-1} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) dP_Y(y) \underset{P_X}{\overset{a.s.}{\equiv}} 1; \text{ and} \qquad (97)$$

$$\int_X \phi'^{-1} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) dP_X(x) \underset{P_Y}{\overset{a.s.}{\equiv}} 1. \qquad (98)$$
*Proof:* The first step consists in proving that the dual functional D in (44) is Gâteaux differentiable on $A_{\phi,c}$. The two first terms of D in (44) are linear, thus the analysis focuses on the third term

$$J(f,g) = \int_{X \times Y} \phi^* \left( \frac{f(x) + g(y) - c(x, y)}{\lambda} \right) dP_X P_Y(x, y). \quad (99)$$

Let $(a, b)$ be a pair of functions in $L^{\infty}(X, P_X) \times L^{\infty}(Y, P_Y)$, $(f, g)$ be a pair of functions in $A_{\phi,c}$ and $t \in \mathbb{R}$ a real number. The variation of $J$ at $(f, g)$ in the direction $(a, b)$ is given by

$$\frac{J(f + ta, g + tb) - J(f, g)}{t} = \int_{X \times Y} \frac{1}{t} \left( \phi^* \left( \frac{f(x) + g(y) + t(a(x) + b(y)) - c(x, y)}{\lambda} \right) - \phi^* \left( \frac{f(x) + g(y) - c(x, y)}{\lambda} \right) \right) dP_X P_Y(x, y), \quad (100)$$

when that quantity is well defined. Since $(f, g) \in A_{\phi,c}$, definition (96) implies the existence of $\gamma < \beta_{\phi}$ such that the range of $(f \oplus g - c)/\lambda$ is almost surely contained in $(-\infty, \gamma]$. Let $\delta = \beta_{\phi} - \gamma > 0$. Since $a$ and $b$ are bounded functions, the norm $\|a \oplus b\|_{\infty} \le \|a\|_{\infty} + \|b\|_{\infty}$ is finite. Choosing $t_0 > 0$ such that $t_0(\|a\|_{\infty} + \|b\|_{\infty})/\lambda < \delta/2$ ensures that for all $|t| \le t_0$, the perturbed argument

$$\frac{f(x) + g(y) - t_0(\|a\|_{\infty} + \|b\|_{\infty}) - c(x, y)}{\lambda} \quad (101)$$

remains bounded away from $\beta_{\phi}$. Specifically, the set

$$B_{t_0} \triangleq \left[ \inf_{x,y \in X \times Y} \frac{f(x) + g(y) - t_0(\|a\|_{\infty} + \|b\|_{\infty}) - c(x, y)}{\lambda}, \sup_{x,y \in X \times Y} \frac{f(x) + g(y) + t_0(\|a\|_{\infty} + \|b\|_{\infty}) - c(x, y)}{\lambda} \right] \quad (102)$$

satisfies the following inclusions

$$B_{t_0} \subset (-\infty, \gamma + \delta/2] \quad (103)$$

$$\subset (-\infty, \beta_{\phi} - \delta/2). \quad (104)$$

For all $|t| \le t_0$, clearly $B_{|t|} \subset B_{t_0}$, where $B_{|t|}$ is defined similarly as $B_{t_0}$ in (102), by replacing $t_0$ with $|t|$. Note that for all $t \in \mathbb{R}$ such that $|t| \le t_0$, the following inclusions hold,

$$\mathrm{Ran} \left( \frac{f \oplus g + t(a \oplus b) - c}{\lambda} \right) \subset B_{|t|} \subset B_{t_0} \subset (-\infty, \beta_{\phi} - \delta/2). \quad (105)$$

Within the interval $(-\infty, \beta_{\phi})$, the function $\phi^*$ is $C^1$ (see Section 1.2). Consequently, according to (105), the expression in (100) is well defined for all $|t| \le t_0$,
and the following pointwise convergence is obtained

$$ \frac{1}{t} \left[ \phi^* \left( \frac{f(x) + g(y) + t(a(x) + b(y)) - c(x, y)}{\lambda} \right) - \phi^* \left( \frac{f(x) + g(y) - c(x, y)}{\lambda} \right) \right] \\ \xrightarrow{t \to 0} (a(x) + b(y)) \phi^{*\prime} \left( \frac{f(x) + g(y) - c(x, y)}{\lambda} \right). \quad (106) $$

Given that $a, b$ and $c$ are bounded functions, the set $B_{t_0}$ is a bounded subset of $(-\infty, \beta_\phi - \delta/2)$. Since the function $\phi^{*\prime}$ is continuous on $(-\infty, \beta_\phi - \delta/2]$ and $B_{t_0}$ is bounded, there exists some $M$ such that $\sup_{B_{t_0}} |\phi^{*\prime}| \le M$. The mean value theorem [32, Corollary 5.11] yields

$$ \frac{1}{t} \left| \phi^* \left( \frac{f(x) + g(y) + t(a(x) + b(y)) - c(x, y)}{\lambda} \right) - \phi^* \left( \frac{f(x) + g(y) - c(x, y)}{\lambda} \right) \right| \\ \le M \left| \frac{1}{t} \frac{t(a(x) + b(y))}{\lambda} \right| \quad (107) \\ \le M \frac{\|a\|_{\infty} + \|b\|_{\infty}}{\lambda}. \quad (108) $$

Since the bound in (108) is finite, the dominated convergence theorem [25, Theorem 1.34] yields

$$ \frac{J((f,g) + t(a,b)) - J(f,g)}{t} \xrightarrow{t \to 0} \int_{X \times Y} \left( \frac{a(x) + b(y)}{\lambda} \right) \phi^{-1} \left( \frac{f(x) + g(y) - c(x,y)}{\lambda} \right) dP_X P_Y(x,y). \quad (109) $$

The Gâteaux differential of $D$ is obtained by combining the differential of the linear terms with the differential of $J$ derived above. Hence, the Gâteaux differential of $D$ at $(f, g)$ in the direction $(a, b)$ is

$$ \langle \partial D(f, g), (a, b) \rangle = \int_X a(x) dP_X(x) + \int_Y b(y) dP_Y(y) \\ - \int_{X \times Y} (a(x) + b(y)) \phi^{-1} \left( \frac{f(x) + g(y) - c(x, y)}{\lambda} \right) dP_X P_Y(x, y), \quad (110) $$

or, equivalently

$$ \langle \partial D(f, g), (a, b) \rangle = \int_X a(x) \left( 1 - \int_Y \phi^{*\prime} \left( \frac{f(x) + g(y) - c(x, y)}{\lambda} \right) dP_Y(y) \right) dP_X(x) \\ + \int_Y b(y) \left( 1 - \int_X \phi^{-1} \left( \frac{f(x) + g(y) - c(x, y)}{\lambda} \right) dP_X(x) \right) dP_Y(y). \quad (111) $$

Suppose that $(\hat{f}, \hat{g}) \in A_{\phi,c}$ satisfies (97) and (98). The equality in (111) directly implies that for all $(a, b) \in L^\infty(X, P_X) \times L^\infty(Y, P_Y)$, the equation
$\langle \partial D(\hat{f}, \hat{g}), (a, b) \rangle = 0$ is satisfied. Since D is concave, [18, Chapter 2.4] implies that D achieves its global maximum at $(\hat{f}, \hat{g})$.

Conversely, suppose D achieves a global maximum at $(\hat{f}, \hat{g})$. From [33, Theorem 1, Chapter 7.4], it holds that for all $(a, b)$, $\langle \partial D(\hat{f}, \hat{g}), (a, b) \rangle = 0$. Setting $b$ to 0 and then $a$ to 0 simplifies this condition to

$$ \int_X a(x) \left( 1 - \int_Y \phi^{*'} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) dP_Y(y) \right) dP_X(x) = 0; \quad (112) $$

$$ \int_Y b(y) \left( 1 - \int_X \phi^{*'} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) dP_X(x) \right) dP_Y(y) = 0, \quad (113) $$

for all $a \in L^\infty(X, P_X)$ and all $b \in L^\infty(Y, P_Y)$. This is equivalent to $(\hat{f}, \hat{g})$ satisfying (97) and (98). ■

## 2.5 Attainment in the primal problem

As demonstrated in [12], an optimal solution to the problem in (13) may fail to exist unless Definition 1.6 is modified to allow measures with a non-zero singular part with respect to the product measure $P_X P_Y$ to yield a finite $\phi$-divergence. In such a case, the optimal coupling may exhibit a non-zero singular part itself. Such a solution presents practical challenges, primarily because the singular component lacks a closed-form expression, and the Sinkhorn algorithm is not available, as it has been shown in [12, Theorem 3]. Ensuring the existence of a computable solution therefore necessitates additional assumptions.

### 2.5.1 Closed-form expression of the optimal coupling

A sufficient condition for the existence of a minimizer for problem (22) is that the optimal potentials of the dual problem (31) belong to the set $A_{\phi,c}$, defined in (96).

**Theorem 2.2.** Suppose that there exists a pair of optimal potentials $(\hat{f}, \hat{g})$ in $A_{\phi,c}$ for the dual problem in (31). Then, strong duality holds between the optimization problems (22) and (31), and the unique minimizer of the primal problem is given by

$$ P_{XY}^{(c, \phi)}(x, y) \triangleq \phi^{*'} \left( \frac{\hat{f} \oplus \hat{g} - c}{\lambda} \right). \quad (114) $$

*Proof:* This proof is adapted from [11, Proposition 3.5]. Given two real numbers $s$ and $t$, and $\lambda$ the regularization parameter set at the beginning of Section 2, the Fenchel-Young inequality, see [18, Theorem 2.3.1], applied to the function $\phi$ implies

$$ \lambda\phi(t) - ts \geq -\lambda\phi^*\left(\frac{s}{\lambda}\right), \quad (115) $$
with equality if and only if $t \in \partial\phi^*(\frac{s}{\lambda})$, where $\partial\phi^*(\frac{s}{\lambda})$ denotes the subgradient of $\phi^*$ at $\frac{s}{\lambda}$. Let $(\hat{f}, \hat{g})$ in $A_{\phi,c}$ be optimal potentials, and fix $(x, y) \in \mathcal{X} \times \mathcal{Y}$. Define the function $r^{(\phi,c)} \triangleq \phi^*(\frac{\hat{f} \oplus \hat{g} - c}{\lambda})$, and the reals $t = r^{(\phi,c)}(x, y)$ and $s = \hat{f}(x) + \hat{g}(y) - c(x, y)$. Since the derivative of a convex function is an element of its subdifferential, it follows that $t \in \partial\phi^*(\frac{s}{\lambda})$; thus, (115) holds with equality, which can be written as

$$ \lambda\phi(r^{(\phi,c)}(x, y)) - r^{(\phi,c)}(x, y)(\hat{f}(x) + \hat{g}(y) - c(x, y)) = -\lambda\phi^*\left(\frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda}\right). \quad (116) $$

Rearranging terms yields

$$ \lambda\phi(r^{(\phi,c)}(x, y)) + c(x, y)r^{(\phi,c)}(x, y) = r^{(\phi,c)}(x, y)\hat{f}(x) + r^{(\phi,c)}(x, y)\hat{g}(y) - \lambda\phi^*\left(\frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda}\right). \quad (117) $$

Integrating both sides of (117) with respect to $P_X P_Y$ leads to

$$ \begin{aligned} & \int_{\mathcal{X} \times \mathcal{Y}} c(x, y) r^{(\phi,c)}(x, y) dP_X P_Y(x, y) + \int_{\mathcal{X} \times \mathcal{Y}} \phi(r^{(\phi,c)}(x, y)) dP_X P_Y(x, y) \\ &= \int_{\mathcal{X} \times \mathcal{Y}} r^{(\phi,c)}(x, y) \hat{f}(x) dP_X P_Y(x, y) + \int_{\mathcal{X} \times \mathcal{Y}} r^{(\phi,c)}(x, y) \hat{g}(y) dP_X P_Y(x, y) \\ & \quad - \lambda \int_{\mathcal{X} \times \mathcal{Y}} \phi^* \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) dP_X P_Y(x, y). \end{aligned} \quad (118) $$

Since $(\hat{f}, \hat{g})$ in $A_{\phi,c}$ are optimal potentials, Lemma 2.11 ensures that $\hat{f}$ satisfies (97) and $\hat{g}$ satisfies (98). This confirms that $r^{(\phi,c)}$ satisfies the constraints of (13), namely $G(A(r^{(\phi,c)})) = 0$, where $G$ is in (20). Since $r^{(\phi,c)}$ satisfies the marginal constraints of the primal problem, specifically $\int_{\mathcal{Y}} r^{(\phi,c)}(\cdot, y)dP_Y(y) = 1$ $P_X$-almost surely, applying Fubini's theorem [25, Theorem 8.8] to the integrable function $(x, y) \mapsto r^{(\phi,c)}(x, y)\hat{f}(x)$ yields the equality

$$ \begin{aligned} \int_{\mathcal{X} \times \mathcal{Y}} r^{(\phi,c)}(x, y) \hat{f}(x) dP_X P_Y(x, y) &= \int_{\mathcal{X}} \hat{f}(x) \left( \int_{\mathcal{Y}} r^{(\phi,c)}(x, y) dP_Y(y) \right) dP_X(x) \quad (119) \\ &= \int_{\mathcal{X}} \hat{f}(x) dP_X(x); \quad \text{and similarly} \end{aligned} \quad (120) $$

$$ \int_{\mathcal{X} \times \mathcal{Y}} r^{(\phi,c)}(x, y) \hat{g}(y) dP_X P_Y(x, y) = \int_{\mathcal{Y}} \hat{g}(y) dP_Y(y). \quad (121) $$

Substitution of (120) and (121) into (118), combined with $G(A(r^{(\phi,c)})) = 0$, shows that the primal functional $P$ defined in (23), when evaluated at $r^{(\phi,c)}$,
satisfies

$$ P(r^{(\phi,c)}) = \int_X \hat{f}(x)dP_X(x) + \int_Y \hat{g}(y)dP_Y(y) \\ -\lambda \int_{X \times Y} \phi^* \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) dP_X P_Y(x, y). \quad (122) $$

Equivalently, $P(r^{(\phi,c)}) = D(\hat{f}, \hat{g}) = \mathcal{D}$, where $D$ is the dual functional in (44) and $\mathcal{D}$ is in (31). Since Lemma 2.5 established that $\mathcal{P} \ge \mathcal{D}$, with $\mathcal{P}$ in (22), it follows that $\mathcal{P} = \mathcal{D}$ and that the minimum in the primal problem is achieved at $r^{(\phi,c)}$. Strict convexity of the primal functional $P$ defined in (23) ensures uniqueness of the solution. ■

## 2.5.2 Sufficient conditions

The validity of Theorem 2.2 rests on the assumption that optimal potentials belong to $\mathcal{A}_{\phi,c}$. Intuitively, this condition requires that the optimal potentials stay uniformly bounded away from the singularity of the dual function $\phi^{*}$. While verifying this property directly for general Polish spaces can be challenging, it holds automatically in several significant scenarios commonly encountered in applications. This section explores sufficient conditions ensuring this requirement is met.

The first sufficient condition concerns divergences with superlinear growth, such as the Kullback-Leibler divergence. In this setting, the domain of the conjugate function $\phi^*$ extends to the entire real line, effectively removing the boundary singularity.

**Lemma 2.12.** If $\beta_{\phi}$ in (7) is equal to $+\infty$, any pair of optimal potentials is in $\mathcal{A}_{\phi,c}$.

*Proof:* Let $(\hat{f}, \hat{g}) \in L^{\infty}(\mathcal{X}, P_X) \times L^{\infty}(\mathcal{Y}, P_Y)$ be a pair of optimal potentials. Since $\hat{f}$ and $\hat{g}$ are essentially bounded and the cost function $c$ is bounded by assumption, the function $\frac{\hat{f} \oplus \hat{g} - c}{\lambda}$ is essentially bounded with respect to the product measure $P_X P_Y$. Therefore, there exists a bounded set $B \subset \mathbb{R}$ such that

$$ \mathrm{Ran} \left( \frac{\hat{f} \oplus \hat{g} - c}{\lambda} \right)_{\substack{\mathrm{a.s.} \\ P_X P_Y}} B. \quad (123) $$

Consequently, this range is contained in some interval $(-\infty, M)$ for some finite $M$. Given that $\beta_{\phi} = +\infty$, the condition $M < \beta_{\phi}$ is satisfied trivially, as the upper bound of the domain is infinite. ■

The second sufficient condition addresses the case of discrete optimal transport, where the underlying spaces are finite sets (e.g., empirical measures, histograms). Here, the finiteness of the space prevents the values of the potentials from accumulating asymptotically against the boundary $\beta_{\phi}$.

**Lemma 2.13.** If the sets $\mathcal{X}$ and $\mathcal{Y}$ are both finite, then any pair of optimal potentials is in the set $\mathcal{A}_{\phi,c}$ in (96).
**Proof:** Let $(\hat{f}, \hat{g})$ be a pair of optimal potentials for the problem in (31). The equalities in (94) and (95) are obtained from Lemma 2.10. From Lemma 2.7, there exists a set $\Gamma \subset \mathcal{X} \times \mathcal{Y}$ satisfying $P_X P_Y(\Gamma) = 1$ and, for all $(x, y) \in \Gamma$,

$$ \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} < \beta_{\phi}. \qquad (124) $$

Since $\mathcal{X}$ and $\mathcal{Y}$ are finite, $\Gamma$ is also finite. Therefore, the following quantity is well defined:

$$ M_{\hat{f}, \hat{g}} \triangleq \max_{(x,y) \in \Gamma} \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda}, \qquad (125) $$

and moreover $M_{\hat{f}, \hat{g}} < \beta_{\phi}$. Since $P_X P_Y(\Gamma) = 1$, it holds that $\frac{\hat{f} \oplus \hat{g} - c}{\lambda} \stackrel{\text{a.s.}}{\le} M_{\hat{f}, \hat{g}}$, with $M_{\hat{f}, \hat{g}} < \beta_{\phi}$. This confirms that the pair of optimal potentials $(\hat{f}, \hat{g})$ belongs to $\mathcal{A}_{\phi, c}$, defined in (96). ■

# 3 Equivalence of the problem to regularization on the family of $\phi$-divergences

## 3.1 Main results

The following theorem constitutes the main result of this report. The result and its proof have been adapted from [16, Theorem 2], which develops in the context of empirical risk minimization with $\phi$-divergence regularization. Although there are fundamental technical differences due to the presence of additional constraints in the case of optimal transport, the two results convey the same meaning in their corresponding scenarios.

**Theorem 3.1.** Assume that the problem in (31) admits a pair of optimal potentials $(\hat{f}, \hat{g}) \in \mathcal{A}_{\phi,c}$, and let $\psi$ be a Legendre-type generator (Definition 1.5). Then, there exists a bounded function $v: \mathcal{X} \times \mathcal{Y} \to \mathbb{R}$ such that the problem in (13) is equivalent, in the sense that both problems admit the same unique minimizer, to the following optimization problem:

$$ \inf_{P_{XY} \in \Pi(P_X, P_Y)} \int_{\mathcal{X} \times \mathcal{Y}} v(x, y) dP_{XY}(x, y) + \lambda D_{\psi}(P_{XY} \| P_X P_Y). \quad (126) $$

**Proof:** Let $(\hat{f}, \hat{g}) \in \mathcal{A}_{\phi,c}$ be maximizing potentials for the dual problem in (31). Define the cost function $v$ as

$$ v(x, y) = -\lambda \psi' \left( \phi'^{-1} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) \right) + 1. \quad (127) $$

Given that the functions $c$, $\hat{f}$ and $\hat{g}$ are bounded and that the range of $(x, y) \mapsto \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda}$ is bounded and contained in some interval of the form $(-\infty, \gamma)$,
with $\gamma < \beta_{\phi}$, the range of $(x, y) \mapsto \phi'^{-1}\left(\frac{\hat{f}(x)+\hat{g}(y)-c(x,y)}{\lambda}\right)$ is bounded and contained in an interval of the form $(\rho, +\infty)$ with $\rho > 0$. Since $\psi'$ is $C^1$ on $(0, \infty)$ (Lemma 1.1), the function

$$ (x, y) \mapsto \psi' \left( \phi'^{-1} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) \right) \quad (128) $$

is also bounded. Thus, the function $v$ is bounded, ensuring that problem (126) fits within the setting established in previous sections. Let $\tilde{f} \in L^{\infty}(\mathcal{X}, P_X)$ and $\tilde{g} \in L^{\infty}(\mathcal{Y}, P_Y)$ be defined by

$$ \tilde{f}(x) = \frac{1}{2}, \quad \forall x \in \mathcal{X}; \quad (129) $$

$$ \tilde{g}(y) = \frac{1}{2}, \quad \forall y \in \mathcal{Y}. \quad (130) $$

The proof proceeds by demonstrating that the functions $\tilde{f}, \tilde{g}$ satisfy (97) and (98) for problem (126). Verification proceeds as follows: let $x \in \mathcal{X}$ be fixed, and compute the integral

$$ \int_{\mathcal{Y}} \psi'^{-1} \left( \frac{\tilde{f}(x) + \tilde{g}(y) - v(x, y)}{\lambda} \right) dP_Y(y) \\ = \int_{\mathcal{Y}} \psi'^{-1} \left( \frac{1 + \lambda \psi' \left( \phi'^{-1} \left( \frac{\hat{f}(x)+\hat{g}(y)-c(x,y)}{\lambda} \right) \right) - 1}{\lambda} \right) dP_Y(y) \quad (131) $$

$$ = \int_{\mathcal{Y}} \phi'^{-1} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) dP_Y(y). \quad (132) $$

Since $(\hat{f}, \hat{g}) \in \mathcal{A}_{\phi,c}$ is the solution to (31), Lemma 2.11 implies

$$ \int_{\mathcal{Y}} \phi'^{-1} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) dP_Y(y) \stackrel{a.s.}{=} 1. \quad (133) $$

Combining (132) and (133), yields

$$ \int_{\mathcal{Y}} \psi'^{-1} \left( \frac{\tilde{f}(x) + \tilde{g}(y) - v(x, y)}{\lambda} \right) dP_Y(y) \stackrel{a.s.}{=} 1. \quad (134) $$

Moreover, applying the same reasoning for a fixed $y \in \mathcal{Y}$ yields

$$ \int_{\mathcal{X}} \psi'^{-1} \left( \frac{\tilde{f}(x) + \tilde{g}(y) - v(x, y)}{\lambda} \right) dP_X(x) \stackrel{a.s.}{=} 1. \quad (135) $$

Lemma 2.11 implies that $(\hat{f}, \hat{g})$ are optimal potentials for the dual of (126), which in turn requires verifying that $(\tilde{f}, \tilde{g}) \in \mathcal{A}_{\psi,v}$. This set is defined analogously to $\mathcal{A}_{\phi,c}$ in (96), by replacing the parameters $\phi$ and $c$ with $\psi$ and $v$ respectively. This condition requires that the range of $l: (x, y) \mapsto \frac{\hat{f}(x)+\hat{g}(y)-v(x,y)}{\lambda}$
satisfies $\mathrm{Ran}(l)_{\hat{P}_X \hat{P}_Y}^{\mathrm{a.s.}}(-\infty, \gamma)$ for some $\gamma < \beta_\psi$. Let $(x, y) \in \mathcal{X} \times \mathcal{Y}$ be fixed, and consider

$$ \frac{\tilde{f}(x) + \tilde{g}(y) - v(x, y)}{\lambda} = \psi' \left( \phi'^{-1} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) \right). \quad (136) $$

Since the pair $(\hat{f}, \hat{g})$ is in the set $\mathcal{A}_{\phi,c}$ defined in (96), and due to continuity and strict monotonicity of $\phi'^{-1}$ on $(-\infty, \beta_\phi)$ (see Section 1.2), the range of $(x, y) \mapsto \phi'^{-1} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right)$ is bounded and satisfies

$$ \mathrm{Ran} \left( \phi'^{-1} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) \right)_{\hat{P}_X \hat{P}_Y}^{\mathrm{a.s.}} (\rho, +\infty), \quad (137) $$

for some $\rho > 0$. Given that $\psi'$ is continuous and strictly increasing on $(0, +\infty)$, the range of $(x, y) \mapsto \psi' \left( \phi'^{-1} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) \right)$ satisfies

$$ \mathrm{Ran} \left( \psi' \left( \phi'^{-1} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) \right) \right)_{\hat{P}_X \hat{P}_Y}^{\mathrm{a.s.}} (-\infty, \gamma), \quad (138) $$

for some $\gamma < \beta_\psi$. Therefore, by (136), $(\tilde{f}, \tilde{g}) \in \mathcal{A}_{\psi,v}$ and satisfies (134) and (135). Lemma 2.11 implies that $(\tilde{f}, \tilde{g})$ are maximizing potentials for the dual problem of (126).

Finally, Theorem 2.2 gives the form of the solution of the problem in (126), denoted by $P_{XY}^{(\nu, \psi)}$ as

$$ \frac{\mathrm{d}P_{XY}^{(\nu, \psi)}}{\mathrm{d}P_X P_Y}(x, y) = \psi'^{-1} \left( \frac{\tilde{f}(x) + \tilde{g}(y) - v(x, y)}{\lambda} \right) \quad (139) $$

$$ = \psi'^{-1} \left( \psi' \left( \phi'^{-1} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) \right) \right) \quad (140) $$

$$ = \phi'^{-1} \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) \quad (141) $$

$$ = \frac{\mathrm{d}P_{XY}^{(c, \phi)}}{\mathrm{d}P_X P_Y}(x, y), \quad (142) $$

where $P_{XY}^{(c, \phi)}$ denotes the solution to the problem in (22). The two problems thus share the same minimizer. ■

## 3.2 Examples

This section illustrates the structural equivalence established in Theorem 3.1 by considering fixed Polish spaces $\mathcal{X}$ and $\mathcal{Y}$, a fixed bounded cost function $c$, and a fixed regularization parameter $\lambda > 0$. Let $\phi$ denote the generator of the
Kullback-Leibler (KL) divergence, defined by $\phi(x) = x \log(x) - x + 1$. For this KL generator, the inverse of its derivative is given by

$$ \phi'^{-1}(y) = \exp(y). \qquad (143) $$

Moreover, the quantity $\beta_{\phi}$ in (7) is equal to $+\infty$. Therefore, it follows from Lemma 2.12 that the corresponding optimal potentials for the dual problem in (31), denoted by $\hat{f}$ and $\hat{g}$, belong to the set $A_{\phi,c}$ in (96). As a consequence, all the assumptions of Theorem 3.1 are satisfied for any target Legendre-type generator $\psi$ satisfying Definition 1.5. In particular, for the transformed cost function defined in (127), the problems in (13) and (126) admit the same unique minimizer. The following examples illustrate the resulting transformed cost function $v$ obtained via (127) for various choices of target regularization.

### 3.2.1 From KL to reverse KL regularization

Consider the target regularization defined by the reverse KL generator $\psi(x) = -\log(x) + x - 1$, whose derivative is $\psi'(x) = 1 - \frac{1}{x}$. From (127) and (143), the transformed cost function $v$ ensuring that the initial problem (13) and the transformed problem (126) share the same unique minimizer is given by

$$ v(x, y) = -\lambda \psi' \left( \exp \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) \right) + 1 \qquad (144) $$

$$ = -\lambda \left( 1 - \exp \left( - \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) \right) + 1 \qquad (145) $$

$$ = \lambda \exp \left( \frac{c(x, y) - \hat{f}(x) - \hat{g}(y)}{\lambda} \right) - \lambda + 1. \qquad (146) $$

### 3.2.2 From KL to Jensen-Shannon regularization

In the case of a target regularization defined by the Jensen-Shannon divergence, the corresponding generator is $\psi(x) = x \log(x) - (x+1) \log\left(\frac{x+1}{2}\right)$, with derivative $\psi'(x) = \log\left(\frac{2x}{x+1}\right)$. The transformed cost function $v$ preserving the minimizer is

$$ v(x, y) = -\lambda \log \left( \frac{2 \exp \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right)}{\exp \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) + 1} \right) + 1 \qquad (147) $$

$$ = \lambda \log \left( \frac{\exp \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right) + 1}{2 \exp \left( \frac{\hat{f}(x) + \hat{g}(y) - c(x, y)}{\lambda} \right)} \right) + 1. \qquad (148) $$

By simplifying the argument within the logarithm, the relation becomes:

$$ v(x, y) = \lambda \log \left( \frac{1 + \exp \left( \frac{c(x, y) - \hat{f}(x) - \hat{g}(y)}{\lambda} \right)}{2} \right) + 1. \qquad (149) $$
These examples confirm that changing the regularizing divergence manifests as
a non-linear adjustment of the cost function, ensuring that the unique minimizer
of the regularized optimal transport problems remain unchanged.

# 4 Conclusions

This work studies $\phi$-divergence-regularized optimal transport problems of the form in (13) under the assumptions that $\mathcal{X}$ and $\mathcal{Y}$ are polish spaces, the cost function $c$ is bounded, and $\phi$ is a Legendre-type generator (Definition 1.5). Section 2.4 establishes the existence of optimal dual potentials for the problem in (31) through the properties of the $(c, \lambda, \phi)$-transform. Section 2.5 establishes the existence of an optimal coupling for the problem in (22), including a closed-form expression, under conditions that Section 2.5.2 shows to follow from sufficient assumptions such as the superlinearity of the generator or the finiteness of the Polish spaces.

The central contribution, stated in Theorem 3.1, highlights a structural equivalence between distinct regularization schemes: replacing the regularizing divergence with another $\psi$-divergence admits compensation through an appropriate modification of the cost function, while preserving the minimizer of the problem. This result offers a unifying theoretical perspective on a broad class of regularized optimal transport formulations and demonstrates that the resulting solutions depend primarily on the interaction between the cost and the regularization term, rather than on the specific choice of divergence alone. The result remains essentially theoretical, since the transformed cost function relies on the solution of the dual problem associated with the original formulation in (127). As a result, the established equivalence does not yield computational improvements in practice. Nevertheless, this perspective may retain practical relevance for sampling purposes, as it provides an alternative characterization of the optimal coupling.
References

[1] I. Tolstikhin, O. Bousquet, S. Gelly, and B. Schölkopf, “Wasserstein auto-encoders,” in *Proceedings of the International Conference on Learning Representations (ICLR)*, May 2018.

[2] M. Arjovsky, S. Chintala, and L. Bottou, “Wasserstein GAN,” in *Proceedings of the International Conference on Machine Learning (ICML)*, Jul. 2017, pp. 214–223.

[3] N. Courty, R. Flamary, D. Tuia, and A. Rakotomamonjy, “Optimal transport for domain adaptation,” *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 39, no. 9, pp. 1852–1865, Sep. 2017.

[4] C. Frogner, C. Zhang, H. Mobahi, M. Araya-Polo, and T. Poggio, “Learning with a wasserstein loss,” in *Proceedings of the 29th International Conference on Neural Information Processing Systems - Volume 2*, Dec. 2015, pp. 2053–2061.

[5] J. Solomon, F. de Goes, G. Peyré, M. Cuturi, A. Butscher, A. Nguyen, T. Du, and L. Guibas, “Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains,” *ACM Transactions on Graphics (TOG)*, vol. 34, no. 4, pp. 1–11, Jul. 2015.

[6] L. M. Pereira and M. H. Amini, “A survey on optimal transport for machine learning: Theory and applications,” *IEEE Access*, vol. 13, pp. 26506–26526, Jan. 2025.

[7] C. Villani, *Optimal Transport: Old and New*, 1st ed. Berlin, Heidelberg: Springer-Verlag, 2009.

[8] G. Peyré and M. Cuturi, *Computational Optimal Transport*, 1st ed. Hanover, MA, USA: Foundations and Trends in Machine Learning, 2020.

[9] M. Cuturi, “Sinkhorn distances: Lightspeed computation of optimal transportation distances,” in *Proceedings of the International Conference on Neural Information Processing Systems (NeurIPS)*, vol. 2, Dec. 2013, pp. 2292–2300.

[10] M. Blondel, V. Seguy, and A. Rolet, “Smooth and sparse optimal transport,” in *Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS)*, Mar. 2018, pp. 314–323.

[11] S. Di Marino and A. Gerolin, “Optimal transport losses and sinkhorn algorithm with general convex regularization,” *arXiv preprint arXiv:2007.00976*, Jul. 2020.

[12] D. Terjék and D. González-Sánchez, “Optimal transport with f-divergence regularization and generalized sinkhorn algorithm,” in *Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS)*, vol. 151, Mar. 2022, pp. 5135–5165.
[13] O. Catoni, *PAC-Bayesian Supervised Classification: The Thermodynamics of Statistical Learning*, 1st ed. Beachwood, OH, USA: Institute of Mathematical Statistics, 2007, vol. 56.

[14] S. M. Perlaza, G. Bisson, I. Esnaola, A. Jean-Marie, and S. Rini, “Empirical risk minimization with relative entropy regularization,” *IEEE Transactions on Information Theory*, vol. 70, no. 7, pp. 5122–5161, Jul. 2024.

[15] F. Daunas, I. Esnaola, S. M. Perlaza, and H. V. Poor, “Empirical risk minimization with relative entropy regularization Type-II,” INRIA, Centre Inria d'Université Côte d'Azur, Sophia Antipolis, France, Tech. Rep. RR-9508, May 2023.

[16] F. Daunas, I. Esnaola, S. M. Perlaza, and H. V. Poor, “Equivalence of empirical risk minimization to regularization on the family of $f$-divergences,” in *Proceedings of the IEEE International Symposium on Information Theory (ISIT)*, Jul. 2024, pp. 759–764.

[17] F. Daunas, I. Esnaola, S. M. Perlaza, and H. V. Poor, “Empirical risk minimization with $f$-divergence regularization in statistical learning,” INRIA, Centre Inria d'Université Côte d'Azur, Sophia Antipolis, France, Tech. Rep. RR-9521, Oct. 2023.

[18] C. Zalinescu, *Convex analysis in general vector spaces*, 1st ed. Singapore: World Scientific, 2002.

[19] J. M. Borwein and A. S. Lewis, “Partially-finite programming in $l_1$ and the existence of maximum entropy estimates,” *SIAM Journal on Optimization*, vol. 3, no. 2, pp. 248–267, May 1993.

[20] E. M. Werner, “$f$-divergence for convex bodies,” in *Asymptotic Geometric Analysis: Proceedings of the Fall 2010 Fields Institute Thematic Program*. New York, NY: Springer, 2013, pp. 381–395.

[21] I. Csiszár, “Information-type measures of difference of probability distributions and indirect observations,” *Studia Scientiarum Mathematicarum Hungarica*, vol. 2, pp. 299–318, 1967.

[22] J. M. Borwein and A. S. Lewis, “Partially-finite programming in $l_1$ and the existence of maximum entropy estimates,” *SIAM Journal on Optimization*, vol. 3, no. 2, pp. 248–267, May 1993.

[23] R. Agrawal and T. Horel, “Optimal bounds between $f$-divergences and integral probability metrics,” *Journal of Machine Learning Research*, vol. 22, no. 1, pp. 5662–5720, Jan. 2021.

[24] Y. Bermudez, G. Bisson, I. Esnaola, and S. M. Perlaza, “Proofs for folklore theorems on the Radon-Nikodym derivative,” INRIA, Centre Inria d'Université Côte d'Azur, Sophia Antipolis, France, Tech. Rep. RR-9591, Jul. 2025.
[25] W. Rudin, *Real and Complex Analysis*, 3rd ed. New York, NY, USA: McGraw-Hill, 1987.

[26] R. E. Megginson, *An Introduction to Banach Space Theory*, 1st ed. New York, NY, USA: Springer, 1998.

[27] H. H. Bauschke and P. L. Combettes, *Convex Analysis and Monotone Operator Theory in Hilbert Spaces*, 1st ed. New York, NY, USA: Springer, 2011.

[28] I. Ekeland and R. Temam, *Convex Analysis and Variational Problems*, 1st ed. Philadelphia, PA, USA: SIAM, 1999.

[29] Y. Wu, “Lecture notes on information-theoretic methods in statistics,” New Haven, CT, Jan. 2020, Department of Statistics and Data Science, Yale University.

[30] R. T. Rockafellar and R. J.-B. Wets, *Variational Analysis*, 1st ed. Berlin, Germany: Springer, 1998.

[31] H. Brézis, *Analyse Fonctionnelle: Théorie et Applications*, 1st ed. Paris, France: Masson, 1983.

[32] W. Rudin, *Principles of Mathematical Analysis*, 3rd ed. New York, NY, USA: McGraw-Hill, 1976.

[33] D. G. Luenberger, *Optimization by Vector Space Methods*, 1st ed. New York, NY, USA: John Wiley & Sons, 1969.
informatics mathematics
Inria

RESEARCH CENTRE
SOPHIA ANTIPOLIS – MÉDITERRANÉE

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-6399