###################################################################### # Quantile regression with generalized lasso penalty and linear constraints # Peng Zeng @ Auburn University # 04-03-2020 ###################################################################### require(gurobi); ###################################################################### # min sum(rho_{tau}(y - x * b)) + lambda * |D * b|_1 # subject to C * b >= d, E * b = f # solve it using Linear Programming by Gurobi # # Input: # x --- n-by-p matrix, predictors # y --- n-by-1 vector, response # lambda --- scalar, nonnegative, tuning parameter in the penalty # tau --- scalar, value in (0, 1), median corresponds to 0.5 # Cmat --- q-by-p matrix, inequality constraints # dvec --- q-by-1 vector, inequality constraints # Emat --- s-by-p matrix, equality constraints # fvec --- s-by-1 vector, equality constraints # Dmat --- m-by-p matrix, default: identity matrix # ... --- more parameters for gurobi # # Output: # b --- p-by-1 vector, estimated coefficients # r --- n-by-1 vector, r = y - x %*% b # z --- m-by-1 vector, z = D %*% b # Qvalue --- scalar, value of the objective function # model --- list, the LP problem to be solved by Gurobi # gurobi --- list, results of LP problem by Gurobi ###################################################################### quantreg.genlasso.fit.gurobi = function( x, y, lambda = 0.0, tau = 0.5, Cmat = NULL, dvec = NULL, Emat = NULL, fvec = NULL, Dmat = diag(rep(1.0, NCOL(x))), ...) { stopifnot(NROW(x) == length(y)); stopifnot(lambda >= 0.0, tau > 0.0, tau < 1.0); stopifnot(NCOL(Dmat) == NCOL(x)); n = NROW(x); p = NCOL(x); m = NROW(Dmat); ###################################################################### # prepare for Linear Programming ###################################################################### model = list(modelsense = 'min'); ###################################################################### # min tau * sum(r+) + (1 - tau) * sum(r-) + lambda * sum(z+ + z-) # arguments: # b -- length p, # r+ -- length n r+ >= 0 # r- -- length n, r- >= 0 # z+ -- length nD, z+ >= 0 # z- -- length nD, z- >= 0 ###################################################################### model$obj = c(rep(0.0, p), rep(tau, n), rep(1.0 - tau, n), rep(lambda, m + m)); model$lb = c(rep(-Inf, p), rep(0.0, n + n + m + m)); model$ub = rep(Inf, p + n + n + m + m); ###################################################################### # contraints: 1 to n # x * b + r+ - r- = y ###################################################################### model$A = cbind(x, eye(n, 1.0), eye(n, -1.0), zeros(n, m + m)); model$rhs = y; model$sense = rep('=', n); ###################################################################### # constraints: n + (1 to nD) # D * b - z+ + z- = 0 or D * b = z ###################################################################### model$A = rbind(model$A, cbind(Dmat, zeros(m, n + n), eye(m, -1.0), eye(m, 1.0))); model$rhs = c(model$rhs, rep(0.0, m)); model$sense = c(model$sense, rep('=', m)); ###################################################################### # constraints: n + nD + (1 to nC) # C * b >= d ###################################################################### if(!is.null(Cmat)) { if(is.null(dvec)) dvec = rep(0.0, NROW(Cmat)); stopifnot(NROW(Cmat) == length(dvec), NCOL(Cmat) == NCOL(x)); q = NROW(Cmat); # update A, rhs, sense model$A = rbind(model$A, cbind(Cmat, zeros(q, n + n + m + m))); model$rhs = c(model$rhs, dvec); model$sense = c(model$sense, rep('>=', q)); } ###################################################################### # constraints: n + nD + nC + (1 to nE) # E * b >= f ###################################################################### if(!is.null(Emat)) { if(is.null(fvec)) fvec = rep(0.0, NROW(Emat)); stopifnot(NROW(Emat) == length(fvec), NCOL(Emat) == NCOL(x)); s = NROW(Emat); # update A, rhs, sense model$A = rbind(model$A, cbind(Emat, zeros(s, n + n + m + m))); model$rhs = c(model$rhs, fvec); model$sense = c(model$sense, rep('=', s)); } ###################################################################### # call gurobi() to solve the Quadratic Programming problem ###################################################################### params = list(OutputFlag = 0, ...); LPfit = gurobi(model, params); if(LPfit$status != 'OPTIMAL') warning('\n\nThe status of gurobi(model, params) is ', LPfit$status, '\n===== Something is wrong! =====\n'); ###################################################################### # get the solution from gurobi() # r = r+ - r- = y - x * b # z = z+ - z- = D * b ###################################################################### b = LPfit$x[1:p]; r = LPfit$x[p + (1:n)] - LPfit$x[p + n + (1:n)]; z = LPfit$x[p + n + n + (1:m)] - LPfit$x[p + n + n + m + (1:m)]; list(x = x, y = y, lambda = lambda, tau = tau, Cmat = Cmat, dvec = dvec, Emat = Emat, fvec = fvec, Dmat = Dmat, b = b, r = r, z = z, Qvalue = LPfit$objval, model = model, gurobi = LPfit); } ###################################################################### # n-by-p matrix of zeros ###################################################################### zeros = function(n, p) { matrix(0.0, nrow = n, ncol = p); } ###################################################################### # p-by-p diagonal matrix with value as its diagonal elements ###################################################################### eye = function(p, value = 1.0) { if(p == 1) matrix(value, nrow = 1) else diag(rep(value, p)); } ###################################################################### # THE END ######################################################################