Linear Regression & the LMS Algorithm

Chapter 1: The Hypothesis — A Line Is Just Two Knobs

In the last chapter you dragged a line around with two handles. Those two handles are the entire model. Let's name them.

A straight line that predicts price from area has the form

Read it aloud: "the prediction h at input x equals theta-zero plus theta-one times x." The symbol h stands for hypothesis — historical jargon for "the function we're proposing as the rule." The two numbers θ₀ and θ₁ are the parameters (also called weights). They are the knobs. Everything about machine learning is choosing good values for knobs like these.

Each knob has a concrete meaning you can see on the plot:

• θ₀ — the intercept. The predicted price when area is zero. It slides the whole line up and down. Think of it as the "base price" the model starts from before it has looked at the size.
• θ₁ — the slope. How many extra thousands of dollars each additional square foot is worth. It tilts the line. A slope of 0.14 means "$140 more per square foot."

Watch what each knob does on its own:

More than one feature

Living area isn't the only thing that sets a price. Suppose we also know the number of bedrooms. Now each house is described by two numbers, and our line becomes a tilted plane:

Here x₁ is living area, x₂ is bedrooms, and we've gained a third knob θ₂ — "dollars per bedroom." There's nothing special about stopping at two. With d features we'd write h(x) = θ₀ + θ₁x₁ + … + θ_dx_d. Writing that sum out every time is tedious, so let's find a cleaner notation.

The x₀ = 1 trick

That lonely θ₀ with no x attached is annoying — it breaks the pattern. Here's a slick fix that the whole field uses: invent a fake feature x₀ that is always equal to 1. Then

and suddenly every term has the same shape θ_jx_j. The hypothesis collapses into one tidy sum:

That last form, θ^Tx, is the dot product of the parameter vector θ = (θ₀, θ₁, …, θ_d) and the feature vector x = (x₀, x₁, …, x_d) = (1, x₁, …, x_d). The little ^T means "transpose" — it just lines the two vectors up so we can multiply them entry-by-entry and add. The dot product is "multiply each feature by its weight and sum the results." That's the beating heart of a linear model, and you'll meet it again in logistic regression, neural networks, and attention layers. Learn it once here.

A worked prediction, by hand

Let's compute one prediction with actual numbers so the dot product stops being abstract. Suppose after training we found

(These are close to the real CS229 fitted values for area + bedrooms.) A house has 2,500 ft² and 4 bedrooms. With the x₀=1 convention, its feature vector is x = (1, 2500, 4). The prediction is the dot product:

# h(x) = θ₀·x₀ + θ₁·x₁ + θ₂·x₂
h(x) = 89.6 × 1  +  0.139 × 2500  +  (−8.74) × 4
     = 89.6      +  347.5         +  (−34.96)
     = 402.14

So the model predicts about $402k. Notice every step: we multiplied each feature by its weight (89.6×1, 0.139×2500, −8.74×4), then added the three products. The negative weight on bedrooms is real and counterintuitive — given a fixed area, more bedrooms means smaller rooms, which the data says is slightly less desirable. We'll revisit that surprise. For now, just internalize: prediction = weighted sum of features.

python
import numpy as np

theta = np.array([89.6, 0.139, -8.74])   # θ₀ (bias), θ₁, θ₂
x     = np.array([1.0, 2500, 4])         # x₀=1 (bias trick), area, beds

# Three equivalent ways to compute the same prediction:
pred = theta[0]*x[0] + theta[1]*x[1] + theta[2]*x[2]  # by hand
pred = np.sum(theta * x)                            # elementwise multiply, then sum
pred = theta @ x                                    # the dot product θᵀx
print(pred)   # 402.14

Chapter 2: Measuring Wrongness — The Cost Function

We have a model with knobs. We have a vague feeling that some knob settings are better than others. To turn that feeling into something a machine can optimize, we need to score any choice of θ with a single number: small when the line fits well, large when it fits badly. Then "find the best line" becomes "find the θ that makes the score smallest" — a math problem, not an opinion.

Start with one house. The line predicts h(x⁽ⁱ⁾); the house actually sold for y⁽ⁱ⁾. The gap between them,

is called the residual — how far off we were on house i. (The superscript (i) is just an index labelling the i-th training example; it is not an exponent. h(x⁽²⁾) means "the prediction for house number 2.") A positive residual means we overpredicted; negative means we underpredicted.

We want a total score across all n houses. Three natural ideas, and why two of them fail:

Squaring does two jobs at once. First, it kills the sign — (−50)² = 2500, same as (+50)² — so errors can't cancel. Second, it is a smooth, bowl-shaped function with a derivative everywhere, which (next chapter) lets us roll downhill on it. There's also a third, deeper reason rooted in probability that we'll earn in Chapter 8. For now, accept squaring as the natural choice and define the cost function (also called the least-squares or mean-squared-error objective):

Two questions you should be asking. Why the ½ out front? Purely for convenience: when we differentiate the square next chapter, the 2 from the power rule cancels the ½, leaving a clean expression. It doesn't change where the minimum is — scaling a function by a constant doesn't move its lowest point — so we're free to put it there. And why is J a function of θ, not x? Because the data (the x's and y's) is fixed — it's the spreadsheet, it doesn't move. The only things we can change are the knobs θ. J(θ) asks: "for this setting of the knobs, how much total squared error do we rack up across all the houses?"

Computing J by hand on a tiny dataset

Abstractions lie; numbers don't. Let's compute J for a 3-point dataset, every step shown. Take these three points and the candidate line h(x) = 1 + 1·x (so θ₀=1, θ₁=1):

Now sum the last column and halve it:

# J(θ) = ½ Σ (h(xᵢ) − yᵢ)²
J = ½ × (1 + 0 + 0)
  = ½ × 1
  = 0.5

So this line scores J = 0.5. Is that good? On its own the number is meaningless — cost is only useful for comparison. Let's try a worse line, h(x) = 0 + 1·x (drop the intercept to 0):

# Line h(x) = x:  predictions are 1, 2, 3
residuals = (1−1), (2−3), (3−4)   =   0, −1, −1
squared   =    0,    1,    1
J = ½ × (0 + 1 + 1) = ½ × 2 = 1.0

This line scores J = 1.0 — twice as costly as the first. So the machine now has a definite, defensible verdict: the line h(x) = 1 + x fits these three points better than h(x) = x. No eyeballing, no opinion — just a smaller number. That comparison is the whole game. The rest of the lesson is finding the θ with the very smallest J.

python
import numpy as np

X = np.array([[1,1],[1,2],[1,3]])  # columns: x₀=1, x₁
y = np.array([1,3,4])

def cost(theta):
    pred = X @ theta              # h(xᵢ) for every row at once
    resid = pred - y
    return 0.5 * np.sum(resid**2)   # ½ Σ residual²

print(cost(np.array([1,1])))  # 0.5  — the good line
print(cost(np.array([0,1])))  # 1.0  — the worse line

Chapter 3: Rolling Downhill — Gradient Descent

We have a score J(θ) and a goal: find the θ that makes it smallest. We could try every possible (θ₀, θ₁) pair — but there are infinitely many, and with d features the search space explodes. We need a smarter method. Here's one of the most important algorithms in all of computing, and it's based on a child's intuition: if you want to reach the bottom of a valley in thick fog, just keep stepping downhill.

Picture the cost J as a landscape. The ground position is the knob setting (θ₀, θ₁); the height at that position is the cost J(θ). For least-squares linear regression this landscape is a single smooth bowl — a convex paraboloid — with exactly one lowest point and no other dips to get trapped in. (We'll prove the "one minimum" claim in Ch. 6.) You can't see the whole bowl, but standing at any point you can feel which way is downhill. So: step downhill, feel again, step again. Repeat until the ground is flat. That flat spot is the minimum.

The update rule

"Feel the slope in the θ_j direction" is precisely the partial derivative ∂J/∂θ_j — the rate at which cost changes as you nudge knob j. If that derivative is positive, increasing θ_j increases cost, so we should decrease θ_j. If negative, the opposite. Either way, we move against the derivative. That gives the gradient descent update:

The := means "assign" — overwrite the old θ_j with the new value (like theta_j = theta_j - ... in code), not an equation claiming the two sides are already equal. The number α (alpha) is the learning rate: how big a step to take. Small α = timid baby steps; large α = bold leaps. We'll see both fail in their own way.

Working out the derivative (every step)

To run the algorithm we need an actual formula for ∂J/∂θ_j. Let's derive it for the simplest case — a single training example (x, y), so the sum in J disappears and J = ½(h(x) − y)². Watch the chain rule peel it apart:

Look at how clean that is. This is why we put the ½ in J: the 2 from squaring and the ½ out front annihilate each other, leaving no stray constant. And the last step is the magic — differentiating the sum θ₀x₀ + θ₁x₁ + … with respect to θ_j, every term is a constant (w.r.t. θ_j) except the one containing θ_j, whose derivative is just its coefficient x_j. So:

Plug that into the update and you get gradient descent for linear regression. We'll unpack what it means — the famous LMS rule — in the very next chapter. First, let's see the descent happen and feel how α controls everything.

The Goldilocks problem of α

The learning rate is the single most temperamental knob in machine learning. Three regimes, all visible in the widget above:

• α too small (e.g. 0.01): every step is a sliver. It converges, but you might wait thousands of iterations. Correct, but glacial.
• α just right: each step makes strong progress and the path curves smoothly into the bullseye. This is what you tune for.
• α too large (push it past ~1.0): a step overshoots the minimum and lands higher up the opposite wall than it started. The next step overshoots worse. Cost grows each iteration and the parameters fly off to infinity — the dreaded divergence.

There's no universal best α; it depends on the data's scale. A standard remedy (Ch. 9) is to rescale features so the bowl is round rather than a long thin canyon — round bowls forgive larger steps.

python
def gradient_descent(X, y, alpha=0.1, iters=1000):
    theta = np.zeros(X.shape[1])          # start at the origin
    n = len(y)
    for _ in range(iters):
        resid = X @ theta - y                # h(xᵢ) − yᵢ for all i
        grad  = X.T @ resid                  # Σ (h−y)·xⱼ  →  ∂J/∂θ for every j at once
        theta = theta - alpha * grad         # the update, all knobs together
    return theta

Chapter 4: The LMS Update Rule

Substitute the derivative we just derived, ∂J/∂θ_j = (h(x) − y)x_j, back into the gradient-descent update θ_j := θ_j − α·∂J/∂θ_j. The two minus signs (one in the update, one if we flip the residual) let us write it in the form everyone quotes. For a single training example (x⁽ⁱ⁾, y⁽ⁱ⁾):

This is the LMS update rule — "LMS" stands for Least Mean Squares — also known as the Widrow–Hoff learning rule, after the engineers who used it in adaptive filters in 1960. It is one of the most important three-term expressions in machine learning, and once you understand it you understand the skeleton of how almost every model on Earth learns, including deep neural networks. So let's read it slowly.

Notice we flipped the residual: instead of (h − y) we wrote the error (y − h), "target minus prediction," and the sign in front became a plus. The piece (y⁽ⁱ⁾ − h_θ(x⁽ⁱ⁾)) is how far below our prediction the truth was — positive when we underpredicted, negative when we overpredicted. Three behaviors fall out of this one term:

• The update is proportional to the error. If our prediction is nearly perfect, (y − h) ≈ 0, so we barely change θ — we don't fix what isn't broken. If we're wildly off, the error is large and we make a big correction. The algorithm spends its effort exactly where it's wrong.
• The sign of the error sets the direction. Underpredicted (y > h, error positive)? Then θ_j increases — push the prediction up next time. Overpredicted? θ_j decreases. The rule self-corrects toward the truth.
• The feature x_j⁽ⁱ⁾ gates and scales the blame. Knob θ_j only moves in proportion to how "active" its feature was on this example. If x_j = 0, that feature played no role in the prediction, so it receives no correction. If x_j is large, that feature drove the prediction hard, so it absorbs more of the fix. Credit and blame flow to the features that were actually responsible.

One LMS step, fully by hand

Let's execute exactly one update with real numbers so nothing is hidden. Setup: one training example with area x₁ = 2 (in thousands of ft², so the numbers stay readable) and true price y = 5. We use the bias trick, so x = (x₀, x₁) = (1, 2). Current knobs θ = (θ₀, θ₁) = (0, 1). Learning rate α = 0.1.

Step 1 — predict. h(x) = θ^Tx = θ₀·1 + θ₁·2 = 0·1 + 1·2 = 2.

Step 2 — error. y − h(x) = 5 − 2 = 3. We underpredicted by 3, so both knobs should rise.

Step 3 — update each knob using θ_j := θ_j + α(y − h)x_j:

# θ₀:  feature x₀ = 1
θ₀ := 0 + 0.1 × 3 × 1  =  0 + 0.3  =  0.3
# θ₁:  feature x₁ = 2  (so it gets twice the correction of θ₀)
θ₁ := 1 + 0.1 × 3 × 2  =  1 + 0.6  =  1.6

New knobs: θ = (0.3, 1.6). Step 4 — verify we improved. Re-predict: h(x) = 0.3·1 + 1.6·2 = 0.3 + 3.2 = 3.5. The new error is 5 − 3.5 = 1.5 — we cut it in half, from 3 down to 1.5, in a single step. And notice θ₁ moved twice as far as θ₀ (0.6 vs 0.3), exactly because its feature x₁=2 was twice as large as x₀=1. The feature gated the blame, just as promised.

python
# One LMS update for a single example (x, y). x includes the bias x₀=1.
def lms_step(theta, x, y, alpha):
    h     = theta @ x                 # prediction θᵀx
    error = y - h                     # target − prediction
    return theta + alpha * error * x     # θⱼ += α·(y−h)·xⱼ  for all j

theta = np.array([0.0, 1.0])
x, y  = np.array([1.0, 2.0]), 5.0
print(lms_step(theta, x, y, 0.1))   # [0.3 1.6] — matches our hand calc

Chapter 5: Batch vs. Stochastic Gradient Descent

The LMS rule of Chapter 4 updated θ from a single example. But we have n houses, not one. How do we use the whole training set? There are two answers, and the difference between them — "look at everything before each step" vs. "step after every example" — is one of the most consequential design choices in large-scale machine learning. Let's build both.

Batch gradient descent: average the wisdom of the crowd

The true cost J(θ) = ½∑_i(h(x⁽ⁱ⁾) − y⁽ⁱ⁾)² sums over all examples, so its true derivative also sums over all examples. We worked out the single-example derivative; the full one is just the sum of those:

Plug it in and you get batch gradient descent: every step scans the entire training set, accumulates the error contribution from all n houses, and takes one well-informed step:

"Batch" because one step consumes the whole batch of data. Each step points in the exact direction of steepest descent on the true cost, so the path is smooth and the convergence is reliable. The cost: every single step requires touching all n examples. With n in the millions, one step is expensive, and you need many steps.

Stochastic gradient descent: step on every example

The alternative throws out the "wait until you've seen everything" rule. Stochastic gradient descent (SGD, also called incremental gradient descent) loops over the examples and applies the single-example LMS update immediately after each one:

Each step uses one example's error as a noisy estimate of the true gradient. Any single example might point slightly wrong — but on average they point downhill, so the parameters drift toward the minimum while jittering around the path. The payoff: SGD starts improving after the very first example, before batch GD has even finished its first scan. On a huge dataset this is transformative.

The price: SGD never quite settles. Near the minimum, each example still yanks θ a little in its own direction, so the parameters keep oscillating in a small cloud around the true optimum rather than landing dead center. In practice that cloud is close enough — and you can shrink it by slowly decreasing α toward zero as training proceeds, which lets the jitter die down and the parameters converge. For large datasets, SGD's head start usually wins, which is why essentially all deep learning uses SGD (or a mini-batch hybrid of the two).

Why the batch gradient really is the sum of LMS terms

Let's verify the "sum of single-example gradients" claim with a 2-point dataset, because it's the conceptual hinge. Points: (x⁽¹⁾, y⁽¹⁾) = ((1,1), 3) and (x⁽²⁾, y⁽²⁾) = ((1,2), 5), with current θ = (0, 1).

# Predictions
h⁽¹⁾ = θᵀx⁽¹⁾ = 0·1 + 1·1 = 1     →  error y−h = 3−1 = 2
h⁽²⁾ = θᵀx⁽²⁾ = 0·1 + 1·2 = 2     →  error y−h = 5−2 = 3

# Per-example contribution to ∂J/∂θⱼ is (h−y)·xⱼ = −(error)·xⱼ
# θ₀ direction (x₀=1 for both):
  example 1:  −2 × 1 = −2
  example 2:  −3 × 1 = −3        sum = −5
# θ₁ direction (x₁ = 1 then 2):
  example 1:  −2 × 1 = −2
  example 2:  −3 × 2 = −6        sum = −8

# Batch gradient ∇J = (−5, −8). One batch step with α=0.1:
θ := (0,1) − 0.1·(−5,−8) = (0.5, 1.8)

The batch step moved both knobs using the combined evidence of both houses at once. SGD, by contrast, would have stepped on house 1, changed θ, then stepped on house 2 from the new θ — two smaller, sequential nudges instead of one summed leap. Same destination, different road.

Chapter 6: The Closed Form — Normal Equations

Gradient descent creeps toward the minimum over many iterations. But for linear regression there's something almost magical: we can jump straight to the answer in one shot, solving an equation instead of iterating. The cost is a perfect bowl, and calculus says the bottom of a bowl is exactly where the slope is zero in every direction. So instead of rolling downhill, set the gradient to zero and solve.

To do that for all parameters at once without drowning in subscripts, we switch to matrix notation. Stack the training inputs as rows of the design matrix X (n rows, one per house; d+1 columns, including the x₀=1 bias column), and stack the targets in a vector y:

Now the magic of matrix multiplication: the matrix–vector product Xθ computes all n predictions at once. Row i of X is house i's features, and (row i)·θ = h_θ(x⁽ⁱ⁾). So Xθ is the vector of every prediction, and Xθ − y is the vector of every residual. Stacking residuals into a vector r, the cost is half the sum of their squares — and "sum of squares of a vector's entries" is exactly r^Tr (a vector dotted with itself). Therefore:

This is the same J as before — same ½, same sum of squared residuals — just written compactly so we can differentiate it in one move instead of d+1 separate ones.

Setting the gradient to zero (every step)

Expand the product and differentiate with respect to the vector θ. Three matrix-calculus facts do all the work: (a) a^Tb = b^Ta (dot products commute), (b) ∇_θ(b^Tθ) = b, and (c) ∇_θ(θ^TAθ) = 2Aθ for symmetric A. Here A = X^TX, which is always symmetric. Watch:

Set that gradient to zero — the defining condition for the bottom of the bowl — and you get the famous normal equations:

This is a plain linear system: a known matrix (X^TX) times the unknown θ equals a known vector (X^Ty). Multiply both sides by the inverse of X^TX and you have the parameters in closed form — no iteration, no learning rate, no convergence to wait for:

Solving the normal equations by hand

Let's grind one out completely on the 3-point set from Chapter 2: points (x,y) = (1,1), (2,3), (3,4). With the bias column, the design matrix and target are

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

Step 1 — form X^TX (a 2×2 matrix). X^T is 2×3; multiply by X (3×2):

# (XᵀX)₀₀ = Σ x₀·x₀ = 1+1+1 = 3
# (XᵀX)₀₁ = Σ x₀·x₁ = 1+2+3 = 6
# (XᵀX)₁₁ = Σ x₁·x₁ = 1+4+9 = 14
XᵀX = [ 3   6 ]
      [ 6  14 ]

Step 2 — form X^Ty (a length-2 vector):

# (Xᵀy)₀ = Σ x₀·y = 1·1 + 1·3 + 1·4 = 8
# (Xᵀy)₁ = Σ x₁·y = 1·1 + 2·3 + 3·4 = 1+6+12 = 19
Xᵀy = [ 8 ]
      [19 ]

Step 3 — invert X^TX. For a 2×2 matrix [[a,b],[c,d]] the inverse is 1/(ad−bc)·[[d,−b],[−c,a]]. The determinant is ad−bc = 3·14 − 6·6 = 42 − 36 = 6.

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

Step 4 — multiply θ = (X^TX)⁻¹X^Ty:

θ = (1/6) [ 14  −6 ] [ 8 ]
          [ −6   3 ] [19 ]

θ₀ = (1/6)(14·8  + (−6)·19) = (1/6)(112 − 114) = (1/6)(−2) = −0.333
θ₁ = (1/6)((−6)·8 + 3·19)   = (1/6)(−48 + 57)  = (1/6)(9)  =  1.5

So the exact best-fit line is h(x) = −0.333 + 1.5x. No iterations, no learning rate — just arithmetic. Sanity check at x=2: h = −0.333 + 3.0 = 2.667, sitting neatly between the y-values 1, 3, 4. This is the provably optimal line; gradient descent from Chapter 3, run long enough, would converge to these same two numbers.

python
X = np.array([[1,1],[1,2],[1,3]], dtype=float)
y = np.array([1,3,4], dtype=float)

# Closed form θ = (XᵀX)⁻¹ Xᵀy — but prefer the numerically stable solver:
theta = np.linalg.solve(X.T @ X, X.T @ y)
print(theta)                         # [-0.333  1.5]

# Or skip the normal equations entirely — lstsq is more robust:
theta, *_ = np.linalg.lstsq(X, y, rcond=None)
print(theta)                         # [-0.333  1.5]

When the inverse doesn't exist

The formula needs (X^TX)⁻¹, which only exists when X^TX is invertible (non-singular). Two situations break it: (1) more features than examples (d+1 > n) — not enough data to pin down all the knobs; and (2) redundant features — e.g. you include area in both square feet and square meters; one is a constant multiple of the other, so they carry no independent information and X^TX collapses. Both are fixable (drop the redundant feature, gather more data, or add a regularization term that makes the matrix invertible again — ridge regression), but the bare formula above silently fails. This is one practical reason gradient descent stays useful: it just keeps creeping downhill and never has to invert anything.

Chapter 7: The Living Regression Lab

Everything so far — the hypothesis, the cost, gradient descent, the LMS rule, batch vs. SGD, the normal equations — now collides in one playground. This is the chapter where it all becomes one moving picture. Take your time here; if the lesson clicks anywhere, it clicks here.

The left panel is the data and the fitted line. The right panel is the cost landscape — the same bowl from Chapter 3, seen from above as contour rings — with a moving dot showing where the current (θ₀, θ₁) sits on it. As the line on the left improves, the dot on the right slides downhill toward the bullseye. Watch the two panels together: a better line is literally a lower point on the bowl. That correspondence — "fitting the data" = "descending the cost surface" — is the single most important idea in the whole lesson, and here you can see both halves at once.

Sit with one observation before moving on: gradient descent and the normal equations land on the same spot. Run GD to convergence, note the line; hit "Normal equations," and the line doesn't budge. Two utterly different procedures — one a thousand tiny downhill steps, the other a single matrix inversion — agree exactly, because both are solving the identical problem: find the bottom of this one convex bowl. When two independent methods agree, you've usually found something real.

(No quiz this chapter — the lab is the test. If you can predict, before clicking, what dragging a point will do to the line and to the dot on the bowl, you understand linear regression.)

Chapter 8: Why Squares? The Probabilistic Story

Back in Chapter 2 we chose to square the residuals, and the reasons were practical: squaring kills the sign and gives a smooth derivative. Honest reasons — but they leave a nagging question. Why squares and not, say, fourth powers, or absolute values? Those also kill the sign. Is least squares just a convenient convention, or is it somehow the right thing to do? This chapter delivers the deep answer, and it's one of the most satisfying results in all of statistics: if the world is a clean line corrupted by Gaussian noise, then minimizing squared error is exactly equivalent to finding the most probable parameters. Least squares isn't arbitrary. It's maximum likelihood in disguise.

A model of where the data comes from

Let's stop treating the y-values as arbitrary and propose a generative story: each true price is given by the perfect line plus some random scatter,

where ε⁽ⁱ⁾ (epsilon) is an error term — everything the line doesn't capture: features we left out (the house's view, the neighbourhood school), measurement noise, plain randomness. We now make a specific, famous assumption about that noise: the ε⁽ⁱ⁾ are independent and identically distributed (IID) and follow a Gaussian (Normal) distribution with mean 0 and some variance σ², written ε⁽ⁱ⁾ ~ 𝒩(0, σ²). Mean 0 says the noise doesn't systematically push prices up or down; the bell shape says small deviations are common and big ones are rare. Its density is the bell curve:

Since y⁽ⁱ⁾ = θ^Tx⁽ⁱ⁾ + ε⁽ⁱ⁾, and θ^Tx⁽ⁱ⁾ is a fixed number once we pick θ, the randomness in y comes entirely from ε. So y⁽ⁱ⁾ is also bell-shaped, centered on the line's prediction. Substituting ε⁽ⁱ⁾ = y⁽ⁱ⁾ − θ^Tx⁽ⁱ⁾ into the density:

Read this as: "given a house's features x and a choice of line θ, here's the probability the house sold for price y." A price right on the line is most probable; prices far from the line are exponentially less probable — and notice that exponent already contains a squared residual. The squares are creeping in, not by our choice this time, but from the shape of the bell curve itself.

Maximum likelihood: pick the θ that makes the data most probable

Here's the principle. Different lines θ assign different probabilities to the prices we actually observed. The maximum likelihood principle says: choose the θ under which the observed data was most probable. It's common sense — prefer the explanation that makes what you saw least surprising.

Because the examples are independent, the probability of the whole dataset is the product of the individual probabilities. Viewed as a function of θ, this product is called the likelihood L(θ):

Maximizing a product of exponentials is awkward. The standard trick: take the logarithm. The log is monotonically increasing, so whatever θ maximizes L also maximizes log L — but the log turns the product into a sum and cancels the exp. Call it the log-likelihood ℓ(θ):

The punchline

Now look at that last line as a function of θ. The first term, n·log(1/√(2π)σ), has no θ in it — it's a constant we can ignore. The 1/σ² is a positive constant multiplier. So maximizing ℓ(θ) means maximizing −½∑(y − θ^Tx)², which — because of the minus sign — means minimizing:

That is exactly our least-squares cost J from Chapter 2. We didn't aim for it; it fell out. Maximizing the probability of the data under Gaussian-noise assumptions is identical to minimizing squared error. The squares were never an arbitrary choice — they are the fingerprint of the Gaussian's exp(−something²) shape. Choose a different noise distribution and you'd get a different cost: Laplacian (double-exponential) noise gives absolute error (the robust loss of Ch. 9). Least squares is the right answer to a specific, named assumption about the world.

One more elegant detail: notice σ² multiplied the whole sum but dropped out of the answer — the best-fit θ doesn't depend on how noisy we think the data is. That's why we never had to know or estimate σ to do linear regression. The noise level changes how confident we should be in the fit, but not the fit itself.

python
# Maximizing Gaussian log-likelihood == minimizing ½Σ(y − ŷ)². Demo:
def neg_log_likelihood(theta, X, y, sigma):
    resid = y - X @ theta
    n = len(y)
    # −ℓ(θ) up to constants that don't depend on θ:
    return n*np.log(np.sqrt(2*np.pi)*sigma) + (1/sigma**2)*0.5*np.sum(resid**2)

# The θ that minimizes −ℓ is the SAME θ that minimizes ½Σresid², for ANY σ>0.
# → least squares IS maximum likelihood under Gaussian noise.

Chapter 9: When It Breaks — Underfitting, Overfitting, Outliers

You now own a complete, working algorithm. A good engineer's next question is always: how does it fail? Knowing the failure modes is what separates "I can run LinearRegression().fit()" from "I understand linear regression." There are four ways it breaks, and each one points to a whole branch of the field.

1. Underfitting — the model is too simple

A straight line can only ever be straight. If the true relationship curves — prices that rise fast then plateau for huge mansions — no choice of θ₀, θ₁ will fit it well. The line systematically misses, leaving a pattern in the residuals. This is underfitting: the model lacks the flexibility to capture the structure that's plainly in the data. The fix is to give it more expressive features.

And here's the trick that surprises people: you can fit a curve with linear regression. Just invent new features that are powers of x — feed the model x, x², x³ as separate columns. The hypothesis h(x) = θ₀ + θ₁x + θ₂x² + θ₃x³ is a wiggly cubic in the input x, yet it's still perfectly linear in the parameters θ, so every formula we derived still applies unchanged. "Linear" regression was never about straight lines — it's about being linear in the knobs.

2. Overfitting — the model is too flexible

If a little curvature helps, why not a lot? Fit a 5th-order polynomial and the curve can pass through every training point exactly, driving J all the way to zero. Victory? No — disaster. Between the points the curve thrashes wildly up and down, and on a new house it predicts nonsense. The model memorized the training data, including its noise, instead of learning the trend. This is overfitting: low training cost, terrible real-world predictions. Remember the Chapter 2 warning — low cost is fit, not wisdom. Here it bites.

The cure is regularization — add a penalty on large θ values to the cost, so the model is rewarded for staying simple unless the data really demands complexity. (That's ridge and lasso regression, and conveniently the ridge penalty also fixes the non-invertibility problem below.) The bias–variance tradeoff that governs this balance gets its own lesson.

3. Outliers — the Achilles' heel of squaring

The very thing that made squared error convenient now backfires. Because errors are squared, a single wild point — a data-entry typo, a mansion sold to a relative for $1 — contributes a huge penalty (its residual squared), and the line bends hard to reduce it, abandoning the other 46 honest points. Least squares is not robust: one bad point can wreck the whole fit. This is the Gaussian assumption from Chapter 8 failing — real noise had a heavy tail the bell curve didn't expect. The remedy is a loss that grows more slowly than the square (absolute error, Huber loss), which is exactly the robust estimation story.

4. Singular X^TX and unscaled features

Two more practical traps. First, as we saw in Chapter 6, if features are redundant (one is a linear combination of others) or you have fewer examples than features, X^TX has no inverse and the normal equations simply don't define an answer. Second, for gradient descent, features on wildly different scales (area in the thousands, bedrooms in single digits) stretch the cost bowl into a long thin canyon, so any α that's safe for one direction crawls hopelessly in the other. The fix is feature scaling — standardize each feature to mean 0, variance 1 — which rounds the bowl and lets gradient descent take confident steps. Always scale your features before running GD.

Chapter 10: Connections & Cheat Sheet

You started this lesson unable to defend any particular line through a scatter of houses. You can now: define a hypothesis, score it with a principled cost, minimize that cost two different ways, and — the deep part — explain why squared error is the right cost when noise is Gaussian. That's a real, complete piece of machine learning, mastered from zero. More importantly, you've met the template that nearly every other model reuses.

The whole lesson on one page

Symbol glossary

Production code — the whole thing, both ways

python
import numpy as np

def fit_normal_equations(X, y):
    """Exact closed-form solution. Best for few features."""
    return np.linalg.lstsq(X, y, rcond=None)[0]   # stabler than (XᵀX)⁻¹Xᵀy

def fit_gradient_descent(X, y, alpha=0.01, iters=5000):
    """Iterative. Scales to huge data. Remember to feature-scale X first."""
    theta = np.zeros(X.shape[1])
    for _ in range(iters):
        theta -= alpha * (X.T @ (X @ theta - y))   # θ −= α·∇J
    return theta

# sklearn one-liner — does the normal equations under the hood:
from sklearn.linear_model import LinearRegression
model = LinearRegression().fit(X_raw, y)     # handles the bias for you
print(model.coef_, model.intercept_)

Where to go next

• Logistic regression — the very next CS229 chapter. Swap the continuous target for a 0/1 label, pass θ^Tx through a sigmoid, and you have classification. The update rule looks identical to LMS — that's not a coincidence, and GLMs explain why.
• Generalized Linear Models (GLMs) — the grand unification: linear regression (Gaussian) and logistic regression (Bernoulli) are two members of one exponential-family framework. The shared shape of their update rules falls right out.
• Locally weighted regression — refit a line near each query point, weighting nearby data more. Fits curves without choosing polynomial degrees.
• Robust estimation — what to do when outliers wreck least squares (Huber, RANSAC, M-estimators).
• Gaussian processes — take the probabilistic view of Chapter 8 to its limit: a full distribution over functions, with error bars on every prediction.
• Regression (EE269) — the same ideas from a signal-processing angle.
• Loss functions — MSE is one entry in a whole zoo; see how cross-entropy and contrastive losses generalize the "score how wrong" idea.

"With four parameters I can fit an elephant, and with five I can make him wiggle his trunk." — John von Neumann, on the danger of overfitting. You now know exactly what he meant.