Normal equation closed-form solution

Overview

The normal equation provides an analytical closed-form solution to linear regression by directly computing the optimal parameters without iterative optimization. Unlike gradient descent, it solves for $\boldsymbol{\theta}$ in a single step using matrix calculus.

Think of it like finding the lowest point in a valley: instead of taking many small steps downhill (gradient descent), we use calculus to jump straight to the bottom in one leap.

The Derivation (From First Principles)

Step 1: Set Up the Cost Function

For linear regression with $m$ training examples and $n$ features:

The cost function (Mean Squared Error) is:

Why $\frac{1}{2m}$? The $\frac{1}{2}$ cancels with the derivative of the square (makes calculus cleaner), and $m$ normalizes across different dataset sizes.

Step 2: Vectorize Using Matrix Notation

Let's construct the design matrix $\mathbf{X}$ (dimensions: $m \times (n+1)$):

Why the column of 1's? That's the $x_0 = 1$ feature that multiplies $\theta_0$ (the intercept/bias term).

The parameter vector $\boldsymbol{\theta} \in \mathbb{R}^{(n+1)}$ and target vector $\mathbf{y} \in \mathbb{R}^m$:

Now all predictions in one shot: $\mathbf{X}\boldsymbol{\theta}$ gives a vector of $m$ predictions.

The cost function becomes:

Why this form? The squared error sum $\sum (h - y)^2$ is exactly the dot product of the error vector with itself: $\mathbf{e}^T \mathbf{e}$.

Step 3: Take the Gradient and Set to Zero

To minimize $J$, we need $\nabla_{\boldsymbol{\theta}} J = \mathbf{0}$.

Expanding the quadratic:

Why this expansion? We use $(a-b)^T(a-b) = a^Ta - 2a^Tb + b^Tb$ with $a = \mathbf{X}\boldsymbol{\theta}$ and $b = \mathbf{y}$.

Taking the gradient with respect to $\boldsymbol{\theta}$ (using matrix calculus rules):

Why these terms?

• $\nabla_{\boldsymbol{\theta}} (\boldsymbol{\theta}^T \mathbf{A} \boldsymbol{\theta}) = 2\mathbf{A}\boldsymbol{\theta}$ when $\mathbf{A}$ is symmetric
• $\nabla_{\boldsymbol{\theta}} (\boldsymbol{\theta}^T \mathbf{b}) = \mathbf{b}$
• Constants vanish

Setting the gradient to zero:

Step 4: Solve for $\boldsymbol{\theta}$

Why is it called "normal"? In linear algebra, "normal" means perpendicular. At the optimum, the residual vector $(\mathbf{X}\boldsymbol{\theta} - \mathbf{y})$ is orthogonal (normal) to the column space of $\mathbf{X}$. This is the least squares solution.

Important caveat on the pseudoinverse: When $\mathbf{X}$ has full column rank (i.e., $\mathbf{X}^T \mathbf{X}$ is invertible), the matrix $(\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T$ equals the Moore-Penrose pseudoinverse $\mathbf{X}^+$, giving the unique least squares solution. When $\mathbf{X}$ is rank-deficient (singular $\mathbf{X}^T \mathbf{X}$), this formula breaks down, and the pseudoinverse must be computed via a more general method such as SVD (Singular Value Decomposition). In that case there are infinitely many least squares solutions, and $\mathbf{X}^+$ returns the one with minimum norm.

Worked Examples

Step 1: Construct $\mathbf{X}$ and $\mathbf{y}$

Why this structure? First column is all 1's for $\theta_0$, second column is the $x$ values.

Step 2: Compute $\mathbf{X}^T \mathbf{X}$

Why this multiplication? Each element $(i,j)$ is the dot product of row $i$ from $\mathbf{X}^T$ with column $j$ from $\mathbf{X}$.

Step 3: Compute $\mathbf{X}^T \mathbf{y}$

Step 4: Invert $\mathbf{X}^T \mathbf{X}$

For a 2×2 matrix $\begin{bmatrix} a & b \\ c & d \end{bmatrix}$, the inverse is $\frac{1}{ad-bc}\begin{bmatrix} d & -b \\ -c & a \end{bmatrix}$.

Why compute the determinant? It's $3(14) - 6(6) = 42 - 36 = 6$, which ensures the matrix is invertible.

Step 5: Compute $\boldsymbol{\theta}$

Result: $\theta_0 = 2/3 \approx 0.67$, $\theta_1 = 3/2 = 1.5$

The fitted line is $y = 0.67 + 1.5x$.

Step 1: Design matrix

Why three columns? Column 0 is for $\theta_0$, column 1 for $\theta_1$ (feature $x_1$), column 2 for $\theta_2$ (feature $x_2$).

Step 2: $\mathbf{X}^T \mathbf{X}$

Step 3: $\mathbf{X}^T \mathbf{y}$

Why these numbers? Row 1: $7 + 10 + 13 = 30$. Row 2: $1(7) + 2(10) + 3(13) = 7 + 20 + 39 = 66$. Row 3: $3(7) + 4(10) + 5(13) = 21 + 40 + 65 = 126$.

Step 4: STOP — check invertibility first!

Look carefully at the feature columns: notice that $x_2 = x_1 + 2$ for every row. This means column 2 of $\mathbf{X}$ equals column 1 plus $2 \times$ (column of 1's):

The columns are linearly dependent, so $\mathbf{X}$ is rank-deficient ($\text{rank} = 2 < 3$). Therefore $\mathbf{X}^T \mathbf{X}$ is singular — its determinant is $0$ and it cannot be inverted. The plain normal equation formula does not apply here.

Why this matters: Because the features are redundant, there is no unique $\boldsymbol{\theta}$. There are infinitely many exact solutions.

Step 5: Solve correctly. The relationship is exactly linear ($y = 2 + 2x_1 + x_2$ is one valid fit), so a perfect (zero-error) fit exists. For example:

But because $x_2 = x_1 + 2$, we can trade weight between the parameters. Another exact solution:

wait — let's keep the constraint honest. Any $\boldsymbol{\theta} = (\theta_0, \theta_1, \theta_2)$ satisfying $\theta_0 + 3\theta_1 + \text{(consistent shift)}$ that reproduces the targets works. The clean way to see it: substitute $x_2 = x_1 + 2$ into the model:

Only the combinations $(\theta_0 + 2\theta_2)$ and $(\theta_1 + \theta_2)$ are determined; the individual values are free. The data requires $(\theta_0 + 2\theta_2) = 4$ and $(\theta_1 + \theta_2) = 3$ (matching $y = 4 + 3x_1$), giving infinitely many $\boldsymbol{\theta}$.

The correct tool: Use the SVD-based pseudoinverse ($\mathbf{X}^+$), which returns the minimum-norm solution among these infinitely many, or remove the redundant feature $x_2$ and solve a well-posed 2-parameter problem.

Lesson: Always check the rank of $\mathbf{X}$ before blindly inverting $\mathbf{X}^T \mathbf{X}$.

When to Use (and When Not to Use)

When to use normal equation:

• $n \leq 10{,}000$ features (matrix inversion is tractable)
• Small to medium datasets ($m < 100{,}000$)
• Need exact solution without tuning hyperparameters
• No regularization needed

When to use gradient descent:

• $n > 10{,}000$ (e.g., image features, text features)
• Very large $m$ (big data scenarios)
• Online learning (data arrives in streams)
• Need regularization (Lasso/Ridge) or other constraints

The mistake: Applying the normal equation when $\mathbf{X}^T \mathbf{X}$ is singular (non-invertible)—exactly what happens in Example 2 above.

Why it feels right: The formula looks simple—just plug in your data!

What goes wrong:

• When $n \geq m$ (more features than examples): system is under-determined
• When features are linearly dependent (multicollinearity): $\text{rank}(\mathbf{X}) < n+1$
• Example: including both "price in dollars" and "price in cents" as separate features

The fix:

1. Remove redundant features (check correlation matrix / rank)
2. Use regularization (Ridge regression adds $\lambda \mathbf{I}$ to make $\mathbf{X}^T \mathbf{X} + \lambda \mathbf{I}$ invertible)
3. Use the SVD-based pseudoinverse (numpy's pinv or lstsq handles singular matrices gracefully and returns the minimum-norm solution)
4. Collect more data if $m < n$

The mistake: Computing the normal equation without feature normalization.

Why it feels right: The normal equation is analytical—scaling shouldn't matter mathematically.

What actually happens:

• Numerical instability in matrix inversion when features have vastly different scales
• $\mathbf{X}^T \mathbf{X}$ becomes ill-conditioned (small changes cause huge errors)
• Example: feature 1 is "age" (20-80), feature 2 is "income" (20,000-200,000)

The fix: Standardize features: $x_j := \frac{x_j - \mu_j}{\sigma_j}$ before applying the normal equation.

Implementation Considerations

Instead of explicitly computing $(\mathbf{X}^T \mathbf{X})^{-1}$, use:

1. QR decomposition: $\mathbf{X} = \mathbf{Q}\mathbf{R}$, then solve $\mathbf{R}\boldsymbol{\theta} = \mathbf{Q}^T \mathbf{y}$
2. SVD (Singular Value Decomposition): Handles near-singular and rank-deficient matrices robustly (and defines the true pseudoinverse)
3. Cholesky decomposition: Faster for positive definite $\mathbf{X}^T \mathbf{X}$

Why? Direct matrix inversion is numerically unstable and computationally expensive. These methods avoid inversion while solving the same system.

import numpy as np
 
# Bad: explicit inversion (fails if X^T X is singular)
theta_bad = np.linalg.inv(X.T @ X) @ X.T @ y
 
# Good: use solve (uses LU decomposition internally)
theta_good = np.linalg.solve(X.T @ X, X.T @ y)
 
# Best: use lstsq (uses SVD, handles singular/rank-deficient matrices)
theta_best = np.linalg.lstsq(X, y, rcond=None)[0]

Imagine you're trying to draw the best straight line through some points on a graph. You could guess a line, measure how far each point is from it, adjust a little bit, and keep doing this until you can't improve anymore (that's gradient descent—lots of little adjustments).

But there's a clever shortcut! Using math, you can calculate the PERFECT line in one go without any guessing. It's like having a magic formula that looks at ALL your points at once and figures out exactly where the line should be.

The "normal equation" is that magic formula. You organize all your data points into a special grid (matrix), do some multiplying and dividing with these grids, and out pops the answer—the exact best line!

The catch? This magic formula gets really slow when you have TONS of data or LOTS of different measurements (features). For small problems, it's the fastest way. For huge problems, you're better off with the "guess-and-adjust" method. And watch out: if two of your measurements are secretly the same thing in disguise (like measuring height in inches AND in centimeters), the magic formula gets confused and can't give one answer—there are infinitely many! You have to throw away the duplicate first.

"The Transpose Times, Inverse Times, Transpose Y"

$\boldsymbol{\theta} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}$

Or think: "X-Transpose-X gets inverted, then X-Transpose-y gets inserted"

Visual: $(XtX)^{-1} \cdot Xt \cdot y$ → the "sandwiched inverse" pattern.

Connections

• Linear Regression Fundamentals — the problem this solves
• Gradient Descent — the iterative alternative
• Ridge Regression — adds $\lambda \mathbf{I}$ to handle non-invertibility
• Moore-Penrose Pseudoinverse — generalized inverse for singular matrices (via SVD)
• QR Decomposition — numerically stable way to solve least squares
• Feature Scaling — preprocessing to avoid numerical issues
• Polynomial Regression — still uses normal equation with transformed features
• Underdetermined vs Overdetermined Systems — when unique solutions exist

#flashcards/

Concept Map

Hinglish (regional understanding)

Hinglish (regional understanding)