Fast non-negative quantile and Huber regression in R

Question

I am looking for a quick way to do a non-negative quantile and Huber regression in R (i.e. with the restriction that all coefficients are> 0). I tried using CVXR for quantile and Huber regression and quantreg for quantile regression, but CVXR very slow and quantreg seems to be wrong when I use non-negativity constraints. Does anyone know a good and quick solution in R, for example. using the Rcplex package or the gurobi R API , thereby using the faster CPLEX or gurobi optimizers?

Note that I need to run the problem size, for example, below 80,000 times, as a result of which I need to update the vector y at each iteration, but still use the same X predictor matrix. In this sense, I consider it inefficient that in CVXR I need to do obj <- sum(quant_loss(y - X %*% beta, tau=0.01)); prob <- Problem(Minimize(obj), constraints = list(beta >= 0)) obj <- sum(quant_loss(y - X %*% beta, tau=0.01)); prob <- Problem(Minimize(obj), constraints = list(beta >= 0)) in each iteration, when the problem actually remains unchanged, and all I want to update is y . Any thoughts on making it all better / faster?

Minimal example:

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

Non-negative regression of quantiles using 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 

The syntax for Huber's non-negative regression is the same, but will use

 M <- 1 ## Huber threshold obj <- sum(CVXR::huber(y - X %*% beta, M)) 

Non-negative regression of quantiles using the quantreg package:

 ### 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... 

Answer 1

To speed up CVXR, you can get the problem data once at the beginning, and then change it in a loop and pass it directly to the solver interface R. The code for this

 prob_data <- get_problem_data(prob, solver = "SCS") 

Then we analyze the arguments and pass them to scs from the scs library. (See Solver.solve in solver.R). You will need to delve into the details of canonization, but I expect that if you just change y at each iteration, it should be a simple modification.