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