Lasso paths, refitting, and choosing λ
LassoPath solves the lasso
minimize ‖y - β₀ - Xβ‖² / (2n) + λ ∑ⱼ ωⱼ |βⱼ|along a decreasing sequence of λ, each from the previous solution; refitlassopath refits least squares on the selected variables. This example shows the options of both and chooses λ by cross-validation.
using HighDimensionalModels
using ProximalBase: MCP
using LinearAlgebra, Printf, Random, Statistics
Random.seed!(3)Random.TaskLocalRNG()more predictors than observations: n = 100, p = 500, ten nonzero coefficients between 0.5 and 2 in absolute value, and an intercept
n, p, s = 100, 500, 10
X = randn(n, p)
β = zeros(p)
β[1:s] = range(2, 0.5; length = s) .* rand((-1, 1), s)
y = 1.0 .+ X * β + randn(n)
nonzeros(b) = count(!iszero, b)nonzeros (generic function with 1 method)The default path
100 values of λ, evenly spaced on the log scale from λmax, the smallest λ with a zero solution, down to λminratio ⋅ λmax, with λminratio = 0.01 when p > n and 1e-4 otherwise; loadings ωⱼ = standard deviation of X[:, j] (standardize_x = true); no intercept unless asked for; and screening with the sequential strong rule
path = LassoPath(X, y; intercept = true)
println(" λ nonzeros β̂₀ passes KKT checks")
for k in 1:11:length(path.λpath)
@printf(" %7.4f %8d %5.2f %6d %10d\n", path.λpath[k], nonzeros(path.βpath[k]),
path.β0path[k], path.iterations[k], path.kkt_checks[k])
end
println("all converged: ", all(path.converged)) λ nonzeros β̂₀ passes KKT checks
2.4330 0 1.06 1 1
1.4586 2 0.97 5 1
0.8744 7 0.95 6 1
0.5242 10 0.92 5 1
0.3142 19 0.89 9 1
0.1884 39 0.93 12 1
0.1129 59 0.93 20 1
0.0677 76 0.96 34 1
0.0406 88 0.99 36 1
0.0243 94 1.00 48 1
all converged: truepath.βpath[k] is a SparseVector, so a path takes memory in proportion to the number of nonzero coefficients
println("βpath[50] is a ", typeof(path.βpath[50]), " with ", nonzeros(path.βpath[50]), " nonzeros")βpath[50] is a SparseArrays.SparseVector{Float64, Int64} with 27 nonzerosOther grids
fewer values on a shorter range ...
short = LassoPath(X, y; nλ = 20, λminratio = 0.1, intercept = true)
@printf("nλ = 20, λminratio = 0.1: λ from %.3f down to %.3f\n", first(short.λpath), last(short.λpath))nλ = 20, λminratio = 0.1: λ from 2.433 down to 0.243... a grid of your own, which should decrease ...
λs = exp.(range(log(1.0), log(0.05); length = 30))
own = LassoPath(X, y, λs; intercept = true)
@printf("own grid: %d values, %d nonzeros at λ = %.2f\n", length(own.λpath),
nonzeros(own.βpath[end]), own.λpath[end])own grid: 30 values, 87 nonzeros at λ = 0.05... or a path that stops after the first solution with more than maxhats nonzero coefficients; λpath then holds only the values solved for
early = LassoPath(X, y; intercept = true, max_hat_s = 20)
@printf("max_hat_s = 20: stops after %d of 100 values, with %d nonzeros\n",
length(early.λpath), nonzeros(early.βpath[end]))max_hat_s = 20: stops after 47 of 100 values, with 22 nonzerosUnstandardized columns
With standardizex = false the loadings are 1, and columns with larger standard deviations are penalized less (see 01lasso.jl). Here the columns have standard deviations mostly between 0.1 and 10, and the variables that enter first are those with large columns rather than those that matter.
Xs = X .* exp.(1.2 * randn(1, p))
function entry_order(path) # variables in the order they become nonzero
order = Int[]
for b in path.βpath, j in findall(!iszero, b)
j in order || push!(order, j)
end
order
end
for standardize_x in (true, false)
path_s = LassoPath(Xs, y .- 1.0; standardize_x)
println("standardize_x = ", rpad(standardize_x, 5), ": first ten variables to enter ",
first(entry_order(path_s), 10))
endstandardize_x = true : first ten variables to enter [1, 2, 3, 5, 4, 10, 6, 7, 322, 8]
standardize_x = false: first ten variables to enter [378, 468, 79, 280, 473, 5, 55, 435, 347, 95]Screening and solver options
With screening = true (the default), each λ is solved first over the variables kept by the strong rule, and the optimality conditions of the others are checked afterwards; violators are added and the problem solved again. kkt_checks[k] counts the checks, 1 when the rule missed nothing. screening = false visits all p variables in every full pass. The solutions are the same up to the accuracy of the solver, which is lowest at the end of the path, where almost n variables are selected. Screening saves time in the full passes when most of many variables stay at zero, as with n = 100 and p = 5000 here.
Xw = randn(100, 5000)
yw = Xw[:, 1:s] * β[1:s] + randn(100)
time_path(; kwargs...) = minimum(@elapsed(LassoPath(Xw, yw; kwargs...)) for _ in 1:3)
p1 = LassoPath(Xw, yw)
p0 = LassoPath(Xw, yw; screening = false)
Δ = [norm(p1.βpath[k] - p0.βpath[k], Inf) for k in eachindex(p1.λpath)]
@printf("largest difference between the two paths: %.1e over the first 80 values of λ, %.1e over the last 20\n",
maximum(Δ[1:80]), maximum(Δ[81:end]))
@printf("time: %.0f ms with screening, %.0f ms without\n", 1000time_path(), 1000time_path(; screening = false))
println("values of λ at which the strong rule missed a variable: ", count(>(1), p1.kkt_checks))largest difference between the two paths: 9.8e-04 over the first 80 values of λ, 3.9e-03 over the last 20
time: 29 ms with screening, 88 ms without
values of λ at which the strong rule missed a variable: 0Options of coordinate descent go in the last positional argument, here a tighter tolerance; screening needs warm starts, so it is off with warm_start = false
tight = LassoPath(Xw, yw, CDOptions(; opt_tol = 1e-9))
println("with opt_tol = 1e-9: ", sum(tight.iterations), " passes along the path instead of ", sum(p1.iterations))with opt_tol = 1e-9: 3823 passes along the path instead of 3801Choosing λ by cross-validation
K-fold cross-validation on the λ grid of the full data: the path is fit without one fold, and its predictions are scored on that fold
function cvlasso(X, y, λpath; K = 5, kwargs...)
n = size(X, 1)
fold = shuffle!(mod1.(1:n, K))
mse = zeros(length(λpath), K)
for k in 1:K
train, test = fold .!= k, fold .== k
pathk = fit(LassoPath, X[train, :], y[train], λpath; intercept = true, kwargs...)
mse[:, k] = vec(mean(abs2, y[test] .- predict(pathk, X[test, :]); dims = 1))
end
vec(mean(mse; dims = 2)), vec(std(mse; dims = 2)) / sqrt(K)
end
cvm, cvse = cvlasso(X, y, path.λpath)
imin = argmin(cvm)
i1se = findfirst(cvm .<= cvm[imin] + cvse[imin]) # largest λ within one standard error
println(" λ nonzeros CV error ‖β̂ - β‖")
for (name, i) in (("minimum CV error", imin), ("one standard error", i1se))
@printf(" %-18s %6.4f %8d %8.3f %8.3f\n", name, path.λpath[i], nonzeros(path.βpath[i]),
cvm[i], norm(path.βpath[i] - β))
end λ nonzeros CV error ‖β̂ - β‖
minimum CV error 0.1183 55 3.327 0.956
one standard error 0.2863 22 3.991 1.293The elastic net and MCP
LassoPath(X, y; alpha) is glmnet's elastic-net path, whose λmax is that of the lasso divided by α, and LassoPath(X, y; penalty = MCP(1, 3)) the path of ncvreg, each solution the local one reached from the previous. MCP stops shrinking coefficients larger than aλ, so its estimates at the same λ are less biased; cvlasso works on any path.
for (name, kw) in (("elastic net, α = 0.5", (alpha = 0.5,)), ("MCP, a = 3", (penalty = MCP(1, 3),)))
pk = LassoPath(X, y; intercept = true, kw...)
ck, _ = cvlasso(X, y, pk.λpath; kw...)
i = argmin(ck)
@printf(" %-20s λmax %6.3f at the minimum CV error: λ = %6.4f, %3d nonzeros, CV error %6.3f, ‖β̂ - β‖ = %.3f\n",
name, pk.λpath[1], pk.λpath[i], nonzeros(pk.βpath[i]), ck[i], norm(pk.βpath[i] - β))
end elastic net, α = 0.5 λmax 4.866 at the minimum CV error: λ = 0.1875, 95 nonzeros, CV error 6.946, ‖β̂ - β‖ = 2.229
MCP, a = 3 λmax 2.433 at the minimum CV error: λ = 0.2067, 16 nonzeros, CV error 1.097, ‖β̂ - β‖ = 0.278Least squares refits
refitlassopath maps every distinct support S along the path to the least squares fit on X[:, S], with the intercept first when the path has one. Refitting removes the shrinkage of the lasso from the selected coefficients.
refits = refit_lasso_path(path, X, y)
println(length(refits), " distinct supports along the path")
S = findall(!iszero, path.βpath[i1se])
βS = refits[S]
βrefit = zeros(p)
βrefit[S] = βS[2:end]
@printf("on the %d variables selected with one standard error: ‖β̂ - β‖ = %.3f for the lasso, %.3f refit\n",
length(S), norm(path.βpath[i1se] - β), norm(βrefit - β))
println("true ", round.(β[1:s]; digits = 2))
println("lasso ", round.(Vector(path.βpath[i1se][1:s]); digits = 2))
println("refit ", round.(βrefit[1:s]; digits = 2))65 distinct supports along the path
on the 22 variables selected with one standard error: ‖β̂ - β‖ = 1.293 for the lasso, 0.659 refit
true [2.0, -1.83, -1.67, 1.5, -1.33, -1.17, 1.0, 0.83, 0.67, -0.5]
lasso [1.69, -1.48, -1.2, 1.04, -0.97, -0.75, 0.6, 0.35, 0.12, -0.52]
refit [1.86, -1.77, -1.49, 1.45, -1.33, -1.02, 0.88, 0.79, 0.62, -0.59]