# Matching, Optimal Transport and Statistical Tests

To explain the “optimal transport” problem, we usually start with Gaspard Monge’s “Mémoire sur la théorie des déblais et des remblais“, where the the problem of transporting a given distribution of matter (a pile of sand for instance) into another (an excavation for instance). This problem is usually formulated using distributions, and we seek the “optimal” transport from one distribution to the other one. The formulation, in the context of distributions has been formulated in the 40’s by Leonid Kantorovich, e.g. from the distribution on the left to the distribution on the right.

Consider now the context of finite sets of points. We want to transport mass from points $\{A_1,\cdots,A_4\}$ to points $\{B_1,\cdots,B_4\}$. It is a complicated combinatorial problem. For 4 points, there are only 24 possible transfer to consider, but it exceeds 20 billions with 15 points (on each side). For instance, the following one is usually seen as inefficient

while the following is usually seen as much better

Of course, it depends on the cost of the transport, which depends on the distance between the origin and the destination. That cost is usually either linear or quadratic.

There are many application of optimal transport in economics, see eg Alfred’s book Optimal Transport Methods in Economics. And there are also applications in statistics, that what I’ve seen while I was discussing with Pierre while I was in Boston, in June. For instance if we want to test whether some sample were drawn from the same distribution,

set.seed(13)
npoints <- 25
mu1 <- c(1,1)
mu2 <- c(0,2)
Sigma1 <- diag(1, 2, 2)
Sigma2 <- diag(1, 2, 2)
Sigma2[2,1] <- Sigma2[1,2] <- -0.5
Sigma1 <- 0.4 * Sigma1
Sigma2 <- 0.4 *Sigma2
library(mnormt)
X1 <- rmnorm(npoints, mean = mu1, Sigma1)
X2 <- rmnorm(npoints, mean = mu2, Sigma2)
plot(X1[,1], X1[,2], ,col="blue")
points(X2[,1], X2[,2], col = "red")

Here we use a parametric model to generate our sample (as always), and we might think of a parametric test (testing whether mean and variance parameters of the two distributions are equal).

or we might prefer a nonparametric test. The idea Pierre mentioned was based on optimal transport. Consider some quadratic loss

ground_p <- 2
p <- 1
w1 <- rep(1/npoints, npoints)
w2 <- rep(1/npoints, npoints)
C <- cost_matrix_Lp(t(X1), t(X2), ground_p)
library(transport)
library(winference)
a <- transport(w1, w2, costm = C^p, method = "shortsimplex")

then it is possible to match points in the two samples

nonzero <- which(a\$mass != 0)
from_indices <- a\$from[nonzero]
to_indices <- a\$to[nonzero]
for (i in from_indices){
segments(X1[from_indices[i],1], X1[from_indices[i],2], X2[to_indices[i], 1], X2[to_indices[i],2])
}

Here we can observe two things. The total cost can be seen as rather large

> cost=function(a,X1,X2){
nonzero <- which(a\$mass != 0)
naa=a[nonzero,]
d=function(i) (X1[naa\$from[i],1]-X2[naa\$to[i],1])^2+(X1[naa\$from[i],2]-X2[naa\$to[i],2])^2
sum(Vectorize(d)(1:npoints))
}
> cost(a,X1,X2)
[1] 9.372472

and the angle of the transport direction is alway in the same direction (more or less)

> angle=function(a,X1,X2){
nonzero <- which(a\$mass != 0)
naa=a[nonzero,]
d=function(i) (X1[naa\$from[i],2]-X2[naa\$to[i],2])/(X1[naa\$from[i],1]-X2[naa\$to[i],1])
atan(Vectorize(d)(1:npoints))
}
> mean(angle(a,X1,X2))
[1] -0.3266797

> library(plotrix)
> ag=(angle(a,X1,X2)/pi)*180
> ag[ag<0]=ag[ag<0]+360
> dag=hist(ag,breaks=seq(0,361,by=1)-.5)
> polar.plot(dag\$counts,seq(0,360,by=1),main=”Test Polar Plot”,lwd=3,line.col=4)

(actually, the following plot has been obtain by generating a thousand of sample of size 25)

In order to have a decent test, we need to see what happens under the null assumption (when drawing samples from the same distribution), see

Here is the optimal matching

Here is the distribution of the total cost, when drawing a thousand samples,

VC=rep(NA,1000)
VA=rep(NA,1000*npoints)
for(s in 1:1000){
X1a <- rmnorm(npoints, mean = mu1, Sigma1)
X1b <- rmnorm(npoints, mean = mu1, Sigma2)
ground_p <- 2
p <- 1
w1 <- rep(1/npoints, npoints)
w2 <- rep(1/npoints, npoints)
C <- cost_matrix_Lp(t(X1a), t(X1b), ground_p)
ab <- transport(w1, w2, costm = C^p, method = "shortsimplex")
VC[s]=cout(ab,X1a,X1b)
VA[s*npoints-(0:(npoints-1))]=angle(ab,X1a,X1b)
}
plot(density(VC)

So our cost of 9 obtained initially was not that high. Observe that when drawing from the same distribution, there is now no pattern in the optimal transport

ag=(VA/pi)*180
ag[ag<0]=ag[ag<0]+360
dag=hist(ag,breaks=seq(0,361,by=1)-.5)
polar.plot(dag\$counts,seq(0,360,by=1),main="Test Polar Plot",lwd=3,line.col=4)

Nice isn’t it? I guess I will spend some time next year working on those transport algorithm, since we have great R packages, and hundreds of applications in economics…

# Risk Uncertainty and Decision Conference

Marc Henry will present our paper Local Utility and Multivariate Risk Aversion (writen with Alfred Galichon, still available online on http://papers.ssrn.com/at the Risk Uncertainty and Decision Conference in Evanston, IL, at Northwestern University.

I will try to upload the slides soon…

# Talk on multivariate comonotonicity and risk measures, JDS in Brussels

Today is the last day of the Journées de Statistique, in Brussels, http://jds2012.ulb.ac.be/. Alfred gave a survey on “Multivariate comonotonicity, stochastic orders and risk measures” in plenary session, this morning (as invited speaker). I have uploaded the slides.

# Local utility and multivariate risk aversion

Marc will give a talk today at the European Center for Advanced Research in Economics and Statistics (ECARES) today, at ULB in Brussels, based on some joint work with also Alfred (the paper can be found online on http://papers.ssrn.com/).

# Multivariate comonotonicity and multivariate dispersion orderings

Marc will give a talk at the Kyoto Symposium in 京都市 on November 29th, on “Multivariate comonotonicity and multivariate dispersion orderings” based on some joint work with also Alfred (the paper can be found online on http://papers.ssrn.com/).

# Dynamic and Multivariate Risk measures

Exposé dans le cadre du groupe de travail Dynamic and multivariate risk measures, organisé par Rose-Anna Dana (Univ. Paris Dauphine) et Alfred Galichon (Ecole Polytechnique). La conférence a lieu du jeudi 16 au vendredi 17, à l’Institut Henri Poincare.

Les orateurs sont Imen Bentahar, Arthur Charpentier, Rama Cont, Nicole El Karoui, Paul Embrechts, Damir Filipovic, Jean-Charles Rochet, Ludger Ruschendorf, Marco Scarsini, et Walter Schachermayer, et le programme est en ligne.