---
title: Using Julia to compute the Kantorovich distance
date : 2014-04-09
--- &lead
```{r setup0, echo=FALSE}
opts_chunk$set(fig.path="assets/fig/KantorovichWithJulia", tidy=FALSE)
```
In the article ['Using R to compute the Kantorovich distance'](http://stla.github.io/stlapblog/posts/KantorovichWithR.html) I have shown how to use the [cddlibb C library](http://web.mit.edu/sage/export/tmp/cddlib-094b/doc/cddlibman.pdf) through R with the help of the [rccd R package](http://cran.r-project.org/web/packages/rcdd/vignettes/vinny.pdf) to compute the Kantorovich distance between two probability measures (on a finite set). In the present article I show how to do so using three
different ways with [Julia](http://julialang.org/):
- **GLPK**: Similarly to the R way, this Julia way uses a C library, the [GLPK (GNU Linear Programming Kit)](http://www.gnu.org/software/glpk/) library, wrapped in a Julia package, named [GLPK](https://gplkjl.readthedocs.org/en/latest/glpk.html) too.
- **JuMP**: Using the [JuMP](https://github.com/JuliaOpt/JuMP.jl) Julia library, a [JuliaOpt](http://juliaopt.org/) project.
- **RationalSimplex**: Using Iain Dunning's [RationalSimplex](https://github.com/IainNZ/RationalSimplex.jl) module:
*Pure Julia implementation of the simplex algorithm for rational numbers*. This way is the only one allowing to get
the exact value when dealing with rational numbers.
In the current version of this article, I will only detail the `GLPK` method.
I express my gratitude to all
guys on [julia-users](https://groups.google.com/forum/#!forum/julia-users)
who kindly provided me their unvaluable help for this article.
# GLPK library
I will try to explain in details the approach using `GLPK`,
without assuming the reader has some knowledge about Julia.
## Setting the data of the problem
First, we define the probability measures $\mu$ and $\nu$ between which we want the Kantorovich distance to be computed.
```{r, eval=FALSE}
mu = [1/7, 2/7, 4/7]
nu = [1/4, 1/4, 1/2]
```
Recall that the Kantorovich distance is defined from an initial distance which we take to be the $0-1$ distance, and is represented in the $D$ matrix defined as follows with Julia:
```{r, eval=FALSE}
n = length(mu)
D = 1.-eye(n);
```
```{r, eval=FALSE}
julia> D
3x3 Array{Float64,2}:
0.0 1.0 1.0
1.0 0.0 1.0
1.0 1.0 0.0
```
Thus, the Julia `eye` function is similar to the R `diag` function.
Note that we have defined $D$ by typing "`1.-eye(n)`" and not "`1-eye(n)`" which
is ambiguous because "`1`" and "`eye(n)`" have not the same size...
I'm afraid you are puzzled because you don't know whether "`1.-eye(n)`" is
```{r, eval=FALSE}
1. - eye(n)
```
or
```{r, eval=FALSE}
1 .- eye(n)
```
? Don't worry, this is very easy to know with Julia, you can get the structure of "`1.-eye(n)`" as an
expression:
```{r, eval=FALSE}
julia> :(1.-eye(n))
:(.-(1,eye(n)))
```
That means the operator "`.-`" acts on the integer "`1`" and the matrix "`eye(n)`", whereas "`1. - eye(n)`" is the expression
for the operator "`-`" acting on the float "`1.`" and "`eye(n)`":
```{r, eval=FALSE}
julia> :(1. - eye(n))
:(-(1.0,eye(n)))
```
## Constraint matrix
The constraints matrix is
$$
A=\begin{pmatrix} M_1 \\ M_2 \end{pmatrix}
$$
where
$$ M_1 = \begin{pmatrix}
1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1
\end{pmatrix} $$
defines the linear equality constraints corresponding to $\mu$
and
$$ M_2 = \begin{pmatrix}
1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \\
0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 \\
0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1
\end{pmatrix}. $$
defines the linear equality constraints corresponding to $\nu$.
I define these matrices as follows in Julia:
```{r, eval=FALSE}
M1=zeros(n,n*n)
for i in 1:n
M1[i,(1:n)+n*(i-1)]=1
end
M2=repmat(eye(n),1,n)
A=vcat(M1,M2);
```
```{r, eval=FALSE}
julia> A
6x9 Array{Float64,2}:
1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 1.0 1.0 1.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 1.0
1.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0
0.0 1.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0
0.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 1.0
```
Recall that the constraints of our problem are the linear equality constraints
$$M_1 P = \begin{pmatrix} \mu(a_1) \\ \mu(a_2) \\ \mu(a_3) \end{pmatrix}
\qquad \text{and} \qquad
M_2 P = \begin{pmatrix} \nu(a_1) \\ \nu(a_2) \\ \nu(a_3) \end{pmatrix}$$
where $P$ is the vector formed by the variables $p_{ij} \geq 0$.
## GLPK in action
First load the package, initialize the problem, and give it a name:
```{r, eval=FALSE}
import GLPK
lp = GLPK.Prob()
GLPK.set_prob_name(lp, "kanto")
```
Computing the Kantorovich distance is a minimization problem, declared as follows:
```{r, eval=FALSE}
GLPK.set_obj_dir(lp, GLPK.MIN)
```
(`obj` refers to *objective function*, the function to be optimized).
Now we use the `GLPK.set_row_bnds` function to set the hand side vector (the *bounds*)
of the linear constraints and to specify the type of
our constraints. Here these are *equality* constraints and this type is specified by `GLPK.FX`:
```{r, eval=FALSE}
GLPK.add_rows(lp, size(A,1))
for i in 1:n
GLPK.set_row_bnds(lp, i, GLPK.FX, mu[i], mu[i])
GLPK.set_row_bnds(lp, i+n, GLPK.FX, nu[i], nu[i])
end
```
Now we specify the positivity constraints $p_{ij} \geq 0$ about the variables $p_{ij}$
corresponding to the columns of the constraint matrix, and we attach to each column the
corresponding coefficient of the objective function, given here by the matrix $D$:
```{r, eval=FALSE}
GLPK.add_cols(lp, size(A,2))
k=0
for i in 1:n, j in 1:n
k = k+1
GLPK.set_col_bnds(lp, k, GLPK.LO, 0.0, 0.0)
GLPK.set_obj_coef(lp, k, D[i,j])
end
```
We are ready ! Load the matrix, run the algorithm :
```{r, eval=FALSE}
GLPK.load_matrix(lp, sparse(A))
GLPK.simplex(lp);
```
and get the solution:
```{r, eval=FALSE}
julia> GLPK.get_obj_val(lp)
0.10714285714285715
```
As we have seen in the [previous article](http://stla.github.io/stlapblog/posts/KantorovichWithR.html), the exact Kantorovich distance between $\mu$ and $\nu$ is $\dfrac{3}{28}$:
```{r, eval=FALSE}
julia> 3/28
0.10714285714285714
```
Have you noticed the results are *not* exactly the same:
```{r, eval=FALSE}
julia> GLPK.get_obj_val(lp) - 3/28
1.3877787807814457e-17
```
Thus, the `GLPK.simplex` method does not achieve
the best approximation of $3/28$ in the 64 bit precision. A better
precision is achieved by the `GLPK.exact` function:
```{r, eval=FALSE}
GLPK.exact(lp);
```
```{r, eval=FALSE}
julia> GLPK.get_obj_val(lp) - 3/28
0.0
```
However, unfortunately, it is not possible to get the rational number $3/28$
with `GLPK`.
# JuMP library
The `JuMP` library allows a very convenient syntax. As compared to `GLPK`, no matrix gymnastics is needed.
Watch this concision:
```{r, eval=FALSE}
using JuMP
mu = [1/7, 2/7, 4/7]
nu = [1/4, 1/4, 1/2]
n = length(mu)
m = Model()
@defVar(m, p[1:n,1:n] >= 0)
@setObjective(m, Min, sum{p[i,j], i in 1:n, j in 1:n; i != j})
for k in 1:n
@addConstraint(m, sum(p[k,:]) == mu[k])
@addConstraint(m, sum(p[:,k]) == nu[k])
end
solve(m)
```
```{r, eval=FALSE}
julia> println("Optimal objective value is:", getObjectiveValue(m))
Optimal objective value is:0.10714285714285715
julia> 3/28
0.10714285714285714
```
As you can see, the best 64-bit precision approximation is not achieved.
But it is possible to get it.
`JuMP` uses a solver, and we have not specified any solver.
Then `JuMP` (actually `MathProgBase)` will search for an available solver and pick one by default.
We can use `GLPK.exact` by calling the `GLPKMathProgInterface` library and specifying the solver
in the `Model` function:
```{r, eval=FALSE}
using GLPKMathProgInterface
m = Model(solver=GLPKSolverLP(method=:Exact))
```
Then re-run the rest of the code, and you'll get:
```{r, eval=FALSE}
julia> getObjectiveValue(m)
0.10714285714285714
```
# RationalSimplex
The `RationalSimplex` module allows to get the exact Kantorovich distance as a rational number
as long as $\mu$ and $\nu$ have rational weights.
Rational numbers are specified in Julia with a double slash:
```{r, eval=FALSE}
mu= [1//7, 2//7, 4//7]
nu = [1//4, 1//4, 1//2]
```
We will not use matrix gymnastics to construct the constraint matrix $A$ with rational entries, we define it
in Julia with our bare hands below.
```{r, eval=FALSE}
include("myfolder/RationalSimplex.jl")
using RationalSimplex
using Base.Test
b = [mu, nu] # 'row bounds'
c = [0//1, 1//1, 1//1, 1//1, 0//1, 1//1, 1//1, 1//1, 0//1] # objective coefficients attached to the columns (D matrix in stacked form)
A = [1//1 1//1 1//1 0//1 0//1 0//1 0//1 0//1 0//1;
0//1 0//1 0//1 1//1 1//1 1//1 0//1 0//1 0//1;
0//1 0//1 0//1 0//1 0//1 0//1 1//1 1//1 1//1;
1//1 0//1 0//1 1//1 0//1 0//1 1//1 0//1 0//1;
0//1 1//1 0//1 0//1 1//1 0//1 0//1 1//1 0//1;
0//1 0//1 1//1 0//1 0//1 1//1 0//1 0//1 1//1]
x = status, simplex(c, :Min, A, b, ['=','=','=','=','=','=']);
```
The `simplex` function provides the solution of the linear programming problem, that is, the values of
$p_{ij}$ achieving the Kantorovich distance:
```{r, eval=FALSE}
julia> x
9-element Array{Rational{Int64},1}:
1//7
0//1
0//1
1//28
1//4
0//1
1//14
0//1
1//2
```
The Kantorovich distance is then obtained as the scalar product of `c` with the solution:
```{r, eval=FALSE}
julia> dot(c,x)
3//28
```