Linear Regression and Least Squares

Definition

With ten features instead of one, what is the 'best' line, and can you find it in one step?

Linear regression models the target as a weighted sum of the features plus an intercept. Ordinary least squares chooses the weights that minimise the sum of squared residuals. With the feature matrix X (including a column of ones), the minimiser satisfies the normal equation XᵀXw = Xᵀy, giving w = (XᵀX)⁻¹Xᵀy when XᵀX is invertible.

In simple words

Each prediction is ŷ = w₀ + w₁x₁ + … + wₚxₚ. The error of the whole fit is the sum of squared gaps between predictions and true values: a bowl-shaped function of the weights. The normal equation finds the bottom of the bowl directly, by setting the gradient to zero.

Why it matters

Least squares is the foundation of regression analysis, the default baseline for numeric prediction, and the starting point for ridge, lasso, logistic regression and neural networks. Its coefficients are directly interpretable.

How it works

Write the loss J(w) = ‖y − Xw‖², the squared distance between targets and predictions.
Its gradient is ∇J = −2Xᵀ(y − Xw).
Setting the gradient to zero gives XᵀXw = Xᵀy.
Solve it with a stable method (np.linalg.lstsq or QR) rather than an explicit inverse.

The maths

The model

\hat{y} = X w

$X$: The n × (p + 1) feature matrix, with a first column of ones
$w$: The weights, including the intercept

The least-squares loss

J (w) = ∥ y - X w ∥^{2} = i = 1 \sum n (y_{i} - x_{i}^{⊤} w)^{2}

The normal equation and its solution

X^{⊤} X w = X^{⊤} y \Rightarrow w = (X^{⊤} X)^{- 1} X^{⊤} y

In code

The normal equation by hand, a numerically stable solver, and scikit-learn all give the same weights.

least_squares.py

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(1)
n = 200
X_raw = rng.normal(size=(n, 2))
y = 4.0 + 3.0 * X_raw[:, 0] - 2.0 * X_raw[:, 1] + rng.normal(scale=0.5, size=n)

X = np.column_stack([np.ones(n), X_raw])           # add the intercept column
w_normal = np.linalg.inv(X.T @ X) @ X.T @ y         # textbook formula
w_lstsq = np.linalg.lstsq(X, y, rcond=None)[0]      # what you should use in practice
model = LinearRegression().fit(X_raw, y)

print("True weights:      [4.0, 3.0, -2.0]")
print("Normal equation:  ", np.round(w_normal, 3))
print("lstsq:            ", np.round(w_lstsq, 3))
print("scikit-learn:     ", np.round([model.intercept_, *model.coef_], 3))
residuals = y - X @ w_lstsq
print("Mean squared residual:", round(float(np.mean(residuals ** 2)), 3))

Output

True weights:      [4.0, 3.0, -2.0]
Normal equation:   [ 4.018  3.013 -1.994]
lstsq:             [ 4.018  3.013 -1.994]
scikit-learn:      [ 4.018  3.013 -1.994]
Mean squared residual: 0.259

Try it yourself

Derive the normal equation

Put the steps of the derivation in order.

1Solve for w = (XᵀX)⁻¹Xᵀy
2Differentiate: ∇J = −2Xᵀy + 2XᵀXw
3Write the loss J(w) = (y − Xw)ᵀ(y − Xw)
4Set ∇J = 0 to get XᵀXw = Xᵀy
5Expand: J = yᵀy − 2wᵀXᵀy + wᵀXᵀXw

Common mistakes

Forgetting the column of ones, which forces the line through the origin.
Computing an explicit inverse for large or ill-conditioned problems.
Reading a coefficient as causal. It is an association, holding other features fixed.

Check yourself

1. What does least squares minimise?

2. When does the normal equation have a unique solution?

3. Why prefer np.linalg.lstsq to inv(XᵀX)?

Summary

Linear regression predicts ŷ = Xw; least squares minimises ‖y − Xw‖².
The normal equation XᵀXw = Xᵀy gives the solution in closed form.
Solve it stably, add an intercept column, and watch for collinearity.

Finished reading?

Sources: James et al., An Introduction to Statistical Learning · Hastie et al., The Elements of Statistical Learning