What you’ll learn

• the matrix form of the model, and how to derive the Normal Equation by setting a gradient to zero
• how to solve a three-parameter fit by hand when the design is orthogonal
• the exact meaning of “holding all other features constant”, and when that phrase is a lie
• why coefficients are only comparable after standardising
• what multicollinearity breaks — and, just as importantly, what it does not
• why scikit-learn uses lstsqlstsq rather than inverting a matrix

Intuition

One feature gives you a line. Two features give you a plane. Ten features give you a hyperplane nobody can picture — but the arithmetic never changes: multiply each feature by its weight, add them all up, add the bias.

What does change is interpretation. With one feature, the slope is simply “how $y$ moves with $x$”. With several correlated features, each coefficient answers a much narrower question: how $y$ moves with this feature once the others have already had their say. Miss that and you will misread almost every model you fit.

  flowchart LR
  X1["x1  size"] --> M["Weighted sum
theta transpose x"]
  X2["x2  bedrooms"] --> M
  X3["x3  age"] --> M
  XN["x_n  ..."] --> M
  M --> Y["Prediction ŷ"]

The math

The model in matrix form

Stack $m$ instances as rows of a design matrix $\mathbf{X} \in \mathbb{R}^{m\times(n+1)}$, whose first column is all ones:

Deriving the Normal Equation

Expand the squared norm as an inner product and drop the constant $1/m$ — it cannot move the minimiser:

Differentiate with respect to $\boldsymbol{\theta}$ using two standard identities — $\nabla_{\boldsymbol{\theta}}(\boldsymbol{\theta}^\top\mathbf{A}\boldsymbol{\theta}) = 2\mathbf{A}\boldsymbol{\theta}$ for symmetric $\mathbf{A}$, and $\nabla_{\boldsymbol{\theta}}(\boldsymbol{\theta}^\top\mathbf{b}) = \mathbf{b}$:

That is the Normal Equation: one formula, no iteration, no learning rate.

The residual is orthogonal to every feature

Rearranged, the Normal Equation says $\mathbf{X}^\top(\mathbf{y} - \mathbf{X}\boldsymbol{\theta}) = \mathbf{0}$: the residual vector is perpendicular to every column of $\mathbf{X}$. Two consequences you can check on any fit:

• because the first column is all ones, the residuals sum to zero
• because every other column is a feature, the residuals are uncorrelated with each feature

If either fails, you did not fit least squares.

Worked example by hand

An orthogonal design — a small experiment where two on/off factors were varied deliberately. Both are coded $-1$ (off) and $+1$ (on), which keeps the arithmetic exact.

Step 1 — build $\mathbf{X}^\top\mathbf{X}$. Every column holds four entries of $\pm 1$, and any two different columns are orthogonal, so all off-diagonal entries vanish:

Step 2 — build $\mathbf{X}^\top\mathbf{y}$.

Step 3 — multiply.

So $\hat{y} = 4.25 + 2.25\,x_1 + 1.25\,x_2$.

Step 4 — check the fit.

Residuals sum to zero, and they are orthogonal to both feature columns — $\sum x_1 e = -0.25 + 0.25 - 0.25 + 0.25 = 0$, and the same for $x_2$. Exactly what the Normal Equation promised.

Reading the coefficients. Because the features are coded $\pm 1$, switching the promotion on moves sales by $2 \times 2.25 = 4.5$ units, and a weekend adds $2 \times 1.25 = 2.5$. The promotion is the stronger lever, and here that comparison is legitimate — both features share a scale and are uncorrelated.

See it move

The Normal Equation solves for the plane in one shot. Below, a flat baseline plane tilts into the least-squares fit over two features; the red stems are residuals shrinking as it locks in.

From scratch

Three ways to solve the same system, in increasing order of how much you should trust them:

import numpy as np
 
X = np.array([[1, -1, -1],
              [1, -1,  1],
              [1,  1, -1],
              [1,  1,  1]], dtype=float)
y = np.array([1.0, 3.0, 5.0, 8.0])
 
# 1. Textbook form — never use this in production
theta_inv = np.linalg.inv(X.T @ X) @ X.T @ y
 
# 2. Solve the linear system instead of inverting: faster, better conditioned
theta_solve = np.linalg.solve(X.T @ X, X.T @ y)
 
# 3. What scikit-learn actually does: least squares via SVD.
#    Still returns an answer when X.T @ X is singular.
theta_lstsq, *_ = np.linalg.lstsq(X, y, rcond=None)
 
print(theta_inv.round(4))     # [4.25 2.25 1.25]
print(theta_solve.round(4))   # [4.25 2.25 1.25]
print(theta_lstsq.round(4))   # [4.25 2.25 1.25]
 
resid = y - X @ theta_lstsq
print(resid.round(4))                 # [ 0.25 -0.25 -0.25  0.25]
print(round(resid.sum(), 10))         # 0.0
print((X.T @ resid).round(10))        # [0. 0. 0.]  <- orthogonal to every column

import numpy as np
 
X = np.array([[1, -1, -1],
              [1, -1,  1],
              [1,  1, -1],
              [1,  1,  1]], dtype=float)
y = np.array([1.0, 3.0, 5.0, 8.0])
 
# 1. Textbook form — never use this in production
theta_inv = np.linalg.inv(X.T @ X) @ X.T @ y
 
# 2. Solve the linear system instead of inverting: faster, better conditioned
theta_solve = np.linalg.solve(X.T @ X, X.T @ y)
 
# 3. What scikit-learn actually does: least squares via SVD.
#    Still returns an answer when X.T @ X is singular.
theta_lstsq, *_ = np.linalg.lstsq(X, y, rcond=None)
 
print(theta_inv.round(4))     # [4.25 2.25 1.25]
print(theta_solve.round(4))   # [4.25 2.25 1.25]
print(theta_lstsq.round(4))   # [4.25 2.25 1.25]
 
resid = y - X @ theta_lstsq
print(resid.round(4))                 # [ 0.25 -0.25 -0.25  0.25]
print(round(resid.sum(), 10))         # 0.0
print((X.T @ resid).round(10))        # [0. 0. 0.]  <- orthogonal to every column

With scikit-learn

from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
 
data = load_diabetes()
X_train, X_test, y_train, y_test = train_test_split(
    data.data, data.target, test_size=0.3, random_state=42
)
 
model = LinearRegression().fit(X_train, y_train)
 
print(f"train R^2 = {model.score(X_train, y_train):.4f}")   # train R^2 = 0.5244
print(f"test  R^2 = {model.score(X_test, y_test):.4f}")     # test  R^2 = 0.4773
 
rmse = mean_squared_error(y_test, model.predict(X_test)) ** 0.5
print(f"test RMSE = {rmse:.2f}")                            # test RMSE = 53.12

from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
 
data = load_diabetes()
X_train, X_test, y_train, y_test = train_test_split(
    data.data, data.target, test_size=0.3, random_state=42
)
 
model = LinearRegression().fit(X_train, y_train)
 
print(f"train R^2 = {model.score(X_train, y_train):.4f}")   # train R^2 = 0.5244
print(f"test  R^2 = {model.score(X_test, y_test):.4f}")     # test  R^2 = 0.4773
 
rmse = mean_squared_error(y_test, model.predict(X_test)) ** 0.5
print(f"test RMSE = {rmse:.2f}")                            # test RMSE = 53.12

Ten features lift test $R^2$ from the 0.34 that BMI alone managed to 0.48. The gap between train (0.52) and test (0.48) is small, which is what a well-behaved linear model looks like.

figurePredicted against actual on held-out datamatplotlib

Reading the plot

• No systematic curve. Points sit around the diagonal rather than bowing above or below it, so a linear model is a reasonable shape here.
• Compression at the extremes. Low actual values are over-predicted and high ones under-predicted — the classic signature of a model that has captured the trend but not the tails.
• Vertical spread is the error. At any actual value the scatter is roughly ±50, which matches the RMSE of 53.

Interpreting coefficients

Raw coefficients are not comparable

A coefficient is expressed in target units per one unit of that feature. A coefficient of 949 on BMI and 3 on age says nothing about relative importance until both features share a scale. Standardise first, and each coefficient becomes “target change per one standard deviation”:

from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
 
data = load_diabetes()
X = StandardScaler().fit_transform(data.data)
model = LinearRegression().fit(X, data.target)
 
for name, coef in sorted(zip(data.feature_names, model.coef_),
                         key=lambda pair: -abs(pair[1])):
    print(f"{name:>5} {coef:8.2f}")
 
#    s1   -37.68
#    s5    35.73
#   bmi    24.73
#    s2    22.68
#    bp    15.43
#   sex   -11.41
#    s4     8.42
#    s3     4.81
#    s6     3.22
#   age    -0.48

from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
 
data = load_diabetes()
X = StandardScaler().fit_transform(data.data)
model = LinearRegression().fit(X, data.target)
 
for name, coef in sorted(zip(data.feature_names, model.coef_),
                         key=lambda pair: -abs(pair[1])):
    print(f"{name:>5} {coef:8.2f}")
 
#    s1   -37.68
#    s5    35.73
#   bmi    24.73
#    s2    22.68
#    bp    15.43
#   sex   -11.41
#    s4     8.42
#    s3     4.81
#    s6     3.22
#   age    -0.48

figureStandardised coefficients, sortedmatplotlib

“Holding everything else constant”

The textbook reading of $\theta_j$ is: the change in $\hat{y}$ for a one-unit increase in $x_j$, with all other features held fixed. That clause is doing enormous work. If two features always move together in the real world, no observation ever shows one changing while the other stays put — so the coefficient describes a scenario your data never contained.

That is why ageage scores −0.48 here despite age plainly mattering to disease progression: its information is already carried by blood pressure and the serum measurements. The coefficient does not say “age does not matter”, it says “age adds nothing beyond the other nine columns”.

Multicollinearity

figure120 refits on fresh samples from the same processmatplotlib

Detect it with the variance inflation factor. For feature $j$, regress it on all the others and take:

A VIF of 1 means no collinearity, above 5 is worth a look, and above 10 means that coefficient is close to meaningless on its own.

import numpy as np
from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
 
data = load_diabetes()
X = data.data
 
for j, name in enumerate(data.feature_names):
    others = np.delete(X, j, axis=1)
    r2 = LinearRegression().fit(others, X[:, j]).score(others, X[:, j])
    vif = 1 / (1 - r2)
    flag = "  <-- high" if vif > 10 else ""
    print(f"{name:>5}  VIF = {vif:7.2f}{flag}")
 
#   age  VIF =    1.22
#   sex  VIF =    1.28
#   bmi  VIF =    1.51
#    bp  VIF =    1.46
#    s1  VIF =   59.20  <-- high
#    s2  VIF =   39.19  <-- high
#    s3  VIF =   15.40  <-- high
#    s4  VIF =    8.89
#    s5  VIF =   10.08  <-- high
#    s6  VIF =    1.48

import numpy as np
from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
 
data = load_diabetes()
X = data.data
 
for j, name in enumerate(data.feature_names):
    others = np.delete(X, j, axis=1)
    r2 = LinearRegression().fit(others, X[:, j]).score(others, X[:, j])
    vif = 1 / (1 - r2)
    flag = "  <-- high" if vif > 10 else ""
    print(f"{name:>5}  VIF = {vif:7.2f}{flag}")
 
#   age  VIF =    1.22
#   sex  VIF =    1.28
#   bmi  VIF =    1.51
#    bp  VIF =    1.46
#    s1  VIF =   59.20  <-- high
#    s2  VIF =   39.19  <-- high
#    s3  VIF =   15.40  <-- high
#    s4  VIF =    8.89
#    s5  VIF =   10.08  <-- high
#    s6  VIF =    1.48

The serum measurements are heavily entangled, which is exactly why s1s1 and s5s5 carry large coefficients of opposite sign: the model is taking a difference between two nearly identical columns.

Pitfalls

Compare

🧪 Try It Yourself

Exercise 1 – Add the bias column

Exercise 2 – Solve the Normal Equation

Exercise 3 – Residuals are orthogonal to the features

Exercise 4 – Fit ten features and score on held-out data

Exercise 5 – Measure a variance inflation factor

Recap

• In matrix form the model is $\hat{\mathbf{y}} = \mathbf{X}\boldsymbol{\theta}$, and setting $\nabla J = 0$ gives $\boldsymbol{\theta} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}$.
• The residual is orthogonal to every column of $\mathbf{X}$ — hence residuals summing to zero.
• On an orthogonal design the whole fit is one division: $\boldsymbol{\theta} = \frac{1}{4}[17, 9, 5]^\top$.
• Ten diabetes features reach test $R^2 = 0.48$ against 0.34 for BMI alone.
• Standardise before comparing coefficients; check VIF before interpreting any of them.
• Collinearity destabilises coefficients, not predictions.

Exercise 6 – When the Normal Equation has no inverse

Next

Continue to Polynomial Regression — bend the straight line into a curve without leaving the linear model behind.