Regularisation: ridge and lasso

In lesson 3.1 a model with many parameters learned the noise: the training error fell while the test error rose. Instead of taking features away from a model, you can charge it for them. Regularisation is a penalty added to the training error: the larger the weights, the larger the penalty. A number α sets its strength. Today you will try two penalties on 10 features of patients with diabetes and see that one of them picks the features by itself.

Try it

• Move α in lasso up from the smallest value. In what order do the weight lines reach zero? How many features are left at α = 15.8, and what does that cost on the test?
• Watch s3. When does it first touch zero, and what happens to it after that?
• Switch to ridge and move α all the way up. Does any weight become exactly zero?

Once you have checked it in the lab, let's name what you see.

The idea in plain words

Ordinary linear regression picks the weights that make the error on the training data as small as possible. Regularisation adds a penalty for the size of the weights to that error. The model now has to choose: fit the data better or keep the weights smaller. A feature keeps a large weight only if it really improves the predictions.

Ridge penalises the sum of the squared weights. It pulls every weight towards zero but sets none to zero: in the lab at α = 398 the test R² falls to 0.405 and all 10 features still have a weight. Lasso penalises the sum of the absolute values of the weights and sets them to zero one by one. At α = 15.8 three features are left: BMI, blood pressure and s5. Without a penalty the test R² is 0.470, with three features 0.386. The model is simpler and cheaper to use, but it predicts less well.

The features are expressed in standard deviations. That matters: the penalty counts the size of the weights, so a feature measured in large units would get a small weight and escape the penalty.

Jak to policzono

• Data: Diabetes, Efron et al. (2004), Annals of Statistics. 442 patients, 10 features measured at the start of the study, target: progression of the disease a year later. The authors gave no licence; scikit-learn redistributes the data set (BSD-3-Clause licence).
• Random split: 309 patients to learn from, 133 for the test. Means and standard deviations of the features computed on the training part only.
• Lasso: 11 values of α from 0.01 to 100. Ridge: 11 values of α from 0.1 to 100000 (ridge in scikit-learn counts the penalty on a different scale from lasso, hence the different numbers).
• Test R²: the share of the variation in the test patients' outcome that the model explains (1 = all of it, 0 = no better than the mean).

Why it works

The lasso penalty grows linearly: shrinking a weight by one unit always saves the same amount of penalty, whether the weight is large or small. If a feature improves the fit by less than the penalty costs, it pays to give it a weight of exactly zero. The ridge penalty grows with the square of the weight, so small weights cost almost nothing. Ridge keeps shrinking a weight but never has to set it to zero.

The penalty also calms the weights of features that say almost the same thing. s1 and s2 are two blood test results strongly linked with each other (in the code below their correlation is 0.90). At the weakest lasso penalty s1 has a weight of −39.0 and s2 of +26.2: two large numbers of opposite sign that almost cancel out. Weights like that change a lot from sample to sample. The penalty pulls them in, and soon lasso removes s2 altogether, because s1 carries almost the same information.

This is not a quirk of one split. In the code below the same experiment runs on 50 random splits: BMI is left last in 37 of them, s5 in the other 13. The feature that goes first changes even more from split to split: age in 21 splits, s2 in 16.

Common mistake

True or false: "lasso gave a feature a weight of zero, so that feature does not matter". False. s2 goes first even though it is strongly linked with s1: lasso picks one feature of the pair, and the other gets zero because it adds nothing new next to s1. Without s1, s2 could be important. The lab shows it with s3 too: at α = 0.158 its weight is zero, at α = 0.398 it is back, because it is linked with s4 (r = −0.73) and lasso moves the weight between them.

The choice between correlated features can be a matter of chance. In the code below the first to go is sometimes s2, sometimes age, and even the last one left is sometimes s5 instead of BMI, depending on which patients end up in the training part. A zero in lasso means "this model did not need this feature, with these data and this α". It does not mean "this feature has nothing to do with the disease".

import numpy as np
from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression, Lasso, lars_path
from sklearn.metrics import r2_score

X, y = load_diabetes(return_X_y=True)
names = np.array(["age", "sex", "bmi", "bp", "s1", "s2", "s3", "s4", "s5", "s6"])
alphas = 10 ** np.linspace(-2, 2, 41)  # the lab's range of penalties, finer
last, first, r2 = [], [], []
for s in range(50):
    tr, te = train_test_split(np.arange(len(y)), test_size=0.3, random_state=s)
    sc = StandardScaler().fit(X[tr])
    Xtr, Xte = sc.transform(X[tr]), sc.transform(X[te])
    path_alphas, _, coefs = lars_path(Xtr, y[tr] - y[tr].mean(), method="lasso")  # the exact lasso path, from the strongest penalty down
    enters = [np.flatnonzero(coefs[j])[0] if coefs[j].any() else len(path_alphas) for j in range(10)]  # where each feature first gets a weight
    last.append(names[np.argmin(enters)]); first.append(names[np.argmax(enters)])  # first to get a weight = last to lose it as the penalty grows
    W = np.array([Lasso(alpha=a, max_iter=100000).fit(Xtr, y[tr]).coef_ for a in alphas])
    k3 = min(k for k in range(len(alphas)) if (W[k] != 0).sum() <= 3)  # the weakest penalty that leaves at most 3 features
    m3 = Lasso(alpha=alphas[k3], max_iter=100000).fit(Xtr, y[tr])
    r2.append((r2_score(y[te], LinearRegression().fit(Xtr, y[tr]).predict(Xte)), r2_score(y[te], m3.predict(Xte))))
r2 = np.array(r2)
print(f"last feature lasso keeps, 50 splits: {', '.join(f'{n} {last.count(n)}' for n in sorted(set(last), key=lambda n: (-last.count(n), n)))}")
print(f"first feature it drops for good: {', '.join(f'{n} {first.count(n)}' for n in sorted(set(first), key=lambda n: (-first.count(n), n)))}")
print(f"test R2: no penalty {r2[:, 0].mean():.3f}, at most 3 features {r2[:, 1].mean():.3f}; lower with 3 in {(r2[:, 1] < r2[:, 0]).sum()} of 50")
print(f"correlation of s1 and s2: {np.corrcoef(X[:, 4], X[:, 5])[0, 1]:.2f}")

last feature lasso keeps, 50 splits: bmi 37, s5 13
first feature it drops for good: age 21, s2 16, s4 9, s1 2, s6 2
test R2: no penalty 0.471, at most 3 features 0.379; lower with 3 in 50 of 50
correlation of s1 and s2: 0.90