Lasso and Square-Root Lasso by QP and SOCP

Dr. Hua Zhou, Mar 15, 2015

This code snippet demonstrates solving lasso and square-root lasso via Quadratic Programming (QP) and Second Order Cone Programming (SOCP) using the Convex.jl package in Julia.

Display system information.

versioninfo()

Julia Version 0.3.6
Commit 0c24dca (2015-02-17 22:12 UTC)
Platform Info:
  System: Darwin (x86_64-apple-darwin13.4.0)
  CPU: Intel(R) Core(TM) i7-3720QM CPU @ 2.60GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

Simulate a data set

using Gadfly

# random seed
srand(2015790)
# n = # obs, p = # predictors excluding intercept
n = 100; p = 10
# Design matrix
X = [ones(n, 1) randn(n, p)]
# True regression coefficients (first 5 are non-zero)
β = [1.0, randn(5, 1), zeros(p - 5, 1)]
# Responses
y = X * β + randn(n, 1)

# plot the true β
plot(x=1:p, y=β[2:end], Geom.point, Guide.title("True β"), 
    Guide.xlabel("j"), Guide.ylabel("βj"))

Solution path of lasso

Load the Convex.jl package.

using Convex

# Use Gurobi solver
using Gurobi
solver = GurobiSolver(OutputFlag=0)
set_default_solver(solver)

## Use Mosek solver
#using Mosek
#solver = MosekSolver(LOG=1)
#set_default_solver(solver)

## Use SCS solver
#using SCS
#solver = SCSSolver(verbose=1)
#set_default_solver(solver)

GurobiSolver({(:OutputFlag,0)})

Set up a grid of \(\lambda\) and solve the lasso problem

at each value of \(\lambda\).

λlist = 0:10:200
βpath = zeros(p + 1, length(λlist))
# βhat are optimization variables
βhat = Variable(size(X, 2))
@time for i = 1:length(λlist)
    λ = λlist[i]
    problem = minimize(0.5 * sum_squares(y - X * βhat) + λ * norm(βhat[2:end], 1))
    solve!(problem)
    βpath[:, i] = βhat.value
end

elapsed time: 0.297374913 seconds (16303756 bytes allocated)

Display solution path.

using DataFrames

lassopath = DataFrame(Predictor = rep(1:p, length(λlist)), 
    Coeff = reshape(βpath[2:end, :], (p * length(λlist))), 
    λ = rep(λlist, each=p))

plot(lassopath, x=:λ, y=:Coeff, color=:Predictor, Geom.line,
    Guide.xlabel("λ"), Guide.ylabel("βhat(λ)"), Guide.title("Lasso"),
    Scale.color_discrete_hue)

Solution path of square-root lasso

Set up a grid of \(\lambda\) and solve the lasso problem

at each value of \(\lambda\).

λlist = 0:0.5:10
βpath = zeros(p + 1, length(λlist))
# βhat are optimization variables
βhat = Variable(size(X, 2))
@time for i = 1:length(λlist)
    λ = λlist[i]
    problem = minimize(norm(y - X * βhat) + λ * norm(βhat[2:end], 1))
    solve!(problem)
    βpath[:, i] = βhat.value
end

elapsed time: 0.169515327 seconds (13903876 bytes allocated)

Display solution path.

sqrtlassopath = DataFrame(Predictor = rep(1:p, length(λlist)), 
    Coeff = reshape(βpath[2:end, :], (p * length(λlist))), 
    λ = rep(λlist, each=p))
plot(sqrtlassopath, x=:λ, y=:Coeff, color=:Predictor, Geom.line,
    Guide.xlabel("λ"), Guide.ylabel("βhat(λ)"), Guide.title("Sqrt-Root Lasso"),
    Scale.color_discrete_hue)