42 - Multivariate and polynomial regression, and feature scaling

The price of a house does not depend only on square metres: the number of rooms, the age of the building, the floor also matter. Imagine the price as the output of a mixing console with many sliders: each feature is a slider, and each slider has its own "sensitivity", that is how much it moves the price when you raise it by one notch. The multivariate linear model is the sum of all the sliders, each multiplied by its own sensitivity. In this lesson you extend linear regression to many features, discover what the coefficients mean (and when they do not mean what you think), see why feature scaling is not a detail, and then use an elegant trick, polynomial regression, to make a model that stays linear learn curves.

From the number to the vector: the multivariate model

With n features x₁, ..., x_n the hypothesis becomes:

h(x) = θ₀ + θ₁x₁ + θ₂x₂ + ... + θ_n x_n = θ · x      (with x₀ = 1)

By adding to each example's row a dummy value x₀ = 1, the intercept becomes a coefficient like the others and the whole prediction is a dot product (sum of term-by-term products) between the parameter vector θ and the feature vector. For all m examples together you write h = Xθ, where X has shape (m, n+1) (m rows, one per example; n+1 columns: the one and the n features) and θ has shape (n+1,). The result h has shape (m,): one prediction per example. Checking shapes at every step is the habit that avoids half of the errors.

The nice thing is that nothing changes in the math seen in the previous lessons. The cost function is the same, J = (1/2m)‖Xθ − y‖², the gradient is the same, ∇J = (1/m)Xᵀ(Xθ − y), and the normal equation is the same, θ = (XᵀX)⁻¹Xᵀy. Only the dimension of θ changes.

A concrete case with synthetic data

We build an artificial housing market where we know the truth: the price is 20,000 + 1,500 per square metre + 8,000 per room − 600 per year of age, plus noise with a standard deviation of 15,000 euros. Rooms grow with square metres (bigger houses have more rooms): this introduces a realistic correlation between features. The numbers are invented for the exercise, they do not come from any real market.

import numpy as np

rng = np.random.default_rng(0)
m = 300
area = rng.uniform(30, 200, m)
rooms = np.round(area / 35 + rng.normal(0, 0.6, m)).clip(1, 8)        # correlated with the area
age = rng.uniform(0, 60, m)
price = 20000 + 1500 * area + 8000 * rooms - 600 * age + rng.normal(0, 15000, m)

X = np.column_stack([np.ones(m), area, rooms, age])                   # shape (300, 4)
theta = np.linalg.solve(X.T @ X, X.T @ price)

print("shape of X:", X.shape, " shape of theta:", theta.shape)
print("estimated coefficients [intercept, area, rooms, age]:", np.round(theta, 1))
print("true coefficients                                   : [20000, 1500, 8000, -600]")
print("correlation between area and rooms: %.3f" % np.corrcoef(area, rooms)[0, 1])

shape of X: (300, 4)  shape of theta: (4,)
estimated coefficients [intercept, area, rooms, age]: [19831.5  1489.2  7907.1  -569.6]
true coefficients                                   : [20000, 1500, 8000, -600]
correlation between area and rooms: 0.913

What the code does. The first lines generate the features and the target by applying the true formula plus Gaussian noise. np.column_stack places the column of ones and the three features side by side, giving the (300, 4) matrix. The solution is the same normal equation as in the previous lesson. The estimated coefficients are close to the true ones but not equal, and the largest relative gap is on age (−570 against −600, 5%).

How to interpret a coefficient: θ_area ≈ 1,489 means "holding rooms and age fixed, one more square metre is associated with about 1,489 more euros". The clause "holding all the others fixed" is essential: each coefficient is the effect of that feature with the others kept still. If the features are correlated (here area and rooms are, at 91%), that clause becomes fragile, and we will see it shortly.

Feature scaling

Features have different units and scales: area ranges from 30 to 200, rooms from 1 to 8, age from 0 to 60. As you saw in the lesson on gradient descent, different scales make the cost surface elongated and the descent slow or unstable. The two most common transformations:

Standardization (z-score):   x' = (x − μ) / σ           mean 0, standard deviation 1
Min-max:                     x' = (x − min) / (max − min)   values in [0, 1]

Standardization is the default choice: it is not distorted by a few extreme values as much as min-max is (a single outlier changes max and compresses all the other values). Remember the hygiene rule: μ and σ (or min and max) are computed on the training set and applied identically to the test set; computing them on the whole dataset is data leakage.

mu, sd = X[:, 1:].mean(axis=0), X[:, 1:].std(axis=0)
Z = np.column_stack([np.ones(m), (X[:, 1:] - mu) / sd])              # standardized features (the column of ones stays)
theta_z = np.linalg.solve(Z.T @ Z, Z.T @ price)

print("condition number of X^T X (raw):           %.0f" % np.linalg.cond(X.T @ X))
print("condition number of Z^T Z (standardized):  %.1f" % np.linalg.cond(Z.T @ Z))
print("standardized coefficients [intercept, area, rooms, age]:", np.round(theta_z, 0))
print("in original units (theta_z[j]/sd[j]):", np.round(theta_z[1:] / sd, 1), " intercept:", np.round(theta_z[0] - (theta_z[1:] * mu / sd).sum(), 1))

condition number of X^T X (raw):           179421
condition number of Z^T Z (standardized):  22.0
standardized coefficients [intercept, area, rooms, age]: [211544.  75330.  12258.  -9719.]
in original units (theta_z[j]/sd[j]): [1489.2 7907.1 -569.6]  intercept: 19831.5

What the code does. Standardization is applied column by column (axis=0) to the features only, keeping the column of ones. The condition number of a matrix measures how "unbalanced" it is (ratio between its largest and smallest eigenvalue): a huge value means that some directions are much steeper than others, so the problem is numerically delicate and gradient descent slow. With raw features it is about 179,000; after standardization it drops to 22.

The coefficients obtained on standardized features have an advantage: they are comparable. The "area" coefficient (about 75,000) indicates how much the price changes per standard deviation of area; the age one (about −9,700) per standard deviation of age. Without standardizing, 1,489 and −570 cannot be compared, because square metres and years are different units. Converting back (last line) the coefficients in original units are recovered: the transformation does not change the model, only its representation.

The multicollinearity trap

Two features are collinear if one is (almost) a linear combination of the other. The model cannot tell their contributions apart: it can attribute the effect to one or the other, in different proportions depending on small changes in the data. The prediction stays good, but the coefficients become unstable and uninterpretable. To measure it we resample the data with the bootstrap (drawing m rows with replacement) 300 times, redo the regression each time and look at the standard deviation of each coefficient. Then we add a nearly duplicated feature, "area in another unit" (the area plus a little noise), and repeat:

def bootstrap_coefficients(X, y, B=300, seed=1):
    r = np.random.default_rng(seed)
    out = []
    for _ in range(B):
        i = r.integers(0, len(y), len(y))
        out.append(np.linalg.lstsq(X[i], y[i], rcond=None)[0])
    return np.array(out)

b_normal = bootstrap_coefficients(X, price)
X_dup = np.column_stack([X, area + np.random.default_rng(2).normal(0, 0.5, m)])   # nearly duplicated area
b_dup = bootstrap_coefficients(X_dup, price)

print("normal model    mean   :", np.round(b_normal.mean(axis=0), 0))
print("normal model    std dev:", np.round(b_normal.std(axis=0), 0))
print("duplicated area mean   :", np.round(b_dup.mean(axis=0), 0))
print("duplicated area std dev:", np.round(b_dup.std(axis=0), 0))
print("sum of the two area coefficients: %.0f (in the normal model: %.0f)" % (b_dup.mean(axis=0)[1] + b_dup.mean(axis=0)[4], b_normal.mean(axis=0)[1]))

normal model    mean   : [19654.  1491.  7855.  -567.]
normal model    std dev: [2709.   49. 1542.   60.]
duplicated area mean   : [19846.   758.  7884.  -567.   730.]
duplicated area std dev: [2755. 1835. 1532.   60. 1837.]
sum of the two area coefficients: 1489 (in the normal model: 1491)

What the code does. bootstrap_coefficients draws indices with replacement, redoes the regression with np.linalg.lstsq (NumPy's least-squares solver, documentation, which also accepts matrices that are not of full rank, where the normal equation breaks down) and collects the coefficients. The comparison is eloquent: with the duplicated feature, the standard deviation of the area coefficient goes from 49 to 1,835, about thirty-seven times more uncertain, and the same holds for the copy (1,837).

The means share the work arbitrarily (758 and 730), but their sum (1,489) stays close to the value of the normal model (1,491): the model knows how much "area" weighs, not to which of the two to attribute it. Even rooms, which are correlated with area at 91%, have an uncertain coefficient (standard deviation 1,542 on a value of about 7,900).

Remedies: remove or combine redundant features, use regularization (the topic of the lesson on overfitting and regularization), or reduce dimensionality with PCA. And an interpretation rule: do not read the coefficients of strongly correlated features one by one.

Polynomial regression: curves with a linear model

A curved relation, like a price that grows more than proportionally with area, does not fit a line well. The trick is surprising in its simplicity: if the curvature is missing, create it as a new feature. From a single feature x you build x, x², x³, ..., x^d and give them to the same linear regression:

h(x) = θ₀ + θ₁x + θ₂x² + θ₃x³ + ... + θ_d x^d

The model stays linear in the parameters θ: the function is a curve in x, but for the learning algorithm the "data" are simply the columns [x, x², x³, ...]. Everything you have seen (normal equation, gradient, descent) works without changes. The degree d is the hyperparameter that controls flexibility: d = 1 is a line, a high d a winding curve.

The picture shows the same cloud of points (a noisy sine wave) fitted with three degrees. Degree 1 is underfitted: too rigid, it cannot represent the trend. Degree 4 follows the curve. Degree 15 is overfitted: it passes close to every training point chasing the noise, and makes wild oscillations where data are sparse. Let us see it with numbers, using a train/test split and implementing the power matrix by hand:

from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.linear_model import LinearRegression

r = np.random.default_rng(1)
x = r.uniform(0, 6, 60)
y = np.sin(x) + r.normal(0, 0.25, 60)                                  # curved relation + noise (sigma 0.25, variance 0.0625)
x_tr, x_te, y_tr, y_te = train_test_split(x, y, test_size=0.5, random_state=1)

def power_matrix(x, degree):
    return np.column_stack([x ** k for k in range(1, degree + 1)])      # columns x, x^2, ..., x^degree

def fit_by_hand(x, y, degree):
    P = power_matrix(x, degree)
    mu, sd = P.mean(axis=0), P.std(axis=0)
    A = np.column_stack([np.ones(len(x)), (P - mu) / sd])               # standardize the powers
    theta = np.linalg.lstsq(A, y, rcond=None)[0]
    predict = lambda xn: np.column_stack([np.ones(len(xn)), (power_matrix(xn, degree) - mu) / sd]) @ theta
    return theta, predict

# check: by hand and with the scikit-learn pipeline give the same predictions (degree 3)
_, pred_hand = fit_by_hand(x_tr, y_tr, 3)
pipe3 = make_pipeline(PolynomialFeatures(3, include_bias=False), StandardScaler(), LinearRegression()).fit(x_tr[:, None], y_tr)
print("max difference between by-hand and sklearn (degree 3): %.2e" % np.abs(pred_hand(x_te) - pipe3.predict(x_te[:, None])).max())

print("degree   train MSE   test MSE")
for degree in (1, 2, 3, 4, 5, 7, 10, 15):
    _, pred = fit_by_hand(x_tr, y_tr, degree)
    print("%6d   %9.4f   %8.4f" % (degree, ((pred(x_tr) - y_tr) ** 2).mean(), ((pred(x_te) - y_te) ** 2).mean()))

max difference between by-hand and sklearn (degree 3): 3.18e-14
degree   train MSE   test MSE
     1      0.1757     0.2357
     2      0.1670     0.2067
     3      0.0318     0.0358
     4      0.0318     0.0359
     5      0.0318     0.0358
     7      0.0271     0.0340
    10      0.0232     0.1026
    15      0.0196     3.4167

What the code does. power_matrix builds the columns x¹, ..., x^d. fit_by_hand standardizes them (high powers have enormous scales: 6¹⁵ is of the order of 5·10¹¹, and without standardization the matrix is practically singular), solves by least squares and returns a prediction function that reapplies the same standardization to new points. The comparison with the scikit-learn pipeline (PolynomialFeatures makes the powers, StandardScaler scales them) confirms the same predictions up to rounding.

Read the table carefully. The error on the training set falls monotonically as the degree grows (from 0.176 to 0.020): more flexibility, the more the model fits the data it has seen. The error on the test set instead falls down to degree 3 (from 0.236 to 0.036), stays flat up to degree 7 (0.034), and then explodes: 0.103 at degree 10 and 3.42 at degree 15, over a hundred times the training error. A comparison based only on the training error would always pick the highest degree; the test says the opposite, and the "U" shape of the test error is the signature of overfitting. A caution: the test set has 30 points, so small differences (0.034 versus 0.036) are not significant. Even the fact that the test error is below the noise variance (0.0625) is an effect of the small sample, not a violation of the theoretical limit.

Extrapolating with polynomials is dangerous

model9 = make_pipeline(PolynomialFeatures(9, include_bias=False), StandardScaler(), LinearRegression()).fit(x_tr[:, None], y_tr)
for xn in (6.5, 7.0, 8.0):
    print("x = %.1f   degree-9 prediction: %10.2f   true value sin(x) = %6.3f" % (xn, model9.predict([[xn]])[0], np.sin(xn)))
print("largest coefficient (absolute value) of degree 9:", np.round(np.abs(model9[-1].coef_).max(), 0))

x = 6.5   degree-9 prediction:     -13.61   true value sin(x) =  0.215
x = 7.0   degree-9 prediction:     -67.92   true value sin(x) =  0.657
x = 8.0   degree-9 prediction:    -663.09   true value sin(x) =  0.989
largest coefficient (absolute value) of degree 9: 10946.0

What the code does: it trains a degree-9 polynomial and queries it just outside the data range (0-6). A high-degree polynomial has a dominant x⁹ term that explodes outside the training range: already at x = 6.5 the prediction is −13.6 against a true value of 0.215; at x = 7 it is −67.9 and at x = 8 it is −663, while the sine never leaves [−1, 1]. The coefficients are also enormous (the largest exceeds 10,000), another symptom of instability. The degree-1 line extrapolates in a boring but bounded way; the high-degree polynomial does not.

Trade-offs

• More features: more information, but more risk of collinearity, overfitting and computational cost (the matrix XᵀX is (n+1)×(n+1)).
• Polynomial degree: more degree = more flexibility = more variance. It is chosen with cross-validation, not with the training error.
• Normal equation or gradient descent: the first is exact but costs roughly cubic time in the number of features (inverting an n×n matrix) and does not work for non-linear models; the second scales to many features and huge datasets, but needs standardized features and a carefully chosen learning rate.

The three items are the same compromise between the capacity of the model and the reliability of the estimate, seen from three sides. More features and more degrees increase flexibility and, with it, variance; the normal equation or gradient descent do not change what the model can learn, but how much it costs to get there. The practical rule that follows is to proceed by measured additions: one more feature or one more degree at a time, each judged by validation and not by the training error.

Real-world traps

• Interpreting a coefficient without the clause "holding the other features fixed", especially with correlated features.
• Scaler fitted on the whole dataset (data leakage). It goes in a Pipeline or is computed on the training set only.
• Forgetting to apply the same scaler to the test set and in production: a model trained on standardized data receives raw data and produces absurd numbers with no error at all.
• Choosing the degree by looking at the training set (or at the test set, used several times).
• Extrapolating with high-degree polynomials.
• Categorical variables treated as numbers (for example the neighbourhood coded 1, 2, 3): they invent an order that does not exist; one-hot encoding is used (a 0/1 column per category; in scikit-learn it is done by OneHotEncoder, which with the drop parameter removes the reference category).
• Shape mismatch: PolynomialFeatures and scikit-learn want X of shape (m, n); an (m,) vector must be converted with x[:, None] or reshape(-1, 1).

These seven traps have a common root: the model never complains. A misread coefficient, a forgotten scaler, a degree chosen on the training set produce plausible numbers and no error. This is why the experiments in this lesson look at the coefficients and at validation, not only at the R², which in the zones example differed by less than half a point between the right coding and the wrong one.

Categorical variables: why numeric coding misleads and how one-hot works

A linear model multiplies each feature by a coefficient: it understands only quantities. A variable like the "zone" of a house (centre, semi-centre, outskirts) is a category, not a quantity. The temptation is to code it 0, 1, 2. But this imposes a structure the data do not have: it says that going from zone 0 to zone 1 is worth exactly as much as going from zone 1 to zone 2, and that zone 2 is "twice" zone 1. The solution is one-hot encoding: a 0/1 column for each category ("is the house in this zone?"). The model then learns a separate contribution for each. There is a subtlety: if you keep all the columns and also the intercept, the summed columns are always 1, exactly the intercept column, and this exact dependency is perfect multicollinearity (the trap seen before): for this reason one is removed (the "reference category"). We run the experiment with a price generated as 50 + 2·square metres + zone effect (0, 30 and 80) + noise of deviation 15.

from sklearn.model_selection import KFold, cross_val_score
from sklearn.linear_model import Ridge

rng_p = np.random.default_rng(0)
n_p = 300
zone = rng_p.integers(0, 3, n_p)
size = rng_p.uniform(40, 140, n_p)
zone_effect = np.array([0, 30, 80])[zone]
price = 50 + 2.0 * size + zone_effect + rng_p.normal(0, 15, n_p)

coded = LinearRegression().fit(np.column_stack([size, zone]), price)
one_hot = np.eye(3)[zone]
dummies = LinearRegression().fit(np.column_stack([size, one_hot[:, 1:]]), price)

print("numeric coding 0/1/2:  square-metre coefficient %.3f   effect per zone «step» %.2f   R² %.4f" % (coded.coef_[0], coded.coef_[1], coded.score(np.column_stack([size, zone]), price)))
print("one-hot (zone 0 as reference): square metres %.3f   zone 1 %.2f   zone 2 %.2f   intercept %.2f   R² %.4f" % (dummies.coef_[0], dummies.coef_[1], dummies.coef_[2], dummies.intercept_, dummies.score(np.column_stack([size, one_hot[:, 1:]]), price)))
print("true values: square metres 2.000   zone 1 30   zone 2 80   intercept 50")

full = LinearRegression().fit(np.column_stack([size, one_hot]), price)
rank = np.linalg.matrix_rank(np.column_stack([np.ones(n_p), size, one_hot]))
print("\nall three columns + intercept: rank %d out of 5 columns; coefficients of the three zones %s (none has a meaning on its own)" % (rank, np.round(full.coef_[1:], 2)))

numeric coding 0/1/2:  square-metre coefficient 1.987   effect per zone «step» 40.06   R² 0.9483
one-hot (zone 0 as reference): square metres 1.989   zone 1 30.98   zone 2 79.44   intercept 50.58   R² 0.9520
true values: square metres 2.000   zone 1 30   zone 2 80   intercept 50

all three columns + intercept: rank 4 out of 5 columns; coefficients of the three zones [-36.8   -5.83  42.63] (none has a meaning on its own)

What the code does. It generates 300 houses with three zones and fits two models. With the numeric coding 0/1/2 the square-metre coefficient is 1.987 (true: 2) but the zone effect is estimated as a single value, 40.06 per "step": this means zone 1 is worth 40.06 more than zone 0 (true: 30) and zone 2 is worth 80.12 more (true: 80). The model has to fit two different effects onto a line, and it gets the intermediate zone wrong by about 10 units. The one-hot, instead, recovers 30.98 and 79.44 (true: 30 and 80) with an intercept of 50.58 (true: 50).

Look at the R²: 0.9483 versus 0.9520, a difference of only 0.4 points. The R² would never have flagged the error: the numeric coding "works" in the statistical sense, but attributes to the intermediate zone a price wrong by almost a third. It is an example of an error that shows in the coefficients and not in the aggregate metrics.

The last line shows the dummy-variable trap: with all three columns plus the intercept the rank is 4 out of 5 columns, the system is singular, and the returned coefficients (−36.8, −5.83, 42.63) are one of the infinitely many equivalent solutions: none has a meaning on its own, but their differences do (−5.83 − (−36.8) = 30.97, the same estimate of zone 1 relative to zone 0).

Interactions: when the effect of one feature depends on another

An additive linear model assumes that the effect of each feature is the same whatever the value of the others. Sometimes it is not so: the effect of square metres on price may be different in the centre and on the outskirts. It is captured by adding as a feature the product of the two (interaction term): the model remains linear in the coefficients, but the slope with respect to x now depends on the group. We generate data with two groups having a different slope (+0.5 for group 0, and 0.5 − 1.5 = −1 for group 1) and compare the two models with cross-validation.

x_i = rng_p.uniform(0, 10, n_p)
group = rng_p.integers(0, 2, n_p)
y_i = 1 + 0.5 * x_i + group * (2 - 1.5 * x_i) + rng_p.normal(0, 1, n_p)    # group 1 has intercept +2 and slope −1.5 extra

folds = KFold(5, shuffle=True, random_state=0)
additive = np.column_stack([x_i, group])
with_interaction = np.column_stack([x_i, group, x_i * group])
print("cross-validated R², additive model:                   %.3f" % cross_val_score(LinearRegression(), additive, y_i, cv=folds, scoring="r2").mean())
print("cross-validated R², with the x·group term:            %.3f" % cross_val_score(LinearRegression(), with_interaction, y_i, cv=folds, scoring="r2").mean())
fit_i = LinearRegression().fit(with_interaction, y_i)
print("coefficients with interaction: intercept %.3f   x %.3f   group %.3f   x·group %.3f   (true: 1, 0.5, 2, −1.5)" % (fit_i.intercept_, *fit_i.coef_))

cross-validated R², additive model:                   0.585
cross-validated R², with the x·group term:            0.913
coefficients with interaction: intercept 1.146   x 0.486   group 1.693   x·group -1.450   (true: 1, 0.5, 2, −1.5)

What the code does. It compares, with 5-fold cross-validation, an additive model (x and group) with one that adds the product x·group. The R² in validation goes from 0.585 to 0.913: the additive model cannot represent two lines with different slopes and fits an average one, while the interaction allows it. The coefficients of the interaction model (intercept 1.146; x 0.486; group 1.693; x·group −1.450) are close to the true ones (1; 0.5; 2; −1.5).

To read them: the slope of x in group 0 is 0.486; in group 1 it is 0.486 − 1.450 = −0.964; and this is the meaning of the interaction term, it is the difference between the slopes of the two groups.

Beware of two things. First, with an interaction the coefficient of "x" alone is no longer "the effect of x", but the effect of x when the group is 0: if you do not keep this in mind you will misread the model. Second, interactions are many (among p features there are p(p − 1)/2 pairs), and each one added without a reason is one more degree of freedom for overfitting. Add them guided by domain knowledge or by an exploration of the residuals, and always verify the gain in validation, as here.

Choosing the polynomial degree: training, validation and truth

x_p = rng_p.uniform(-3, 3, 60)
y_p = np.sin(x_p) + 0.3 * rng_p.normal(size=60)
x_grid = np.linspace(-3, 3, 400)                                            # dense grid to measure the error against the true curve sin(x)

print("degree   training MSE   cross-validated MSE   MSE against the true curve")
for degree in (1, 2, 3, 5, 7, 9, 12, 15):
    pipe = make_pipeline(PolynomialFeatures(degree), StandardScaler(), LinearRegression())
    cv_mse = -cross_val_score(pipe, x_p[:, None], y_p, cv=folds, scoring="neg_mean_squared_error").mean()
    pipe.fit(x_p[:, None], y_p)
    print("%-6d   %-14.3f   %-28.3f   %.3f" % (degree, np.mean((pipe.predict(x_p[:, None]) - y_p) ** 2), cv_mse, np.mean((pipe.predict(x_grid[:, None]) - np.sin(x_grid)) ** 2)))

degree   training MSE   cross-validated MSE   MSE against the true curve
1        0.225            0.242                          0.176
2        0.215            0.261                          0.182
3        0.082            0.093                          0.014
5        0.080            0.104                          0.013
7        0.076            0.103                          0.018
9        0.071            0.135                          0.028
12       0.064            0.499                          0.045
15       0.058            2.559                          0.069

What the code does. It generates 60 points from the sin(x) curve with noise 0.3 (variance 0.09), fits polynomials of increasing degree and reports three errors: on training, in cross-validation (5 folds) and against the true curve on a dense grid, a quantity that in practice you do not have but that here we can compute because we know the truth. The training error always falls, from 0.225 with degree 1 to 0.058 with degree 15: more flexibility, more fit to the points, by construction. The validation error has a minimum at degree 3 (0.093), rises slowly to 0.104 with degree 5 and 0.103 with degree 7, then explodes: 0.135 with degree 9, 0.499 with degree 12, 2.559 with degree 15. The gap between training and validation is overfitting opening up.

The error against the true curve gives the same indication (minimum 0.013–0.014 for degrees 3 and 5), but with an important difference: at degree 15 it is only 0.069, not 2.559. Cross-validation exaggerates the damage, and there are two plausible explanations, which I have not isolated here with an experiment: each validation fold contains points at the edges of the range that the polynomial, trained on the other folds, must extrapolate (where high-degree polynomials are terrible, as seen), and the MSE estimate is dominated by a few points with enormous errors. It would not be a flaw of the method, but an indication that with little data and very flexible models the CV estimate has high variance.

The operating conclusion remains valid: the right degree is 3, chosen without looking at the true curve, and degrees between 3 and 7 are nearly equivalent; choose the simplest within the flat zone.

Measuring multicollinearity: the VIF and its cost on the coefficients

The variance inflation factor (VIF) of a feature is 1 / (1 − R²_j), where R²_j is the R² of the regression of that feature on all the others. If the feature is predictable from the others, R²_j is close to 1 and the VIF explodes: the coefficient has a variance equal to VIF times what it would have if the feature were independent. We create three features: x₁, x₂ ≈ x₁ with a small noise (deviation 0.1), and an independent x₃, with target y = x₁ + x₃ + noise.

x1 = rng_p.normal(size=500)
x2 = x1 + rng_p.normal(0, 0.1, 500)                                         # almost a copy of x1
x3 = rng_p.normal(size=500)
X_v = np.column_stack([x1, x2, x3])
y_v = x1 + x3 + rng_p.normal(0, 1, 500)

for j in range(3):
    others = np.delete(X_v, j, axis=1)
    r2_j = LinearRegression().fit(others, X_v[:, j]).score(others, X_v[:, j])
    print("VIF of feature x%d: %.1f" % (j + 1, 1 / (1 - r2_j)))

boot = np.array([LinearRegression().fit(X_v[i], y_v[i]).coef_ for i in (np.random.default_rng(s).integers(0, 500, 500) for s in range(200))])
print("coefficients (mean over 200 resamples):   %s   (truth: 1, 0, 1)" % np.round(boot.mean(axis=0), 3))
print("standard deviation of the coefficients:           %s" % np.round(boot.std(axis=0), 3))
print("standard deviation of the sum of the coefficients of x1 and x2: %.3f" % (boot[:, 0] + boot[:, 1]).std())

VIF of feature x1: 112.8
VIF of feature x2: 112.8
VIF of feature x3: 1.0
coefficients (mean over 200 resamples):   [ 1.048 -0.072  1.032]   (truth: 1, 0, 1)
standard deviation of the coefficients:           [0.481 0.477 0.043]
standard deviation of the sum of the coefficients of x1 and x2: 0.046

What the code does. For each feature it computes the R² of the regression on the other two and derives the VIF. x₁ and x₂ each have a VIF of 112.8 (their correlation is practically total), x₃ has 1.0, the minimum possible: it is independent of the others. As a rule of thumb a VIF above 5–10 is considered worrying; here we are more than twenty times the strictest threshold of 5.

The cost can be read in the last part, which resamples the data 200 times and refits the model. The mean coefficients are 1.048, −0.072 and 1.032 (truth: 1, 0, 1): on average the model is not wrong, and this is why multicollinearity is insidious, because the predictions are good.

But the standard deviation of the coefficients of x₁ and x₂ is 0.481 and 0.477, while for x₃ it is 0.043: eleven times larger, and √112.8 ≈ 10.6, exactly as the theory predicts (the standard deviation grows with the square root of the VIF). The line added at the end confirms it: the standard deviation of the sum of the two coefficients is 0.046, ten times smaller than that of each one (about 0.48). The model knows how much y depends on the pair, but cannot distribute it between the two, and in a single sample the two coefficients compensate each other. If you only care about predicting, you can ignore the problem; if you interpret the coefficients ("how much does x₁ weigh?"), the conclusions are unreliable, and it pays to remove one of the two features, combine them, or use regularization, as in overfitting and regularization.

Three exercises

Exercise 1. In the zone experiment, change the reference category (use zone 2 instead of zone 0): how do the intercept and coefficients change, and what remains unchanged? Why does an arbitrary choice of reference not change the predictions? Exercise 2. Replace linear regression with Ridge (alpha = 1) in the VIF experiment: does the standard deviation of the coefficients of x₁ and x₂ change? What compromise do you pay? Exercise 3. Repeat the degree selection with 600 points instead of 60: does the best degree in cross-validation remain 3? What happens to the gap between training and validation error with high degrees?

Other lessons in the series that link to this one: Linear Algebra 1: Vectors, Matrices, Systems.

Sources

• scikit-learn, Linear Models - the description of multicollinearity (least-squares estimate sensitive to error on the target) and the section on polynomial regression as a linear model on non-linear functions of the data, consulted on 2026-11-12
• scikit-learn, Preprocessing data - standardization (removing the mean and dividing by the standard deviation), consulted on 2026-11-12
• scikit-learn, OneHotEncoder - one-hot encoding and the drop parameter, consulted on 2026-11-12
• NumPy, numpy.linalg.lstsq - the least-squares solver used in the bootstrap, consulted on 2026-11-12

The figures throughout the lesson are the author's measurements: Python 3.14.4, NumPy 2.3.5, scikit-learn 1.8.0, every listing run twice with a fixed seed and identical output. The VIF thresholds (5-10) are a convention, not a standard.

In the next lesson the same machinery is used to predict a category instead of a number: logistic regression, from the line to the probability.