###################################################################### # Least squares estimate wih lasso penalty # Peng Zeng @ Auburn University # 11-17-2021 ###################################################################### require(gurobi); ###################################################################### # min 0.5 * sum((y - x * b)^2) + lambda * sum(abs(b)) # solve it using by Gurobi # # Input: # x --- n-by-p matrix, predictors # y --- n-by-1 vector, response # lambda --- scalar, nonnegative, tuning parameter in the penalty # ... --- more parameters for solve() in CVXR # # Output: # b --- p-by-1 vector, estimated coefficients # Qvalue --- scalar, value of the objective function # gurobi --- list, results of QP problem by Gurobi ###################################################################### lasso.fit.gurobi = function(x, y, lambda = 0.0, ...) { stopifnot(NROW(x) == length(y)); stopifnot(lambda >= 0.0); n = NROW(x); p = NCOL(x); xtx = t(x) %*% x; xty = t(x) %*% y; ###################################################################### # prepare for Quadratic Programming ###################################################################### model = list(modelsense = "min"); ############################################################ # min 0.5 * b' * S * b - a' * b + sum(abs(lambda * b)) # where S = xtx and a = xty # arguments: b+, b- (length: p + p) ############################################################ model$Q = 0.5 * rbind(cbind(xtx, -xtx), cbind(-xtx, xtx)); model$obj = c(lambda - xty, lambda +xty); model$objcon = 0.5 * sum(y * y); model$lb = rep(0.0, p + p); model$ub = rep(Inf, p + p); ############################################################ # no constraints needed, nrow = 0 ############################################################ model$A = matrix(0.0, nrow = 0, ncol = p + p); model$sense = rep(">", 0); model$rhs = rep(0.0, 0); ###################################################################### # call gurobi() to solve the Quadratic Programming problem ###################################################################### params = list(OutputFlag = 0, ...); QPfit = gurobi(model, params); if(QPfit$status != "OPTIMAL") warning("\n\nThe status of gurobi(model, params) is ", QPfit$status, "\n===== Something is wrong! =====\n"); b = QPfit$x[1:p] - QPfit$x[(p+1):(p+p)]; list(x = x, y = y, lambda = lambda, b = b, Qvalue = QPfit$objval, gurobi = QPfit); } ###################################################################### # THE END ######################################################################