Linear Regression

Linear regression predicts a real-valued target as a weighted sum of input features. Murphy (§11.1) states the model’s key property as the assumption that the expected output is a linear function of the input, $\mathbb{E}[y \mid \mathbf{x}] = \mathbf{w}^\top \mathbf{x}$ , which makes it easy to fit and to interpret. In computer vision the features may be pixels of a patch, filter responses or a network’s embedding, and the target a brightness, a coordinate offset or a quality score: facial landmark detection describes early fitting algorithms whose parameter updates were learned by linear regression, and super-resolution and correlation filter tracking fit linear maps with a ridge penalty.

The algebra of fitting, the normal equations, QR and SVD solvers, conditioning, weighted least squares and the Gauss–Markov theorem, belongs to Least Squares. This article treats linear regression as a statistical estimator: what ordinary least squares (OLS) estimates and how its fit is judged, how ridge regression shrinks coefficients along the singular directions of the data, how the lasso sets coefficients exactly to zero, how the elastic net combines the two, and how scikit-learn scales each penalty. The example fits all four estimators to texture patches whose target depends on five known pixels.

Problem

Given $N$ training pairs $(\mathbf{x}_n, y_n)$ with $\mathbf{x}_n \in \mathbb{R}^D$ and $y_n \in \mathbb{R}$ , choose weights $\mathbf{w}$ and an offset $w_0$ so that $w_0 + \mathbf{w}^\top \mathbf{x}$ predicts $y$ for new inputs drawn from the same source. Murphy (§11.2.1) writes the probabilistic model as

p(y \mid \mathbf{x}, \boldsymbol{\theta}) = \mathcal{N}\big(y \mid w_0 + \mathbf{w}^\top \mathbf{x}, \sigma^2\big),

with the noise variance $\sigma^2$ as a third parameter. With one input the model is simple linear regression, with several inputs multiple linear regression, and with a vector of outputs multivariate linear regression. Replacing $\mathbf{x}$ by a fixed feature map $\phi(\mathbf{x})$ , such as polynomial terms, keeps the model linear in its parameters even though it is no longer linear in the inputs (§11.2.1).

Two questions decide which estimator to use. The first is whether the data determine $\mathbf{w}$ at all: with fewer examples than features ( $N < D$ ), as in the example below, or with nearly collinear features, many weight vectors fit the training data equally well. The second is what the weights are for: prediction alone, or reading off which features matter.

Inputs and Outputs

Input: a design matrix $X \in \mathbb{R}^{N \times D}$ with one example per row, a target vector $\mathbf{y} \in \mathbb{R}^N$ , and for the penalized estimators a regularization strength ( $\lambda$ or scikit-learn’s alpha) and, for the elastic net, a mixing parameter.
Output: the coefficient vector $\hat{\mathbf{w}}$ and offset $\hat{w}_0$ . Predictions are $\hat{y} = \hat{w}_0 + \hat{\mathbf{w}}^\top \mathbf{x}$ . OLS and ridge return a dense $\hat{\mathbf{w}}$ ; the lasso and elastic net return a sparse one, whose nonzero entries are the selected features.

Intuition

Each coefficient $w_d$ is the change in the predicted output per unit change of feature $d$ with all other features held fixed: Murphy (§11.2.3.3) shows that the least-squares coefficients equal the partial regression coefficients. When two features move together, as neighbouring pixels of a smooth image do, the data barely constrain how weight is split between them, and OLS can put a large positive weight on one and a large negative weight on the other with almost no change to the fit.

Ridge regression adds a cost on the squared size of the weights. It leaves the directions the data determine well almost untouched and pulls poorly determined directions towards zero, so a correlated pair shares its weight. The lasso charges the absolute size of each weight instead. That cost has a corner at zero, so a feature whose correlation with the remaining residual is below a threshold receives exactly zero weight, and the fitted model uses only a subset of the features. The elastic net charges both.

Algorithm

OLS and ridge have closed forms; the lasso and elastic net are solved iteratively. Coordinate descent, the method scikit-learn’s Lasso and ElasticNet use, runs as follows (Murphy, Algorithm 11.1, which calls it the shooting algorithm):

Initialize $\mathbf{w}$ (Murphy starts from the ridge solution; scikit-learn starts from zero, or from the previous solution when warm-starting along a path).
For each coordinate $d = 1, \dots, D$ in turn: compute the correlation $c_d$ of feature $d$ with the residual left by all the other features, and set $w_d$ to the soft-thresholded value of $c_d$ divided by the feature’s squared norm.
Repeat the sweep until the coefficients stop changing.

To fit a whole regularization path, solve for a decreasing sequence of penalties, starting each fit from the previous solution (warm starting, Murphy §11.4.9.4).

Mathematical Formulation

Ordinary least squares as an estimator

Maximizing the Gaussian likelihood over $\mathbf{w}$ minimizes the residual sum of squares; Murphy (Eq. 11.6) defines it with a factor $\tfrac12$ , $\mathrm{RSS}(\mathbf{w}) = \tfrac12\lVert X\mathbf{w} - \mathbf{y} \rVert_2^2$ . Setting its gradient to zero gives the normal equations and, when $X$ has full column rank, the OLS solution $\hat{\mathbf{w}} = (X^\top X)^{-1}X^\top\mathbf{y}$ (Eq. 11.9). Least Squares derives this, shows why solvers factor $X$ instead of inverting $X^\top X$ , and gives the estimate’s covariance $\sigma^2 (X^\top X)^{-1}$ and the Gauss–Markov theorem. The maximum-likelihood estimate of the noise variance is the mean squared residual, $\hat{\sigma}^2 = \frac1N \sum_n (y_n - \hat{\mathbf{w}}^\top\mathbf{x}_n)^2$ (Eq. 11.49).

It follows from the model that OLS is unbiased when $X$ has full column rank: substituting $\mathbf{y} = X\mathbf{w} + \boldsymbol{\epsilon}$ gives $\hat{\mathbf{w}} = \mathbf{w} + (X^\top X)^{-1}X^\top\boldsymbol{\epsilon}$ , whose mean is $\mathbf{w}$ . All the error is variance, and it is large in the directions where $X^\top X$ has small eigenvalues. When $N < D$ the matrix $X^\top X$ is singular and the normal equations have infinitely many solutions; the pseudoinverse picks the one of smallest norm (Singular Value Decomposition). That solution lies in the row space of $X$ , so the component of the true $\mathbf{w}$ outside it is lost however small the noise; Inductive Bias shows the resulting peak of test error near the interpolation threshold.

The offset. Murphy (§11.2.3.1) shows that $\mathbf{w}$ can be fitted on centred data and the offset recovered afterwards as $\hat{w}_0 = \bar{y} - \bar{\mathbf{x}}^\top\hat{\mathbf{w}}$ . scikit-learn does exactly this: with fit_intercept=True, its preprocessing centres dense $X$ and $\mathbf{y}$ , stores the means, and sets intercept_ to y_offset - X_offset @ coef_ (_base.py). With fit_intercept=False the data are taken as already centred. Because the offset is fitted separately, the penalties below never apply to it; Murphy (§11.3) notes that the offset only affects the global mean of the output and does not contribute to overfitting.

A generative reading. Fitting a joint Gaussian to $(\mathbf{x}, y)$ and conditioning on $\mathbf{x}$ gives $\mathbf{w} = \Sigma_{xx}^{-1}\Sigma_{xy}$ , the same coefficients as least squares; Murphy (§11.2.3.5) adds that this equivalence holds only for Gaussian models.

Goodness of fit. Residuals $r_n = y_n - \hat{y}_n$ should look like structureless noise; a curved trend in a residual plot indicates a missing term, and with many inputs one plots $\hat{y}_n$ against $y_n$ instead (§11.2.4.1). The coefficient of determination compares the residual sum of squares with the total sum of squares about the mean,

R^2 = 1 - \frac{\sum_n (y_n - \hat{y}_n)^2}{\sum_n (y_n - \bar{y})^2},

so a model no better than predicting $\bar{y}$ has $R^2 = 0$ and a perfect fit $R^2 = 1$ (Eq. 11.51). The root of the mean squared error is the other usual summary. On training data with $N \le D$ and rows in general position, OLS reaches $R^2 = 1$ whatever the target, which is why fit is judged on held-out data.

Ridge regression: shrinkage along singular directions

Ridge regression is the MAP estimate under a zero-mean Gaussian prior on the weights, and minimizes the penalized sum of squares (Murphy, Eq. 11.55)

J(\mathbf{w}) = \lVert \mathbf{y} - X\mathbf{w} \rVert_2^2 + \lambda \lVert \mathbf{w} \rVert_2^2, \qquad \hat{\mathbf{w}}_{\text{ridge}} = (X^\top X + \lambda I)^{-1} X^\top \mathbf{y},

with $\lambda = \sigma^2/\tau^2$ for a prior variance $\tau^2$ . Least Squares presents this as Tikhonov regularization with $\Gamma = I$ , and Regularization relates it to weight decay. The matrix $X^\top X + \lambda I$ is invertible for any $\lambda > 0$ , so the estimate is unique even when $N < D$ .

The singular value decomposition $X = USV^\top$ shows what the penalty does. Murphy (§11.3.2) writes the training predictions as

\hat{\mathbf{y}} = X\hat{\mathbf{w}}_{\text{ridge}} = \sum_j \mathbf{u}_j \, \frac{\sigma_j^2}{\sigma_j^2 + \lambda} \, \mathbf{u}_j^\top \mathbf{y},

where $\sigma_j$ are the singular values of $X$ , whereas least squares keeps every factor at 1. A direction with $\sigma_j^2 \gg \lambda$ passes almost unchanged; a direction with $\sigma_j^2 \ll \lambda$ is suppressed. The small singular values belong to the directions in which the weights are most uncertain, so ridge shrinks exactly the poorly determined combinations of features towards zero; Murphy calls this shrinkage. The effective number of degrees of freedom

\mathrm{dof}(\lambda) = \sum_j \frac{\sigma_j^2}{\sigma_j^2 + \lambda}

falls from the number of nonzero singular values at $\lambda = 0$ towards 0 as $\lambda \to \infty$ (Eq. 11.71); it is the trace of the hat matrix $X(X^\top X + \lambda I)^{-1}X^\top$ , which the example checks. Murphy contrasts this soft weighting with principal components regression, which keeps the first $K$ directions and discards the rest, and reports that the latter predicts less well.

Two consequences follow from the formula. A ridge coefficient is zero only by coincidence, so ridge never selects features. And since $\hat{\mathbf{w}}_{\text{ridge}} = V \operatorname{diag}\!\big(\sigma_j/(\sigma_j^2 + \lambda)\big) U^\top \mathbf{y}$ is a combination of the right singular vectors, it lies in the row space of $X$ , like the minimum-norm OLS solution: when $N < D$ , ridge cannot recover the part of $\mathbf{w}$ that the training inputs never excite. It trades variance for bias, a trade-off the Bias–Variance Tradeoff article covers. The penalty is not invariant to rescaling a feature, so features are usually standardized first (Murphy §11.2.2.3); Evaluating Representations tabulates which feature transformations ridge is and is not invariant to.

The lasso: soft thresholding and sparsity

The lasso, “least absolute shrinkage and selection operator”, replaces the squared penalty by the $\ell_1$ norm. Murphy (§11.4.1) introduces it as MAP estimation under a Laplace prior, citing Tibshirani (1996), with the objective (Eq. 11.75)

\lVert X\mathbf{w} - \mathbf{y} \rVert_2^2 + \lambda \lVert \mathbf{w} \rVert_1 .

The $\ell_1$ norm is the tightest convex relaxation of the count of nonzero weights (§11.4.1). Written as a constrained problem, $\min \mathrm{NLL}(\mathbf{w})$ subject to $\lVert\mathbf{w}\rVert_1 \le B$ , the feasible set is a cross-polytope whose corners lie on the coordinate axes; Murphy’s geometric argument (§11.4.2) is that the likelihood’s contours tend to touch it first at a corner, especially in high dimensions, and corners are sparse solutions. The $\ell_2$ ball has no corners.

The coordinate-wise solution makes this exact. With all weights but $w_d$ fixed, define $a_d = \sum_n x_{nd}^2$ and $c_d = \sum_n x_{nd}\,(y_n - \mathbf{w}_{-d}^\top \mathbf{x}_{n,-d})$ , the correlation of feature $d$ with the residual of the others. Setting the subgradient to zero gives (Eqs. 11.87–11.89)

\hat{w}_d = \mathrm{SoftThreshold}\!\left(\frac{c_d}{a_d}, \frac{\lambda}{a_d}\right), \qquad \mathrm{SoftThreshold}(x, \delta) = \operatorname{sign}(x)\,\big(\lvert x\rvert - \delta\big)_+ .

A feature whose residual correlation satisfies $\lvert c_d\rvert \le \lambda$ gets exactly zero weight; any other is moved towards zero by a fixed amount. In our reading, Murphy’s Eq. 11.75 and this update differ by a factor of 2 in $\lambda$ : the update follows from the gradient $a_d w_d - c_d$ of the halved sum of squares of Eq. 11.6, while Eq. 11.75 writes the full sum of squares. Constants of this kind are why penalty values do not carry between texts and libraries (see Parameters).

Three properties follow (Murphy §11.4.3–§11.4.6):

Shrinkage bias. Even large coefficients are shrunk by the threshold, so the lasso is biased. A two-stage debiasing fits the lasso to find the support and then refits the selected features by least squares.
A piecewise-linear path. All weights are zero above $\lambda_{\max} = \lVert X^\top\mathbf{y}\rVert_\infty$ (Eq. 11.90). As $\lambda$ decreases, features enter the model, and Murphy states, citing Efron et al. (2004), that the solution path is piecewise linear in $\lambda$ between critical values where the active set changes. With $D > N$ , Murphy states that the lasso can select at most $N$ variables before reaching the complete set.
Selection is not prediction. A method that recovers the true support as $N \to \infty$ is called model-selection consistent. Murphy notes that cross-validation chooses $\lambda$ for predictive accuracy, which tends to keep it small enough to admit irrelevant variables, and reports that Meinshausen and Bühlmann (2006) proved the prediction-optimal $\lambda$ does not give selection consistency. In Murphy’s sparse-signal illustration ( $D = 4096$ , 160 spikes, $N = 1024$ , $\lambda = 0.1\lambda_{\max}$ ), the lasso finds the spikes but with shrunken heights, debiasing recovers them, and the minimum-norm solution is dense.

The elastic net

The elastic net adds both penalties (Murphy, Eq. 11.101):

\lVert \mathbf{y} - X\mathbf{w} \rVert_2^2 + \lambda_2 \lVert \mathbf{w} \rVert_2^2 + \lambda_1 \lVert \mathbf{w} \rVert_1 .

Murphy (§11.4.8), citing Zou and Hastie (2005), states that for $\lambda_2 > 0$ the objective is strictly convex, so it has a unique minimum even when $X$ is not full rank; that strictly convex penalties give a grouping effect, so identical features receive identical weights where the lasso may zero one and keep the other; and that with $D > N$ the elastic net can select more than $N$ variables. The scikit-learn guide puts the practical consequence this way: with several correlated features, the lasso “is likely to pick one of these at random, while elastic-net is likely to pick both”. In coordinate descent the $\ell_2$ term only adds a constant to the denominator of the soft-thresholded update, as the scikit-learn kernel does (_cd_fast.pyx).

Parameters

The same name means different things in different places, because each source scales the loss and penalty differently. scikit-learn 1.7.2’s docstrings give:

Estimator	Objective minimized	Defaults
`LinearRegression`	$\lVert \mathbf{y} - X\mathbf{w} \rVert_2^2$	`fit_intercept=True`; dense input is solved with `scipy.linalg.lstsq`
`Ridge`	$\lVert \mathbf{y} - X\mathbf{w} \rVert_2^2 + \alpha \lVert \mathbf{w} \rVert_2^2$	`alpha=1.0`, `solver='auto'`
`Lasso`	$\frac{1}{2N}\lVert \mathbf{y} - X\mathbf{w} \rVert_2^2 + \alpha \lVert \mathbf{w} \rVert_1$	`alpha=1.0`, `max_iter=1000`, `tol=1e-4`, `selection='cyclic'`
`ElasticNet`	$\frac{1}{2N}\lVert \mathbf{y} - X\mathbf{w} \rVert_2^2 + \alpha\rho\lVert \mathbf{w} \rVert_1 + \frac{\alpha(1-\rho)}{2}\lVert \mathbf{w} \rVert_2^2$	`alpha=1.0`, `l1_ratio` $\rho$ `=0.5`

Penalty strength. Ridge’s alpha has no $1/N$ , while Lasso’s and ElasticNet’s divide the loss by $2N$ , so the two alphas are on different scales and neither is Murphy’s $\lambda$ . Inside the solver, enet_path multiplies alpha by the number of samples before calling the coordinate-descent kernel, which minimizes $\tfrac12\lVert\mathbf{y} - X\mathbf{w}\rVert^2 + \alpha\lVert\mathbf{w}\rVert_1 + \tfrac{\beta}{2}\lVert\mathbf{w}\rVert^2$ . In this scaling the all-zero solution starts at $\alpha_{\max} = \lVert X^\top\mathbf{y}\rVert_\infty / (N\rho)$ . The strength is chosen by cross-validation (RidgeCV, LassoCV, ElasticNetCV) or, for the lasso, by AIC or BIC in LassoLarsIC; Murphy (§11.3.3) describes the empirical-Bayes alternative of maximizing the marginal likelihood; uncertainty estimation treats Bayesian linear regression itself.
l1_ratio. $\rho = 1$ is the lasso and $\rho = 0$ ridge; the docstring states that l1_ratio corresponds to glmnet’s alpha and alpha to glmnet’s lambda, and that l1_ratio <= 0.01 is not reliable unless one supplies one’s own sequence of alphas.
fit_intercept. Centres the data and fits an unpenalized offset (see above). Set it to False only for data that are centred by construction.
Solver choices. Ridge with solver='auto' resolves to cholesky for dense input without positivity constraints. When there are more features than samples, that branch solves the dual system with the $N \times N$ kernel matrix $XX^\top$ and maps back with $X^\top$ ; if that solve fails as singular, the kernel routine warns and uses a least-squares solution instead (_ridge.py). LinearRegression passes lstsq a relative cutoff for small singular values of max(X.shape) times machine epsilon.
tol and max_iter (lasso and elastic net). The coordinate-descent kernel checks the duality gap once the largest coefficient update in a sweep is small relative to the largest coefficient, and stops when the gap is below tol times $\mathbf{y}^\top\mathbf{y}$ ; if max_iter sweeps end without that, it warns with a ConvergenceWarning.
selection. 'cyclic' visits the coordinates in order; 'random' draws a coordinate at random at each step.

Complexity

The scikit-learn guide gives OLS, computed by the SVD of $X$ , a cost of $O(ND^2)$ for $N \ge D$ , and gives ridge the same order. When $D > N$ , Murphy (§11.3.1.2) shows that ridge can be solved in $N$ dimensions by replacing each $\mathbf{x}_i$ with its coordinates in the row space, in $O(DN^2)$ operations; scikit-learn’s dual branch also works in $N$ dimensions, through the kernel matrix. One coordinate-descent sweep for the lasso costs $O(ND)$ with residual updates, or $O(D^2)$ with a precomputed Gram matrix $X^\top X$ as in the example; the number of sweeps varies with the problem and the penalty. Murphy (§11.4.9) also describes projected and proximal gradient methods (ISTA and FISTA) and LARS, which computes the whole lasso path by moving from one change of the active set to the next.

Implementation

Sparse pixel weights from texture patches

The example asks what each estimator recovers when the target depends on a few pixels of an image patch, and how that changes when neighbouring pixels are correlated, as they are in natural images.

Data. Each input is a $16 \times 16$ patch, so $D = 256$ . A patch is cut from the centre of its own $32 \times 32$ field of independent standard normal noise, either left as it is (white) or filtered with a Gaussian of standard deviation 2 pixels (corr, scipy.ndimage.gaussian_filter with its default truncation at 4 standard deviations). The 8-pixel margin is the radius of the truncated filter, so every kept pixel sees the whole kernel and the texture is stationary within the patch. Each filtered field is divided by the norm of the filter’s impulse response, which gives every pixel variance 1; the exact correlation of horizontal neighbours, computed from the same impulse response, is 0.9394. The target is $y = \mathbf{w}^\top\mathbf{x} + \epsilon$ with $\epsilon \sim \mathcal{N}(0, 0.5^2)$ , so the best possible test MSE is 0.25. The true $\mathbf{w}$ has five nonzero entries, $+1, -1, +1, -1, +1$ at (row, column) $(3,3), (3,12), (8,8), (12,4), (12,12)$ . Each seed draws 150 training, 150 validation and 2,000 test patches; the two textures use the same noise draws, so they differ only by the filter. Inputs and targets have zero mean by construction, so every fit omits the offset (fit_intercept=False).

Estimators.

OLS: the minimum-norm solution pinv(X) @ y.
Ridge: computed from the SVD shrinkage formula, with $\lambda$ chosen by validation MSE from 29 values between $10^{-3}$ and $10^4$ (scikit-learn’s scaling).
Lasso and elastic net ( $\rho = 0.5$ ): fitted by a NumPy cyclic coordinate descent in scikit-learn’s scaling. It is warm-started along the path $\alpha_{\max}\,10^{-2i/20}$ , $i = 0, \dots, 20$ , and stops when no coefficient changes by more than $10^{-10}$ in a sweep, or after 20,000 sweeps. The sources leave the choice of $\alpha$ open, so both rules are printed: 0.1max is $\alpha = 0.1\,\alpha_{\max}$ , the choice in Murphy’s illustration, and val is the path value with the lowest validation MSE.
Oracle: least squares on the five true pixels, which knows the support. It is a reference, not a competitor.

Measures. The coefficient error $\lVert\hat{\mathbf{w}} - \mathbf{w}\rVert_2$ and the test MSE on noisy targets (Mean Squared Error) are printed first. The support table counts:

nonzero: coefficients that are not exactly zero;
true hits: how many of the five true pixels are among them;
near: nonzero coefficients within two pixels (Chebyshev distance) of a true pixel, true pixels included. The five near regions cover 123 of the 256 pixels, so a selection spread evenly over the patch would have a near fraction of 0.480;
precision and near frac: true hits and near as fractions of the nonzero count, averaged over seeds. Recall is true hits divided by 5 (Precision and Recall);
exact: the number of seeds on which the nonzero set is exactly the true set.

Each row prints the mean with the minimum and maximum over seeds 0–9.

Checks.

C1: OLS agrees with scikit-learn’s LinearRegression (largest difference at most $10^{-7}$ ) and with $X^\top(XX^\top)^{-1}\mathbf{y}$ (at most $10^{-3}$ , a loose bar for an explicit solve with $XX^\top$ ), and interpolates the training targets (training MSE at most $10^{-20}$ ).
C2: the ridge SVD form agrees with the closed form and with Ridge, and $\mathrm{dof}(\lambda)$ equals the trace of the hat matrix (each to within $10^{-8}$ ); every shrinkage factor lies in $(0, 1)$ .
C3 and C4: the coordinate-descent coefficients agree with Lasso and ElasticNet run at a tolerance of $10^{-12}$ (largest difference at most $10^{-6}$ ), with the same zero pattern, and violate the optimality (KKT) conditions by at most $10^{-6}$ relative to the $\ell_1$ weight.
C5: the pixel variance is within 0.02 of 1 and the neighbour correlation within 0.01 of its exact value.
C6: on noise-free data, OLS with 400 white patches recovers $\mathbf{w}$ to within $10^{-10}$ . The lasso at $10^{-3}\alpha_{\max}$ with 150 white patches selects exactly the true pixels, with the coefficient error of the closed-form on-support solution to within $10^{-6}$ . The elastic net there agrees with scikit-learn and satisfies KKT.

The script uses NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.7.2 with one BLAS thread. It took about two minutes on a shared machine; two runs gave identical output.

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
    os.environ[v] = "1"                          # one BLAS thread: reproducible sums
import sys
import numpy as np
from scipy.ndimage import gaussian_filter
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet

SEEDS = [int(a) for a in sys.argv[1:]] or list(range(10))
P, PAD, D = 16, 8, 256                           # 16x16 patches, 256 pixel features
TRUE = {(3, 3): 1.0, (3, 12): -1.0, (8, 8): 1.0, (12, 4): -1.0, (12, 12): 1.0}
W = np.zeros((P, P))
for rc, v in TRUE.items():
    W[rc] = v
W = W.ravel()
SUPP = np.flatnonzero(W)
NOISE, N_TR, N_VA, N_TE = 0.5, 150, 150, 2000
rows, cols = np.divmod(np.arange(D), P)
NEAR = np.zeros(D, bool)                         # within 2 px (Chebyshev) of a true pixel
for r, c in TRUE:
    NEAR |= np.maximum(abs(rows - r), abs(cols - c)) <= 2


def kernel(s):                                   # the filter's response to a unit impulse
    imp = np.zeros((2 * PAD + 1, 2 * PAD + 1)); imp[PAD, PAD] = 1.0
    return gaussian_filter(imp, sigma=s, mode="reflect", truncate=4.0)


def make(seed, s, n_tr=N_TR, noise=NOISE):
    """Gaussian texture patches with filter sigma s (0 = white), unit pixel variance."""
    rng = np.random.default_rng(seed)
    n = n_tr + N_VA + N_TE
    z = rng.standard_normal((n, P + 2 * PAD, P + 2 * PAD))
    e = rng.standard_normal(n) * noise
    if s == 0:
        X = z[:, PAD:PAD + P, PAD:PAD + P]
    else:
        X = np.stack([gaussian_filter(f, sigma=s, mode="reflect", truncate=4.0) for f in z])
        X = X[:, PAD:PAD + P, PAD:PAD + P] / np.sqrt((kernel(s) ** 2).sum())
    X = X.reshape(n, D)
    y = X @ W + e
    a, b = n_tr, n_tr + N_VA
    return (X[:a], y[:a]), (X[a:b], y[a:b]), (X[b:], y[b:])


def cd(X, y, alpha, rho=1.0, w0=None, tol=1e-10, max_sweeps=20000):
    """Cyclic coordinate descent for (1/2n)||y-Xw||^2 + alpha*rho*||w||_1 + alpha*(1-rho)/2*||w||^2."""
    n, d = X.shape
    G, b = X.T @ X / n, X.T @ y / n              # covariance form
    w = np.zeros(d) if w0 is None else w0.copy()
    g = G @ w
    l1, l2 = alpha * rho, alpha * (1 - rho)
    for sweep in range(max_sweeps):
        dmax = 0.0
        for j in range(d):
            c = b[j] - g[j] + G[j, j] * w[j]     # correlation of pixel j with the partial residual
            new = np.sign(c) * max(abs(c) - l1, 0.0) / (G[j, j] + l2)   # soft thresholding
            if new != w[j]:
                g += (new - w[j]) * G[:, j]
                dmax = max(dmax, abs(new - w[j]))
                w[j] = new
        if dmax < tol:
            return w, True
    return w, False


def kkt(X, y, w, alpha, rho=1.0):
    """Largest violation of the optimality conditions, relative to alpha*rho."""
    grad = X.T @ (y - X @ w) / len(y) - alpha * (1 - rho) * w
    nz, l1 = w != 0, alpha * rho
    v1 = np.abs(grad[nz] - l1 * np.sign(w[nz])).max() if nz.any() else 0.0
    v0 = max(np.abs(grad[~nz]).max() - l1, 0.0)
    return max(v1, v0) / l1


def path(X, y, Xv, yv, rho):
    """Warm-started path alpha_max * 10^(-2i/20), i = 0..20; returns fits at i = 10 and at the best validation MSE."""
    amax = np.abs(X.T @ y).max() / (len(y) * rho)
    w, fits = None, []
    for i in range(21):
        a = amax * 10 ** (-2 * i / 20)
        w, ok = cd(X, y, a, rho, w0=w)
        fits.append((np.mean((Xv @ w - yv) ** 2), a, w.copy(), ok, i))
    return {"0.1max": fits[10], "val": min(fits, key=lambda f: f[0])}


LAMS = 10 ** np.linspace(-3, 4, 29)
res, chk = {}, {k: [] for k in ("C1lr", "C1mn", "C1tr", "C2", "C2dof", "C2fac",
                                "C3", "C3z", "C3kkt", "C4", "C4z", "C4kkt", "C5var", "C5lag",
                                "C6a", "C6b_set", "C6b_err", "C6c", "C6c_kkt")}
info = {k: [] for k in ("lam", "dof", "edge", "noconv")}
lag_exact = {0: 0.0, 2: float((kernel(2)[:, :-1] * kernel(2)[:, 1:]).sum() / (kernel(2) ** 2).sum())}

for seed in SEEDS:
    for s, tex in ((0, "white"), (2, "corr")):
        (X, y), (Xv, yv), (Xt, yt) = make(seed, s)
        A = np.vstack([X, Xv, Xt]).reshape(-1, P, P)  # C5: the generator's statistics
        chk["C5var"].append(abs(A.var() - 1))
        chk["C5lag"].append(abs((A[:, :, :-1] * A[:, :, 1:]).mean() / A.var() - lag_exact[s]))
        fit = {}
        # OLS, minimum norm (n < d: interpolates the training targets)
        w = np.linalg.pinv(X) @ y
        chk["C1lr"].append(np.abs(w - LinearRegression(fit_intercept=False).fit(X, y).coef_).max())
        chk["C1mn"].append(np.abs(w - X.T @ np.linalg.solve(X @ X.T, y)).max())
        chk["C1tr"].append(np.mean((X @ w - y) ** 2))
        fit["OLS", "-"] = w
        # ridge: SVD shrinkage, lambda chosen on the validation set
        U, sv, Vt = np.linalg.svd(X, full_matrices=False)
        UTy = U.T @ y
        lam = min(LAMS, key=lambda l: np.mean((Xv @ (Vt.T @ (sv / (sv**2 + l) * UTy)) - yv) ** 2))
        w = Vt.T @ (sv / (sv**2 + lam) * UTy)
        w_closed = np.linalg.solve(X.T @ X + lam * np.eye(D), X.T @ y)
        w_sk = Ridge(alpha=lam, fit_intercept=False).fit(X, y).coef_
        chk["C2"].append(max(np.abs(w - w_closed).max(), np.abs(w - w_sk).max()))
        fac = sv**2 / (sv**2 + lam)                   # shrinkage factor of each singular direction
        H = X @ np.linalg.solve(X.T @ X + lam * np.eye(D), X.T)
        chk["C2dof"].append(abs(fac.sum() - np.trace(H)))
        chk["C2fac"].append(bool(fac.min() > 0 and fac.max() < 1))
        info["lam"].append((tex, lam)); info["dof"].append((tex, fac.sum()))
        fit["ridge", "val"] = w
        # lasso (rho = 1) and elastic net (rho = 0.5) by NumPy coordinate descent
        for rho, name, c in ((1.0, "lasso", "C3"), (0.5, "enet", "C4")):
            for rule, (_, a, w, ok, i) in path(X, y, Xv, yv, rho).items():
                if rho == 1.0:
                    sk = Lasso(alpha=a, fit_intercept=False, tol=1e-12, max_iter=100000)
                else:
                    sk = ElasticNet(alpha=a, l1_ratio=rho, fit_intercept=False, tol=1e-12, max_iter=100000)
                w_sk = sk.fit(X, y).coef_
                chk[c].append(np.abs(w - w_sk).max())
                chk[c + "z"].append(bool(np.array_equal(w == 0, w_sk == 0)))
                chk[c + "kkt"].append(kkt(X, y, w, a, rho))
                info["edge"].append(rule == "val" and i in (0, 20)); info["noconv"].append(not ok)
                fit[name, rule] = w
        w = np.zeros(D); w[SUPP] = np.linalg.lstsq(X[:, SUPP], y, rcond=None)[0]
        fit["oracle", "-"] = w                        # reference: OLS on the true support
        for key, w in fit.items():
            nz = w != 0
            nnz, tp, near = int(nz.sum()), int(nz[SUPP].sum()), int((nz & NEAR).sum())
            res.setdefault((tex,) + key, []).append(
                (np.linalg.norm(w - W), np.mean((Xt @ w - yt) ** 2), nnz, tp,
                 tp / nnz, near / nnz, int(np.array_equal(np.flatnonzero(nz), SUPP)), near))
    # C6: noise-free base checks
    (X4, y4), _, _ = make(seed, 0, n_tr=400, noise=0.0)
    chk["C6a"].append(np.abs(np.linalg.pinv(X4) @ y4 - W).max())
    (X, y), _, _ = make(seed, 0, noise=0.0)
    G, b = X.T @ X / len(y), X.T @ y / len(y)
    a = np.abs(b).max() * 1e-3
    w, _ = cd(X, y, a, 1.0)
    w_cf = np.zeros(D)                                # closed form on the true support
    w_cf[SUPP] = np.linalg.solve(G[np.ix_(SUPP, SUPP)], b[SUPP] - a * np.sign(W[SUPP]))
    chk["C6b_set"].append(bool(np.array_equal(np.flatnonzero(w), SUPP)))
    chk["C6b_err"].append(abs(np.linalg.norm(w - W) - np.linalg.norm(w_cf - W)))
    a2 = np.abs(b).max() / 0.5 * 1e-3
    w, _ = cd(X, y, a2, 0.5)
    w_sk = ElasticNet(alpha=a2, l1_ratio=0.5, fit_intercept=False, tol=1e-12, max_iter=100000).fit(X, y).coef_
    chk["C6c"].append(np.abs(w - w_sk).max()); chk["C6c_kkt"].append(kkt(X, y, w, a2, 0.5))

print(f"seeds {SEEDS[0]}-{SEEDS[-1]}; mean [min, max] over seeds")
print(f"{'texture':7} {'method':6} {'rule':6} {'coef error':>22} {'test MSE':>25}")
for (tex, m, rule), v in res.items():
    v = np.array(v)
    print(f"{tex:7} {m:6} {rule:6} {v[:,0].mean():7.3f} [{v[:,0].min():.3f}, {v[:,0].max():.3f}]"
          f" {v[:,1].mean():8.4f} [{v[:,1].min():.4f}, {v[:,1].max():.4f}]")
print(f"\n{'texture':7} {'method':6} {'rule':6} {'nonzero':>17} {'true hits':>13} {'near':>17} {'precision':>9} {'near frac':>9} {'exact':>6}")
for (tex, m, rule), v in res.items():
    v = np.array(v)
    print(f"{tex:7} {m:6} {rule:6} {v[:,2].mean():7.2f} [{v[:,2].min():3.0f}, {v[:,2].max():3.0f}]"
          f" {v[:,3].mean():4.2f} [{v[:,3].min():.0f}, {v[:,3].max():.0f}]"
          f" {v[:,7].mean():7.2f} [{v[:,7].min():3.0f}, {v[:,7].max():3.0f}] {v[:,4].mean():9.3f}"
          f" {v[:,5].mean():9.3f} {int(v[:,6].sum()):3d}/{len(v)}")
for tex in ("white", "corr"):
    lam = [l for t, l in info["lam"] if t == tex]; dof = [d for t, d in info["dof"] if t == tex]
    print(f"ridge {tex}: lambda [{min(lam):.3g}, {max(lam):.3g}], dof(lambda) [{min(dof):.1f}, {max(dof):.1f}] of n = {N_TR}")
print(f"validation choice at a grid end: {sum(info['edge'])}; unconverged CD fits: {sum(info['noconv'])}")
print("\nchecks (largest value over seeds and textures)")
print(f"C1 OLS   |pinv-LinearRegression| {max(chk['C1lr']):.1e}  |pinv-X^T(XX^T)^-1 y| {max(chk['C1mn']):.1e}  train MSE {max(chk['C1tr']):.1e}")
print(f"C2 ridge |svd-closed, svd-sklearn| {max(chk['C2']):.1e}  |dof-trace(H)| {max(chk['C2dof']):.1e}  factors in (0,1): {all(chk['C2fac'])}")
for c, n in (("C3", "lasso"), ("C4", "enet ")):
    print(f"{c} {n} |cd-sklearn| {max(chk[c]):.1e}  same zeros: {all(chk[c + 'z'])}  KKT violation {max(chk[c + 'kkt']):.1e}")
print(f"C5 data  |pixel var - 1| {max(chk['C5var']):.4f}  |lag-1 corr - exact| {max(chk['C5lag']):.4f}  (exact: white 0, corr {lag_exact[2]:.4f})")
print(f"C6 noise-free  OLS n=400 max error {max(chk['C6a']):.1e};  lasso exact support: {all(chk['C6b_set'])},"
      f" |error - closed form| {max(chk['C6b_err']):.1e};  enet |cd-sklearn| {max(chk['C6c']):.1e}, KKT {max(chk['C6c_kkt']):.1e}")

Output:

seeds 0-9; mean [min, max] over seeds
texture method rule               coef error                  test MSE
white   OLS    -        1.553 [1.434, 1.672]   2.6691 [2.2848, 3.1547]
white   ridge  val      1.539 [1.420, 1.627]   2.6336 [2.2662, 2.9881]
white   lasso  0.1max   0.294 [0.230, 0.421]   0.3419 [0.3123, 0.4473]
white   lasso  val      0.226 [0.153, 0.320]   0.3062 [0.2760, 0.3640]
white   enet   0.1max   0.533 [0.479, 0.662]   0.5393 [0.4590, 0.7248]
white   enet   val      0.352 [0.288, 0.449]   0.3798 [0.3338, 0.4686]
white   oracle -        0.088 [0.045, 0.149]   0.2610 [0.2509, 0.2768]
corr    OLS    -      1294.613 [875.473, 1858.195]   4.5893 [2.2325, 7.6397]
corr    ridge  val      2.110 [2.088, 2.165]   0.3427 [0.3026, 0.4291]
corr    lasso  0.1max   0.768 [0.542, 1.140]   0.3864 [0.3210, 0.4668]
corr    lasso  val      0.742 [0.494, 1.406]   0.2980 [0.2671, 0.3375]
corr    enet   0.1max   1.862 [1.809, 1.914]   0.4622 [0.3736, 0.5728]
corr    enet   val      1.739 [1.664, 1.826]   0.3109 [0.2784, 0.3631]
corr    oracle -        0.081 [0.043, 0.133]   0.2596 [0.2518, 0.2723]

texture method rule             nonzero     true hits              near precision near frac  exact
white   OLS    -       256.00 [256, 256] 5.00 [5, 5]  123.00 [123, 123]     0.020     0.480   0/10
white   ridge  val     256.00 [256, 256] 5.00 [5, 5]  123.00 [123, 123]     0.020     0.480   0/10
white   lasso  0.1max    8.00 [  6,  11] 5.00 [5, 5]    6.60 [  5,   9]     0.648     0.832   0/10
white   lasso  val      25.70 [ 10,  39] 5.00 [5, 5]   14.50 [  8,  21]     0.228     0.600   0/10
white   enet   0.1max   13.80 [ 10,  20] 5.00 [5, 5]   10.00 [  7,  17]     0.383     0.733   0/10
white   enet   val      55.90 [ 45,  63] 5.00 [5, 5]   29.70 [ 23,  34]     0.090     0.533   0/10
white   oracle -         5.00 [  5,   5] 5.00 [5, 5]    5.00 [  5,   5]     1.000     1.000  10/10
corr    OLS    -       256.00 [256, 256] 5.00 [5, 5]  123.00 [123, 123]     0.020     0.480   0/10
corr    ridge  val     256.00 [256, 256] 5.00 [5, 5]  123.00 [123, 123]     0.020     0.480   0/10
corr    lasso  0.1max   12.20 [ 10,  15] 5.00 [5, 5]   12.00 [ 10,  15]     0.416     0.985   0/10
corr    lasso  val      21.10 [ 17,  26] 5.00 [5, 5]   15.80 [ 13,  18]     0.241     0.757   0/10
corr    enet   0.1max   45.60 [ 42,  50] 5.00 [5, 5]   45.20 [ 42,  49]     0.110     0.992   0/10
corr    enet   val      49.30 [ 42,  62] 5.00 [5, 5]   40.40 [ 35,  44]     0.103     0.824   0/10
corr    oracle -         5.00 [  5,   5] 5.00 [5, 5]    5.00 [  5,   5]     1.000     1.000  10/10
ridge white: lambda [0.001, 31.6], dof(lambda) [121.1, 150.0] of n = 150
ridge corr: lambda [3.16, 31.6], dof(lambda) [34.0, 50.6] of n = 150
validation choice at a grid end: 1; unconverged CD fits: 1

checks (largest value over seeds and textures)
C1 OLS   |pinv-LinearRegression| 9.5e-09  |pinv-X^T(XX^T)^-1 y| 6.1e-05  train MSE 1.2e-23
C2 ridge |svd-closed, svd-sklearn| 6.6e-11  |dof-trace(H)| 5.7e-14  factors in (0,1): True
C3 lasso |cd-sklearn| 7.9e-06  same zeros: True  KKT violation 8.0e-08
C4 enet  |cd-sklearn| 1.0e-08  same zeros: True  KKT violation 1.3e-08
C5 data  |pixel var - 1| 0.0111  |lag-1 corr - exact| 0.0014  (exact: white 0, corr 0.9394)
C6 noise-free  OLS n=400 max error 2.0e-15;  lasso exact support: True, |error - closed form| 1.0e-12;  enet |cd-sklearn| 1.6e-11, KKT 1.0e-08

Checks.

C1, C2, C4, C5 and C6 held.
C3’s coefficient-agreement part failed. The largest difference between the NumPy coordinate descent and Lasso was $7.9\times10^{-6}$ , above the $10^{-6}$ tolerance. Its other two parts held: the zero patterns were identical and the largest KKT violation was $8.0\times10^{-8}$ .
Two flags. The script printed one coordinate-descent fit that hit the 20,000-sweep cap and one validation choice at an end of the $\alpha$ path. The output does not show which fits these are; both are included in the tables.

The supports compared below are identical to scikit-learn’s.

OLS and ridge. On the white texture the two are nearly the same. The coefficient error is 1.553 for OLS and 1.539 for ridge, and the test MSE is about ten times the noise floor. Ridge’s validation choice of $\lambda$ ranged down to $10^{-3}$ , the low end of its grid, which leaves $\mathrm{dof}(\lambda)$ at 150.0 as printed, so ridge is then almost the same as OLS. The construction limits what shrinkage can do here: both solutions lie in the 150-dimensional row space of $X$ , and the true weights’ component outside it is lost to both.

The correlated texture separates them. The minimum-norm OLS solution has coefficient errors around 1,300; in our reading, the smoothing leaves $X$ with many small singular values, which the pseudoinverse divides by, and ridge’s printed $\mathrm{dof}(\lambda)$ of 34–51 out of 150 points the same way. Ridge keeps 34–51 effective degrees of freedom and predicts almost as well as the lasso (test MSE 0.343 against 0.298), yet its coefficient error of 2.110 is the worst of the penalized fits. In our reading, it predicts by spreading weight over the smooth combinations of pixels that the texture actually contains, not by finding the five pixels. Every OLS and ridge coefficient is nonzero, and their near fraction is the uniform 0.480.

Lasso. All five true pixels were selected on every seed, for both textures and both rules. No lasso fit selected only the true pixels (exact 0/10 in every lasso row), so the errors are false positives. The validation rule selected more pixels than 0.1max, as Murphy’s account of prediction-optimal penalties leads one to expect: on the white texture, 25.70 pixels against 8.00. In exchange it predicted better (0.3062 against 0.3419) and had a smaller coefficient error (0.226 against 0.294).

On the correlated texture at 0.1max, the extra pixels were almost all neighbours of true ones (near fraction 0.985). At the validation rule the near fraction fell to 0.757. Even with all its false positives, the lasso’s coefficient error on the white texture stayed far below ridge’s, and at the validation rule its test MSE was within 0.05 of the oracle’s.

Elastic net. On the correlated texture it kept whole neighbourhoods: at 0.1max, 45.60 nonzero pixels of which 45.20 were near true ones, against 12.20 and 12.00 for the lasso. Its test MSE stayed close to the lasso’s at the validation rule (0.3109 against 0.2980), but its coefficient error was more than twice as large (1.739 against 0.742). In our reading, the weight of each true pixel is shared out over its correlated neighbours, which is the grouping effect at work. On the white texture, where there are no correlated groups, the elastic net selected more pixels than the lasso and predicted worse under both rules.

Hypotheses. Each was planned before any run.

The lasso recovers the support on white texture. It finds all true pixels on every seed under both rules: held. The validation rule selects more pixels than the 0.1max rule (25.70 against 8.00): held.
With correlated pixels, precision drops. Under 0.1max it did (0.416 against white’s 0.648): held. Under the validation rule it did not (0.241 against white’s 0.228): failed.
The selected pixels lie near true ones (mean near fraction at least 0.8, a bar the plan fixed without computing what a correct lasso reaches). Held under 0.1max (0.985); failed under the validation rule (0.757).
Test MSE stays low with correlated pixels. The lasso’s 0.2980 at the validation rule is below ridge’s 0.3427 and OLS’s 4.5893, and below 1.5 times the oracle’s 0.2596 (0.3894): held.
Ridge and OLS zero no coefficient. 256 nonzero on every seed and both textures: held.
On correlated pixels, the elastic net selects more pixels, and more near pixels, than the lasso. Nonzero 45.60 against 12.20 (0.1max) and 49.30 against 21.10 (val); near 45.20 against 12.00 and 40.40 against 15.80: held under both rules.
On white texture, coefficient error and test MSE fall from OLS to ridge to the lasso (validation rule). Coefficient error 1.553 > 1.539 > 0.226 and test MSE 2.6691 > 2.6336 > 0.3062: both held.

Process note. The example was planned before it was run, and the plan lists five departures from the batch plan’s description of this example:

The elastic net was added, because this article owns it and the correlated texture is the case that motivates it.
Two $\alpha$ rules and a validation set were added, because the sources leave the choice open.
The oracle row was added as the reachable value and the bar for the test-MSE hypothesis.
The “near” measure was made concrete as Chebyshev distance 2, and correlation length as the filter’s standard deviation.
Noise-free base checks were added.

Three further changes were made:

At planning, after the reachable-value and timing scripts and before the plan’s checksum was recorded, the plan records four changes. Two changed planned items, below; the noise-free lasso check became a single fit at the same $\alpha$ $α$ instead of a path point; and the ridge $\lambda$ $λ$ and $\mathrm{dof}(\lambda)$ $dof (λ)$ ranges and the counts of grid-end choices and unconverged fits were added to the printout, deciding no check.
- The $\alpha$ path was shortened. The plan records the reason: one lasso path on the correlated texture took 29.3 s, with one $\alpha$ “stopped at the 20000-sweep cap (not converged)”. The first grid ran to $10^{-3}\alpha_{\max}$ in 31 steps; the shortened one runs to $10^{-2}\alpha_{\max}$ in 21.
- The noise-free elastic-net check was redefined. As first planned it asked for exactly the true pixels, which no correct method can meet on every seed: on two seeds “the closed-form on-support elastic-net solution violates the off-support optimality condition”. It now asks for agreement with scikit-learn and KKT.
After the pilot, the near-count column was added as a print-only repair. The plan lists that measure, but the pilot script printed only the near fraction. The editor approved the change, the repaired pilot was run once as final, and the run shown here followed with no other change.
The plan’s own time stamps are estimates. The working log gives the true order from file times: plan, reachable-value script, timing script, checksum, pilot, run.

Properties and Behavior

OLS is unbiased when $X$ has full column rank, and the Gauss–Markov theorem makes it the minimum-variance linear unbiased estimator (Least Squares). With collinear features its variance is large. The scikit-learn guide notes that the estimates “rely on the independence of the features”, and that a nearly singular design makes them highly sensitive to random errors in the target. With $N < D$ it interpolates the training data.
Ridge is always unique for $\lambda > 0$ and shrinks the least-determined directions most. It never zeroes a coefficient, and its predictions can be good while its coefficients are not, as on the correlated texture above.
Lasso produces exact zeros and so selects features. It shrinks the coefficients it keeps, and at a prediction-optimal penalty it tends to keep irrelevant features (Murphy §11.4.6). The example’s validation rule selected four to five times as many pixels as are relevant. With correlated features its choice among them is unstable; Murphy (§11.4.8) notes that with identical features the lasso may keep either one.
Elastic net keeps groups of correlated features together and has a unique solution. Over-selection is its price: on the correlated texture its selections were nearly ten times the true support.
Prediction and interpretation are different goals. In the example ridge and the lasso predicted the correlated texture almost equally well with very different coefficients. A low test error says nothing about whether the coefficients identify the features that generate the target.

Limitations

Linearity in the parameters. Nonlinear relations need a feature map chosen in advance; learning the features is the subject of neural networks, and kernel methods use an implicit map.
Squared loss. Squared error penalizes large residuals quadratically, so outliers can give a poor fit (Murphy §11.6). Least Squares describes robust penalties and RANSAC; Murphy (§11.6) describes Laplace and Student-t likelihoods, the Huber loss and RANSAC.
Scale dependence of the penalties. Ridge, lasso and elastic net change with the units of each feature, so inputs are standardized before fitting, and penalty values do not transfer between scalings or libraries.
Selection is fragile with correlated inputs. Pixels of natural images, filter responses at nearby scales and embedding dimensions can be correlated. With identical features the lasso may keep either one (Murphy §11.4.8), and model-selection consistency is a theoretical notion that assumes the data come from the model (Murphy §11.4.6).
Uncertainty. The estimators above return point estimates. Uncertainty in the weights and predictions comes from the sampling distribution of OLS or from Bayesian linear regression (Uncertainty Estimation).

Variants

Weighted least squares for heteroscedastic noise (Murphy §11.2.2.4; Least Squares).
Polynomial and spline regression: linear regression on fixed basis functions (Murphy §11.2.1, §11.5).
Principal components regression: regression on the leading principal components (Principal Component Analysis), which Murphy (§11.3.2) contrasts with ridge’s soft weighting.
Group lasso: penalizes the $\ell_2$ norm of predefined groups of weights, so a whole group enters or leaves together (Murphy §11.4.7). The $\ell_\infty$ variant also tends to give equal magnitudes within a group.
Least angle regression (LARS): computes the lasso path by jumping between changes of the active set; without variable removal it gives a slightly different algorithm, least angle regression or LAR (Murphy §11.4.9.4, citing Efron et al., 2004). scikit-learn provides it as Lars, LassoLars and LassoLarsIC.
Recursive (online) estimation: Murphy (§11.2.3.4) updates the estimate as data arrive for one input; for $D$ -dimensional inputs he points to stochastic gradient descent, which gives the least mean squares algorithm.
Robust regression: Laplace or Student-t likelihoods, the Huber loss and RANSAC (Murphy §11.6).
Bayesian linear regression: a posterior over weights in place of a point estimate (Uncertainty Estimation).
Classification: replacing the Gaussian likelihood by a Bernoulli one gives logistic regression.
Sparse coding elsewhere on the site: super-resolution describes coding patches with an $\ell_1$ -regularized least-squares problem, and motion segmentation regresses each trajectory on the others with a lasso penalty in sparse subspace clustering. With orthonormal features the lasso reduces to soft thresholding each least-squares coefficient (Murphy §11.4.5, Eq. 11.95); thresholding the coefficients of an orthogonal wavelet basis, as in image denoising, is in our reading the image counterpart.