Ch.6: Multiple Linear Regression - The Math

Inspired by: YouTube

Ch.1 and Ch.2 introduced simple linear regression and derived its closed-form OLS solution for a single input feature. Ch.3 established how to measure regression model performance using metrics like MAE, MSE, and R². Ch.4 covered the assumptions a fitted line needs to roughly satisfy to be trustworthy. Ch.5 extended the intuition to multiple dimensions, fitting a plane or hyperplane through data with multiple input columns.

This post completes the picture for multiple linear regression. Trying to calculate partial derivatives scalar-by-scalar for $k$ features quickly becomes unmanageable. Matrix algebra provides a clean way to solve for all coefficients at once. Below is the full derivation of the closed-form normal equation, an implementation built from scratch in NumPy, and a side-by-side comparison against scikit-learn on the diabetes dataset.

Matrix Form of the Equation

In simple linear regression, the model equation is $y = m x + b$. When moving to multiple input features $x_1, x_2, \dots, x_k$, the prediction for a single row $i$ becomes:

ŷᵢ = β₀ + β₁·xᵢ₁ + β₂·xᵢ₂ + ... + βₖ·xᵢₖ

For a dataset with $n$ rows and $k$ features, writing $n$ individual equations is tedious. Matrix notation packs all of them into a single line.

First, prepend a column of 1s to the feature matrix $X$. Multiplying that all-ones column by $\beta_0$ contributes exactly $\beta_0$ to every single row's prediction, no matter what the row's other features are, which is exactly what an intercept is supposed to do. It's the matrix equivalent of writing $y = mx + b \cdot 1$ instead of $y = mx + b$: same equation, but now every term is "some coefficient times some column," which is what lets the whole thing collapse into one matrix multiplication. The augmented matrix $X$ has shape $(n, k+1)$:

Next, represent all target values in a vector $y$ of shape $(n, 1)$, and all parameters in a coefficient vector $\beta$ of shape $(k+1, 1)$:

With these definitions, the predictions for every single row in the dataset are computed simultaneously via matrix multiplication:

ŷ = X · β

Where $\hat{y}$ is a column vector of shape $(n, 1)$.

Formulating the Cost Function

Just as in simple linear regression, the goal is to minimize the Sum of Squared Errors across all $n$ samples. Define the error vector $e$ as the difference between actual targets and predictions:

e = y - ŷ = y - X · β

The sum of squared errors $E(\beta)$ is equal to the inner product of the error vector with itself, $e^T e$:

Why $e^Te$ gives the sum of squares: Each entry of $e$ is just a residual, $e_i = y_i - \hat{y}_i$. Take a 2-row example: $e^T$ is the row $\begin{bmatrix} y_1 - \hat{y}_1 & y_2 - \hat{y}_2 \end{bmatrix}$, and $e$ is the column $\begin{bmatrix} y_1 - \hat{y}_1 \\ y_2 - \hat{y}_2 \end{bmatrix}$. Standard matrix multiplication pairs up each matching entry and adds the products: $e^T e = (y_1 - \hat{y}_1)(y_1 - \hat{y}_1) + (y_2 - \hat{y}_2)(y_2 - \hat{y}_2) = (y_1 - \hat{y}_1)^2 + (y_2 - \hat{y}_2)^2$. Since every entry is being paired with itself, the product comes out to a sum of squared residuals for free, no extra summation notation needed.

E(β) = eᵀ · e = (y - X·β)ᵀ · (y - X·β)

Expanding the transpose $(A - B)^T = A^T - B^T$:

E(β) = (yᵀ - βᵀ·Xᵀ) · (y - X·β)

Distributing terms across the product:

E(β) = yᵀ·y - yᵀ·X·β - βᵀ·Xᵀ·y + βᵀ·Xᵀ·X·β

Look closely at the middle two terms: $y^T X \beta$ and $\beta^T X^T y$.

The vector dimensions are $y^T$ of shape $(1, n)$, $X$ of shape $(n, k+1)$, and $\beta$ of shape $(k+1, 1)$. Their product $y^T X \beta$ is a $1 \times 1$ scalar.

Since the transpose of a scalar is itself, taking the transpose of $y^T X \beta$ yields:

(yᵀ · X · β)ᵀ = βᵀ · Xᵀ · (yᵀ)ᵀ = βᵀ · Xᵀ · y

Because both middle terms represent the exact same scalar value, they can be combined:

E(β) = yᵀ·y - 2·βᵀ·Xᵀ·y + βᵀ·Xᵀ·X·β

This single expression is the cost function in matrix form.

Minimizing the Cost Function

$E(\beta)$ is a quadratic function of $\beta$: expand it out and every term is either constant, linear in $\beta$, or $\beta$ multiplied by itself. In one dimension, a quadratic like $E(\beta) = a\beta^2 - 2b\beta + c$ traces out a parabola, a single bowl shape with one lowest point, and calculus finds that point by setting $\frac{dE}{d\beta} = 0$. With $k+1$ coefficients instead of one, $E(\beta)$ traces out a $(k+2)$-dimensional bowl (a paraboloid) instead of a 2D parabola, but it is still convex: one global minimum, no other flat spots to get stuck at. Setting every partial derivative to zero at once, i.e. the full gradient $\frac{\partial E}{\partial \beta} = 0$, finds that single lowest point directly, the same way ordinary calculus finds the bottom of a parabola.

Differentiating each term in the cost function with respect to $\beta$, term by term:

1. $\frac{\partial}{\partial \beta}(y^T y) = 0$, since $y^T y$ does not depend on $\beta$ at all, it's a constant with respect to what we're differentiating, same as how $\frac{d}{dx}(5) = 0$ in ordinary calculus.
2. $\frac{\partial}{\partial \beta}(-2 \beta^T X^T y) = -2 X^T y$. This term is linear in $\beta$ (no $\beta$ multiplied by itself), so it behaves like differentiating $-2a\beta$ with respect to $\beta$ in scalar calculus, which gives the constant $-2a$ back. Here that "constant" is the vector $X^T y$. Formally this uses the matrix derivative identity $\frac{\partial}{\partial x}(x^T A) = A^T$.
3. $\frac{\partial}{\partial \beta}(\beta^T X^T X \beta) = 2 X^T X \beta$. This term is quadratic in $\beta$ (it appears twice), the matrix analog of $\frac{d}{dx}(ax^2) = 2ax$ in scalar calculus. Formally this uses the identity $\frac{\partial}{\partial x}(x^T A x) = 2 A x$ for symmetric matrices $A = X^T X$.

Putting these derivatives together:

Dividing by 2 and moving terms gives:

Xᵀ·X·β = Xᵀ·y

This matrix equation is known as the Normal Equation.

Assuming $X^T X$ is invertible (non-singular), multiply both sides by $(X^T X)^{-1}$:

β = (Xᵀ · X)⁻¹ · Xᵀ · y

That is the complete closed-form solution. A single matrix expression solves for every coefficient, including the intercept, in one step.

Computational Complexity and Limitations

The closed-form OLS formula is elegant, but it comes with a practical computational trade-off. Computing $(X^T X)^{-1}$ requires inverting a matrix of shape $(k+1, k+1)$, where $k$ is the number of features. Using Gaussian elimination, matrix inversion has a computational complexity of roughly $O(k^3)$. With 10 or even 100 features this is instantaneous, but for image data, text data, or any dataset with a very large number of columns, computing $k^3$ operations becomes extremely slow.

That's the motivation for Gradient Descent: an approximation technique that arrives at the coefficients iteratively instead of via a single formula. It doesn't produce the exact same answer as the normal equation, but it lands close, and it stays fast regardless of feature count. Most of the time sklearn's LinearRegression class (which uses the closed-form solution) is enough, but SGDRegressor is available when gradient descent is the better fit.

Coding It From Scratch

The normal equation $\beta = (X^T X)^{-1} X^T y$ translates into a compact Python class using NumPy:

import numpy as np
 
class MeraLR:
 
    def __init__(self):
        self.coef_ = None
        self.intercept_ = None
 
    def fit(self, X_train, y_train):
        # Step 1: Prepend a column of 1s for the intercept
        X_train = np.insert(X_train, 0, 1, axis=1)
        
        # Step 2: Calculate betas using the Normal Equation
        # beta = (X^T * X)^-1 * X^T * y
        betas = np.linalg.inv(np.dot(X_train.T, X_train)).dot(X_train.T).dot(y_train)
        
        # Step 3: Extract intercept (beta_0) and feature coefficients (beta_1 to beta_k)
        self.intercept_ = betas[0]
        self.coef_ = betas[1:]
 
    def predict(self, X_test):
        y_pred = np.dot(X_test, self.coef_) + self.intercept_
        return y_pred

To verify the implementation, compare MeraLR against scikit-learn's LinearRegression using the standard diabetes dataset (10 input features):

from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
 
# Load diabetes dataset (10 features)
X, y = load_diabetes(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=2)
 
# Benchmark 1: sklearn LinearRegression
sk_lr = LinearRegression()
sk_lr.fit(X_train, y_train)
y_pred_sk = sk_lr.predict(X_test)
 
print("--- Sklearn LinearRegression ---")
print("Intercept:", sk_lr.intercept_)
print("Coefs:", sk_lr.coef_)
print("R2 score:", r2_score(y_test, y_pred_sk))
 
# Benchmark 2: Custom MeraLR
mera_lr = MeraLR()
mera_lr.fit(X_train, y_train)
y_pred_mera = mera_lr.predict(X_test)
 
print("\n--- Custom MeraLR ---")
print("Intercept:", mera_lr.intercept_)
print("Coefs:", mera_lr.coef_)
print("R2 score:", r2_score(y_test, y_pred_mera))

Running this code produces the following empirical output:

--- Sklearn LinearRegression ---
Intercept: 151.88331005254165
Coefs: [  -9.15865318 -205.45432163  516.69374454  340.61999905 -895.5520019
  561.22067904  153.89310954  126.73139688  861.12700152   52.42112238]
R2 score: 0.439933866156897

--- Custom MeraLR ---
Intercept: 151.8833100525417
Coefs: [  -9.15865318 -205.45432163  516.69374454  340.61999905 -895.5520019
  561.22067904  153.89310954  126.73139688  861.12700152   52.42112238]
R2 score: 0.4399338661568968

Both implementations match to 14 decimal places. The custom class derived strictly from matrix calculus produces the exact same predictions, intercept, feature coefficients, and $R^2$ score ($0.4399$) as scikit-learn's built-in model.

What's Next

The closed-form normal equation works well when the number of features is manageable, but its $O(k^3)$ computational cost becomes a bottleneck on large, high-dimensional datasets. The next topic in optimization addresses this limitation: Gradient Descent, an iterative approach that finds optimal regression parameters step-by-step without computing matrix inverses.