Tag Archives: distance

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))
[1] 1.166631
[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
[1] 0.5167714
[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


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}}}}, defineU(\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\rbraceFor 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})^kso 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\rbracewhere\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\rbraceFor 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,

x = (1:6)/7
y = runif(9)
[1] 0.14 0.29 0.43 0.57 0.71 0.86
[1] 0.29 0.79 0.41 0.88 0.94 0.05
Wxy = wasserstein(x,y[1:6])

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)

here with a very small pernalty

or a larger one

Sxy = sinkhorn(x, y[1:6], p = 2, lambda = 0.05)

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 examplex_{{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))
[,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,

D’un abus de langage…

Dans mon cours d’analyse des données, je commençais par présenter les diverses distances et normes, en particulier, celles qui dérivent d’un produit scalaire, comme la norme Euclidienne \|\boldsymbol{X}\|_{\mathbb{I}}^2=\boldsymbol{X}^\top\boldsymbol{X}, ou plus intéressante, la norme que l’on peut construire à partir d’une matrice \boldsymbol{M} symétrique définie positive, \|\boldsymbol{X}\|_{\boldsymbol{M}}^2=\boldsymbol{X}^\top\boldsymbol{M}\boldsymbol{X}. Un exemple bien connu étant la norme de Mahalanobis, ou la distance associée, où on prend pour \boldsymbol{M}  la matrice \boldsymbol{\Sigma}^{-1}\boldsymbol{\Sigma}=\text{Var}[\boldsymbol{X}]. Sauf que, par abus de langage, j’appelle aussi distance de Mahalanobis le distance obtenue en prenant pour matrice \boldsymbol{M} la matrice diagonale \boldsymbol{\Delta}_{\boldsymbol{\Sigma}}^{-1} définie par \boldsymbol{\Delta}_{\boldsymbol{\Sigma},ii}=\boldsymbol{\Sigma}_{ii}. En fait, c’est cette matrice que l’on utilise principale en analyse de données, par exemple en analyse des correspondances principales… Prenons un petit exemple, avec quelques données centrées

for(j in 1:ncol(X)) X[,j]=X[,j]-mean(X[,j])

et seulement 4 variables. Par défaut, quand on fait une analyse en composantes principales, on va réduire les données

acp = PCA(X,scale.unit = TRUE,graph = FALSE)
fviz_pca_var(acp, col.var = "cos2")

Le fait de centrer et réduire les données fait que l’on travaille avec la matrice de corrélation (\boldsymbol{X}^\top\boldsymbol{X}, à un facteur n, ou n-1, près). D’où le cercle de corrélation qui apparaît sur la figure. Et abusivement, quand je présente l’ACP, je dis qu’on a deux manière de voir : (1) utiliser le produit scalaire et la distance usuelle (Euclidienne) sur les données centrées et réduites, (2) utiliser la distance de Mahalanobis sur les données centrées. Mon abus de langage tient du fait qu’on utilise la version avec la matrice diagonale \boldsymbol{\Delta}_{\boldsymbol{\Sigma}}^{-1} , et pas la vraie distance de Mahalanobis… En discutant avec Gabriel après mon cours de ce matin, je me suis rendu compte de la confusion possible (même si je me faisais la réflexion lors de mes enregistrements de capsules, les jours passés). Mais en réalité, faire une analyse en composantes principales avec la vraie distance de Mahalanobis n’aurait aucun sens.

Pour s’en convaincre, utilisons les données centrées et réduites,

S = var(X)
Sd = diag(diag(S))
X1 = X%*%solve(chol(Sd))
A 1.0000000 0.10113165 -0.88021475 0.4733021
B 0.1011316 1.00000000 0.07885109 -0.7400999
C -0.8802147 0.07885109 1.00000000 -0.5743660
D 0.4733021 -0.74009993 -0.57436600 1.0000000


acp1 = PCA(X1,scale.unit = FALSE,graph = FALSE)
fviz_pca_var(acp1, col.var = "cos2")

On n’a plus ici, le cercle des corrélations car la fonction ne sait pas que mes données ont été centrées et réduite au préalable… mais on a la même sortie qu’auparavant à part ça…

Si maintenant on utilise la norme de Mahalanobis, on peut là aussi tranformer les données, en notant que \boldsymbol{X}^\top\boldsymbol{M}\boldsymbol{X} peut aussi s’écrire \boldsymbol{X}^\top\boldsymbol{M}^{1/2}\boldsymbol{M}^{1/2}\boldsymbol{X} soit encore la distance euclidienne standard \tilde{\boldsymbol{X}}^\top\tilde{\boldsymbol{X}} sur les données transformées \tilde{\boldsymbol{X}}={\boldsymbol{X}}{\boldsymbol{M}^{1/2}},

S = var(X)
X2 = X%*%solve(chol(S))
> var(X2)
              A             B             C             D
A  1.000000e+00 -3.582611e-17  5.307644e-16 -1.004784e-16
B -3.582611e-17  1.000000e+00 -3.021333e-16  4.164692e-17
C  5.307644e-16 -3.021333e-16  1.000000e+00 -3.649967e-16
D -1.004784e-16  4.164692e-17 -3.649967e-16  1.000000e+00

ce qui donne des données… orthogonales (pour ne pas dire indépendantes, on retrouve ici le processus de whitening, utilisé en ICA, analyse des composantes indépendantes. Et comme on a rendu les données indépendantes… l’analyse des composantes principales n’aura plus grand chose à raconter…

acp2 = PCA(X2,scale.unit = FALSE,graph = FALSE)
fviz_pca_var(acp2, col.var = "cos2")

Chaque axe explique ici pour 25% de l’inertie… bref, on ne va jamais faire une ACP avec la distance de Mahalanobis (la vraie). Elle ne va pas normaliser les variables, mais les orthonormaliser. Et pourtant, par abus de langage, je vais continuer à utiliser ce terme… Désolé…

INF7100, “rappels” de maths

En parallèle des trois parties de mon intervention sur la science des données, dans le cadre du cours INF7100, j’ai mis en ligne une série des capsules de “rappels” de maths.

  • 401: Concepts généraux (logarithme, exponentielle, dérivée, intégrale, approximation) pdf video (38:57)
  • 411: Optimisation (condition du premier ordre, descente de gradient, Lagrangien) pdf video (31:33)
  • 412: Ensembles (intersection, unions, fusion de bases, convexité, centroïdes) pdf video (10:13)
  • 421: Vecteurs et Matrices (orthogonalité, distance, norme, tenseurs) pdf video (43:20)
  • 422: Algèbre Linéaire (déterminant, valeurs propres, vecteurs propres, formes quadratiques, approximations) pdf
  • 423: Coordonnées (données circulaires, coordonnées polaires, projection, spatial) pdf
  • 431: Loi Gaussienne (loi normale, cas multivariée, limite et approximation, loi des grands nombres) pdf video (17:37)
  • 432: Lois Puissance (changement d’échelle, loi puissance, Zipf, Benford) pdf
  • 433: Loi de Durée (loi géométrique, loi exponentielle, loi de Poisson, loi des petits nombres) pdf
  • 434: Modèle de Markov (données séquentielles) pdf

Where hiding if you don’t want to get wet ?

Following the previous post, two additional remarks. Following a comment by@cosi, I have investigated quickly a binomial fit to the distribution of the number of people not getting wet, with a fixed number of players on the field. It looks like it should be a binomial distribution with a fixed probability (2/3) and with size parameter affine in the number of players. @guigui suggested some connexion with with “birds on a wire” problem (see e.g. http://www.cut-the-knot.org/)

for(i in 1:20){

for the implied size parameter above, and below the implied probability parameter.


(as functions of the number of players). I’d be glad to get more details on that 2/3 probability.

Now, let us investigate another question sent by email: “Where should you hide if you don’t want to get wet ?” A first idea could be the following: given that some players are already on the field, where should I go if I do not want to get wet ? Below are some simulations for 7 or 25 players (already on the field). The red area is the area so that I will become someone’s target (perhaps even the target of two players…). The green area is the safe zone.

(with 7 players above, and 25 below)

It looks like, on the border, it might be safer than in the middle of the field. But we have to confirm that intuition… or at least see if that intuition is valid.

Based on what was done the other day, it is possible to look where people that got wet were located (instead of counting dry players as done in the previous function). So here, we simply look where non wet players were standing

(d=as.matrix(dist(cbind(x,y), method = "euclidean",upper=TRUE)))
whonotwet=( (1:n) %notin% names(table(dmin)) )
for(i in 1:p){
for(j in 1:p){

based on function

"%notin%" <- function(x, y) x[!x %in% y]

On a given grid, we count people playing the game that ended dry (with might avoid boundary bias on nonparametric smooth estimator of distribution, as we’ll see later on). For instance with 11 players,

for(s in 1:100000){

Then we can plot the distribution, on the field,

COL=rev(heat.colors(101)); p=25
for(i in 1:p){
for(j in 1:p){

Red means a lot of non-wet people (i.e safer zones). Graphs below are with 7 and 11 players respectively (from the left to the right)

with the following distribution on the diagonal: corners are almost 4 times safer than the middle of the field, with 7 players,

Below are plotted distributions of locations of non-wet players when the total number of players was either 25 (on the left) and 101 (on the right)

with again on the diagonal

Hence, the border is rather safe, but next to the border, it is no safe any longer: is someone is standing right on the border, he will probably shoot at you: there is no one behind him ! This explains the stange behavior on the borders (and corners, thanks JP for the intuitive explanation).
But would it be completely different with a field shaped as a disk ?

using the previous technique of working on a fixed grid (or correcting for boundary bias, since the disk might cover only a fraction of the grid-square), or keeping coordinates of non-wet players, and using standard kernel-based estimator of the distribution (the light yellow circle outside the disk is simply due to bias of the kernel estimator on the border)

method = "euclidean",upper=TRUE)))
whonotwet=( (1:n) %notin% names(table(dmin)) )

for(s in 1:10000){

K=kde(x=M11, H=HP)


And note that the distribution of the number of players ending the game dry is the same, for a square field, or a disk,

method = "euclidean",upper=TRUE)))

method = "euclidean",upper=TRUE)))


But so far, it was still simple… I wonder what it might become if we consider a non-convex place, with walls, where player might hide…. Next time, a post on indoor paint-ball !

Vitesse et comportement dangereux au volant (partie 3)

Après plusieurs billets sur les accidents de la route (ici et ), parlons un peu de vitesse des véhicules. J’ai pu récupérer des données via l’inrets, d’observations près de Saclay en région parisienne, sur une 2×2 limitée à 110 km/h.


J’ai ainsi le passage de plusieurs millions de véhicules, sur quelques jours, incluant tous les véhicules. Mais si j’ai signé des papiers quant à la confidentialité des données, je dois pouvoir mettre des dessins construits à partir de ces données…

  • A quelle vitesse roulent les conducteurs parisiens ?

Comme pour les accidents, distinguons semaine et week end. En semaine, on retrouve les pics du matin et du soir, qui poussent les gens à ralentir. La courbe en noir est la vitesse moyenne des véhicules qui passent sur ce tronçon, et la courbe en bleu, la vitesse des 5% des véhicules les plus rapides,

Ces pics sont légèrement atténués le week end,

avec ici un trafic très dense tout l’après midi, entre 15 et 20 heures.

  • Quelle proportion de véhicules dépasse la limite autorisée ?

Comme pour les accidents, distinguons semaine et week end. En semaine, on retrouve les pics du matin et du soir, avec en gras la proportion des véhicules qui dépasse 110 km/h, et en trait fin, la proportion qui dépasse 115 km/h (soit 5% de plus),

avec des choses finalement assez proches le week end,

Autrement dit, aux heures de pointe, peu de monde dépasse la vitesse autorisée, mais vers midi, ou surtout la nuit, environ un véhicule sur quatre dépasse la vitesse autorisée (mais seulement un sur six la dépasse de plus de 5%).

  • Vitesse et respect des distances de sécurité

La base est très complète, puisque j’ai tous les véhicules qui sont passés près de Saclay, avec l’heure d’observation, la vitesse, ce qui permet d’obtenir la distance entre deux véhicules qui se suivent (en supposant que le premier véhicule maintienne sa vitesse constante).
Rappelons qu’il existe un lien théorique entre la vitesse d’un véhicule et la distance de freinage (expliqué ici ou ). Pour faire simple, la distance de freinage https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein02.png est la somme d’une distance parcourue avant de réagir https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein04.png, et d’une distance de freinage proprement dite https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein03.png,


Si on suppose que la personne met une seconde à réagir, la distance de réaction est simplement la vitesse https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein05.png(exprimée en m/h), i.e.


Quant au freinage à proprement parler, rappelons que les équations qui caractérisent le mouvement d’un objet sont de la forme


Aussi, en mettant la deuxième au carré, et en combinant avec la troisième, on a


avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein15.png (on cherche à arrêter le véhicule) et https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein14.png, aussi, avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/frein16.png, on obtient quelque chose qui peut s’écrire


à condition d’avoir des vitesses en km/h. Aussi en agglomération, il faut 23 m pour s’arrêter. Et sur notre route à Saclay, pour une personne qui roule à 110 km/h, il faut compter 110 m (ce qui correspond au message véhiculé par la sécurité routière).
Sur nos données, on peut représenter en rouge les individus qui ne respectent pas la distance de sécurité.

La courbe bleue est la distance moyenne des véhicule en fonction de la distance qui les sépare du véhicule qui les précède,
Aussi, en faisant une régression logistique (lissée), on obtient que 60% des conducteurs qui roulent entre 50 et 110 km/h ne respectent pas les distances de sécurité.