Fast Simplex in R

Question

I am not sure whether this is the right place to ask, but I need to efficiently implement a very large linear programming problem in R (and there is no way around R here). I have tried out some packages like lpSolve but the results seem unsatisfactory. I would be glad for any advice on packages, or alternatively a better place to ask this question.

Here is the problem:

N <- 10^4 # number of products
K <- 10^4 # number of scenarios

### Get expectation and covariance matrix

mu <- rep(100,N)
A <- matrix(rnorm(N^2,0,1), nrow=N, ncol=N)
Sigma <- t(A) %*% A 

R <- mvrnorm(K, mu, Sigma) # create scenarios
means <- apply(R, 2, mean) # computes mean for each product

### The LP

# There are some additional constraints to pure expectation maximization
# This leads to additional variables 

c <- c(-means,0,rep(0,K))
A1 <- c(rep(1,N),0,rep(0,K))
A2 <- c(rep(0,N),1,-rep((0.05*K)^(-1),K))
A3 <- cbind(R,rep(-1,K),diag(1,K))
A <- rbind(A1,A2,A3)
b <- c(1,98,rep(0,K))

system.time(lp <- lp(direction = "min", objective.in = c, const.mat = A, 
                     const.dir = c("=", ">=", rep(">=",K)), const.rhs = b))

Provide a sample data and code you have tried so far. Please visit [How to make a great R reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). — UseR10085, Oct 28 '20 at 08:19
Gurobi comes with an R interface. There are also R interfaces for Cplex. — Erwin Kalvelagen, Oct 28 '20 at 08:20
`CVXR` if your problem is convex. It is fast. I second @BappaDas, please give more information, the kind of problem you want to solve and how large is the data. — Stéphane Laurent, Oct 28 '20 at 08:57
The ROI optimization package provides access to a suite of alternative solvers including free (no license required): http://roi.r-forge.r-project.org/ The related ompr package https://dirkschumacher.github.io/ompr/ provides a more intuitive modeling platform that can call ROI functions for execution. — SteveM, Oct 28 '20 at 11:01

score 2 · Answer 1 · edited Nov 04 '20 at 11:31

You can try the Rglpk package. In the kantorovich package, there are two functions which solve the same linear programming problem: one using Rglpk, the other one using lpSolve. Benchmarks show that the first one is faster when the data is large.

library(kantorovich)
library(microbenchmark)

mu <- rpois(50, 20)
mu <- mu/sum(mu)
nu <- rpois(50, 20)
nu <- nu/sum(nu)

microbenchmark(
  Rglpk = kantorovich_glpk(mu, nu),
  lpSolve = kantorovich_lp(mu, nu),
  times = 3
)
# Unit: milliseconds
#    expr       min        lq     mean    median        uq       max neval cld
#   Rglpk  402.0955  552.3754  605.159  702.6553  706.6908  710.7264     3  a 
# lpSolve 1092.7131 1184.7517 1327.208 1276.7904 1444.4552 1612.1200     3   b

Edit: better with 'CVXR'

I tried and CVXR with the ECOS solver is faster:

Unit: milliseconds
      expr       min        lq      mean    median        uq       max neval cld
     Rglpk  364.2730  378.3198  383.0741  392.3666  392.4746  392.5827     3  b 
   lpSolve 1012.4465 1087.0728 1128.5777 1161.6990 1186.6433 1211.5875     3   c
      CVXR  370.5944  386.2229  392.8022  401.8515  403.9061  405.9607     3  b 
 CVXR_GLPK  483.2246  488.4495  508.6683  493.6744  521.3902  549.1060     3  b 
 CVXR_ECOS  219.0252  222.9561  224.3361  226.8871  226.9915  227.0960     3 a

Fast Simplex in R

1 Answers1