我正在寻找在R中进行非负分位数和Huber回归的快速方法(即所有系数都> 0的约束条件)。我尝试使用
请注意,我需要运行像下面80,000次这样的问题规模,每次迭代中只需更新
最小示例:
使用CVXR进行非负分位数回归:
CVXR
软件包进行分位数和Huber回归以及quantreg
软件包进行分位数回归,但是CVXR
非常慢,当我使用非负性约束时,quantreg
似乎有缺陷。是否有人知道在R中的一个好的快速解决方案,例如使用Rcplex
软件包或R gurobi API,从而使用更快的CPLEX或gurobi优化器?请注意,我需要运行像下面80,000次这样的问题规模,每次迭代中只需更新
y
向量,但仍然使用相同的预测矩阵X
。在这个意义上,我觉得在CVXR
中,我现在必须在每次迭代中执行obj <- sum(quant_loss(y - X %*% beta, tau=0.01)); prob <- Problem(Minimize(obj), constraints = list(beta >= 0))
,当事实上问题保持不变,我想要更新的只是y
。有什么想法可以更好/更快地解决所有这些问题吗?最小示例:
## Generate problem data
n <- 7 # n predictor vars
m <- 518 # n cases
set.seed(1289)
beta_true <- 5 * matrix(stats::rnorm(n), nrow = n)+20
X <- matrix(stats::rnorm(m * n), nrow = m, ncol = n)
y_true <- X %*% beta_true
eps <- matrix(stats::rnorm(m), nrow = m)
y <- y_true + eps
使用CVXR进行非负分位数回归:
## Solve nonnegative quantile regression problem using CVX
require(CVXR)
beta <- Variable(n)
quant_loss <- function(u, tau) { 0.5*abs(u) + (tau - 0.5)*u }
obj <- sum(quant_loss(y - X %*% beta, tau=0.01))
prob <- Problem(Minimize(obj), constraints = list(beta >= 0))
system.time(beta_cvx <- pmax(solve(prob, solver="SCS")$getValue(beta), 0)) # estimated coefficients, note that they ocasionally can go - though and I had to clip at 0
# 0.47s
cor(beta_true,beta_cvx) # correlation=0.99985, OK but very slow
非负 Huber 回归的语法相同,但需要使用以下内容:
M <- 1 ## Huber threshold
obj <- sum(CVXR::huber(y - X %*% beta, M))
使用 quantreg
包进行非负分位数回归:
### Solve nonnegative quantile regression problem using quantreg package with method="fnc"
require(quantreg)
R <- rbind(diag(n),-diag(n))
r <- c(rep(0,n),-rep(1E10,n)) # specify bounds of coefficients, I want them to be nonnegative, and 1E10 should ideally be Inf
system.time(beta_rq <- coef(rq(y~0+X, R=R, r=r, tau=0.5, method="fnc"))) # estimated coefficients
# 0.12s
cor(beta_true,beta_rq) # correlation=-0.477, no good, and even worse with tau=0.01...