The Matrix Formulation

Topic — Linear regression has a closed-form solution. Not an approximation, not an iterative search — one line of linear algebra that returns the exact best-fit coefficients. This page derives it, verifies it against scikit-learn digit for digit, and then shows the three cases where it collapses.

This is also where the linear algebra earns its place: matrices, transposes, inverses and rank all appear for a concrete reason rather than as background.

---

From equations to matrices

Multiple linear regression with two features and five samples is five equations:

$8 = b_0 + b_1(1) + b_2(2)$

$7 = b_0 + b_1(2) + b_2(1)$

$15 = b_0 + b_1(3) + b_2(4)$

$14 = b_0 + b_1(4) + b_2(3)$

$20 = b_0 + b_1(5) + b_2(5)$

Five equations, three unknowns. Writing them out one by one doesn't scale to 30 features and 20,000 rows, so they are stacked into matrices:

Output

   X (5 samples x 2 features) =
[[1. 2.]
 [2. 1.]
 [3. 4.]
 [4. 3.]
 [5. 5.]]
   y = [ 8.  7. 15. 14. 20.]

The whole system becomes one statement:

$y = X \cdot \beta$

where $X$ holds the features, $\beta$ (beta) is the column of coefficients we want, and $y$ is the target.

---

The column of ones

There is a problem with $y = X \cdot \beta$ as written: the intercept $b_0$ has no feature to multiply. Every other coefficient is attached to a column of $X$; $b_0$ stands alone.

The fix is to give it a column — one filled entirely with $1$s:

Output

   prepend a column of ones -> the design matrix:
[[1. 1. 2.]
 [1. 2. 1.]
 [1. 3. 4.]
 [1. 4. 3.]
 [1. 5. 5.]]
   shape (5, 3): 3 unknowns to solve for (b0, b1, b2)

Now $b_0$ multiplies a column that is $1$ in every row, so it contributes $b_0 \times 1 = b_0$ to every prediction — exactly what an intercept does. The result is called the design matrix, and the intercept has become an ordinary coefficient.

tip

Why this trick matters

It is not cosmetic. Once the intercept is inside $\beta$, there is no special case — one matrix equation solves for the intercept and all the slopes simultaneously. Every formula below works unchanged whether or not there is an intercept.

The same device appears throughout machine learning; in a neural network the equivalent is called the bias term, and it is implemented the same way.

Note the shape bookkeeping: $X$ is $(5, 3)$ and $\beta$ must be $(3,)$ for $X \cdot \beta$ to produce the $(5,)$ predictions. The number of coefficients equals the number of columns, including the ones column.

---

The normal equation

We want the $\beta$ that minimises the sum of squared errors. Because squared error is convex and differentiable, the minimum is where the derivative is zero — and solving that condition gives a formula directly:

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

This is the normal equation. Reading it right to left:

Piece	Shape here	What it is
$X^T$	$(3, 5)$	The transpose — rows become columns
$X^T y$	$(3,)$	Each feature's correlation with the target, unnormalised
$X^T X$	$(3, 3)$	Every feature against every other — the Gram matrix
$(X^T X)^{-1}$	$(3, 3)$	Its inverse
$\beta$	$(3,)$	The coefficients

Worked on our data:

Output

   X^T X =
[[ 5. 15. 15.]
 [15. 55. 53.]
 [15. 53. 55.]]
   X^T y = [ 64. 223. 225.]

   beta = [2.1333 1.2778 2.2778]
   -> b0 = 2.1333, b1 = 1.2778, b2 = 2.2778

Two entries in $X^T X$ are worth naming. The top-left $5$ is the number of samples — it is the ones column dotted with itself. And $15$ appears three times: the sum of each feature column, which is that column dotted with the ones column.

$X^T X$ is always square ($\text{features} \times \text{features}$) and symmetric, since the correlation of feature $i$ with $j$ equals that of $j$ with $i$. That symmetry is what makes the inverse tractable.

Verified against scikit-learn

Output

   scikit-learn: intercept = 2.1333, coef = [1.2778 2.2778]
   identical: True

Digit for digit. LinearRegression is not doing anything mysterious — it is solving this equation, and np.allclose confirms the agreement exactly.

note

What the model is

After fitting, a linear regression is three numbers. $b_0 = 2.1333$, $b_1 = 1.2778$, $b_2 = 2.2778$. Prediction is one dot product.

That is the entire trained model — no stored training data, no tree structure, no iteration state. It's worth contrasting with 1-nearest-neighbour, which stores the whole training set and does all its work at prediction time. Linear regression is the opposite extreme: expensive to fit, trivial to deploy.

---

The Gram Matrix & Collinearity Workbench

Experiment with the interactive matrix calculator below. We blend Column 2 ($x_2$) from our original independent text data ([2, 1, 4, 3, 5]) to perfectly collinear data (2 * x1 = [2, 4, 6, 8, 10]). Watch how the elements of the Gram matrix ($X^T X$), the determinant, and the matrix rank shift in real-time.

Gram Matrix & Collinearity Workbench

Slide from left to right to watch Column 2 transition from perfect collinearity to full independence.

Column Independence Blending: 100% Independent (Original)

Fitted Design Matrix X:

[1, x₁, x₂]

[1.0, 1.0, 2.00]

[1.0, 2.0, 1.00]

[1.0, 3.0, 4.00]

[1.0, 4.0, 3.00]

[1.0, 5.0, 5.00]

Gram Matrix (XᵀX):

[[5.0, 15.0, 15.0],

[15.0, 55.0, 53.0],

[15.0, 53.0, 55.0]]

Determinant: det(XᵀX)

180.0000

Measures close-to-singular state.

Design Matrix Rank:

rank(X) = 3

Must be equal to column count (3) for unique solution.

---

Three ways it breaks

The formula requires $(X^T X)^{-1}$ to exist. A matrix has an inverse only if it is full rank, and there are three ordinary situations in which it isn't.

1. Perfectly correlated features

The previous page showed coefficients destabilising as correlation approached $1.0$. At exactly $1.0$, they stop existing:

Output

   feature 2 = 2 x feature 1, exactly
   det(X^T X) = 0.000e+00
   rank(X) = 2, but X has 3 columns
   np.linalg.inv -> LinAlgError: Singular matrix

The determinant is exactly zero and the matrix is singular. rank(X) = 2 against 3 columns means one column carries no information the others don't already have — so the system has infinitely many solutions, all fitting equally well, and no basis for choosing among them.

This is the endpoint of the collinearity story: at $\text{corr} = 0.999$ the coefficient split was arbitrary but unique; at $\text{corr} = 1.0$ it isn't even unique.

What scikit-learn does instead

Output

   pseudo-inverse   = [-0.0015  0.5991  1.1982]
   scikit-learn     = [-0.0015 0.5991 1.1982]
   no exception — sklearn solves with lstsq/SVD, not an explicit inverse
   the true effect of 3.0 is split across the identical columns:
     0.5991 + 2 x 1.1982 = 2.9954

LinearRegression returned an answer where np.linalg.inv raised. It does not use an explicit inverse — it calls scipy.linalg.lstsq, which uses SVD and returns the minimum-norm solution when there are many: of all the equally-good answers, the one with the smallest coefficients.

And look at the last line. The true effect was $3.0$, and the pseudo-inverse split it as $0.5991 + 2 \times 1.1982 = 2.9954$. The total is recovered exactly even though the individual coefficients are arbitrary — precisely the invariant the previous page measured.

warning

No exception is not the same as no problem

LinearRegression will fit perfectly collinear data and report plausible-looking coefficients. Nothing warns you. The predictions are fine; the coefficients are one arbitrary choice from an infinite set, and any interpretation of them is meaningless.

Check np.linalg.matrix_rank(X) against X.shape[1], or look at the correlation matrix, before reading coefficients from any model.

2. More features than samples

Output

   X shape (3, 5) -> 5 unknowns from 3 equations
   rank(X) = 3: infinitely many exact solutions
   train R2 = 1.0000   <- meaningless

Three equations cannot pin down five unknowns. $X^T X$ is $(5, 5)$ with rank at most $3$, so it is singular for the same reason as above — and the model achieves a perfect $R^2 = 1.0000$ on its training data by passing exactly through every point.

That perfect score carries no information whatsoever. It is the memorisation problem in its purest algebraic form: with $p \ge n$ an exact fit always exists.

This case is common in practice — genomics, text with one column per word, and any wide one-hot encoding. It is one of the main reasons regularisation exists, since a penalty term makes $X^T X + \lambda I$ invertible even when $X^T X$ is not.

3. A constant column

A column with the same value in every row is a multiple of the ones column, so it too destroys independence. That is precisely why the pre-flight audit checked for constant columns — the Region column that was EU in every row would have made the design matrix singular.

Failure	Cause	Detect with
Perfect collinearity	A column is a combination of others	matrix_rank(X) < X.shape[1]
$p \ge n$	Fewer samples than coefficients	X.shape[0] <= X.shape[1]
Constant column	It duplicates the ones column	df[c].nunique() <= 1

---

Why anyone bothers with gradient descent

If there is an exact formula, why does the next section spend its time on an iterative search?

Normal Equation vs. Gradient Descent Cost Comparison

Adjust the slider below to change the number of features $p$ (assuming $n = 100,000$ samples and $100$ gradient descent epochs) to watch the computational operations (FLOPs) of matrix inversion explode.

Normal Equation vs. Gradient Descent Cost

See how the cubic scaling O(p³) of matrix inversion completely crushes the normal equation as dimensionality grows.

Number of Features (p): 100

Normal Equation (O(p³ + np²)):~1.00e+9 FLOPS

Gradient Descent (O(n · p · epochs)):~1.00e+9 FLOPS

At p = 100 features, Normal equation is highly optimal and executes near-instantly.

Two rows decide it. The $p^3$ cost makes the closed form impractical once features run to tens of thousands. And more fundamentally, the closed form exists only for squared error — logistic regression's log loss has no closed-form solution at all, which is why the very next algorithm needs an iterative method.

So the normal equation is the right tool for the problem it solves, and gradient descent is what you use for everything else.

---

2. Implementation Lab

tip

Lab Exercise: The Matrix Formulation

An interactive, fully-functional Google Colab / Jupyter Notebook is available to experiment with this concept live in your browser.

How to run the lab:

1. Click the "Open In Colab" badge above to launch the interactive notebook.
2. Click "Run all" or execute cells individually using Shift + Enter.

Things to try inside the script:

• Experiment 1 (Transpose Math): Adjust the elements of the raw matrix and re-calculate $(X^T X)^{-1} X^T y$ manually to ensure floating point agreement with scikit-learn.
• Experiment 2 (Singularity analysis): Watch OLS crash during a perfect $r=1.0$ linear dependency, and observe how scikit-learn's internal SVD solvers bypass the exception.
• Experiment 3 (Dimension trap): Fit a model with fewer samples than features ($p \ge n$) and observe the empty $R^2 = 1.0000$ training score.

---

3. Common Mistakes

1. Forgetting Ones Column

• The Misconception: "I can just pass raw features X directly to the normal equation."
• Wrong Thinking: Believing the equation automatically fits an intercept without a ones column.
• The Right Principle: Without the ones column, the intercept is forced to $0$, meaning the fitted plane is mathematically constrained to pass through the origin $(0,0,\dots)$. Always prepend a column of $1$s.

2. No error = No problem

• The Misconception: "Scikit-Learn fitted the model without raising any error, so the coefficients are completely trustworthy."
• Wrong Thinking: Believing a completed fit implies there are no collinearities or singular matrices.
• The Right Principle: Scikit-learn does not solve the equation via inv(XtX). It uses SVD minimum-norm solutions, meaning it will happily fit perfectly collinear data and output arbitrary coefficients without warning you. Always check matrix rank beforehand.

3. Computing inv(XtX) explicitly

• The Misconception: "I must write inv(X.T @ X) in my OLS numpy code."
• Wrong Thinking: Assuming explicit matrix inversion is numerically stable.
• The Right Principle: Explicit matrix inversion is highly unstable. In production or scratch code, always use np.linalg.solve(XtX, Xty) or SVD solvers (pinv / lstsq) for superior numerical stability.

4. Closed form universal applicability

• The Misconception: "The closed-form normal equation can solve any linear regression variant, like Logistic Regression."
• Wrong Thinking: Attempting to solve classification cost functions using matrix inversions.
• The Right Principle: The closed-form normal equation exists only for squared error (because its loss surface is perfectly quadratic/convex and has a single global minimum). Logistic regression uses logarithmic loss, which has no analytical closed-form solution and must be solved iteratively.

---

4. Summary

Core Theory

• Design Matrix: The raw feature matrix $X$ with a column of $1$s prepended, merging the intercept as an ordinary coefficient.
• Normal Equation: $\beta = (X^T X)^{-1} X^T y$ is the unique analytical solution derived by setting the SSE derivative to zero.
• Gram Matrix: $X^T X$ is a square, symmetric matrix showing every feature against every other. Its top-left entry counts rows.
• Full Rank Rule: OLS requires $X$ to be full rank. Rank deficiency destroys matrix invertibility.

Engineering Realities

• LinearRegression internal: Uses SVD/lstsq solvers matching the pseudo-inverse (pinv) rather than standard OLS inversion (inv).
• Unwarned collinearity: Scikit-learn bypasses singular matrix errors, outputting arbitrary weights under perfect collinearity without alerting the user.
• $p \ge n$ Trap: Having fewer samples than features guarantees a meaningless $R^2 = 1.0000$ training score due to memorization.
• Cubic Inversion scaling: Inverting $X^T X$ scales as $O(p^3)$, completely ruling out the normal equation for wide datasets where gradient descent scales as $O(p)$.

info

Key Takeaways

• Prediction overhead: Once fit, a linear regression is merely a vector of coefficients. Prediction requires a single dot product—zero stored training data, zero iterative overhead.
• Regularisation fix: Under the $p \ge n$ genomics/text wide-data trap, L2 regularization (adding $\lambda I$ to $X^T X$) forces the singular matrix to become invertible.
• Solver choice: Use the Normal Equation for exact, rapid solutions on narrow datasets. Use Gradient Descent for wide, massive, or streaming datasets, or for non-squared-error losses.

Next in this section: Loss Functions — the closed form exists only for squared error, so the next section begins by asking what else you might minimise, and what that choice costs.

---

5. Active Recall

Attempt each question first, then click to reveal detailed answers, mathematical derivations, and sample explanations.

Interactive Checkpoint Quiz

Interactive Checkpoint: The Pseudo-Inverse Partition

Under perfect collinearity (Column 3 = 2 × Column 2), the true coefficient effect of the features is exactly 3.0. Standard inversion crashes. What does Scikit-Learn's SVD solver return?

Review Flashcards

1. [THEORY] Why do we prepend a column of ones to $X$, and what does this do algebraically to $b_0$?

• The Intercept Trick: The linear equation $y = b_0 + b_1 x_1 + \dots$ has $b_0$ standing alone.
• By prepending a column of $1$s to the feature matrix $X$ (forming the design matrix), $b_0$ now multiplies a feature value of $1.0$ in every single row.
• This incorporates the intercept as an ordinary coefficient inside the parameter vector $\beta$: $\beta = \begin{bmatrix} b_0 & b_1 & \dots & b_k \end{bmatrix}^T$
• The Benefit: This eliminates the intercept as an algebraic special case. The same single matrix equation solves for both the intercept and all slope coefficients simultaneously.

2. [THEORY] Why does the closed-form normal equation exist ONLY for squared error loss?

• The Derivation: The Sum of Squared Errors loss is: $L = \sum (y_i - X_i \beta)^2 = (y - X\beta)^T (y - X\beta)$
• Because squared error is a quadratic function, its derivative with respect to $\beta$ is perfectly linear: $\frac{\partial L}{\partial \beta} = -2 X^T(y - X\beta) = 0 \implies X^T X \beta = X^T y$
• Because this first-order condition is perfectly linear, we can solve for $\beta$ analytically via a single matrix inversion: $\beta = (X^T X)^{-1} X^T y$
• Other Losses: For other losses (like absolute error or logistic log-loss), the loss derivative is not linear, meaning there is no linear system to invert and a closed-form analytical solution cannot exist.

3. [ANALYZE] Why is a training score of $R^2 = 1.0000$ mathematically guaranteed when $p \ge n$?

• The Geometry: In $n$-dimensional space, $n$ data points define a specific set of coordinates.
• If the number of features (dimension parameters) $p$ is greater than or equal to $n$, we have more degrees of freedom (unknown variables) than coordinate constraints (equations).
• In linear algebra, a system of linear equations with more unknowns than equations has infinitely many exact solutions.
• This guarantees that a hyperplane can always be positioned to pass exactly through every single training point, yielding a perfect, yet completely overfitted and meaningless, $R^2 = 1.0000$.

4. [ANALYZE] Why can't we use standard OLS matrix inversion on perfectly collinear data, and how does Scikit-Learn resolve it?

• Standard Inversion Fail: Under perfect collinearity, one column is a scalar multiple of another, meaning the design matrix columns are not linearly independent. This makes $X^T X$ singular, with a determinant of exactly $0.00$, making standard OLS inversion $(X^T X)^{-1}$ fail with a Singular matrix division error.
• The Scikit-Learn Solver: Scikit-learn does not call standard matrix inversion inv. It uses Singular Value Decomposition (SVD) and calculates the Moore-Penrose pseudo-inverse ($X^+$), outputting the minimum-norm solution. While individual coefficients are arbitrary, they are numerically stable and the combined total effect is recovered perfectly.

5. [ANALYZE] At $p = 10,000$ features, why does the normal equation completely collapse in terms of computational scaling compared to Gradient Descent?

• Matrix Inversion Scaling: Computing $(X^T X)^{-1}$ requires a matrix inversion that scales cubically as $O(p^3)$.
• At $p = 100$, $p^3 = 10^6$ operations (instantaneous).
• At $p = 10,000$, $p^3 = 10^{12}$ (1 trillion operations), requiring extreme processing power and causing severe RAM/memory overheads.
• Gradient Descent Scaling: Gradient descent requires calculating the gradient vector which scales linearly as $O(n \cdot p)$ per iteration. For high-dimensional datasets ($p \ge 10,000$), gradient descent's linear complexity is the only computationally feasible approach.