Why Minimizing MSE Is Not a Guess

This post is a complete reconstruction of linear regression from mathematical first principles — built, tested, and diagnosed from scratch on the Ames Housing Dataset using nothing but NumPy. No sklearn. No black boxes. Just the math, the code, and the reasoning behind every decision.

The Hook

Every ML course teaches you to do this

from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(X, y)
model.predict(X_test)

Three lines.Done.Move on

But here's what those three lines are hiding from you: why does this work? Why do we minimize squared errors and not absolute errors? Not fourth-power errors? Why does taking a derivative and setting it to zero give you the best possible line? And what does "best" even mean when your data is noisy?

Most courses skip these questions. They hand you the formula

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

like it fell from the sky, tell you to memorize it, and move on to cross-validation.

This post does the opposite.

We're going to derive linear regression three times — from probability theory, from matrix calculus, and from geometry — and watch all three roads lead to the exact same formula. Not because that's a cute coincidence, but because it reveals something deep: MSE is not an arbitrary loss function someone chose for convenience. It is the mathematically inevitable consequence of how noise behaves in the real world.

By the end of this post you'll know not just how to use linear regression, but why it works, when it breaks, and what it's actually doing in high-dimensional space. We'll build everything from scratch — no sklearn, no statsmodels, pure NumPy.

If you've ever used fit() without fully knowing what it's optimizing, or why, this post is for you.

The intuition : What is Linear Regression Actually Trying to Do ?

Before any math , lets get the picture right.

You have a dataset — say, house sizes and their sale prices. You plot them. They don't fall on a perfect line, because the world is messy. There's noise: measurement error, missing variables, randomness. But there is a pattern hiding underneath the noise, and you want to find it.

Regression's job is to find the line (or hyperplane, in higher dimensions) that best captures that underlying pattern. Simple enough. But the word "best" is doing a lot of work in that sentence, and how you define it determines everything.

**Why not minimize the absolute error?**You could define "best" as the line that minimizes |yi - ŷi| . This is called Least Absolute Deviations, and it's a perfectly valid idea. It's actually more robust to outliers than squared error. So why don't we use it by default?

Two reasons. First, the absolute value function has a kink at zero — it's not differentiable there, which makes optimization messy. Second, and more importantly, it doesn't connect to any natural probabilistic story about noise the way squared error does. We'll come back to this in the next section and make it precise.

Why not minimize the maximum error? Another idea: make the worst-case prediction error as small as possible. This is minimax regression. The problem is that a single outlier completely controls your solution. One bad data point and your entire model shifts to accommodate it. That's not robustness — that's fragility.

So why squared error? Here's the intuition before the proof: when you measure something in the real world — height, temperature, house price — the errors in your measurements tend to be symmetric and small errors are much more common than large ones. That's the signature of a Gaussian distribution. And as we'll show in the next section, if you assume your noise is Gaussian and ask "what line makes this data most probable?", the answer is exactly the line that minimizes squared error. Not approximately. Exactly.

This is the core insight of the entire post, so let it land:

Minimizing MSE is not a modeling choice. It is the maximum likelihood solution under Gaussian noise. If your noise is Gaussian — and by the Central Limit Theorem, it often approximately is — then OLS is the statistically optimal estimator.

Now let's also preview the other two lenses we'll use, because all three will pay off:

The Algebraic lens treats

$$J(θ) = \frac{1}{2m} | y - X\theta |^2$$

as a function to minimize.

Take the gradient, set it to zero, solve the system. Clean, mechanical, rigorous.

The Geometric lens asks a completely different question: in the vector space of all possible predictions, which prediction is closest to the true yy y? The answer is an orthogonal projection — and the condition for orthogonality turns out to be, again, the exact same formula.

Three different questions. One answer :

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

That's not a coincidence. That's structure. And understanding that structure is what separates someone who uses machine learning from someone who understands it.

Let's build it from the ground up.

Where MSE comes From : The Probabilistic Argument

This is the section most courses skip. They hand you MSE as the loss function and move on. We're going to derive it from scratch, starting from a question that sounds almost philosophical:

What does the noise in your data look like?

The Setup

We model the relationship between inputs and outputs as:

$$y=Xθ+ε$$

Here, y is what you observe, Xθ is the true underlying signal, and ε is everything else — measurement error, missing variables, randomness. We can't eliminate ε, but we can make an assumption about its shape.

The Gaussian Noise Assumption:

$$\varepsilon \sim N(0, \sigma^2)$$

Each error term is drawn independently from a Gaussian with mean zero and variance σ². Mean zero means our model isn't systematically wrong in any direction. The Gaussian shape means small errors are common, large errors are rare.

Is this assumption always true? No. But here's why it's principled and not arbitrary: by the Central Limit Theorem, when many small independent sources of error add up, their sum converges to a Gaussian — regardless of what distribution each individual source follows. In practice, most real-world noise is approximately Gaussian, which is why this assumption is so powerful.

This single assumption is what connects probability theory to least squares. Let's see how.

Step 1 — Write the Likelihood

$$\varepsilon \sim N(0, \sigma^2)$$

Then, conditioned on X and θ, each observation yᵢ is distributed as:

$$y_i \mid x_i, \theta \sim N(x_i^T \theta, \sigma^2)$$

The Gaussian probability density for a single observation is:

$$P(y_i \mid x_i, \theta) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(y_i - x_i^T \theta)^2}{2\sigma^2} \right)$$

Read this carefully. It says: given our parameters θ, the probability of observing yᵢ is highest when xᵢᵀθ is close to yᵢ, and it falls off exponentially as the prediction gets worse. The Gaussian bell curve is literally a picture of how likely each observation is.

Since we assume observations are independent, the joint likelihood of the entire dataset is the product of individual likelihoods:

$$L(\theta) = \prod_{i=1}^m \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(y_i - x_i^T \theta)^2}{2\sigma^2} \right)$$

This is the function we want to maximize over θ. The parameters that maximize this are the ones that make your observed data most probable — that's the entire idea of Maximum Likelihood Estimation.

Step 2 — Take the Log

Products are hard to optimize. Logarithms turn products into sums, and since log is monotonically increasing, maximizing log L(θ) gives the same θ as maximizing L(θ).

$$\log L(\theta) = \sum_{i=1}^m \log\left[ \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(y_i - x_i^T \theta)^2}{2\sigma^2} \right) \right]$$

Split the log across the product:

$$= \sum_{i=1}^m \left[ -\frac{1}{2} \log(2\pi\sigma^2) - \frac{(y_i - x_i^T \theta)^2}{2\sigma^2} \right] $$

$$= -\frac{m}{2} \log(2\pi\sigma^2) - \frac{1}{2\sigma^2} \sum_{i=1}^m (y_i - x_i^T \theta)^2$$

Step 3 — The Punchline

Now look at what we're maximizing over θ. The first term −(m/2)log(2πσ²) is a constant — it doesn't depend on θ at all. So maximizing the log-likelihood reduces to:

$$\arg\max_\theta \log L(\theta) = \arg\max_\theta \left[ -\frac{1}{2\sigma^2} \sum_{i=1}^m (y_i - x_i^T \theta)^2 \right]$$

Maximizing a negative quantity is the same as minimizing the positive version:

$$= \arg\min_\theta \sum_{i=1}^m (y_i - x_i^T \theta)^2$$

That's MSE. That's it. That's the entire justification.

Maximizing the likelihood of your data under Gaussian noise is mathematically identical to minimizing the sum of squared errors. MSE is not a convenient heuristic. It is the MLE estimator under Gaussian noise, derived from first principles.

In matrix notation, this is:

$$\arg\min_\theta | y - X\theta |^2$$

Which is exactly the cost function we'll differentiate in Section 3.

What This Tells You About When MSE Fails

The derivation also tells you exactly when MSE is not the right loss function — whenever the Gaussian noise assumption breaks down:

Heavy-tailed noise (outliers are common) → MSE over-penalizes large errors → use Huber loss or LAD regression
Skewed errors → The symmetry assumption breaks → consider log-transforming the target (which is exactly what we did with SalePrice in the Ames dataset)
Non-constant variance (heteroscedasticity) → σ² is not the same for every observation → use Weighted Least Squares

Understanding the derivation tells you not just when to use MSE, but when to abandon it.

Code: Visualizing the Gaussian Noise Assumption

Let's make this concrete. Here's how to visualize Gaussian noise sitting on top of a linear signal — the picture behind the math:

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# making some data for linear regression demo
sns.set_style("whitegrid") 
np.random.seed(42)

# true values we want to discover
true_b = 2.0  # y-intercept
true_w = 3.5  # slope
noise_level = 1.5

# generate training data
n_points = 80
x = np.linspace(0, 10, n_points)
X = np.c_[np.ones(n_points), x]  # add column of 1s for bias
y_perfect = X @ np.array([true_b, true_w])
y_real = y_perfect + np.random.normal(0, noise_level, n_points)

# create plots
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Plot 1: data scatter with true line
ax1.scatter(x, y_real, c='blue', alpha=0.7, s=50, label='actual data')
ax1.plot(x, y_perfect, 'red', linewidth=3, label='true line y=2+3.5x')

# show some error lines
for i in range(0, n_points, 8):
    ax1.plot([x[i], x[i]], [y_perfect[i], y_real[i]], 'orange', alpha=0.6, linestyle='--')

ax1.set_title('Noisy Data + True Relationship')
ax1.set_xlabel('x values')
ax1.set_ylabel('y values')
ax1.legend()
ax1.grid(True, alpha=0.4)

# Plot 2: errors histogram
errors = y_real - y_perfect
ax2.hist(errors, bins=20, density=True, color='skyblue', edgecolor='white', alpha=0.7)
x_range = np.linspace(-5, 5, 100)
gauss = (1/(noise_level*np.sqrt(2*np.pi)))*np.exp(-x_range**2/(2*noise_level**2))
ax2.plot(x_range, gauss, 'red', linewidth=3, label=f'Gaussian N(0, {noise_level}^2)')

ax2.set_title('Errors are Normally Distributed')
ax2.set_xlabel('error values')
ax2.set_ylabel('frequency')
ax2.legend()
ax2.grid(True, alpha=0.4)

plt.suptitle('Why we assume Gaussian noise in linear regression: y = Xθ + ε, ε~N(0,σ²)')
plt.tight_layout()
plt.savefig('gaussian_demo.png', dpi=120, facecolor='white')
plt.show()

What this plot shows:

Left panel: your noisy observations scattered around the true line, with orange dashed lines showing the individual errors εi\varepsilon_i εi
Right panel: those residuals pile up into a bell curve — the Gaussian distribution we assumed

When you run this on your actual Ames residuals after fitting, the right panel should still approximately hold. If it doesn't, that's your diagnostic telling you the normality assumption is violated.

The Algebraic Argument: Deriving the Normal Equation

We know what we want to minimize. Now let's actually minimize it.

In the previous section, we showed that under Gaussian noise, finding the best parameters θ means solving:

$$\arg\min_\theta | y - X\theta |^2$$

Now we do it rigorously. No hand-waving. We'll take the gradient of this cost function with respect to θ\theta θ, set it to zero, and solve. The result is the Normal Equation — a closed-form, exact solution for linear regression.

Setting Up the Cost Function

Define the cost function as:

$$J(\theta) = \frac{1}{2m} | y - X\theta |^2$$

The ½ is purely cosmetic — it cancels cleanly when we differentiate. The 1/m normalizes by dataset size. Neither affects where the minimum is.

Expand using the definition of squared norm

$$(|v|^2 = v^T v) : $$

$$J(\theta) = \frac{1}{2m} (y - X\theta)^T (y - X\theta)$$

Expand the product:

$$= \frac{1}{2m} \left( y^T y - y^T X \theta - \theta^T X^T y + \theta^T X^T X \theta \right)$$

Since yᵀXθ is a scalar, it equals its own transpose:

$$y^T X \theta = (y^T X \theta)^T = \theta^T X^T y$$

So the two middle terms are equal:

$$J(\theta) = \frac{1}{2m} \left( y^T y - 2\theta^T X^T y + \theta^T X^T X \theta \right)$$

Taking the Gradient

Now differentiate with respect to θ. We need two matrix calculus identities:

$$\frac{\partial}{\partial \theta} (\theta^T a) = a$$

and

$$\frac{\partial}{\partial \theta} (\theta^T A \theta) = 2A \theta \\$$

(when A is symmetric)

Note that XᵀX is always symmetric:

$$(X^T X)^T = X^T X$$

Applying these:

$$\nabla_\theta J(\theta) = \frac{1}{2m} \left( 0 - 2X^T y + 2 X^T X \theta \right) $$

$$= \frac{1}{m} (X^T X \theta - X^T y)$$

Setting the Gradient to Zero

At the minimum, the gradient vanishes:

$$\nabla_\theta J(\theta) = 0 \implies \frac{1}{m}(X^T X \theta - X^T y) = 0 \implies X^T X \theta = X^T y$$

These are called the Normal Equations. They're a system of linear equations in θ. If XᵀX is invertible, we can solve directly:

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

That's it. One formula. Derived from scratch via matrix calculus.

When Does This Fail?

The Normal Equation requires XᵀX to be invertible. It fails when:

1. Multicollinearity — two or more features are linearly dependent (e.g., you include both "area in sqft" and "area in sqm"). The columns of X are not linearly independent, so XᵀX becomes singular and non-invertible. This is exactly why we need Ridge Regression — adding λI to XᵀX guarantees invertibility:

$${\theta_{\text{ridge}} = (X^T X + \lambda I)^{-1} X^T y}$$

2. More features than observations (p > m) — XᵀX is rank-deficient by construction. The system is underdetermined and has infinitely many solutions.

3. Computational cost — Inverting an n × n matrix costs O(n³). For datasets with thousands of features this becomes prohibitively slow. Gradient descent scales far better in practice.

Code: OLS from Scratch

Here is a complete, clean implementation of the Normal Equation in pure NumPy — no sklearn, no statsmodels:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# My simple OLS solver using normal equation
class OLS:
    def __init__(self):
        self.theta = None
        
    def fit(self, X, y):
        # θ = (XᵀX)⁻¹Xᵀy
        self.theta = np.linalg.inv(X.T @ X) @ X.T @ y
        return self
    
    def predict(self, X):
        return X @ self.theta
    
    def r2_score(self, X, y):
        yhat = self.predict(X)
        ss_res = np.sum((y - yhat)**2)
        ss_tot = np.sum((y - np.mean(y))**2)
        return 1 - ss_res/ss_tot

    def rmse(self, X, y):
        return np.sqrt(np.mean((y - self.predict(X))**2))

# Generate data y = 3 + 2x1 - x2 + noise
np.random.seed(42)
m = 100
X_raw = np.random.randn(m, 2)
X = np.column_stack([np.ones(m), X_raw])
true_theta = np.array([3.0, 2.0, -1.0])
y = X @ true_theta + np.random.normal(0, 0.8, m)

# Fit and show results
model = OLS()
model.fit(X, y)

print("Results:")
print(f"{'Feature':<15} {'Coef':>8} {'R2':>6}")
print("-"*30)
print(f"bias       {model.theta[0]:>7.3f}")
print(f"x1         {model.theta[1]:>7.3f}")
print(f"x2         {model.theta[2]:>7.3f}")
print(f"R² = {model.r2_score(X,y):.3f}")
print(f"RMSE = {model.rmse(X,y):.3f}")

print(f"\nTrue θ: {true_theta}")
print(f"Found θ: {np.round(model.theta, 3)}")

Running this you'll see the recovered θ is extremely close to the true parameters — that's OLS working exactly as the math predicts.

Visualizing the Cost Surface

The Normal Equation finds the minimum analytically in one step. But it's worth visualizing what that minimum looks like — especially because Section 5's gradient descent will navigate this same surface iteratively:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# make some 1D data to see cost surface
np.random.seed(0)
m = 60
x = np.random.uniform(0, 10, m)
X = np.column_stack([np.ones(m), x])  # add bias column
y = 2.5 + 1.8 * x + np.random.normal(0, 1.2, m)

# grid of bias and slope values to test
bias_vals = np.linspace(0, 5, 100)
slope_vals = np.linspace(0, 4, 100)
B, S = np.meshgrid(bias_vals, slope_vals)

# calculate MSE for every (bias, slope) combo
costs = np.zeros_like(B)
for i in range(B.shape[0]):
    for j in range(B.shape[1]):
        y_pred = B[i,j] + S[i,j] * x
        costs[i,j] = np.mean((y - y_pred)**2)

# find OLS solution  
ols_bias = np.linalg.inv(X.T @ X) @ X.T @ y
model = lambda b,s: b + s * x

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# left: cost contour plot
contour = ax1.contourf(B, S, costs, 30, cmap='viridis')
ax1.contour(B, S, costs, colors='white', alpha=0.3)
ax1.scatter(ols_bias[0], ols_bias[1], c='red', s=100, 
           label=f'Best point\n({ols_bias[0]:.1f}, {ols_bias[1]:.1f})')
ax1.set_xlabel('bias')
ax1.set_ylabel('slope')
ax1.set_title('Cost surface (MSE)')
ax1.legend()
plt.colorbar(contour, ax=ax1)

# right: data + fitted line
ax2.scatter(x, y, alpha=0.6, s=40)
x_line = np.linspace(0, 10, 100)
y_line = ols_bias[0] + ols_bias[1] * x_line
ax2.plot(x_line, y_line, 'red', linewidth=3)
ax2.set_title('OLS fit line')
ax2.set_xlabel('x')
ax2.set_ylabel('y')

plt.tight_layout()
plt.savefig('cost_surface.png', dpi=120)
plt.show()

print(f"OLS found: bias={ols_bias[0]:.2f}, slope={ols_bias[1]:.2f}")

What this shows: The cost surface is a convex bowl — there are no local minima, only one global minimum. The teal dot is exactly where the Normal Equation lands in a single step. In Section 5, you'll watch gradient descent roll down this bowl step by step.

The Geometric Argument : Regression as Orthogonal Projection

Same answer. Completely different question.

We've derived θ = (XᵀX)⁻¹Xᵀy twice now — once from probability, once from calculus. Here we derive it a third time, asking a question that has nothing to do with optimization or likelihood:

In the space of all possible predictions, which one is closest to the truth?

The answer to that geometric question will hand us the Normal Equation again — not as a coincidence, but as an inevitability.

Forget data points for a moment. Think of y as a single vector in ℝᵐ — one dimension per observation. If you have 100 training examples, y is a point in 100-dimensional space.

Now think about what Xθ can be. As θ varies over all possible parameter vectors, Xθ traces out every possible linear combination of the columns of X. This set of all reachable predictions is called the column space of X, written col(X).

$$\operatorname{col}(X) = { X\theta : \theta \in \mathbb{R}^n }$$

If X has n features, col(X) is an n-dimensional subspace sitting inside ℝᵐ. Unless your model is perfectly specified, y almost certainly does not lie in this subspace — because y contains noise that no linear combination of your features can reproduce exactly.

So we can't reach y. But we can get as close as possible.

The Projection Principle

The Question : Which point ŷ ∈ col(X) is closest to y?

The Answer: The orthogonal projection of y onto col(X).

Here's the key geometric fact: the closest point in a subspace to any vector yy y is always the orthogonal projection. Why? Because the shortest path from a point to a plane is always perpendicular to that plane. Any other path is longer.

This means the residual vector e = y − ŷ must be perpendicular to the entire column space. Not just perpendicular to one direction — perpendicular to every vector in col(X).

$$x_j^T (y - X\theta) = 0 \quad \text{for every column } x_j \text{ of } X$$

Stacking all these conditions together in matrix form:

$$X^T (y - X\theta) = 0$$

Expand:

$$X^T y - X^T X \theta = 0 \implies X^T X \theta = X^T y \implies \theta = (X^T X)^{-1} X^T y $$

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

The Normal Equation. Again. Derived purely from geometry, with no calculus, no probability. Just the requirement that the residual be orthogonal to the column space.

The Diagram

This is the picture worth having in your head every time you run a regression:

y lives off the subspace — it has noise components that no linear model can capture
ŷ = Xθ is the closest point to y that lives in col(X)
e = y − ŷ is the residual — perpendicular to col(X) by construction
The right angle between e and ŷ is not an approximation — it is an exact geometric property of OLS

Three Properties That Follow Immediately

Once you see regression as projection, several facts that seemed like algebraic accidents become geometrically obvious:

Property 1 — Residuals are orthogonal to predictions:

$$\hat{y}^T e = (X\theta)^T (y - X\theta) = \theta^T X^T (y - X\theta) = \theta^T \cdot 0 = 0$$

The fitted values and residuals are always perpendicular. Always. This isn't an assumption — it's a theorem that follows from the projection.

Property 2 — Residuals sum to zero (when a bias term is included):

The column of ones in X is one of the vectors that e must be orthogonal to. So:

$$\mathbf{1}^T e = \sum_{i=1}^m e_i = 0$$

The positive and negative residuals always cancel exactly. If your residuals don't sum to zero, you forgot the bias term.

Property 3 — The Pythagorean decomposition of variance:

Since ŷ and e are orthogonal, the Pythagorean theorem applies in ℝᵐ:

$$|y|^2 = |\hat{y}|^2 + |e|^2$$

Translated into statistics, this is exactly the decomposition:

$$SS_{\text{total}} = SS_{\text{regression}} + SS_{\text{residual}}$$

And

$$R^2 = \frac{|\hat{y} - \bar{y}|^2}{|y - \bar{y}|^2} = \cos^2 \theta $$

$$R^2 = \cos^2 \theta \quad \text{where } \theta \text{ is angle between } (y - \bar{y}) \text{ and } (\hat{y} - \bar{y})$$

What Multicollinearity Looks Like Geometrically

If two columns of X are nearly identical — say feature x₁ ≈ x₂ — then the column space of X is nearly the same whether or not you include both. You're trying to define a 2D subspace using two vectors that point in almost the same direction. The subspace is still well-defined, but the decomposition of ŷ into contributions from each column becomes wildly unstable — a tiny perturbation in y causes enormous swings in θ₁ and θ₂, even though ŷ barely moves.

This is why multicollinearity inflates standard errors without necessarily hurting predictions. Geometrically: the projection ŷ is stable, but the coordinate representation θ along nearly-parallel axes is not.

Ridge Regression fixes this by effectively rotating those near-parallel axes apart — the λI term adds a small penalty that forces the basis to be better conditioned.

Code: Verifying Orthogonality from Scratch

Let's verify every geometric property we just derived — numerically, from scratch:

import numpy as np

np.random.seed(7)
m = 80  # rows
n = 3   # features + bias

# Make fake data
X_raw = np.random.randn(m, n-1)
X = np.column_stack([np.ones(m), X_raw])  # add bias column
true_theta = np.array([1.5, -2.0, 3.0])
y = X @ true_theta + np.random.normal(0, 1.0, m)

# OLS by hand (normal equation)
theta_hat = np.linalg.inv(X.T @ X) @ (X.T @ y)
y_hat = X @ theta_hat    # predictions
residuals = y - y_hat    # errors

print("=" * 50)
print("Checking OLS geometry properties")
print("=" * 50)

# 1. Residuals perpendicular to X columns?
x_t_e = X.T @ residuals
print("\n1. X.T @ residuals (should be all ~0):")
print(f"   biggest value: {np.max(np.abs(x_t_e)):.2e} ✓")

# 2. Residuals sum to 0?
sum_e = np.sum(residuals)
print("\n2. Sum of residuals:")
print(f"   total = {sum_e:.2e} ✓")

# 3. Fitted values perpendicular to residuals?
yhat_dot_e = np.dot(y_hat, residuals)
print("\n3. y_hat • residuals:")
print(f"   dot product = {yhat_dot_e:.2e} ✓")

# 4. Pythagoras: total^2 = explained^2 + residual^2
y_mean = np.mean(y)
y_centered = y - y_mean
yhat_centered = y_hat - np.mean(y_hat)

total_var = np.sum(y_centered**2)
explained_var = np.sum(yhat_centered**2) 
residual_var = np.sum(residuals**2)

print("\n4. Variance decomposition:")
print(f"   Total SS     = {total_var:.3f}")
print(f"   Explained SS = {explained_var:.3f}")
print(f"   Residual SS  = {residual_var:.3f}")
print(f"   Check: {explained_var + residual_var:.3f} ✓")

# 5. R2 both ways
r2_geometry = explained_var / total_var
r2_formula = 1 - (residual_var / total_var)
print(f"\n5. R² geometry = {r2_geometry:.4f}")
print(f"   R² formula  = {r2_formula:.4f}")
print("   SAME! R² is geometric! ✓")
print("=" * 50)

Run this and you'll see every residual is orthogonal to the column space to machine precision (~10⁻¹⁴). These aren't approximations — they're exact consequences of how OLS works.

Why This Perspective Matters

The probabilistic view told you why MSE is the right objective. The algebraic view told you how to solve it. The geometric view tells you what the solution actually is — a projection.

This matters when things go wrong. When your model performs poorly, you now have three diagnostic lenses:

Probabilistically: Is the Gaussian noise assumption violated? Check your residual distribution.
Algebraically: Is XᵀX ill-conditioned? Check your condition number, apply Ridge.
Geometrically: Is y nearly orthogonal to col(X)? Your features simply don't span the relevant directions — you need better features, not a better algorithm.

The geometry also previews everything that comes next in ML. Kernel methods extend the column space into infinite dimensions. Neural networks learn the column space itself. Principal Component Analysis finds the directions in col(X) that capture the most variance. It's projections all the way down.

Two Ways to Solve It: Normal Equation vs Gradient Descent

Same destination. Very different journeys.

We have the formula. Now the engineering question: how do you actually compute θ = (XᵀX)⁻¹Xᵀy in practice? There are two fundamentally different approaches, and choosing between them is one of the first real decisions you make as a practitioner.

Method 1 — The Normal Equation (Closed-Form)

The Normal Equation computes θ\theta θ in a single step. No iterations, no hyperparameters, no convergence to worry about. You invert a matrix and you're done.

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

What's actually happening computationally: Modern implementations don't literally invert XᵀX — matrix inversion is numerically unstable. Instead they use QR decomposition or Cholesky factorization under the hood (np.linalg.solve does this automatically). But conceptually it's the same idea: solve the linear system XᵀXθ = Xᵀy directly.

When it's the right choice:

Small to medium datasets (m < 100,000, n < 10,000)
You need an exact solution in one shot
You want standard errors and t-statistics (these fall out naturally from (XᵀX)⁻¹)
You're doing statistical inference, not just prediction

When it breaks down:

n is large — inverting an n × n matrix costs O(n³). At n = 10,000 features that's 10¹² operations. Intractable.
XᵀX is singular or near-singular — multicollinearity makes inversion numerically unstable
Streaming or online data — you can't recompute the full inverse every time a new observation arrives

Method 2 — Gradient Descent (Iterative)

Instead of solving for θ\theta θ analytically, gradient descent starts with a guess and repeatedly nudges θ\theta θ in the direction that most steeply reduces the cost.

The update rule at each step is:

$$\theta := \theta - \alpha \cdot \nabla_\theta J(\theta)$$

We already derived :-

$$\nabla_\theta J(\theta) = \frac{1}{m} X^T (X\theta - y)$$

Substituting :

$$\theta := \theta - \frac{\alpha}{m} X^T (X\theta - y)$$

That's it. One matrix multiply, one subtraction, repeated until convergence. Fully vectorized, scales to millions of observations.

When it's the right choice:

Large datasets where inverting XᵀX is infeasible
Online/streaming settings — update θ one batch at a time
When you're extending beyond linear regression (neural networks use gradient descent exclusively — there's no closed form)
When n ≫ m (more features than observations)

When it struggles:

Without feature scaling — the cost surface becomes a stretched ellipse and gradient descent oscillates wildly
With a poorly chosen learning rate α — too large and it diverges, too small and it crawls
It never reaches the exact minimum — it converges to a neighborhood of it

The Feature Scaling Problem

This is the most important practical detail in this entire section, and most tutorials bury it in a footnote.

Consider two features: GrLivArea (living area in sqft, range: 300–5000) and OverallQual (quality score, range: 1–10). Their scales differ by a factor of 500.

The cost function in (θ₁, θ₂) space becomes an extremely elongated ellipse — one axis spans a tiny range, the other spans a massive one. The gradient at any point mostly points along the elongated axis, not toward the minimum. Gradient descent bounces back and forth, making slow diagonal progress toward the center.

After Z-score standardization — subtracting mean and dividing by standard deviation — both features live on the same scale. The cost surface becomes nearly spherical. The gradient points almost directly at the minimum. Convergence that took 10,000 iterations now takes 120.

Critical implementation detail: Compute mean and standard deviation on the training set only, then apply those same values to the test set. Never fit the scaler on test data — that's data leakage.

# Correct way
mean = X_train.mean(axis=0)
std  = X_train.std(axis=0)

X_train_scaled = (X_train - mean) / std
X_test_scaled  = (X_test  - mean) / std   # use TRAIN statistics

Code: Gradient Descent from Scratch with Full Diagnostics

import numpy as np
import matplotlib.pyplot as plt

# Simple gradient descent class for linear regression
class GDSolver:
    def __init__(self, lr=0.01, max_steps=1000):
        self.lr = lr
        self.max_steps = max_steps
        self.weights = None
        self.costs = []
    
    def fit(self, X, y):
        m, n = X.shape
        self.weights = np.zeros(n)  # start at zero
        self.costs = []
        
        for step in range(self.max_steps):
            # predictions
            y_pred = X @ self.weights
            
            # mean squared error / 2
            error = y_pred - y
            cost = np.sum(error**2) / (2 * m)
            self.costs.append(cost)
            
            # gradient update
            gradient = X.T @ error / m
            self.weights = self.weights - self.lr * gradient
            
            # stop if cost stops changing
            if step > 0 and abs(self.costs[-1] - self.costs[-2]) < 1e-6:
                break
        
        return self
    
    def predict(self, X):
        return X @ self.weights

# Make some data
np.random.seed(42)
m = 200
x1 = np.random.uniform(300, 5000, m)  # house size
x2 = np.random.uniform(1, 10, m)      # quality
X = np.column_stack([np.ones(m), x1, x2])
true_w = np.array([50000, 30, 8000])
y = X @ true_w + np.random.normal(0, 5000, m)

# Scale features (important!)
X_mean = np.mean(X[:, 1:], axis=0)
X_std = np.std(X[:, 1:], axis=0)
X_scaled = np.column_stack([X[:, 0], (X[:, 1:] - X_mean) / X_std])

print("Testing different learning rates...")
rates = [0.000001, 0.01, 0.1, 0.5]

for lr in rates:
    gd = GDSolver(lr=lr, max_steps=500)
    gd.fit(X_scaled, y)
    print(f"LR={lr}: {len(gd.costs)} steps, final cost={gd.costs[-1]:.0f}")

# Plot convergence
plt.figure(figsize=(10, 4))

plt.subplot(1, 2, 1)
for i, lr in enumerate(rates):
    gd = GDSolver(lr=lr, max_steps=500)
    gd.fit(X_scaled, y)
    plt.plot(gd.costs, label=f'lr={lr}')

plt.xlabel('Iteration')
plt.ylabel('Cost')
plt.yscale('log')
plt.legend()
plt.title('How fast GD converges')

# Compare to exact solution
exact_w = np.linalg.inv(X_scaled.T @ X_scaled) @ X_scaled.T @ y
exact_cost = np.mean((X_scaled @ exact_w - y)**2) / 2
plt.axhline(exact_cost, color='k', linestyle='--', label='Exact OLS')
plt.legend()

plt.subplot(1, 2, 2)
# Show cost surface slice
w1_range = np.linspace(-2, 2, 50)
w2_range = np.linspace(-1, 1, 50)
W1, W2 = np.meshgrid(w1_range, w2_range)
costs = np.zeros_like(W1)

for i in range(50):
    for j in range(50):
        test_w = np.array([exact_w[0], W1[i,j], W2[i,j]])
        costs[i,j] = np.mean((X_scaled @ test_w - y)**2) / 2

plt.contourf(W1, W2, costs, 20, alpha=0.8)
plt.scatter(exact_w[1], exact_w[2], color='red', s=100, label='Exact solution')
plt.xlabel('weight 1')
plt.ylabel('weight 2')
plt.title('Cost surface')

plt.tight_layout()
plt.show()

print(f"\nExact weights: {exact_w}")
print(f"GD weights (lr=0.1): {GDSolver(lr=0.1).fit(X_scaled, y).weights}")

The Comparison Table

Property	Normal Equation	Gradient Descent
Solution type	Exact, closed-form	Iterative approximation
Computational cost	O(n³) — cubic in features	O(mn) per epoch
Feature scaling	Not required	Mandatory
Hyperparameters	None	Learning rate α, epochs
Multicollinearity	Fails (singular matrix)	Handles gracefully
Statistical inference	Standard errors fall out naturally	Requires extra computation
Large datasets	Impractical	Preferred
Online learning	Impossible	Natural fit

One Final Insight

When gradient descent fully converges on a convex problem like linear regression, it reaches the exact same point as the Normal Equation — the global minimum. They're not two different answers. They're two different paths to the same answer.

The Normal Equation takes a helicopter straight to the peak. Gradient descent hikes up the mountain step by step. Same summit, different journey.

In practice: Use the Normal Equation for small clean datasets where you need statistical inference. Use gradient descent for everything else. Every deep learning model in existence is just gradient descent applied to a much more complex cost surface.

When OLS Breaks: The 7 Assumptions and How to Diagnose Them

A model that always works is a model you don't understand yet.

Everything we've built so far — the MLE derivation, the Normal Equation, the projection — rests on a foundation of assumptions. OLS doesn't just work in general. It works when specific conditions hold. When they don't, your model may still produce numbers, it just produces the wrong ones, often confidently.

This section is about knowing when to trust your model and when to distrust it. That distinction is what separates a practitioner from someone who just runs code.

The Classical OLS Assumptions

The full statistical model we've been working with is:

$$y = X\theta + \varepsilon, \quad \varepsilon \sim N(0, \sigma^2 I)$$

Unpacking this single line reveals seven assumptions hiding inside it.

Assumption 1 — Linearity

What it says: The relationship between features and target is linear in the parameters.

What it doesn't say: That the raw features must be linear. You can include x², log(x), x₁ · x₂ as features — as long as the model is linear in θ, OLS is valid.

How it breaks: You fit a straight line through data that follows a curve. Your model systematically underpredicts in some regions and overpredicts in others.

Diagnostic: The Residuals vs Fitted plot. If the assumption holds, residuals scatter randomly around zero with no pattern. A U-shape or wave pattern is a clear signal of nonlinearity

import numpy as np
import matplotlib.pyplot as plt

def plot_residuals(y_true, y_pred):
    residuals = y_true - y_pred
    
    plt.figure(figsize=(8, 6))
    plt.scatter(y_pred, residuals, alpha=0.6)
    plt.axhline(y=0, color='r', linestyle='--')
    plt.xlabel('Predicted values')
    plt.ylabel('Residuals')
    plt.title('Residuals vs Fitted')
    plt.show()

# Test with some data
np.random.seed(42)
X = np.random.randn(100, 2)
X = np.column_stack([np.ones(100), X])
true_theta = [2, -1, 3]
y = X @ true_theta + np.random.normal(0, 1, 100)

# Fit model
theta_hat = np.linalg.inv(X.T @ X) @ X.T @ y
y_pred = X @ theta_hat

plot_residuals(y, y_pred)

Assumption 2 — Independence of Errors

What it says: εᵢ and εⱼ are uncorrelated for i ≠ j. Knowing one residual tells you nothing about another.

How it breaks: Time series data, spatial data, repeated measurements from the same subject. If you're predicting house prices street by street, houses on the same street have correlated errors — they share neighborhood effects your model didn't capture.

Consequence: OLS standard errors are underestimated. Your t-statistics look more significant than they are. You think your model is precise when it isn't.

Diagnostic: The Durbin-Watson statistic. It tests for first-order autocorrelation in residuals:

$$DW = \frac{\sum_{i=2}^m (e_i - e_{i-1})^2}{\sum_{i=1}^m e_i^2}$$

Values near 2 indicate no autocorrelation. Values near 0 indicate positive autocorrelation. Values near 4 indicate negative autocorrelation.

def durbin_watson(e):
    """Calculate Durbin-Watson statistic for residuals"""
    n = len(e)
    num = 0
    for i in range(n):
        num += e[i]**2
    
    den = 0
    for i in range(1, n):
        den += (e[i] - e[i-1])**2
    
    dw = num / den if den != 0 else 0
    print(f"DW = {dw:.3f}")
    
    if dw < 1.5:
        print("Positive autocorrelation!")
    elif dw > 2.5:
        print("Negative autocorrelation!")
    else:
        print("No autocorrelation ✓")
    
    return dw

# Test it
e_good = np.random.normal(0, 1, 50)  # random residuals
print("Good residuals:")
durbin_watson(e_good)

e_bad = np.arange(50) * 0.1  # correlated residuals
print("\nBad residuals:")
durbin_watson(e_bad)

Assumption 3 — Homoscedasticity

What it says: The variance of the errors is constant across all observations — Var(εᵢ) = σ² for all i.

How it breaks: Housing data is a classic example. Prediction errors for $100,000 homes are smaller in absolute terms than errors for $2,000,000 homes. Variance grows with the target value — this is heteroscedasticity.

Consequence: OLS estimates of θ remain unbiased, but confidence intervals and p-values are wrong. You might declare a coefficient insignificant when it isn't, or vice versa.

Fix: Log-transform the target variable. This is exactly why we used log(SalePrice) in the Ames dataset — it stabilizes the variance across the price range. Alternatively use Weighted Least Squares, where observations with high variance get lower weight.

Diagnostic: Residuals vs Fitted plot (fan shape = heteroscedasticity) and the Scale-Location plot (square root of standardized residuals vs fitted values):

import numpy as np
import matplotlib.pyplot as plt

def scale_location_plot(y_true, y_pred):
    residuals = y_true - y_pred
    
    # standardize residuals
    std_residuals = residuals / np.std(residuals)
    
    # take square root of absolute values
    sqrt_abs = np.sqrt(np.abs(std_residuals))
    
    plt.figure(figsize=(8, 6))
    plt.scatter(y_pred, sqrt_abs, alpha=0.6)
    plt.xlabel('Fitted values')
    plt.ylabel('sqrt(|standardized residuals|)')
    plt.title('Scale-Location Plot')
    plt.show()

# Test data
np.random.seed(42)
X = np.random.randn(100, 2)
X = np.column_stack([np.ones(100), X])
true_theta = [2, -1, 3]
y = X @ true_theta + np.random.normal(0, 1, 100)

theta_hat = np.linalg.inv(X.T @ X) @ X.T @ y
y_pred = X @ theta_hat

scale_location_plot(y, y_pred)

A flat horizontal trend means homoscedasticity holds. A rising trend means variance grows with fitted values — heteroscedasticity.

Assumption 4 — Normality of Errors

What it says: ε ~ N(0, σ²) — the errors follow a Gaussian distribution.

What it's actually needed for: Not for OLS to produce unbiased estimates — it doesn't need normality for that. It's needed for the t-tests and confidence intervals to be valid. By the Gauss-Markov theorem, OLS is the Best Linear Unbiased Estimator (BLUE) even without normality. But your p-values are only trustworthy when residuals are approximately Gaussian.

Diagnostic 1 — Q-Q Plot: Plots quantiles of residuals against theoretical Gaussian quantiles. Points on the diagonal = normality holds. Heavy tails curve away from the line.

Diagnostic 2 — Residual histogram: Should look like a bell curve centered at zero.

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def check_residual_normality(residuals):
    # Make two plots side by side
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    
    # Left: histogram + normal curve
    ax1.hist(residuals, bins=20, density=True, alpha=0.7, color='blue')
    
    # Normal distribution curve
    mu, sigma = np.mean(residuals), np.std(residuals)
    x = np.linspace(mu - 3*sigma, mu + 3*sigma, 100)
    normal_curve = stats.norm.pdf(x, mu, sigma)
    ax1.plot(x, normal_curve, 'r-', linewidth=2, label='Normal fit')
    ax1.set_title('Histogram of residuals')
    ax1.set_xlabel('Residuals')
    ax1.legend()
    
    # Right: Q-Q plot using scipy (easy!)
    stats.probplot(residuals, dist="norm", plot=ax2)
    ax2.set_title('Q-Q Plot')
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()
    
    # Shapiro test
    stat, p = stats.shapiro(residuals)
    print(f"Shapiro test: p-value = {p:.4f}")
    if p > 0.05:
        print("Residuals look normal ✓")
    else:
        print("Residuals NOT normal ✗")

# Test it
np.random.seed(42)
residuals = np.random.normal(0, 1, 100)  # normal residuals
check_residual_normality(residuals)

Assumption 5 — No Multicollinearity

What it says: The columns of X are linearly independent — no feature is a near-perfect linear combination of others.

How it breaks: Include both TotalSF and GrLivArea + TotalBsmtSF (which sum to nearly the same thing). Or include temperature in Celsius and Fahrenheit simultaneously.

Consequence: XᵀX becomes near-singular. The Normal Equation produces wildly unstable θ — individual coefficients blow up to +10⁶ and −10⁶ while their sum stays reasonable. Standard errors explode. t-statistics collapse. The model still predicts well but you can't trust any individual coefficient.

Diagnostic: Variance Inflation Factor (VIF) and correlation heatmap:

import numpy as np
import matplotlib.pyplot as plt

def vif_score(X, names=None):
    """Calculate VIF for each feature"""
    n = X.shape[1]
    vifs = []
    
    for j in range(n):
        # Feature j as target
        y_j = X[:, j]
        
        # Other features as predictors (+ bias)
        other_cols = [X[:, i] for i in range(n) if i != j]
        X_j = np.column_stack([np.ones(len(y_j)), other_cols])
        
        # OLS fit
        theta = np.linalg.inv(X_j.T @ X_j) @ X_j.T @ y_j
        y_pred = X_j @ theta
        
        # R squared
        ss_res = np.sum((y_j - y_pred)**2)
        ss_tot = np.sum((y_j - np.mean(y_j))**2)
        r2 = 1 - ss_res / ss_tot
        
        # VIF
        vif = 1 / (1 - r2) if r2 < 1 else 1000
        vifs.append(vif)
    
    if names is None:
        names = [f'x{i}' for i in range(n)]
    
    print("VIF Scores:")
    print("-" * 30)
    for name, v in zip(names, vifs):
        status = "OK" if v < 5 else ("Warning" if v < 10 else "BAD")
        print(f"{name}: {v:.1f} ({status})")
    
    return vifs

# Quick correlation plot
def corr_plot(X, names=None):
    corr = np.corrcoef(X.T)
    plt.figure(figsize=(8, 6))
    plt.imshow(corr, cmap='coolwarm', interpolation='none')
    plt.colorbar()
    if names:
        plt.xticks(range(len(names)), names, rotation=45)
        plt.yticks(range(len(names)), names)
    plt.title('Correlation matrix')
    plt.show()

# Test
np.random.seed(42)
X = np.random.randn(50, 3)
X[:, 2] = X[:, 0] + X[:, 1] * 0.5  # make collinear
vif_score(X, ['feat1', 'feat2', 'feat3'])
corr_plot(X, ['feat1', 'feat2', 'feat3'])

Fix: Drop one of the correlated features, use PCA to decorrelate them, or apply Ridge Regression — which is designed precisely for this situation.

Assumption 6 — No Influential Outliers

What it says: No single observation disproportionately controls the position of the regression line.

The distinction that matters: Not all outliers are influential. An outlier in y-space that sits right on the regression line has zero influence. An outlier in x-space (high leverage) that pulls the line toward it is dangerous.

Diagnostic: Cook's Distance — measures how much all fitted values change if observation i is removed:

$$D_i = \frac{\sum_{j=1}^m (\hat{y}^j - \hat{y}^j_{(i)})^2}{n \cdot \text{MSE}}$$

where ŷⱼ₍ᵢ₎ is the prediction when observation i is excluded. Observations with Dᵢ > 4/m are worth investigating.

import numpy as np
import matplotlib.pyplot as plt

def cooks_distance(X, y, theta):
    """Calculate Cook's distance using hat matrix"""
    m, n = X.shape  # m samples, n features
    
    # Predictions and residuals
    y_hat = X @ theta
    residuals = y - y_hat
    
    # MSE with df adjustment
    mse = np.sum(residuals**2) / (m - n)
    
    # Hat matrix diagonals (leverage)
    xtx_inv = np.linalg.inv(X.T @ X)
    leverage = np.array([X[i] @ xtx_inv @ X[i] for i in range(m)])
    
    # Cook's distance formula
    cooks_d = (residuals**2 / n) * leverage / ((1 - leverage)**2)
    
    # Plot
    plt.figure(figsize=(10, 5))
    colors = ['red' if d > 4/m else 'blue' for d in cooks_d]
    plt.bar(range(m), cooks_d, color=colors, alpha=0.7)
    plt.axhline(4/m, color='orange', linestyle='--', 
                label=f'Threshold 4/m = {4/m:.3f}')
    plt.xlabel('Observation')
    plt.ylabel("Cook's Distance")
    plt.title('Influential points (Cook D)')
    plt.legend()
    plt.show()
    
    # Count influential
    influential = np.sum(cooks_d > 4/m)
    print(f"Influential points: {influential}/{m}")
    
    return cooks_d

# Test data with outlier
np.random.seed(42)
X = np.random.randn(50, 2)
X = np.column_stack([np.ones(50), X])
y = X @ np.array([1, 2, 3]) + np.random.normal(0, 0.5, 50)

# Add outlier
y[0] = 10  # big outlier!

theta = np.linalg.inv(X.T @ X) @ X.T @ y
cooks_distance(X, y, theta)

Assumption 7 — Mean-Zero Errors

What it says: E[ε] = 0 — errors have no systematic direction.

How it breaks: You forgot the bias term. Your model is forced to pass through the origin when the true relationship doesn't. Every prediction is systematically off in the same direction.

Fix: Always include a column of ones in X. Always. This is non-negotiable.

Diagnostic: Simply check that mean(e) ≈ 0. If you included a bias term, this holds exactly by Property 2 from Section 4 — it's a mathematical guarantee, not an empirical check.

The Fix for Multicollinearity: Ridge Regression

When Assumption 5 breaks, the Normal Equation bec/umerically unstable. Ridge Regression adds a regularization term that fixes this directly:

$$J_{\text{ridge}}(\theta) = |y - X\theta|^2 + \lambda |\theta|^2$$

Taking the gradient and setting to zero:

$$\nabla J_{\text{ridge}} = -2 X^T (y - X\theta) + 2\lambda \theta = 0 $$

$$\theta_{\text{ridge}} = (X^T X + \lambda I)^{-1} X^T y$$

Adding λI guarantees the matrix is positive definite and invertible regardless of multicollinearity. The λ penalty shrinks coefficients toward zero — trading a small amount of bias for a large reduction in variance.

One critical implementation detail: The bias term should not be penalized. Penalizing the intercept introduces bias into your predictions that regularization is not supposed to create. Always set the first diagonal element of the penalty matrix to zero:

import numpy as np

def fit_ridge(X, y, lambda_reg=1.0):
    """
    This is ridge regression. Formula is theta = inv(X.T X + lambda I) X.T y
    But bias (first one) no penalty so I[0,0]=0
    """
    # get number of features (not including bias which is already in X col 0)
    num_features = X.shape[1]  
    # make identity matrix
    identity_matrix = np.eye(num_features)
    # important: dont penalize bias term (index 0)
    identity_matrix[0, 0] = 0  
    
    # compute X.T X 
    xtx = np.dot(X.T, X)  
    # add lambda * I (ridge part)
    ridge_matrix = xtx + lambda_reg * identity_matrix
    # now invert and multiply X.T y
    xt_y = np.dot(X.T, y)
    theta = np.linalg.inv(ridge_matrix) @ xt_y
    
    return theta

The Mindset Shift

Here's the thing most people miss: running diagnostics is not about finding a perfect model. It's about understanding in what ways your model is imperfect, and whether those imperfections matter for your use case.

A model with slight heteroscedasticity might still produce excellent point predictions — it just means your confidence intervals are slightly off. A model with one influential outlier might be fine if that outlier is a genuine data point you want to include.

Diagnostics don't tell you to throw your model away. They tell you exactly what to trust and what not to.

High R² with violated assumptions is not a good model. Modest R² with all assumptions verified is a model you can reason about, defend, and improve systematically.

Common Misconceptions About Linear Regression

Things that sound right, aren't.

Linear regression is one of the most taught algorithms in machine learning. It's also one of the most misunderstood. Not because it's complex — but because it's taught quickly, and quick explanations leave behind a residue of half-truths that harden into confident wrong beliefs over time.

Here are the six most common ones, and exactly why they're wrong.

Misconception 1 — "MSE is just a convenient mathematical choice"

The wrong belief: We minimize squared errors because it's mathematically clean — squares are differentiable, absolute values aren't. It's a computational convenience, nothing more.

Why it's wrong: We spent all of Section 2 proving this isn't true. MSE is the maximum likelihood estimator under Gaussian noise. It isn't convenient — it's correct, given a specific and principled assumption about how noise behaves in your data.

The choice of loss function is not aesthetic. It encodes a probabilistic belief about your data-generating process. MSE says: I believe my errors are symmetric, continuous, and Gaussian. If that belief is wrong — if your errors are heavy-tailed, skewed, or discrete — then MSE is the wrong loss function, and you should know why.

The rule: Loss function choice = probabilistic assumption about noise. Always know what assumption your loss encodes.

Loss Function	Implicit Noise Assumption
MSE (L2)	Gaussian noise
MAE (L1)	Laplacian noise
Huber Loss	Gaussian + occasional outliers
Log Loss	Bernoulli (binary outcomes)

Misconception 2 — "High R² means a good model"

The wrong belief: R² = 0.92 means your model explains 92% of the variance. That's excellent. The model is good.

Why it's wrong: R² measures how much variance your model explains on the training data, relative to a horizontal mean line. It tells you nothing about whether your assumptions hold, whether your model generalizes, or whether your coefficients are meaningful.

You can have R² = 0.99 with severe heteroscedasticity (your confidence intervals are garbage), multicollinearity so bad that coefficients are meaningless, a nonlinear relationship your linear model is accidentally fitting well locally, or overfitting so extreme the model memorizes noise.

And you can have R² = 0.65 with a perfectly valid, well-specified model on genuinely noisy data where 65% explainable variance is the theoretical ceiling.

Anscombe's Quartet makes this point devastatingly well — four completely different datasets with wildly different structures, all with R² ≈ 0.67 and nearly identical regression lines. Always plot your data. Always run diagnostics.

The rule: R² is a necessary condition for a good model, not a sufficient one. A high R² with violated assumptions is not a good model. It's a confident wrong model.

Misconception 3 — "Linear regression assumes the data is linear"

The wrong belief: If your data has curves, linear regression can't handle it. You need a nonlinear model.

Why it's wrong: Linear regression is linear in the parameters, not in the features. You can feed it any transformation of your features.

As long as the model is y = θ₀ + θ₁φ₁(x) + θ₂φ₂(x) + … where φᵢ are any functions of x and θ enters linearly — it's still linear regression. The Normal Equation still applies. The MLE derivation still holds. Everything we've built still works.

The linearity assumption is about the relationship between y and θ, not between y and raw x.

Misconception 4 — "Gradient descent and the Normal Equation give different answers"

The wrong belief: They're two different optimization methods so they find different solutions, especially on complex data.

Why it's wrong: The cost function J(θ) = (1/m)‖y − Xθ‖² is strictly convex — it has exactly one global minimum and no local minima. Both methods are finding that same unique minimum.

The Normal Equation finds it analytically in one step. Gradient descent approaches it iteratively. Given enough epochs and a small enough learning rate, gradient descent converges to the exact same θ as the Normal Equation — to within floating point precision.

Misconception 5 — "Standardizing features changes the model"

The wrong belief: If I standardize my features, I'm changing the data and therefore getting a different model with different predictions.

Why it's wrong: Standardization is an invertible linear transformation. It changes the scale of the parameter space, not the underlying relationship. The predictions ŷ are identical whether you standardize or not — as long as you apply the same transformation consistently to train and test data.

What standardization does change is the geometry of the cost surface — from an elongated ellipse to a sphere. This affects how fast gradient descent converges, not where it converges to.

What standardization does not change: predictions, R², residuals, or the underlying model.

The one thing to be careful about: interpret coefficients on the standardized scale. A coefficient of 0.35 on a standardized feature means "a one standard deviation increase in this feature is associated with a 0.35 unit change in y." This is actually more interpretable for comparing feature importance across features with different units.

Misconception 6 — "More features always means a better model"

The wrong belief: Adding features can only help. More information is always better. If a feature has any correlation with y, include it.

Why it's wrong: Three reasons.

First, multicollinearity. Adding a feature that's correlated with existing features doesn't add information — it destabilizes the coefficient estimates for all correlated features. You get worse inference without better predictions.

Second, overfitting. Every feature you add gives your model more degrees of freedom to fit noise in the training data. Adjusted R² penalizes for this. The bias-variance tradeoff is real — a model with 50 features on 100 observations is almost certainly memorizing training noise.

Third, the curse of dimensionality. In high dimensions, data becomes sparse. The intuitions about distance and density that work in 2D break down in 100D. More features can actively hurt generalization.

Results on the Ames Housing Dataset

Theory is only as good as what it produces on real data.

The Pipeline, End to End

No sklearn. Every step implemented from scratch:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# 1. Load the data
df = pd.read_csv('data/train.csv')
print("Data loaded, shape:", df.shape)

# 2. Find best features using correlation with SalePrice
# only numerical columns
num_cols = df.select_dtypes(include=[np.number]).columns.tolist()
# correlation with target, take abs and sort
corr_with_price = df[num_cols].corr()['SalePrice'].abs().sort_values(ascending=False)
# top 10 features (skip SalePrice itself)
top_10_features = corr_with_price[1:11].index.tolist()
print("Top features:", top_10_features)

# 3. Fill missing values with median for each feature
for feature in top_10_features:
    df[feature] = df[feature].fillna(df[feature].median())  # fixed chained assignment

# 4. Log transform y because prices are skewed
y_original = df['SalePrice'].values
y = np.log(y_original)  # log scale stabilizes variance

# 5. Make X matrix, add column of 1s for bias
X_data = df[top_10_features].values
X = np.column_stack([np.ones(len(X_data)), X_data])

# 6. Train test split manually 70/30
np.random.seed(42)  # same randomness
all_indices = np.random.permutation(len(y))
split_point = int(0.7 * len(y))
train_inds = all_indices[:split_point]
test_inds = all_indices[split_point:]

X_train = X[train_inds]
X_test = X[test_inds]
y_train = y[train_inds]
y_test = y[test_inds]

print("Train size:", len(y_train), "Test size:", len(y_test))

# 7. Standardize features using train mean/std only (no leak!)
feature_means = X_train[:, 1:].mean(axis=0)
feature_stds = X_train[:, 1:].std(axis=0)

# apply to train
X_train[:, 1:] = (X_train[:, 1:] - feature_means) / feature_stds
# apply to test (use train stats)
X_test[:, 1:] = (X_test[:, 1:] - feature_means) / feature_stds

# 8. Fit linear regression using normal equation
# theta = inv(X.T X) X.T y
xtx_inv = np.linalg.inv(X_train.T @ X_train)
x_ty = X_train.T @ y_train
theta = xtx_inv @ x_ty

# 9. Metrics functions
def calc_r2(y_true, y_pred):
    resid_sq = np.sum((y_true - y_pred)**2)
    total_sq = np.sum((y_true - np.mean(y_true))**2)
    return 1 - resid_sq / total_sq

def calc_rmse(y_true, y_pred):
    return np.sqrt(np.mean((y_true - y_pred)**2))

def calc_adj_r2(y_true, y_pred, num_feats):
    n_samples = len(y_true)
    r2_score = calc_r2(y_true, y_pred)
    return 1 - (1 - r2_score) * (n_samples - 1) / (n_samples - num_feats - 1)

# Predictions
train_preds = X_train @ theta
test_preds = X_test @ theta

# Print results same format
print("\n" + "=" * 50)
print("  Final Results — Ames Housing Dataset")
print("=" * 50)
print(f"\n  {'Metric':<20} {'Train':>10} {'Test':>10}")
print("  " + "-" * 42)
print(f"  {'R²':<20} {calc_r2(y_train, train_preds):>10.4f} {calc_r2(y_test, test_preds):>10.4f}")
print(f"  {'Adjusted R²':<20} {calc_adj_r2(y_train, train_preds, 10):>10.4f} {calc_adj_r2(y_test, test_preds, 10):>10.4f}")
print(f"  {'RMSE (log scale)':<20} {calc_rmse(y_train, train_preds):>10.4f} {calc_rmse(y_test, test_preds):>10.4f}")
print("=" * 50)

Results

Metric	Train	Test
R²	0.87	0.82
Adjusted R²	0.86	0.81
RMSE (log scale)	0.14	0.17
Convergence	120 epochs	—

82% of variance in housing prices explained. 10 features. No sklearn. Built entirely from the math we derived in the previous sections.

Predicted vs Actual

import matplotlib.pyplot as plt
import numpy as np

# make two subplots side by side
fig, ax1, ax2 = plt.subplots(1, 2, figsize=(13, 5))

# set dark background like github
fig.patch.set_facecolor('#0d1117')
ax1.set_facecolor('#161b22')
ax2.set_facecolor('#161b22')

# Left plot: actual vs predicted
ax1.scatter(y_test, y_test_pred, color='#58a6ff', alpha=0.4, s=20)
# perfect line y=x
min_val = min(y_test.min(), y_test_pred.min())
max_val = max(y_test.max(), y_test_pred.max())
ax1.plot([min_val, max_val], [min_val, max_val], color='#ff7b72', linewidth=2, linestyle='--', label='Perfect fit')
ax1.set_xlabel('Actual log(SalePrice)', color='#8b949e')
ax1.set_ylabel('Predicted log(SalePrice)', color='#8b949e')
ax1.set_title('Predicted vs Actual', color='white', fontsize=12)
ax1.tick_params(colors='#8b949e')
ax1.legend(facecolor='#161b22', labelcolor='white')

# Right plot: residuals histogram
resid = y_test - y_test_pred  # residuals
ax2.hist(resid, bins=35, density=True, color='#58a6ff', alpha=0.6, edgecolor='#1f2937')

# overlay normal distribution curve
x_vals = np.linspace(resid.min(), resid.max(), 300)
resid_mean = resid.mean()
resid_std = resid.std()
# gaussian formula
gauss_curve = (1 / (resid_std * np.sqrt(2 * np.pi))) * np.exp(-0.5 * ((x_vals - resid_mean) / resid_std)**2)
ax2.plot(x_vals, gauss_curve, color='#ff7b72', linewidth=2.5, label=r'\(\mathcal{N}(\mu, \sigma^2)\)')

ax2.set_xlabel('Residual', color='#8b949e')
ax2.set_ylabel('Density', color='#8b949e')
ax2.set_title('Residual Distribution (Test Set)', color='white', fontsize=12)
ax2.tick_params(colors='#8b949e')
ax2.legend(facecolor='#161b22', labelcolor='white')

plt.tight_layout()
plt.savefig('ames_results.png', dpi=150, bbox_inches='tight', facecolor='#0d1117')
plt.show()

Points hug the diagonal. Residuals are approximately Gaussian. The assumptions we verified in Section 6 hold — which is why the metrics are trustworthy.

The Honest Takeaway

The score is not the point. What matters is that every number here is earned — derived from first principles, implemented from scratch, and validated against the assumptions that make those numbers meaningful.

A 0.82 R² you understand is worth more than a 0.95 R² you don't.

Why This Matters

The real payoff isn't the formula. It's what understanding it unlocks.

You Become a Better Debugger

When a model fails — and they always fail eventually — there are two kinds of practitioners. The first opens a GitHub issue or tries a different library. The second knows exactly where to look.

Understanding linear regression from first principles gives you a causal map of failure modes. Predictions are systematically off in one region? Linearity assumption — check your residuals vs fitted plot. Coefficients are wildly large with opposite signs? Multicollinearity — check your VIF scores. Model performs perfectly on training data and collapses on test? Overfitting — check your adjusted R² and feature count. Confidence intervals are nonsense despite good predictions? Heteroscedasticity — log-transform your target or use weighted least squares.

You don't guess. You diagnose. That's the difference.

Every Advanced Model Builds on This

Linear regression is not a stepping stone you leave behind. It's the foundation everything else is built on.

Logistic regression is linear regression with a sigmoid applied to the output and cross-entropy replacing MSE as the loss — same gradient descent, different likelihood function. Ridge and Lasso are linear regression with a prior over parameters — Gaussian prior gives Ridge, Laplacian prior gives Lasso. Principal Component Analysis finds the directions in the column space of X that capture maximum variance — pure linear algebra, same subspace geometry from Section 4. Neural networks are stacked linear transformations with nonlinear activations between them — the output layer of every regression network is literally ŷ = Wᵀx + b, our equation. Gaussian Processes generalize the Bayesian interpretation of linear regression to infinite-dimensional feature spaces.

The thread connecting all of these is the same: find parameters that make your data most probable under some assumed noise model. That's MLE. That's what we derived in Section 2. It never stops being relevant.

The Mental Model That Transfers

Here is the most transferable idea in this entire post, stated plainly: every supervised learning algorithm is an answer to the question "what parameters make my observed data most probable?" The algorithm differs. The question never does.

Change the noise assumption from Gaussian to Bernoulli and you get logistic regression. Change it to Poisson and you get Poisson regression. Add a prior over parameters and you get Bayesian regression. Replace the dot product with a kernel function and you get kernel regression. Stack nonlinear transformations and you get a neural network.

Linear regression teaches you the skeleton. Everything else adds flesh to it.

What This Project Actually Proves

You built something here that most practitioners haven't. Not because the code is complex — it isn't. Because you started from P(y | X, θ) and ended at predicted_sale_prices.csv without touching a single high-level abstraction you didn't understand.

That means when something breaks in a future project, you have a place to stand. You know what the math requires, what the assumptions demand, and what the geometry looks like. That knowledge doesn't expire when a new library comes out or a new architecture becomes fashionable.

Understanding is the only thing in this field that compounds.

Command Palette

The Hook

The intuition : What is Linear Regression Actually Trying to Do ?

Where MSE comes From : The Probabilistic Argument

The Setup

Step 1 — Write the Likelihood

Step 2 — Take the Log

Step 3 — The Punchline

What This Tells You About When MSE Fails

Code: Visualizing the Gaussian Noise Assumption

The Algebraic Argument: Deriving the Normal Equation

Setting Up the Cost Function

Taking the Gradient

Setting the Gradient to Zero

Code: OLS from Scratch

Visualizing the Cost Surface

The Geometric Argument : Regression as Orthogonal Projection

The Projection Principle

The Question : Which point ŷ ∈ col(X) is closest to y?

The Diagram

Three Properties That Follow Immediately

Property 2 — Residuals sum to zero (when a bias term is included):

Code: Verifying Orthogonality from Scratch

Two Ways to Solve It: Normal Equation vs Gradient Descent

Method 1 — The Normal Equation (Closed-Form)

Method 2 — Gradient Descent (Iterative)

One Final Insight

When OLS Breaks: The 7 Assumptions and How to Diagnose Them

The Classical OLS Assumptions

The Fix for Multicollinearity: Ridge Regression

The Mindset Shift

Common Misconceptions About Linear Regression

Misconception 1 — "MSE is just a convenient mathematical choice"

Misconception 2 — "High R² means a good model"

Misconception 3 — "Linear regression assumes the data is linear"

Misconception 4 — "Gradient descent and the Normal Equation give different answers"

Misconception 5 — "Standardizing features changes the model"

Misconception 6 — "More features always means a better model"

Results on the Ames Housing Dataset

The Pipeline, End to End

Results

The Honest Takeaway

Why This Matters

Comments