Linear regression

Problem

You have examples with numeric features and a numeric target: hours studied and an exam score, square footage and a sale price, ad spend and sign-ups. You want a rule that predicts the target for new examples, and you want to read off how much each feature matters.

Linear regression is the simplest model that does both. It predicts with a weighted sum, ŷ = w·x + b, so each weight says how much the prediction moves per unit of its feature. It has a closed-form solution, so there is no learning rate to tune. And it is the baseline every fancier model should beat: a neural network's last layer is a linear regression on learned features, and many "deep" wins vanish once a linear model gets the same features.

Intuition

Draw a line through a cloud of points. For each point, the vertical gap between the point and the line is its residual, y − ŷ. Least squares picks the line that makes the sum of squared residuals as small as possible. Squaring makes every gap count as positive, punishes big misses hardest, and (as the loss-functions page shows) is exactly what you get by assuming Gaussian noise around the line.

Two facts fall out of that choice. The best line always passes through the point of means (x̄, ȳ). And its slope is the covariance of x and y divided by the variance of x: how much y moves with x, per unit of x's own spread.

Left: the fit and its residuals. Right: the residual plot. Residuals sum to zero and show no pattern against the fitted values, which is what a well-specified linear model leaves behind.

Mechanics

The model and the loss

Stack n examples with d features into a matrix X (n × d) and the targets into a vector y. Append a column of ones to X so the intercept becomes one more weight; then the model is a single matrix-vector product:

ŷ = Xw    L(w) = (1/n)‖y − Xw‖² = (1/n) Σ_i (y_i − x_i·w)²

"Linear" means linear in the weights, not in the raw inputs. Columns of X can be x², log x, one-hot categories or products of features; the solver doesn't care.

The normal equations, derived

L is a quadratic bowl in w, so its minimum is where the gradient is zero. Expand ‖y − Xw‖² = yᵀy − 2wᵀXᵀy + wᵀXᵀXw and differentiate term by term:

∇L = (2/n)(XᵀXw − Xᵀy) = 0  ⟹  XᵀX w = Xᵀy

These are the normal equations: d equations in d unknowns. When the columns of X are linearly independent, XᵀX is invertible and the solution is unique. The name is geometric: rearranged, they say Xᵀ(y − Xw) = 0, so the residual vector is normal (perpendicular) to every column of X. Two consequences you can check on any fit: with an intercept column, the residuals sum to zero; and they are uncorrelated with every feature.

The Hessian, (2/n)XᵀX, is constant and positive semidefinite, so the loss is convex: there are no local minima to get stuck in, only one bowl.

How it is actually solved

Writing w = (XᵀX)⁻¹Xᵀy is fine on paper and bad in floating point. Forming XᵀX squares the condition number of X, so a problem that loses three digits of precision loses six. Real solvers factor X itself: QR, or the SVD X = UΣVᵀ, which gives w = VΣ⁺Uᵀy and handles rank-deficient X by zeroing tiny singular values (the minimum-norm solution).

scikit-learn's LinearRegression.fit does three things. First, _preprocess_data subtracts the column means from X and the mean from y, so the intercept drops out of the problem. Then, for dense X, it calls scipy.linalg.lstsq on the centered data, an SVD-based least-squares solve (sparse X goes to the iterative lsqr instead, and positive=True to non-negative least squares). Finally _set_intercept recovers b = ȳ − x̄·w, the "line passes through the means" fact from above.

Gradient descent on the same problem

When d is huge or the data doesn't fit in memory, you skip the factorization and walk downhill instead: w ← w − η∇L. Because the loss is a fixed quadratic, how fast that walk goes is set entirely by the Hessian's eigenvalues. Gradient descent is stable only for η < 2/λ_max, and with that step size the flat direction, with eigenvalue λ_min, shrinks by only about a factor (1 − 2λ_min/λ_max) per step. The ratio λ_max/λ_min, the condition number, is the number to watch.

For the five students with raw x, the Hessian over (w, b) is [[22, 6], [6, 2]], with eigenvalues 23.66 and 0.338: a condition number of 70. The bowl is a long, narrow valley because slope and intercept trade off against each other: raise the slope and you must lower the intercept to keep the line through the data. Center x first and the cross term vanishes: the Hessian becomes diag(4, 2), condition number 2. That is one reason _preprocess_data centers, and why you standardize features before training anything by gradient descent.

Same data, same loss, same algorithm. With raw x the first step overshoots into the valley and the next 299 creep along it; centering removes the slope-intercept coupling and the path goes almost straight to the minimum.

Ridge and lasso

Least squares trusts the data completely. With many features, few examples, or features that move together, that trust is misplaced: tiny changes in the data swing the weights wildly, and XᵀX can be singular. Regularization adds a penalty on the size of w.

Ridge (L2) adds α‖w‖². The normal equations become (XᵀX + αI)w = Xᵀy, which is invertible for any α > 0. In SVD terms, the component of the solution along each singular direction is scaled by s/(s² + α) instead of 1/s: directions with small singular values, the unstable ones, are shrunk the most. Ridge shrinks every weight toward zero but never exactly to zero, and when two features are correlated it splits the credit between them.

Lasso (L1) adds α‖w‖₁. There is no closed form, because |w| has a kink at zero. Coordinate descent solves it: cycle through the weights, and for each one compute the least-squares value ρ given all the others, then apply the soft threshold w_j = sign(ρ)·max(|ρ| − α, 0). Any weight whose correlation with the remaining residual is below α is set to exactly zero, so lasso selects features.

Both have a probabilistic reading: ridge is the maximum a posteriori estimate under a Gaussian prior on w, lasso under a Laplace prior. Elastic net mixes the two penalties to get lasso's sparsity with ridge's stable handling of correlated groups.

Two scaling details bite in practice. The penalty treats all weights alike, so features must be on comparable scales (standardize first, or the penalty mostly punishes whichever feature has the smallest units). And scikit-learn's two objectives are scaled differently: Ridge minimizes ‖y − Xw‖² + α‖w‖², while Lasso minimizes (1/(2n))‖y − Xw‖² + α‖w‖₁, so the same α means very different things in each.

Worked example

Ridge bends every path smoothly toward zero and pulls the correlated pair together; lasso's paths are piecewise straight and hit zero one feature at a time, weakest first.

Interpreting coefficients

A coefficient is the change in prediction per unit change in its feature, holding the other features fixed. Three cautions follow.

• Units. A raw coefficient's size depends on the feature's units, so compare standardized coefficients (per standard deviation) when asking which feature matters more. To turn hours studied's 9.04 per standard deviation back into points per hour, divide by its standard deviation, about 2.29: roughly 3.95 points per hour.
• Correlation. "Holding the others fixed" is a counterfactual the data may never show. When two features rise together, the fit can trade weight between them at almost no cost in error, so the split 9.04 versus 4.51 above is fragile. Their sum is well determined; the individual values aren't.
• Causation. A coefficient measures association conditional on the included features. Leave out a confounder and the coefficient absorbs its effect.

Reading residual plots

Plot residuals against fitted values (or against each feature) after every fit. A well-specified model leaves a structureless band around zero. A curve means a missing nonlinear term: add x² or a spline. A funnel that widens with ŷ means the noise grows with the prediction (heteroscedasticity): model log y, or use weights. A few points far outside the band are outliers that squared error lets dominate: check them, or switch to a Huber loss. A trend against time or row order means the examples aren't independent.

Implementation

import torch
import torch.nn.functional as F

x = torch.tensor([[1.0], [2.0], [3.0], [4.0], [5.0]])
y = torch.tensor([[52.0], [60.0], [61.0], [75.0], [82.0]])

# closed form: append a column of ones and solve least squares
A = torch.cat([x, torch.ones_like(x)], dim=1)
torch.linalg.lstsq(A, y).solution            # tensor([[ 7.5000], [43.5000]])

# gradient descent: center x first so the bowl is round
model = torch.nn.Linear(1, 1)
opt = torch.optim.SGD(model.parameters(), lr=0.2)
xc = x - x.mean()
for _ in range(100):
    opt.zero_grad()
    loss = F.mse_loss(model(xc), y)
    loss.backward()
    opt.step()
model.weight, model.bias                     # ≈ 7.5 and 66.0 (the intercept at x̄ = 3)

from sklearn.linear_model import Lasso, LinearRegression, Ridge
from sklearn.preprocessing import StandardScaler

X = [[1, 2, 7], [2, 1, 6], [3, 4, 5], [4, 3, 6], [5, 6, 9], [6, 5, 9], [7, 8, 4], [8, 7, 7]]
y = [49, 49, 59, 61, 71, 76, 83, 85]
Z = StandardScaler().fit_transform(X)

LinearRegression().fit(Z, y).coef_           # [9.04, 4.51, 0.60]
Ridge(alpha=10).fit(Z, y).coef_              # [4.26, 3.93, 0.44]
Lasso(alpha=1.0).fit(Z, y).coef_             # [8.84, 3.69, 0.  ]

The pieces in plain Python:

def fit_line(xs, ys):
    n = len(xs)
    mx, my = sum(xs) / n, sum(ys) / n
    slope = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / sum((x - mx) ** 2 for x in xs)
    return slope, my - slope * mx

def soft_threshold(rho, alpha):
    return max(abs(rho) - alpha, 0.0) * (1 if rho > 0 else -1)

def lasso(Z, y, alpha, sweeps=500):
    """Coordinate descent for (1/2n)||y - Zw||^2 + alpha*||w||_1; Z standardized, y centered."""
    n, d = len(Z), len(Z[0])
    w, r = [0.0] * d, list(y)
    for _ in range(sweeps):
        for j in range(d):
            rho = sum(row[j] * ri for row, ri in zip(Z, r)) / n + w[j]
            new = soft_threshold(rho, alpha)
            r = [ri - row[j] * (new - w[j]) for row, ri in zip(Z, r)]
            w[j] = new
    return w

fit_line([1, 2, 3, 4, 5], [52, 60, 61, 75, 82])   # (7.5, 43.5)

Tradeoffs

• Closed form versus iterative. A direct solve costs about nd² + d³ and needs X in memory; it is exact and has no learning rate. Gradient descent costs nd per step and streams data, so it wins when d is large or n doesn't fit, but its speed depends on conditioning and you must pick a step size.
• Ridge versus lasso. Ridge is stable, has a closed form, and keeps every feature. Lasso gives sparse, readable models but picks one feature from a correlated group somewhat arbitrarily and can be unstable across resamples. Elastic net sits between them. Choose α by cross-validation (RidgeCV, LassoCV), never on the training error, which always prefers α = 0.
• Squared error's assumptions. Least squares is optimal when the noise is Gaussian with constant variance. Heavy tails let outliers steer the fit; switch the loss to Huber or absolute error.
• Linearity. The model is linear in its features. Feature engineering (polynomials, interactions, splines) buys flexibility but each added column costs variance; past a point, trees or neural networks learn the features instead.
• Extrapolation. A line keeps going. Predictions far outside the range of the training features are guesses, however confident R² looks.

Check yourself

Connections

References

• Hastie, Tibshirani and Friedman. The Elements of Statistical Learning, 2nd ed., chapter 3. Springer, 2009.
• Hoerl and Kennard. Ridge Regression: Biased Estimation for Nonorthogonal Problems. Technometrics, 1970.
• Tibshirani. Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical Society B, 1996.
• Friedman, Hastie and Tibshirani. Regularization Paths for Generalized Linear Models via Coordinate Descent. Journal of Statistical Software, 2010.
• Zou and Hastie. Regularization and Variable Selection via the Elastic Net. Journal of the Royal Statistical Society B, 2005.
• Trefethen and Bau. Numerical Linear Algebra, lectures 11 and 18 (least squares and its conditioning). SIAM, 1997.
• scikit-learn 1.9.1: _base.py, _ridge.py, _coordinate_descent.py and _cd_fast.pyx.
• PyTorch v2.15.0: torch/linalg/__init__.py.