Linear Regression with Normal Equations and Gradient Descent

By Better Wise

This is the first line of code you will write in this course, and I am not going to let you write it yet. Before np.linalg.solve touches your keyboard, the gradient has to be derived on paper, because every later model in this course (logistic regression, softmax, the MLP) is the exact same procedure with a different loss. Recall that you already have the linear algebra (SVD, eigenvalues) and the calculus (chain rule) for this, so the only new work is bookkeeping: get $X^\top(Xw - y)$ right once and you will recognize it everywhere.

By the end of this lesson you can:

• Write the least-squares objective and derive the normal equations by setting the gradient to zero.
• Argue why the normal equations can be unreliable when $X^\top X$ is ill-conditioned or $n \gg d$.
• Implement batch gradient descent in pure NumPy with a correct per-feature gradient.
• State the cost of each method ($O(nd^2 + d^3)$ versus $O(nd)$ per step) and pick the right one for a given problem size.

The Least-Squares Objective and Its Gradient

Model: $\hat{y} = Xw + b$. Fold the intercept into $w$ by appending a column of ones to $X$, so $\hat{y} = Xw$ with $w_0 = b$. Objective:

$$ J(w) = \frac{1}{2n}\lVert Xw - y \rVert_2^2 $$

Expand: $J(w) = \frac{1}{2n}(w^\top X^\top X w - 2 w^\top X^\top y + y^\top y)$. Differentiate with respect to $w$ (the middle term gives two identical contributions, hence the factor 2):

$$ \nabla_w J = \frac{1}{n} X^\top (Xw - y) $$

Read it in English: the residual projected back onto each feature column. That is exactly what gradient descent will compute every step, and it is the same shape as the softmax gradient later. Note the $1/n$: it makes the learning rate independent of dataset size, which matters when you compare runs on different folds.

Normal Equations: Set the Gradient to Zero

Setting $\nabla_w J = 0$ gives $X^\top X w = X^\top y$, so:

$$ \hat{w} = (X^\top X)^{-1} X^\top y $$

Never compute the inverse. Use np.linalg.solve on the $d \times d$ system, or np.linalg.lstsq directly on $X, y$, which handles rank deficiency via the pseudo-inverse.

import numpy as np

def fit_normal_equations(X, y):
    # X already has a ones column for the intercept
    return np.linalg.lstsq(X, y, rcond=None)[0]

def predict(X, w):
    return X @ w

Cost: forming $X^\top X$ is $O(nd^2)$, factorizing it is $O(d^3)$. Cheap when $d$ is small (hundreds), brutal when $d$ is 50,000.

Warning: If two feature columns are nearly collinear, $X^\top X$ is nearly singular. solve will not raise an error, it will return huge, sign-flipping coefficients that fit training data well and generalize terribly. Check np.linalg.cond(X.T @ X); a condition number above roughly $10^{8}$ (in float32 even lower) is your signal that centering, dropping a column, or adding a ridge penalty is required.

Practice

1. You have 500 rows, 3 features, plus intercept. Your colleague has 100,000 rows, 40,000 sparse features. Which one should use the normal equations? Guidance: the first, $d = 4$, so $O(d^3)$ is nothing. The second needs an iterative method, because $X^\top X$ alone is 40,000 squared floats, about 12.8 GB in float64.

2. Derive $\nabla_w J$ for the unnormalized loss $J = \lVert Xw - y \rVert^2$ (no $1/2n$). Answer: $2X^\top(Xw - y)$. The scale factor changes nothing about the optimum, only the learning rate you need.

Batch Gradient Descent

Update rule, applied to the full batch each iteration:

$$ w \leftarrow w - \eta \cdot \frac{1}{n} X^\top (Xw - y) $$

def fit_gd(X, y, lr=0.1, n_iters=1000, tol=1e-10):
    n, d = X.shape
    w = np.zeros(d)
    prev = np.inf
    for it in range(n_iters):
        resid = X @ w - y                 # (n,)
        grad = X.T @ resid / n            # (d,)
        w -= lr * grad
        loss = 0.5 * np.mean(resid ** 2)
        if abs(prev - loss) < tol:
            break
        prev = loss
    return w, it

Three things that decide whether this converges:

The curvature that governs stability is the largest eigenvalue of $\frac{1}{n}X^\top X$: gradient descent converges safely when $\eta < 2/\lambda_{\max}$. That is where your linear algebra pays off: after standardizing, $\lambda_{\max} \le d$ roughly, so $\eta \approx 0.1$ is a sane first guess.

Note: Gradient descent and the normal equations reach the same optimum for a full-rank problem. If your GD loss plateaus well above the closed-form loss, the bug is almost never the objective, it is the learning rate or unscaled features.

Practice

3. Add a standardize(X) helper that returns standardized features plus the mean and std, and a destandardize(w, mu, sigma) that maps fitted weights back to the original scale. Guidance: after centering, the intercept becomes $\bar{y} - \mu^\top w_{\text{scaled}}$; verify by predicting on the raw test set and comparing to the scaled pipeline, the predictions must match to ~1e-10.

4. Benchmark both fitters on sklearn.datasets.make_regression(n_samples=2000, n_features=500) and log wall-clock time and final MSE. Guidance: the normal equations should win on time here, but GD with 200 iterations is competitive and much lighter on memory. This measurement is exactly the comparison you will need for the capstone.

Choosing Between Them

Rule of thumb: under a few thousand features, close the laptop and use lstsq. Above that, or when data arrives in chunks, use gradient descent and standardize first. And always report the loss of both, because they should agree.

Field Notes Before You Write the Code

• Derive before you debug. If you cannot write $\frac{1}{n}X^\top(Xw - y)$ from memory, you will not be able to tell a transposed matrix from a wrong sign when the loss diverges.
• Always fit an intercept by appending a ones column, never by hoping the data is centered.
• Report RMSE, not raw SSE, so numbers stay comparable across dataset sizes and folds.
• Do not skip the closed-form check. Fit both methods once, compare losses, and only then trust your GD loop. This habit becomes the finite-difference check in the next lesson.
• Center before you inspect condition numbers. Centering alone often drops the condition number by orders of magnitude.

What to Take Away

• The gradient is $\frac{1}{n}X^\top(Xw - y)$, the normal equations are its zero, and both methods find the same optimum on full-rank data.
• Use the normal equations (via lstsq) for small $d$; use standardized batch gradient descent for large, sparse, or streaming data.
• Conditioning and learning rate are the two failure modes, and both are diagnosable before any model training happens.

Next up: a finite-difference harness that checks every gradient you ever write, so you will have a machine, not just your memory, confirming the derivations.## Recall Drill for the Two Solvers

Fill these in from memory before the next lesson. If you miss one, re-derive it rather than looking it up.

The gradient of mean squared error with respect to the weights is _____. Setting that gradient to zero yields the normal equations _____. Factorizing the d by d Gram matrix costs _____, while one gradient descent step costs _____. To fit an intercept you append a column of _____ to the design matrix. Calling lstsq avoids forming the inverse of X^T X because it solves via _____. Features should be _____ before gradient descent so one learning rate works for every coordinate.

Front: What does the linear regression model assume about the target? Back: That it is a linear function of the features plus noise: y = Xw + b.

Front: Write the MSE gradient with respect to w. Back: (1/n) X^T (Xw - y).

Front: Write the normal equations. Back: X^T X w = X^T y.

Front: What dominates the cost of the closed-form solution? Back: O(d^3) for the factorization of the d by d Gram matrix, plus O(nd^2) to build it.

Front: What is the cost of one batch gradient descent iteration? Back: O(nd): one pass over the data for the residual and one for the gradient.

Front: Why should you call lstsq instead of inverting X^T X? Back: Forming the inverse squares the condition number and amplifies numerical error; QR or SVD does not.

Front: What does appending a column of ones accomplish? Back: It gives the model an intercept term instead of forcing the fit through the origin.

Front: How does feature scale hurt gradient descent? Back: Unscaled features make the loss surface ill-conditioned, so the step size that is stable on one coordinate explodes on another.

Front: What does the loss curve look like when the learning rate is too large? Back: It oscillates or grows instead of decreasing monotonically.

Front: What does the loss curve look like when the learning rate is too small? Back: It decreases smoothly but far too slowly, and the run hits the iteration cap.

Front: Which metric should you report to compare runs across dataset sizes? Back: RMSE, the square root of the mean squared error, not the raw sum of squares.

Front: When does gradient descent beat the normal equations? Back: When n or d is large, when data arrives in chunks or does not fit in memory, and when you later swap in a regularized or mini-batch loss.

Front: On full-rank data, what should the two methods agree on? Back: The same optimal weights and therefore the same loss, up to numerical tolerance.

Two habits to carry into the rest of the course: keep these cards warm, and keep the closed-form solution as your reference answer whenever a training loop misbehaves.