######################################################################
# Least squares estimate wih lasso penalty 
# Peng Zeng @ Auburn University 
# 11-17-2021
######################################################################

require(CVXR); 

######################################################################
# min 0.5 * sum((y - x * b)^2) + lambda * sum(abs(b))
# solve it using by CVXR
#
# 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 
#    cvx    --- list, results of LP problem by CVXR 
######################################################################

lasso.fit.cvxr = function(
    x, y, lambda = 0.0, ...)
{
    stopifnot(NROW(x) == length(y)); 
    stopifnot(lambda >= 0.0); 

    n = NROW(x); p = NCOL(x); 

    ######################################################################
    # specify model
    ######################################################################

    bvar = Variable(p);     
    obj = 0.5 * sum_squares(y - x %*% bvar) + lambda * p_norm(bvar, 1); 

    ######################################################################
    # define contraints
    ######################################################################

    prob = Problem(Minimize(obj)); 
    result = solve(prob, ...); 

    b = drop( result$getValue(bvar) ); 

    list(x = x, y = y, lambda = lambda,          
         b = b, Qvalue = result$value, cvx = result);
}

######################################################################
# THE END
######################################################################