In the article 'Using R to compute the Kantorovich distance' I have shown how to use the cddlibb C library through R with the help of the rccd R package 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:

• GLPK: Similarly to the R way, this Julia way uses a C library, the GLPK (GNU Linear Programming Kit) library, wrapped in a Julia package, named GLPK too.

• JuMP: Using the JuMP Julia library, a JuliaOpt project.

• RationalSimplex: Using Iain Dunning's RationalSimplex 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 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.

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:

n = length(mu)
D = 1.-eye(n);

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

1. - eye(n)

or

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:

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)":

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:

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);

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:

import GLPK 
lp = GLPK.Prob()
GLPK.set_prob_name(lp, "kanto")

Computing the Kantorovich distance is a minimization problem, declared as follows:

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:

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\):

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 :

GLPK.load_matrix(lp, sparse(A))
GLPK.simplex(lp);

and get the solution:

julia> GLPK.get_obj_val(lp)
0.10714285714285715

As we have seen in the previous article, the exact Kantorovich distance between \(\mu\) and \(\nu\) is \(\dfrac{3}{28}\):

julia> 3/28
0.10714285714285714

Have you noticed the results are not exactly the same:

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:

GLPK.exact(lp);

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:

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)

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:

using GLPKMathProgInterface
m = Model(solver=GLPKSolverLP(method=:Exact))

Then re-run the rest of the code, and you'll get:

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:

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.

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:

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:

julia> dot(c,x)
3//28