Regularization

The optimizers note finished the machinery: least squares can be fitted for any number of examples and features. Fitting it is not the problem any more; trusting the fit is. The least-squares note showed two nearly identical features whose coefficients swung between $-9.90$ and $11.88$ from one draw of the noise to the next, and a training $R^2$ that climbed to $0.92$ with 40 useless features while the validation $R^2$ fell below zero. This note changes the objective rather than the optimizer: it adds a penalty that makes large weights expensive. The penalty trades a little bias for a large reduction of variance, has a Bayesian reading as a prior, and in its L1 form selects features by setting weights exactly to zero.

Notation. As in the least-squares note, $X \in \mathbb R^{n \times d}$ is the design matrix and $y \in \mathbb R^n$ the targets. Throughout this note the features are standardized (mean 0, variance 1) and the target is centred, so no intercept is needed: it is $\bar y$ and is never penalized. $\lambda \ge 0$ is the strength of the penalty, $\lVert w\rVert_2^2 = \sum_j w_j^2$ and $\lVert w\rVert_1 = \sum_j \lvert w_j\rvert$. $p = \operatorname{rank} X \le \min(n, d)$ is the rank of the design matrix. Its singular value decomposition is $X = U\Sigma V^\top$, where $U \in \mathbb R^{n \times p}$ and $V \in \mathbb R^{d \times p}$ have orthonormal columns $u_j$ and $v_j$ and $\Sigma$ is diagonal with the singular values $\sigma_1 \ge \dots \ge \sigma_p > 0$ (always with an index, unlike the noise level $\sigma$), so the nonzero eigenvalues of $X^\top X$ are $\sigma_j^2$. $S(z, a) = \operatorname{sign}(z)\max(\lvert z\rvert - a, 0)$ is the soft-thresholding function.

Ridge: a penalty on the size of the weights#

Ridge regression (Hoerl and Kennard, 1970) minimizes the squared error plus the squared length of the weight vector:

$$\hat w_\lambda = \arg\min_w \; \lVert y - Xw\rVert^2 + \lambda\lVert w\rVert_2^2 .$$

The gradient is $2X^\top(Xw - y) + 2\lambda w$, and setting it to zero gives

$$\hat w_\lambda = (X^\top X + \lambda I)^{-1}X^\top y .$$

The matrix $X^\top X + \lambda I$ has eigenvalues $\sigma_j^2 + \lambda$ (and $\lambda$ along the directions with no data), all at least $\lambda$. For any $\lambda > 0$ it is positive definite and invertible: the solution exists and is unique even with collinear features or with more features than examples, the cases where least squares broke down. Scaling matters: the penalty treats a weight of 1 on a feature measured in metres the same as on one measured in millimetres, which is why the features are standardized first, with statistics computed on the training data only.

Worked example. One feature $x = (-1, 0, 1)$ and centred target $y = (-2, 0, 2)$: $x^\top x = 2$ and $x^\top y = 4$. Least squares gives $w = 4/2 = 2$; ridge with $\lambda = 2$ gives $4/(2 + 2) = 1$, half the size. Now duplicate the feature, $X = [x, x]$. Least squares has infinitely many solutions, every $(w_1, w_2)$ with $w_1 + w_2 = 2$. Ridge has one: by symmetry $w_1 = w_2 = a$, and the first equation of $(X^\top X + 2I)w = X^\top y$ reads $4a + 2a = 4$, so $a = 2/3$. The penalty splits the weight equally between identical features, the split with the smallest norm.

On the nearly identical features of the least-squares note (correlation $0.9999$, true weights $(1, 1)$), three fresh draws of the noise give least-squares weights $(4.81, -2.86)$, $(-15.84, 17.95)$ and $(-0.11, 2.09)$, and ridge with $\lambda = 1$ gives $(0.99, 0.95)$, $(0.97, 1.13)$ and $(0.98, 0.99)$ (the Code section).

Which directions it shrinks#

Substitute the SVD. Since $X^\top X = V\Sigma^2V^\top$ and $X^\top y = V\Sigma U^\top y$, the fitted values are

$$X\hat w_\lambda = \sum_{j=1}^p u_j\, \frac{\sigma_j^2}{\sigma_j^2 + \lambda}\, u_j^\top y .$$

Least squares projects $y$ onto each principal direction $u_j$ of the data and keeps the projection in full. Ridge multiplies it by $\sigma_j^2/(\sigma_j^2 + \lambda)$: directions in which the features vary a lot ($\sigma_j^2 \gg \lambda$) are kept almost intact, directions of small variance ($\sigma_j^2 \ll \lambda$) are shrunk almost to zero (Hastie, Tibshirani and Friedman, 2009, section 3.4.1). These are exactly the directions that made least squares unstable, because the variance of the least-squares coefficients along $v_j$ is $\sigma^2/\sigma_j^2$. The sum of the factors, $\operatorname{df}(\lambda) = \sum_j \sigma_j^2/(\sigma_j^2 + \lambda)$, is the trace of the ridge hat matrix and is called the effective degrees of freedom: $p$ at $\lambda = 0$, falling to 0 as $\lambda \to \infty$.

On the diabetes data (442 patients, 10 standardized features), the eigenvalues of $X^\top X$ range from $1778.7$ down to $34.6$ and $3.8$. The smallest belongs to a combination dominated by the blood-serum measurements s1 and s2, whose correlation is $0.90$. With $\lambda = 10$ the largest direction is kept at $0.994$, the two smallest at $0.776$ and $0.275$, and $\operatorname{df} = 8.83$; with $\lambda = 100$, $\operatorname{df} = 6.59$; with $\lambda = 1000$, $2.51$.

The shrinkage introduces bias: $E[\hat w_\lambda] = (X^\top X + \lambda I)^{-1}X^\top X\, w \ne w$. It removes much more variance, which is why the mean squared error can fall, the trade-off of the estimators note that the bias–variance note applies to every model of the course.

The Bayesian reading#

The penalty is a prior. With Gaussian noise of variance $\sigma^2$ and independent priors $w_j \sim N(0, \tau^2)$, the negative log-posterior is the ridge objective divided by $2\sigma^2$, with $\lambda = \sigma^2/\tau^2$; a Laplace prior $p(w_j) \propto e^{-\lvert w_j\rvert/b}$ gives the L1 penalty of the lasso below, with $\lambda = 2\sigma^2/b$. The derivation is in the MAP section of the Bayesian inference note. Two consequences are worth stating here:

• A strong penalty is a confident prior that the weights are small; $\lambda \to 0$ is a flat prior and gives back least squares.
• For the Gaussian prior the posterior is Gaussian, so $\hat w_\lambda$ is both its mode and its mean, and the posterior covariance $\sigma^2(X^\top X + \lambda I)^{-1}$ gives uncertainties. For the Laplace prior the lasso estimate is only the mode; the posterior puts zero probability on any weight being exactly 0, so the sparsity below is a property of the optimization, not of the posterior.

The constrained form: why L1 gives zeros#

For each $\lambda$ there is a radius $t$ such that the penalized problem has the same solution as the constrained one,

$$\min_w \lVert y - Xw\rVert^2 \quad \text{subject to} \quad \lVert w\rVert_2 \le t \quad (\text{ridge}) \qquad \text{or} \quad \lVert w\rVert_1 \le t \quad (\text{lasso}) ,$$

by the Lagrangian duality of convex problems: $\lambda$ is the multiplier of the constraint, and a smaller $t$ corresponds to a larger $\lambda$. The constrained form explains the different behaviour of the two penalties.

Figure

Why the lasso zeroes weights and ridge does not

The contours of the squared error are ellipses around the least-squares solution $\hat w$. The constrained solution is the point of the feasible set on the lowest contour, the point where a growing ellipse first touches the set. The L2 ball is round, so the touching point is wherever the ellipse’s shape puts it, generically with all coordinates nonzero. The L1 ball is a diamond (a cross-polytope in $d$ dimensions) whose corners lie on the axes, and its corners stick out: an ellipse approaching from a wide range of directions hits a corner first, and at a corner some coordinates are exactly zero. In the figure the condition is visible as a computation: at the corner $(0, 1)$ the negative gradient of the loss is $(1.4, 1.6)$, which lies strictly between the normals $(1, 1)$ and $(-1, 1)$ of the two edges that meet there, so no move along the diamond’s boundary lowers the loss.

The lasso#

The lasso (Tibshirani, 1996) is the L1-penalized problem,

$$\hat w = \arg\min_w \; \lVert y - Xw\rVert^2 + \lambda\lVert w\rVert_1 .$$

It is convex but not differentiable where a weight is zero, and it has no closed form in general. With orthonormal features, $X^\top X = I$, it does. Write $z = X^\top y$, the least-squares solution in that case; then $\lVert y - Xw\rVert^2 = \lVert y\rVert^2 - 2w^\top z + \lVert w\rVert^2$, and the objective splits into one problem per coordinate, $\min_{w_j} w_j^2 - 2z_jw_j + \lambda\lvert w_j\rvert$. For $w_j > 0$ the derivative $2w_j - 2z_j + \lambda$ vanishes at $w_j = z_j - \lambda/2$, which is positive only if $z_j > \lambda/2$; symmetrically for $w_j < 0$; otherwise the minimum is at 0. Together:

$$\hat w_j^{\text{lasso}} = S(z_j, \lambda/2), \qquad \hat w_j^{\text{ridge}} = \frac{z_j}{1 + \lambda} .$$

Ridge shrinks every coefficient by the same factor; the lasso subtracts the same amount from every coefficient’s absolute value and sets those smaller than $\lambda/2$ to zero.

Worked example. Orthonormal features with $z = (3, 1, -0.5)$ and $\lambda = 2$. The lasso gives $(S(3, 1), S(1, 1), S(-0.5, 1)) = (2, 0, 0)$: two features are dropped. Ridge gives $(3, 1, -0.5)/3 = (1, 0.333, -0.167)$: all three are kept, smaller. scikit-learn’s Lasso minimizes $\frac{1}{2n}\lVert y - Xw\rVert^2 + \alpha\lVert w\rVert_1$, so its $\alpha$ is $\lambda/(2n)$ in this note’s units; with that conversion it returns $(2, 0, 0)$ (the Code section).

Coordinate descent. For general $X$, minimize over one weight at a time with the others fixed. The one-dimensional problem is the orthonormal one with $z$ replaced by the correlation of feature $j$ with the partial residual $r^{(j)} = y - \sum_{k \ne j}x_k w_k$:

$$w_j \leftarrow \frac{S\big(x_j^\top r^{(j)},\ \lambda/2\big)}{x_j^\top x_j} .$$

Cycling over the coordinates converges for the lasso (Friedman et al., 2007), and computing the solutions for a decreasing sequence of $\lambda$, each started from the previous one, gives the whole regularization path cheaply; this is what scikit-learn does. A related method takes a gradient step on the squared error and then soft-thresholds all coordinates (the proximal gradient method); the least angle regression algorithm computes the path exactly, because it is piecewise linear in $\lambda$ (Efron et al., 2004).

Where the path starts. At $w = 0$ the subgradient condition for optimality is $\lvert 2x_j^\top y\rvert \le \lambda$ for every $j$. So all weights are zero for $\lambda \ge \lambda_{\max} = 2\max_j\lvert x_j^\top y\rvert$, and the first feature to enter is the one most correlated with the target. On the diabetes data $\lambda_{\max} = 39{,}921$, attained by bmi.

Figure

Coefficient paths of ridge and the lasso

The two paths in the figure show the difference on real data. Ridge pulls the opposite coefficients of the correlated pair s1 and s2 ($-37.7$ and $22.7$ in least squares) towards each other and then both towards zero; no coefficient is ever exactly zero. The lasso adds features one at a time as $\lambda$ decreases, in the order bmi, s5, bp, s3, sex, s6, s1, s4, s2, age, and at the value chosen by cross-validation it keeps seven of the ten.

Correlated features: the elastic net#

The lasso’s selection is unstable within a group of correlated features. When two features carry the same signal, the L1 penalty is the same for any split of the weight between them, $\lvert a\rvert + \lvert b\rvert = \lvert a + b\rvert$ for weights of one sign, and the squared error barely distinguishes the splits, so the lasso tends to keep one feature and drop the other, and which one depends on the noise. On two copies of one signal with correlation $0.999$ plus a useless feature, three draws of the data give lasso weights $(1.46, 0, 0)$, $(0, 1.40, 0)$ and $(0.60, 0.76, 0)$ (scikit-learn’s $\alpha = 0.5$).

The elastic net (Zou and Hastie, 2005) adds a ridge term: $\lambda_1\lVert w\rVert_1 + \lambda_2\lVert w\rVert_2^2$. The squared term is strictly convex, so it prefers the equal split, and Theorem 1 of the paper bounds the difference between the weights of two features by a multiple of $\sqrt{1 - \rho}$, where $\rho$ is their correlation: the grouping effect. The same three draws give $(0.76, 0.76, 0)$, $(0.72, 0.74, 0)$ and $(0.71, 0.72, 0)$ (scikit-learn’s $\alpha = 0.5$ and l1_ratio=0.5), and the useless feature is still removed.

Choosing the penalty by validation#

The penalty strength is a hyperparameter, and the validation note gives the procedure: fit for a grid of $\lambda$ on a log scale, estimate the error of each by cross-validation, keep the best, and report the final performance with a nested scheme or a test set, since the best of many cross-validated scores is optimistic. Ridge has a shortcut: its fitted values are linear in $y$, $\hat y = H_\lambda y$, and the leave-one-out residual is $(y_i - \hat y_i)/(1 - h_{ii})$, so leave-one-out cross-validation costs one fit per $\lambda$ (Hastie, Tibshirani and Friedman, 2009, section 7.10.1, eq. 7.51); scikit-learn’s RidgeCV uses it.

How much regularization helps depends on how many features there are per example. On the diabetes data, 5-fold cross-validated $R^2$:

Features	Count	Least squares	Ridge, best $\lambda$	Lasso, best $\alpha$
original	10	$0.489$	$0.490$ ($\lambda = 1$)	$0.490$ ($\alpha = 1$)
with all products of degree 2	65	$0.412$	$0.475$ ($\lambda = 100$)	$0.490$ ($\alpha = 2$)
with all products up to degree 3	285	$-41.7$	$0.441$ ($\lambda = 300$)	$0.470$ ($\alpha = 3$)

With ten features and 442 patients, the variance of least squares is already small and the penalty changes nothing. With 285 features on about 350 training patients per fold, least squares fits the noise so badly that its validation $R^2$ is $-41.7$; ridge and the lasso bring it back to $0.44$ and $0.47$. On this data the added features carry little signal, and the best regularized models only match the plain linear model on the original ten; regularization makes many features safe, not useful.

Other ways of limiting a model act like penalties and appear later in the course: stopping gradient descent early (roughly an L2 penalty with $\lambda \approx 1/(\eta\tau)$ after $\tau$ steps), dropout and data augmentation, all in the training note, and pruning in the decision-tree note. Weight decay, shrinking every weight by a fixed fraction at each step, is the same as an L2 penalty for plain SGD; with adaptive optimizers such as Adam the two differ, and the decoupled form, AdamW, is the one used in practice (weight decay under Adam). Weight decay regularizes the fine-tuning of language models, and an L1 penalty on activations makes the sparse autoencoder sparse.

What changed#

	Least squares	Ridge	Lasso
Objective	$\lVert y - Xw\rVert^2$	$+\,\lambda\lVert w\rVert_2^2$	$+\,\lambda\lVert w\rVert_1$
Solution	$(X^\top X)^{-1}X^\top y$, if it exists	$(X^\top X + \lambda I)^{-1}X^\top y$, always	no closed form; coordinate descent
Collinear features or $d > n$	not unique, unstable	unique, weight shared	unique fit, selection unstable within groups
Zeros	not enforced	not enforced (only by coincidence)	exact zeros: feature selection
Prior	flat	Gaussian	Laplace (as MAP)
Hyperparameter	none	$\lambda$, by cross-validation	$\lambda$, by cross-validation
Bias and variance	unbiased, high variance	biased, lower variance	biased, lower variance

Limitations#

• The penalty depends on the units. Features must be standardized, and the result changes with the scaling chosen.
• Bias. Both penalties shrink large true coefficients too; the lasso’s estimates of the features it keeps are biased towards zero, and refitting least squares on the selected features is a common remedy.
• The lasso keeps at most $n$ features when $d > n$ and the features are in general position, for example drawn from a continuous distribution: the solution is then unique with probability one and has at most $n$ nonzero weights (Tibshirani, 2013). With exactly duplicated columns the fit stays unique but the weights do not: any split of a weight between two copies, with the same sign, has the same fit and penalty. Within strongly correlated groups the choice is arbitrary; the elastic net softens both limits.
• Selection is not inference. A feature dropped by the lasso may still matter, and standard errors computed after selection are too small.
• One more hyperparameter. $\lambda$ costs a cross-validation loop and makes an honest error estimate need a nested one.

What you should remember#

• Ridge: $\hat w_\lambda = (X^\top X + \lambda I)^{-1}X^\top y$, unique for every $\lambda > 0$; it shrinks each principal direction by $\sigma_j^2/(\sigma_j^2 + \lambda)$, the small-variance directions most; $\operatorname{df}(\lambda) = \sum_j \sigma_j^2/(\sigma_j^2 + \lambda)$.
• The penalties are priors: Gaussian gives ridge with $\lambda = \sigma^2/\tau^2$, Laplace gives the lasso as a MAP estimate.
• The L1 ball has corners on the axes, so the lasso sets weights exactly to zero; with orthonormal features it soft-thresholds the least-squares coefficients at $\lambda/2$, while ridge divides them by $1 + \lambda$.
• All lasso weights are zero for $\lambda \ge 2\max_j\lvert x_j^\top y\rvert$; the path is computed by coordinate descent with warm starts.
• The elastic net adds an L2 term so that correlated features share their weight.
• Standardize features, do not penalize the intercept, choose $\lambda$ by cross-validation and report with nested cross-validation.

Code#

The worked examples, ridge on collinear features, the shrinkage factors and the lasso path on diabetes, cross-validation with polynomial features, and the elastic net. The printed numbers come from NumPy 2.1.3 and scikit-learn 1.9.1, and with max_iter=50_000 no lasso fit raises a convergence warning:

import numpy as np, warnings
from sklearn.datasets import load_diabetes
from sklearn.linear_model import Ridge, Lasso, ElasticNet, LinearRegression, lasso_path
from sklearn.model_selection import KFold, cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
warnings.filterwarnings('ignore', category=FutureWarning)   # only scikit-learn 1.9's notice about n_alphas; convergence warnings stay on
rng = np.random.default_rng(0)
ridge = lambda X, y, lam: np.linalg.solve(X.T @ X + lam * np.eye(X.shape[1]), X.T @ y)

# worked example: x = (-1, 0, 1), centred y = (-2, 0, 2); then the same feature twice
x = np.array([[-1.], [0], [1]]); y = np.array([-2., 0, 2])
print(ridge(x, y, 1e-12), ridge(x, y, 2.0), ridge(np.hstack([x, x]), y, 2.0).round(3))   # [2.] [1.] [0.667 0.667]

# orthonormal features: the lasso soft-thresholds z = X^T y at lambda/2, ridge divides it by 1 + lambda
n = 50
Xo = np.linalg.qr(rng.normal(size=(n, 3)))[0]            # orthonormal columns
z = np.array([3.0, 1.0, -0.5]); yo = Xo @ z
lam = 2.0                                                # objective |y - Xw|^2 + lam |w|_1; scikit-learn's alpha = lam / (2n)
print(Lasso(alpha=lam / (2 * n), fit_intercept=False, tol=1e-12).fit(Xo, yo).coef_.round(3) + 0, ridge(Xo, yo, lam).round(3))
# [2. 0. 0.] [ 1.     0.333 -0.167]

# two nearly identical features (the least-squares note's example): ridge with lambda = 1 is stable
x1 = rng.normal(size=100); X2 = np.column_stack([x1, x1 + 0.01 * rng.normal(size=100)])
for _ in range(3):
    t = X2.sum(1) + rng.normal(size=100)
    print(np.linalg.lstsq(X2, t, rcond=None)[0].round(2), ridge(X2, t, 1.0).round(2))
# [ 4.81 -2.86] [0.99 0.95] | [-15.84  17.95] [0.97 1.13] | [-0.11  2.09] [0.98 0.99]

# diabetes, standardized features and centred target: closed form = scikit-learn, shrinkage per principal direction
X, yd = load_diabetes(return_X_y=True)
Xs = (X - X.mean(0)) / X.std(0); yc = yd - yd.mean()
print(np.abs(ridge(Xs, yc, 10.0) - Ridge(alpha=10.0).fit(Xs, yc).coef_).max() < 1e-9)   # True
d2 = np.linalg.eigvalsh(Xs.T @ Xs)[::-1]
print(d2.round(1))
for lam in (1, 10, 100, 1000):
    f = d2 / (d2 + lam)
    print(lam, f[[0, 8, 9]].round(3), f.sum().round(2))  # shrinkage of the largest and two smallest directions; effective df
# [1778.7  659.6  533.   422.3  292.7  266.4  237.2  191.7   34.6    3.8]
# 1 [0.999 0.972 0.791] 9.74 | 10 [0.994 0.776 0.275] 8.83 | 100 [0.947 0.257 0.036] 6.59 | 1000 [0.64  0.033 0.004] 2.51

# the lasso path: the largest useful lambda, and the order in which features enter
names = load_diabetes().feature_names
print(round(2 * np.abs(Xs.T @ yc).max(), 0), names[np.abs(Xs.T @ yc).argmax()])         # 39921.0 bmi
alphas, coefs, _ = lasso_path(Xs, yc, n_alphas=2000, eps=1e-4)
order = []
for k in range(coefs.shape[1]):
    order += [names[j] for j in np.flatnonzero(coefs[:, k]) if names[j] not in order]
print(order)   # ['bmi', 's5', 'bp', 's3', 'sex', 's6', 's1', 's4', 's2', 'age']

# choosing lambda by 5-fold cross-validation, on the 10 features and on all their products up to degree 2 and 3
cv = KFold(5, shuffle=True, random_state=0)
for deg in (1, 2, 3):
    pipe = lambda m: make_pipeline(PolynomialFeatures(deg, include_bias=False), StandardScaler(), m)
    score = lambda m: round(cross_val_score(pipe(m), X, yd, cv=cv, scoring='r2').mean(), 3)
    p = PolynomialFeatures(deg, include_bias=False).fit(X).n_output_features_
    r = {a: score(Ridge(alpha=a)) for a in (1, 10, 100, 300, 1000)}
    l = {a: score(Lasso(alpha=a, max_iter=50_000)) for a in (0.3, 1, 2, 3, 5)}
    print(deg, p, score(LinearRegression()), max(r, key=r.get), max(r.values()), max(l, key=l.get), max(l.values()))
# 1 10 0.489 1 0.49 1 0.49 | 2 65 0.412 100 0.475 2 0.49 | 3 285 -41.699 300 0.441 3 0.47

# two copies of one signal (correlation 0.999) and a useless feature, three draws: the lasso picks a copy, the elastic net shares
for seed in (0, 1, 2):
    r = np.random.default_rng(seed); zz = r.normal(size=200)
    Xe = np.column_stack([zz + 0.03 * r.normal(size=200), zz + 0.03 * r.normal(size=200), r.normal(size=200)])
    Xe = (Xe - Xe.mean(0)) / Xe.std(0); ye = 2 * zz + r.normal(size=200); ye -= ye.mean()
    print(Lasso(alpha=0.5).fit(Xe, ye).coef_.round(2) + 0, ElasticNet(alpha=0.5, l1_ratio=0.5).fit(Xe, ye).coef_.round(2) + 0)
# [1.46 0.   0.  ] [0.76 0.76 0.  ] | [0.  1.4 0. ] [0.72 0.74 0.  ] | [0.6  0.76 0.  ] [0.71 0.72 0.  ]

Next#

The linear regression section ends with a model that predicts a number, is fitted by formula or by descent, and does not chase the noise. Many targets are not numbers but classes: spam or not, default or repayment. The next note asks what a linear model should minimize when the answer is a side of a hyperplane, why squared error is the wrong loss for that, and which losses replace it. The ridge penalty comes along: it is the second term of the support vector machine objective and the heart of kernel ridge regression.

References#

• E. Sokolov et al. Машинное обучение (HSE), lecture notes, lecture 4: regularization, hyperparameters, sparse models and the proximal step, feature transformations.
• A. E. Hoerl, R. W. Kennard. Ridge regression: biased estimation for nonorthogonal problems. Technometrics, 12(1), 1970.
• R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society B, 58(1), 1996.
• H. Zou, T. Hastie. Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society B, 67(2), 2005. Section 2.3, Theorem 1 (the grouping effect).
• B. Efron, T. Hastie, I. Johnstone, R. Tibshirani. Least angle regression. Annals of Statistics, 32(2), 2004.
• J. Friedman, T. Hastie, H. Höfling, R. Tibshirani. Pathwise coordinate optimization. Annals of Applied Statistics, 1(2), 2007.
• R. J. Tibshirani. The lasso problem and uniqueness. Electronic Journal of Statistics, 7, 2013. Uniqueness for features from a continuous distribution and the bound of $n$ nonzero weights.
• T. Hastie, R. Tibshirani, J. Friedman. The Elements of Statistical Learning, 2nd ed. Springer, 2009. Section 3.4.1 (ridge regression, the SVD view and effective degrees of freedom), Section 3.4.2 (the lasso), Section 3.4.3 (the constraint regions and soft thresholding), Section 3.8.6 (pathwise coordinate optimization), Section 7.10.1 (cross-validation and the leave-one-out shortcut).
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Section 3.1.4 (regularized least squares), Section 3.3 (Bayesian linear regression).