# On some generalisations of Optimal Transport problem

Luca Nenna

► To cite this version:

Luca Nenna. On some generalisations of Optimal Transport problem. Optimization and Control [math.OC]. Paris Saclay University, 2024. tel-04554531

HAL Id: tel-04554531

https://inria.hal.science/tel-04554531

Submitted on 22 Apr 2024

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 Attribution 4.0 International License
# Sur quelques généralisations du problème du transport optimal

Habilitation à diriger des recherches de l'Université Paris-Saclay

présentée et soutenue à Orsay, le 06/03/2024, par

Luca NENNA

## Composition du jury

**Yann BRENIER**
DR, CNRS et Université Paris-Saclay

Rapporteur

**Young-Heon KIM**
Professeur des universités, UBC

Rapporteur

**Gabriel PEYRÉ**
DR, CNRS et ENS Paris

Rapporteur

**Julie DELON**
Professeure des universités, Université Paris-Cité

Examinatrice

**Jean-Marie MIREBEAU**
DR, ENS Paris-Saclay

Examineur

**Simona ROTA-NODARI**
Professeure des universités, Université Côte d'Azur

Examinatrice
per aspera ad astra
On pains au chocolat, electrons, entropy and
the theory (and numerics) to rule them all.
Remerciements

Le temps passe vite, même trop vite, et l'écriture de ce manuscrit pour l'Habilitation à Diriger les Recherches a été l'occasion pour moi de mesurer une fois de plus l'immense chance que j'ai de faire un travail qui me passionne. Je souhaiterais profiter de cette occasion pour adresser mes plus profonds remerciements aux différentes personnes qui m'ont permis d'avoir cette chance.

Tout d'abord, je voudrais remercier très chaleureusement Jean-David Benamou et Guillaume Carlier pour avoir accepté, il y a longtemps désormais d'encadrer ma thèse de doctorat. Je pense que si j'avais fait un choix différent il y a presque 11 ans, je ne serais pas si heureux (et fier) du travail accompli et plein d'énergie pour tout le travail qui m'attend. Encore merci pour tout ce que vous m'avez donné.

I am deeply thankful to Gabriel Peyré and Young-Heon Kim for accepting to be referees of this thesis. I am glad that they appreciated this work and I have found their comments very inspiring.

Je tiens également à remercier Yann Brenier (merci pour les discussions sur les mathématiques, la politique et les trains...ça a été et ça sera toujours un vrai plaisir!!), Julie Delon, Jean-Marie Mirebeau et Simona Rota-Nodari qui m'ont fait l'honneur d'être les membres de mon jury.

Un moka-merci à l'équipe MOKAPLAN (Vincent, JD, Antonin, Irène, Paul, Flavien, Derya et tous les moka-thésards et aussi les moka-anciens) qui m'a accueilli pendant ma thèse et en-core aujourd'hui me donne abri au moins une fois par semaine dans ses locaux sur Paris! Vous devez être fiers de cette communauté que vous avez construit et que vous continuez à construire. Toujours côté Inria, je voudrais remercier l'équipe MATHERIALS pour m'avoir accueilli en délégation ces deux dernières années (délégation qui m'a permis de rédiger ce manuscrit). Bien sûr je peux pas oublier le 4ème étage du bâtiment A

J'ai rencontré beaucoup trop de personnes passionnantes ces dernières années (une pensée spéciale pour Rodrigue Lelotte et Mathias Dus avec son alter ego J.-C. ); parmi elles merci encore à mes collaborateur.trice.s : Brendan Pass, Mathieu Lewin, Simone Di Marino, Paul Pagon, Virginie Ehrlacher, Luigi De Pascale, Quentin Mérigot, Léonard Monsaigeon, Jean-Baptiste Casteras, Hamza Ennaji, Augusto Gerolin...merci à vous tou.te.s car à chaque discussion, à chaque pause café ou match de tennis vous m'apprenez avec patience (et il faut beaucoup de patience avec moi) énormément de choses sur la vie et sur les maths. Je sais qu'il n'est pas toujours facile de travailler avec moi, mais j'espère que vous avez apprécié nos résultats et le temps que nous avons passé (et que nous passerons encore) ensemble.
Merci à Maxime, Paul, Quentin, Thomas et Léonard pour leur relectures, pas seulement de ce document, (oui Thomas j'ai copié tes remerciements ;) mais c'est factuel), pour GOTA, pour
ParMA, les matchs de tennis et les discussions de maths autour d'une bière fraîche sous le ciel bleu de Lisbonne.

Merci au Laboratoire de Mathématiques d'Orsay, et surtout à l'équipe ANEDP, pour m'avoir fait confiance et recruté il y a quelques années comme McF. J'espère que ce manuscrit puisse témoigner que la confiance n'était pas mal placée. Merci à Jean-Marie, Maya et Cyril pour avoir partagé et partager encore le bureau 2F1 avec la vue (plus ou moins) sur la bucolique vallée de l'Yvette. Et surtout merci à Thomas Galloüet, Quentin Mérigot, Yann Brenier, Bertrand Maury, Bruno Lévy, Hugo Leclerc et Sylvain Faure avec qui on commence cette nouvelle aventure avec l'équipe ParMA.

Merci à vous tou.te.s avec qui je partage ma vie à Paris (Ci, Simon, Max, Ricky, JLO, Dario, Franchie, Vituzzo, Meggie, André, Pierre, Damian, Paul, Samra...et les plus petits Lupo et VL...et encore les anciens parisiens Vale, Paolo, Anna et Davide...et j'ai peut-être oublié queluq'un), mes vacances, mes séjours à Lisbonne ou en Italie, à mes amis de toujours (Laura, Claudio, Monica, Daniele, Alby, Paolo, Carlotta, Eleonora×2), à mes parents, à ma famille et à Antonietta et Roger qui m'ont accueilli il y a presque 11 ans et m'ont permis de faire de Paris ma maison. À vous tou.te.s, je veux dédier ce manuscrit même si je sais que c'est rien par rapport au support, à l'amitié et à l'amour que vous me donnez chaque jour.
# Contents

<table><tr><td>Introduction (en français)</td><td>2</td></tr><tr><td>Introduction (in English)</td><td>10</td></tr><tr><td>1 Multi-Marginal Optimal Transport</td><td>21</td></tr><tr><td>1.1 Introduction to Multi-Marginal Optimal Transport</td><td>21</td></tr><tr><td>1.2 From Risk measures to multi-marginal optimal transport</td><td>24</td></tr><tr><td>1.2.1 Equivalence between maximizing spectral risk measures and multi-marginal optimal transport</td><td>26</td></tr><tr><td>1.2.2 Solutions for one-dimensional assets and compatible payouts</td><td>28</td></tr><tr><td>1.2.3 Higher dimensional assets</td><td>29</td></tr><tr><td>1.2.4 Multidimensional measures of risk</td><td>30</td></tr><tr><td>1.3 Entropic Multi-Marginal Optimal Transport</td><td>31</td></tr><tr><td>1.4 Rate of convergences for entropic multi-marginal optimal transport</td><td>33</td></tr><tr><td>1.4.1 Upper bound</td><td>34</td></tr><tr><td>1.4.2 Lower bound</td><td>35</td></tr><tr><td>1.4.3 Matching bound</td><td>36</td></tr><tr><td>1.5 An ODE characterisation of discrete entropic MMOT</td><td>37</td></tr><tr><td>1.5.1 Discrete Entropic Multi-Marginal Optimal Transport problem</td><td>37</td></tr><tr><td>1.5.2 An ODE characterisation of multi-marginal optimal transport</td><td>38</td></tr><tr><td>1.5.3 The algorithm</td><td>40</td></tr><tr><td>1.6 Perspectives</td><td>43</td></tr><tr><td>2 Generalized Optimal Transport in Mathematical Physics</td><td>46</td></tr><tr><td>2.1 Optimal Transport in Density Functional Theory and beyond</td><td>46</td></tr><tr><td>2.2 Grand Canonical Optimal Transport</td><td>49</td></tr><tr><td>2.2.1 The theory of Grand Canonical Optimal Transport</td><td>51</td></tr><tr><td>2.2.2 Entropic Grand Canonical Optimal Transport</td><td>57</td></tr><tr><td>2.3 Quantum Optimal Transport and Lieb Functional</td><td>60</td></tr><tr><td>2.3.1 Moment approximation of Lieb functional: the primal problem</td><td>61</td></tr><tr><td>2.3.2 Convergence to the exact Lieb functional</td><td>62</td></tr><tr><td>2.3.3 Convergence rate of the ground state energy in the bounded domain case</td><td>63</td></tr><tr><td>2.3.4 The dual problem</td><td>65</td></tr><tr><td>2.4 Perspectives</td><td>66</td></tr></table>
CONTENTS

<table><tr><td>3 Unequal dimensional Optimal Transport</td><td>68</td></tr><tr><td>3.1 Some variational problems involving Unequal dimensional Optimal Transport</td><td>70</td></tr><tr><td>3.1.1 A sufficient condition for nestedness</td><td>72</td></tr><tr><td>3.1.2 The congestion case</td><td>73</td></tr><tr><td>3.1.3 The interaction case and fixed point characterization</td><td>76</td></tr><tr><td>3.1.4 The double minimization</td><td>78</td></tr><tr><td>3.1.5 The hedonic pricing problem</td><td>79</td></tr><tr><td>3.2 A class of metrics involving singular measures</td><td>80</td></tr><tr><td>3.2.1 The v-based Wasserstein metric</td><td>82</td></tr><tr><td>3.2.2 Structure of geodesics for the v-based Wasserstein metric</td><td>84</td></tr><tr><td>3.2.3 Geodesic convexity of functionals</td><td>85</td></tr><tr><td>3.2.4 The μ-based Wasserstein metric</td><td>86</td></tr><tr><td>3.2.5 Fixed point characterization of solutions to variational problems</td><td>87</td></tr><tr><td>3.3 Perspectives</td><td>90</td></tr></table>
Notations

In all this manuscript, $d \in \mathbb{N}^*$ denotes the dimension of the ambient space $\mathbb{R}^d$ and $m \in \mathbb{N}$ is an integer such that $m \ge 2$.

<table>
  <tr>
    <td>X</td>
    <td>a subset of R<sup>d</sup></td>
  </tr>
  <tr>
    <td>X<sub>i</sub></td>
    <td>a subset of R<sup>d</sup> for any index i &isin; {1, ..., m}</td>
  </tr>
  <tr>
    <td>X</td>
    <td>product X<sub>1</sub> &times; ... &times; X<sub>m</sub> whenever (X<sub>i</sub>)<sub>1&le;i&le;m</sub> is a family of m subsets of R<sup>d</sup>;</td>
  </tr>
  <tr>
    <td>X&#773;<sub>-i</sub></td>
    <td>product &Pi;<sub>1&le;j&le;m, j&ne;i</sub> X<sub>j</sub></td>
  </tr>
  <tr>
    <td>x<sub>i</sub>, x, x</td>
    <td>a point in X&#773;<sub>i</sub>, in some X, and in X respectively;</td>
  </tr>
  <tr>
    <td>x<sub>q</sub></td>
    <td>(x<sub>i</sub>)<sub>i&isin;q</sub> if q &sube; {1, ..., m} and x &isin; X;</td>
  </tr>
  <tr>
    <td>e<sub>i</sub></td>
    <td>i-th coordinate map e<sub>i</sub>: x = (x<sub>1</sub>, ..., x<sub>m</sub>) &mapsto; x<sub>i</sub>;</td>
  </tr>
  <tr>
    <td>A &isin; B</td>
    <td>A &sube; K &sube; B for some compact set K;</td>
  </tr>
  <tr>
    <td>&lang;&middot;, &middot;&lang;</td>
    <td>Euclidean scalar product;</td>
  </tr>
  <tr>
    <td>|&middot;|</td>
    <td>Euclidean norm on R<sup>d</sup>;</td>
  </tr>
  <tr>
    <td>||&middot;||</td>
    <td>norm on R<sup>d</sup> &times; ... &times; R<sup>d</sup> defined by ||x|| = max<sub>1&le;i&le;m</sub> |x<sub>i</sub>| if x = (x<sub>1</sub>, ..., x<sub>m</sub>);</td>
  </tr>
  <tr>
    <td>B<sub>r</sub>(x)</td>
    <td>open ball of radius r centered at x &isin; (R<sup>d</sup>)<sup>m</sup> for the above norm;</td>
  </tr>
  <tr>
    <td>C<sup>0,1</sup><sub>loc</sub>(X)</td>
    <td>space of real-valued locally Lipschitz functions on X which is a sub-manifold of R<sup>d</sup> or (R<sup>d</sup>)<sup>m</sup>;</td>
  </tr>
  <tr>
    <td>[f]C<sup>0,1</sup><sub>loc</sub>(X)</td>
    <td>Lipschitz constant of f: X &rarr; R where X is a subset of R<sup>d</sup> or (R<sup>d</sup>)<sup>m</sup> for the above norms;</td>
  </tr>
  <tr>
    <td>C<sup>1,1</sup><sub>loc</sub>(X)</td>
    <td>space of differentiable real-valued functions on X, a sub-manifold of R<sup>d</sup> or (R<sup>d</sup>)<sup>m</sup>, with locally Lipschitz differential;</td>
  </tr>
  <tr>
    <td>P(X)</td>
    <td>space of probability measures on a metric space X;</td>
  </tr>
  <tr>
    <td>&Pi;(&mu;<sub>1</sub>, ..., &mu;<sub>m</sub>)</td>
    <td>space of probability measures on X having &mu;<sub>1</sub>, ..., &mu;<sub>m</sub> as marginals;</td>
  </tr>
  <tr>
    <td>&Pi;(&rho;)</td>
    <td>space of probability measures on X having all marginals equals to &rho;;</td>
  </tr>
  <tr>
    <td>||&middot;||<sub>L<sup>p</sup>(&mu;)</sub></td>
    <td>L<sup>p</sup> norm induced by a measure &mu;, where p &isin; [1, +&infin;];</td>
  </tr>
  <tr>
    <td>spt &mu;</td>
    <td>support of the measure &mu;;</td>
  </tr>
  <tr>
    <td>H<sup>s</sup><sub>X</sub></td>
    <td>s-dimensional Hausdorff measure on the metric space X endowed with the Borel &sigma;-algebra (the subscript X will often be dropped);</td>
  </tr>
  <tr>
    <td>M<sub>d</sub>(R)</td>
    <td>space of real matrices of size d &times; d, endowed with the Frobenius norm induced by the scalar product A &middot; B := Tr(A<sup>T</sup>B), for A, B &isin; M<sub>d</sub>(R);</td>
  </tr>
  <tr>
    <td>S<sub>d</sub>(R)</td>
    <td>subspace of real symmetric matrices of size d &times; d;</td>
  </tr>
  <tr>
    <td>&Delta;<sub>P</sub></td>
    <td>simplex of P-uples t = (t<sub>p</sub>)<sub>p&isin;P</sub> such that t<sub>p</sub> &ge; 0 for all p &isin; P and &Sigma;<sub>p&isin;P</sub> t<sub>p</sub> = 1;</td>
  </tr>
  <tr>
    <td>H<sup>m</sup><sub>1</sub></td>
    <td>is the set of admissible electronic wavefunctions for a system of m electrons with finite kinetic energy, that is H<sup>m</sup><sub>1</sub> := &Lambda;<sup>m</sup><sub>i=1</sub> H<sup>1</sup>(R<sup>3</sup>). In the same fashion we define H<sup>m</sup><sub>0</sub> := &Lambda;<sup>m</sup><sub>i=1</sub> L<sup>2</sup>(R<sup>3</sup>) and H<sup>m</sup><sub>0</sub> := &Lambda;<sup>m</sup><sub>i=1</sub> H<sup>2</sup>(R<sup>3</sup>);</td>
  </tr>
  <tr>
    <td>G<sup>+</sup><sub>1</sub>(H<sup>m</sup><sub>0</sub>)</td>
    <td>denotes the set of trace-class self-adjoint non-negative operators on H<sup>m</sup><sub>0</sub>.</td>
  </tr>
</table>
Introduction (en français)

Le transport optimal (OT) est un domaine de recherche très dynamique et il y a tant d'exemples et d'applications qu'il est difficile de choisir l'un d'entre eux pour introduire ce sujet fascinant; nous mentionnons les manuels "old and new" comme [AGS04; Vil09; Vil03; San15; Gal16; PC19]. Nous avons décidé d'introduire le problème du transport optimal classique en utilisant un exemple donné par Villani dans [Vil03].

Situons notre contexte à Paris et prenons en considération l'ensemble des boulangeries qui produisent des pains au chocolat/chocolatines (le lecteur peut bien sûr choisir la viennoiserie qui lui plaît le plus). Ces derniers doivent être livrés chaque matin aux cafés dans lesquels les clients pourront les goûter. La production et la consommation de pains au chocolat sont décrites par $\mu_1$ et $\mu_2$, respectivement. Nous supposons que la quantité total de la production et de la consommation est la même. Le problème du transport optimal consiste alors à transporter la quantité de pains au chocolat/chocolatines produite par la boulangerie $x_1 \in X_1$ ($X_1 \subset \mathbb{R}^{d_1}$ est l'ensemble des boulangeries caractérisées par $d_1$ propriétés e.g., emplacement précis, taille, etc.) à un café $x_2 \in X_2$ ($X_2 \subset \mathbb{R}^{d_2}$ est l'ensemble des cafés caractérisés par $d_2$ propriétés) de sorte que le coût de transport $c(x_1, x_2)$, par exemple la distance entre la boulangerie et le café, soit minimal. En transport optimal classique, nous avons $d_1 = d_2$, mais dans ce qui suit, nous allons essayer d'aborder le problème lorsque le nombre de caractéristiques (alias les dimensions) diffère. Le problème introduit par Monge dans [Mon81] consiste à chercher un transport $x_2 = T(x_1)$ qui nous indique le café $x_2$ où l'intégralité des pains au chocolat produits par $x_1$ sera délivrés. Toutefois, ce problème est difficile à traiter et on nécessite une relaxation (connue sous le nom de problème de Monge-Kantorovich) : on permet, par exemple, qu'une partie des pains au chocolat produite en $x_1$ soit envoyée au café $x_2$ et l'autre à $x'_2$. Dans ce cas, nous cherchons un couplage optimal $\gamma(x_1, x_2)$ qui nous indique comment la masse en $x_1$ est répartie sur chaque $x_2 \in X_2$. Le problème peut alors être formulé comme suit : étant donné deux mesures de probabilité $\mu_1 \in \mathcal{P}(X_1)$ et $\mu_2 \in \mathcal{P}(X_2)$ et une fonction de coût continue $c: X_1 \times X_2 \to \mathbb{R}$ le problème de Monge-Kantorovich est définie comme suit

$$ \mathrm{OT}^c(\mu_1, \mu_2) := \inf_{\gamma \in \Pi(\mu_1, \mu_2)} \int_{X_1 \times X_2} c(x_1, x_2) d\gamma(x_1, x_2), \quad (1) $$

où $\Pi(\mu_1, \mu_2)$ désigne l'ensemble des mesures de probabilité sur $X_1 \times X_2$ ayant $\mu_1$ et $\mu_2$ comme marginales, c'est-à-dire qu'étant donnée la projection canonique $e_i: X_1 \times X_2 \to X_i$, nous avons $e_i, \sharp \gamma = \mu_i$ pour $i=1,2$.

Si le couplage optimal $\gamma(x_1, x_2)$ consiste à assigner un unique café à chaque boulangerie $x_1$, alors on dit que le couplage est déterministe et qu'un transport optimal $T$ existe. Lorsque le coût $c(x_1, x_2)$ est la distance au carré, l'existence d'un couplage déterministe a été prouvée par Yann Brenier dans [Bre91]; de plus, dans ce cas, (1) est une distance au carré sur l'ensemble
des mesures de probabilité.

De nombreuses questions se posent à ce stade : que se passe-t-il si le nombre de caractéristiques des boulangeries et des cafés est différent (e.g., $d_1 \neq d_2$)? Et si l'une des distributions est inconnue et que nous souhaitons l'obtenir en minimisant le OT$^c$ plus une autre fonctionnelle ? De plus, nous pourrions ajouter une troisième catégorie, les hôtels, de sorte que nous devons maintenant minimiser le coût d'envoi des pains au chocolat/chocolatines à la fois aux cafés et aux hôtels à Paris. Ce problème peut-il encore être interprété comme un problème d'OT ? Ce manuscrit (qui résume en partie mes recherches après la thèse) est, en fait, consacré à répondre à ces questions et à quelques autres concernant les modèles de transport optimal apparaissant en physique, en particulier en physique quantique.

Plus en détail, ce manuscrit est divisé en trois chapitres organisés, selon la thématique, comme suit :

* Le chapitre 1 est consacré (1) au transport optimal multi-marginales et à une appli-
cation aux mesures de risque, (2) aux taux de convergence pour le transport optimal
multi-marginales entropique et, enfin, (3) à une caractérisation du transport optimal multi-
marginales entropique (discret) par une EDO. Mes co-auteurs sur ces sujets sont H. Ennaji,
Q. Mérigot, B. Pass et P. Pegon et les résultats des articles suivants [Enn+22; NP23c;
NP23a] y sont résumés ;

* Dans le **chapitre 2** nous décrivons une généralisation du transport optimal (1) au cas où nous laissons le nombre de marginales varier et (2) au cadre quantique dans lequel on propose une approximation parcimonieuse de la fonctionnelle de Lieb. Mes co-auteurs sur ces sujets sont S. Di Marino, M. Lewin et V. Ehrlacher et nous nous concentrons sur les articles suivants [DLN22; EN23]

* Dans le **chapitre 3** nous considérons le transport optimal lorsque les dimensions des marginales (à savoir le nombre de caractéristiques dans l'exemple ci-dessus) sont différentes et certains problèmes variationnels connexes découlant de la théorie des jeux (e.g., l'équilibre de Cournot-Nash). De plus, nous proposons une classe de métriques sur l'ensemble des mesures de probabilité, impliquant des mesures singulières, qui nous aident à prouver la convergence d'une méthode numérique pour résoudre les problèmes variationnels mentionnés ci-dessus. Mon co-auteur sur ces sujets est B. Pass et nous considérons les résultats obtenus dans [NP20; NP23b].

Des perspectives sont présentées à la fin de chaque chapitre. Nous allons maintenant présenter plus précisément le contenu des différents chapitres.

Transport optimal multi-marginales et régularisation entropique

Au début de cette introduction, nous avons présenté le transport optimal classique dans le cas où nous devons comparer deux populations, à savoir les boulangeries et les cafés. Cependant, il est très naturel de traiter plusieurs populations (boulangeries, cafés, hôtels, supermarchés, etc.), de sorte qu'une manière simple de les comparer est de définir le problème de transport optimal multi-marginales suivant : étant donné $m$ distributions $\mu_i \in \mathcal{P}(X_i)$, avec $X_i \subset \mathbb{R}^{d_i}$ (pour simplifier, nous pouvons prendre $d_i = d$ pour tout $i = 1, \dots, m$), et une fonction de coût continue
$c: X \to \mathbb{R}$, le problème d'optimisation se lit comme suit

$$ \text{MOT}^c(\mu_1, \dots, \mu_m) := \inf_{\gamma \in \Pi(\mu_1, \dots, \mu_m)} \int_X c(x_1, \dots, x_m) d\gamma(x_1, \dots, x_m), \quad (2) $$

où $\Pi(\mu_1, \dots, \mu_m)$ désigne l'ensemble des mesures de probabilité sur $X := \prod_{i=1}^m X_i$ ayant $\mu_i$ comme marginales. Il s'avère que ce type de problème se pose dans de nombreux domaines d'application, y compris l'économie, [CE10], les mathématiques financières [BHP13; DS14a; DS14b; Enn+22], les statistiques [BK18; CCG16], le traitement d'images [Rab+11], la tomographie [Abr+17], l'apprentissage automatique [Haa+21; TJK22], la dynamique des fluides [Bre89] et la physique et la chimie quantique, dans le cadre de la théorie de la fonctionnelle de la densité [BDG12; CFK13]. La structure des solutions du problème de transport optimal multi-marginales est une question notoirement délicate, et n'est pas encore bien comprise, malgré les efforts considérables déployés par de nombreux chercheurs [GŚ98; Car03; CN08; Hei02; Pas11b; Pas12; KP14b; KP15; CS16; PV21; CDD15; PV22; MP17] ; voir aussi [Pas15] et [DGN17].

La classe des problèmes de transport optimal multi-marginales unidimensionnels qui se posent dans le contexte de l'estimation du risque présente un intérêt particulier. Dans ce cadre, pour des raisons de sécurité, il est nécessaire d'avoir une estimation pessimiste du risque. De tels problèmes se posent, par exemple, pour évaluer le risque d'une installation industrielle proche d'une rivière et protégée par une digue (voir [IL15]). En particulier, on sait comment évaluer un risque (tel que le niveau d'eau dans la rivière, comparé à la hauteur de la digue) en fonction de quelques variables (par exemple, la largeur de la rivière, le débit annuel maximal, etc. Ce que l'on ne sait pas, en revanche, c'est le **couplage** entre ces variables qui permet d'obtenir le scénario le plus défavorable. En gros, nous voulons résoudre le problème de maximisation suivant

$$ \max_{\gamma \in \Pi(\mu_1, \dots, \mu_m)} \mathcal{R}_{\alpha}(b_{\#\gamma}), \quad (3) $$

où $\mathcal{R}_{\alpha}(\mu) := \int_{0}^{1} F_{\mu}^{-1}(t)\alpha(t)dt$ est la mesure spectrale du risque pour une fonction non négative et non décroissante $\alpha: [0, 1] \to \mathbb{R}_{+}$ avec $\int_{0}^{1} \alpha(t)dt = 1$, $F_{\mu}^{-1}$ est la pseudo-inverse de la fonction de répartition de $\mu$ et $b$ est la fonction output (e.g., elle représente la hauteur simulée d'une rivière à risque d'inondation). Dans [Enn+22], nous parvenons (1) à montrer que pour tout problème de mesures de risque (3) peut être reformulé comme un problème MOT$^c$ avec $m+1$ marginales ; (2) à prouver que la solution peut être explicitement caractérisée en utilisant un raffinement de la théorie existante du transport multi-marginales sur un espace ambiant unidimensionnel ; (3) à étudier le cas dans lequel les variables sont en dimension supérieure ; (4) à étendre au cadre dans lequel la fonction output est multivariée.

Une façon naturelle de résoudre numériquement le problème multi-marginales présenté ci-dessus est de recourir à la régularisation entropique, ce qui nous amène au deuxième type de généralisation du transport optimal présenté dans ce manuscrit : **le transport optimal multi-marginales entropique**. Étant données comme avant $m$ mesures de probabilité $\mu_i$ et une fonction coût $c$, le problème de minimisation est défini comme suit

$$ \text{MOT}_{\varepsilon}^{c}(\mu_{1}, \dots, \mu_{m}) := \inf_{\gamma \in \Pi(\mu_{1}, \dots, \mu_{m})} \int_{X} cd\gamma + \varepsilon \text{Ent}(\gamma | \otimes_{i=1}^{m} \mu_{i}), \quad (4) $$

où $\text{Ent}(\cdot | \otimes_{i=1}^{m} \mu_i)$ est l'entropie relative de Boltzmann-Shannon (ou la divergence de Kullback-Leibler) par rapport à la mesure produit $\otimes_{i=1}^{m} \mu_i$. Grâce aux propriétés algébriques de l'entropie
la méthode des projections alternées sur les $m$ contraintes de marginal correspond, dans ce cas, au célèbre algorithme de Sinkhorn [Sin67] qui est aujourd'hui largement utilisé pour résoudre numériquement les problèmes de transport optimal et les problèmes variationnels connexes (voir par exemple [PC19]). Puisque $\text{MOT}_\varepsilon^c$ peut être considéré comme une perturbation de $\text{MOT}^c$, il est naturel d'en étudier le comportement lorsque $\varepsilon$ tend vers 0. Cela nous a motivés, voir [NP23c], à étudier le taux de convergence du coût entropique $\text{MOT}_\varepsilon^c$ vers $\text{MOT}_0^c$ sous certaines hypothèses sur les fonctions coût et les marginales. Les résultats que nous avons obtenus avec P. Pegon, et détaillés dans la section 1.4, peuvent être résumés comme suit : (1) Nous établissons deux **bornes supérieures**, l'une valable pour des coûts localement Lipschitz et une autre plus fine valable pour des coûts localement semi-concaves. Pour cette dernière, nous améliorons la borne supérieure (obtenue pour les coûts localement Lipschitz) d'un facteur de 1/2, obtenant l'inégalité suivante pour une certaine $C^* \in \mathbb{R}_+$ :

$$ \text{MOT}_{\varepsilon}^{c} \leq \text{MOT}_{0}^{c} + \frac{1}{2} \left( \sum_{i=1}^{m} d_i - \max_{1 \leq i \leq m} d_i \right) \varepsilon \log(1/\varepsilon) + C^{*} \varepsilon. \qquad (5) $$

(2) Pour la **borne inférieure**, en fonction d'un $\kappa$ dépendant d'une condition de signature sur les dérivées secondes mixtes du coût, nous avons

$$ \text{MOT}_{\varepsilon}^{c} \geq \text{MOT}_{0}^{c} + \frac{\kappa}{2} \varepsilon \log(1/\varepsilon) - C_{*} \varepsilon. \qquad (6) $$

pour une certaine $C_* \in \mathbb{R}_+$. À remarquer que quand $\kappa = \left( \sum_{i=1}^m d_i - \max_{1 \leq i \leq m} d_i \right)$ on obtient un développement à l'ordre 1 du coût entropique. Ces résultats nous ont également permis d'étendre les résultats de [CPT23; EN22] aux fonctions de coût dégénérées dans le cas des deux marginales.

La dernière section de ce chapitre 1.5 traite d'une caractérisation du problème de transport optimal multi-marginales entropique discret pour le cas d'une fonction coût pairwise, c'est-à-dire $c(x_1, \dots, x_m) = \sum_{i<j} w(x_i, x_j)$ avec $w$ symétrique, via une équation différentielle ordinaire. Des exemples de fonctions coût de ce type sont le coût Gangbo-Świech où $w(x_i, x_j) = |x_i - x_j|^2$, qui est lié au problème du barycentre de Wasserstein [AC11], et le coût de Coulomb $w(x_i, x_j) = 1/|x_i - x_j|$, qui apparaît dans le cadre de la théorie fonctionnelle de la densité [CFK13; BDG12]. Supposons maintenant que $X_i = X$, pour tout $i$, sont des ensembles finis et que les marginales sont des sommes de diracs. Dans le cas où les marginales sont toutes égales $\mu_i = \rho = \sum_{x \in X} \rho_x \delta_x$ (e.g., en chimie quantique car elles représentent des électrons indiscernables), on peut montrer que le problème discret dual de (4) prend la forme suivante (pour simplifier, nous conservons la notation du cas continu, mais les intégrales doivent être comprises comme des sommes) :

$$ \inf_{\varphi: X \to \mathbb{R}} \{ \tilde{\Phi}(\varphi, \eta) \}, \qquad (7) $$

où $\eta > 0$ est un paramètre qui permet d'ajouter/enlever certaines interactions entre les marginales,

$$ \tilde{\Phi}(\varphi, \eta) := -(m-1) \int_X \varphi d\rho + \varepsilon \int_X \log \left( \int_{X^{m-1}} \exp \left( \frac{\sum_{i=2}^m \varphi(x_i) - c_\eta(\boldsymbol{x})}{\varepsilon} \right) d \otimes^{m-1} \rho \right) d\rho, $$

et

$$ c_{\eta}(x_1, \dots, x_m) := \eta \sum_{i=2}^{m} \sum_{j=i+1}^{m} w(x_i, x_j) + \sum_{i=2}^{m} w(x_1, x_i). $$
Notez que lorsque $\eta = 1$, nous retrouvons le coût pair-wise défini ci-dessus et nous devons minimiser un seul potentiel car nous exploitons l'indiscernabilité des marginales. Ensuite, en différenciant la condition d'optimalité de (20) par rapport à $\eta$, nous obtenons le problème de Cauchy suivant

$$ \left\{ \begin{array}{l} \frac{d\varphi}{d\eta}(\eta) = -[D_{\varphi,\varphi}^2 \tilde{\Phi}(\varphi(\eta), \eta)]^{-1} \frac{\partial}{\partial\eta} \nabla_{\varphi} \tilde{\Phi}(\varphi(\eta), \eta), \\ \varphi(0) = \varphi_w, \end{array} \right. \qquad (8) $$

où la valeur initiale $\varphi(0)$ de $\varphi$ lorsque $\eta = 0$ coïncide avec le potentiel optimal $\varphi_w$ pour le problème de transport optimal à deux marginales avec le coût $w$. Après avoir prouvé que l'EDO dans (21) est bien posée, on peut l'utiliser pour caractériser le problème multi-marginales entropique original à $\eta = 1$. En exploitant ensuite la régularité du membre de droite, nous pouvons résoudre numériquement (21) en utilisant des méthodes d'ordre élevé, montrant ainsi que c'est (parfois) une approche préférable (en termes de nombre d'itérations pour converger et de temps CPU) à celle de Sinkhorn. Nous concluons le chapitre en utilisant cette approche EDO également pour le cas où la fonction coût est associée à la discrétisation temporelle du principe de Brenier pour les équations d'Euler incompressibles [Bre89].

## Transport optimal en physique mathématique

Nous continuons à explorer certaines généralisations du transport optimal, en particulier le transport multi-marginales, en nous concentrant sur l'application du transport optimal dans la théorie de la fonctionnelle de densité (DFT). L'une des quantités clés de la DFT est celle que l'on appelle la fonctionnelle de Levy-Lieb, qui est l'énergie minimale (cinétique plus interaction électron-électron) d'un système quantique ayant $\rho$ comme densité électronique prescrite. En d'autres termes, nous devons résoudre le problème de minimisation suivant

$$ F_{LL}[\rho] := \inf_{\substack{\Psi \in \mathcal{H}_1^m \\ \rho_{\Psi} = \rho}} \frac{1}{2} \int_{\mathbb{R}^{3m}} |\nabla \Psi|^2 + \int_{\mathbb{R}^{3m}} V |\Psi|^2, \qquad (9) $$

où $\mathcal{H}_1^m$ est l'ensemble des fonctions d'onde antisymétriques $\Psi$ pour un système de $m$ électrons avec une énergie cinétique finie, la fonction $V$ est le potentiel d'interaction coulombienne électron-électron, et $\rho_{\Psi}$ est la densité électronique associée à la fonction d'onde $\Psi$, c'est-à-dire la fonction réelle définie sur $\mathbb{R}^3$ comme suit

$$ \forall x \in \mathbb{R}^3, \quad \rho_{\Psi}(x) := m \int_{(\mathbb{R}^3)^{m-1}} |\Psi(x, x_2, \dots, x_m)|^2 dx_2 \dots dx_m. $$

Il se trouve que $F_{LL}$ n'est pas convexe, il est donc pratique de considérer une convexification proposée par Lieb [Lie83], où la minimisation est effectuée sur l'ensemble des états mixtes plutôt que sur l'ensemble des états purs comme dans (9): nn gros, les fonctions d'onde sont maintenant remplacées par des opérateurs et l'énergie dans (9) est exprimée en termes de traces. Plus précisément, nous considérons le problème de minimisation sur l'ensemble $\mathcal{S}_1^+$ des opérateurs auto-adjoints non négatifs de trace 1 ayant $\rho$ comme marge, c'est-à-dire

$$ F_L[\rho] := \inf_{\substack{\Gamma \in \mathcal{S}_1^+ \\ \rho_{\Gamma} = \rho}} \mathrm{Tr}(H_m \Gamma), \qquad (10) $$
où $H_m := -\frac{1}{2}\Delta + V$. Remarquez également que l'opérateur optimal $\Gamma$ peut être diagonalisé, c'est-à-dire qu'il existe $(\Psi_i)_{i \in \mathbb{N}^*}$ et $(\alpha_i)_{i \in \mathbb{N}^*}$ tels que $\Gamma = \sum_i \alpha_i \Psi_i$. On considère maintenant un paramètre semi-classique $\epsilon = \hbar^2$, en faisant un rescaling de la densité $\rho$, on peut étudier la limite lorsque $\epsilon \to 0$ à la fois de (22) et de (23). Les fonctionnelles ci-dessus convergent (voir [Lew18] pour la fonctionnelle de Lieb et [CFK13; CFK18; BD17] pour celle de Lévy-Lieb) vers le problème de transport optimal multi-marginales avec un coût de Coulomb

$$ \inf_{\gamma \in \Pi(\rho)} \int_{\mathbb{R}^{3m}} \sum_{i<j} \frac{1}{|x_i - x_j|} d\gamma \quad (11) $$

où maintenant les marginales sont toutes égales à la densité électronique donnée $\rho$. Le MOT$^c$ avec un coût de Coulomb est donc une manière rigoureuse d'approcher le terme de répulsion électron-électron dans la fonctionnelle de Levy-Lieb/Lieb. Pour cette raison, (11) joue un rôle central dans la DFT computationnelle, et trouver des moyens de la calculer (ou de l'approcher) de manière efficace est crucial. Ce chapitre est en quelque sorte consacré à trouver des généralisations de MOT pour approcher (11) d'une manière plus facile à traiter (numériquement).

La première partie de ce chapitre 2.2 est consacrée à une généralisation du MOT$^c$ qui approxime la répulsion électron-électron dans la fonctionnelle de Levy-Lieb en donnant des informations sur l'importance des fluctuations des corrélations dans le système. En particulier, nous permettons maintenant au nombre de marginales $m$ de varier. Ce modèle est inspiré de la physique statistique, où il est généralement connu sous le nom de **grand-canonique** [Rue99]. Le problème du Transport Optimal Grand Canonique (GC-OT) peut être formulé sous la forme suivante :

$$ \mathrm{GCOT}_{0}^{c}(\rho) := \inf \left\{ \sum_{n=0}^{\infty} \int_{X^n} c_n \, d\gamma_n : \sum_{n=0}^{\infty} \gamma_n(X^n) = 1, \quad \sum_{n=1}^{\infty} n \gamma_n(\cdot, X^{n-1}) = \rho \right\} \quad (12) $$

où chaque $\gamma_n$ est une mesure symétrique sur $X^n$, $c_n$ est le coût symétrique pour le problème $n$-marginal (par exemple, le potentiel de Coulomb) et $\gamma_n(\cdot, X^{n-1})$ est la première marge de $\gamma_n$. Remarquez le facteur $n$ multipliant la marge dans la contrainte impliquant $\rho$, ce qui tient compte du fait qu'il y a $n$ marges égales de ce genre puisque toutes les $\gamma_n$ sont symétriques. La famille $\gamma = (\gamma_n)_{n \ge 0}$ forme une probabilité qui décrit le comportement de certains agents dont le nombre est inconnu ou peut varier. Dans cette interprétation, $\gamma_n(X^n)$ est la probabilité qu'il y ait $n$ agents et $\gamma_0 \in [0, 1]$ est celle qu'il n'y ait aucun agent du tout. Dans le problème GC-OT (12), seule la quantité moyenne $\rho$ est fixée et des fluctuations du nombre d'agents sont autorisées. Résoudre le problème de minimisation $\mathrm{GCOT}_0^c(\rho)$ nécessite en particulier de déterminer la meilleure façon de répartir le nombre d'agents à travers les mesures $\gamma_n$, afin de reproduire la moyenne donnée $\rho$, en fonction des coûts correspondants $c = (c_n)n \ge 0$. Il est à noter que, dans le cas où nous considérons le potentiel de Coulomb comme $c_n$, résoudre ce problème est censé fournir des informations sur les corrélations dans le système, en fonction de la densité $\rho$, qui pourraient en suite être utilisées en DFT computationnelle. À remarquer que la contrainte de symétrie sur les mesures $\gamma_n$ ne nous empêche pas de traiter des systèmes avec différents types d'agents. Revenons à l'exemple des boulangeries/cafés : si nous voulons transporter une certaine quantité de pains au chocolat/chocolatines, alors $X$ sera un ensemble fini contenant les propriétés (emplacement précis, taille, etc.) à la fois des boulangeries et des cafés. La symétrie signifie simplement que les pains au chocolat/chocolatines sont tous identiques et que nous ne voulons pas distinguer lequel est envoyé où. La principale différence avec l'approche standard est que maintenant nous
pouvons définir un problème de transport sans connaître exactement ni la quantité de pains au chocolat/chocolatines produits ni la demande. Nous discutons plusieurs propriétés mathématiques du problème GT-OT (12), dans le cadre de la théorie du transport optimal. (1) Nous formulons le problème et montrons l'existence d'un minimiseur $\gamma = (\gamma_n)_{n \ge 0}$, sous des hypothèses sur les coûts. (2) Nous déduisons certaines propriétés sur le support de $\gamma$, c'est-à-dire sur le nombre de $\gamma_n$ qui sont non nuls (nous discutons en particulier le cas du potentiel de Coulomb). (3) Nous étudions également le problème tronqué où toutes les $\gamma_n$ sont supposées s'annuler après un certain $N_{max}$, et la convergence vers le problème réel lorsque $N_{max} \to \infty$. Cela est utile pour effectuer des calculs numériques. (4) Nous nous concentrons sur la théorie de la dualité et l'existence du potentiel dual. (5) Nous étudions le cas de la dimension $d=1$, confirmant ainsi certaines conjectures sur la forme du plan optimal formulées dans [MSG13]. (6) Enfin, nous étudions la régularisation entropique du problème GC-OT.

Dans la dernière partie du chapitre 2.3, nous nous concentrons sur la fonctionnelle de Lieb (23), qui peut être considérée en quelque sorte comme une version quantique d'un problème standard de MOT. À partir de ces observations, nous visons à utiliser certains résultats récents sur l'approximation par moments du MOT. Plus précisément, l'approche d'abord considérée dans [Alf+21] consiste à introduire une approximation, lorsque le nombre de moments $M$ tend vers l'infini, des problèmes de transport multi-marginaux où les contraintes de marginales sont remplacées par des contraintes de moment associées à $M$ "fonctions de moment" qui sont des fonctions à valeurs réelles définies sur $\mathbb{R}^3$.

La solution de ce problème de transport optimal avec contrainte des moments est toujours une mesure de probabilité définie sur $\mathbb{R}^{3m}$, mais elle a également une structure parcimonieuse au sens où elle peut être écrite sous la forme d'une mesure discrète concentrée sur un certain nombre de points appartenant à $\mathbb{R}^{3m}$. Ce nombre de points évolue au plus de manière linéaire avec le nombre de contraintes de moments. Dans notre cadre, cela se traduit par la démonstration que les solutions de l'approximation de (10) par contraintes de moments peuvent être écrites sous la forme $\Gamma = \sum_{k=1}^K \alpha_k |\Psi_k\rangle\langle\Psi_k|$, où $K \in \mathbb{N}^*$ évolue au plus de manière linéaire avec le nombre de contraintes de moments et $\Psi_k$ sont des fonctions d'onde. En particulier, nous sommes en mesure de prouver (1) l'existence de minimiseurs pour le problème approché ainsi que la structure parcimonieuse. (2) Nous montrons la convergence de l'approximation vers la fonctionnelle de Lieb exacte. De plus, nous prouvons également le taux de convergence de l'approximation de l'énergie de l'état fondamental associée vers la valeur exacte. (3) Enfin, nous présentons quelques résultats sur la formulation duale du problème approché, en remarquant en particulier que le problème dual peut être reformulé comme un problème de programmation semi-définie positive.

## Transport optimal entre dimensions différentes

Le dernier chapitre de ce manuscrit traite du cas où le nombre de propriétés des boulangeries et des cafés diffère, c'est-à-dire $d_1 \neq d_2$, ce qui conduit à une autre généralisation du transport optimal, bien étudiée dans [CMP17; MP20], connue sous le nom de **transport optimal entre dimensions différentes**. Nous nous concentrons notamment sur l'étude de certains problèmes variationnels impliquant ce type de terme de transport optimal :

$$ \inf_{\mu_1, \mu_2} \mathcal{J}(\mu_1, \mu_2) := \mathrm{OT}^c(\mu_1, \mu_2) + \mathcal{F}(\mu_2) + \mathcal{G}(\mu_1). \qquad (13) $$

On s'intéresse à la caractérisation des minimiseurs de $(\mu_1, \mu_2) \mapsto \mathcal{J}(\mu_1, \mu_2)$, ainsi que les
minimiseurs des sous-problèmes obtenus lorsque soit $\mu_1$ soit $\mu_2$ est fixé : $\mu_2 \mapsto \mathcal{J}(\mu_1, \mu_2)$ et $\mu_1 \mapsto \mathcal{J}(\mu_1, \mu_2)$. Des problèmes de ces formes, pour divers choix des fonctionnelles $\mathcal{F}$ et $\mathcal{G}$, se posent dans une grande variété d'applications, notamment : les flots de gradient sur l'espace de Wasserstein (où la minimisation $\mu_2 \mapsto \mathcal{J}(\mu_1, \mu_2)$ représente une étape dans un flot de gradient discret en temps, JKO [JKO98]). Ou encore l'interpolation par déplacement (lorsque $\mathcal{F}$ est la distance de Wasserstein par rapport à une deuxième mesure de probabilité), les équilibres de Cournot-Nash en théorie des jeux, la planification urbaine, etc. Pour clarifier la raison de la prise en compte d'un terme OT de dimensions inégales, considérons le sous-problème où la distribution $\mu_1$, disons la distribution des boulangeries, est fixée et où une personne très riche, ayant le désir incontrôlable d'ouvrir de nombreux cafés le long de la ligne de métro la plus fréquentée de Paris, veut trouver la distribution des cafés qui minimise $\mu_2 \mapsto \text{OT}^c(\mu_1, \mu_2) + \mathcal{F}(\mu_2)$ : le coût de transport des pains au chocolat/chocolatines plus une fonctionnelle qui, par exemple, évite une concentration trop élevée de cafés au même endroit. Maintenant, en modélisant le fait que $\mu_2$ est concentré sur une courbe (la ligne de métro), on assume que l'espace $X_2$ a une dimension plus basse que $X_1$, par exemple $d_2 = 1$.

La première section du chapitre 3.1 est entièrement consacrée à l'étude des minimiseurs de tels problèmes variationnels. Remarquez que le problème OT entre dimensions différentes est traitable (c'est-à-dire que nous pouvons construire une solution presque explicite) lorsqu'une certaine condition, dite de **nestedness** (emboîtement), concernant le coût $c$ et les marginales $\mu_1$ et $\mu_2$ est satisfaite. Malheureusement, dans le contexte actuel, seulement le coût et **une** des marginales (et **aucune** dans le cas de doubles minimisations) sont connus. Nous prouvons ici que, sous diverses hypothèses sur $c$, $\mu_1$ et $X_2$, $(c, \mu_1, \mu_2)$ est nested chaque fois que $\mu_2$ minimise $\mu_2 \mapsto \mathcal{J}(\mu_1, \mu_2)$, pour différentes choix de la fonctionnelle $\mathcal{F}$ ; des résultats analogues sont également établis pour certaines formes spécifiques de $\mathcal{G}$ pour les minimisations dans l'espace de dimension supérieure, $\mu_1 \mapsto \mathcal{J}(\mu_1, \mu_2)$, et pour les doubles minimisations $(\mu_1, \mu_2) \mapsto \mathcal{J}(\mu_1, \mu_2)$. Nous montrons ensuite que cette garantie a priori de nestedness rend le problème de caractérisation ou d'identification des minimiseurs beaucoup plus traitable ; dans différents contextes, en fonction de la forme précise de $\mathcal{F}$, nous établissons que les solutions peuvent être caractérisées par des équations différentielles ( locales), peuvent être calculées numériquement par un schéma itératif convergent, ou peuvent être obtenues sous une forme presque explicite.

La deuxième partie de ce chapitre 3.2 traite d'une classe de métriques sur l'espace des mesures de probabilité impliquant des mesures singulières. Plus précisément, nous développons la théorie d'une métrique, que nous appelons la métrique de Wasserstein avec base $\nu$ et que nous notons $W_\nu$, sur l'ensemble des mesures de probabilité $\mathcal{P}(X)$ sur un domaine $X \subseteq \mathbb{R}^d$. Cette métrique est basée sur un léger affinement de la notion de **géodésiques généralisées**, voir [AGS04], par rapport à une mesure de base $\nu$ et est pertinente en particulier dans le cas où $\nu$ est singulière par rapport à la mesure de Lebesgue $d$-dimensionnelle ; elle est également liée au concept de transport optimal linéarisé [Wan+13]. La métrique $W_\nu$ est définie en termes d'un problème variationnel impliquant le transport optimal vers $\nu$ ; nous la caractérisons également en termes d'intégrations de la métrique de Wasserstein classique (transport optimal avec le coût quadra-tique) entre les probabilités conditionnelles lorsque les mesures sont désintégrées par rapport au transport optimal vers $\nu$, et à travers la limite de certains problèmes de transport optimal multi-marginales. En revenant au problème d'appariement entre les boulangeries et les cafés. L'idée du $W_\nu$ consiste à apparier les deux populations en tenant compte des informations fournies par la mesure de base $\nu$. Cela est particulièrement clair lorsqu'on considère la caractérisation en termes d'intégrations de la métrique de Wasserstein classique entre les probabilités conditionnelles :
considérez $\nu$ comme une mesure de probabilité unidimensionnelle telle que, par exemple, pour un
quartier/arrondissement $y$ de Paris (elle est unidimensionnelle car elle pourrait simplement être
un ensemble d'étiquettes), $\nu(y)$ vous indique à quel point il est important d'apparier les boulan-
geries et les cafés à proximité. Ensuite, l'appariement via la métrique $W_\nu$ fonctionne comme suit
: (1) identifiez les boulangeries et les cafés proches d'un certain $y$ (la désintégration par rapport
au transport optimal vers $\nu$), (2) appariez-les en utilisant la distance de Wasserstein classique.
À mesure que nous faisons varier la mesure de base $\nu$, la métrique $W_\nu$ interpole entre la
métrique de Wasserstein quadratique habituelle (obtenue lorsque $\nu$ est une mesure de Dirac) et
une métrique associée aux uniques géodésiques généralisées définies lorsque $\nu$ est suffisamment
régulière (par exemple, absolument continue par rapport à la mesure de Lebesgue). Lorsque $\nu$ se
concentre sur une sous-variété de dimension inférieure de $\mathbb{R}^d$, nous démontrons que le problème
variationnel dans la définition de la métrique $W_\nu$ a une solution unique et nous établissons la
convexité géodésique de la classe habituelle de fonctionnelles. Nous introduisons également une
classe de métriques définies par rapport à une mesure de base fixée $\mu$ de dimension supérieure,
sur l'ensemble des mesures qui sont absolument continues par rapport à une deuxième mesure
de base fixée $\sigma$. À l'aide de cette métrique, nous prouvons la convergence d'un schéma itératif
pour résoudre un problème variationnel étudié dans la première section du chapitre.

Permettez-moi de conclure cette introduction en mentionnant quelques travaux qui ne sont pas
résumés dans ce manuscrit : [Ben+18] ( étude des jeux à champ moyen variationnels du second
ordre et lien avec la minimisation d'une entropie relative), [DGN23] (analyse des minimiseurs de la
fonctionnelle de Levy-Lieb dans le cas bosonique) et [NP22] (méthodes numériques pour résoudre
des problèmes variationnels impliquant le transport optimal entre dimensions différentes).
# Introduction (in English)

Optimal Transportation is a very dynamic research field and there are so many examples and applications that it is actually difficult to choose one in order to introduce this fascinating subject; we refer the reader to some “old and new” textbooks such as [AGS04; Vil09; Vil03; San15; Gal16]. Here, we have decided to introduce the classical optimal transport problem by using an example given by Villani in [Vil03].

Imagine that we are in Paris and consider the bakeries, producing pains au chocolat/chocolatines (the reader is free to choose his favorite French pastry), which should be shipped each morning to the local cafés where customers will enjoy eating them. The production and the consumption of pains au chocolat/chocolatines are described by a distribution $\mu_1$ and $\mu_2$, respectively; we assume that the total amount of the production and the consumption is the same. Then, the Optimal Transportation problem consists in transporting the amount of pains au chocolat/chocolatines produced by the bakery $x_1 \in X_1$ ($X_1 \subset \mathbb{R}^{d_1}$ is the set of bakeries differentiated by $d_1$ properties e.g., precise location, size, etc) to a café $x_2 \in X_2$ ($X_2 \subset \mathbb{R}^{d_2}$ is the set of cafés differentiated by $d_2$ properties) such that the transport cost $c(x_1, x_2)$, for instance the distance between the bakery and the café, is minimized. In classical optimal transport we have that $d_1 = d_2$, but in the following we will try to tackle the problem when the number of characteristics (a.k.a. the dimensions) differs. The problem introduced by Monge in [Mon81] consists in searching for a transport map $x_2 = T(x_1)$ which tells us the café $x_2$ where all the amount of pains au chocolat/chocolatines produced by $x_1$ is sent. However this problem is difficult to treat and a relaxation (known as the **Monge-Kantorovich problem**) is needed: we allow to split the “mass” at $x_1$ such that, as an example, a fraction is sent to $x_2$ and the other one to $x'_2$. So, in this second case we are looking for an optimal coupling $\gamma(x_1, x_2)$ which tells us how mass at $x_1$ is distributed to all $x_2 \in X_2$. The problem can be then formulated as follows: given two probability measures $\mu_1 \in \mathcal{P}(X_1)$ and $\mu_2 \in \mathcal{P}(X_2)$ and a continuous cost function $c: X_1 \times X_2 \to \mathbb{R}$ the **Monge-Kantorovich problem** reads as

$$ \mathrm{OT}^c(\mu_1, \mu_2) := \inf_{\gamma \in \Pi(\mu_1, \mu_2)} \int_{X_1 \times X_2} c(x_1, x_2) d\gamma(x_1, x_2), \quad (14) $$

where $\Pi(\mu_1, \mu_2)$ denotes the set of probability measures on $X_1 \times X_2$ having $\mu_1$ and $\mu_2$ as marginals, that is, given the canonical projection $e_i: X_1 \times X_2 \to X_i$, we have $e_{i,\#} \gamma = \mu_i$ for $i=1,2$. In the case in which the “cheapest” coupling $\gamma(x_1, x_2)$ is to assign a unique café to each bakery $x_1$, then we say that the optimal coupling is *deterministic* and an optimal map exists. When cost $c(x_1, x_2)$ is the squared distance, the existence of such a deterministic coupling was proved by Yann Brenier in his seminal work [Bre91]; moreover in this case (14) is a squared distance over the set of probability measures.

Many questions arise at this point: what happens if the number of characteristics of bakeries and cafés is different (e.g., $d_1 \neq d_2$)? and if one of the distributions is unknown and we wish to
obtain it via minimizing the OT^c plus another functional? Moreover, we can also add a third category, e.g., the hotels, so that now we have to minimize the cost of sending pains au chocolat/chocolatines to both the cafés and the hotel in Paris. Can this problem still be interpreted as an OT problem? This manuscript (which partially summarizes my research after the Ph.D.) is, actually, devoted to answering these questions and a few other ones dealing with optimal transport models arising in Physics, especially in Quantum Physics.

More precisely, this manuscript is divided in three chapters, thematically, organized as follows:

* **Chapter 1** is devoted to (1) multi-marginal optimal transport and an application to risk measures, (2) convergence rates for entropic multi-marginal optimal transport and, finally, (3) an ODE characterization of (discrete) entropic multi-marginal optimal transport. My co-authors on these topics are H. Ennaji, Q. Mérigot, B. Pass and P. Pegon and the results of the following papers [Enn+22; NP23c; NP23a] are summarized;

* In **Chapter 2** we describe a generalization of optimal transport to (1) the case in which we let the number of marginals vary and (2) to the quantum framework deriving a sparse approximation of the Lieb Functional. My co-authors on these topics are S. Di Marino, M. Lewin and V. Ehrlacher and we focus on the following papers [DLN22; EN23]

* In **Chapter 3** we consider the optimal transport when the dimensions of the marginals (namely the number of characteristics in the example above) are different and some related variational problems arising in Game Theory (e.g., the Cournot-Nash equilibrium). Moreover, we propose a class of metrics on the set of probability measures, involving singular measures, which help us to prove convergence of a numerical method to solve the above mentioned variational problems. My co-author on these topics is B. Pass and we consider the results obtained in [NP20; NP23b].

At the end of each chapter we present some perspectives. We now introduce more precisely the content of the different chapters.

## Multi-Marginal Optimal Transport and Entropic Regularization

At the beginning of the introduction we presented the standard optimal transport in the case in which we have to match/compare two populations, namely the bakeries and the cafés. However, it is very natural to deal with many populations (e.g., bakeries, cafés, hotels, supermarkets, etc) so that a simple way to compare them is to define the following **multi-marginal optimal transport problem**: given $m$ distributions $\mu_i \in \mathcal{P}(X_i)$, with $X_i \subset \mathbb{R}^{d_i}$ (for simplicity we can take $d_i = d$ for all $i = 1, \dots, m$), and a continuous cost function $c: \mathbf{X} \to \mathbb{R}$ the optimization problem reads as

$$ \text{MOT}^c(\mu_1, \dots, \mu_m) := \inf_{\gamma \in \Pi(\mu_1, \dots, \mu_m)} \int_{\mathbf{X}} c(x_1, \dots, x_m) d\gamma(x_1, \dots, x_m), \quad (15) $$

where $\Pi(\mu_1, \dots, \mu_m)$ denotes the set of probability measures on $\mathbf{X} := \prod_{i=1}^m X_i$ having $\mu_i$ as marginals. It turns out that this kind of problem arises in many areas of applications including economics [CE10], financial mathematics [BHP13; DS14a; DS14b; Enn+22], statistics [BK18; CCG16], image processing [Rab+11], tomography [Abr+17], machine learning [Haa+21; TJK22], fluid dynamics [Bre89] and quantum physics and chemistry, in the framework of density functional
theory [BDG12; CFK13]. The structure of solutions to the multi-marginal optimal transport problem is a notoriously delicate issue, and is still not well understood, despite substantial efforts by many researchers [GŚ98; Car03; CN08; Hei02; Pas11b; Pas12; KP14b; KP15; CS16; PV21; CDD15; PV22; MP17]; see also the surveys [Pas15] and [DGN17].

Of particular interest is the class of one dimensional multi-marginal optimal transport problems arising in the context of **risk estimation**. In this framework, for safety reasons, one needs to have a pessimistic estimation of the risk. Such problems arise for evaluating the risk of an industrial facility close to a river and protected by a dyke (see [IL15]). In particular, one knows how to evaluate a risk (such as the level of water in the river, compared to the height of the dyke) depending on a few variables (e.g. river width, maximal annual flowrate, etc.) whose probability distribution is known. What is unknown however, is the **coupling** between these variables yielding the worst-case scenario. Roughly speaking we want to solve the following maximization problem

$$ \max_{\gamma \in \Pi(\mu_1, \dots, \mu_m)} \mathcal{R}_{\alpha}(b_{\#}\gamma), \qquad (16) $$

where $\mathcal{R}_{\alpha}(\mu) := \int_{0}^{1} F_{\mu}^{-1}(t)\alpha(t)dt$ is the spectral risk measure for a non-negative and non-decreasing function $\alpha: [0,1] \to \mathbb{R}_{+}$ with $\int_{0}^{1} \alpha(t)dt = 1$, $F_{\mu}^{-1}$ is the pseudo-inverse function of $\mu$ and $b$ is the output function (e.g., it represents the simulated height of a river at risk of flooding). In [Enn+22] we (1) show that for any risk measures problem (16) can be reformulated as a MOT$^c$ problem with $m+1$ marginals; (2) prove that the solution can be explicitly characterized by using a careful refinement of the existing theory of multi-marginal transport over one dimensional ambient space; (3) study the case in which the underline variables are in higher dimension; (4) extend to the framework in which the output function is multivariate.

A natural way to solve numerically the multi-marginal problem presented above is by means of the well known entropic regularization which brings us to the second kind of generalization of optimal transport presented in this manuscript: **the entropic multi-marginal optimal transport problem**. Given as before $m$ probability measures $\mu_i$ and a cost function $c$ the minimization problem is defined as

$$ \mathrm{MOT}_{\varepsilon}^{c}(\mu_{1}, \dots, \mu_{m}) := \inf_{\gamma \in \Pi(\mu_{1}, \dots, \mu_{m})} \int_{\mathcal{X}} cd\gamma + \varepsilon \mathrm{Ent}(\gamma | \otimes_{i=1}^{m} \mu_{i}), \qquad (17) $$

where $\mathrm{Ent}(\cdot | \otimes_{i=1}^{m} \mu_i)$ is the Boltzmann-Shannon relative entropy (or Kullback-Leibler divergence) w.r.t. the product measure $\otimes_{i=1}^{m} \mu_i$. Notice that when it comes to solve the regularized problem by using the alternating Kullback-Leibler projections on the $m$ marginals constraints, by the algebraic properties of the entropy such an algorithm corresponds to the celebrated Sinkhorn's algorithm [Sin67] which is nowadays widely used to numerically solve optimal transport and related variational problems (see for instance [PC19]). Now, since $\mathrm{MOT}_{\varepsilon}^{c}$ can be seen as perturbation of $\mathrm{MOT}^{c}$ it is natural to study the behaviour as $\varepsilon$ vanishes. This motivated us, see [NP23c], to investigate the rate of convergence of the entropic cost $\mathrm{MOT}_{\varepsilon}^{c}$ to $\mathrm{MOT}_{0}^{c}$ under some mild assumptions on the cost functions and the marginals. The results we obtained together with P. Pegon, and detailed in Section 1.4, can be summarized as follows:

(1) we establish two **upper bounds**, one valid for locally Lipschitz costs and a finer one valid for locally semi-concave costs. For the latter, we improve the upper bound by a $1/2$ factor,
obtaining the following inequality for some $C^* \in \mathbb{R}_+$

$$ \mathrm{MOT}_{\varepsilon}^{c} \leq \mathrm{MOT}_{0}^{c} + \frac{1}{2} \left( \sum_{i=1}^{m} d_{i} - \max_{1 \leq i \leq m} d_{i} \right) \varepsilon \log(1/\varepsilon) + C^{*} \varepsilon. \quad (18) $$

(2) For the lower bound, given a $\kappa$ depending on a signature condition on the second mixed derivatives of the cost, we have

$$ \mathrm{MOT}_{\varepsilon}^{c} \geq \mathrm{MOT}_{0}^{c} + \frac{\kappa}{2} \varepsilon \log(1/\varepsilon) - C_{*} \varepsilon. \quad (19) $$

for some $C_* \in \mathbb{R}_+$.

These also helped us to extend the results in [CPT23; EN22] to degenerate cost function in the two marginals case. Notice that when $\kappa = \left( \sum_{i=1}^{m} d_i - \max_{1 \le i \le m} d_i \right)$ we obtain a first order development of the entropic cost.

The last section of this chapter 1.5 deals with a characterisation of the discrete entropic multi-marginal optimal transport problem for the case of pair-wise cost function, that is $c(x_1, \dots, x_m) = \sum_{i<j} w(x_i, x_j)$ with $w$ is symmetric, via an ordinary differential equation. Typical examples of pair-wise cost functions are the Gangbo-Święch cost where $w(x_i, x_j) = |x_i - x_j|^2$, which is related to the Wasserstein barycenter problem [AC11], and the Coulomb cost $w(x_i, x_j) = 1/|x_i - x_j|$, which appears in the framework of Density Functional theory [CFK13; BDG12]. Assume now that $X_i = X$, for all $i$, are finite sets and the marginals are sum of diracs, one can show that if the marginals are all equal $\mu_i = \rho = \sum_{x \in X} \rho_x \delta_x$ (this is the case of the Coulomb cost since they represent electrons which are indistinguishable) then the discrete problem dual to (17) takes the following simple form (for simplicity we keep the continuous case notation but integrals must be intended as sums)

$$ \inf_{\varphi: X \to \mathbb{R}} \{ \tilde{\Phi}(\varphi, \eta) \}, \quad (20) $$

where

$$ \tilde{\Phi}(\varphi, \eta) := -(m-1) \int_X \varphi d\rho + \varepsilon \int_X \log \left( \int_{X^{m-1}} \exp \left( \frac{\sum_{i=2}^m \varphi(x_i) - c_\eta(\boldsymbol{x})}{\varepsilon} \right) d \otimes^{m-1} \rho \right) d\rho, $$

$$ c_{\eta}(x_1, \dots, x_m) := \eta \sum_{i=2}^{m} \sum_{j=i+1}^{m} w(x_i, x_j) + \sum_{i=2}^{m} w(x_1, x_i), $$

and $\eta > 0$ is a parameter which allows to interpolate between the original and a special cost function for which the problem's complexity scales linearly in $m$. Notice that when $\eta = 1$ we retrieve the pair-wise cost defined above and we have to minimize just over one potential since we exploit the undistinguishability of the marginals. Now, differentiating the optimality condition of (20) with respect to $\eta$ we obtain the following Cauchy problem:

$$ \left\{ \begin{array}{l} \frac{d\varphi}{d\eta}(\eta) = -[D_{\varphi,\varphi}^2 \tilde{\Phi}(\varphi(\eta), \eta)]^{-1} \frac{\partial}{\partial\eta} \nabla_{\varphi} \tilde{\Phi}(\varphi(\eta), \eta), \\ \varphi(0) = \varphi_w, \end{array} \right. \quad (21) $$
where, the initial value $\varphi(0)$ of $\varphi$ when $\eta = 0$ coincides with the optimal potential $\varphi_w$ for the two marginal optimal transport problem with cost $w$. After proving the well-posedness of ODE in (21) one can use it in order to characterize the original entropic multi-marginal problem obtained for $\eta = 1$. By exploiting then the regularity of the rhs, we can solve numerically (21) by using some high order methods showing that it is (sometimes) a preferable approach (in terms of number of iterations to converge and CPU times) to Sinkhorn's. We conclude the chapter by using this ODE approach also for the case in which the cost function associated to the time discretization of Brenier's principle for incompressible Euler equations [Bre89].

## Optimal Transport in mathematical physics

We study some other generalizations of optimal transport, and in particular multi-marginal transport, by focusing on the application of OT in Density Functional Theory (DFT). One of the key quantities in DFT is the so-called Levy-Lieb functional which is the least energy (kinetic plus electron-electron interaction) of a quantum system having the prescribed electron density $\rho$. Namely we have to solve the following minimization problem

$$F_{LL}[\rho] := \inf_{\substack{\Psi \in \mathcal{H}_1^m \\ \rho_{\Psi} = \rho}} \frac{1}{2} \int_{\mathbb{R}^{3m}} |\nabla \Psi|^2 + \int_{\mathbb{R}^{3m}} V |\Psi|^2, \quad (22)$$

where $\mathcal{H}_1^m$ is the set of antysymmetric wavefunctions $\Psi$ for a system of $m$ electrons with finite kinetic energy, the function $V$ is the electron-electron Coulomb interaction potential, $\rho_{\Psi}$ is the electronic density associated to the wavefunction $\Psi$, namely the real-valued function defined over $\mathbb{R}^3$ as follows:

$$\forall x \in \mathbb{R}^3, \quad \rho_{\Psi}(x) := m \int_{(\mathbb{R}^3)^{m-1}} |\Psi(x, x_2, \dots, x_m)|^2 dx_2 \dots dx_m.$$

It happens that $F_{LL}$ is not convex so it is convenient to look at a convexification proposed by Lieb [Lie83] where the minimization is performed over the set of mixed states instead of the set of pure ones as in (22): roughly speaking the wavefunctions are now replaced by operators and the energy in (22) is expressed in terms of traces. More precisely, we consider the minimization problem over the set $\mathcal{S}_1^+$ of trace-class self-adjoint non-negative operators having $\rho$ as marginal, that is

$$F_L[\rho] := \inf_{\substack{\Gamma \in \mathcal{S}_1^+ \\ \rho_{\Gamma} = \rho}} \mathrm{Tr}(H_m \Gamma), \quad (23)$$

where $H_m := -\frac{1}{2}\Delta + V$. Notice also that the optimal $\Gamma$ can be diagonalized, that is there exists $(\Psi_i)_{i \in \mathbb{N}^*}$ and $(\alpha_i)_{i \in \mathbb{N}^*}$ such that $\Gamma = \sum_i \alpha_i \Psi_i$. Consider now an effective semi-classical parameter $\epsilon = \hbar^2$, by scaling the density $\rho$, and look at the limit for $\epsilon \to 0$ of both (22) and (23). The above functionals converge (see [Lew18] for the Lieb functional and [CFK13; CFK18; BD17] for the Lévy-Lieb one) to the multi-marginal optimal transport problem with Coulomb cost

$$\inf_{\gamma \in \Pi(\rho)} \int_{\mathbb{R}^{3m}} \sum_{i<j} \frac{1}{|x_i - x_j|} d\gamma \quad (24)$$

where now the marginals are all equal to the given electron density $\rho$. The MOT$^c$ with Coulomb cost is then a rigorous way to approximate the electron-electron repulsion term in Levy-Lieb/Lieb
functional. For this reason, (24) plays a central role in computational DFT and so finding ways to efficiently compute (or approximate) is crucial. This chapter is in some sense devoted to find generalizations of MOT to approximate (24) in a more (numerical) tractable way.

The first part of this chapter 2.2 is devoted to a generalization of MOT<sup>c</sup> which approximates the electron-electron repulsion in Levy-Lieb functional by giving some information on the amount of fluctuations of correlations in the system. In particular we now allow the number of marginals m to vary. This model takes its root in statistical physics, where it usually goes under the name **grand-canonical** [Rue99]. In short, the Grand-Canonical Optimal Transport (GC-OT) problem can be formulated in the Kantorovich form as follows

$$ \mathrm{GCOT}_{0}^{c}(\rho) := \inf \left\{ \sum_{n=0}^{\infty} \int_{X^n} c_n \, d\gamma_n : \sum_{n=0}^{\infty} \gamma_n(X^n) = 1, \quad \sum_{n=1}^{\infty} n \gamma_n(\cdot, X^{n-1}) = \rho \right\} \quad (25) $$

where each $\gamma_n$ is a symmetric measure on $X^n$, $c_n$ is the symmetric cost for the n-marginal problem (for instance the Coulomb potential) and $\gamma_n(\cdot, X^{n-1})$ is the first marginal of $\gamma_n$. Notice the factor n multiplying the marginal in the constraint involving $\rho$, which accounts for the fact that there are n equal such marginals since all the $\gamma_n$ are symmetric. The family $\gamma = (\gamma_n)_{n \ge 0}$ forms a probability which describes the behavior of some agents whose number is unknown or can vary. In this interpretation $\gamma_n(X^n)$ is the probability that there are n agents and $\gamma_0 \in [0, 1]$ is the one that there is no agent at all. In the GC-OT problem (25) only the *average* quantity $\rho$ is fixed and fluctuations of the number of agents are allowed. Solving the minimization problem $\mathrm{GCOT}_0^c(\rho)$ requires in particular to determine the best way to distribute the number of agents through the measures $\gamma_n$, in order to reproduce the given average $\rho$, depending on the corresponding costs $c = (c_n)_{n \ge 0}$. It is worth noticing that, in the case in which we consider the Coulomb potential as $c_n$, solving this problem is believed to give some information on the amount of correlations in the system, depending on the density $\rho$, which could then be used in the approximation of the more precise quantum problem. Let us mention that the symmetric constraint on the measures $\gamma_n$ does not prevent us from treating systems with different kinds of agents. Going back to the bakeries/cafés example, if we want to transport a certain amount of pains au chocolat/chocolatines then X will be a finite set containing the properties (precise location, size, etc) of both bakeries and cafés. The symmetry just means that the pains au chocolat/chocolatines are all the same and we do not want to distinguish which one is sent where. The main difference with the standard approach is that now we can define an OT problem without knowing exactly either the quantity of pains au chocolat/chocolatines produced or the demand. We discuss several mathematical properties of the GT-OT problem (25), within the framework of optimal transport theory. (1) We formulate the problem and show the existence of a minimizer $\gamma = (\gamma_n)_{n \ge 0}$, under appropriate assumptions on the costs. (2) We derive some properties on the *support* of $\gamma$, that is, on how many of the $\gamma_n$'s are non zero (we particularly discuss the Coulomb case). (3) We also study the truncated problem where all the $\gamma_n$ are assumed to vanish after some $N_{\max}$ and the convergence to the true problem when $N_{\max} \to \infty$. This is useful for numerical purposes. (4) We focus on duality theory and the existence of the dual potential. (5) We study the 1D case, confirming thereby some conjectures on the shape of the optimal plan made in [MSG13]. (6) Finally, we study the entropic regularization of the GC-OT problem.

In the last part of the chapter 2.3 we focus on the Lieb functional (23) which can be considered in some sense a quantum version of a standard MOT problem. Starting from these observations we aim at using some recent results about the moment approximation of MOT. More precisely, the
approach first considered in [Alf+21] consists in introducing an approximation, as the number of moments $M$ tends to infinity, of the exact multi-marginal transport problems where the marginal constraints are replaced by moment constraints associated to a $M$ "moment functions" which are real-valued functions defined on $\mathbb{R}^3$. The solution of this moment-constrained optimal transport problem is still a probability measure defined on $\mathbb{R}^{3m}$ but is also a sparse object in the sense that it can be written as a discrete measure charging a number of points belonging to $\mathbb{R}^{3m}$ which scales at most linearly with the number of moment constraints. Back to our setting, this translates into proving that the solutions of moment constrained approximations of (23) can be written under the form $\Gamma = \sum_{k=1}^K \alpha_k |\Psi_k\rangle\langle\Psi_k|$, where $K \in \mathbb{N}^*$ scales at most linearly with the number of moment constraints and $\Psi_k$ are some wave-functions. In particular we are able to prove (1) existence of minimizers for the approximate problem as well as the sparse structure. (2) We show convergence of the moment-constrained approximation towards the exact Lieb functional. Moreover, we also prove some rates of convergence of the associated approximation of the ground state energy to the exact one. (3) We finally present some results about the dual formulation of the moment-constrained problem noticing in particular that the dual problem can be re-written as a semi-definite positive programming problem.

## Unequal dimensional Optimal Transport

The final chapter of this manuscript deals with the case in which the number of properties of bakeries and cafés differs, that is $d_1 \neq d_2$ which yields to another generalization of optimal transport, well studied in [CMP17; MP20], known as **unequal dimensional optimal transport**. We especially focus on studying some variational problems involving this kind of optimal transport term, that is

$$ \inf_{\mu_1, \mu_2} \mathcal{J}(\mu_1, \mu_2) := \mathrm{OT}^c(\mu_1, \mu_2) + \mathcal{F}(\mu_2) + \mathcal{G}(\mu_1). \qquad (26) $$

We are interested in characterizing the minimizers of $(\mu_1, \mu_2) \mapsto \mathcal{J}(\mu_1, \mu_2)$, as well as the minimizers of the subproblems obtained when either $\mu_1$ or $\mu_2$ is fixed: $\mu_2 \mapsto \mathcal{J}(\mu_1, \mu_2)$ and $\mu_1 \mapsto \mathcal{J}(\mu_1, \mu_2)$. Problems of these general forms, for various choices of the functionals $\mathcal{F}$ and $\mathcal{G}$, arise in a wide variety of applications, including: gradient flows on Wasserstein space (where the minimization $\mu_2 \mapsto \mathcal{J}(\mu_1, \mu_2)$ represents one step in a discrete in time gradient flow, the celebrated JKO scheme [JKO98]), displacement interpolation (when $\mathcal{F}$ is the Wasserstein distance to a second probability measure), Cournot-Nash equilibria in game theory, city planning problems, etc. To clarify the reason for considering an unequal dimensional OT term, consider the subproblem where the population $\mu_1$, say the distribution of bakeries, is fixed and one, say a very rich man/woman with the uncontrollable desire of opening many cafés along the busiest metro line in Paris, wants to find the distribution of cafés which minimizes $\mu_2 \mapsto \mathrm{OT}^c(\mu_1, \mu_2) + \mathcal{F}(\mu_2)$: the cost of transporting the pains au chocolat/chocolatines plus a functional which, for example, avoids an excessive concentration of cafés at the same place. Now modeling the fact that $\mu_2$ is concentrated on a curve (the metro line) amounts to saying that the space $X_2$ is lower dimensional than $X_1$, for example $d_2 = 1$.

The first section of the chapter 3.1 is exactly devoted to studying the minimizers of such variational problems. Notice that when a joint condition, the **nestedness**, on the cost $c$ and marginals $\mu_1$ and $\mu_2$ is satisfied, the unequal dimensional OT problem is tractable (namely we can build an almost explicit solution). Unfortunately, in the present context, only the cost and **one** of the marginals (**neither** in the case of double minimizations) is prescribed. Here we prove that, under various conditions on $c$, $\mu_1$ and $X_2$, $(c, \mu_1, \mu_2)$ is nested whenever $\mu_2$ minimizes $\mu_2 \mapsto \mathcal{J}(\mu_1, \mu_2)$,
for a variety of different choices of the functional $\mathcal{F}$; analogous results for certain specific forms of $\mathcal{G}$ are also established for minimizations on the higher dimensional space, $\mu_1 \mapsto \mathcal{J}(\mu_1, \mu_2)$, and for double minimizations $(\mu_1, \mu_2) \mapsto \mathcal{J}(\mu_1, \mu_2)$. We go on to demonstrate that this a priori guarantee of nestedness makes the problem of characterizing or identifying the minimizers much more tractable; in different contexts, depending on the precise form of $\mathcal{F}$, we establish that solutions can be characterized by (local) differential equations, can be computed numerically by a convergent iterative scheme, or can be derived in almost closed form.

The second part of this chapter 3.2 deals with a class of metrics on the space of probability measures involving singular measures. More precisely, we develop the theory of a metric, which we call the $\nu$-based Wasserstein metric and denote by $W_\nu$, on the set of probability measures $\mathcal{P}(X)$ on a domain $X \subseteq \mathbb{R}^d$. This metric is based on a slight refinement of the well-known notion of *generalized geodesics*, see [AGS04], with respect to a base measure $\nu$ and is relevant in particular for the case when $\nu$ is singular with respect to $d$-dimensional Lebesgue measure; it is also closely related to the concept of linearized optimal transport [Wan+13]. The $\nu$-based Wasserstein metric is defined in terms of an iterated variational problem involving optimal transport to $\nu$; we also characterize it in terms of integration of classical Wasserstein metric ($\text{OT}^c$ with the quadratic cost) between the conditional probabilities when measures are disintegrated with respect to optimal transport to $\nu$, and through limits of certain multi-marginal optimal transport problems. Going back to the matching problem between bakeries and cafés, the idea underlying $W_\nu$ consists in matching the two populations by taking into account the information given by the base measure $\nu$. This is particularly clear when considering the characterization in terms of integration of classical Wasserstein metric between the conditional probabilities: consider $\nu$ as a one dimensional probability measure such that, for example, for a neighborhood/arrondissement $y$ of Paris (it is one dimensional because it could be just a set of labels), $\nu(y)$ tells you how important is matching bakeries and cafés near to it. Then the matching via the $W_\nu$ metrics works then as follows: (1) identify the bakeries and cafés close to some $y$ (the disintegration with respect to optimal transport to $\nu$), (2) match them using the classic Wasserstein distance.

As we vary the base measure $\nu$, the $\nu$-based Wasserstein metric interpolates between the usual quadratic Wasserstein metric (obtained when $\nu$ is a Dirac mass) and a metric associated with the uniquely defined generalized geodesics obtained when $\nu$ is sufficiently regular (eg, absolutely continuous with respect to Lebesgue). When $\nu$ concentrates on a lower dimensional submanifold of $\mathbb{R}^d$, we prove that the variational problem in the definition of the $\nu$-based Wasserstein metric has a unique solution. We also establish geodesic convexity of the usual class of functionals. We also introduce a class of metrics defined relative to a fixed higher dimensional base measure $\mu$, on the set of measures which are absolutely continuous with respect to a second fixed based measure $\sigma$. By means of this metric we prove convergence of an iterative scheme to solve a variational problem studied in the first section of the chapter.

Let me end this introduction by mentioning some works which are not summarized in this manuscript: [Ben+18] (study of second order variational mean-field games and equivalence with the minimization of a relative entropy), [DGN23] (analysis of the minimizers of the Levy-Lieb functional in the bosonic case) and [NP22] (numerical methods in order to solve variational problems involving unequal dimensional optimal transport).

**List of papers:**

* J.-D. Benamou, G. Carlier, M. Cuturi, L. Nenna, and G. Peyré. “Iterative Bregman
Projections for Regularized Transportation Problems". In: *SIAM Journal on Scientific Computing* 37.2 (2015), A1111–A1138. ISSN: 1064-8275. DOI: 10.1137/141000439

* J.-D. Benamou, G. Carlier, and L. Nenna. "A Numerical Method to Solve Multi-Marginal Optimal Transport Problems with Coulomb Cost". In: *Splitting Methods in Communication, Imaging, Science, and Engineering*. Springer International Publishing, 2016, pp. 577–601

* S. Di Marino, A. Gerolin, and L. Nenna. "Optimal Transportation Theory with Repulsive Costs". In: *Topological Optimization and Optimal Transport*. Radon Series on Computational and Applied Mathematics. https://www.degruyter.com/view/books/9783110430417/9783110430417_010/9783110430417-010.xml; De Gruyter, 2017

* M. Seidl, S. Di Marino, A. Gerolin, L. Nenna, K. J. H. Giesbertz, and P. Gori-Giorgi. "The Strictly-Correlated Electron Functional for Spherically Symmetric Systems Revisited". In: *Physical Review A: Atomic, Molecular, and Optical Physics* (2017)

* A. Blanchet, G. Carlier, and L. Nenna. "Computation of Cournot-Nash Equilibria by Entropic Regularization". In: *Vietnam Journal of Mathematics* (2017). ISSN: 2305-2228. DOI: 10.1007/s10013-017-0255-x

* J.-D. Benamou, G. Carlier, and L. Nenna. "Generalized Incompressible Flows, Multi-Marginal Transport and Sinkhorn Algorithm". In: *Numerische Mathematik* 142.1 (2019), pp. 33–54. ISSN: 0945-3245. DOI: 10.1007/s00211-018-0995-x

* J.-D. Benamou, G. Carlier, S. Di Marino, and L. Nenna. "An Entropy Minimization Approach to Second-Order Variational Mean-Field Games". 2018

* L. Nenna and B. Pass. "Variational Problems Involving Unequal Dimensional Optimal Transport". In: *Journal de Mathématiques Pures et Appliquées* 139 (2020), pp. 83–108

* S. Di Marino, M. Lewin, and L. Nenna. *Grand-Canonical Optimal Transport*. 2022. URL: http://arxiv.org/abs/2201.06859. preprint

* H. Ennaji, Q. Mérigot, L. Nenna, and B. Pass. "Robust Risk Management via Multi-Marginal Optimal Transport". In: *to appear on Journal of Optimization Theory and Applications* (2022). DOI: 10.48550/ARXIV.2211.07694

* L. Nenna and B. Pass. "A Note on Cournot-Nash Equilibria and Optimal Transport between Unequal Dimensions". In: *Optimal Transport Statistics for Economics and Related Topics*. Vol. 483. Studies in Systems, Decision and Control. arXiv, 2022. DOI: 10.48550/arXiv.2209.14888

* S. Di Marino, A. Gerolin, and L. Nenna. "Universal Diagonal Estimates for Minimizers of the Levy-Lieb Functional". In: *Letters in Mathematical Physics* (2023). DOI: 10.48550/arXiv.2303.00496

* L. Nenna and B. Pass. "Transport Type Metrics on the Space of Probability Measures Involving Singular Base Measures". In: *Applied Mathematics & Optimization* 87.2 (2023), p. 28
* L. Nenna and B. Pass. *An ODE Characterisation of Multi-Marginal Optimal Transport with Pairwise Cost Functions*. 2023. URL: http://arxiv.org/abs/2212.12492. preprint

* L. Nenna and P. Pegon. “Convergence Rate of Entropy-Regularized Multi-Marginal Optimal Transport Costs”. In: *to appear on Canadian Journal of Mathematics* arXiv:2307.03023 (2023). DOI: 10.48550/arXiv.2307.03023. preprint

* V. Ehrlacher and L. Nenna. *A Sparse Approximation of the Lieb Functional with Moment Constraints*. 2023. URL: http://arxiv.org/abs/2306.00806. preprint

List of papers in preparation:

* J.B. Casteras, L. Monsaingeon, **L. Nenna**, *Large Deviation Principle and Gamma-convergence for the Sticky-Schödinger problem*.

* V. Ehrlacher, **L. Nenna**, *A numerical method for moments constraint Lieb functional*.

* V. Ehrlacher, **L. Nenna**, *Reduced-order modeling for parametrized optimal transport problems*.

* L. De Pascale, **L. Nenna**, *Variational multi-populations mean field games*.

* Joshua Zoen-Git Hiew, **L. Nenna**, B. Pass, *An ordinary differential equation for entropic optimal transport and its linearly constrained variants*.

* **L. Nenna**, B. Pass, *Regularity results for some multi-marginal optimal transport problems*.

* M. Seidl, S. Di Marino, A. Gerolin, **L. Nenna**, K. Giesbertz, P. Gori-Giorgi, *The strictly-correlated electron functional for spherically symmetric systems revisited II: SGS CONJECTURE*.
# Chapter 1

# Multi-Marginal Optimal Transport and Entropic Regularization

This chapter summarizes the main results on Multi-Marginal Optimal Transport with a particular focus on the risk measure applications and the entropic counterpart.

All the results are contained in the following papers

Section 1.2 H. Ennaji, Q. Mérigot, L. Nenna, and B. Pass. “Robust Risk Management via Multi-Marginal Optimal Transport”. In: to appear on *Journal of Optimization Theory and Applications* (2022). DOI: 10.48550/ARXIV.2211.07694

Section 1.4 L. Nenna and P. Pegon. “Convergence Rate of Entropy-Regularized Multi-Marginal Optimal Transport Costs”. In: to appear on *Canadian Journal of Mathematics* arXiv:2307.03023 (2023). DOI: 10.48550/arXiv.2307.03023. preprint

Section 1.5 L. Nenna and B. Pass. *An ODE Characterisation of Multi-Marginal Optimal Transport with Pairwise Cost Functions*. 2023. URL: http://arxiv.org/abs/2212.12492. preprint

## 1.1 Introduction to Multi-Marginal Optimal Transport

In this section we briefly introduce a generalization of standard Optimal Transport problem to the case in which many marginals are involved instead of two.

Given $m$ compactly supported probability measures $\mu_i$ on sub-manifolds $X_i$ of dimension $d_i$ in $\mathbb{R}^d$ for $i \in \{1, \dots, m\}$ and a continuous cost function $c: X_1 \times X_2 \times \dots \times X_m \to \mathbb{R}_+$, the multi-marginal optimal transport problem consists in solving the following optimization problem

$$ \text{MOT}^c(\mu_1, \dots, \mu_m) := \inf_{\gamma \in \Pi(\mu_1, \dots, \mu_m)} \int_X c(x_1, \dots, x_m) d\gamma \quad (\text{MOT}) $$

where $X := X_1 \times X_2 \times \dots \times X_m$ and $\Pi(\mu_1, \dots, \mu_m)$ denotes the set of all probability measures $\gamma$ having $\mu_i$ as $i$-th marginal, i.e. $(e_i)_\gamma = \mu_i$ where $e_i: (x_1, \dots, x_m) \mapsto x_i$, for every $i \in \{1, \dots, m\}$. The formulation above is also known as the *Kantorovich problem* and it amounts to a linear minimization problem over a convex, weakly compact set; it is then not difficult to prove the existence of a solution by the direct method of calculus of variations. Much of the attention in the optimal transport community is rather focused on uniqueness and the structure of the minimizers. In particular one is mainly interested in determining if the solution is concentrated on the graph
of a function $(T_2, \dots, T_m)$ over the first marginal, where $(T_i)_{\#} \mu_1 = \mu_i$ for $i \in \{1, \dots, m\}$, in which case this function induces a solution à la Monge, that is $\gamma = (\mathrm{id}, T_2, \dots, T_m)_{\#} \mu_1$.

In the two marginals setting, the theory is fairly well understood and it is well-known that under mild conditions on the cost function (e.g. twist condition) and marginals (e.g. being absolutely continuous with respect to Lebesgue), the solution to (MOT) is unique and is concentrated on the graph of a function; we refer the reader to [San15] to have a glimpse of it. The extension to the multi-marginal case is still not well understood, but it has attracted recently a lot of attention due to a diverse variety of applications.

In particular in his seminal works [Pas11b; Pas12] Pass established some conditions, more restrictive than in the two marginals case, to ensure the existence of a solution concentrated on a graph. In this chapter we rely, for the high dimensional case, on the following (local) result in [Pas12] giving an upper bound on the dimension of the support of the solution to (MOT). Let $P$ be the set of partitions of $\{1, \dots, m\}$ into two non-empty disjoint subsets: $p = \{p_-, p_+\} \in P$ if $p_- \cup p_+ = \{1, \dots, m\}$, $p_- \cap p_+ = \emptyset$ and $p_-, p_+ \neq \emptyset$. Then for each $p \in P$ we denote by $g_p$ the bilinear form on the tangent bundle $TX$

$$g_p := D_{p-p_+}^2 c + D_{p_+p_-}^2 c \quad \text{where} \quad D_{pq}^2 c := \sum_{i \in p, j \in q} D_{x_i x_j}^2 c$$

for every $p, q \subseteq \{1, \dots, m\}$, and $D_{x_i x_j}^2 c := \sum_{\alpha_i, \alpha_j} \frac{\partial^2 c}{x_i^{\alpha_i} x_k^{\alpha_k}} dx_i^{\alpha_i} \otimes dx_j^{\alpha_j}$, defined for every $i, j$ on the whole tangent bundle $TX$. Define

$$G_c := \left\{ \sum_{p \in P} t_p g_p \mid (t_p)_{p \in P} \in \Delta_P \right\} \qquad (1.1)$$

to be the convex hull generated by the $g_p$, then it is easy to verify that each $g \in G_c$ is symmetric and therefore its signature, denoted by $(d^+(g), d^-(g), d^0(g))$, is well defined. Then, the following result from [Pas12] gives a control on the dimension of the support of the optimizer(s) in terms of these signatures.

**Theorem 1.1.1** (Part of [Pas12, Theorem 2.3]). Let $\gamma$ a solution to (MOT) and suppose that the signature of some $g \in G_c$ at a point $x \in X$ is $(d^+, d^-, d^0)$, that is the number of positive, negative and zero eigenvalues. Then, there exists a neighbourhood $N_x$ of $x$ such that $N_x \cap \text{spt } \gamma$ is contained in a Lipschitz sub-manifold of $X$ with dimension no greater than $\sum_{i=1}^m d_i - d^+$.

**Remark 1.1.2.** For the following it is important to notice that by standard linear algebra arguments we have for each $g \in G_c$ that $d^+(g) \le \sum_{i=1}^m d_i - \max_i d_i$. This implies that the smallest bound on the dimension of $\text{spt } \gamma$ which Theorem 1.1.1 can provide is $\max_i d_i$.

**Remark 1.1.3** (Two marginals case). When $m = 2$, the only $g \in G_c$ coincides precisely with the pseudo-metric introduced by Kim and McCann in [KM10]. Assuming for simplicity that $d_1 = d_2 = d$, they noted that $g$ has signature $(d, d, 0)$ whenever $c$ is non-degenerate so Theorem 1.1.1 generalizes their result since it applies even when non-degeneracy fails providing new information in the two marginals case: the signature of $g$ is $(r, r, 2d - 2r)$ where $r$ is the rank of $D_{x_1 x_2}^2 c$. Notice that this will help us in Section 1.4 to generalize the results established in [CPT23; EN22], concerning the rate of convergence of entropic optimal transport, to the case of a degenerate cost function.
It is well known that under some mild assumptions the Kantorovich problem (MOT) is dual to the following

$$ \sup \left\{ \sum_{i=1}^{m} \int_{X_i} \varphi_i(x_i) d\mu_i \mid \varphi_i \in \mathcal{C}_b(X_i), \sum_{i=1}^{m} \varphi_i(x_i) \le c(x_1, \dots, x_m) \right\}. \quad (\text{MD}) $$

Besides, it admits solutions $(\varphi_i)_{1 \le i \le m}$, called Kantorovich potentials, when $c$ is continuous and all the $X_i$'s are compact, and these solutions may be assumed $c$-conjugate, in the sense that for every $i \in \{1, \dots, m\}$

$$ \forall x \in X_i, \quad \varphi_i(x) = \inf_{(x_j)_{j \neq i} \in X_{-i}} c(x_1, \dots, x_{i-1}, x, x_{i+1}, \dots, x_m) - \sum_{1 \le j \le m, j \neq i} \varphi_j(x_j). \quad (1.2) $$

**One dimensional case** The most tractable of Multi-Marginal Optimal Transport problems occurs when the underlying spaces $X_i$ are one dimensional, and the mixed second derivatives of $c$ interact in a certain way.

**Definition 1.1.4 (Compatibility).** Suppose each $X_i \subset \mathbb{R}$ is a bounded real interval. Assume that $c \in \mathcal{C}^2(X_1 \times \dots \times X_m)$. We say that $s$ is strictly compatible if for each three distinct indices $i, j, k \in \{1, 1, \dots, m\}$ and each $(x_1, \dots, x_m) \in X_1 \times \dots \times X_m$ we have

$$ \frac{\partial^2 c}{\partial x_i \partial x_j} \left[ \frac{\partial^2 c}{\partial x_k \partial x_j} \right] \frac{\partial^2 c}{\partial x_k \partial x_i} (x_1, \dots, x_m) < 0. \quad (1.3) $$

We will say that $c$ is weakly compatible (or simply compatible) if for each $i \neq j$, we have either $\frac{\partial^2 c}{\partial x_i \partial x_j} \ge 0$ throughout $X_1 \times \dots \times X_m$ or $\frac{\partial^2 c}{\partial x_i \partial x_j} \le 0$ throughout $X_1 \times \dots \times X_m$, and

$$ \frac{\partial^2 c}{\partial x_i \partial x_j} \left[ \frac{\partial^2 c}{\partial x_k \partial x_j} \right] \frac{\partial^2 c}{\partial x_k \partial x_i} (x_1, \dots, x_m) \le 0, \quad (1.4) $$

for each distinct $i, j, k \in \{1, \dots, m\}$ and each $(x_1, \dots, x_m) \in X_1 \times \dots \times X_m$.

Note that if $c$ is strictly compatible, (1.3) implies that each $\frac{\partial^2 c}{\partial x_i \partial x_j} \neq 0$ throughout the connected domain $X_1 \times \dots \times X_m$; continuity then yields that each $\frac{\partial^2 c}{\partial x_i \partial x_j} \neq 0$ is either always positive or always negative. As (1.3) clearly implies (1.4), then strict compatibility implies compatibility.

The fact that the mixed partials $\frac{\partial^2 c}{\partial x_i \partial x_j}$ do not change signs under the compatibility condition allows us to partition the set $\{1, 2, \dots, m\} = P_+ \cup P_-$ of indices into disjoint subsets $P_+$ and $P_-$ such that $1 \in P_-$ and for each $i \neq j$, $\frac{\partial^2 c}{\partial x_i \partial x_j} \le 0$ throughout $X_1 \times \dots \times X_m$ if either both $i$ and $j$ are in $P_-$ or if both are in $P_+$, and $\frac{\partial^2 c}{\partial x_i \partial x_j} \ge 0$ throughout $X_1 \times \dots \times X_m$ otherwise. The inequalities are strict if strict compatibility holds.

**Definition 1.1.5 (Pseudo-inverse).** Given a probability measure $\mu$ on the real line $\mathbb{R}$, the generalized or pseudo-inverse of the cumulative function $F_\mu$ is defined as

$$ F_{\mu}^{-1}(t) = \inf\{x \in \mathbb{R} \mid F_{\mu}(x) \ge t\}. $$
**Definition 1.1.6.** For a compatible c, we define the c - comonotone coupling by:

$$
\gamma = (G_1, \dots, G_m) \# \mathcal{L}_{[0,1]}, \qquad (1.5)
$$

where $G_1 = F_{\mu_1}^{-1}$ and for each $i = 2, \dots, m$

$$
G_i(t) = \begin{cases} F_{\mu_i}^{-1}(t) & \text{if } i \in P_{-}, \\ F_{\mu_i}^{-1}(1-t) & \text{if } i \in P_{+}. \end{cases} \qquad (1.6)
$$

The following result can be found in [Car03] for strictly submodular functions (all the second mixed derivatives are strictly negative) as well as in [Pas11a] for compatible functions (where its proof appears together with the formulation of the compatibility conditions and the observation that it is equivalent to strict submodularity after changing coordinates).

**Theorem 1.1.7 ([Car03; Pas11a])**. *Suppose that c is strictly compatible. Then the c-comonotone coupling (1.5) is the unique optimizer in (MOT).*

Though the following result, we have obtained in [Enn+22], does not seem to be available in the literature, it is a straightforward consequence of Theorem 1.1.7 and the well known stability of optimal transport with respect to perturbations of the cost function. It asserts that the optimality of the measure constructed in Theorem 1.1.7 still holds if the strong compatibility assumption is relaxed to weak compatibility, although the uniqueness assertion may not.

**Proposition 1.1.8 ([Enn+22])**. *Suppose that c is weakly compatible. Then the c-comonotone coupling (1.5) is optimal in (MOT).*

1.2 From Risk measures to multi-marginal optimal transport

**Remark 1.2.1 (Notations).** Throughout all this section we adapt the language and the notations to the case in which optimal transport is applied to financial mathematics; that is instead of defining problem (MOT) as a minimization problem, we consider it as a maximization one and the cost function is now denoted as a surplus function b, namely

$$
\sup_{\gamma \in \Pi(\mu_1, \dots, \mu_m)} \int_{\mathbf{X}} b(x_1, \dots, x_m) d\gamma(x_1, \dots, x_m). \qquad (1.7)
$$

Notice that the two formulations are clearly equivalent by taking $b = -c$. All the theorems stated in the previous section still hold and will be used in the following.

In a variety of problems in operations research, a variable of interest $b = b(x_1, x_2, \dots, x_m)$ depends on several underlying random variables, whose individual distributions are known (or can be estimated) but whose joint distribution is not. A natural example arises when the variables represent parameters in a physical system whose individual distributions can be estimated empirically or through modeling (or a combination of both) but whose dependence structure cannot; an example of this flavour, originating in [IL15], in which the output $b$ represents the simulated height of a river at risk of flooding, and the underlying variables include various design parameters and climate dependent factors can be found below. Methods used in risk management to evaluate the resulting aggregate level of risk naturally depend on the distribution of the output variable $b$, and therefore, in turn, on the joint distribution of the $x_i$. A natural problem
is therefore to determine bounds, or worst-case scenarios, on these risk measures; that is, to maximize a given risk measure over all possible joint distributions of the $x_i$ with known marginal distributions.

Problems such as this have received extensive attention within the financial risk management community; see, for example, the monographs [Rüs13] and [MFE15] and the references within them. In the simplest of these problems, the underlying variables $x_i$ are typically real valued and the output variable $b$ is often assumed to have a particular structure (in many cases, it is a weighted sum of the $x_i$, reflecting the value of a portfolio built out of underlying assets with values $x_i$, or a function of this weighted sum). In these cases, explicit solutions for the maximizing couplings of the $x_i$ can sometimes be obtained (see, for example, [EWW15; EP06; PR13; Puc13; WW11]).

For more general output functions $b$, and possibly multi-dimensional underlying variables $x_i$, much less is known about the dependence structure of the maximizing joint distributions. Motivated by applications arising in the analysis of counter-party credit risk, a recent paper of Ghossoub-Hall-Saunders [GHS20] studies this more general setting systematically, allowing the underlying variables to take values in very general spaces and the output $b$ to take a very general form. They observe that the resulting problem is a generalization of the optimal transport problem; in fact, in the simplest case, when the spectral function is identically equal to one and the number $m$ of underlying variables is 2, the problem is exactly a classical optimal transport problem with surplus function $b = -c$. Although the analysis in [GHS20] focused on the $m = 2$ case, they note that their results can be extended to the $m \ge 3$ setting, in which case maximizing spectral risks becomes a generalization of the multi-marginal optimal transport problem.

In a joint work with H. Ennaji, Q. Mérigot and B. Pass we obtain the following results which will be more detailed below:

* We show that for any spectral risk measure, the maximization can in fact be formulated as a traditional multi-marginal optimal transport problem with $m+1$ marginals: the given marginals distributions of the $x_i$ as well as another distribution arising from the particular form of the spectral function (see Theorem 1.2.7 below). Although this formulation slightly increases the underlying dimension, it makes the maximization problem linear and much more tractable – indeed, the results and techniques, both theoretical and computational, in the substantial literature on multi-marginal optimal transport become applicable. Moreover, we derive an interesting equivalence between partial optimal transport and multi-marginal optimal transport which, up to our knowledge, has not been observed before in the literature, see Remark 1.2.4 and 1.2.8.

* We derive an explicit characterization of solutions for very general spectral risk measures and a substantial class of output functions $b$, through a careful refinement of the existing theory of multi-marginal optimal transport on one dimensional ambient spaces (see Theorem 1.2.10 and the examples following it). For $b$ falling outside this class, explicit solutions are likely generally unattainable; however, our formulation of the problem can potentially facilitate the use of a very broad range of computational methods for optimal transport to approximate solutions numerically (see Section 8.1 in [Ben21] and the references therein).

* For underlying variables in higher dimensional spaces, explicit solutions are generally not possible; however, there are known conditions on the cost/surplus function under which the solution to the multi-marginal optimal transport problem concentrates on a graph over the first variable, see [KP14a]. These conditions were not developed with our application in
mind; in particular, most previous work has focused on the case when each of the variables lives in set of the same dimension, whereas in our setting, the extra variable $x_0$ arising from the spectral function is one dimensional while the underlying variables $x_1, ..., x_m$ are higher dimensional. When $m = 2$, we provide a condition on $b$ under which the solution is of Monge type (Theorem 1.2.13).

* We also extend these ideas to the setting where the output function $b$ is multi-variate valued, using the maximal correlation risk measures developed in [Rüs06; EGH12]. We are again able to identify conditions under which the solution concentrates on a graph (Propositions 1.2.15 and 1.2.16); once more, these results do not seem readily attainable without using our optimal transport reformulation.

### 1.2.1 Equivalence between maximizing spectral risk measures and multi-marginal optimal transport

Let the real line describe a certain level of risk (e.g., the level of radiation in the nuclear power plant) and $\mu$ be a probability measure on $\mathbb{R}$ which can be interpreted as the distribution of risk. We will consider the following form of quantifier of the risk associated with $\mu$.

**Definition 1.2.2 (spectral risk measure).** A functional $\mathcal{R}: \mathcal{P}(\mathbb{R}) \to \mathbb{R}$ is called a spectral risk measure if it takes the form $\mathcal{R} = \mathcal{R}_\alpha$, where

$$ \mathcal{R}_{\alpha}(\mu) := \int_{0}^{1} F_{\mu}^{-1}(t)\alpha(t)dt. \qquad (1.8) $$

for a non-negative and nondecreasing function $\alpha: [0, 1] \to \mathbb{R}_+$ with $\int_{0}^{1} \alpha(t)dt = 1$.

**Remark 1.2.3.** Note that by taking $\alpha$ to be defined and real valued on $[0, 1]$, we are tacitly assuming that it is bounded, since $\alpha(0) \le \alpha(t) \le \alpha(1)$ by monotonicity; the same assumption was made in [GHS20]. This is mostly for technical convenience; we expect that most of our results can be extended to the more general case $\alpha: [0, 1) \to \mathbb{R}_+$, where, as above, $\alpha$ is monotone and $\int_{0}^{1} \alpha(t)dt = 1$, but possibly $\lim_{t \to 1} \alpha(t) = \infty$, under suitable hypothesis (for instance, decay conditions on $\alpha\#\mathcal{L}_{[0,1]}$).

For a given $\alpha$, we will often refer to the spectral risk measure $\mathcal{R}_\alpha$ as the $\alpha$-risk, and $\alpha$ as the spectral function.

It is well known in literature (see, for example, [Rüs06] and [EGH12], in which this fact is the basis for a multi-variate extension) the following variational characterization of spectral risk measures, that is given a probability measure $\mu$ with $\int_{\mathbb{R}} x d\mu(x) > -\infty$,

$$ \mathcal{R}_{\alpha}(\mu) = \max_{\pi \in \Pi(\alpha_{\#}\mathcal{L}_{[0,1]}, \mu)} \int_{\mathbb{R} \times \mathbb{R}} xy d\pi(x, y), \qquad (1.9) $$

where we have denoted $\Pi(\alpha_{\#}\mathcal{L}_{[0,1]}, \mu)$ the space of probability measures on $\mathbb{R}^2$ with marginals $\alpha_{\#}\mathcal{L}_{[0,1]}$ and $\mu$. Notice that thanks to this characterization the spectral risk measure is indeed a standard optimal transport problem with cost function $c(x, y) = -xy$. One can easily prove this by exploiting the fact that in one dimension the optimal transport plan is obtained by using the pseudo inverse of the cumulative function, as (1.5). Return now to the case in which there are $m$ parameters which enter into the estimation of risk through the output function $b$, and
we wish to evaluate the worst case scenario for the α-risk $\mathcal{R}_\alpha(b_\#\gamma)$ of $b_\#\gamma$ among couplings $\gamma \in \Pi(\mu_1, \mu_2, \dots, \mu_m)$ of the marginal distributions $\mu_i$. That is, we want to maximize:

$$ \max_{\gamma \in \Pi(\mu_1, \mu_2, \dots, \mu_m)} \mathcal{R}_\alpha(b_\#\gamma). \qquad (1.10) $$

**Remark 1.2.4** (Maximizing Expected Shortfall is optimal partial transport). A special case of particular importance in risk management applications occurs when $\alpha = \alpha_{m_0} := \frac{1}{m_0} \mathbf{1}_{[1-m_0,1]}$, in which case (1.8) is also known as the Expected Shortfall, or Conditional Value at Risk. In this setting, the maximization problem (1.10) is actually equivalent to a well known variant of the optimal transport problem; indeed, it can be reformulated into

$$ \max_{\gamma \in \Pi_{m_0}(\mu_1, \dots, \mu_m)} \frac{1}{m_0} \int b(x_1, \dots, x_m) d\gamma(x_1, \dots, x_m) $$

where $\Pi_{m_0}(\mu_1, \dots, \mu_m)$ denotes the set of non-negative measures $\gamma$ on $X_1 \times \dots \times X_m$ with total mass $m_0$, such that its $i$-marginal $\gamma_i$ is dominated by $\mu_i$ for each $i$; that is $\int_{X_i} \varphi d\gamma_i \le \int_{X_i} \varphi d\mu_i$ for all non-negative test functions $\varphi \in \mathcal{C}^0(X_i)$. This is known as the optimal partial transport problem when $m=2$ [CM10; Fig10], and the multi-marginal optimal partial transport problem [KP14c] when $m > 2$. As we will show below, the problem is in fact equivalent to an ordinary multi-marginal optimal transport problem with an additional marginal; see Theorem 1.2.7 and Remark 1.2.8 below.

Our first main contribution is to show that the spectral risk maximization problem (1.10) is equivalent to the multi-marginal optimal transport problem (1.7) with $X_0 = [\alpha(0), \alpha(1)] \subseteq \mathbb{R}$, $\mu_0 = \alpha_{\#} \mathcal{L}_{[0,1]}$, the other $X_i$ and $\mu_i$ representing the domains and distributions of the underlying variables, respectively, and

$$ s(x_0, x_1, \dots, x_m) = x_0 b(x_1, \dots, x_m). \qquad (1.11) $$

The intuition behind this is simple: if the variables $(x_0, \dots, x_m)$ are optimally coupled for $s$, $x_0$ and the variable $b$ must be coupled in an increasing way, this is captured precisely by Lemma 1.2.5. This monotonicity leads exactly to the spectral risk measure evaluated on the distribution of $b$.

**Lemma 1.2.5** ([Enn+22]). Suppose that $\pi \in \Pi(\mu_0, \mu_1, \dots, \mu_m)$ and let $\gamma$ be the $(1, \dots, m)$-marginal of $\pi$, that is the projection of $\pi$ on $X_1 \times \dots \times X_m$. Then, for $s$ given by (1.11),

$$ \int_{X_0 \times \mathbf{X}} s(x_0, x_1, \dots, x_m) d\pi(x_0, x_1, \dots, x_m) \le \mathcal{R}_\alpha(b_\#\gamma) $$

Furthermore, we have equality if and only if the support of

$$ \tau_\pi = \left( (x_0, x_1, x_2, \dots, x_d) \mapsto (x_0, b(x_1, x_2, \dots, x_d)) \right)_{\#} \pi \in \mathcal{P}(\mathbb{R}^2) $$

is monotone increasing.

**Lemma 1.2.6** ([Enn+22]). Given any measure $\gamma \in \Pi(\mu_1, \dots, \mu_m)$, there exists a $\pi \in \Pi(\mu_0, \mu_1, \dots, \mu_m)$ whose $(1, \dots, m)$-marginal is $\gamma$, such that $((x_0, x_1, \dots, x_m) \mapsto (x_0, b(x_1, x_2, \dots, x_m))_{\#} \pi$ has monotone increasing support.
The proof of Lemma 1.2.5 relies on the variational characterization of the spectral risk (1.9) where as for Lemma 1.2.6 one uses comonotonic transport plans and measure disintegration. We are now ready to establish the equivalence between the maximal $\alpha$-risk problem (1.10) and the multi-marginal optimal transport problem (MOT).

**Theorem 1.2.7** ([Enn+22]). A probability measure $\pi$ in $\Pi(\mu_0, \mu_1, \dots, \mu_m)$ is optimal for (1.7) with cost function (1.11) if and only if its $(1, \dots, d)$-marginal is optimal in (1.10), and

$$ \tau_{\pi} = \left( (x_0, x_1, x_2, \dots, x_d) \mapsto (x_0, b(x_1, x_2, \dots, x_d)) \right)_{\#} \pi $$

has monotone increasing support.

**Remark 1.2.8.** Returning to the Expected Shortfall Case, $\alpha = \alpha_{m_0} = \frac{1}{m_0}1_{[1-m_0,1]}$, in view of Remark 1.2.4, Theorem 1.2.7 shows that the optimal partial transport problem for the surplus $b(x_1, \dots, x_m)$ can be transformed into a multi-marginal transport problem for the surplus $s = x_0b$ and additional marginal $\mu_0 = \alpha_{\#}L_{[0,1]}$. To the best of our knowledge, this equivalence between these two well studied mathematical problems has not been observed before in the optimal transport literature. In addition, this perspective, together with results in [CM10; Fig10; KP14c] allows us to immediately identify conditions under which the active (that is, the part that couples to $\alpha > 0$) part of the optimal $\gamma$ in (1.10) is uniquely determined and concentrates on a graph over $x_1$. Furthermore, algorithms to compute the solution are readily available [Ben+15; IN18].

### 1.2.2 Solutions for one-dimensional assets and compatible payouts

We now turn our attention to the structure of maximizers in (1.10) when the underlying variables are one dimensional. Assume that each $x_i \in \mathbb{R}$ and each $\mu_i$ is supported on an interval, $X_i = [\underline{x}_i, \bar{x}_i]$, for all $i = 1, \dots, m$. We note that in our setting, the first marginal, $\mu_0$ is supported always on the interval $[\underline{x}_0, \bar{x}_0] := [\alpha(0), \alpha(1)]$ with $\alpha(0) \ge 0$.

**Lemma 1.2.9** ([Enn+22]). Suppose that $b$ is compatible and monotone increasing in each $x_i \in P_+$ and monotone decreasing for each $x_i \in P_-$. Then the $s$-comonotone coupling (1.5) is optimal for (1.7) with surplus function given by (1.11). The maximal value is

$$ \int_0^1 \alpha(t)b(G_1(t), G_2(t), \dots, G_m(t))dt. \qquad (1.12) $$

Furthermore, if in addition the monotonicity of $b$ with respect to each argument is strict, the solution is unique on the support of $\alpha$; that is, any other solution $\bar{\pi}$ coincides with $\pi = (G_\alpha, G_1, \dots, G_m)_{\#}L_{[0,1]}$ on $(0, \alpha(1)] \times [\underline{x}_1, \bar{x}_1] \times \dots \times [\underline{x}_m, \bar{x}_m]$.

Clearly full uniqueness of the optimal $\pi$ cannot hold if $\alpha = 0$ on a set of positive measure, or, equivalently, $\mu_0(\{0\}) > 0$; in this case, we can rearrange the part of the $(1, \dots, d)$-marginal $\gamma$ of $\pi$ that couples to $\alpha = 0$ in any way without affecting the value of $\int sd\pi$. The preceding Lemma identifies conditions under which this is the only source of non-uniqueness, so that the optimal $\pi$ is uniquely determined on $(0, \alpha(1)] \times [\underline{x}_1, \bar{x}_1] \times \dots \times [\underline{x}_m, \bar{x}_m]$.

Note also that we are able to obtain uniqueness off the set where $\alpha = 0$ despite the fact that the surplus function may not be strongly compatible. This is done by exploiting an idea developed in [PV22] (see in particular Proposition 4.1 there); namely that one only needs strongly compatible interactions between certain key pairs of variables (rather than all of them).

Thanks to Theorem 1.2.7, the preceding result then easily yields the following characterization of solutions to (1.10).
**Theorem 1.2.10 ([Enn+22])**. Assume that b is compatible, monotone increasing in each $x_i \in P_+$ and monotone decreasing in each $x_i \in P_-$. Then $(G_1, \dots, G_m)_{\#\mathcal{L}_{[0,1]}}$ maximizes (1.10) and the maximal value is given by (1.12). If in addition the monotonicity is strict, then letting $a = \mu_0\{0\}$ any other maximizer $\gamma$ must couple the regions defined by $G_i([0, a])$, and must coincide with the b-comonotone coupling $(G_1, \dots, G_m)_{\#\mathcal{L}_{[0,1]}}$ on $G_1([a, 1]) \times G_2([a, 1]) \times \dots \times G_m([a, 1])$.

**Example 1.2.11.** If $b(x_1, \dots, x_m) = -\sum_{1 \le i < j \le m} (x_i - x_j)^2$ and if $\alpha = 1$, problem (1.10) is equivalent to the computation of Wasserstein barycenters [AC11], while with $\alpha = \alpha_{m_0}$ we get the partial Wasserstein barycenter problem [KP14c]. In both cases, since b is compatible the answer can be calculated explicitly by Theorem 1.2.10.

**Example 1.2.12 (Sensitivity analysis and maximal river flow).** We briefly describe a recurring example from [IL15], used throughout that paper to illustrate issues in sensitivity analysis. In that setting, one wants to understand the influence of the dependence structure (among other factors) between several contributing inputs on an output behavior. We consider a simple model which involves the height of river at risk of flooding and compares it to the height of a dyke which protects industrial facilities. The maximal annual overflow $S$ of a river is modelled by

$$S = Z_{\nu} + \left( \frac{Q}{BK_{s}\sqrt{\frac{Z_{m}-Z_{\nu}}{L}}} \right)^{0.6} - H_{d} - C_{b},$$

where $\left( \frac{Q}{BK_s\sqrt{\frac{Z_m-Z_\nu}{L}}} \right)^{0.6}$ is the maximal annual height of the river and the variables $Q, K_s, Z_\nu, Z_m, H_d, C_b, L, B$ are physical quantities whose values are modelled as random variables due to their variation in time and space, measurement inaccuracies, or uncertainty of their true values. Their individual distributions are modelled in [IL15] (see Table 1 on p.4), and are all absolutely continuous with respect to Lebesgue measure. The $\alpha$-risk (1.8), where the $x_i$ are the variables $Q, K_s, Z_\nu, Z_m, H_d, C_b, L, B$ and $b = S$ is the overflow, quantifies the risk of overflow. In [IL15], the variables were assumed to be independent, although other dependence structures are certainly possible; (1.10) asks what is the maximal risk over all possible dependence structures. Notice then that the surplus function $b$ is compatible and satisfies the strict monotonicity with respect to each variable required in Theorem 1.2.10. The $(1, \dots, m)$ marginal $\gamma$ of the s-comonotone solution $\pi$ defined in (1.5) is optimal, and the unique optimizer on the support of $\alpha$, $\{\alpha > 0\}$. In this case, the $G_i$ corresponding to $Z_\nu, Q$ and $L$ are monotone increasing, while the $G_i$ corresponding to the other variables are decreasing. They can be retrieved explicitly by computing the quantile functions associated to each marginals as described above.

### 1.2.3 Higher dimensional assets

As explained at the beginning of the chapter, there is an analogous theory of multi-marginal optimal transport when the underlying variables lie in more general spaces. Although in these cases it is generally not possible to derive explicit solutions as we did above, it is possible to prove analogous structural properties of optimal couplings, namely, that solutions are of Monge type (that is, concentrated on graphs over $x_1$), for certain surplus functions $s$ (see [G$\check{S}$98] for an early result in this direction, for a particular $s$, and [KP14a] and [Pas11b] for general, sufficient conditions on $s$).

We present one such result below, when each $X_i \subseteq \mathbb{R}^d$, to illustrate how the theory can be adapted to the present setting. For simplicity, we restrict our attention to the $m = 2$ case.
The conditions we impose on $b$ are related to the local differential condition found in [Pas11b]; however, the proof in [Pas11b] was restricted to the case when each marginal was supported on a space of the same dimension, whereas here our first marginal $\mu_0$, corresponding to $x_0 = \alpha$, has one dimensional support while the other marginals are supported in $\mathbb{R}^d$. Our proof is inspired by the argument in [Pas11b], but modified to fit the present setting here. Similar results can be established for larger $m$, using similar arguments. However, the conditions imposed on $b$ become more stringent for larger $m$, and more complicated to state. In what follows, we will assume that $b$ is $\mathbb{C}^2$; $D_{x_i}b(x_1, x_2) = (\frac{\partial b}{\partial x_i^1}, \frac{\partial b}{\partial x_i^2}, \dots, \frac{\partial b}{\partial x_i^d})$ represents the gradient of the function $b: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$ with respect to the variable $x_i \in \mathbb{R}^d$, for $i = 1, 2$. Similarly, $D_{x_1x_2}^2 b(x_1, x_2) = (\frac{\partial^2 b}{\partial x_1^1 \partial x_2^1})_{rl}$ is the $d$ by $d$ matrix of mixed second order derivatives (each entry is a derivative with respect to one coordinate from each of $x_1$ and $x_2$).

**Theorem 1.2.13 (Non-linear Brenier Theorem, [Enn+22]).** Suppose that $m = 2$, that the domains $X_1, X_2 \subseteq \mathbb{R}^d$ are compact and that $\mu_1$ is absolutely continuous with respect to Lebesgue measure. Assume that $x_2 \mapsto D_{x_1}b(x_1, x_2)$ is injective for each fixed $x_1 \in X_1$, and that for each $(x_1, x_2) \in X_1 \times X_2$, $\det(D_{x_1x_2}^2b(x_1, x_2)) \neq 0$ and

$$ \langle D_{x_2}b(x_1, x_2), [D_{x_1x_2}^2 b(x_1, x_2)]^{-1} D_{x_1}b(x_1, x_2) \rangle > 0. \quad (1.13) $$

Then the part of the solution to (1.10) away from $\alpha = 0$ concentrates on the graph of a function over $x_1$. Furthermore, if $|\{\alpha = 0\}| = 0$, the solution is unique.

**Remark 1.2.14.** The injectivity of $x_2 \mapsto D_{x_1}b(x_1, x_2)$ is known as the twist condition in the optimal transport literature; it is well known that, together with the absolute continuity of $\mu_1$, it guarantees the Monge structure of the solution to the two marginal optimal transport problem with surplus $b$ (see, for example, [San15]). The invertibility of $D_{x_1x_2}^2b(x_1, x_2)$ is frequently referred to as non-degeneracy, and can be seen as a linearized version of the twist. Here we require additional hypotheses as we are dealing with the more sophisticated 3 marginal problem with surplus $s(\alpha, x_1, x_2) = \alpha b(x_1, x_2)$, which is equivalent to the spectral risk maximization (1.10) by Theorem 1.2.7.

### 1.2.4 Multidimensional measures of risk

We end this section by showing a final result obtained in [Enn+22] concerning the framework proposed in [Rüs06] in which risk is measured in a multi-dimensional way: instead of a single output variable, we now have several, depending on the same underlying random variables; this corresponds to a vector valued output function $b: X_1 \times \dots \times X_m \to \mathbb{R}^d$. A natural form of multi-variate risk measures on the distribution $b_{\#}\gamma$ of output variables is then the maximal correlation measure from [Rüs06], which is defined by

$$ \mathcal{R}_{\nu}(b_{\#}\gamma) = \max_{\eta \in \Pi(b_{\#}\gamma, \nu)} \int_{\mathbb{R}^d \times \mathbb{R}^d} \langle z, y \rangle \, d\eta $$

for some probability measure $\nu \in \mathcal{P}(\mathbb{R}^d)$. We note that it was proven in [EGH12] that any strongly coherent multi-variate risk measure takes this form. We consider the problem of maximizing $\mathcal{R}_{\nu}(b_{\#}\gamma)$ over all $\gamma \in \Pi(\mu_1, \dots, \mu_m)$, where the $\mu_i$ as before represent the distributions of the underlying variables. Exactly as in Theorem 1.2.7, one can show that this problem is equivalent to the multi-marginal problem

$$ \max_{\pi \in \Pi(\nu, \mu_1, \dots, \mu_d)} \int \langle b(x_1, \dots, x_m), y \rangle \, d\pi. \quad (1.14) $$
This problem is more challenging than the case of a scalar valued b; nonetheless, we are able to obtain the following results in particular cases.

**Proposition 1.2.15** ([Enn+22]). Suppose that the underlying variables $x_i$ are one dimensional. Assume that $\nu$ is concentrated on a smooth curve, that is $\nu = f_{\#L_{[0,1]}}$ with $f : [0, 1] \to \mathbb{R}^d$, where each component of $f$ is positive and monotone increasing, and each component $b_j$ of $b = (b_1, b_2, \dots, b_n)$ has positive second mixed derivatives and is monotone increasing in each $x_i$. Then $((t, x_1, \dots, x_m) \mapsto (f(t), x_1, \dots, x_m))_{\#\pi}$ is optimal in (1.14), where $\pi = (\text{id}, F_{\mu_1}^{-1}, F_{\mu_2}^{-1}, \dots, F_{\mu_m}^{-1})_{\#L_{[0,1]}}$ is the comonotone coupling of $\mu_0, \mu_1, \dots, \mu_m$, where $\mu_0 = L_{[0,1]}$.

Notice that the positivity of mixed second derivatives of each components $b_j$ implies that they are compatibles so that one can again exploit the fact that in one dimension the optimal plan can be obtained by using the pseudo-inverse of the cumulative functions. In particular it turns out that the cost function $\langle b(x_1, \dots, x_m), f(t) \rangle = \sum_{j=1}^d b_j(x_1, \dots, x_m)f_j(t)$ is compatible, and so Proposition 1.1.8 implies the desired result.

For more general, diffuse $\nu$, we are able to prove that the solution is of Monge form and unique, provided that $b$ is invertible.

**Proposition 1.2.16** ([Enn+22]). Assume that $\nu$ is absolutely continuous with respect to $d$ dimensional Lebesgue measure and $b : X_1 \times \dots \times X_m \to \mathbb{R}^d$ is invertible. Then there exists a unique solution to (1.14). Furthermore, it concentrates on a graph over $y$.

Note that the invertibility assumption on $b$ implies that the sum of the dimensions of the $X_i$ must be less than or equal to the dimension $d$ of $\nu$.

## 1.3 Entropic Multi-Marginal Optimal Transport

We consider a $m$-uple of probability measures $\mu_i$ compactly supported on sub-manifolds $X_i \subseteq \mathbb{R}^d$ of dimension $d_i$ and a cost function $c : X_1 \times \dots \times X_m \to \mathbb{R}_+$. The Entropic Multi-Marginal Optimal Transport problem is defined as :

$$ \text{MOT}_{\varepsilon}^{c} := \inf \left\{ \int_{X_1 \times \dots \times X_m} c d\gamma + \varepsilon \text{Ent}(\gamma | \otimes_{i=1}^{m} \mu_i) \mid \gamma \in \Pi(\mu_1, \dots, \mu_m) \right\}, \quad (\text{MOT}_{\varepsilon}) $$

where $\text{Ent}(\cdot | \otimes_{i=1}^{m} \mu_i)$ is the Boltzmann-Shannon relative entropy (or Kullback-Leibler divergence) w.r.t. the product measure $\otimes_{i=1}^{m} \mu_i$, defined for general probability measures $p, q$ as

$$ \text{Ent}(p|q) = \begin{cases} \int_{\mathbb{R}^d} \rho \log(\rho) dq & \text{if } p = \rho q, \\ +\infty & \text{otherwise.} \end{cases} $$

Notice that the multi-marginal optimal transport problem we have seen before corresponds to the case where $\varepsilon = 0$.

Let us highlight that Entropic optimal transport (EOT) has found applications and proved to be an efficient way to approximate Optimal Transport (OT) problems, especially from a computational viewpoint. Indeed, when it comes to solving EOT by alternating Kullback-Leibler projections on the two marginal constraints, by the algebraic properties of the entropy such iterative projections correspond to the celebrated Sinkhorn's algorithm [Sin64], applied in this framework in the pioneering works [Cut13; Ben+15]. The simplicity and the good convergence
guarantees (see [FL89; MG20]) of this method compared to the algorithms used for the OT problems, then determined the success of EOT for applications in machine learning, statistics, image processing, language processing and other areas (see the monograph [PC19] and references therein). As for multi-marginal optimal transport one can derive the dual problem of (MOT<sub>ε</sub>) which reads as

$$ \mathrm{MOT}_{\varepsilon}^{c} = \varepsilon + \sup \left\{ \sum_{i=1}^{m} \int_{X_i} \varphi_i(x_i) \mathrm{d}\mu_i - \varepsilon \int_{\mathbf{X}} e^{\frac{\sum_{i=1}^{m} \varphi_i(x_i) - c(\mathbf{x})}{\varepsilon}} \mathrm{d}\otimes_{i=1}^{m} \mu_i \mid \varphi_i \in \mathcal{C}_b(X_i) \right\}, \quad (\mathrm{MD}_{\varepsilon}) $$

which is invariant by $(\varphi_1, \dots, \varphi_m) \mapsto (\varphi_1 + \lambda_1, \dots, \varphi_m + \lambda_m)$ where $(\lambda_1, \dots, \lambda_m) \in \mathbb{R}^m$ and $\sum_{i=1}^m \lambda_i = 0$, see [Léo14; NW22; MG20] for some recent presentations. It admits an equivalent log-sum-exp form:

$$ \mathrm{MOT}_{\varepsilon}^{c} = \sup \left\{ \sum_{i=1}^{m} \int_{X_i} \varphi_i(x_i) \mathrm{d}\mu_i - \varepsilon \log \left( \int_{\mathbf{X}} e^{\frac{\sum_{i=1}^{m} \varphi_i(x_i) - c(\mathbf{x})}{\varepsilon}} \mathrm{d}\otimes_{i=1}^{m} \mu_i \right) \mid \varphi_i \in \mathcal{C}_b(X_i) \right\}, \quad (\mathrm{MD}_{\varepsilon}') $$

which is invariant by the same transformations without assuming $\sum_{i=1}^m \lambda_i = 0$.

From $(\mathrm{MOT}_\varepsilon)$ and $(\mathrm{MD}_\varepsilon)$ we recover, as $\varepsilon \to 0$, the unregularized multi-marginal optimal transport (MOT) and its dual (MD) we have introduced above. The link between multi-marginal optimal transport and its entropic regularization is very strong and a consequence of the $\Gamma$-convergence of $(\mathrm{MOT}_\varepsilon)$ towards (MOT) (one can adapt the proof in [Car+17] or see [BCN19; GKR20] for $\Gamma$-convergence in some specific cases) is that

$$ \lim_{\varepsilon \to 0} \mathrm{MOT}_{\varepsilon}^{c} = \mathrm{MOT}_{0}^{c}. $$

By the direct method in the calculus of variations and strict convexity of the entropy, one can show that $(\mathrm{MOT}_\varepsilon)$ admits a unique solution $\gamma_\varepsilon$, called *optimal entropic plan*. Moreover, there exist $m$ real-valued Borel functions $\varphi_i^\varepsilon$ such that

$$ \gamma_{\varepsilon} = \exp \left( \frac{\bigoplus_{i=1}^{m} \varphi_{i}^{\varepsilon} - c}{\varepsilon} \right) \otimes_{i=1}^{m} \mu_{i}, \quad (1.15) $$

where $\bigoplus_{i=1}^m \varphi_i^\varepsilon := (x_1, \dots, x_m) \mapsto \sum_{i=1}^m \varphi_i^\varepsilon(x_i)$, and in particular we have that

$$ \mathrm{MOT}_{\varepsilon}^{c} = \sum_{i=1}^{m} \int_{X_i} \varphi_i^{\varepsilon} \mathrm{d}\mu_i \quad (1.16) $$

and these functions have continuous representatives and are uniquely determined up a.e. to additive constants. The reader is referred to the analysis of [MG20], to [Nen16] for the extension to the multi-marginal setting, and to [BL92; BLN94; Csi75; FG97; RT98] for earlier references on the two marginals framework.

The functions $\varphi_i^\varepsilon$ in (1.15) are called *Schrödinger potentials*, the terminology being motivated by the fact that they solve the dual problem $(\mathrm{MD}_\varepsilon)$ and are as such the (unique up to invariance) solutions to the so-called *Schrödinger system*: for all $i \in \{1, \dots, m\}$,

$$ \varphi_i(x_i) = -\varepsilon \log \int_{\mathbf{X}_{-i}} e^{\frac{\bigoplus_{1 \le j \le m, j \ne i} \varphi_j^\varepsilon - c(\mathbf{x})}{\varepsilon}} \mathrm{d}\bigoplus_{1 \le j \le m, j \ne i} \mu_j \quad \text{for } \mu_i\text{-a.e. } x_i, \quad (1.17) $$

where $\mathbf{X}_{-i} = \prod_{1 \le j \le m, j \ne i}^m X_j$. Note that (1.17) is a *softmin* version of the multi-marginal c-conjugacy relation for Kantorovich potentials.
## 1.4 Rate of convergences for entropic multi-marginal optimal transport

Since $MOT_{\varepsilon}^c$ can be seen as a perturbation of $MOT_0^c$, it is natural to study the behaviour as $\varepsilon$ vanishes. In this section we summarize the joint work with P. Pegon [NP23c] where we were mainly interested in investigating the rate of convergence of the entropic cost $MOT_{\varepsilon}^c$ to $MOT_0^c$ under some mild assumptions on the cost functions and marginals. In particular we have extended the techniques introduced in [CPT23] to the multi-marginal case which will also let us generalize the bounds in [CPT23] to the case of degenerate cost functions. For the two marginals and non-degenerate case we also refer the reader to a very recent (and elegant) paper [MS23] where the authors push a little further the analysis of the convergence rate by disentangling the roles of $\int cd\gamma$ and the relative entropy in the total cost and deriving convergence rate for both these terms. Notice that concerning the convergence rate of the entropic multi-marginal optimal transport an upper bound has been already established in [EN22], which depends on the number of marginals and the quantization dimension of the optimal solutions to $(MOT_{\varepsilon})$ with $\varepsilon = 0$.

Our main findings can be summarized as follows:

* we establish two **upper bounds**, one valid for locally Lipschitz costs and a finer one valid for locally semi-concave costs. For the latter, we improve the upper bound by a 1/2 factor, obtaining the following inequality for some $C^* \in \mathbb{R}_+$

$$MOT_{\varepsilon}^{c} \leq MOT_{0}^{c} + \frac{1}{2} \left( \sum_{i=1}^{m} d_i - \max_{1 \leq i \leq m} d_i \right) \varepsilon \log(1/\varepsilon) + C^{*} \varepsilon. \quad (1.18)$$

We stress that this upper bound is smaller or equal than the one provided in [EN22, Theorem 3.8], which is of the form $\frac{1}{2}(m-1)D\varepsilon \log(1/\varepsilon) + O(\varepsilon)$ where $D$ is a quantization dimension of the support of an optimal transport plan. Thus $D$ must be greater or equal than the maximum dimension of the support of the marginals, and of course $\sum_{i=1}^{m} d_i - \max_{1 \leq i \leq m} d_i \leq (m-1)\max_{1 \leq i \leq m} d_i$. The inequality may be strict for example in the two marginals case with unequal dimension.

* For the **lower bound**, from the dual formulation of $(MOT_{\varepsilon})$ we have

$$MOT_{\varepsilon}^{c} \geq MOT_{0}^{c} - \varepsilon \log \int_{X} e^{-\frac{E(x_{1}, \dots, x_{m})}{\varepsilon}} d \otimes_{i=1}^{m} \mu_{i}(x_{i}),$$

where $E(x_1, \dots, x_m) = c(x_1, \dots, x_m) - \bigoplus_{i=1}^{m} \varphi_i(x_i)$ is the duality gap and $(\varphi_1, \dots, \varphi_m)$ are Kantorovich potentials for the un-regularized problem $(MOT_{\varepsilon})$ with $\varepsilon = 0$. By proving that $E$ detaches quadratically from the set $\{E = 0\}$ allows us to estimate the previous integral so that we get the following: given a $\kappa$ depending on a signature condition (see $(\text{PS}(\kappa))$) on the second mixed derivatives of the cost, the lower bound can be summarized as follows

$$MOT_{\varepsilon}^{c} \geq MOT_{0}^{c} + \frac{\kappa}{2} \varepsilon \log(1/\varepsilon) - C_{*} \varepsilon. \quad (1.19)$$

for some $C_* \in \mathbb{R}_+$.
1.4.1 Upper bound

In this section we detail a little more our main results concerning the upper bounds and we give some elements of proof. We first consider the case of locally Lipschitz costs for which we obtain the following theorem.

**Theorem 1.4.1 ([NP23c])**. Assume that for $i \in \{1, \dots, m\}$, $\mu_i \in \mathcal{P}(X_i)$ is a compactly supported measure on a Lipschitz sub-manifold $X_i$ of dimension $d_i$ and $c \in \mathbb{C}_{\text{loc}}^{0,1}(\mathbf{X})$, then

$$ \mathrm{MOT}_{\varepsilon}^{c} \leq \mathrm{MOT}_{0}^{c} + \left( \sum_{i=1}^{m} d_i - \max_{j \in \{1, \dots, m\}} d_j \right) \varepsilon \log(1/\varepsilon) + O(\varepsilon). \quad (1.20) $$

The proofs rely, as in [CPT23], on a multi-marginal variant of the block approximation introduced in [Car+17]. In particular, given the optimal solution $\gamma_0$ to the un-regularized problem one can build a competitor $\gamma_\delta$ for the entropic problem by using the block approximation. Then we have

$$ \mathrm{MOT}_{\varepsilon}^{c} - \mathrm{MOT}_{0}^{c} \leq \int cd(\gamma_{\delta} - \gamma_{0}) + \varepsilon \mathrm{Ent}(\gamma_{\delta} | \otimes_{i=1}^{m} \mu_{i}), \quad (1.21) $$

where the first term can be controlled by using the fact that the cost function is Lipschitz and the second term is smaller equal than $\sum_{i=1}^{m-1} H_\delta(\mu_i)$ where $H_\delta(\mu_i)$ is Renyi's dimension of the measure $\mu_i$. Now, by exploiting the fact we deal with a relative entropy and the block approximation, we can get rid of the Renyi's dimension of one measure which can be the one concentrated on the higher dimensional sub-manifold (in this case $\mu_m$ in order to simplify notation). Thus, since the $\mu_i$ are concentrated on submanifold of dimension $d_i$ we have that $H_\delta(\mu_i) \leq d_i \log(1/\delta) + C$, with $C > 0$, and by choosing $\delta = \varepsilon$ we obtain the desired result. Notice that in this case the bound will depend only on the dimension of the support of the marginals **and not on the optimal transport plan**. Besides, notice that by taking $m = 2$ and $d_1 = d_2 = d$, one easily retrieves [CPT23, Proposition 3.1]. More involved, it is the case in which we consider a locally semi-concave cost functions. In particular we can use an integral variant of Alexandrov Theorem which is proved in [CPT23].

**Lemma 1.4.2 ([CPT23, Lemma 3.6])**. Let $f: \Omega \to \mathbb{R}$ be a $\lambda$-concave function defined on a convex open set $\Omega \subseteq \mathbb{R}^d$, for some $\lambda \geq 0$. There exists a constant $C \geq 0$ depending only on $d$ such that:

$$ \int_{\Omega} \sup_{y \in B_r(x) \cap \Omega} |f(y) - (f(x) + \langle \nabla f(x), (y-x) \rangle)| dx \leq Cr^2 \mathcal{H}^{d-1}(\partial \Omega)([f]_{\mathbb{C}^{0,1}(\Omega)} + \lambda \operatorname{diam}(\Omega)), \quad (1.22) $$

where $[f]_{\mathbb{C}^{0,1}(\Omega)}$ denotes the Lipschitz constant of $f$.

$$ H_{\delta}(\mu) = \inf \left\{ \sum_{n \in \mathbb{N}} \mu(A_n) \log(1/\mu(A_n)) \mid \forall n, \operatorname{diam}(A_n) \leq \delta, \text{ and } X = \bigcup_{n \in \mathbb{N}} A_n \right\}, $$

where the infimum is taken over countable partitions $(A_n)_{n \in \mathbb{N}}$ of $X$ by Borel subsets of diameter less than $\delta$, and we define the lower and upper entropy dimension of $\mu$ respectively by:

$$ \underline{\dim}_R(\mu) := \liminf_{\delta \to 0^+} \frac{H_\delta(\mu)}{\log(1/\delta)}, \quad \overline{\dim}_R(\mu) := \limsup_{\delta \to 0^+} \frac{H_\delta(\mu)}{\log(1/\delta)}. $$

<footnote n="1">¹Rényi dimension (following [You82]): If $\mu$ is a probability measure over a metric space $X$, we set for every $\delta > 0$,</footnote>
## 1.4. RATE OF CONVERGENCES FOR ENTROPIC MULTI-MARGINAL OPTIMAL TRANSPORT

However, since we are considering (and exploiting) probability measures concentrated on smooth sub-manifolds of dimension $d_i$, we have to carefully define local semi-concavity in order to apply the previous Lemma:

**Definition 1.4.3.** A function $f: X \to \mathbb{R}$ defined on a $\mathbb{C}^2$ sub-manifold $X \subseteq \mathbb{R}^N$ of dimension $d$ is locally semiconcave if for every $x \in X$ there exists a local chart (i.e. a $\mathbb{C}^2$ diffeomorphism) $\psi: U \to \Omega$ where $U \subseteq X$ is an open neighborhood of $x$ and $\Omega$ is an open convex subset of $\mathbb{R}^d$, such that $f \circ \psi^{-1}$ is $\lambda$-concave for some $\lambda \in \mathbb{R}$, meaning $f \circ \psi^{-1} - \lambda \frac{|\partial^2}{2}$ is concave on $\Omega$.

We obtain then the following theorem

**Theorem 1.4.4.** Let $c$ be locally semiconcave and assume that for every $i \in \{1, \dots, m\}$, $X_i \subseteq \mathbb{R}^N$ is a $\mathbb{C}^2$ sub-manifold of dimension $d_i$ and $\mu_i \in L^\infty(\mathcal{H}_{X_i}^{d_i})$ is a probability measure compactly supported in $X_i$. Then there exists constants $\varepsilon_0, C^* \ge 0$ such that for $\varepsilon \in (0, \varepsilon_0]$

$$ \mathrm{MOT}_{\varepsilon}^{c} \le \mathrm{MOT}_{0}^{c} + \frac{1}{2} \left( \sum_{i=1}^{m} d_i - \max_{1 \le i \le m} d_i \right) \varepsilon \log(1/\varepsilon) + C^* \varepsilon. \quad (1.23) $$

The proof consists in better estimating the term $\int cd(\gamma_{\delta} - \gamma_{0})$ in (1.21). Indeed we have that

$$ \int cd(\gamma_{\delta} - \gamma_{0}) = \int Ed\gamma_{\delta}, $$

where $E = c - \bigoplus_{i=1}^{m} \varphi_i$ with $\varphi_i$ being the optimal potentials for the un-regularized problem. By exploiting the fact that $c \in \mathbb{C}^2$, the local semiconcavity of the potentials (inherited by the cost function), the fact that $\mu_i \in L^\infty$ and Lemma 1.4.2 we can improve the upper bound of a factor 1/2.

### 1.4.2 Lower bound

We are going to establish a lower bound in the same form as the fine upper bound of 1.4.4, the dimensional constant being this time related to the signature of some bilinear forms by using Theorem 1.1.1 and following ideas from the proof of it proposed in [Pas12]. Let us consider the following signature condition:

$$ \text{for every } \boldsymbol{x} \in \mathbf{X}, \quad d_c^{+}(\boldsymbol{x}) \ge \kappa \quad \text{where} \quad d_c^{\pm}(\boldsymbol{x}) := \max \{d^{\pm}(g)(\boldsymbol{x}) \mid g \in G_c\}. \quad (\text{PS}(\kappa)) $$

Then we can state the main result concerning the lower bound.

**Theorem 1.4.5.** Let $c \in \mathbb{C}^2(\mathbf{X})$ and assume that for every $i \in \{1, \dots, m\}$, $X_i \subseteq \mathbb{R}^d$ is a $\mathbb{C}^2$ sub-manifold of dimension $d_i$ and $\mu_i \in L^\infty(\mathcal{H}_{X_i}^{d_i})$ be a probability measure compactly supported in $X_i$. If $(\text{PS}(\kappa))$ is satisfied, then there exists a constant $C_* \in [0, \infty)$ such that for every $\varepsilon > 0$,

$$ \mathrm{MOT}_{\varepsilon}^{c} \ge \mathrm{MOT}_{0}^{c} + \frac{\kappa}{2} \varepsilon \log(1/\varepsilon) - C_{*} \varepsilon. \quad (1.24) $$

As already mentioned at the beginning of this section the proof of this result relies on the dual formulation of $(\mathrm{MOT}_{\varepsilon})$

$$ \mathrm{MOT}_{\varepsilon}^{c} \ge \mathrm{MOT}_{0}^{c} - \varepsilon \log \int_{\mathbf{X}} e^{-\frac{E(x_1, \dots, x_m)}{\varepsilon}} d \otimes_{i=1}^{m} \mu_i(x_i), $$
where $E(x_1, \dots, x_m) = c(x_1, \dots, x_m) - \bigoplus_{i=1}^m \varphi_i(x_i)$ is the duality gap and $(\varphi_1, \dots, \varphi_m)$ are Kantorovich potentials for the un-regularized problem (MOT$_\varepsilon$) with $\varepsilon = 0$. By using the singular values decomposition of a well-chosen bilinear form obtained as an average of mixed second derivatives of the cost and the signature condition introduced above we are able to prove that E detaches quadratically from the set $\{E = 0\}$ that is

$$ \frac{E(\mathbf{x}') + E(\mathbf{x})}{2} \geq |u^+(\mathbf{x}') - u^+(\mathbf{x})|^2 - |u^-(\mathbf{x}') - u^-(\mathbf{x})|^2 - \eta(r) |u(\mathbf{x}') - u(\mathbf{x})|^2 \quad (1.25) $$

for $\eta(r) \ge 0$, $\eta \to 0$, and where $u = (u^0, u^+, u^-)$ are the local coordinate associate to the singular values decomposition of a suitable bilinear form in $G_c$, as defined in (1.1). This allows us to estimate the previous integral, in a similar way as in [CPT23], namely for some constant $C > 0$ and for every $\varepsilon > 0$, we get

$$ \int e^{-E/\varepsilon} d \otimes_{i=1}^{m} \mu_i \leq C \varepsilon^{\kappa/2}. $$

**Remark 1.4.6.** Theorem 1.4.5 improves the result in [EN22] where only an upper bound depending on the quantization dimension of the solution to the un-regularized problem is provided. Moreover, this slightly more flexible use of Minty's trick compared to [CPT23] allows us to obtain a lower bound also for degenerate cost functions in the two marginals setting.

### 1.4.3 Matching bound

We present here some examples where the upper and lower bounds coincide giving us a development of the entropic cost.

**Example 1.4.7 (Two marginals case).** In previous works [CPT23; EN22] concerning the rate of convergence for the two marginals problem, it was assumed that the cost function must satisfy a non degeneracy condition, that is $D_{x_1x_2}^2 c$ must be of full rank. A direct consequence of our analysis is that we can provide a lower bound (the upper bound does not depend on such a condition) for costs for which the non-degeneracy condition fails. Let $r$ be the rank of $D_{x_1x_2}^2 c$ at the point where the non-degeneracy condition fails, then the signature of $\bar{g}$ at this point is given by $(r, r, 2d - 2r)$ meaning that locally the support of the optimal $\gamma_0$ is at most $2d - r$ dimensional. Thus, the bounds become

$$ \frac{r}{2} \varepsilon \log(1/\varepsilon) - C_* \varepsilon \leq \mathrm{OT}_{\varepsilon}^c - \mathrm{OT}_0^c \leq \frac{d}{2} \varepsilon \log(1/\varepsilon) + C^* \varepsilon, $$

for some constants $C_*, C^* > 0$. Notice that if $D_{x_1,x_2}^2 c$ has full rank then $r = d$ and we retrieve the matching bound results of [CPT23; EN22].

**Example 1.4.8 (Two marginals case and unequal dimension).** Consider now the two marginals case but unequal dimensional (the unequal dimensional OT will be discussed more precisely in chapter 3), that is for example $d_1 > d_2$. Then, if $D_{x_1,x_2}^2 c$ has full rank, that is $r = d_2$, we obtain a matching bound depending only on the lower dimensional marginal

$$ \frac{d_2}{2} \varepsilon \log(1/\varepsilon) - C_* \varepsilon \leq \mathrm{OT}_{\varepsilon}^c - \mathrm{OT}_0^c \leq \frac{d_2}{2} \varepsilon \log(1/\varepsilon) + C^* \varepsilon, $$

for some constants $C_*, C^* > 0$. If $\mu_1$ is absolutely continuous with respect to $H^{d_1}$ on some smooth sub-manifold of dimension $d_1$, then any optimal transport plan would be concentrated on a set of Hausdorff dimension no less than $d_1$, and thus the upper bound given in [EN22, Theorem 3.8] would be $\frac{d_1}{2} \varepsilon \log(1/\varepsilon) + O(\varepsilon)$, which is strictly worse than our estimate.
**Example 1.4.9 (Gangbo-Święch cost and Wasserstein barycenter).** Suppose that $c(x_1, \dots, x_m) = \sum_{i<j} |x_i - x_j|^2$, known as the Gangbo-Święch cost [GS98]. Notice that the cost is equivalent to $c(x_1, \dots, x_m) = h(\sum_{i<j} x_i)$ where $h$ is $\mathbb{C}^2$ and $D^2h < 0$, then the signature of $\bar{g}$ is $((m-1)d, d, 0)$ and we have a matching bound

$$ \frac{1}{2}((m-1)d)\varepsilon \log(1/\varepsilon) - C_*\varepsilon \leq \mathrm{MOT}_{\varepsilon}^c - \mathrm{MOT}_{0}^c \leq \frac{1}{2}((m-1)d)\varepsilon \log(1/\varepsilon) + C^*\varepsilon. $$

Notice now that considering the $\mathrm{MOT}_0^c$ problem with a cost $c(x_1, \dots, x_m) = \sum_i |x_i - T(x_1, \dots, x_m)|^2$, where $T(x_1, \dots, x_m) = \sum_{i=1}^m \lambda_i x_i$ is the Euclidean barycenter, is equivalent to the $\mathrm{MOT}_0^c$ with the Gangbo-Święch cost and the matching bound above still holds. Moreover, the multi-marginal problem with this particular cost has been shown [AC11] to be equivalent to the Wasserstein barycenter, that is $T_{\gamma_0} = \nu$ is the barycenter of $\mu_1, \dots, \mu_m$.

## 1.5 An ODE characterisation of discrete entropic MMOT

In many cases of interest, the cost function $c(x_1, \dots, x_m) = \sum_{i,j=1, i \neq j}^m w(x_i, x_j)$ is given by a sum of two marginal cost functions; when $w(x_i, x_j) = |x_i - x_j|^2$, for instance the multi-marginal problem is equivalent to the well known Wasserstein barycenter problem (see Proposition 4.2 in [AC11]), while the Coulomb cost $w(x_i, x_j) = \frac{1}{|x_i - x_j|}$ plays a central role in the quantum chemistry applications pioneered in [CFK13] and [BDG12].

In a joint work with B. Pass [NP23a] for such pairwise interaction costs, we develop a continuation method which, by introducing a suitable one parameter family of cost functions, establishes a link between the original multi-marginal problem and a simpler one whose complexity scales linearly in the number of marginals. For discrete marginals, we show that, after the addition of an entropic regularization term, the solution of the original multi-marginal problem can be recovered by solving an ordinary differential equation (ODE) whose initial condition is the solution to the simpler problem.

This method is actually inspired by the one introduced in [CGS09] to compute the Monge solution of the two marginal problem, starting from the Knothe-Rosenblatt rearrangement; note, however, note that since we apply this strategy to a regularized problem, our resulting ODE enjoys better regularity than the one in [CGS09], which, in turn, makes it amenable to higher order numerical schemes. The above mentioned differential equation will be derived by differentiating the optimality conditions of the dual problem; in particular by penalizing the constraints with the soft-max function we will obtain a well defined ODE for which existence and uniqueness of a solution can be established.

### 1.5.1 Discrete Entropic Multi-Marginal Optimal Transport problem

From now on we take finite sets $X_i = \{x_i^1, \dots, x_i^N\}$ and discrete measures $\mu_i = \sum_{j=1}^N \mu_i^j \delta_{x_i^j}$. For sake of simplicity we will keep the continuous in space notation, but integral must be understood as finite sums. Consider then the discrete entropic multi-marginal optimal transport problem defined as

$$ \mathrm{MOT}_{\varepsilon}^{c,M}(\mu_1, \dots, \mu_m) := \inf \left\{ \int_X cd\gamma + \varepsilon \mathrm{Ent}(\gamma | \otimes_{i=1}^m \mu_i) \mid \gamma \in \Pi(\mu_1, \dots, \mu_m) \right\}, \quad (\mathrm{MOT}_{\varepsilon}^M) $$
as well as the associated dual problem (written in a slightly different way with respect to the previous section)

$$ \sup \left\{ \sum_{i=1}^{m} \int_{X_i} \varphi_i d\mu_i - \varepsilon \int_{\mathbf{X}} \exp \left( \frac{\sum_i \varphi_i - c(\mathbf{x})}{\varepsilon} \right) d \otimes_{i=1}^{m} \mu_i \right\}, \quad (1.26) $$

where now the $\varphi_i \in \mathbb{R}^{X_i}$ are functions over finite sets. We note in particular that (1.26) is an *un-constrained* finite dimensional concave maximization problem. Solutions may be computed using a multi-marginal version of the Sinkhorn algorithm [Ben+15; Cut13; PC19; Gal16; Chi+16], and one can then recover the optimal $\gamma$ in $(\text{MOT}_{\varepsilon}^M)$ from the solutions $\varphi_1, \dots, \varphi_m$ to (1.26) via the well known formula:

$$ \gamma_x = \exp \left( \frac{\sum_{i=1}^{m} \varphi_i(x_i) - c(\mathbf{x})}{\varepsilon} \right) \mu_{x_1}^1 \mu_{x_2}^2 \dots \mu_{x_m}^m. $$

Here, we are especially interested in cost functions $c(x_1, \dots, x_m)$ involving pair-wise interactions, that is

$$ c(x_1, \dots, x_m) = \sum_{i<j}^{m} w(x_i, x_j). $$

Such costs are ubiquitous in applications: for example, for systems of interacting classical particles in [CFK13; BDG12], c is a pair-wise cost, with $w(x - y) = \frac{1}{|x - y|}$, known as the Coulomb cost. The case where $w(x, y) = |x - y|^2$ is the quadratic distance is well known to be equivalent to the Wasserstein barycenter problem (see Proposition 4.2 in [AC11]), which has a wide variety of applications in statistics, machine learning and image processing, among other areas. Pairwise costs with pair-dependent interactions, that is, costs of the form $c(x_1, \dots, x_m) = \sum_{i<j} w_{ij}(x_i, x_j)$ where the $w_{ij}$ are not necessarily the same for each choice of i and j, also arise in the time discretization of Brenier's relaxation [Bre89] of Arnold's variational formulation of incompressible Euler equation, and in inference problems for probabilistic graphical models [Haa+21]. Our results adapt immediately to such costs. Let us now consider costs $c_{\eta}$ of the form

$$ c_{\eta}(x_1, \dots, x_m) := \eta \sum_{i=2}^{m} \sum_{j=i+1}^{m} w(x_i, x_j) + \sum_{i=2}^{m} w(x_1, x_i). \quad (1.27) $$

It is clear that when $\eta = 1$ we retrieve a pair-wise cost as defined above whereas in the limit $\eta \to 0$ we obtain a cost involving only the interactions between $x_1$ and the other $x_i$ individually. Later on, we will develop an ordinary differential equation that governs the evolution with $\eta$ of the solutions to the regularized dual problem (1.26). In particular Proposition 1 in [NP23a] asserts that the initial condition for that equation (that is, the solutions when $\eta = 0$) can be recovered by solving each of the individual two marginal problems between $\mu_1$ and $\mu_i$.

### 1.5.2 An ODE characterisation of multi-marginal optimal transport

We now turn our attention to developing an ODE for the Kantorovich potentials, working with the regularized discrete problem $(\text{MOT}_{\varepsilon}^M)$ and its dual (1.26) with pairwise cost (1.27), we make the following, standing assumptions throughout this section:

1. (Equal marginals) All the marginals are equal $\mu_i = \rho = \sum_{x \in X} \rho_x \delta_x$, where X is a finite subset.
1.5. AN ODE CHARACTERISATION OF DISCRETE ENTROPIC MMOT

2. (Symmetric cost) The two body cost $w$ is symmetric $w(x, y) = w(x, y)$.

3. (Finite cost) The two body cost function $w: X \times X \to \mathbb{R}$ is everywhere real-valued.

Notice now that although the cost (1.27) at $\eta = 1$ is symmetric in the variables $x_1, x_2, \dots, x_m$, the one at $\eta < 1$ is not. It is, however, symmetric in the variables $x_2, \dots, x_m$; this means that the optimal $\varphi_i$ in (1.26) satisfy $\varphi_i = \varphi_j = \varphi$ for $i, j \ge 2$ and so, setting $\varphi_1 = \psi$, we can rewrite (1.26) as

$$ \inf_{\varphi, \psi \in \mathbb{R}^X} \{\Phi(\varphi, \psi, \eta)\}, \qquad (1.28) $$

where

$$ \Phi(\varphi, \psi, \eta) := -(m-1) \int_X \varphi d\rho - \int_X \psi d\rho + \varepsilon \int_{X^m} \exp \left( \frac{\sum_{i=2}^m \varphi(x_i) + \psi(x_1) - c_\eta(\mathbf{x})}{\varepsilon} \right) d \otimes^m \rho. $$

Since the functional $\Phi(\varphi, \psi, \eta)$ is convex on the set $\{\varphi, \psi: X \to \mathbb{R}\} \approx \mathbb{R}^{2|X|}$, as the sum of a linear and an exponential function, optimal solutions $(\varphi^*, \psi^*)$ can be characterized by the first order optimality conditions $\nabla_{\varphi} \Phi = \nabla_{\psi} \Phi = 0$ which allow to rewrite the optimal potential $\psi^*$ as a function of $\psi^*$. By exploiting this problem (1.28) can be re-written in a more compact form

$$ \inf_{\varphi: X \to \mathbb{R}} \{\tilde{\Phi}(\varphi, \eta)\}, \qquad (1.29) $$

where

$$ \tilde{\Phi}(\varphi, \eta) := -(m-1) \int_X \varphi d\rho + \varepsilon \int_X \log \left( \int_{X^{m-1}} \exp \left( \frac{\sum_{i=2}^m \varphi(x_i) - c_\eta(\mathbf{x})}{\varepsilon} \right) d \otimes^{m-1} \rho \right) d\rho. $$

It is well known that the solution to (1.26) is unique up to the addition of constants $\varphi_i \mapsto \varphi_i + C_i$ adding to 0, $\sum_{i=1}^m C_i = 0$; thus, solutions to (1.29) are unique up to the addition of a single constant, $\varphi \mapsto \varphi + C$. We therefore impose the normalization

$$ \varphi(x_0) = 0 \qquad (1.30) $$

for all $\eta \in [0, 1]$ and a fixed $x_0 \in X$. The problem (1.29), restricted to $\varphi$'s satisfying (1.30) then has a unique solution; the function $\tilde{\Phi}(\cdot, \eta)$ is strictly convex when restricted to this set, and the solution $\varphi_\eta^*$ can be characterized by the optimality condition $\nabla_{\varphi} \tilde{\Phi}(\varphi^*, \eta) = 0$. Our numerical method consists then in solving an ODE for the evolution of $\varphi(\eta) = \varphi_\eta^*$ (we consider now a curve of optimal potential $\varphi_\eta^*$) obtained by differentiating the above mentioned optimality condition with respect to $\eta$. It turns out that $\varphi$ can be fully characterized as the solution to the following Cauchy problem:

$$ \left\{ \begin{array}{l} \displaystyle \frac{d\varphi}{d\eta}(\eta) = -[D_{\varphi, \varphi}^2 \tilde{\Phi}(\varphi(\eta), \eta)]^{-1} \frac{\partial}{\partial \eta} \nabla_{\varphi} \tilde{\Phi}(\varphi(\eta), \eta), \\ \varphi(0) = \varphi_w, \end{array} \right. \qquad (1.31) $$

where, the initial value $\varphi(0)$ of $\varphi$ when $\eta = 0$ coincides with the optimal potential $\varphi_w$ for the two marginal optimal transport problem with cost $w$.

The following well-posedness theorem then holds.
**Theorem 1.5.1.** Let $\varphi(\eta)$ be the solution to (1.29) for all $\eta \in [0, 1]$. Then $\eta \mapsto \varphi(\eta)$ is $\mathbb{C}^1$ and is the unique solution to the Cauchy problem (1.31).

The main step of the proof consists in showing that the pure second derivatives with respect to $\varphi$ as well as the mixed second derivatives with respect to $\varphi$ and $\varepsilon$ exist and are Lipschitz, and the Hessian with respect to $\varphi$ is invertible. In particular thanks to the boundedness of the cost function, namely $||c_{\eta}||_{\infty} \le M$, allows to prove that the potentials (the techniques we used were largely by inspired by Carlier's beautiful work on the linear convergence of multi-marginal Sinkhorn [Car22]) are also bounded and one can restrict the study of the well-posedness of the ODE on the set

$$U := \{\varphi \in \mathbb{R}^X \mid \varphi_{x_0} = 0, \|\varphi\|_{\infty} \le 4M\}. \qquad (1.32)$$

On this set the functional $\tilde{\Phi}$ is now strongly concave and this ensures that the Hessian is invertible.

### 1.5.3 The algorithm

The algorithm consists in discretizing (1.31) by an explicit Euler scheme (notice that one could also use some high order method for the ODEs). Let $h$ be the step size and set $\varphi(0) = \varphi_w$ the solution of a 2 marginal problem with cost $w$, then the $\varphi$ can be defined inductively as detailed in 1.

<caption>**Algorithm 1** Algorithm to compute the $\varphi$ via explicit Euler method</caption>

Require: $\varphi(0) = \varphi_w$
1: **while** $||\varphi^{(k+1)} - \varphi^{(k)}|| < \text{tol}$ **do**
2: &nbsp;&nbsp;&nbsp;&nbsp;$D^{(k)} := D_{\varphi, \varphi}^2 \tilde{\Phi}(\varphi^{(k)}, kh)$
3: &nbsp;&nbsp;&nbsp;&nbsp;$b^{(k)} := -\frac{\partial}{\partial \eta} \nabla_{\varphi} \tilde{\Phi}(\varphi^{(k)}, kh)$
4: &nbsp;&nbsp;&nbsp;&nbsp;Solve $D^{(k)}z = b^{(k)}$
5: &nbsp;&nbsp;&nbsp;&nbsp;$\varphi^{(k+1)} = \varphi^{(k)} + hz$
6: **end while**

**Remark 1.5.2 (Complexity in the marginals).** We highlight that this ODE approach does not overpass the NP hardness (see [AB21]) of some multi-marginal problems, as the one with the Coulomb cost, and it still suffers the exponential complexity in the number of marginals. However, since the ODE is smooth, one can try apply an high order method (we will do a careful analysis of it later in the section) which can converge quickly in the number of iterations and have a computational time competitive with respect to Sinkhorn. We also note that, if the symmetry assumptions are dropped, as we have brifley discussed above, one needs to solve a system of $(N-1)m$ ODEs. The matrix $D^{(k)}$ in Algorithm 1 is then be $(N-1)m \times (N-1)m$ instead of $m \times m$, increasing the leading term in the complexity by the relatively manageable multiplicative factor $(N-1)^2$. Numerical results for one such case are presented below.

In Figures 1.2 and 1.1 we show the numerical result obtained with $\varepsilon = 0.006$, $h = 1/1000$, $m = 3$, the uniform measure on $[0, 1]$ uniformly discretized with 100 gridpoints and the pairwise interaction $w(x, y) = -\log(0.1 + |x - y|)$. Notice that since we have developed our continuation method by the entropic regularization of optimal transport, we can easily reconstruct the optimal (regularized) plan at each time $k$ by using the potential $\varphi^{(k)}$. Moreover, it is interesting to notice
that the optimal plan at each time step of the ODE stays deterministic (taking into account the blurring effect of the entropic regularization); that is, it is concentrated on a low dimensional structure.

<figure><img src="image_2.png" /><figcaption>Figure 1.1: (Log cost) potential $\varphi(\eta)$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 1.2: (Log cost) support of the coupling $\gamma_{1,2}^\eta$.</figcaption></figure>

**Remark 1.5.3** (Intermediate step of the ODE approach). It is interesting to highlight that each *ith-step* of the ODE approach *actually returns the optimal solution to the multi-marginal problem with the cost* $c_{kh}$. This, in particular, implies that one can also use this approach in order to *retrieve the solutions of many multi-marginal problems with different costs, by choosing a suitable* $c_\eta$, *instead of solving these problems individually using Sinkhorn*.

We are finally interested in comparing the ODE approach and the Sinkhorn algorithm in terms of performance. In order to do this we consider the optimal solution of the regularized dual problem obtained with a gradient descent with backtracking and take it as the reference solution to compute the relative error $\frac{\|\varphi - \varphi_{ref}\|_{\infty}}{\|\varphi_{ref}\|_{\infty}}$. Concerning the ODE, since we have already remarked above that it is smooth, we consider different high order methods such as 3rd, 5th and 8th order Runge-Kutta methods. By looking at Table 1.1, it is clear that all the methods achieve almost the same relative error, but the number of iterations to reach it as well as the CPU time in seconds slightly differs. In particular we remark that a 3rd order RK is faster than a Sinkhorn in terms of computational time and iterations but less precise. The other RK methods achieve comparable precision to Sinkhorn with less iterations but the computational cost at the step of the ODE becomes now quite onerous implying a significant increase (especially 8th RK) in terms of CPU time.
<table><caption>Table 1.1: Comparison between the ODE approach and Sinkhorn for the uniform density and 400 gridpoints</caption><thead><tr><th></th><th>3rd RK</th><th>5th RK</th><th>8th RK</th><th>Sinkhorn</th></tr></thead><tbody><tr><td>relative error</td><td>1.47 × 10<sup>-5</sup></td><td>7.8 × 10<sup>-6</sup></td><td>7.62 × 10<sup>-6</sup></td><td>5.46 × 10<sup>-6</sup></td></tr><tr><td>iterations</td><td>87</td><td>87</td><td>87</td><td>820</td></tr><tr><td>CPU time (sec)</td><td>72.39</td><td>158.9</td><td>385.1</td><td>102.8</td></tr></tbody></table>

**The non-symmetric ODE: Euler case** In a series of papers [Bre89; Bre93; Bre99] Brenier proposed a relaxation of the incompressible Euler equation with constrained initial and final data interpreted as a geodesic on the group of measure preserving diffeomorphisms. Brenier's relaxed formulation consists in finding a probability measure over absolutely continuous paths which minimizes the average kinetic energy. In this framework the incompressibility is encoded by an additional constraint that at each time t, the distribution of position need be uniform. If we consider a uniform discretization of [0, T] (where T is the final time) with m steps in time, we recover a multi-marginal formulation of the Brenier principle with the specific cost function

$$c(x_1, \dots, x_m) = \frac{m^2}{2T^2} \sum_{i=1}^{m-1} |x_{i+1} - x_i|^2$$

representing the discretized, in time, kinetic energy (|·| denotes the standard euclidean norm). The coupling γ is the probability to find a *generalized particle* on the discrete path x₁, ..., x_m. Because the fluid is incompressible, the ith marginal μᵢ of γ is the uniform on the d-dimensional cube [0, 1]ᵈ. For each i ∈ {1, ..., m} the transition probability from time 1 to time i is given by the coupling

$$\gamma_{1,i} = \left( (x_1, \dots, x_m) \mapsto (x_1, x_i) \right)_{\#} \gamma,$$

which represents the probability of finding a *generalized particle* initially at x₁ to be at xᵢ at time i. In order to impose the initial and final constraint we include by penalization, that is by adding a term to the cost function which now reads as

$$c(x_1, \dots, x_m) = \frac{m^2}{2T^2} \sum_{i=1}^{m-1} |x_{i+1} - x_i|^2 + \beta |F(x_1) - x_m|^2,$$

where β > 0 is a penalization parameter in order to enforce the initial-final constraint. F(x₁) represents the prescribed final position of the particle initially at position x₁, and the coupling γ can be interpreted as a (generalization of) a discrete time geodesic between the identity mapping and F on the space of measure preserving maps. For an efficient implementation of Sinkhorn in this case we refer the reader to [BCN19; Nen16; Ben+15]. If we consider now the ODE setting, we have now to deal with a non symmetric case (the cost is not symmetric in the marginals anymore) and so to solve a system, still well posed, of ODEs. In particular we consider the following c_η cost

$$c_{\eta}(x_1, \dots, x_m) = \frac{m^2}{2T^2} |x_2 - x_1|^2 + \eta \left( \frac{m^2}{2T^2} \sum_{i=2}^{m-1} |x_{i+1} - x_i|^2 \right) + \beta |F(x_1) - x_m|^2.$$

For the numerical simulations we took 100 gridpoints discretization of [0, 1], m = 9 time marginals, β = 20 and ε = 0.002. In this case, we solved the ODE system by using a 5th
order Runge-Kutta method and $h = 1/100$. In Figures 1.3 and 1.4 we plot the transition couplings for the initial-final for two different initial-final configurations: $F(x) = 1 - x$ and $F(x) = (x + 1/2) \mod 1$.

<figure><img src="image_4.png" /><figcaption>Figure 1.3: Transition plan $\gamma_{1,i}$ for $F(x) = 1 - x$ of Brenier relaxed formulation of incompressible Euler equations.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 1.4: Transition plan $\gamma_{1,i}$ for $F(x) = (x + 1/2) \mod 1$ of Brenier relaxed formulation of incompressible Euler equations.</figcaption></figure>

## 1.6 Perspectives

Let us mention in this section some research perspectives:

* **Particle discretization of MOT<sup>c</sup>:** we want study the behaviour of a simple **particle discretization** of the problem, where the solution is a sum of a finite number of Dirac masses whose positions are free to move, that is $\gamma = \sum_{i=1}^{N} \delta_{X_i}$ where $X_i \in \mathbb{R}^{dm}$. The constraints could be enforced via penalization, using optimal transport data attachment terms. The problem we want to solve then takes the following form

$$ \min_{X^1, \dots, X^m} \sum_{i=1}^{N} c(X_i^1, \dots, X_i^m) + \lambda \left( \sum_{j=1}^{m} \mathrm{OT}^{d^2}(\mu_j, e_{j\#}\gamma) \right), \quad (1.33) $$

where we take the quadratic cost $d^2(x, y) = |x - y|^2$ in the penalization terms and $X_i^j$ is the $j$-th component of the vector $X_i$. This discretization makes the discrete problem non-convex, but one could hope to study the convergence of global minimizers of MOT<sup>c</sup>
depending on the dimension of the support of the solution, on the size of the discretization and the strength of the penalization. Establishing convergence results for this discretization in the spirit of Γ-convergence should then be within hand, using similar arguments as [Sar22], and will extend the result in [BGV19] for the case of Wasserstein splines. Notice that a coarser discretization, thanks to a better understanding of the structure of the solution, will make the numerical method competitive with respect to the actual ones: in the entropic case the dimension of the support of the solution is naturally high dimensional since the regularization effect, meaning that a finer discretization is needed. A recent result of [MSS21] makes us optimistic that despite the non-convexity of this particle discretization, it would be possible to prove (at least in some cases) that a simple gradient descent algorithm is sufficient to reach a low energy configuration. This is a difficult question, that would be of prime interest to many communities, from optimal transport to machine-learning. Moreover, we also plan to study the discretization in the case in which the OT penalization terms are replaced by the entropic counterpart $OT_{\epsilon}^{c}(\mu_1, \mu_2)$. Finally we would like to implement this algorithm and experiment it numerically on some applications as the risk estimation problem presented in [IL15]. A second step will consist in testing (and comparing performances) other numerical methods such as the **Sinkhorn algorithm** based on the entropic regularization of Optimal Transport [BCN19; Ben+15; Ben+18] or the ones arising from the approximation of **optimal transport problems with marginal moments constraints** [Alf+21; ACE22]. In particular we are also interested in studying a recent approach based on genetic column generation method proposed in [FP22; FSV22]: convergence for multi-marginal problem is still an open question as well as a possible link with Frank-Wolfe method.

* **Entropic optimal transport and Mean Field Games:** Establish convergence rates for optimizers of the $\epsilon$-regularized primal problem $OT_{\epsilon}^{c}$ and its dual as $\epsilon \to 0$. While convergence (up to subsequence) is easily obtained for compactly supported measures by simple compactness arguments, obtaining quantitative rates is substantially more difficult. The fully discrete case is covered by the general results of [CM94; Wee18] and the semi-discrete case by [ANS22; Del22]. We propose to focus on the continuous case, which is not so well understood, and establish sharp convergence rates for the primal problem (e.g. in terms of Wasserstein distance) and for the dual problem, with assumptions on the cost $(\mathbb{C}^2$ regularity, twist condition in particular) and on the marginals but as few assumptions as possible on the unregularized optimizers (in contrast to what is usually done in the literature). Study the asymptotics as $\epsilon \to 0$ in more general situations including the multi-marginal case and non-twisted or degenerate costs. A first direction is to generalize the convergence rates of the entropic cost established in [CPT23] to the multi-marginal case with possibly degenerate costs. A second direction, as another development from [CPT23], is to tackle the question of the convergence of entropic plans.

Furthermore, consider the dual problem to $OT_{\epsilon}^{c}$, that is

$$ \sup_{\varphi_1, \varphi_2} \int \varphi_1(x_1) d\mu_1 + \int \varphi_2(x_2) d\mu_2 - \epsilon \int e^{\frac{\varphi_1(x_1) + \varphi_2(x_2) - c(x_1, x_2)}{\epsilon}} d\mu_1 d\mu_2, \quad (1.34) $$

by optimality condition one can derive the Sinkhorn's algorithm

$$ \varphi_i^{n+1}(x_i) = -\log \left( \int e^{\frac{\varphi_j(x_j) - c(x_1, x_2)}{\epsilon}} d\mu_j \right), \quad \text{for } i=1,2 \text{ and } i \neq j. \quad (1.35) $$
Linear convergence for the multi-marginal case has been recently established in [Car22]
however the same result for a generalized Sinkhorn in order to solve variational second
order Mean Field Games (see [Ben+18] for more details) is still an open question. In this
case we would like to establish it by starting with simple congestion functionals and then
relax the assumptions in order to treat more general case. Moreover, we plan to extend the
variational formulation to the case of a **multi-population** systems in order to study the
dynamic of several populations. We also want to tackle the convergence for the variational
formulation of second order mean field games [Ben+18] when the noise vanishes.

Finally, for the entropic transport we have developed an ODE approach for the case of
pair-wise cost function, we would like to extend it to a more general cost functions. Notice
that in this case one can easily take $c_{\eta}(x_1, \dots, x_m) = \eta c(x_1, \dots, x_m)$ and having the trivial
null solutions as initial condition for the Cauchy problem.
# Chapter 2

# Generalized Optimal Transport in
Mathematical Physics

This chapter summarizes the main results obtained for some applications of Optimal Transport in Mathematical Physics. In particular we will focus on (1) a Statistical Mechanics extension of multi-marginal transport, namely the Grand Canonical OT, which allows to deal with the case of a varying number of marginals; (2) an approximation of Lieb's functional which can be seen as a quantum regularization of optimal transport.

All the results are contained in the following papers

Section 2.2 S. Di Marino, M. Lewin, and L. Nenna. *Grand-Canonical Optimal Transport*. 2022. URL: http://arxiv.org/abs/2201.06859. preprint

Section 2.3 V. Ehrlacher and L. Nenna. *A Sparse Approximation of the Lieb Functional with Moment Constraints*. 2023. URL: http://arxiv.org/abs/2306.00806. preprint

## 2.1 Optimal Transport in Density Functional Theory and beyond

Density Functional Theory (see the seminal paper by Hohenberg and Kohn [HK64]) attempts to describe all the relevant information about a many-body quantum system at ground state in terms of the one electron density $\rho$. Following Levy and Lieb's approach [Lev79; Lie83], the ground state energy can be rephrased as the following variational principle involving only the electron density

$$ \varepsilon_0[v] = \inf_{\substack{\rho \in A^m \\ \int_{\mathbb{R}^3} v(x) d\rho < +\infty}} \left\{ F_{LL}[\rho] + \int_{\mathbb{R}^3} v(x) d\rho \right\}, $$

where $A^m = \{\rho \in L^1(\mathbb{R}^3) : \rho \ge 0, \sqrt{\rho} \in H^1, \rho(\mathbb{R}^3) = m\}$ is the set of admissible densities, $v$ is an external potential and the Levy-Lieb functional $F_{LL}$ is defined as

$$ F_{LL}[\rho] := \inf_{\substack{\Psi \in \mathcal{H}_1^m \\ \rho_{\Psi} = \rho}} \frac{1}{2} \int_{\mathbb{R}^{3m}} |\nabla \Psi|^2 + \int_{\mathbb{R}^{3m}} V |\Psi|^2, \qquad (2.1) $$

where
(i) $\mathcal{H}_1^m := \bigwedge_{i=1}^m H^1(\mathbb{R}^3)$ is the set of admissible electronic wavefunctions for a system of $m$ electrons with finite kinetic energy, that is the set of antisymmetric functions of $H^1(\mathbb{R}^{3m})$;

(ii) for any $\Psi \in \mathcal{H}_1^m$ and $x \in \mathbb{R}^3$, $\rho_\Psi$ is the electronic density associated to the wavefunction $\Psi$, namely the real-valued function defined over $\mathbb{R}^3$ as follows:

$$ \forall x \in \mathbb{R}^3, \quad \rho_{\Psi}(x) := m \int_{(\mathbb{R}^3)^{m-1}} |\Psi(x, x_2, \dots, x_m)|^2 dx_2 \dots dx_m. $$

(iii) the function $V : (\mathbb{R}^3)^m \to \mathbb{R}_+ \cup \{+\infty\}$ is the electron-electron Coulomb interaction potential

$$ V(x_1, \dots, x_m) := \sum_{i<j} \frac{1}{|x_i - x_j|}. $$

The Levy-Lieb functional is indeed the lowest possible (kinetic plus interaction) energy of a quantum system having the prescribed density $\rho$. This universal functional is the central object of Density Functional Theory, since knowing it would allow one to compute the ground state energy of a system with any external potential $v$. For a complete review on it we refer the reader to [LLS23]. However, it turns out that $F_{LL}$ is not convex, it is therefore convenient to look at a convexification proposed by Lieb [Lie83] where the minimization is performed over the set of mixed states instead of the set of pure ones as in (2.1). More precisely, we consider here the alternative minimization problem

$$ F_L[\rho] := \inf_{\substack{\Gamma \in \mathcal{S}_1^+(\mathcal{H}_0^m) \\ \rho_\Gamma = \rho}} \mathrm{Tr}(H_m \Gamma), \qquad (2.2) $$

where $H_m := -\frac{1}{2}\Delta + V$ is a self-adjoint operator on $\mathcal{H}_0^m := \bigwedge_{i=1}^m L^2(\mathbb{R}^3)$ with domain $D(H_m) = \mathcal{H}_2^m := \bigwedge_{i=1}^m H^2(\mathbb{R}^3)$, $\mathcal{S}_1^+(\mathcal{H}_0^m)$ denotes the set of trace-class self-adjoint non-negative operators on $\mathcal{H}_0^m$. For all $\Gamma \in \mathcal{S}_1^+(\mathcal{H}_0^m)$, there exists an orthonormal basis $(\Psi_i)_{i \in \mathbb{N}^*}$ of $\mathcal{H}_0^m$ and a non-increasing sequence $(\alpha_i)_{i \in \mathbb{N}^*}$ of non-negative numbers such that

$$ \Gamma = \sum_{i=1}^{+\infty} \alpha_i |\Psi_i\rangle\langle\Psi_i|, \qquad (2.3) $$

using so-called bra-ket notation. Then, the associated electronic density $\rho_\Gamma$ is defined as follows: for all $x \in \mathbb{R}^3$,

$$ \rho_\Gamma(x) := m \sum_{i=1}^{+\infty} \alpha_i \int_{(\mathbb{R}^3)^{m-1}} |\Psi_i(x, x_2, \dots, x_m)|^2 dx_2 \dots dx_m = \sum_{i=1}^{+\infty} \alpha_i \rho_{\Psi_i}(x). $$

We know that there exist positive constants $\varepsilon, D > 0$ such that $H_m + D \ge \varepsilon(-\Delta + \mathrm{Id})$ (in the sense of self-adjoint operators on $\mathcal{H}_0^m$). We also denote by $\mathcal{S}_{1,1}(\mathcal{H}_0^m)$ the set of self-adjoint operators $\Gamma$ on $\mathcal{H}_0^m$ with finite kinetic energy, i.e. such that $\mathrm{Tr}(|H_m + D|^{1/2}\Gamma|H_m + D|^{1/2}) < +\infty$.

**Remark 2.1.1.** It can then be easily checked that, $\Gamma \in \mathcal{S}_{1,1}(\mathcal{H}_0^m)$ if and only if $\Gamma \in \mathcal{S}_1(\mathcal{H}_0^m)$ and $\mathrm{Tr}(H_m \Gamma) < +\infty$. Then, if $\Gamma$ admits an eigendecomposition of the form (2.3), necessarily $\Psi_i \in \mathcal{H}_1^m$ as soon as $\alpha_i > 0$.
It is well-known, see for instance [Lie83], that the infimum in (2.1) and (2.2) is attained.

**Remark 2.1.2** (Convexification). It is worth highlighting that $F_L$ is indeed the convexification of $F_{LL}$ in the sense that

$$F_L[\rho] = \inf_{\substack{\forall i \ge 1, \alpha_i \ge 0, \rho_i \in A^m \\ \sum_{i=1}^{+\infty} \alpha_i = 1 \\ \sum_{i=1}^{+\infty} \alpha_i \rho_i = \rho}} \sum_{i=1}^{+\infty} \alpha_i F_{LL}[\rho_i]$$

It is useful noticing that $F_L$ admits a dual problem.

**Theorem 2.1.3** ([Lie83]). Duality holds in the sense that

$$F_L[\rho] = \sup_{\substack{v \in L^{3/2}(\mathbb{R}^3) + L^\infty(\mathbb{R}^3) \\ H_m^v \ge 0}} \left\{ \int_{\mathbb{R}^3} v(x) \rho(x) dx \right\}, \quad (2.4)$$

where

$$H_m^v = H_m - \sum_{i=1}^{m} v(x_i).$$

The constraint in (2.4) has to be understood in the sense of self-adjoint operators, namely for all $\Psi \in \mathcal{H}_1^m$, $\langle \Psi | H_m^v | \Psi \rangle \ge 0$.

**Remark 2.1.4.** Notice that (2.4) is very similar to the dual problem introduced defined in (MD).

It is possible to introduce an effective semi-classical parameter $\epsilon = \hbar^2$ by scaling the density $\rho$. Namely, for $\rho_\epsilon = \epsilon^3 \rho(\epsilon x)$ we have

$$\begin{align}
\frac{F_L[\rho_\epsilon]}{\epsilon} &= \inf_{\substack{\Gamma \in \mathcal{S}_1^+(\mathcal{H}_0^m) \\ \rho_\Gamma = \rho}} \mathrm{Tr} \left( -\frac{\epsilon}{2} \Delta + V \right) \Gamma, \nonumber \\
\frac{F_{LL}[\rho_\epsilon]}{\epsilon} &= \inf_{\substack{\Psi \in \mathcal{H}_1^m \\ \rho_\Psi = \rho}} \frac{\epsilon}{2} \int_{\mathbb{R}^{3m}} |\nabla \Psi|^2 + \int_{\mathbb{R}^{3m}} V |\Psi|^2 \tag{2.5}
\end{align}$$

In the limit $\epsilon \to 0$, the above functionals converge (see [Lew18] for Lieb functional and [CFK13; CFK18; BD17] for the Lévy-Lieb ones) to the multi-marginal optimal transport problem with Coulomb cost

$$\inf_{\gamma \in \Pi(\rho)} \int_{\mathbb{R}^{3m}} \sum_{i<j} \frac{1}{|x_i - x_j|} d\gamma(x_1, \dots, x_m) \quad (2.6)$$

where now the marginals are all equal to the given electron density $\rho$.

**Remark 2.1.5.** Notice that now the mass of $\rho$ is not equal to 1 anymore even if by a slight abuse of notation we will still define the OT problems on the set of probability measures. However the fact that the mass is equal to the number of particles/electrons m of the system does not affect the minimizers but only the value of the energy.

In Theoretical Chemistry literature the 0-limit regime is also known as the Strictly Correlated Electrons (SCE) state and it has been firstly introduce in [SGS07] and it is used in order to approximate the electron-electron repulsion term in Lévy-Lieb/Lieb functionals. The literature about the SCE approximation is growing considerably. Recent developments include results on
the existence and non-existence of Monge-type solutions (e.g., [CDD15; CS16; CFK13; BDG12; BDK21; CD15]), structural properties of Kantorovich potentials (e.g., [CDS19; BCD18]), and design of efficient computational algorithms (e.g., [BCN16; Nen16; KY19; Alf+21; Kho+20; ACE22; FP22; FSV22; FP23]).

In this chapter we mainly focus on two different generalization of multi-marginal optimal transport:

* **Grand Canonical Optimal Transport** This model takes its roots in Statistical Physics, where it usually goes under the name *Grand-Canonical* [Rue99]. It has recently been used for Coulomb and Riesz costs where it naturally occurs in the large-m limit of the multi-marginal problem [LLS19; LLS18; LLS23; CP19]. A truncated version has been studied in [Bou+21a; Bou+21b]. Its entropic regularization is well known in the literature [CCL84; CC84; JKT23] and plays a central role in the density functional theory of inhomogeneous classical fluids at positive temperature [Eva79]. Here we will define both Grand Canonical OT and entropic grand canonical OT, study the existence and characterization of solutions as well as develop the duality theory.

* **Quantum Optimal Transport** We provide a sparse approximation of the Lieb functional (2.2) by using some techniques originally developed to solve multi-marginal optimal transport [Alf+21]. In particular we can look at this functional as a (linear) regularization of multi-marginal OT in the quantum framework.

## 2.2 Grand Canonical Optimal Transport

The Grand-Canonical Optimal Transport (GC-OT) problem can be formulated in the Kantorovich form as follows

$$ \mathrm{GCOT}_{0}^{c}(\rho) := \inf_{\gamma \in \Pi_{\mathrm{GC}}} \sum_{n=0}^{\infty} \int_{X^n} c_n \, d\gamma_n, \qquad (2.7) $$

where

$$ \Pi_{\mathrm{GC}}(\rho) := \left\{ \gamma = (\gamma_n)_{n \ge 0} : \gamma_0 \in [0,1], \quad \gamma_n \in \mathcal{M}_{\mathrm{sym}}(X^n) \right. \\ \left. \qquad \gamma_0 + \sum_{n \ge 1} \gamma_n(X^n) = 1, \quad \rho_{\gamma} = \rho \right\}, \quad (2.8) $$

where $\mathcal{M}_{\mathrm{sym}}(X^n)$ denotes the set of symmetric positive measures on $X^n$, $c_n$ is the symmetric cost for the n-marginal problem and $\rho_{\gamma_n}$ is the first marginal of $\gamma_n$. Notice the factor $n$ multiplying the marginal in the constraint involving $\rho$, which accounts for the fact that there are $n$ equal such marginals since all the $\gamma_n$ are symmetric. The family $\gamma = (\gamma_n)_{n \ge 0}$ forms a probability which describes the behavior of some agents whose number is unknown or can vary. In this interpretation $\gamma_n(X^n)$ is the probability that there are $n$ agents and $\gamma_0 \in [0, 1]$ is the one that there is no agent at all. In the GC-OT problem (2.7) only the *average* quantity $\rho$ is fixed and fluctuations of the number of agents are allowed. Solving the minimization problem $\mathrm{GCOT}_0^c(\rho)$ requires in particular to determine the best way to distribute the number of agents through the measures $\gamma_n$, in order to reproduce the given average $\rho$, depending on the corresponding costs $c = (c_n)_{n \ge 0}$.

In order to motivate the problem, let us first explain why it is natural to let the number of marginals vary, even if we start with a system which has a well defined number of agents
m. Consider a symmetric Borel probability measure $\mathcal{P}$ over $X^m$, where $X \subset \mathbb{R}^n$ is any given Borel and $m \ge 2$ is typically a large number. In applications, $\mathcal{P}$ describes $d$ properties of $m$ agents, taking their values in $X$ (e.g. their space location along the $d$ coordinates in $\mathbb{R}^d$), with $\mathcal{P}(A_1 \times \cdots \times A_m)$ being the probability that agent 1 is in $A_1$, agent 2 is in $A_2$, etc. The symmetry of $\mathcal{P}$ simply means that our agents are all identical and indistinguishable from one another. When $m$ is very large it seems natural to allow it to vary a little, for instance to account for the fact that it can probably not be known exactly. However, grand-canonical states also occur naturally at fixed $m$ when we look at subsystems. Let us explain this now. Let us fix a subset $A \subset X$. Imagine that we would like to ignore what is happening outside of $A$ and only concentrate on what the agents in $A$ are doing. Of course, although there are exactly $m$ agents in total, the number of agents in the subset $A$ can vary from 0 to $m$. A simple counting argument leads us to introducing the measures $\gamma_n$ on $A^n$ defined by

$$ \left\{ \begin{array}{ll} \gamma_0 = \mathcal{P}((X \setminus A)^m) & \text{for } n=0, \\ \gamma_n(B_n) = \binom{m}{n} \mathcal{P}(B_n \times (X \setminus A)^{m-n}) & \text{for } 1 \le n \le m-1, \\ \gamma_m(B_m) = \mathcal{P}(B_m) & \text{for } n=m, \\ 0 & \text{for } n \ge m+1, \end{array} \right. \qquad (2.9) $$

for any $B_n \subset A^n$. In other words, the conditional measure $\gamma_n$ describes what $n$ agents are doing in $A$, in the situation that there are exactly $n$ agents in $A$ and $m-n$ outside. The combinatorial factor in the definition of $\gamma_n$ is because in the expression it was assumed, by symmetry, that the $n$ first are in $A$ whereas the $m-n$ last ones are outside. Then $\gamma_n(A^n)$ is the probability that there are exactly $n$ agents in $A$ and thus

$$ \gamma_0 + \sum_{n=1}^{m} \gamma_n(A^n) = 1. $$

On the other hand, the average number of agents in $A$ is given by

$$ \sum_{n=1}^{m} n \gamma_n(A^n) \in [0, m]. $$

The corresponding density of agents in $A$ is the positive measure defined by

$$ \rho_{\gamma}(B) := \gamma_{1}(B) + \sum_{n=2}^{m} n \gamma_{n}(B \times A^{n-1}) $$

for every $B \subset A$. It is the sum of the marginals with the multiplication factor $n$ seen before. A calculation shows that

$$ \rho_{\gamma}(B) := m \mathcal{P}(B \times X^{m-1}), $$

that is, $\rho_\gamma$ is just the restriction of $m$ times the first marginal of $\mathcal{P}$ to the set $A$. Our conclusion is that although the total number of agents $m$ can be fixed, it is never fixed as soon as we look at a subsystem, which is always represented by a grand-canonical probability $\gamma$. Here we have $\gamma_n \equiv 0$ for $n \ge m+1$ because we start with $m$ agents in total. But if we let $m$ be arbitrarily large with $A$ fixed, then $\gamma$ can in principle have infinitely many non-trivial $\gamma_n$'s, as in (2.7).

Note that there is definitely some loss of information when we look at this $\gamma$ instead of the big probability $\mathcal{P}$ since everything which is happening outside of $A$ has been averaged over. Even
if we would cover $X$ with several $A$'s and look at the corresponding localized probabilities, we would in general not be able to reconstruct $\mathcal{P}$, since the correlation between the different domains is discarded. We believe that approximating a large multi-marginal problem by a collection of smaller local grand-canonical ones is a strategy which might be helpful in practice.

### 2.2.1 The theory of Grand Canonical Optimal Transport

#### Existence

We start to state some results concerning the existence of minimizers for (2.7) as well as the relation with a classical multi-marginal optimal transport problem. Although one can allow any possible cost $c_n$ for a mathematical exercise, in practice the $c_n$ are often related with one another. A large part of our results is stated by using costs describing pairwise interactions between agents/particles, that is all the $c_n$ are expressed in terms of two-agent cost $c_2$ only.

**Definition 2.2.1 (Pairwise costs).** A pairwise grand-canonical cost takes the form

$$c_0 = c_1 = 0, \quad c_n(x_1, \dots, x_n) = \sum_{1 \le j < k \le n} c_2(x_j, x_k) \quad \text{for } n \ge 3. \qquad (2.10)$$

The well-posedness of the grand-canonical problem (2.7) requires specific assumptions on $c_n$ to avoid a collapse due to the possibility of having infinitely many agents in the system. The first assumption we need is the following one

**Definition 2.2.2 (Stability).** Let $X \subset \mathbb{R}^d$. We say that the family of symmetric costs $c = (c_n)_{n \ge 0}$ is stable whenever there exist two constants $A, B \ge 0$ such that $c_n \ge -A - Bn$ on $X^n$ for all $n \ge 0$.

Stability implies that $\text{GCOT}_0^c(\rho) > -\infty$ and so our problem is a well defined minimization problem, which is manifestly convex in $\rho$. Although stability is a good condition for $\text{GCOT}_0^c(\rho)$ to be well defined for all finite measures $\rho$, it is not sufficient to obtain the existence of optimizers. Consider for instance costs which are positive but not large enough for $n \gg 1$, e.g., $c_0 = 0$ and $c_n > 0$ with $||c_n||_{L^\infty(X^n)} = o(n)$. Now using the following trial state

$$\gamma_0 = 1 - \frac{\rho(X)}{m}, \quad \gamma_m = \frac{\rho(X)}{m} \otimes^m \left( \frac{\rho}{\rho(X)} \right), \quad \gamma_n \equiv 0 \quad \forall n \notin \{0, m\} \qquad (2.11)$$

we have that $0 < \text{GCOT}_0^c(\rho) \le \rho(X) ||c_m||_{L^\infty(X^m)}/m \to 0$. Hence $\text{GCOT}_0^c(\rho) = 0$ is never attained for $\rho \ne 0$. This example shows that the cost $c_n$ should be at least as large as $n$ for some configurations to hope to have minimizers. Notice that for a pairwise cost as in the definition above with $c_2(0, 0) > 0$ then we have $c_n(0, \dots, 0) = \frac{n(n-1)}{2}c_2(0, 0)$ which blows up like $n^2$. The following remedy is the adaptation of another classical concept in statistical mechanics [Rue99].

**Definition 2.2.3 (Super-stability).** We say that the family of costs $c = (c_n)_{n \ge 0}$ is super-stable if it is stable and if for any compact set $K \subset \mathbb{R}^d$, there exists $\varepsilon_K > 0$ and $n_K \in \mathbb{N}$ such that

$$c_n(x_1, \dots, x_n) \ge -\frac{n}{\varepsilon_K} + \varepsilon_K \left( \sum_{j=1}^{n} 1_{X \cap K}(x_j) \right)^2 \quad \text{on } X^n \text{ for all } n \ge n_K. \qquad (2.12)$$

Note that such a cost blows up quadratically in the number of agents in any fixed domain $K$. Then the following existence theorem holds
**Theorem 2.2.4 (Existence of optimizers [DLN22])**. Let $X \subset \mathbb{R}^d$ be any Borel set. Let $c = (c_n)_{n \ge 0}$ be a superstable family of lower semi-continuous costs. Then any finite $\text{GCOT}_0^c(\rho)$ admits a minimizer $\gamma^*$. Moreover $\rho \mapsto \text{GCOT}_0^c(\rho)$ is convex and lower semi-continuous for the tight convergence of measures.

The proof relies on the standard direct method of calculus of variations. However let us remark that the super-stability plays a central role in showing that the candidate $\gamma^*$ to be the optimizer has the right average marginal: it is clear that $\rho_{\gamma^*} \le \rho$ and we also have that for any compact $K$

$$ \rho_{\gamma^*}(K) \ge \lim_{k \to +\infty} \left( \rho(K) - \sum_{n \ge M} \rho_{\gamma_n^k} \right), $$

where $\gamma^k$ is the minimizing (sub)sequence. Hence thanks to superstability we get that the tail of the sum in the above equation behaves as $1/\sqrt{M}$. By letting $M \to \infty$ we can conclude that $\rho_{\gamma^*}(K) \ge \rho(K)$ obtaining existence for $\text{GCOT}_0^c(\rho)$.

If we decompose $\rho$ into any possible convex combination of densities with integer masses, then as an immediate consequence of the existence theorem we have that the $\text{GCOT}^c$ is the the convex hull of multi-marginal optimal transport problem that is

$$ \text{GCOT}_0^c(\rho) = \min_{\substack{\rho=\sum_{n \ge 1} \alpha_n \rho_n \\ \rho_n(\bar{X})=n \\ \sum_{n \ge 0} \alpha_n=1}} \sum_{n \ge 0} \alpha_n \text{MOT}_0^c(\rho_n), \quad (2.13) $$

where $\text{MOT}_0^c(\rho_n)$ is the standard $n$-marginal problem with cost $c_n$. Under an extra slightly technical hypothesis on the family of cost

$$ \lim_{\substack{\min_k |y_k| \to \infty \\ \min_{k \ne \ell} |y_k - y_\ell| \to \infty}} c_{n+m}(x_1, \dots, x_n, y_1, \dots, y_m) = c_n(x_1, \dots, x_n), \quad (2.14) $$

one can also prove that $\text{GCOT}^c$ is also the weak closure of $\text{MOT}_0^c$.

**Theorem 2.2.5 (Weak lower semi-continuous envelope [DLN22])**. Let $\rho$ be any finite measure so that $\text{GCOT}_0^c(\rho) < +\infty$. Then there exists a sequence $\rho^k$ such that $N_k := \rho^k(\mathbb{R}^d) \in \mathbb{N}$,

$$ \rho^k \to \rho \quad \text{locally and} \quad \lim_{k \to \infty} \text{MOT}_0^c(\rho^k) = \text{GCOT}_0^c(\rho). $$

The fact that we only have local convergence is due to the possible unboundedness of the number $N_k$ of agents in $\rho^k$. This is because we need infinitely many agents to be able to reproduce a $\gamma = (\gamma_n)_{n \ge 0}$ with $\gamma_n \ne 0$ for arbitrarily large $n$. If we know that there exists a minimizer which satisfies $\gamma_n \equiv 0$ for $n$ large, then the local convergence can be replaced by weak convergence in duality with $C_0(\mathbb{R}^d)$.

### Truncation of the support

When $\text{GCOT}_0^c(\rho)$ has a minimizer $\gamma$, a natural question is to ask how many agents are necessary to minimize the given grand-canonical cost $c$, that is, how many of the $\gamma_n$'s are non zero. Let us firstly introduce the following definition of support for a grand canonical probability measure.
**Definition 2.2.6 (Support).** We call

$$ \mathrm{Supp}(\gamma) = \{n \ge 0 : \gamma_n \ne 0\} $$

the support in $n$ of a grand-canonical probability $\gamma = (\gamma_n)_{n \ge 0}$ and we say that $\gamma$ has a compact support whenever $\mathrm{Supp}(\gamma)$ is bounded.

Notice now that in numerical simulations we have to truncate the support of the $\gamma$ in order to make $\mathrm{GCOT}^c$ computable that is

$$ \mathrm{GCOT}^{c, \le N}(\rho) := \inf_{\substack{\gamma \in \Pi_{\mathrm{GC}}(\rho) \\ \mathrm{Supp}(\gamma) \subset [0, N]}} \sum_{n \ge 0} \int_{X^n} c_n d\gamma_n. \quad (2.15) $$

Thus, knowing the size of the support of exact minimizers plays a central role in order to suppress or diminish the approximation due to the truncation. Good quantitative estimates are then useful and they will often depend on $\rho$.

**Remark 2.2.7 (Truncated Grand Canonical).** Notice that $\mathrm{GCOT}^{c, \le N}(\rho)$ has been deeply studied in a beautiful paper [Bou+21b] by Bouchitté and coauthors where it was obtained as the weak closure of $\mathrm{MOT}_0^c$ in a similar spirit has in our work. Moreover, in [DLN22][Thm 3.1] we can provide that $\mathrm{GCOT}^{c, \le N}(\rho)$ convergences toward $\mathrm{GCOT}^c(\rho)$ when $N \to +\infty$.

Notice that in our work [DLN22] we discuss several possible conditions on $c_2 > 0$ which imply that the support of minimizers is always compact. There are some overlap between these conditions but we have not found a general simple theory which covers everything. To be more precise, we give a quantitative estimate on the support of minimizers in the following cases

(i) when $X$ is a bounded domain,

(ii) when $1/c_2$ satisfies a triangle-type inequality, which is for instance the case of Riesz interactions $c_2(x, y) = |x - y|^{-s}$ with $s > 0$ in $\mathbb{R}^d$,

(iii) when $c_2$ is asymptotically doubling.

In the Coulomb case $c_2(x, y) = |x - y|^{-1}$ we will show that the length of the support is controlled by $\sqrt{\rho(X)}$ but it can in principle grow with $\rho(X)$. Let us summarise the results concerning the support of minimisers:

* [DLN22][Thm 3.2]: Consider the case where the cost is pairwise and such that $0 < m \le c_2 \le M \le \infty$ and the domain is bounded. Then $\mathrm{GCOT}^c(\rho)$ is finite for any non-negative measure $\rho$ and when $\rho(X) > 1$ a minimizer $\gamma^*$ satisfies

$$ \mathrm{Supp}(\gamma^*) \subset \left[ \frac{m}{M} |\rho(X)|, 1 + \frac{M}{m} (|\rho(X)| - 1) \right]. $$

When $\rho(x) \le 1$ then $\mathrm{GCOT}^c(\rho) = 0$ and the unique minimizers is of the form $\gamma^* = (1 - \rho(X), \rho, 0, \dots)$.

* [DLN22][Thm 3.3]: Consider as before a bounded domain and a pairwise cost with $c_2$ lower semi-continuous and such that $1/c_2$ satisfies the triangle-type inequality

$$ \frac{1}{c_2(x, y)} \le Z \left( \frac{1}{c_2(x, z)} + \frac{1}{c_2(z, y)} \right), \quad \forall x, y, z \in X, $$
for a some positive constant $Z$ then for any non-negative measures $\rho$ such that $\text{GCOT}^c(\rho) < +\infty$ and $\rho(X) > 1$ a minimizer has a compact support

$$ \text{Supp}(\gamma^*) \subset \left[ \frac{\lfloor \rho(X) \rfloor + 1}{2Z + 1}, (2Z + 1)(\lfloor \rho(X) \rfloor - 1) \right]. $$

The above results is indeed inspired by a celebrated paper by Lieb [Lie84] on the maximum number of electrons that a molecule can bind. Moreover, notice that our theorem can be applied to the 3D Coulomb potential or more generally to Riesz costs where $c_2(x, y) = |x - y|^{-s}$, $s > 0$ in any dimension. In these cases we have $Z = \max(1, 2^{s-1})$.

* [DLN22][Thm 3.4]: We also consider the case of an asymptotically doubling cost, namely a condition to quantify how fast $c_2$ can converge to zero. For a cost like this there exists a positive constant $C$ such that for some large enough $R$ we have $m(2r) \ge CM(r)$ for all $r \ge R$, where

$$ m(r) = \inf\{c_2(x, y) : |x - y| \le r\}, \quad M(r) = \sup\{c_2(x, y) : |x - y| \ge r\}. $$

A typical example is the Riesz class of costs where $C = 2^{-s}$. In this case, if we assume that $\rho(X) > 1$ and it is not too concentrate, that is $\kappa := \sup_{x \in X} \rho(B_r(x)) < 1$ (for instance $\rho$ has no atom). The support of the minimiser is characterized as follows

$$ \text{Supp}(\gamma^*) \subset \left[0, 1 + \rho(X) \max \left( \frac{2}{C^2}, \frac{M(r)}{(1 - \kappa)m(2R_0)} \right) \right] \quad (2.16) $$

where $R_0$ is the smallest radius such that $\rho(X \setminus B_{R_0}) \le 1/2$ and $R_0 \ge R$. Notice that the Coulomb cost is indeed asymptotically doublig but the support estimate (2.16) is actually less precise then in our previous result.

* **The 3D Coulomb cost:** the case of the Coulomb cost deserves a special attention since it plays a central role in Density Functional Theory and the approximation of Lévy-Lieb/Lieb functionals. First, notice that the bound of theorem [DLN22][3.3] gives

$$ \text{Supp}(\gamma^*) \subset \left[ \frac{N_0 + 1}{3}, 3N_0 - 3 \right], \quad \rho(X) = N_0 \in \mathbb{N} \setminus \{1\}. $$

A natural question is then to ask whether the length of the support is really of the order $N_0$ or smaller. We remind that the ionization conjecture is a celebrated problem in quantum mathematical physics [Sim00; LS10] which states that a nucleus with charge $Z$ can bind at most $m \le Z + M$ electrons for some $M$. A natural similar conjecture in our case would be that the support of any minimizers for $\text{GCOT}^c(\rho)$ is included in $[\rho(\mathbb{R}^d) - M, \rho(\mathbb{R}^d) + M]$ for a universal constant $M$. This happens to be true in dimension $d = 1$ with $M = 1$, but wrong in dimension $d \ge 2$. We can however prove that the length of the support is much smaller that the average $\rho(\mathbb{R}^d)$, at most of the order of $\sqrt{\rho(\mathbb{R}^d)}$ in any dimension. In particular [DLN22][Thm 3.5] state that for $X = \mathbb{R}^d$ and $c_2(x, y) = |x - y|^{-1}$ for any non-negative measures such that the grand canonical cost is finite and $\rho(\mathbb{R}^d) > 1$ any minimizers satisfies

$$ \text{supp}(\gamma^*) \subset \left[ \lfloor \rho(X) \rfloor - \frac{1}{2} \sqrt{8 \lfloor \rho(X) \rfloor + 9} + \frac{3}{2}, \lfloor \rho(X) \rfloor + \frac{1}{2} \sqrt{8 \lfloor \rho(X) \rfloor - 7} - \frac{1}{2} \right]. \quad (2.17) $$
Moreover, when $\rho(\mathbb{R}^d) = 2$ then $\text{Supp}(\gamma^*) = \{2\}$. Let us remark that this result has been obtained by exploiting a method introduced in [FNV18; FKN16] and an extension of classical c-monotonicity to the grand canonical case, see [DLN22][Lemma 3.1]. We end the discussion on the high dimensional case with an example which has support of the length of the order $\rho(X)^\alpha$ with $\alpha \sim 0.38$. Let us consider 6 points $x_1, \dots, x_6$ in the plane $\mathbb{R}^2$ and the uniform measure $\rho = \frac{1}{2} \sum_{i=1}^6 \delta_{x_i}$ such that $\text{GCOT}^c(\rho)$ admits

$$ \gamma_2^* = \frac{1}{2} \delta_{x_1} \otimes_s \delta_{x_2}, \qquad \gamma_4^* = \frac{1}{2} \delta_{x_3} \otimes_s \delta_{x_4} \otimes_s \delta_{x_5} \otimes_s \delta_{x_6} \qquad (2.18) $$

as unique minimizer (in [DLN22][section 3.4.2] we explain in detail how this is possible). Then, form these points it is possible to construct inductively a sequence $(y_j^{(k)})_{j=1}^{6k} \subset \mathbb{R}^d$ so that,

$$ \rho^{(k)} = \frac{1}{2} \sum_{j=1}^{6k} \delta_{y_j^{(k)}}, \qquad \rho^{(k)}(\mathbb{R}^d) = \frac{6k}{2}, $$

the grand-canonical problem $\text{GCOT}_0^c(\rho^{(k)})$ admits a unique minimizer $\gamma^{(k)}$, which satisfies

$$ \text{Supp}(\gamma^{(k)}) = \left\{ \frac{6k - 2^k}{2}, \frac{6k + 2^k}{2} \right\}. $$

Notice that this example disproves the ionization-kind conjecture for high dimensional Coulomb cost.

* **The 1D problem:** as in the classical case of (multi-marginal) optimal transport the one dimensional case happens to be more tractable by exploiting the structure of $\mathbb{R}$. In this case we are able to extend the results of [CDD15] to the grand canonical framework, confirming in particular the shape of the optimal plans considered [MSG13]. We prove, see [DLN22][Thm 5.1], that for $m < \rho(\mathbb{R}) < m + 1$ we have that $\text{Supp}(\gamma) = \{m, m + 1\}$, whereas for $\rho(\mathbb{R}) \in \mathbb{N}$ the grand-canonical optimal solution is actually the canonical one. In addition, we show that the points are strictly correlated on the support of the optimal plan, that is, is given by a Monge state. When the mass of $\rho$ is not an integer, only one point is removed in a some region of the support. The running hypothesis is that $c = (c_n)_{n \ge 0}$ is a pairwise cost with $c_2(x, y)$ satisfying

$$ c_2(x, y) = w(|x - y|), \quad w: \mathbb{R} \to [0, +\infty] \text{ convex and decreasing.} \qquad (2.19) $$

Let us just clarify a little bit what a Monge state means in this case. For any atomless $\rho \in \mathcal{M}(\mathbb{R})$ let us consider the unique $m \in \mathbb{N}$ and $\eta \in [0, 1)$ such that $\rho(\mathbb{R}) = m + \eta$. We then split $\mathbb{R}$ into consecutive intervals where $\rho$ has the alternating masses $\eta$ and $1 - \eta$: $x_0 = -\infty \le x_1 \le x_2 \le \dots \le x_{2m} \le x_{2m+1} = +\infty$ where

$$ \rho((x_{2i}, x_{2i+1})) = \eta, \qquad \rho((x_{2i+1}, x_{2i+2})) = 1 - \eta, \quad \forall i = 0, \dots, m-1. $$

Of course, we automatically obtain $\rho((x_{2n}, x_{2n+1})) = \eta$ too. We now define $T: (x_0, x_{2m-1}) \to (x_2, x_{2m+1})$ as the $\rho$-a.e. unique increasing function such that $\rho((x, T(x))) = 1$. We have $T((x_i, x_{i+1})) \subseteq (x_{i+2}, x_{i+3})$. Iteratively we then define $T_1(x) = T(x)$ and $T_{i+1}(x) = T(T_i(x))$. Next we define $\gamma^*$ by

$$ \gamma_i^* = \begin{cases} 0 & \text{for } i \notin \{m, m+1\}, \\ \text{Sym.}[(\text{id}, T_1, T_2, \dots, T_{n-1})\sharp\rho|_{(x_1, x_2)}] & \text{for } i = m, \\ \text{Sym.}[(\text{id}, T_1, T_2, \dots, T_n)\sharp\rho|_{(x_0, x_1)}] & \text{for } i = m+1. \end{cases} \qquad (2.20) $$
where Sym. means that we are symmetrizing the plan. We have $\gamma^* \in \Pi_{GC}(\rho)$. If $\rho(\mathbb{R}) \in \mathbb{N}$, then $\gamma^*$ is exactly the canonical optimal plan considered in [CDD15]. Note that $\gamma^*$ does not depend on the cost, it is only a function of the density $\rho$.

### Duality

We conclude this section by studying the dual problem to (2.7). We know that the variables dual to the density $\rho$ are one-agent costs of the form $\Phi = (\Phi_n)_{n \ge 0}$ with $\Phi_0 = 0$ and $\Phi_n = \sum_{j=1}^n \varphi(x_j)$ for a given $\varphi \in C_b^0(\mathbb{R}^d)$. We should however not forget the other constraint that $\gamma$ forms a probability, $\gamma_0 + \sum_{n \ge 1} \gamma_n(X^n) = 1$, which requires the introduction of an additional Lagrange multiplier $\beta$. This constraint is independent of the density constraint, on the contrary to the usual multi-marginal problem. This leads us to the following dual problem

$$ \mathcal{D}(\rho) := \sup \left\{ \beta + \int_X \varphi(x) \, d\rho(x) : \beta \le c_0, \quad \varphi \in C_b(\mathbb{R}^d), \right. \\ \left. \beta + \sum_{j=1}^n \varphi(x_j) \le c_n(x_1, \dots, x_n), \quad \forall n \ge 1 \right\}. \quad (2.21) $$

If we take any $\beta$ and $\varphi \in C_b(\mathbb{R}^d)$ satisfying the above constraints and any $\gamma \in \Pi_{GC}(\rho)$, then we have

$$ \langle \gamma, c \rangle := c_0\gamma_0 + \sum_{n \ge 1} \int_{X^n} c_n \, d\gamma_n \ge \beta + \int_X \varphi(x) \, d\rho(x) $$

which proves that $\text{GCOT}_0^c(\rho) \ge \mathcal{D}(\rho)$ in all situations. We would like to have equality. It is possible to rewrite the dual problem in a slightly different manner. Let us introduce the grand-canonical ground state energy in the potential $\varphi$

$$ E(\varphi) := \inf_{\substack{n \ge 0 \\ x_1, \dots, x_n \in X^n}} \left( c_n(x_1, \dots, x_n) - \sum_{j=1}^n \varphi(x_j) \right) = \inf_{\gamma} \langle \gamma, c - \Phi \rangle. $$

The last infimum taken over all grand-canonical probabilities $\gamma$ so that $\langle \gamma, c \rangle < \infty$ and $\rho_\gamma(X) < \infty$ (without any other constraint on $\rho_\gamma$). Like in Density Functional Theory [Lie83; LLS23], we rewrite the infimum over $\gamma$ as an infimum over $\rho$ and then an infimum over all $\gamma$ having this density $\rho$:

$$ E(\varphi) := \inf_{\rho} \left\{ \text{GCOT}_{0}^{c}(\rho) - \int_{X} \varphi(x) \, d\rho(x) \right\} $$

with the infimum taken over all finite non-negative measures. We see that $E$ is nothing but the Legendre-Fenchel transform of $\text{GCOT}_0^c$. On the other hand, in (2.21) the largest possible $\beta$ at fixed $\varphi$ is indeed equal to $E(\varphi)$, and therefore we can rewrite (2.21) in the form

$$ \mathcal{D}(\rho) = \sup_{\varphi \in C_b(\mathbb{R}^d)} \left\{ \int_X \varphi(x) \, d\rho(x) + E(\varphi) \right\}. \quad (2.22) $$

Thus, $\mathcal{D}$ is the Legendre-Fenchel transform of $E$. From Theorem 2.2.4 and the Fenchel duality theorem for convex lower semi-continuous functions [Sim11], we conclude the following, which is an extension of a well known result in the multi-marginal case [Kel].

**Theorem 2.2.8 (Duality).** Let $X \subset \mathbb{R}^d$ be any Borel set. Let $c = (c_n)_{n \ge 0}$ be a superstable family of lower semi-continuous costs. Then we have $\text{GCOT}_0^c(\rho) = \mathcal{D}(\rho)$.