23 - Linear Algebra 1: Vectors, Matrices, Systems

Machine Learning Engineer course · Module 02 · Mathematics and statistics · Lesson 8

Goal. See linear algebra for what it is: the geometry of space, not tables of numbers. At the end you will be able to read a vector as an arrow and as a row of data, compute and interpret the dot product and the cosine similarity, understand a matrix as a transformation of space (with the determinant as an area factor), solve a linear system with Gaussian elimination written from scratch, recognise an ill-conditioned problem, and derive the normal equations of least squares with many variables.

Myth to dispel. "To solve Ax = b you compute the inverse: x = A⁻¹b." It is the slowest and least accurate method. The system is solved directly (we will do it ourselves, by hand, and then with np.linalg.solve).

Introduction

A dataset is a matrix: rows for examples, columns for features. A linear model is a product X @ w. A neural network is a succession of multiplications by matrices and of non-linearities. To understand what these operations do (and when they fail) you need the geometry underneath. This lesson builds it; the next adds eigenvalues, SVD and PCA.

The path starts from vectors and the dot product, moves to matrices as transformations of space and to the determinant, solves linear systems with Gaussian elimination and tackles rank, inverse and conditioning, to arrive at the normal equations of least squares. The matrix product takes up the composition of functions from sets and the chain rule from Calculus 1; vectorisation in NumPy is its working environment.

Vectors

The analogy: the displacement

A vector is a displacement: "three steps east and two north" is the vector (3, 2), regardless of where you start. It is drawn as an arrow from the origin to the point (3, 2). With d components it is a displacement in a d-dimensional space.

The same object is also a row of data: a customer aged 40, income 35,000, 3 purchases is the vector (40, 35000, 3) in the feature space (lesson 1: the Cartesian product ℝ³). Here geometry merges with machine learning: an example is a point, a model is something that acts on points.

Two operations:

• Sum: the components are added, (3, 2) + (1, 4) = (4, 6): two displacements in sequence.
• Multiplication by a scalar c: stretches (or shortens, or reverses if c < 0) the arrow.

A linear combination c₁v₁ + c₂v₂ + … is the recipe for building new vectors with these two operations. The set of all linear combinations of the vectors v₁, …, v_k is their span (the subspace they generate): for one vector it is a line, for two non-parallel vectors a plane.

import numpy as np

v = np.array([3.0, 2.0])
w = np.array([1.0, 4.0])
print("sum                :", v + w)
print("2.5 * v            :", 2.5 * v)
print("combination 2v - w :", 2 * v - w)

# a customer as a feature vector
cliente = np.array([40, 35_000, 3])
print("shape:", cliente.shape, "| space R^", cliente.size)

Output:

sum                : [4. 6.]
2.5 * v            : [7.5 5. ]
combination 2v - w : [5. 0.]
shape: (3,) | space R^ 3

What the code does. It performs the two elementary operations with NumPy arrays (which do them component by component, as in lesson 5) and shows that a customer is a vector in ℝ³. The expression 2 * v - w is a linear combination with coefficients 2 and −1.

The dot product, the norm and the cosine similarity

Definition and meaning

The dot product of a and b is the sum of the products of the components:

a · b = Σ aᵢ·bᵢ.

It has a geometric interpretation: a · b = |a|·|b|·cos θ, with |a| = √(a · a) the length (norm) of a and θ the angle between the two vectors. Let us prove the correspondence in the plane: write a = |a|(cos α, sin α) and b = |b|(cos β, sin β). Then a · b = |a||b|(cos α cos β + sin α sin β) = |a||b| cos(α − β), by the cosine subtraction formula, and α − β is the angle between the vectors. ∎

Consequences:

• Perpendicular ⇔ a · b = 0 (the cosine of 90° is 0).
• Projection: the "shadow" of b on the direction of a has length (a · b)/|a|.
• Cauchy-Schwarz inequality: since |cos θ| ≤ 1, it holds that |a · b| ≤ |a|·|b|.
• Cosine similarity: cos θ = (a · b)/(|a||b|) measures how much two vectors point in the same direction, ignoring their length: value 1 same direction, 0 perpendicular, −1 opposite.

The most used norms: ‖x‖₂ = √Σ xᵢ² (Euclidean), ‖x‖₁ = Σ |xᵢ| (sum of absolute values, "taxi distance"), ‖x‖∞ = max |xᵢ|. The distance between two points is the norm of the difference.

a = np.array([1.0, 2.0, 3.0])
b = np.array([4.0, 5.0, 6.0])
print("a·b by hand  :", sum(x * y for x, y in zip(a, b)), "| with @:", a @ b, "| np.dot:", np.dot(a, b))
print("norms of a   : L2 =", round(float(np.linalg.norm(a)), 4), "| L1 =", np.linalg.norm(a, 1), "| Linf =", np.linalg.norm(a, np.inf))

cos = (a @ b) / (np.linalg.norm(a) * np.linalg.norm(b))
print("cos θ =", round(float(cos), 4), "| angle in degrees =", round(float(np.degrees(np.arccos(cos))), 2))
print("perpendicular (1,2)·(-2,1) =", np.array([1, 2]) @ np.array([-2, 1]))

# Cauchy-Schwarz on 1000 random pairs
rng = np.random.default_rng(0)
A, B = rng.normal(size=(1000, 5)), rng.normal(size=(1000, 5))
lhs = np.abs((A * B).sum(axis=1))
rhs = np.linalg.norm(A, axis=1) * np.linalg.norm(B, axis=1)
print("Cauchy-Schwarz holds for all 1000 pairs:", bool(np.all(lhs <= rhs)))

Output:

a·b by hand  : 32.0 | with @: 32.0 | np.dot: 32.0
norms of a   : L2 = 3.7417 | L1 = 6.0 | Linf = 3.0
cos θ = 0.9746 | angle in degrees = 12.93
perpendicular (1,2)·(-2,1) = 0
Cauchy-Schwarz holds for all 1000 pairs: True

What the code does. It computes the dot product in three equivalent ways (explicit sum, operator @, np.dot), the three norms of the vector a, and the angle between a and b by inverting the cosine. Then it verifies that (1, 2) and (−2, 1) are perpendicular (product 0) and checks the Cauchy-Schwarz inequality on 1000 random pairs: (A * B).sum(axis=1) is the dot product row by row, and np.linalg.norm(..., axis=1) computes the norms of all the rows in one go.

An application: similar documents

We represent three small texts as vectors that count words, and compare them with the cosine similarity. The advantage over distance: a text twice as long, which repeats the same words, must come out identical in direction.

vocabolario = ["cat", "dog", "fish", "eats", "sleeps"]
doc1 = np.array([3, 0, 0, 2, 1], dtype=float)         # talks about a cat
doc2 = doc1 * 2                                        # the same text, repeated twice
doc3 = np.array([0, 3, 1, 1, 2], dtype=float)         # talks about a dog and a fish

def coseno(x, y):
    return float(x @ y / (np.linalg.norm(x) * np.linalg.norm(y)))

print("cosine(doc1, doc2) =", round(coseno(doc1, doc2), 4), "| Euclidean distance =", round(float(np.linalg.norm(doc1 - doc2)), 3))
print("cosine(doc1, doc3) =", round(coseno(doc1, doc3), 4), "| Euclidean distance =", round(float(np.linalg.norm(doc1 - doc3)), 3))

Output:

cosine(doc1, doc2) = 1.0 | Euclidean distance = 3.742
cosine(doc1, doc3) = 0.276 | Euclidean distance = 4.583

What the code does. doc2 is doc1 with all frequencies doubled: it has a different length but the same direction. The cosine similarity recognises it (1.0), the Euclidean distance considers them far apart (it has length |doc1|). This is why in document retrieval and in embeddings (vectors that represent words or sentences) the cosine is used.

Trade-off. The cosine ignores the magnitude: two customers, one who spends 10 and one who spends 1000 in the same proportions, come out identical. If the magnitude matters, the distance (or a different normalisation) is needed.

Matrices

The analogy: the machine that transforms space

A matrix A of dimensions m × n (m rows, n columns) is a machine that takes a vector x with n components and returns one with m components: y = A x. In particular, if A is square (n × n) it transforms the space into itself: rotates, stretches, squashes, tilts.

The rule to compute it: each component of y is the dot product of a row of A with x:

yᵢ = Σⱼ Aᵢⱼ·xⱼ.

But there is a second way of seeing it, much more illuminating, the column point of view:

A x = x₁·(column 1 of A) + x₂·(column 2 of A) + …

that is A x is a linear combination of the columns of A, with the coefficients given by x. In particular the columns of A are the images of the basis vectors: if e₁ = (1, 0) and e₂ = (0, 1), then A e₁ is the first column and A e₂ the second. Knowing where the basis vectors go, you know the entire transformation.

The determinant

The unit square (sides e₁ and e₂, area 1) becomes the parallelogram with sides A e₁ and A e₂. Its area is the determinant (in absolute value). For A = [[a, b], [c, d]] the columns are (a, c) and (b, d), and the area of the parallelogram they generate is |ad − bc|: hence det A = ad − bc. (The sign says whether the transformation flips the orientation.) In the figure: 2·1.5 − 1·0.5 = 2.5.

If det A = 0 the parallelogram collapses onto a line (area zero): the transformation loses an entire dimension, and is not invertible (you cannot reconstruct the input).

A = np.array([[2.0, 1.0],
              [0.5, 1.5]])
x = np.array([3.0, -1.0])

print("A @ x                        :", A @ x)
print("combination of the columns   :", 3 * A[:, 0] + (-1) * A[:, 1])
print("det A (formula ad - bc)      :", A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0], "| np.linalg.det:", round(float(np.linalg.det(A)), 6))

quadrato = np.array([[0, 1, 1, 0], [0, 0, 1, 1]], dtype=float)       # the 4 vertices, one per column
print("vertices of the square after A:\n", A @ quadrato)

B = np.array([[1.0, 2.0], [2.0, 4.0]])                                # second column = 2 × first
print("det of a matrix with dependent columns:", np.linalg.det(B))

Output:

A @ x                        : [5. 0.]
combination of the columns   : [5. 0.]
det A (formula ad - bc)      : 2.5 | np.linalg.det: 2.5
vertices of the square after A:
 [[0.  2.  3.  1. ]
 [0.  0.5 2.  1.5]]
det of a matrix with dependent columns: 0.0

What the code does. A @ x and the combination of the columns (3·column1 − 1·column2) give the same vector: they are two readings of the same operation. The determinant computed with the formula coincides with NumPy's. Multiplying A by the matrix of the four vertices of the square (one per column) transforms them all together: the vertices of the parallelogram are obtained. The last matrix has the second column equal to double the first: the columns are dependent, the span is a single line and the determinant is zero.

The matrix product: composing transformations

Apply first B and then A to a vector: A(Bx). It is still a linear transformation, and its matrix is the product AB. This is why (AB)x = A(Bx): the matrix product is the composition of functions (lesson 1, where we introduced g ∘ f; and the chain rule of lesson 6). The elements are (AB)ᵢⱼ = Σₖ Aᵢₖ Bₖⱼ, and the dimensions must match: (m × n)(n × p) = (m × p).

The product is not commutative: rotating and then stretching is not the same as stretching and then rotating.

theta = np.radians(90)
R = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta),  np.cos(theta)]])          # 90° counter-clockwise rotation
S = np.array([[2.0, 0.0],
              [0.0, 1.0]])                                # stretches the x direction by 2
v = np.array([1.0, 1.0])

print("rotate first, then stretch  S(Rv):", np.round(S @ (R @ v), 4))
print("stretch first, then rotate  R(Sv):", np.round(R @ (S @ v), 4))
print("S@R == R@S ?", np.allclose(S @ R, R @ S))
print("(S R) v = S (R v):", np.allclose((S @ R) @ v, S @ (R @ v)))

M, N = np.arange(6).reshape(2, 3), np.arange(12).reshape(3, 4)
print("shapes: (2,3)@(3,4) ->", (M @ N).shape)
print("(MN)^T = N^T M^T:", np.array_equal((M @ N).T, N.T @ M.T))

Output:

rotate first, then stretch  S(Rv): [-2.  1.]
stretch first, then rotate  R(Sv): [-1.  2.]
S@R == R@S ? False
(S R) v = S (R v): True
shapes: (2,3)@(3,4) -> (2, 4)
(MN)^T = N^T M^T: True

What the code does. R is the 90° rotation matrix (the basis vectors end up on (0, 1) and (−1, 0)), S stretches by 2 along x. Applied in a different order they give different vectors: (−2, 1) against (−1, 2) (up to rounding), and np.allclose(S @ R, R @ S) is false. Composition is associative: (SR)v = S(Rv). The last lines show the rule of dimensions and the property of the transpose: (MN)ᵀ = NᵀMᵀ (the order is reversed).

Linear systems

What they are

A linear system is a set of equations in which the unknowns appear only multiplied by numbers and added. Example with three unknowns:

 2x +  y −  z =   8
−3x −  y + 2z = −11
−2x +  y + 2z =  −3

In compact form: A x = b, with A the matrix of coefficients, x the vector of unknowns, b the known terms. Read with the columns: "find the combination of the columns of A that produces b". Geometrically each equation is a plane in space, and the solution is the point where the three planes meet.

Gaussian elimination

The idea (already used by Gauss and before him by Chinese mathematicians) is to combine the equations to make the unknowns disappear one at a time, until a triangular form is obtained that is solved "from the bottom":

1. The first equation is used to eliminate x from the others: a suitable multiple of the first is subtracted from each.
2. The second equation (already without x) is used to eliminate y from the following ones.
3. One arrives at an upper triangular system (U x = c): the last equation gives the last unknown directly; then one substitutes backwards (back substitution).

The number used to eliminate, in the working row, is called the pivot. If the pivot is zero or very small a problem arises: rows are swapped so as to use as pivot the value largest in modulus in the column (partial pivoting), to avoid divisions by tiny numbers that amplify rounding errors.

def gauss(A, b):
    """Solves A x = b by Gaussian elimination with partial pivoting. Returns x."""
    A = np.array(A, dtype=float)            # copies: we do not modify the inputs
    b = np.array(b, dtype=float)
    n = len(b)
    for k in range(n):
        # pivoting: bring to position k the row with the largest absolute value in column k
        p = k + int(np.argmax(np.abs(A[k:, k])))
        if abs(A[p, k]) < 1e-12:
            raise ValueError("singular (or nearly singular) matrix: no unique solution")
        if p != k:
            A[[k, p]], b[[k, p]] = A[[p, k]], b[[p, k]]
        # elimination: zero out column k below the pivot
        for i in range(k + 1, n):
            m = A[i, k] / A[k, k]           # multiplier
            A[i, k:] -= m * A[k, k:]
            b[i] -= m * b[k]
    # back substitution
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - A[i, i + 1:] @ x[i + 1:]) / A[i, i]
    return x

A = [[2, 1, -1], [-3, -1, 2], [-2, 1, 2]]
b = [8, -11, -3]
x = gauss(A, b)
print("solution with Gauss:", np.round(x, 6))
print("np.linalg.solve    :", np.round(np.linalg.solve(A, b), 6))
print("check A x - b      :", np.round(np.array(A) @ x - np.array(b), 12))

try:
    gauss([[1, 2], [2, 4]], [3, 6])
except ValueError as e:
    print("system with dependent rows:", e)

Output:

solution with Gauss: [ 2.  3. -1.]
np.linalg.solve    : [ 2.  3. -1.]
check A x - b      : [0. 0. 0.]
system with dependent rows: singular (or nearly singular) matrix: no unique solution

What the code does, step by step. The outer loop chooses, for each column k, the row with the largest pivot (np.argmax on the absolute value of the part below the diagonal) and swaps it into position k if needed. Then, for each row below, it computes the multiplier m = A[i,k]/A[k,k] and subtracts m times the pivot row (from both matrices A and b): the column k below the pivot becomes zero. At the end the back substitution computes x from the last unknown to the first. The result (2, 3, −1) coincides with np.linalg.solve and the residual A x − b is zero. The last example has the second row equal to double the first: no usable pivot, and the function flags it instead of returning garbage.

Cost. Elimination requires about n³/3 multiplications: doubling the number of unknowns multiplies the time by 8. For systems with thousands of unknowns it is still feasible; for millions you need to exploit the structure (sparse matrices, iterative methods like the gradient descent of the last lesson).

Rank, inverse, and when a system does not have a unique solution

The rank of a matrix is the number of linearly independent columns (columns that cannot be obtained as a combination of the others); it coincides with the number of independent rows and is the "dimension" of the space the matrix can reach. For a square n × n matrix:

• rank n: invertible, det ≠ 0, Ax = b has exactly one solution for every b.
• rank < n: singular, det = 0, Ax = b has either no solution or infinitely many.

The inverse A⁻¹ is the matrix that "undoes" A: A⁻¹A = I, with I the identity matrix (the transformation that does nothing). It exists only if A is invertible. But computing it to solve a system is a waste: it costs about three times more than elimination and is less accurate. In practice np.linalg.solve(A, b), never np.linalg.inv(A) @ b.

C = np.array([[1.0, 2.0, 3.0],
              [4.0, 5.0, 6.0],
              [7.0, 8.0, 9.0]])
print("rank of C:", np.linalg.matrix_rank(C), "| det C =", round(float(np.linalg.det(C)), 10))
print("column 3 = 2·column 2 − column 1:", np.allclose(C[:, 2], 2 * C[:, 1] - C[:, 0]))

rng = np.random.default_rng(1)
G = rng.normal(size=(300, 300))
g = rng.normal(size=300)
x1 = np.linalg.solve(G, g)
x2 = np.linalg.inv(G) @ g
print("residual with solve:", f"{np.linalg.norm(G @ x1 - g):.2e}", "| residual with inv:", f"{np.linalg.norm(G @ x2 - g):.2e}")

Output:

rank of C: 2 | det C = 0.0
column 3 = 2·column 2 − column 1: True
residual with solve: 1.40e-12 | residual with inv: 2.06e-12

What the code does. The matrix C (the numbers from 1 to 9) looks "generic" but has rank 2: the third column is a combination of the first two (2·col2 − col1). The determinant is zero (up to rounding) and therefore it is not invertible. In the second part it solves a random 300 × 300 system with solve and with the inverse: the two residuals are comparable on this well-conditioned matrix (here solve is slightly better); the difference is felt in computing time (the inverse costs about triple in number of operations, according to the standard counts: time not measured in this listing) and, above all, on ill-conditioned matrices, where the explicit inverse amplifies errors.

Conditioning

A system can be theoretically solvable (det ≠ 0) and practically unsolvable: it is the case of an ill-conditioned matrix, that is almost singular. Small perturbations of the data (even just the rounding of a float) produce enormous variations in the solution.

The analogy: the lever. Moving a boulder with a long lever requires a small force, but a small error of the hand at the end translates into a large displacement of the boulder. An ill-conditioned matrix is a very long lever: it amplifies errors.

The condition number κ(A) is the ratio between the maximum and minimum "stretching" that A applies to vectors (we will compute it with the SVD in the next lesson). Rule of thumb: about log₁₀ κ decimal digits of accuracy are lost in the solution.

# a nearly singular matrix: two almost equal rows
A_buona = np.array([[1.0, 1.0], [1.0, -1.0]])
A_cattiva = np.array([[1.0, 1.0], [1.0, 1.0001]])
b = np.array([2.0, 2.0])

for nome, A in (("well-conditioned", A_buona), ("ill-conditioned", A_cattiva)):
    x = np.linalg.solve(A, b)
    x_pert = np.linalg.solve(A, b + np.array([0.0, 0.0001]))        # perturb b by 1e-4
    print(f"{nome}: κ = {np.linalg.cond(A):.3e} | x = {np.round(x, 4)} | with perturbed b: {np.round(x_pert, 4)}")

# the degree-14 polynomial of the Matplotlib lesson: the Vandermonde matrix
for grado in (3, 5, 9, 14):
    punti = np.linspace(0, 1, 15)
    V = np.vander(punti, grado + 1)                                    # columns: x^degree, ..., x, 1
    print(f"degree {grado:>2}: condition number of the data matrix = {np.linalg.cond(V):.2e}")

Output:

well-conditioned: κ = 1.000e+00 | x = [2. 0.] | with perturbed b: [ 2.e+00 -1.e-04]
ill-conditioned: κ = 4.000e+04 | x = [2. 0.] | with perturbed b: [1. 1.]
degree  3: condition number of the data matrix = 1.05e+02
degree  5: condition number of the data matrix = 3.16e+03
degree  9: condition number of the data matrix = 4.63e+06
degree 14: condition number of the data matrix = 4.03e+11

What the code does. It solves two almost identical systems. For the "good" matrix perturbing b by 10⁻⁴ changes the solution by that order; for the "bad" one (two almost equal rows) the same perturbation moves the solution by orders of magnitude more. The second part rediscovers the mysterious error of the degree-14 polynomial met in Visualising Data with Matplotlib: it builds the Vandermonde matrix (the columns are the powers of x) used to fit a polynomial to data; its condition number grows rapidly with the degree, up to enormous values for degree 14. Fitting that polynomial means solving a system with that matrix, and rounding errors are amplified by the factor κ: this is why the degree-14 polynomial "exploded". It is the link between linear algebra and floating-point arithmetic.

Remedies. Standardise the features; remove nearly collinear variables; use regularisation (adding a term that penalises large weights, which improves conditioning); use numerically stable methods (QR, SVD) instead of the normal equations.

Least squares with many variables

The problem

With n examples and d features, X is an n × d matrix, w a vector of d weights, y the vector of the n values to predict. We want w such that X w ≈ y. Usually n > d (more equations than unknowns): the system X w = y has no exact solution, and we look for the best in a quadratic sense, minimising

L(w) = ‖y − X w‖² = (y − Xw)ᵀ(y − Xw).

The geometric point of view

Xw is a linear combination of the columns of X: it lies in the plane (subspace) generated by the columns. y in general lies outside that plane. The point of the plane closest to y is its orthogonal projection ŷ, and what remains (the residual r = y − ŷ) is perpendicular to the plane, that is perpendicular to each column of X:

Xᵀ (y − X w) = 0 ⇒ XᵀX · w = Xᵀy.

These are the normal equations (normal in the sense of "perpendicular").

The derivation with the gradient

Let us verify with the tools of lesson 7. Expanding: L(w) = yᵀy − 2wᵀXᵀy + wᵀXᵀXw. The gradient with respect to w is ∇L = −2Xᵀy + 2XᵀXw (obtained as for partial derivatives, term by term: the derivative of wᵀa is a, that of wᵀMw for M symmetric is 2Mw). Imposing ∇L = 0 gives again XᵀXw = Xᵀy: two different roads (geometry and calculus) lead to the same equation. If the columns of X are independent, XᵀX is invertible and w = (XᵀX)⁻¹Xᵀy.

rng = np.random.default_rng(5)
n, d = 500, 3
X = rng.normal(size=(n, d))
w_vero = np.array([2.0, -1.0, 0.5])
y = X @ w_vero + rng.normal(0, 0.3, size=n)

# 1) normal equations
w_normali = np.linalg.solve(X.T @ X, X.T @ y)
# 2) stable solver (uses SVD/QR)
w_lstsq = np.linalg.lstsq(X, y, rcond=None)[0]
# 3) scikit-learn
from sklearn.linear_model import LinearRegression
w_sk = LinearRegression(fit_intercept=False).fit(X, y).coef_

print("true weights      :", w_vero)
print("normal equations  :", np.round(w_normali, 4))
print("np.linalg.lstsq   :", np.round(w_lstsq, 4))
print("scikit-learn      :", np.round(w_sk, 4))

residuo = y - X @ w_lstsq
print("Xᵀ r (must be ≈ 0, residual orthogonal to the columns):", np.round(X.T @ residuo, 8))
print("R² =", round(1 - (residuo ** 2).sum() / ((y - y.mean()) ** 2).sum(), 4))

Output:

true weights      : [ 2.  -1.   0.5]
normal equations  : [ 1.9888 -1.0144  0.5059]
np.linalg.lstsq   : [ 1.9888 -1.0144  0.5059]
scikit-learn      : [ 1.9888 -1.0144  0.5059]
Xᵀ r (must be ≈ 0, residual orthogonal to the columns): [-0. -0. -0.]
R² = 0.981

What the code does. It generates 500 examples with 3 features and true weights (2, −1, 0.5), plus noise. It computes the weights in three ways (normal equations, lstsq, scikit-learn) and obtains the same result, close to the true weights. The last part verifies the perpendicularity condition: Xᵀr is almost zero (up to rounding), that is the residual is orthogonal to all the columns, exactly as in the drawing. fit_intercept=False because here the model has no intercept (to include it just add a column of ones to X).

A warning: the normal equations double the conditioning

The conditioning of XᵀX is the square of that of X: κ(XᵀX) = κ(X)². If κ(X) = 10⁶, κ(XᵀX) = 10¹² and almost all the digits of a float64 are lost. This is why, when the columns are almost collinear, lstsq (based on stable factorisations) is preferred.

Xc = np.column_stack([X[:, 0], X[:, 1], X[:, 0] + 1e-6 * rng.normal(size=n)])     # third column ≈ first
print("κ(X)   =", f"{np.linalg.cond(Xc):.2e}", "| κ(XᵀX) =", f"{np.linalg.cond(Xc.T @ Xc):.2e}", "(about the square)")
yc = Xc @ np.array([1.0, 1.0, 1.0]) + rng.normal(0, 0.3, size=n)
print("coefficients (lstsq):", np.round(np.linalg.lstsq(Xc, yc, rcond=None)[0], 3))
print("sum of the two collinear coefficients:", round(float(np.linalg.lstsq(Xc, yc, rcond=None)[0][[0, 2]].sum()), 3), "(the sum is well determined, the individual ones are not)")

Output:

κ(X)   = 1.82e+06 | κ(XᵀX) = 3.33e+12 (about the square)
coefficients (lstsq): [ 1.700241e+03  1.004000e+00 -1.698237e+03]
sum of the two collinear coefficients: 2.004 (the sum is well determined, the individual ones are not)

What the code does. It builds a matrix in which the third column is the first plus a tiny noise (collinearity): the two columns carry almost the same information. The condition number explodes, and that of the normal equations is about the square. The coefficients of lstsq on the two collinear columns turn out unstable (they change a lot with small variations of the data), but their sum is well determined: the data can say how much "the common direction" weighs but not how to split it. Whoever interprets the coefficients as the "importance" of a variable must take this into account.

Typical mistakes (symptom, cause, remedy)

• * instead of @. It multiplies element by element instead of the matrix product (lesson 5). Remedy: @, and check the shapes.
• Vectors (n,) and (n, 1) confused. Broadcasting produces (n, n) matrices without error. Remedy: assert on the shapes.
• Inverting to solve. inv(A) @ b is slow and imprecise. Remedy: solve.
• Testing det == 0 for singularity. A determinant 1e-17 is "zero" but not exactly zero. Remedy: rank (matrix_rank) or conditioning.
• Ignoring collinearity. Enormous and unstable coefficients. Remedy: remove redundant variables, regularise.
• Order of multiplications. AB ≠ BA; (AB)ᵀ = BᵀAᵀ. Remedy: reason about the dimensions.
• Normal equations with ill-conditioned matrices. Remedy: lstsq, QR, SVD.

Three mistakes out of seven concern the shapes of arrays (* versus @, (n,) versus (n, 1), order of the factors) and are prevented by writing the shape next to every variable. The other four (inverse, determinant, collinearity, normal equations) have the same root: treating a nearly singular matrix as if it were well conditioned, and to find out it is enough to compute the condition number before trusting the result.

How to work in practice

Think in shapes: write next to every array its shape. Use the solvers (solve, lstsq) and not the inverses, check the conditioning with np.linalg.cond before trusting a solution and standardise the features to improve it. Verify with the residual (‖Ax − b‖ small and Xᵀr ≈ 0 for least squares) and draw in 2D before moving to high dimensions.

The essential vocabulary. A vector is both a displacement and an ordered list of numbers; a linear combination is a sum of scaled vectors and the span the set of all the possible ones. The dot product is Σ aᵢbᵢ = |a||b|cos θ and the norm the length. The matrix represents a linear transformation, the determinant the factor by which it multiplies areas, the rank the number of independent columns, the inverse the matrix that undoes A. The pivot is the value used to eliminate an unknown, the condition number measures how much the solution amplifies errors on the data, the normal equations are XᵀXw = Xᵀy and collinearity is the near-proportionality between features.

Mini-project

1. Implement cosine_similarity_matrix(X) that returns the n × n matrix of cosine similarities between all the rows, without loops (hint: normalise the rows and use a single matrix product X_norm @ X_norm.T). Verify it on 5 invented documents.
2. Extend gauss so that it also returns the determinant (the product of the pivots, with the sign changed at every row swap) and compare it with np.linalg.det.
3. With scikit-learn's breast-cancer data, solve the linear regression of the first feature index on the others with the normal equations and with lstsq: compare the condition number of X and of XᵀX, with and without standardisation.

The first point is an immediate application of the dot product: normalising the rows turns all the similarities into a single matrix product, and it is the trick used to compare thousands of embeddings. The second extends the Gauss code with the determinant, and the third lets you see with your own eyes that the conditioning of the normal equations is the square of that of X.

Challenge

• Easy. Prove that two non-zero vectors are perpendicular if and only if |a + b|² = |a|² + |b|² (generalised Pythagorean theorem).
• Medium. Prove that det(AB) = det(A)·det(B) for 2 × 2 matrices and interpret it geometrically (areas that multiply).
• Hard. Prove that κ(XᵀX) = κ(X)² using the fact that the "stretching lengths" of XᵀX are the squares of those of X (we will justify them with the SVD).
• Expert. Implement the orthogonal projection of y onto the span of the columns of X with the Gram-Schmidt method (orthogonalisation) and compare it with X @ lstsq(X, y).
• Enterprise. Take a linear model from your organisation with at least 10 features. Compute κ(X) before and after standardisation, identify the pairs of features with correlation > 0.9, and report how much the coefficients change by removing one of them.

The first two challenges are short proofs that can be checked with a numerical case; the third can be closed only with the SVD of the next lesson, and the fourth forces you to write the orthogonal projection instead of calling lstsq. The last is the most useful for work: compare the coefficients with and without a nearly duplicated variable, and you will see how little the "importance" read from the weights is worth.

What you can do now

You can read a vector as an arrow and as a row of data, use the cosine to compare directions, interpret a matrix from its columns and the determinant as an area factor, and solve a system with Gaussian elimination with pivoting instead of with the inverse. You can recognise an ill-conditioned system with np.linalg.cond and derive the normal equations both with geometry and with the gradient. The next lesson, eigenvalues, SVD and PCA, gives the precise meaning of the condition number; least squares with many variables returns in Multivariate regression.

The same series also includes: 22 - Calculus 2: Gradient, Integrals, Gradient Descent (the previous lesson), 25 - Matrix Calculus and Automatic Differentiation: From the Jacobian to an Autograd Written From Scratch (which refers back to this lesson).