Linear Algebra for Machine Learning

If you have ever looked at a spreadsheet, you have already looked at the central object of machine learning. A table with one row per customer and one column per thing you know about them is, in mathematical language, a matrix. Nothing more exotic than that.

The examples begin with multiplying and adding numbers, then introduce the geometry and notation as they become useful. The code uses NumPy arrays; run the body blocks in order so later calculations can reuse earlier variables. The projection and conditioning sections take more work than the opening arithmetic, and can be revisited after the examples.

A data table is a matrix

Suppose we are looking at five houses, and for each one we know three things: floor area in square meters, the number of rooms, and the age of the building in years. Written as a table that is five rows by three columns.

Which way around this goes is a convention rather than a law, and it is worth stating explicitly because the two orientations are easy to confuse. In scikit-learn and throughout this series, a row is one observation and a column is one feature. One row is one house; one column is one measurement taken on every house. Other conventions exist — some statistics texts and some signal-processing code put features in rows — so the first thing to check when borrowing code is which layout it expects.

import numpy as np

# five houses, three columns: floor area (m2), number of rooms, age (years)
X = np.array([
    [ 85.0, 3.0, 12.0],
    [120.0, 4.0,  5.0],
    [ 60.0, 2.0, 30.0],
    [150.0, 5.0,  2.0],
    [ 95.0, 3.0, 18.0],
])
print("shape (rows, columns):", X.shape)
print("row 0 (first house)  :", X[0])
print("column 0 (all areas) :", X[:, 0])
# shape (rows, columns): (5, 3)
# row 0 (first house)  : [85.  3. 12.]
# column 0 (all areas) : [ 85. 120.  60. 150.  95.]

That table has a standard name and standard letters. It is called the design matrix and written \(X\). The number of rows is \(n\) — here five houses — and the number of columns is \(p\) — here three features. So \(X\) has shape \(n \times p\), which people write as \(X \in \mathbb{R}^{n \times p}\). Read that symbol as “\(X\) is a table of real numbers with \(n\) rows and \(p\) columns”; it says nothing more than that.

Many tabular problems have \(n\) much larger than \(p\) — a million customers, fifty things known about each — and that tall, thin shape explains a good deal of what follows. It is not universal. Genomics, text and many scientific datasets routinely have more features than rows, and \(p > n\) changes which of the results below still apply; the covariance section returns to this.

Making a prediction is multiplying and adding

Now suppose someone hands us a rule for pricing a house: charge 3 units per square meter, add 10 units for each room, and subtract 0.5 units for each year of age. To price the first house we multiply each of its numbers by the matching rate and add the results up.

That is \(85 \times 3 + 3 \times 10 + 12 \times (-0.5) = 255 + 30-6 = 279\). Multiply pairwise, then sum. This operation happens so often that it has its own name, the dot product, and its own notation: mathematicians write it \(a^{\top}b\). In NumPy it is computed with @, which is Python’s matrix-multiplication operator rather than a symbol for the dot product itself.

# X is the house table from the previous block
price_rule = np.array([3.0, 10.0, -0.5])   # per m2, per room, per year of age

longhand = 85.0 * 3.0 + 3.0 * 10.0 + 12.0 * (-0.5)
print("written out by hand :", longhand)
print("as a dot product    :", X[0] @ price_rule)
print("all five at once    :", X @ price_rule)
# written out by hand : 279.0
# as a dot product    : 279.0
# all five at once    : [279.  397.5 185.  499.  306. ]

The last line is the point. Writing X @ price_rule priced all five houses in one stroke, because multiplying a whole table by a list of rates just repeats the dot product once per row. In symbols the list of rates is called \(\beta\) and the whole operation is written \(X\beta\).

The operation \(X\beta\) is the core computation inside a linear model. Linear regression uses it directly. Logistic regression passes the result through the logistic sigmoid, the inverse of the logit link. A neural-network layer normally adds a bias term and then a nonlinear activation. For simplicity this pricing example has no intercept — a baseline price applying to every house. Most practical models compute \(X\beta + b\), or equivalently glue a column of ones onto \(X\) so the intercept rides along inside \(\beta\), which is what the temperature example later does. The numbers in \(\beta\) are called weights or coefficients; they are the “rates” in our pricing rule.

What the dot product actually measures

So far the dot product was just bookkeeping. It also has a geometric meaning, and that meaning is why it appears everywhere.

Think of a list of two numbers as an arrow drawn from the origin on a sheet of paper: \([3, 4]\) means go 3 right and 4 up. The dot product of two such arrows combines two things — how long they are, and how closely they point the same way — through the relation \(a^{\top}b = \|a\|\,\|b\|\cos\theta\). With the lengths held fixed it is large and positive when the arrows point the same way, exactly zero when they are at right angles, and negative when they point in opposite directions.

import numpy as np

a = np.array([3.0, 4.0])
print("length of a:", np.linalg.norm(a))

for name, b in [("same direction", np.array([6.0, 8.0])),
                ("perpendicular ", np.array([-4.0, 3.0])),
                ("opposite      ", np.array([-3.0, -4.0]))]:
    dot = a @ b
    cos = dot / (np.linalg.norm(a) * np.linalg.norm(b))
    print(f"{name}  dot {dot:6.1f}   cosine {cos:5.2f}")
# length of a: 5.0
# same direction  dot   50.0   cosine  1.00
# perpendicular   dot    0.0   cosine  0.00
# opposite        dot  -25.0   cosine -1.00

Two further pieces of vocabulary come straight out of that snippet. The length of an arrow — its norm, written \(\|a\|\) — is the ordinary distance from the origin, found by Pythagoras: \(\sqrt{3^2 + 4^2} = 5\). For two nonzero vectors, the cosine is the dot product after both arrows have been rescaled to length one, which strips out size and leaves only direction. It lies between \(-1\) and \(1\). A zero vector has no direction, so this cosine formula is undefined when either norm is zero.

Nothing about this changes when the arrows have fifty numbers instead of two. We cannot draw it, but every formula still works, and “perpendicular” still means the dot product is zero. That is the practical reason the geometric language is worth keeping: it lets you reason about fifty-dimensional data using pictures from a sheet of paper.

Cosine or distance? The choice changes the answer

There are two natural ways to ask whether two rows of a table are similar, and they disagree. Euclidean distance — ordinary straight-line distance — asks “are these two close together?” Cosine asks “do these two point the same way?”

Take two customers who buy the same products in the same proportions, but one spends ten times as much. Their cosine is 1, reflecting identical purchase proportions in this representation. Their Euclidean distance reflects the difference in amounts; its numerical size depends on the units. The useful comparison depends on the task.

• If purchase proportions matter more than total spending, cosine is the reasonable choice.
• If the absolute amounts are the quantity of interest, Euclidean distance keeps the information cosine throws away.
• Cosine is often useful for document vectors when relative word usage matters more than overall magnitude. Document length alone does not establish whether two topics differ.

The metric and preprocessing jointly define similarity. Purchase proportions and spending amounts answer different questions; compare choices against the task instead of treating either metric as universally correct.

Fitting: finding the weights from data

Until now we were handed the pricing rule. In practice nobody hands it to us: we see the houses and the prices they sold for, and we must work backwards to the rates. That is what fitting a model means.

Below we generate prices from the rule we invented, add a little random noise to imitate the messiness of real sales, then throw the rule away and try to recover it from the numbers alone.

import numpy as np

rng = np.random.default_rng(0)
true_rule = np.array([3.0, 10.0, -0.5])
y = X @ true_rule + rng.normal(scale=4.0, size=5)   # the prices they sold for

fitted = np.linalg.lstsq(X, y, rcond=None)[0]
print("rule we invented :", true_rule)
print("rule recovered   :", np.round(fitted, 3))
print("actual prices    :", np.round(y, 2))
print("fitted prices    :", np.round(X @ fitted, 2))
print("leftover         :", np.round(y - X @ fitted, 2))
# rule we invented : [ 3.  10.  -0.5]
# rule recovered   : [ 2.701 18.814 -0.442]
# actual prices    : [279.5  396.97 187.56 499.42 303.86]
# fitted prices    : [280.7  397.14 186.41 498.31 305.06]
# leftover         : [-1.2  -0.17  1.15  1.11 -1.2 ]

The fitted prices sit within about one unit of the actual prices, so the model reproduces these five houses closely. But the recovered room rate is 18.8 where the truth was 10. Note what has and has not been shown: we scored the fit on the same five houses it was fitted to, so this says nothing yet about unseen houses. What it does show is that a close fit and trustworthy coefficients are separate things.

Two reasons. Five houses is almost no data for estimating three rates. And area and room count move together — big houses have more rooms — so the fit can trade one against the other without changing its predictions much. That second reason is important enough that the rest of this article keeps returning to it.

Two symbols to meet before the next section

The next section uses two operations that have appeared only as symbols so far. Both are simpler than the notation suggests, and it is worth spending a minute on each before they arrive in a formula.

Transpose, written \(X^{\top}\), flips a table on its diagonal: rows become columns and columns become rows. Nothing is added or removed, only relabeled.

The inverse, written \(A^{-1}\), is the matrix that undoes what \(A\) does. If multiplying by \(A\) turns one list of numbers into another, multiplying that result by \(A^{-1}\) gets you back where you started.

import numpy as np

small = np.array([[ 85.0, 3.0],                # two houses, two features
                  [120.0, 4.0]])
print("rows are houses:")
print(small)
print("transposed — rows are now features:")
print(small.T)

A = np.array([[2.0, 1.0],
              [1.0, 3.0]])
v = np.array([4.0, 5.0])
print("A @ v            :", A @ v)
print("undone by inv(A) :", np.linalg.inv(A) @ (A @ v))
# rows are houses:
# [[ 85.   3.]
#  [120.   4.]]
# transposed — rows are now features:
# [[ 85. 120.]
#  [  3.   4.]]
# A @ v            : [13. 19.]
# undone by inv(A) : [4. 5.]

Two cautions carry into the next section. Only square tables can have an inverse, and even then it may not exist — when one column is a copy or a combination of the others, no operation can undo the flattening that caused. That case has a name, and it appears later under the heading of collinearity.

Least squares is a shadow

What exactly did lstsq optimize? It chose the weights that make the total of the squared leftovers as small as possible. “Leftover” is the gap between the actual price and the fitted price, normally called the residual. Squaring makes every gap nonnegative, so overshooting by 3 counts the same as undershooting by 3, and penalizes big misses much more than small ones.

In symbols: find the \(\beta\) that minimizes \(\|y-X\beta\|_2^2\). Reading it left to right — take the actual values \(y\), subtract the predictions \(X\beta\), measure the length of what remains, and make that length as small as you can.

We can also build the prediction vector by scaling each column of \(X\) and adding the results. That is exactly what multiplying \(X\) by a weight vector does.

combo = 3.0 * X[:, 0] + 10.0 * X[:, 1] + (-0.5) * X[:, 2]
print("weighted columns added:", combo)
print("X @ price_rule        :", X @ price_rule)
# weighted columns added: [279.  397.5 185.  499.  306. ]
# X @ price_rule        : [279.  397.5 185.  499.  306. ]

For this full-column-rank house table, changing the weights changes the prediction vector.

Whatever weights you choose, \(X\beta\) can only produce prices that are some mixture of the columns of \(X\). All the reachable combinations form a linear subspace. A sheet of paper through the origin in a three-dimensional room is a lower-dimensional analogy; here the three independent columns span a three-dimensional subspace of five-dimensional prediction vectors. This surface has a name, the column space of \(X\).

The vector of actual prices \(y\) is a point somewhere in the room, and in general it is not on that sheet of paper — no set of rates reproduces the real prices exactly. So we settle for the closest point on the sheet. The closest point is found by dropping straight down onto it: it is the shadow of \(y\) on the sheet, cast by a light directly above. Mathematicians call the shadow a projection.

That picture immediately gives the defining property. If you drop straight down, the line from \(y\) to its shadow is perpendicular to the sheet — and therefore perpendicular to every column of \(X\). Writing “perpendicular” as “dot product equals zero” gives \(X^{\top}(y-X\beta) = 0\), which rearranges into the famous normal equations \(X^{\top}X\beta = X^{\top}y\). The word “normal” here is the old geometric term for perpendicular; it has nothing to do with the normal distribution.

When no column of \(X\) can be written as a linear combination of the others, \(X\) has full column rank. In that case, the matrix that casts the shadow is \(H = X(X^{\top}X)^{-1}X^{\top}\). It is nicknamed the hat matrix because it turns \(y\) into \(\hat{y}\), the fitted values. The code below checks the projection using a separate random matrix.

If the columns are linearly dependent, that inverse does not exist. The more general expression is \(H = XX^{+}\), where \(X^{+}\) is the pseudoinverse. NumPy computes it with np.linalg.pinv using the SVD, which we introduce later. Where the inverse exists, the two expressions give the same least-squares projection.

import numpy as np

rng = np.random.default_rng(0)
Xr = rng.normal(size=(50, 3))                 # fresh random data, not the houses
yr = rng.normal(size=50)

br = np.linalg.lstsq(Xr, yr, rcond=None)[0]
yr_hat = Xr @ br
Hr = Xr @ np.linalg.inv(Xr.T @ Xr) @ Xr.T     # the hat matrix

print("lstsq equals projection:", np.allclose(yr_hat, Hr @ yr))
print("residual vs columns:", np.abs(Xr.T @ (yr - yr_hat)).max())
# lstsq equals projection: True
# residual vs columns: 8.722189637211386e-15

The second number is about 8.7e-15, a tiny floating-point residual rather than an exact zero. This numerical check agrees with the orthogonality predicted by the least-squares equations at the scale of this example.

One practical footnote. The line computing the hat matrix — \(H\) in the formula, Hr in the code — explicitly inverts \(X^{\top}X\), which is written that way here only to mirror the formula. Numerical libraries generally avoid forming \((X^{\top}X)^{-1}\) explicitly, for reasons covered at the end of the article. NumPy’s lstsq uses an SVD instead, which is more numerically stable; whether it is also faster depends on the shape of the problem and the solver.

When there is no single right answer: collinearity

Here is a situation that stretches the shadow picture in an instructive way. The shadow itself survives — the projection \(\hat{y}\) stays unique — but the coordinates used to describe it stop being unique. Imagine your data table records the temperature twice — once in Celsius and once in Fahrenheit. The second column carries no information the first did not already have; it is just the first column rescaled and shifted.

In a model that also includes an intercept, the separate coefficients of Celsius and Fahrenheit cannot be identified: the intercept and two temperature columns are linearly dependent. Infinitely many coordinated changes to their coefficients leave the predictions unchanged. This is an identifiability issue; even uniquely fitted regression coefficients would not, by themselves, establish causal effects.

import numpy as np

rng = np.random.default_rng(7)
n = 200
celsius = rng.normal(loc=20.0, scale=8.0, size=n)
fahrenheit = celsius * 9 / 5 + 32.0              # exactly determined by celsius
y = 2.0 * celsius + 5.0 + rng.normal(scale=1.0, size=n)

Xt = np.column_stack([np.ones(n), celsius, fahrenheit])   # intercept, C, F
print("columns:", Xt.shape[1], "  rank:", np.linalg.matrix_rank(Xt))

w0 = np.linalg.lstsq(Xt, y, rcond=None)[0]       # one least-squares solution
z = np.array([-32.0, -1.8, 1.0])                 # F - 1.8C - 32 = 0, so Xt @ z = 0

for t in (-10.0, 0.0, 10.0):
    w = w0 + t * z
    print(f"  t={t:6.1f}  weights {np.round(w, 3)}"
          f"   largest prediction gap {np.max(np.abs(Xt @ w - Xt @ w0)):.1e}"
          f"   rmse {np.sqrt(np.mean((y - Xt @ w) ** 2)):.6f}")
# columns: 3   rank: 2
#   t= -10.0  weights [319.908  19.723  -9.845]   largest prediction gap 1.9e-13   rmse 0.970066
#   t=   0.0  weights [-0.092  1.723  0.155]   largest prediction gap 0.0e+00   rmse 0.970066
#   t=  10.0  weights [-320.092  -16.277   10.155]   largest prediction gap 1.8e-13   rmse 0.970066

The vector \(z\) is chosen so that \(X_t z = 0\): because Fahrenheit is exactly \(1.8\) times Celsius plus \(32\), that particular combination of the three columns cancels to nothing. Adding any multiple of it to a solution therefore changes the weights without changing a single prediction. The middle column confirms it — the predictions differ by around \(10^{-13}\), which is floating-point dust.

The diagnostic is in the first line of output: three columns, but rank 2. The rank of a table is the number of columns that genuinely add something new. Rank below the number of columns means at least one column is redundant and \(X^{\top}X\) cannot be inverted. The fitted values stay unique; it is the coefficient vector that stops being.

Exact redundancy can be detected by a rank check. With nearly dependent but still independent columns, the least-squares coefficient vector is unique but can be sensitive to noise. The size of that sensitivity also depends on the noise level and the data; the next experiment measures it for a particular design.

The experiment below makes that concrete. It builds two features with a controlled correlation \(\rho\), then refits the model 500 times on fresh noise, and reports how much the coefficients move around. It also prints a quantity called the condition number, which gets a section of its own at the end of this article; for now, read it as one number summarising how close to redundant the columns have become.

import numpy as np

rng = np.random.default_rng(1)
n = 100
x1 = rng.normal(size=n)

for rho in (0.0, 0.9, 0.99, 0.999):
    x2 = rho * x1 + np.sqrt(1 - rho**2) * rng.normal(size=n)
    Xp = np.column_stack([x1, x2])
    betas = []
    for _ in range(500):                       # same Xp, fresh noise each time
        y = 1.0 * x1 + 1.0 * x2 + rng.normal(size=n)
        betas.append(np.linalg.lstsq(Xp, y, rcond=None)[0])
    betas = np.array(betas)
    print(f"rho {rho:5.3f}  cond {np.linalg.cond(Xp):8.1f}"
          f"  sd(beta) {np.round(betas.std(0), 3)}"
          f"  sd(sum) {betas.sum(1).std():.3f}")
# rho 0.000  cond      1.2  sd(beta) [0.116 0.102]  sd(sum) 0.149
# rho 0.900  cond      3.6  sd(beta) [0.226 0.205]  sd(sum) 0.118
# rho 0.990  cond     11.3  sd(beta) [0.683 0.692]  sd(sum) 0.121
# rho 0.999  cond     38.0  sd(beta) [2.144 2.139]  sd(sum) 0.124

As the correlation climbs from 0 to 0.999, the standard deviation of the first coefficient grows from 0.116 to 2.144, roughly eighteenfold; the second coefficient shows a similar increase. The variability of their sum stays almost unchanged over the same range: 0.149 to 0.124.

That contrast is what collinearity does. Severe collinearity makes individual coefficients hard to interpret and highly sensitive to sampling noise, while the fitted values remain comparatively stable. A model can fit well and still report coefficients that flip sign between two random samples drawn from the same population.

So when a stakeholder asks which feature matters most, the coefficients of a linear model cannot answer until the correlation structure has been checked. There is a standard response, but it is built from ideas the next two sections introduce, so it waits until the end.

The covariance matrix

To see why correlated features create small, unstable directions — and how those directions can be found — we need one more object. Center each column of the table by subtracting its own average, call the result \(X_c\), and compute \(S = X_c^{\top}X_c / (n-1)\). This is the covariance matrix, and it is a compact summary of how the features vary and move together.

It is a \(p \times p\) grid — one row and one column per feature. The diagonal entries are the variances: how much each feature spreads out on its own. The off-diagonal entry in row \(i\), column \(j\) is the covariance between features \(i\) and \(j\): positive if they tend to rise together, negative if one rises as the other falls.

Two limits are worth knowing before you rely on it. First, it has \(p(p+1)/2\) distinct entries to estimate, so the number of quantities to pin down grows with the square of the feature count while the evidence grows only with \(n\). Unregularized sample covariance estimates can therefore be unstable when the sample size is small relative to the number of features. After centering, the sample covariance matrix is singular whenever \(p \geq n\), since \(X_c\) has at most \(n-1\) independent rows. Shrinkage estimators or structural assumptions are the usual response in high-dimensional settings. Second, it only sees straight-line relationships. Two features can have a covariance of exactly zero while one is a deterministic nonlinear function of the other, as \(x^2\) is of \(x\) on data spread symmetrically around zero. The dependence runs one way here: \(x\) fixes \(x^2\), while \(x^2\) leaves the sign of \(x\) unknown.

Eigenvectors: the directions a cloud stretches

Picture your data as a cloud of points. Most real clouds are not round — they are stretched, like a rugby ball or a flat pancake. It is natural to ask: along which direction is the cloud longest? Which perpendicular direction has the most spread left? These directions are eigenvectors of the covariance matrix, and each corresponding eigenvalue measures variance along that direction.

For a small calculation, take \(S=\begin{pmatrix}9&0\\0&1\end{pmatrix}\). Multiplying gives \(S(1,0)^{\top}=9(1,0)^{\top}\) and \(S(0,1)^{\top}=1(0,1)^{\top}\). Each vector keeps its direction, so they are eigenvectors with eigenvalues 9 and 1. If \(S\) is a covariance matrix, those values are variances; the standard deviations along the two directions are 3 and 1.

Formally, an eigenvector is a nonzero vector satisfying \(Sv=\lambda v\). For a covariance matrix, the eigenvalues are nonnegative and describe variance along unit eigenvectors. An eigenvector with the largest eigenvalue gives a direction of greatest variance, which is what principal component analysis seeks. If eigenvalues tie, the corresponding directions need not be unique.

There is a second, better route to the same answer. The singular value decomposition factors the centered data directly as \(X_c = U \Sigma V^{\top}\), without ever building \(S\). Each piece holds one thing. \(V\) holds the directions in feature space, one per column — the same eigenvectors as before. \(\Sigma\) is diagonal and holds one non-negative number per direction, the singular value, measuring how far the data extends along it. Multiplying \(U\) by \(\Sigma\) gives each observation’s coordinates along the directions in \(V\). So the variances come from the singular values as \(\lambda_i = \sigma_i^2/(n-1)\). The code below computes both ways and compares.

import numpy as np

rng = np.random.default_rng(2)
A = rng.normal(size=(500, 4)) @ rng.normal(size=(4, 4))
Ac = A - A.mean(0)

S = Ac.T @ Ac / (len(Ac) - 1)
w, V = np.linalg.eigh(S)                       # symmetric eigensolver
U, sv, Vt = np.linalg.svd(Ac, full_matrices=False)

print("eig variances", np.round(np.sort(w)[::-1], 4))
print("svd variances", np.round(sv**2 / (len(Ac) - 1), 4))
# eig variances [1.23506e+01 6.36340e+00 1.14850e+00 1.30000e-03]
# svd variances [1.23506e+01 6.36340e+00 1.14850e+00 1.30000e-03]

The two agree to four decimals, as they must. For a real symmetric matrix, use eigh, which is specialized for this structure and returns real eigenvalues in ascending order.

Look at the four variances: 12.35, 6.36, 1.15, and 0.0013. The last direction carries almost nothing. It is tempting to conclude it is noise and drop it, but be careful — this data genuinely occupies all four dimensions, and that thin direction may be exactly the signal you care about. Deciding whether a small eigenvalue is negligible or important is a modeling judgment, not something the arithmetic can settle.

Condition number, and why scaling is not cosmetic

For a full-column-rank matrix, the 2-norm condition number is \(\kappa_2(X)=\sigma_{\max}/\sigma_{\min}\). A small minimum singular value means some coefficient changes produce only small changes in predictions. Recovering those coefficients from perturbed observations can therefore be sensitive.

For a nonsingular linear system, \(\kappa_2(X)\) bounds the amplification of relative right-hand-side error into relative solution error when the matrix is fixed. This motivates the rough rule that a condition number near \(10^k\) can put about \(k\) decimal digits at risk in a stable floating-point solve. It is not an exact digit count or a universal bound for every least-squares perturbation: sensitivity to changes in the design matrix also depends on the residual and problem geometry. Column scaling and near-dependence can both produce large condition numbers.

This is what makes the normal equations a risky route to a least-squares solution. For a full-column-rank \(X\), \(\kappa_2(X^{\top}X) = \kappa_2(X)^2\), so the normal-equation system can be much more sensitive to roundoff than working directly with the design matrix. Forming \(X^{\top}X\) is not wrong in itself — the covariance matrix earlier in this article is exactly that computation — but using it to solve for \(\beta\) spends accuracy you did not have to spend.

Feature scaling can improve conditioning when column magnitudes differ mainly because of units. Standardization removes those magnitude differences, but does not remove linear dependence or change the magnitude of pairwise correlations. Its effect on the overall condition number depends on the design. Scaling and regularization can therefore serve complementary purposes.

Scaling is not only a numerical convenience, though. It changes what a regularization penalty punishes and what a distance-based method considers close — the very choice we discussed under cosine versus distance. That is why scaling belongs inside the cross-validation fold and must be fitted on training data only: it is part of the model, not a tidying step performed beforehand.

Ridge: a response to instability

Now that eigenvalues and the condition number are on the table, the standard answer to collinearity can be stated properly. With the features and response centered and the intercept handled separately, ridge regression with \(\lambda>0\) solves \((X^{\top}X + \lambda I)\beta = X^{\top}y\) in place of the plain normal equations, which adds \(\lambda\) to every eigenvalue of \(X^{\top}X\). Here \(I\) is the identity matrix: ones on the diagonal and zeros elsewhere. The near-zero directions that made \(\sigma_{\min}\) tiny — and the condition number huge — get a floor, so noise amplification along them is reduced rather than removed. The intercept is normally excluded from the penalty, since shrinking it toward zero would bias the overall level of the predictions.

The penalty discourages large slope coefficients, trading some bias for greater stability, so choosing \(\lambda\) means choosing how much of that trade to accept. Note how this differs from scaling: standardization fixes conditioning that came from units, while ridge addresses dependence between the columns themselves. The two solve different problems and are routinely used together. Ridge gets its own article.

What to carry forward

These ideas form one framework rather than a list. Projection explains why least-squares residuals are orthogonal to the features. Rank explains when the coefficients are determined at all. Singular values explain why nearly dependent columns produce numerical instability, and the covariance matrix connects that instability to the directions PCA finds. The same framework reappears in ridge regression, PCA, logistic regression, and the linear layers inside a neural network.

Exercises

1. Cosine versus Euclidean. Build customers that vary both in total spending and in purchase mix, then rank them by cosine similarity and by Euclidean distance to a query customer. Report where the two rankings disagree.

Expected observation: a customer that ties for the highest cosine similarity while being the farthest away by Euclidean distance.

import numpy as np
q = np.array([10.0, 20.0, 30.0])               # the query customer
cust = {"same mix, 10x":  np.array([100.0, 200.0, 300.0]),
        "same mix, 1x":   np.array([10.0, 20.0, 30.0]),
        "different mix":  np.array([30.0, 20.0, 10.0])}
for name, c in cust.items():
    cos = c @ q / (np.linalg.norm(c) * np.linalg.norm(q))
    dist = np.linalg.norm(c - q)
    print(f"{name:16s} cosine {cos:.4f}   euclidean {dist:8.2f}")
# same mix, 10x    cosine 1.0000   euclidean   336.75
# same mix, 1x     cosine 1.0000   euclidean     0.00
# different mix    cosine 0.7143   euclidean    28.28

2. Projection idempotence. The trace of a square matrix is the sum of the entries on its main diagonal. Verify numerically that the hat matrix satisfies \(H^2 = H\) and \(H^{\top} = H\), and report its trace. Explain what the trace counts.

Expected observation: a trace of 3 for this matrix. In general the trace equals \(\operatorname{rank}(X)\), which matches the column count only when the columns are independent.

import numpy as np
rng = np.random.default_rng(0)
X = rng.normal(size=(50, 3))
H = X @ np.linalg.inv(X.T @ X) @ X.T
print("idempotent", np.allclose(H @ H, H))     # True
print("symmetric ", np.allclose(H, H.T))       # True
print("trace     ", round(float(np.trace(H)), 6))   # 3.0
print("rank      ", np.linalg.matrix_rank(H))       # 3

3. Condition number and lost digits. Solve the same least-squares problem through the normal equations and through lstsq, on a matrix whose condition number you control. Report where the two answers diverge.

Expected observation: the normal-equation error growing far faster than the SVD route’s. Depending on your BLAS build it may end in a singular-matrix failure.

import numpy as np
rng = np.random.default_rng(0)
n = 200
for kappa in (1e2, 1e6, 1e10):
    U, _ = np.linalg.qr(rng.normal(size=(n, 3)))
    V, _ = np.linalg.qr(rng.normal(size=(3, 3)))
    s = np.array([1.0, 1.0 / kappa**0.5, 1.0 / kappa])
    X = U * s @ V.T                            # condition number = kappa
    beta_true = np.array([1.0, 2.0, 3.0])
    y = X @ beta_true

    try:
        normal = np.linalg.solve(X.T @ X, X.T @ y)
        normal_err = f"{np.abs(normal - beta_true).max():.3e}"
    except np.linalg.LinAlgError:
        normal_err = "FAILED"
    svd = np.linalg.lstsq(X, y, rcond=None)[0]
    print(f"kappa {kappa:7.0e}  normal err {normal_err:>9s}"
          f"   lstsq err {np.abs(svd - beta_true).max():.3e}")
# kappa   1e+02  normal err 3.602e-12   lstsq err 4.885e-15
# kappa   1e+06  normal err 1.867e-04   lstsq err 7.905e-13
# kappa   1e+10  normal err    FAILED   lstsq err 2.443e-07