R中的快速单形

编程问答 2022-06-06

问题描述

我不确定这是问不对的地方，但是我需要在R中有效地实现一个非常大的线性编程问题（这里R周围没有办法）。我已经尝试了一些lpSolve之类的软件包，但结果似乎并不令人满意。我很乐意提供有关包装的建议，或者是提出这个问题的好地方。

问题出在这里：

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,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,rep(0,K))
A1 <- c(rep(1,N),K))
A2 <- c(rep(0,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,K))

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

解决方法

您可以尝试使用Rglpk软件包。在kantorovich包中，有两个函数可以解决相同的线性编程问题：一个函数使用Rglpk，另一个函数使用lpSolve。基准测试表明，当数据量较大时，第一个比较快。

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

编辑：使用“ CVXR”更好

我尝试过使用CVXR求解器进行ECOS更快： 

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