Linear Regression from First Principles

OLS, normal equation, geometric interpretation, assumptions

Here is a true story from over two hundred years ago. On the very first night of 1801, an astronomer spotted a faint new dot in the sky — a small world we now call Ceres. He watched it for about forty nights, and then it drifted behind the glare of the Sun and vanished. Everyone wanted it back, but to know where to point a telescope, you had to predict its whole orbit from just a handful of shaky, imperfect measurements. Nobody could do it.

Then a 24-year-old named Carl Friedrich Gauss tried something new. He took all those messy measurements and found the single orbit that fit them best — not perfectly, because the measurements had errors in them, but best in a way he could actually prove. He told the astronomers where to look. Ceres was right there.

The idea Gauss leaned on is called least squares, and it is the beating heart of what we now call linear regression. Before it, drawing the best line through noisy data was guesswork — everyone eyeballed it and got a different answer. Least squares turned "fit a line to the dots" from an art into a solved problem with one right answer. Two hundred years later it is still one of the first tools any data scientist reaches for. Let me show you how it works, brick by brick.

The setup

Forget planets for a moment and take something homelier. You sell houses. For every house that sold last month you wrote down a few facts — how big it is, how many bedrooms, how nice the area is — and the price it sold for. Now a new house comes up, and you need to guess its price.

Your instinct is the right one. A bigger house costs more. More bedrooms cost more. A nicer area costs more. So take each fact, give it a weight for how much it matters, and add them up:

price ≈ w₁ × size + w₂ × bedrooms + w₃ × area score

That is the whole model. Pick the three weights and you can price any house. So the entire game comes down to one question: how do we pick good weights?

To make every step ahead checkable by hand, shrink this down to one fact — size — and five real houses. Sizes, in hundreds of square feet: House A = 10, B = 15, C = 20, D = 25, E = 30 (that's 1,000 to 3,000 sqft). Prices, in thousands of dollars: 200, 250, 280, 310, 360. These ten numbers are the running example for the rest of this module — every weight, every residual, every score below is computed from exactly this table, so you can check each one yourself.

How wrong are we?

To find good weights, we first have to measure bad ones. You already know what last month's houses actually sold for, so use your weights to predict those same houses and compare each guess to the truth. How far off you were on one house is called the residual — a fancier word for the miss.

One miss is easy. But we want a single score for how wrong the weights are across *all* the houses. The most natural idea is to average the misses. Careful though — some guesses are too high and some too low, so the misses have opposite signs and cancel out. A terrible model could average to nearly zero just by luck.

The fix is to throw away the signs. Take the size of each miss, ignore whether it was over or under, and average those. That gives you the mean absolute error, or MAE — plainly, "on average, how many dollars off are we?" It is honest, easy to read, and a couple of wild houses barely move it.

But when it comes to actually *finding* the best weights, we usually reach for MAE's cousin. Instead of taking the size of each miss, we square it, then average. That is the mean squared error, or MSE, and it does two things MAE does not. First, squaring makes one big miss count for far more than several small ones, so the model works hardest to avoid the embarrassing, way-off guesses. Second — and this is the part that matters most — squaring makes the whole thing smooth, which, as you are about to see, is exactly what lets us solve for the best weights in one clean shot.

Either way, this single number — how wrong the weights are, all together — is called the loss. And now the whole job fits in one sentence: find the weights that make the loss as small as possible.

The beautiful part: solving it in one step

Here is where least squares earns its fame. Picture the loss as a landscape: the weights are your position, and the loss is the height of the ground under you. Because we squared the misses, that landscape is a smooth bowl — one single lowest point, no other dips to fall into. And the bottom of a bowl is exactly the spot where the ground goes flat.

So we do not have to wander around trying weights and checking the loss. We just ask math for the point where the slope is flat, and it hands back one formula that gives the best weights directly:

$θ̂ = (XᵀX)⁻¹Xᵀy$

Do not worry about the symbols yet. The thing to feel is the magic: for a straight-line model, we never search — we solve, in a single step, and get the provably best weights every time. That is the same trick that let Gauss find a lost planet. The one-step solution is called ordinary least squares, or OLS, and that is the entire engine.

For one feature, that matrix formula collapses into arithmetic you can run by hand — and this is where the five houses stop being a story and become a computation. First find the mean size, x̄ = (10+15+20+25+30)/5 = 20, and the mean price, ȳ = (200+250+280+310+360)/5 = 280. Now measure how size and price move together: multiply each house's distance from the mean size by its distance from the mean price, and add those up. House A: (10−20)×(200−280) = (−10)×(−80) = 800. B: (−5)×(−30) = 150. C: (0)×(0) = 0. D: (5)×(30) = 150. E: (10)×(80) = 800. Total: Sxy = 1900 — call this the *co-movement* sum. Next measure how spread-out size alone is: square each house's distance from the mean size and add those up: 100+25+0+25+100 = Sxx = 250 — the *size-spread* sum. The best-fit slope is just their ratio: slope = Sxy / Sxx = 1900 / 250 = 7.6 — every extra hundred square feet adds \$7,600, or \$76 per square foot. The intercept anchors the line at the means: intercept = ȳ − slope×x̄ = 280 − 7.6×20 = 280 − 152 = 128 — a house with (hypothetically) zero size would price at \$128k. So OLS hands back one exact model: price = 128 + 7.6 × size, no searching required.

Check it against the houses it was fit on. House A (size 10): 128 + 7.6×10 = 204, but it actually sold for 200 — a residual of 200 − 204 = −4. B (15): 128+114=242 vs 250, residual +8. C (20): 128+152=280 vs 280, residual 0 — a perfect hit, because the line passes through the means exactly. D (25): 128+190=318 vs 310, residual −8. E (30): 128+228=356 vs 360, residual +4. Notice the residuals sum to exactly zero (−4+8+0−8+4=0) — no accident, it's a direct consequence of the flat-slope condition that produced this line in the first place.

A warning about trusting the weights

Now a catch that trips up almost everyone. Suppose two of your facts move together — bigger houses usually have more bedrooms, so size and bedroom count rise and fall as a pair. The model wants to hand out credit for the price, some to size and some to bedrooms. But because the two always move together, it cannot tell which one is really doing the work. It might load most of the weight onto size, or most onto bedrooms, and both choices predict prices about equally well.

So the weights turn shaky. Train on a slightly different batch of houses and those two weights can jump around — one run gives size a positive weight, the next run gives it a negative weight, as if a bigger house should cost *less*. Meanwhile the actual price predictions barely change. This is called collinearity. And notice carefully what is and is not moving: the two facts are still exactly as correlated as before — that is a fixed fact about your data and it does not budge. It is the *weights* that have gone unreliable, not the correlation.

Why care? Because people read the weights to decide which fact matters. See a weight near zero and think "this one does nothing, drop it," and you can be badly wrong — the weight may be tiny only because its twin grabbed the credit. The usual fix is a method called Ridge, which keeps the weights from growing wild and shaky. The Regularisation module picks up exactly here — and reuses this module's own Sxy=1900 and Sxx=250 to show, in real numbers, how Ridge and Lasso each tame a weight like our slope of 7.6.

How good is the fit, really?

Say you have trained the model and its total squared error comes out to some big number. Is that good? On its own it means nothing. Fifty million squared-dollars — compared to what?

You need a yardstick, and the fairest one is the dumbest possible model: the one that ignores every fact about the house and just guesses the same number every time — the average price. Call it the lazy model. How wrong is the lazy model? Take each house's price, subtract the average, square it, and add it all up. On the five houses: (200−280)²+(250−280)²+(280−280)²+(310−280)²+(360−280)² = 6400+900+0+900+6400 = 14,600. That total, 14,600, is the sum of squared deviations — add up how far every price sits from the average, squared. Divide it by how many houses you have (5) and you get the *average* squared spread, which has a famous name: the variance: 14,600 ÷ 5 = 2,920. (Take its square root and you are back in real dollars — that is the standard deviation, ≈ $54.0k. And the average you started from is the mean. Mean, then variance, then standard deviation, each one built on the one before.)

Now compare it to your model's own squared error — the one you can compute right now, from the residuals worked out above (−4, +8, 0, −8, +4): 16+64+0+64+16 = 160. Your model is stuck with 160; the lazy model was stuck with 14,600.

Now comes the clever comparison. R² just asks: of all the error the lazy mean-model was stuck with, how much did your model clear away?

R² = 1 − (your model's squared error ÷ the lazy model's squared error) = 1 − (160 ÷ 14,600) = 1 − 0.0110 = 0.989

Read it straight off. R² = 0 means you are no better than guessing the average — useless. R² = 1 means you nailed every house. Our 0.989 means the size-only model wiped out 98.9% of the error the lazy model had — size alone very nearly explains these five prices. Now the number finally means something, because it is measured against a floor.

R² is the first thing to look at when you judge a linear regression — but it hides one trap. Add a new fact to the model, even one that is pure random noise, and R² never goes down. The model can always use junk to shave a sliver off the training error. So plain R² quietly rewards piling on useless features. The fix is adjusted R², which charges a small fee for every fact you add: adjusted R² = 1 − (1−R²)×(n−1)/(n−p−1), where n is the number of houses and p the number of features. With n=5, p=1: 1 − (1−0.989)×(4/3) = 1 − 0.0110×1.333 = 1 − 0.0146 = 0.985 — a little lower than plain R², exactly the fee for spending one feature. Add a second feature that is pure noise and p becomes 2: even if R² ticks up to, say, 0.991, adjusted R² now pays a steeper fee (n−p−1 shrinks from 3 to 2) and can fall below 0.985 even as plain R² rises — that's the whole point of the adjustment. So when you are deciding whether a feature belongs, trust adjusted R², not plain R².

One last habit — the most useful one

Even a high R² can fool you, and this final trick catches it. Our own R² of 0.989 is a great score — but suppose the real relationship actually curves, and you fit a straight line anyway. You can still score an impressive R² and still be wrong about the shape.

To catch it, take your misses and plot them against your predictions. If the straight line is right, the misses should scatter randomly around zero — no pattern at all. Plot our five: prediction 204 → residual −4; 242 → +8; 280 → 0; 318 → −8; 356 → +4. No trend, no U, no fan — small, sign-flipping misses with no relationship to the prediction size, exactly what a genuinely linear fit should look like. But if the truth was a curve, the misses fall into a clear shape — often a U, too low in the middle and too high at the ends. That U is the model quietly telling you the straight line is the wrong shape. R² will never warn you. The plot of the misses always will.

So before you trust any straight-line model: plot the misses and look. Gauss would have.

The fine print: what OLS quietly assumes

That residual plot isn't just a nicety — it's how you check the promises OLS silently makes. There are five, and it helps to split them into two piles.

The first pile is what you need just to *predict* well. One: the relationship really is a straight line (linearity). Two: no feature is a near-copy of another (no perfect multicollinearity) — the collinearity trap from earlier. If these hold, your predictions are trustworthy even if nothing else does.

The second pile is what you additionally need to *trust the weights and their uncertainty* — the inference layer we're about to meet. Three: the errors have zero mean *given the features* (exogeneity) — the model isn't systematically wrong in any region. Four: the errors all have the same spread (homoscedasticity). Five: the errors are independent of each other (no autocorrelation). Break these and your *predictions* can still be fine, but every confidence interval and p-value you quote is off. That split — prediction assumptions versus inference assumptions — is the thing most people fumble in an interview.

When the spread isn't even: heteroscedasticity

Look again at a residual plot that fans out — tight on the left, wide on the right. That's heteroscedasticity: the errors are bigger for pricier houses. Here's the subtle part interviewers love. Your weights are *still unbiased* — on average they're centered on the truth. What breaks is the *uncertainty* around them: the standard errors, confidence intervals, and hypothesis tests all become unreliable, usually too optimistic. So you keep predicting fine but start *believing* things about your weights that aren't warranted. The fixes: robust ("sandwich") standard errors, or a transform like modeling log(price) instead of price.

One weird house can move the whole line: leverage and influence

"What if a single data point changes your slope?" is a classic — and you can watch it happen to our own five-house line. Add a sixth house, F: 8,000 sqft (size 80 in our hundreds-of-sqft units), sold for \$300k. Two different things are going on. Leverage is a point with an extreme *feature* value — F sits at size 80, far out past the other five, which cluster between 10 and 30. That alone gives it the *potential* to swing the line hard, just by being far out — under the original line (price = 128 + 7.6×size), F "should" cost 128+7.6×80 = 736k, but it sold for only 300k, a huge gap. Influence is when a point *actually* does swing the fit — high leverage *and* a price that fights the trend, which is exactly F's situation.

Refit with all six houses and watch it happen: the new means are x̄=30, ȳ≈283.3, and the same co-movement/spread calculation from before now gives Sxy=2,900 and Sxx=3,250, so slope = 2900/3250 ≈ 0.89 — down from 7.6. One added house out of six collapsed the slope by 88%, from "\$76 per square foot" to "\$9 per square foot." Cook's distance measures exactly this: how much every weight would move if you deleted that one house — and it would flag House F immediately. A point with big Cook's distance is one row quietly steering your whole model. You find these by looking, not by trusting the summary metrics — R² won't flinch.

Reading the weights like a statistician: the inference layer

So far we've used the weights to *predict*. But often the real question is "does size actually matter, or did we imagine it?" That needs the inference layer, and our five houses (before House F showed up) have exactly enough numbers to compute it. Around each weight you compute a standard error — how much that weight would jitter across different samples. It's built from the residual variance: MSE_resid = SSE/(n−2) = 160/3 ≈ 53.33 (that "−2" spends one degree of freedom on the slope, one on the intercept), then SE(slope) = √(MSE_resid / Sxx) = √(53.33/250) ≈ 0.462. Divide the weight by its standard error and you get a t-statistic: 7.6/0.462 ≈ 16.4 — with only 3 degrees of freedom left (n−2=3), a t this large is far out in the tail, giving a p-value well under 0.001 — the odds of seeing a slope this big by chance, if size truly had no effect on price, are vanishingly small. Wrap the weight in the 95% critical value for 3 degrees of freedom (t≈3.18) and you get a confidence interval: 7.6 ± 3.18×0.462 ≈ 7.6 ± 1.47, or [6.13, 9.07]. In dollars per square foot: size adds about \$76/sqft, and we're 95% confident it's between \$61/sqft and \$91/sqft — that's the difference between "\$76" and a number you can actually stand behind.

And why is OLS the natural tool for this? The Gauss-Markov theorem: when those assumptions hold, OLS is BLUE — the Best Linear Unbiased Estimator, meaning among all unbiased straight-line methods, none has smaller variance. That's the deep reason least squares earns its place, not just tradition. (Note: scikit-learn gives you the weights but not p-values or intervals — for those you reach for statsmodels.)

How the computer actually solves it

We wrote the one-step answer as $θ̂ = (XᵀX)⁻¹Xᵀy$, but no careful library computes it that literally. Forming XᵀX and inverting it *squares* how sensitive the math is to nearly-collinear features, so it can blow up numerically. Real solvers (scikit-learn included) instead use QR or SVD decompositions — same answer in exact arithmetic, far more stable when features are correlated or poorly scaled. Worth knowing that the textbook formula and the production code disagree on purpose.

Picking the right yardstick

R² tells you how much you beat the lazy mean-model, but for reporting error you'll usually quote one of a few — and our five houses give real numbers for each. MAE is the plain average dollar miss: (4+8+0+8+4)/5 = \$4,800 — easy to read, shrugs off a few wild houses. RMSE squares before averaging, so it punishes big misses harder and stays in dollars: √(SSE/n) = √(160/5) = √32 ≈ \$5,657 — noticeably above the \$4,800 MAE, exactly because RMSE lets the two \$8k misses (D and B) count for more than the two \$4k misses. Use RMSE when large errors are especially costly. MAPE reports the miss as a *percentage* of each price: (4/200 + 8/250 + 0/280 + 8/310 + 4/360)/5 × 100 ≈ 1.78% — this travels across scales but explodes when true values are near zero and punishes over-prediction unevenly. And R² (0.989, computed above) is the unitless "fraction of variance explained" for a quick sense of fit. Match the metric to the question: absolute dollars (MAE/RMSE), relative error (MAPE), or overall fit (R²).

Key points

• In one line: OLS finds the weights that make the total squared miss as small as possible — slope = Sxy/Sxx, intercept = ȳ − slope×x̄. Use linear regression first whenever you are predicting a number and "add up the facts, each with its own weight" is a sensible guess. Its big advantage is that you can read the weights straight off — on our five-house example, slope 7.6 means "+\$7,600 per 100 sqft," computed directly from Sxy=1900 and Sxx=250, no search required. It is fast, it is simple, and it is the baseline every fancier model has to beat. Move to something else when the residual plot curves (the real shape is not a straight line), when facts clearly work together, or when the answer has to stay inside fixed limits — like a probability between 0 and 1, which a straight line cannot respect.
• The trap that catches people: when two facts move together, their weights stop being trustworthy — but the predictions still look fine, so nothing warns you. If size and total rooms almost always rise together, OLS splits the credit between them however it likes, and a different batch of houses splits it differently. The two weights wobble, but their combined effect stays steady — so the predictions look healthy. That is why "this weight is near zero, let us drop the fact" is dangerous: the weight may be small only because its twin took the credit. Add Ridge (a small penalty) before you trust any list of which facts matter.
• The one check to run every time: plot the misses against the predictions. If the misses scatter evenly around zero, the straight line fits. If they form a U or any clear pattern, the shape is wrong — and a curved plot with R² = 0.95 is worse than a messy one with R² = 0.60, because now you are confidently wrong. If the misses fan out wider as predictions grow, the errors get bigger for pricier houses; the predictions can still be fine, but any confidence range you quote is off. R² hides all of this. The plot of the misses shows it.
• Split the assumptions into two piles — predicting well needs less than trusting the weights. To *predict* you mainly need a genuinely linear relationship and no near-duplicate features. To also *trust the standard errors, p-values, and intervals*, you additionally need errors with zero mean given the features (exogeneity), equal spread (homoscedasticity), and independence (no autocorrelation). Heteroscedasticity is the classic gotcha: it leaves your weights unbiased but makes their uncertainty wrong, so predictions stay fine while every confidence interval quietly lies. Fix with robust standard errors or a log transform.
• Know the difference between leverage and influence — and that one row can steer the whole model. Leverage is a point with an extreme *feature* value (it has the potential to swing the line); influence is when it actually does (extreme feature value *and* a target that fights the trend). Add one 8,000 sqft house priced at \$300k to our five-house fit and the slope collapses from 7.6 to about 0.89 — an 88% drop from a single added row; Cook's distance is the number that would flag exactly that row. And this is what the inference layer buys you: on the original five houses, SE(slope)≈0.462, t≈16.4, giving a 95% confidence interval of [6.13, 9.07] — "+\$76/sqft, 95% sure it's \$61–\$91/sqft" instead of a bare "\$76." OLS is the natural tool because Gauss-Markov says it is BLUE — the minimum-variance unbiased linear estimator — when the assumptions hold. (scikit-learn gives weights only; use statsmodels for p-values and intervals.)
• Match the error metric to the question, and don't compute the normal equation literally. On our five houses: MAE = \$4,800 (average dollar miss, robust, readable); RMSE = √32 ≈ \$5,657 (squares first, so it punishes the two \$8k misses harder than the two \$4k ones — quote it when large errors are costly); MAPE ≈ 1.78% (percentage error, travels across scales but explodes near zero); R² = 0.989 (unitless fraction of variance explained). On implementation: nobody careful inverts XᵀX directly — squaring the condition number is numerically dangerous — so real solvers (scikit-learn included) use QR or SVD for the same answer with far more stability on correlated or badly scaled features.

Recap

• OLS = weights that minimise total squared miss. Squaring: signs don't cancel, big misses dominate, loss is a smooth bowl.
• One-step solve: θ̂ = (XᵀX)⁻¹Xᵀy → for one feature, slope = Sxy/Sxx = 1900/250 = 7.6, intercept = ȳ−slope·x̄ = 128 (5-house example).
• R² = fraction of the lazy mean-model's error you cleared → 1−160/14,600 = 0.989; adjusted R² (0.985) charges a fee per feature (junk features never lower plain R²).
• Always plot the misses vs predictions. Our five residuals (−4,+8,0,−8,+4) scatter with no pattern — that's what a genuine linear fit looks like; a U-pattern says wrong shape, and R² can't see it.
• Leverage vs influence: one added extreme point (8,000 sqft, \$300k) collapsed slope 7.6→0.89, an 88% swing — Cook's distance flags exactly this.
• Inference layer: SE(slope)≈0.462, t≈16.4, 95% CI [6.13, 9.07] — from MSE_resid=SSE/(n−2) and SE=√(MSE_resid/Sxx).
• Two assumption piles: predict well (linearity, no near-duplicate features) vs trust the weights (exogeneity, homoscedasticity, independence). Heteroscedasticity → weights unbiased but standard errors lie.
• Real solvers use QR/SVD, not literal (XᵀX)⁻¹ — inverting squares the numerical sensitivity.

Check your understanding