Everything in this Part, simple regression, multiple regression, even the linear layer of a neural network, is the same calculation underneath: stack your data into a matrix, and solve one equation. That equation is β̂ = (XᵀX)⁻¹Xᵀy, and it is worth understanding, not just calling.
This chapter is the theory behind step 2 (fit the model) of the five-step workflow. And the Gauss-Markov conditions that make OLS trustworthy are the same LINE conditions you check in steps 1 and 4, seen from the estimator's side.
From a Line to a Matrix
With one predictor you can write two normal equations by hand. With twenty predictors you need linear algebra. Stack every observation's response into a vector y and its predictors into the design matrix X (that leading column of 1s is what gives you an intercept), and the whole model collapses to y = Xβ + ε.
Minimizing the squared residuals leads to the normal equations XᵀXβ̂ = Xᵀy, and as long as the columns of X are not redundant (so XᵀX can be inverted) the solution is unique:
That single line is the entire fitting step, for one predictor or a thousand. When two columns of X are nearly redundant (multicollinearity), XᵀX is nearly singular and the inverse blows up, which is why collinearity makes coefficients unstable, a problem Chapters 94 and 97 return to.
The Geometry of Least Squares
There is a picture that makes OLS click. Think of the response y as a single point in n-dimensional space. Every possible set of fitted values Xβ traces out a flat subspace, the column space of X. OLS finds the point in that subspace closest to y.
The closest point is the perpendicular projection of y onto the subspace. That is the fitted vector ŷ, and the leftover e = y − ŷ is the residual, which sits at a right angle to the subspace. "Perpendicular" is exactly the algebraic condition Xᵀe = 0: the residuals are uncorrelated with every predictor. The map from y to ŷ is the hat matrix H = X(XᵀX)⁻¹Xᵀ, so ŷ = Hy, it literally puts the hat on y.
Gauss-Markov, BLUE, and When It Fails
Why least squares, and not some other rule? The Gauss-Markov theorem gives the answer: under a short list of conditions, OLS is the Best Linear Unbiased Estimator (BLUE). Among all estimators that are linear in y and unbiased, none has smaller variance. "Best" means most precise.
The conditions are the LINE assumptions from the Simple Linear Regression chapter, wearing their estimator-theory names. And critically, when a condition fails, OLS does not simply break, it loses a specific guarantee, which tells you exactly how to patch it:
| Gauss-Markov condition | LINE name | If it fails | Remedy |
|---|---|---|---|
| Linear in parameters | Linearity | estimates are biased | transform, add terms, re-specify (see Regularization & Flexible Models) |
| Exogeneity E[ε|X]=0 | (correct model) | estimates are biased | add omitted variables, instruments (see Econometrics & Panel Data) |
| Homoscedasticity | Equal variance | still unbiased, but SEs wrong, not BLUE | robust (HC) SEs, or weighted/GLS |
| No autocorrelation | Independence | still unbiased, but SEs wrong, not BLUE | Newey-West SEs, time-series models (see Econometrics & Panel Data) |
The pattern is worth memorizing. Break linearity or exogeneity and the coefficients themselves are wrong. Break equal variance or independence and the coefficients are still on target, but their standard errors, and therefore every p-value and confidence interval, are not to be trusted. The fix for that second case is usually painless: heteroscedasticity-robust standard errors, one argument in statsmodels.
Real-World Example: What Predicts a Wage?
We fit hourly wage on education and experience, first to confirm the matrix
formula reproduces statsmodels exactly, then to catch a Gauss-Markov violation and fix it.
One row per worker with education (years), experience
(years), and hourly wage (dollars).
Fit (step 2). Computing β̂ = (XᵀX)⁻¹Xᵀy by hand gives the identical coefficients statsmodels returns, the formula is not a black box:
| Coefficient | Estimate | Plain-language reading |
|---|---|---|
| Intercept | $2.42 | anchors the plane (no schooling, no experience is extrapolation) |
| education | +$2.38 / year | each additional year of schooling adds about 2.38 dollars/hour |
| experience | +$0.32 / year | each year of experience adds about 0.32 dollars/hour |
| R² | 0.57 | the two predictors explain 57% of wage variation |
Predict. The fitted plane turns into a prediction by plugging values in:
A worker with a 4-year degree (16 years) and 10 years of experience: 2.42 + 2.376×16 + 0.319×10 ≈ $43.60/hour. Across the whole sample the predictions track the actual wages closely:
Check the conditions (step 4). The residual spread grows with education, higher-educated workers
have far more variable pay. A Breusch-Pagan test confirms it (p < 0.001): homoscedasticity fails. By
Gauss-Markov, the coefficients above are still unbiased, but the classical standard errors are not reliable. The
remedy is one keyword, cov_type="HC3":
| Standard error on education | Value | Verdict |
|---|---|---|
| Classical OLS SE | 0.127 | assumes equal variance, which is violated here |
| HC3 robust SE | 0.134 | valid under heteroscedasticity, the one to report |
Here the two are close, so the conclusion (education clearly matters) is unchanged, but the discipline is the point: diagnose first, then report the standard errors your data actually justify. When the gap is large, reporting the classical SE would overstate your certainty.
OLS in Machine Learning & AI
The normal equations are everywhere in ML, sometimes solved exactly, more often approached by gradient descent.
| Idea (this chapter) | In ML / AI it becomes | Example |
|---|---|---|
| β̂ = (XᵀX)⁻¹Xᵀy | The closed-form (normal-equation) solver | LinearRegression in scikit-learn |
| Minimizing squared error | Gradient descent on an MSE loss | how models too big to invert XᵀX are trained |
| (XᵀX)⁻¹ blows up | Ill-conditioning from collinear features | ridge adds λI to stabilize the inverse (see Regularization & Flexible Models) |
| Hat matrix / projection | Leverage & influence diagnostics | finding points that dominate the fit |
| OLS under Normal errors | Maximum likelihood estimation | least squares = MLE when errors are Gaussian |
The linear layer at the output of nearly every deep network is doing regression on learned features, and its weights are found by minimizing squared (or cross-entropy) loss, the same objective as OLS, just solved by gradient descent instead of a matrix inverse. The conditioning of XᵀX is why feature scaling and regularization matter: correlated or badly-scaled inputs make the loss surface a stretched valley that is slow and unstable to descend, exactly the finite-sample face of a near-singular XᵀX. And the fact that OLS equals maximum likelihood under Gaussian noise is the bridge from this chapter to the probabilistic view of learning, where a loss function is a negative log-likelihood and the choice of loss encodes an assumption about the noise.
See the matrix algebra come alive
The companion notebook computes β̂ = (XᵀX)⁻¹Xᵀy with a few NumPy lines and
confirms it matches statsmodels to the last digit, builds the hat matrix and verifies the
residuals are orthogonal to X, runs a Monte-Carlo Gauss-Markov demo showing OLS is the
tightest unbiased estimator, and then diagnoses heteroscedasticity in the-ols-framework--wages.xlsx and
switches to robust standard errors.
View opens the rendered notebook instantly (no setup). Open in Colab runs &
edits it live in your browser. To run locally, install numpy, pandas, scipy,
matplotlib, seaborn, statsmodels, and openpyxl and launch
jupyter notebook.
🎓 Key Takeaways
- ✓Any regression is y = Xβ + ε, fit by the normal equations with the closed form β̂ = (XᵀX)⁻¹Xᵀy.
- ✓Least squares is a projection: ŷ is the perpendicular projection of y onto the column space of X, so the residual satisfies Xᵀe = 0.
- ✓Gauss-Markov → BLUE: under its conditions OLS is the minimum-variance linear unbiased estimator.
- ✓Know the failure modes: broken linearity/exogeneity bias the coefficients; broken equal-variance/independence only spoil the standard errors, fix those with robust SEs or GLS.
- ✓Real data: the wage model's hand-computed β̂ matches statsmodels exactly; Breusch-Pagan flags heteroscedasticity, so HC3 robust SEs are the ones to report.
Practice Challenges
Five short challenges. Use NumPy for the matrix algebra and statsmodels to check yourself.
Fit by matrix, then by library
Build X with a column of 1s and compute β̂ = (XᵀX)⁻¹Xᵀy; confirm it matches ols(...).fit().
np.linalg.solve(X.T@X, X.T@y) is more stable than inverting.Orthogonality of residuals
Show Xᵀe ≈ 0 and that total = explained + residual sum of squares.
X.T @ e and the SST/SSR/SSE identity.The hat matrix
Form H = X(XᵀX)⁻¹Xᵀ; verify ŷ = Hy and that trace(H) equals the number of parameters.
Gauss-Markov by simulation
Across many simulated samples, compare the variance of OLS with a simple alternative unbiased estimator.
Robust standard errors
Fit the wage model, run Breusch-Pagan, then compare classical and HC3 standard errors and CIs.
.fit(cov_type="HC3") versus the default.A fully-worked solutions notebook walks through all five challenges, each verified in code. Try them yourself first, then compare.
Quiz: Test Yourself
Eight quick questions on the OLS framework. Answer them, hit Check Answers, and keep refining until you score 100%. Your progress is saved, so you can hop back to the chapter and return anytime.
You have the engine: one matrix formula, a geometric picture, and the conditions that make it optimal. Multiple Linear Regression puts many predictors to work, and shows how to read a coefficient when the others are "held constant."