Linear regression two ways — closed form and gradient descent

Linear regression is worth doing slowly, because it is the only interesting model whose optimum can be written down exactly. Everything after it is the same objective, solved approximately.

The setup

Predict ŷ = Xw, where X is n × d — one row per sample, with a leading column of ones for the intercept — and w is d × 1. Squared error:

J(w) = (1/n) ‖Xw - y‖²

Squared rather than absolute, for two reasons worth keeping separate: it is differentiable everywhere, and its minimiser is the conditional mean. If you wanted the conditional median you would minimise absolute error — a different model, not a worse solver.

Deriving the normal equation

Expand, then differentiate:

J(w) = (1/n)(Xw - y)ᵀ(Xw - y)
     = (1/n)(wᵀXᵀXw - 2wᵀXᵀy + yᵀy)

∂J/∂w = (2/n)(XᵀXw - Xᵀy)

Set the gradient to zero — legitimate here because J is convex, so a stationary point is the global minimum, not a candidate to be checked:

XᵀXw = Xᵀy            the normal equation
w    = (XᵀX)⁻¹Xᵀy     when XᵀX is invertible

By hand on three points

Data: x = 1, 2, 3 with y = 1, 2, 2.

X = [[1, 1],      y = [1,
     [1, 2],           2,
     [1, 3]]           2]

XᵀX = [[3,  6],       Xᵀy = [ 5,
       [6, 14]]              11]

Those entries are only n = 3, Σx = 6, Σx² = 14, Σy = 5, Σxy = 11 — the normal equation is bookkeeping over sums.

det(XᵀX) = 3(14) - 6(6) = 6

(XᵀX)⁻¹ = (1/6) [[14, -6],
                 [-6,  3]]

w = (1/6) [14(5) - 6(11),  -6(5) + 3(11)]
  = (1/6) [4, 3]
  = [0.667, 0.500]

So ŷ = 0.667 + 0.5x. Check the residuals:

x = 1:  ŷ = 1.167   r = -0.167
x = 2:  ŷ = 1.667   r = +0.333
x = 3:  ŷ = 2.167   r = -0.167
                    Σr =  0.000   ✓

Residuals summing to zero is not luck — it is the first normal equation, the one contributed by the intercept column, restated. If your residuals do not sum to roughly zero, you have not fitted an intercept.

Why anyone bothers with gradient descent

	Closed form	Gradient descent
Cost	O(nd² + d³)	O(nd) per step
Memory	forms XᵀX, d × d	one batch at a time
Requires	XᵀX invertible	nothing
Answer	exact, one shot	approximate, many steps

The d³ is the inversion. At d = 100 it is free. At d = 100,000 — text features, one-hot categoricals — it is about 10¹⁵ operations and a d × d matrix that does not fit in memory. That is the crossover, and it is the only reason the exact solution is not always used.

XᵀX also goes singular the moment two features are perfectly collinear: height in cm beside height in inches, or a complete set of one-hot columns beside an intercept. The exact answer then does not exist. The honest fixes are to drop the redundant column, or to add ridge regularisation, which makes (XᵀX + λI) invertible for any λ > 0.

In code, both ways

import numpy as np

X = np.array([[1.0, 1.0], [1.0, 2.0], [1.0, 3.0]])
y = np.array([1.0, 2.0, 2.0])

# Closed form via lstsq, never inv(): stabler, and it copes with a
# rank-deficient X instead of exploding.
w_exact, *_ = np.linalg.lstsq(X, y, rcond=None)

# Gradient descent
w, eta, n = np.zeros(2), 0.1, len(y)
for _ in range(500):
    grad = (2 / n) * X.T @ (X @ w - y)
    w -= eta * grad

print(w_exact, w)   # [0.667 0.5]  [0.667 0.5]

Use lstsq, not inv(X.T @ X) @ X.T @ y. Explicit inversion is slower and loses precision — it squares the condition number of the problem.

What to take forward

• Convex loss means “gradient zero” is the answer.
• The normal equation is that condition, rearranged.
• Residuals sum to zero when an intercept is fitted — a free sanity check.
• Prefer lstsq to inversion; prefer gradient descent once d is large.

---

— Ishaan SandhwarRegressionLinear algebraOptimisation