# Fairness and discrimination, PhD Course, #4 Wasserstein Distances and Optimal Transport

For the fourth course, we will discuss Wasserstein distance and Optimal Transport. Last week, we mentioned distances, dissimilarity and divergences. But before talking about Wasserstein, we should mention Cramer distance.

#### Cramer and Wasserstein distances

The definition of Cramér distance, for $k\geq1$, is

while Wasserstein will be (also for $k\geq1$)

If we consider cumulative distribution functions, for the first one (Cramer), we consider some sort of “vertical” distance, while for the second one (Wasserstein), we consider some “horizontal” one,

Obviously, when $k=1$, the two distances are identical

c1 = function(x) abs(pnorm(x,0,1)-pnorm(x,1,2)) w1 = function(x) abs(qnorm(x,0,1)-qnorm(x,1,2)) integrate(c1,-Inf,Inf)$value [1] 1.166631 integrate(w1,0,1)$value [1] 1.166636

But when $k>1$, it is no longer the case.

c2 = function(x) (pnorm(x,0,1)-pnorm(x,1,2))^2 w2 = function(u) (qnorm(u,0,1)-qnorm(u,1,2))^2 sqrt(integrate(c2,-Inf,Inf)$value) [1] 0.5167714 sqrt(integrate(w2,0,1)$value) [1] 1.414214

For instance, we can illustrate with a simple multinomial distribution, and the distance with some Binomial one, with some parametric inference based on distance minimization $$\theta^\star=\text{argmin}\{d(p,q_{\theta})\}$$(where here a multinomial distribution with parameters $\boldsymbol{p}=(.5,.1,.4)$, taking values respectively in $\{0,1,10\}$, while the binomial distribution has probabilities $\boldsymbol{q}_{\theta}=(1-\theta,\theta)$, taking values in $\{0,10\}$)

One can prove that

while

When $k=1$, observe that the distance is easy to compute when distributions are ordered

When $k=2$, the two distances are not equal

In the Gaussian (and the Bernoulli) case, we can get an expression for the distance (and much more as we will see later on)

There are several representations for $W_2$

And finally, we can also discuss $W_{\infty}$

#### Wasserstein distances, and optimal transport

Wasserstein distance can also we written using some sort of expected values, when considering random variables instead of distributions, and some best-case scenario, or cheapest transportation cost,

which lead to the so call Kantorovich problem

An alternative way to look at this problem is to consider a transport map, and a push-forward measure

This is simply

Of course such mapping exist

We can then consider Monge problem

And interestingly, those two problems are (somehow) equivalent

#### Discrete case

If $\boldsymbol{a}_{{A}}\in\mathbb{R}_+^{\color{red}{n_{{A}}}}$ and $\boldsymbol{a}_{{B}}\in\mathbb{R}_+^{\color{blue}{n_{{B}}}}$, define$$U(\boldsymbol{a}_{{A}},\boldsymbol{a}_{{B}})=\big\lbrace M\in\mathbb{R}_+^{\color{red}{n_{{A}}}\times\color{blue}{n_{{B}}}}:M\boldsymbol{1}_{\color{blue}{n_{{B}}}}=\boldsymbol{a}_{A}\text{ and }{M}^\top\boldsymbol{1}_{\color{red}{n_{{A}}}}=\boldsymbol{a}_{B}\big\rbrace$$For convenience, let $U_{\color{red}{n_{{A}}},\color{blue}{n_{{B}}}}$ denote $$\displaystyle{U\left(\boldsymbol{1}_{n_{{A}}},\frac{\color{red}{n_{{A}}}}{\color{blue}{n_{{B}}}}\boldsymbol{1}_{n_{{B}}}\right)}$$ (so that $U_{\color{red}{n},\color{blue}{n}}$ is the set of permutation matrices associated with $\mathcal{S}_n$). Let $C_{i,j}=d(x_i,y_{j})^k$so that $$W_k^k(\boldsymbol{x},\boldsymbol{y}) = \underset{P\in U_{\color{red}{n_{{A}}},\color{blue}{n_{{B}}}}}{\text{argmin}} \Big\lbrace \langle P,C\rangle \Big\rbrace$$where$$\langle P,C\rangle = \sum_{i=1}^{\color{red}{n_{{A}}}} \sum_{j=1}^{\color{blue}{n_{{B}}}} P_{i,j}C_{i,j}$$ then consider $$P^* \in \underset{P\in U_{\color{red}{n_A},\color{blue}{n_B}}}{\text{argmin}} \Big\lbrace \langle P,C\rangle \Big\rbrace$$For the slides, if we have the same sample sizes in the two groups, we have

we can illustrate below (with costs, or distances)

And with different group sizes,

i.e., if we consider real datasets

And as usual, we can consider some penalized version. Define $$\mathcal{E}(P) = -\sum_{i=1}^{\color{red}{n_{{A}}}} \sum_{j=1}^{\color{blue}{n_{{B}}}} P_{i,j}\log P_{i,j}$$or$$\mathcal{E}'(P) = -\sum_{i=1}^{\color{red}{n_{{A}}}} \sum_{j=1}^{\color{blue}{n_{{B}}}} P_{i,j}\big[\log P_{i,j}-1\big]$$ or $$\mathcal{E}'(P) = -\sum_{i=1}^{\color{red}{n_{{A}}}} \sum_{j=1}^{\color{blue}{n_{{B}}}} P_{i,j}\big[\log P_{i,j}-1\big]$$ Define $$P^*_\gamma = \underset{P\in U_{\color{red}{n_{{A}}},\color{blue}{n_{{B}}}}}{\text{argmin}} \Big\lbrace \langle P,C\rangle -\gamma \mathcal{E}(P) \Big\rbrace$$ The problem is strictly convex.

#### Sinkhorn relaxation

This idea is related to the following theorem

Consider a simple optimal transportation problem between 6 points to 6 other points,

set.seed(123) x = (1:6)/7 y = runif(9) x [1] 0.14 0.29 0.43 0.57 0.71 0.86 y[1:6] [1] 0.29 0.79 0.41 0.88 0.94 0.05 library(T4transport) Wxy = wasserstein(x,y[1:6]) Wxy$plan that we can visualize below (the first observation of $\boldsymbol{x}$ is matched with the last one of $\boldsymbol{y}$, the second observation of $\boldsymbol{x}$ is matched with the first one of $\boldsymbol{y}$, etc) We observe that we simply match according to ranks. But we can also use a penalized version Sxy = sinkhorn(x, y[1:6], p = 2, lambda = 0.001) Sxy$plan

here with a very small pernalty

or a larger one

Sxy = sinkhorn(x, y[1:6], p = 2, lambda = 0.05) Sxy$plan In the discrete case, optimal transport can be related to Hardy-Littlewood-Polya inequality, that is related to the idea of matching based on ranks (corresponding to a monotone mapping function) We have then In the bivariate dicrete case, we have #### Optimal mapping We have mentioned that, in the univariate setting and clearly, $\mathcal{T}^\star$ is increasing. In the Gaussian case, for example$$x_{{B}}=\mathcal{T}^\star(x_{{A}})= \mu_{{B}}+\sigma_{{B}}\sigma_{{A}}^{-1} (x_A-\mu_{{A}}).$$In the multivariate case, we need a more general concept of increasingness to define an “increasing” mapping $\mathcal{T}^\star:\mathbb{R}^k\to\mathbb{R}^k$. In the Gaussian case, for example, we have a linear mapping,$$\boldsymbol{x}_{{B}} = \mathcal{T}^\star(\boldsymbol{x}_{{A}})=\boldsymbol{\mu}_{{B}} + \boldsymbol{A}(\boldsymbol{x}_{{A}}-\boldsymbol{\mu}_{{A}})$$where $\boldsymbol{A}$ is a symmetric positive matrix that satisfies $$\boldsymbol{A}\boldsymbol{\Sigma}_{{A}}\boldsymbol{A}=\boldsymbol{\Sigma}_{{B}},$$ which has a unique solution given by $$\boldsymbol{A}=\boldsymbol{\Sigma}_{{A}}^{-1/2}\big(\boldsymbol{\Sigma}_{{A}}^{1/2}\boldsymbol{\Sigma}_{{B}}\boldsymbol{\Sigma}_{{A}}^{1/2}\big)^{1/2}\boldsymbol{\Sigma}_{{A}}^{-1/2},$$ where $\boldsymbol{M}^{1/2}$ is the square root of the square (symmetric) positive matrix $\boldsymbol{M}$ based on the Schur decomposition ($\boldsymbol{M}^{1/2}$ is a positive symmetric matrix). In R, for example, use the expm package, M = expm::sqrtm(matrix(c(1,1.2,1.2,2),2,2)) M [,1] [,2] [1,] 0.8244771 0.5658953 [2,] 0.5658953 1.2960565 M %*% M [,1] [,2] [1,] 1.0 1.2 [2,] 1.2 2.0 #### Optimal mapping, on real data To illustrate, it is possible to consider the optimal matching, between the height of $n$ men and $n$ women, Another example (discussed in Optimal Transport for Counterfactual Estimation: A Method for Causal Inference – with a nice R notebook created by Ewen), consider Black and non-Black mothers in the U.S. or the joint mapping, in dimension 2 We will spend more time on those functions (and the related concept) in a few weeks, when discussing barycenters and geodesics… More details in the slides (online) and in the forthcoming textbook, # La normalité ? avant tout une histoire de taille Un billet rapide pour préciser le fond de ma pensée lorsque j’explique aux étudiants (qui travaillent sur leur devoir) que sur les bases de données avec beaucoup d’observations, il est très dur d’accepter l’hypothèse de normalité. Pour ça, simulons des jeux de données: on a une seule variable explicative, un bruit gaussien, homoscédastique. Mais on va mal spécifier le modèle. En particulier, le vrai modèle ne sera pas linéaire. Voilà ce qu’on simule, et que l’on va chercher à ajuster. n=50 X1=runif(n) epsilon=rnorm(n) Y=-2+8*X1-4*(X1<.2)+epsilon base=data.frame(Y,X1) reg=lm(Y~X1,data=base) Autrement dit, les résidus (du modèle linéaire) ne devraient pas être gaussiens, car on n’ajuste pas le bon modèle. Maintenant, mon point est le suivant: avec peu de données, la plupart des tests de normalité ne vont pas rejeter l’hypothèse de normalité. Alors qu’avec beaucoup de données, l’hypothèse de normalité sera presque toujours rejetée. Prenons par exemple le test Lillie, évoqué dans un vieux billet, qui est une adaptation du test de Kolmogorov-Smirnov quand on compare la distribution empirique avec la loi normale dont les paramètres ont été obtenues par maximum de vraisemblance. Sur le graphique ci-dessous, à gauche, on a en abscisse la taille de l’échantillon (qui va de 10 à 1000), et en ordonnée la p-value du test Lillie. Ou plus précisément, la courbe bleue est la région où s’est trouvée la p-value dans 90% des cas (la courbe supérieure étant le quantile à 95% de la p-value sur 10,000 échantillons simulés). La courbe rouge est la valeur moyenne de la p-value. A droite, on a – toujours en fonction de la taille de l’échantillon – la proportion des échantillons simulés pour lesquels la p-value était supérieure à 5%. Autrement dit, plus la taille de l’échantillon est grande, plus on a de chance de rejeter (fort justement d’ailleurs) l’hypothèse de normalité. Et ce, uniquement sur un effet taille. Ou pour être plus rigoureux, sur les petits échantillons, on accepte beaucoup trop souvent (à tort cette fois) l’hypothèse de normalité. On peut aussi regarder le test de Shapiro-Wilk ou d’Anderson-Darling, ou encore de Cramér-von Mises Le code R pour faire ces tests est relativement simple, X1=runif(n) Y=-2+8*X1-4*(X1<.2)+rnorm(n) base=data.frame(Y,X1) reg=lm(Y~X1,data=base) AD[n,i]=ad.test(residuals(reg))$p.value
SW[n,i]=shapiro.test(residuals(reg))$p.value CV[n,i]=cvm.test(residuals(reg))$p.value
LL[n,i]=lillie.test(residuals(reg))\$p.value

Pour ceux qui veulent aller plus loin, se cache derrière la notion de puissance du test…