The lasso

minimize  ‖y - β₀ - Xβ‖² / (2n) + λ ∑ⱼ ωⱼ |βⱼ|

with lasso and its parameters: the penalty level λ, an unpenalized intercept β₀, standardized columns, and penalty loadings ωⱼ.

using HighDimensionalModels
using ProximalBase: MCP, SCAD
using LinearAlgebra, Printf, Random, Statistics

Random.seed!(1)

Random.TaskLocalRNG()

n = 200 observations of p = 50 predictors, of which the first five matter

n, p = 200, 50
X = randn(n, p)
β = [3.0, -2.0, 1.5, -1.0, 0.5, zeros(p - 5)...]
y = X * β + randn(n)

support(x) = findall(!iszero, x)

support (generic function with 1 method)

The penalty level λ

λ = σ √(2 log(p) / n) is the choice of the theory, with noise level σ = 1 here; 03lassopath.jl chooses λ by cross-validation and 04sqrtscaled_feasible.jl without knowing σ.

λ = sqrt(2 * log(p) / n)
sol = lasso(X, y, λ)

@printf("λ = %.3f selects %s\n", λ, support(sol.x))
println("β̂[1:5] = ", round.(sol.x[1:5]; digits = 2), ", shrunk towards zero")
println("σ̂ = ", round(sol.σ; digits = 3), ", the standard deviation of the residuals")
println("$(sol.iterations) passes over the coordinates, converged: $(sol.converged)")

λ = 0.198 selects [1, 2, 3, 4, 5]
β̂[1:5] = [2.81, -1.81, 1.36, -0.73, 0.4], shrunk towards zero
σ̂ = 1.081, the standard deviation of the residuals
10 passes over the coordinates, converged: true

a larger λ gives a sparser solution; all coefficients are zero for λ ≥ λmax = maxⱼ |X[:, j]'y| / n

λmax = maximum(abs, X' * y) / n
println("\n     λ  nonzeros  ‖β̂ - β‖")
for λ in (1.01λmax, 0.5λmax, 0.2, 0.1, 0.05, 0.01)
  s = lasso(X, y, λ)
  @printf("%6.3f  %8d  %7.3f\n", λ, count(!iszero, s.x), norm(s.x - β))
end

     λ  nonzeros  ‖β̂ - β‖
 3.460         0    4.062
 1.713         2    2.893
 0.200         5    0.425
 0.100        15    0.265
 0.050        30    0.300
 0.010        44    0.451

An intercept

With intercept = true, the intercept is not penalized and sol.intercept = ȳ - x̄'β̂. The columns and y are centered on the fly, so X is not copied. Without an intercept, the predictors have to explain the mean of y as well.

y2 = 2.0 .+ y
for intercept in (false, true)
  s = lasso(X, y2, λ; intercept)
  @printf("intercept = %-5s  β̂₀ = %.2f  RMS of the residuals %.2f  selects %s\n",
          intercept, s.intercept, sqrt(mean(abs2, residuals(s))), support(s.x))
end

intercept = false  β̂₀ = 0.00  RMS of the residuals 2.08  selects [1, 2, 3, 4, 5, 11, 14, 19, 26, 27, 29, 31, 34, 36, 39, 42, 43, 45, 49]
intercept = true   β̂₀ = 1.96  RMS of the residuals 1.08  selects [1, 2, 3, 4, 5]

predictions at new points, and the other functions of StatsAPI: coef(sol) is [β̂₀; β̂] with an intercept, dof(sol) the number of nonzero coefficients plus one for the intercept

sol = lasso(X, y2, λ; intercept = true)
Xnew = randn(3, p)
ŷ = predict(sol, Xnew)                         # sol.intercept .+ Xnew * sol.x
println("predictions ", round.(ŷ; digits = 2), ", true means ", round.(2.0 .+ Xnew * β; digits = 2))
println("coef(sol)[1:3] = ", round.(Vector(coef(sol))[1:3]; digits = 2), ", nobs = ", nobs(sol), ", dof = ", dof(sol))

predictions [1.51, 2.2, 0.86], true means [1.52, 2.79, 0.81]
coef(sol)[1:3] = [1.96, 2.81, -1.81], nobs = 200, dof = 6

Standardization

The penalty λ ∑ⱼ |βⱼ| depends on the units of the columns. Measuring predictor 10, which does not matter, in centimeters instead of meters multiplies its column by 100, so that a coefficient a hundred times smaller fits the same, at a hundredth of the penalty. standardize_x = true takes the loadings ωⱼ = standard deviation of X[:, j], which gives the lasso on standardized columns with the coefficients transformed back: the fit then does not depend on the units.

Xcm = copy(X)
Xcm[:, 10] .*= 100
for standardize_x in (false, true)
  a = lasso(X, y, λ; standardize_x)
  b = lasso(Xcm, y, λ; standardize_x)
  @printf("standardize_x = %-5s  selects %s in meters and %s in centimeters\n",
          standardize_x, support(a.x), support(b.x))
end

standardize_x = false  selects [1, 2, 3, 4, 5] in meters and [1, 2, 3, 4, 5, 10] in centimeters
standardize_x = true   selects [1, 2, 3, 4, 5] in meters and [1, 2, 3, 4, 5] in centimeters

Penalty loadings ωⱼ

lasso(X, y, λ, ω) penalizes λ ωⱼ |βⱼ|. A zero loading leaves a coefficient unpenalized, so that its variable is always in the model, for example a treatment whose effect is of interest.

ω = ones(p)
ω[[20, 30]] .= 0
sol = lasso(X, y, λ, ω)
println("with ω₂₀ = ω₃₀ = 0 the lasso selects ", support(sol.x))

with ω₂₀ = ω₃₀ = 0 the lasso selects [1, 2, 3, 4, 5, 20, 30]

The adaptive lasso (Zou, 2006) takes ωⱼ = 1 / |β̃ⱼ| from a first estimate β̃, here least squares (n > p). Large coefficients are penalized less, so they are less biased, and small ones more, so fewer variables that do not matter are selected as λ decreases. adaptive_lasso(X, y, λ; init = β̃) is lasso(X, y, λ, 1 ./ abs.(β̃)).

β̃ = X \ y
println("\n           lasso               adaptive lasso")
println("     λ  nonzeros  ‖β̂ - β‖   nonzeros  ‖β̂ - β‖")
for λ in (0.2, 0.1, 0.05, 0.02)
  a = lasso(X, y, λ)
  b = adaptive_lasso(X, y, λ; init = β̃)
  @printf("%6.2f  %8d  %7.3f   %8d  %7.3f\n",
          λ, count(!iszero, a.x), norm(a.x - β), count(!iszero, b.x), norm(b.x - β))
end

           lasso               adaptive lasso
     λ  nonzeros  ‖β̂ - β‖   nonzeros  ‖β̂ - β‖
  0.20         5    0.425          5    0.472
  0.10        15    0.265          5    0.247
  0.05        30    0.300          5    0.155
  0.02        39    0.405          7    0.130

By default the first estimate is the lasso at the same λ (init = :lasso, or :ridge, at the level initλ): a variable that it drops stays out.

c = adaptive_lasso(X, y, 0.05)
println("from a first lasso at λ = 0.05, the adaptive lasso selects ", support(c.x))

from a first lasso at λ = 0.05, the adaptive lasso selects [1, 2, 3, 4, 5]

Other penalties

alpha < 1 gives the elastic net λ(α‖β‖₁ + (1 - α)‖β‖²/2) of glmnet, which keeps groups of correlated variables together; penalty gives MCP or SCAD, whose level λ replaces. With standardize_x the penalty applies to the coefficients of the standardized columns, as in glmnet and ncvreg. At the level λ of the theory, MCP and SCAD leave the coefficients larger than aλ unshrunk, so their estimates are closer to β.

for (name, kw) in (("lasso", (;)), ("elastic net, α = 0.5", (alpha = 0.5,)),
                   ("MCP, a = 3", (penalty = MCP(1, 3),)), ("SCAD, a = 3.7", (penalty = SCAD(1, 3.7),)))
  est = lasso(X, y, λ; intercept = true, standardize_x = true, kw...)
  @printf("%-22s selects %2d variables, ‖β̂ - β‖ = %.3f\n", name, count(!iszero, est.x), norm(est.x - β))
end

lasso                  selects  5 variables, ‖β̂ - β‖ = 0.422
elastic net, α = 0.5   selects 19 variables, ‖β̂ - β‖ = 0.577
MCP, a = 3             selects  5 variables, ‖β̂ - β‖ = 0.124
SCAD, a = 3.7          selects  5 variables, ‖β̂ - β‖ = 0.122

Bounds on the coefficients

lower and upper bound the coefficients, numbers or one bound per coefficient, as glmnet's lower.limits and upper.limits; zero must be feasible, and lower = upper = 0 keeps a variable out. lower = 0 gives the nonnegative lasso. Here it is the wrong model: β₂ and β₄ are negative, so they are set to zero, and noise variables enter to fit what they explained.

nn = lasso(X, y, λ; intercept = true, standardize_x = true, lower = 0)
println("lower = 0: β̂[1:5] = ", round.(nn.x[1:5]; digits = 2), ", and ",
        count(!iszero, nn.x[6:end]), " noise variables")
cap = lasso(X, y, λ; intercept = true, standardize_x = true, lower = -1, upper = 1)
println("with -1 ≤ βⱼ ≤ 1, β̂[1:5] = ", round.(cap.x[1:5]; digits = 2))

lower = 0: β̂[1:5] = [2.82, 0.0, 1.26, 0.0, 0.54], and 9 noise variables
with -1 ≤ βⱼ ≤ 1, β̂[1:5] = [1.0, -1.0, 1.0, -0.76, 0.46]

Number types and solver options

lasso converts X, y and λ to a common floating point type: integer data give a Float64 solution and Float32 data a Float32 one

sol32 = lasso(Float32.(X), Float32.(y), 0.2)
println("eltype of the solution from Float32 data: ", eltype(sol32.x))

eltype of the solution from Float32 data: Float32

The last positional argument sets the options of coordinate descent (02solveroptions.jl), here a tighter tolerance

sol = lasso(X, y, λ, CDOptions(; opt_tol = 1e-10))
println("with opt_tol = 1e-10: $(sol.iterations) passes")

with opt_tol = 1e-10: 15 passes