The kantorovich package

Stéphane Laurent

2024-04-23

The kantorovich package has two main features:

Exact results

With the help of the rccd and gmp packages, the kantorovich package can return the exact values of the extreme joinings and of the Kantorovich distance.

Quick example

As an example, take \(\mu\) and \(\nu\) the uniform probability measures on a finite set having three elements.

mu <- nu <- c(1/3, 1/3, 1/3)

The ejoinings function returns the extreme joinings of \(\mu\) and \(\nu\). In this case these are the \(6!\) permutation matrices:

library(kantorovich)
ejoinings(mu, nu)
## Message: You should enter mu and nu in rational with the gmp package.
## [[1]]
##           1         2         3
## 1 0.3333333 0.0000000 0.0000000
## 2 0.0000000 0.0000000 0.3333333
## 3 0.0000000 0.3333333 0.0000000
## 
## [[2]]
##           1         2         3
## 1 0.3333333 0.0000000 0.0000000
## 2 0.0000000 0.3333333 0.0000000
## 3 0.0000000 0.0000000 0.3333333
## 
## [[3]]
##           1         2         3
## 1 0.0000000 0.3333333 0.0000000
## 2 0.0000000 0.0000000 0.3333333
## 3 0.3333333 0.0000000 0.0000000
## 
## [[4]]
##           1         2         3
## 1 0.0000000 0.3333333 0.0000000
## 2 0.3333333 0.0000000 0.0000000
## 3 0.0000000 0.0000000 0.3333333
## 
## [[5]]
##           1         2         3
## 1 0.0000000 0.0000000 0.3333333
## 2 0.0000000 0.3333333 0.0000000
## 3 0.3333333 0.0000000 0.0000000
## 
## [[6]]
##           1         2         3
## 1 0.0000000 0.0000000 0.3333333
## 2 0.3333333 0.0000000 0.0000000
## 3 0.0000000 0.3333333 0.0000000

Since mu and nu were unnamed, the vector names c(1,2,3) has been automatically assigned to them. The Kantorovich distance between \(\mu\) and \(\nu\) is relative to a given distance on the state space of \(\mu\) and \(\nu\), represented by their vector names. By default, the kantorovich package takes the discrete \(0\mathrm{-}1\) distance. Obviously the Kantorovich distance is \(0\) on this example, because \(\mu=\nu\).

kantorovich(mu, nu)
## Message: You should enter mu and nu in rational with the gmp package.
## [1] 0

Note the message returned by both the ejoinings and the kantorovich functions. In order to get exact results, use rational numbers with the gmp package:

library(gmp)
mu <- nu <- as.bigq(c(1,1,1), c(3,3,3)) # shorter: as.bigq(c(1,1,1), 3) 
ejoinings(mu, nu)
## [[1]]
##   1     2     3    
## 1 "1/3" "0"   "0"  
## 2 "0"   "0"   "1/3"
## 3 "0"   "1/3" "0"  
## 
## [[2]]
##   1     2     3    
## 1 "1/3" "0"   "0"  
## 2 "0"   "1/3" "0"  
## 3 "0"   "0"   "1/3"
## 
## [[3]]
##   1     2     3    
## 1 "0"   "1/3" "0"  
## 2 "0"   "0"   "1/3"
## 3 "1/3" "0"   "0"  
## 
## [[4]]
##   1     2     3    
## 1 "0"   "1/3" "0"  
## 2 "1/3" "0"   "0"  
## 3 "0"   "0"   "1/3"
## 
## [[5]]
##   1     2     3    
## 1 "0"   "0"   "1/3"
## 2 "0"   "1/3" "0"  
## 3 "1/3" "0"   "0"  
## 
## [[6]]
##   1     2     3    
## 1 "0"   "0"   "1/3"
## 2 "1/3" "0"   "0"  
## 3 "0"   "1/3" "0"

User-specified distance

Let us try an example with a user-specified distance. Let’s say that the state space of \(\mu\) and \(\nu\) is \(\{a, b, c\}\), and then we use c("a","b","c") as the vector names.

mu <- as.bigq(c(1,2,4), 7)
nu <- as.bigq(c(3,1,5), 9)
names(mu) <- names(nu) <- c("a", "b", "c")

Define distance as a matrix

The distance can be specified as a matrix.

Assume the distance \(\rho\) is given by \(\rho(a,b)=1\), \(\rho(a,c)=2\) and \(\rho(b,c)=4\). The bigq matrices offered by the gmp package do not handle dimension names. But, in our example, the distance \(\rho\) takes only integer values, therefore one can use a numerical matrix:

M <- matrix(
  c(
    c(0, 1, 2),
    c(1, 0, 4),
    c(2, 4, 0)
  ), 
  byrow = TRUE, nrow = 3,
  dimnames = list(c("a","b","c"), c("a","b","c")))
kantorovich(mu, nu, dist=M)
## Big Rational ('bigq') :
## [1] 13/63

If the distance takes rational values, one can proceed as before with a character matrix:

M <- matrix(
  c(
    c("0", "3/13", "2/13"),
    c("1/13", "0", "4/13"),
    c("2/13", "4/13", "0")
  ), 
  byrow = TRUE, nrow = 3,
  dimnames = list(c("a","b","c"), c("a","b","c")))
kantorovich(mu, nu, dist=M)
## Big Rational ('bigq') :
## [1] 1/63

Define distance as a function

One can enter the distance as a function. In such an example, this does not sound convenient:

rho <- function(x,y){
  if(x==y) {
    return(0)
  } else {
    if(x=="a" && y=="b") return(1)
    if(x=="a" && y=="c") return(2)
    if(x=="b" && y=="c") return(4)
    return(rho(y,x))
  } 
}
kantorovich(mu, nu, dist=rho)
## Big Rational ('bigq') :
## [1] 13/63

Using a function could be more convenient in the case when the names are numbers:

names(mu) <- names(nu) <- 1:3

But one has to be aware that there are in character mode:

names(mu)
## [1] "1" "2" "3"

Thus, one can define a distance function as follows, for example with \(\rho(x,y)=\frac{|x-y|}{1+|x-y|}\):

rho <- function(x,y){
  x <- as.numeric(x); y <- as.numeric(y)
  return(as.bigq(abs(x-y), 1+abs(x-y)))
}
kantorovich(mu, nu, dist=rho)
## Big Rational ('bigq') :
## [1] 37/378

A non-square example

The kantorovich package also handles the case when mu and nu have different lengths, such as this example:

mu <- as.bigq(c(1,2,4), 7)
nu <- as.bigq(c(3,1), 4)
names(mu) <- c("a", "b", "c") 
names(nu) <- c("b", "c")
ejoinings(mu, nu)
## Caution: some names of mu and/or nu were missing or not compatible - automatic change
## [[1]]
##   b      c    
## a "1/7"  "0"  
## b "1/28" "1/4"
## c "4/7"  "0"  
## 
## [[2]]
##   b      c    
## a "1/7"  "0"  
## b "2/7"  "0"  
## c "9/28" "1/4"
## 
## [[3]]
##   b      c     
## a "0"    "1/7" 
## b "5/28" "3/28"
## c "4/7"  "0"   
## 
## [[4]]
##   b       c     
## a "0"     "1/7" 
## b "2/7"   "0"   
## c "13/28" "3/28"
kantorovich(mu, nu)
## Caution: some names of mu and/or nu were missing or not compatible - automatic change
## Big Rational ('bigq') :
## [1] 13/28

Note the caution message. The kantorovich package has to handle the fact that mu is given on the set \(\{a,b,c\}\) while nu is given on the set \(\{b,c\}\). It detects that the second set is included in the first one, and then treats nu on the bigger set \(\{a,b,c\}\) by assigning \(\nu(a)=0\). To avoid this caution message, the user has to enter this \(0\) value:

nu <- as.bigq(c(0,3,1), 4)
names(nu) <- c("a", "b", "c")

Other solvers

The kantorovich package provides three other functions to compute the Kantorovich distance:

Contrary to the kantorovich function, these two functions do not take care of the names of the two vectors mu and nu representing the two probability measures \(\mu\) and \(\nu\), and the distance to be minimized on average must be given as a matrix only, not a function.

A better precision is achieved by kantorovich_glpk. For instance, take the previous example for which we found \(13/63\) as the exact Kantorovich distance:

mu <- c(1,2,4)/7
nu <- c(3,1,5)/9
M <- matrix(
  c(
    c(0, 1, 2),
    c(1, 0, 4),
    c(2, 4, 0)
  ), 
  byrow = TRUE, nrow = 3)
kanto_lp <- kantorovich_lp(mu, nu, dist=M)
kanto_glpk <- kantorovich_glpk(mu, nu, dist=M)
kanto_CVX <- kantorovich_CVX(mu, nu, dist=M)

Then kantorovich_lp and kantorovich_CVX do not return the better decimal approximation of \(13/63\):

print(kanto_lp, digits=22)
## [1] 0.2063492063492062544849
print(kanto_glpk, digits=22)
## [1] 0.2063492063492063377517
print(kanto_CVX, digits=22)
## [1] 0.2063492063214846794494
print(13/63, digits=22)
## [1] 0.2063492063492063377517

But kantorovich_CVX is the fastest one, and it handles the case when the marginal probability measures mu and nu have a large support.