Lasso by Coordinate Descent

Fit lasso regression with cyclic coordinate descent.

Lasso adds an L1 penalty to least squares:

loss(w) = (1 / (2n)) * ||y - X w||^2 + lam * sum_j |w_j|

Unlike ridge, the L1 penalty sets some weights to exactly zero, so lasso also selects features. The penalty has no derivative at zero, so plain gradient descent handles it badly. Coordinate descent instead solves for one weight at a time with the others held fixed, and that one-dimensional problem has an exact answer given by soft thresholding:

r_j   = y - X w + X[:, j] * w_j          # residual without feature j
rho_j = (1 / n) * X[:, j] . r_j
z_j   = (1 / n) * X[:, j] . X[:, j]
w_j   = soft(rho_j, lam) / z_j
soft(a, lam) = sign(a) * max(|a| - lam, 0)

Start from w = 0. One sweep updates j = 0, 1, ..., d - 1 in order, each update using the latest values of the other weights. Run n_iters sweeps.

Function signature

def lasso_cd(X, y, lam, n_iters=100):

Example

Input

X = np.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0], [0.0, 1.0]])
y = np.array([3.0, 1.0, 1.0, -1.0])
lasso_cd(X, y, lam=0.25)

Output

array([1.5, 0. ])

Reasoning

The two columns are orthogonal, so each weight is solved exactly in the first sweep. For feature 0, rho = (3 + 1) / 4 = 1 and z = 2 / 4 = 0.5. Shrinking 1 by 0.25 gives 0.75, and 0.75 / 0.5 = 1.5. Least squares alone would give 2. For feature 1, rho = (1 - 1) / 4 = 0, which is inside the threshold, so the weight is exactly 0.

Constraints and notes

• There is no intercept. Centre the data yourself if you need one.
• No column of X is all zeros.
• With lam = 0 the method converges to ordinary least squares.