
Using Scala to compute the Kantorovich distance
20140609
SourceRecall that we previously provided an example of calculation of a Kantorovich distance using R and using Julia. Now it's time to do it with Scala and its breeze library, a ScalaNLP projet.
Our example is a linear programing problem given as follows. The unknown variables \(p_{ij}\) have a tabular structure \[P=\begin{pmatrix} p_{11} & p_{12} & p_{13} \\ p_{21} & p_{22} & p_{23} \\ p_{31} & p_{32} & p_{33} \end{pmatrix},\] are constrained by the following linear equality/inequality constraints: \[\begin{cases} {\rm (1a) } \quad \sum_j p_{ij} = \mu(a_i) & \forall i \\ {\rm (1b) } \quad \sum_i p_{ij} = \nu(a_j) & \forall j \\ {\rm (2) } \quad p_{ij} \geq 0 & \forall i,j \\ \end{cases}\] where \(\mu\) and \(\nu\) are two probability vectors, taken to be \[ \mu=\left(\frac{1}{7},\frac{1}{7},\frac{1}{7}\right) \qquad \text{and } \qquad \mu=\left(\frac{1}{4},\frac{1}{4},\frac{1}{2}\right) \] in our example, and the desideratum is to minimize the linear combination \[ \sum_{i,j} D_{ij}p_{ij}\] where \(D\) is a distance matrix, taken to be \[ D = \begin{pmatrix} 0 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \end{pmatrix}\] in our example.
Breeze expressions
The linear programming solver provided by
Breeze
is based on Apache's Simplex Solver. Most linear programming solvers are conceived to take as input the matrix of linear constraints, and this is quite annoying to construct this matrix in our situation: for the linear equality constraints this 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\). See how I defined these matrices using R and using Julia (with theGLPK
library), it was a bit tricky or a bit painful.As we have seen, the Julia
JuMP
package provides a very convenient way to write, and maybe more importantly, to read, the linear programing problem, using expressions to define the unknown variables \(p_{ij}\), the objective function and the constraints, which avoids in particular to construct the \(A\) matrix :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)
Using the Scala
breeze
library we also set the problem by using expressions. This is not as concise as the Julia JuMP library and this requires a bit of gymnastics with the Scala language, but essentially the spirit is the same: we write the problem, using expressions. Another Java library, WVLPsolver, allows to encode the problem with a syntax closer to the one ofJuMP
, but I have never tried it yet. Anyway there is no difficuly withbreeze
once you get a basic knowledge in the Scala language: mainly you only have to use a new type,Expression
, a member of thelp
object.Scala breeze in action
First, we define the probability measures \(\mu\) and \(\nu\) between which we want the Kantorovich distance to be computed, and the distance matrix \(D\):
var Mu = Array(1.0/7, 2.0/7, 4.0/7) var Nu = Array(1.0/4, 1.0/4, 1.0/2) val n = Mu.length val D = Array.fill(n,n)(0) for( i < 0 to n1 ) for( j < 0 to n1 ) if( !(i==j) ) D(i)(j) = 1
To show an array in Scala, do:
scala> println(D.deep.mkString("\n")) Array(0, 1, 1) Array(1, 0, 1) Array(1, 1, 0)
Before introducing the \(P\) matrix of variables, we need to create a new linear programing object:
val lp = new breeze.optimize.linear.LinearProgram()
and then we can define \(P\) as an array of
Real
variables:val P = Array.fill(n,n)(lp.Real())
Here is our matrix of unknown variables:
scala> println(P.deep.mkString("\n")) Array(x_0, x_1, x_2) Array(x_3, x_4, x_5) Array(x_6, x_7, x_8)
The following code defines the objective function as an expression:
val Objective = ( for( i < 0 to 2 ) yield ( for( (x, a) < P(i) zip D(i) ) yield (x * a) ).reduce(_ + _) ).reduce(_ + _)
scala> println(Objective) (x_0) * 0.0 + (x_1) * 1.0 + (x_2) * 1.0 + (x_3) * 1.0 + (x_4) * 0.0 + (x_5) * 1.0 + (x_6) * 1.0 + (x_7) * 1.0 + (x_8) * 0.0
Let us explain the above code. The following command returns the scalar product of two vectors
V1
andV2
:for( (x, a) < V1 zip V2 ) yield (x * a) ).reduce(_ + _)
Then our code firstly generates the scalar product (as an expression) of the \(i\)th line
P(i)
of \(P\) and the \(i\)th lineD(i)
of \(D\), for every \(i=0,1,2\) (indexation starts at \(0\) in Scala, not at \(1\)), and then it generates the sum of these three scalar products, which is the expression for our objective function.Now we define the constraints in a
Constraint
array. There are \(2n\) equality constraints and \(n^2\) positivity (inequality) constraints, we will store them in this array:var Constraints = new Array[lp.Constraint](2*n+n*n)
The equality constraints are about the sum of the variables in each row and each column of the \(P\) matrix. We set them as follows:
val tP=P.transpose for (i < 0 to n1) { Constraints(i) = P(i).reduce[lp.Expression](_ + _) =:= Mu(i) Constraints(i+n) = tP(i).reduce[lp.Expression](_ + _) =:= Nu(i) }
And finally we set the positivity constraints as follows:
for (i < 0 to n1) { var k = 2*n + n*i for (j < 0 to n1) Constraints(k+j) = P(i)(j) >=0 }
There's quite bit of gymnastics in the previous code... I'll possibly simplify it in a future version, as I'm at my first steps with Scala
breeze
.Since the library only allows to maximize the objective function, whereas we aim to minimize it, we take its opposite in the following
lpp
object which fully represents our linear programing problem :val lpp = ( Objective subjectTo( Constraints:_* ) ) val result = lp.maximize(lpp)
And then we get the Kantorovich distance:
scala> println(result.value) 0.1071428571428571
which almost achieves the best 64bit precision approximation of the exact value \(3/28\):
scala> println(3.0/28) 0.10714285714285714