Wrong solution with the default solver

[stla]

Hello,

I have a function which runs psolve and the user can choose the solver. With the default solver, the result is not correct.

> kantorovich_CVX(mu, nu) # wrong
[1] 0.1573778
> kantorovich_CVX(mu, nu, solver = "ECOS") # correct
[1] 0.158657

I will post a reprex later.

> xfun::session_info("CVXR")
R version 3.6.3 (2020-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS, RStudio 1.3.1093

Locale:
  LC_CTYPE=fr_BE.UTF-8       LC_NUMERIC=C               LC_TIME=fr_BE.UTF-8       
  LC_COLLATE=fr_BE.UTF-8     LC_MONETARY=fr_BE.UTF-8    LC_MESSAGES=fr_BE.UTF-8   
  LC_PAPER=fr_BE.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
  LC_TELEPHONE=C             LC_MEASUREMENT=fr_BE.UTF-8 LC_IDENTIFICATION=C       

Package version:
  bit_4.0.4           bit64_4.0.5         CVXR_1.0            ECOSolveR_0.5.3    
  gmp_0.6.0           graphics_3.6.3      grDevices_3.6.3     grid_3.6.3         
  lattice_0.20.41     Matrix_1.2.18       methods_3.6.3       osqp_0.6.0.3       
  R6_2.4.1            Rcpp_1.0.0          RcppEigen_0.3.3.7.0 Rmpfr_0.8.1        
  scs_1.3.2           stats_3.6.3         utils_3.6.3       

[stla]

And with ECOS or the default solver, I get a wrong solution of the problem. Only GLPK is correct.

[stla]

Here is the reprex.

kantorovich_CVX <- function(
  mu, nu, dist=NULL, solution=FALSE, stop_if_fail=TRUE, solver = "ECOS", ...
){
  m <- length(mu)
  n <- length(nu)
  # checks
  if(m != n) stop("mu and nu do not have the same length")
  if(!is.null(dist)){
    if(!is(dist, "matrix") || mode(dist) != "numeric")
      stop("dist must be a numeric matrix")
    if(nrow(dist)!=m || ncol(dist)!=m)
      stop("invalid dimensions of the dist matrix")
  }
  if(sum(mu)!=1 || sum(nu)!=1 || any(mu<0) || any(nu<0)){
    message("Warning: mu and/or nu are not probability measures")
  }
  #
  if(is.null(dist)) dist <- 1-diag(m)

  obj <- c(t(dist))
  A <- rbind(t(model.matrix(~0+gl(m,n)))[,],
             t(model.matrix(~0+factor(rep(1:n,m))))[,])
  x <- Variable(m*n)
  objective   <- Minimize(t(obj) %*% x)
  constraints <- list(x >= 0, A%*%x == c(mu,nu))
  problem     <- Problem(objective, constraints)

  kanto <- psolve(problem, solver = solver, ...)
  # status
  if(kanto$status != "optimal"){
    if(stop_if_fail){
      stop(sprintf("No optimal solution found: status %s \n", kanto$status))
    }else{
      warning(sprintf("No optimal solution found: status %s \n", kanto$status))
      return(kanto)
    }
  }
  # output
  out <- kanto$value
  if(solution) attr(out, "solution") <-
    matrix(kanto$getValue(x), nrow=m, byrow=TRUE)
  out
}

mu <- c(1/4, 3/4, 0, 0)
nu <- c(0, 1/2, 1/2, 0)
dist <- structure(c(0, 1/3, 2/3, 1, 1/3, 0, 1/3, 2/3,
                       2/3, 1/3, 0, 1/3, 1, 2/3, 1/3, 0), .Dim = c(4L, 4L))
kanto <- kantorovich_CVX(mu, nu, dist = dist, solution = TRUE, solver = NA) # default solver
kanto_ECOS <- kantorovich_CVX(mu, nu, dist = dist, solution = TRUE, solver = "ECOS") 
kanto_GLPK <- kantorovich_CVX(mu, nu, dist = dist, solution = TRUE, solver = "GLPK") 

> kanto
[1] 0.249998
attr(,"solution")
              [,1]          [,2]         [,3]          [,4]
[1,]  1.681841e-06  1.067960e-01 1.432006e-01  1.682059e-06
[2,] -1.675999e-06  3.932054e-01 3.567980e-01 -1.676071e-06
[3,] -2.919593e-09 -6.841739e-07 6.900872e-07 -2.992153e-09
[4,] -2.921088e-09 -6.841754e-07 6.900917e-07 -2.993648e-09
> kanto_ECOS
[1] 0.25
attr(,"solution")
              [,1]          [,2]         [,3]          [,4]
[1,]  8.329396e-12  1.219023e-01 1.280977e-01  1.913124e-12
[2,] -1.679298e-12  3.780977e-01 3.719023e-01  1.946880e-12
[3,] -3.018023e-12  1.153635e-12 4.083524e-12 -2.217068e-12
[4,] -3.618322e-12 -1.842714e-13 5.443933e-12 -1.631836e-12
> kanto_GLPK
[1] 0.25
attr(,"solution") # correct
     [,1] [,2] [,3] [,4]
[1,]    0 0.25  0.0    0
[2,]    0 0.25  0.5    0
[3,]    0 0.00  0.0    0
[4,]    0 0.00  0.0    0

[stla]

Hmm... in fact no, the solution is correct. But there are multiple solutions. Looks like GLPK searches among the extremal solutions.

I'm understand the impression that warm_start does not work as expected: verbose = TRUE shows that it is "on" when you set this option to FALSE.