Linear algebra in NumPy, by doing

Linear algebra is arithmetic on whole lists of numbers at once. Since a model represents everything — a word, an image, a user — as a list of numbers, it is the arithmetic of the entire field. This section covers the handful of operations that actually recur: combining two lists into a single number that says how alike they are, using a grid of numbers to move points around, and finding the line that best fits a cloud of them. Each is done by hand on small numbers first and then in NumPy, so the meaning and the code arrive together.

A vector, and the dot product#

import numpy as np

a = np.array([3.0, 4.0])
b = np.array([4.0, 3.0])

print(a @ b)                  # 3*4 + 4*3 = 24.0   the dot product
print(np.linalg.norm(a))      # sqrt(3^2 + 4^2) = 5.0   the length
cos = (a @ b) / (np.linalg.norm(a) * np.linalg.norm(b))
print(cos)                    # 0.96: pointing in nearly the same direction

The dot product multiplies matching entries and adds them up. Divided by both lengths, it gives the cosine of the angle between the vectors: 1 for the same direction, 0 for unrelated, -1 for opposite. This is exactly how embedding similarity works.

Check the arithmetic against the output before going further. a is [3, 4] and b is [4, 3], so the dot product is 3 times 4 plus 4 times 3, which is 24. The length of a is the square root of 9 plus 16, which is the square root of 25, which is 5, and b has the same length by symmetry. So the cosine is 24 divided by 25, which is 0.96 — close to 1, because the two vectors are nearly the same direction but not quite.

The reason to write this as a @ b rather than a Python loop is not elegance. NumPy hands the work to compiled code that operates on the whole array at once, and the difference is not subtle.

import time

rng = np.random.default_rng(0)
u = rng.standard_normal(1_000_000)
v = rng.standard_normal(1_000_000)

start = time.perf_counter()
total = 0.0
for i in range(len(u)):
    total += u[i] * v[i]
loop_seconds = time.perf_counter() - start

start = time.perf_counter()
fast = u @ v
numpy_seconds = time.perf_counter() - start

print(round(loop_seconds, 3), round(numpy_seconds, 6))   # 0.219 0.000359
print(round(loop_seconds / numpy_seconds))               # 611, and 464 on a second run
print(total, fast)   # -204.34915522872583 -204.34915522874348

A matrix moves points#

theta = np.pi / 2
rotate = np.array([[np.cos(theta), -np.sin(theta)],
                   [np.sin(theta),  np.cos(theta)]])
stretch = np.array([[2.0, 0.0],
                    [0.0, 0.5]])

p = np.array([1.0, 0.0])
print(rotate @ p)             # [0, 1]: turned a quarter circle
print(stretch @ p)            # [2, 0]: doubled along x
print(stretch @ rotate @ p)   # rotate first, then stretch: [0, 0.5]

Multiplying a matrix by a vector transforms it: rotates, stretches, projects. Multiplying two matrices makes one transformation that does both, and order matters. A neural network layer is a learned transformation of this kind followed by a non-linearity.

Two details in that output are worth stopping on. The first entry of the rotation result prints as 6.123234e-17 rather than 0, because the cosine of pi over 2 is not exactly zero in floating point; treat anything of that size as zero and never test a float with an equality. And read the last two lines as composition applied right to left: stretch @ rotate @ p rotates first and then stretches, giving [0, 0.5], whereas rotate @ stretch @ p stretches first and gives [0, 2]. Same two operations, results four times apart. This is the same fact as A @ B not equalling B @ A, seen geometrically.

print(rotate @ stretch @ p)      # [1.2246468e-16 2.0000000e+00]   stretch first
print(stretch @ rotate @ p)      # [1.2246468e-16 5.0000000e-01]   rotate first

One more thing about that stretch matrix is worth noticing, because it is the bridge to everything in the next section. It leaves the direction [1, 0] pointing exactly where it was and only doubles its length, and it does the same to [0, 1], halving it. Those two directions are the eigenvectors of that matrix, and 2 and 0.5 are its eigenvalues. Eigenvectors are not an abstraction added later; they are the directions a transformation does not turn.

Matrix multiplication, one entry at a time#

A = np.array([[1, 2], [3, 4], [5, 6]])    # 3 x 2
B = np.array([[7, 8, 9], [10, 11, 12]])   # 2 x 3

C = np.zeros((3, 3))
for i in range(3):
    for j in range(3):
        C[i, j] = A[i, :] @ B[:, j]       # row i of A dotted with column j of B

print(np.allclose(C, A @ B))              # True

Solving equations#

# 2x + y = 5  and  x - y = 1
A = np.array([[2.0, 1.0], [1.0, -1.0]])
b = np.array([5.0, 1.0])
print(np.linalg.solve(A, b))     # [2, 1]

Substitute to check: 2 times 2 plus 1 is 5, and 2 minus 1 is 1. Both hold, so [2, 1] is the answer. Writing the system as a matrix and a right-hand side is what turns two lines of school algebra into one call that scales to ten thousand unknowns.

It fails in exactly one interesting way, and the failure is worth meeting deliberately. If one equation is a multiple of another, they say the same thing and there is no unique answer.

singular = np.array([[1.0, 2.0], [2.0, 4.0]])    # row 2 is row 1 doubled
try:
    np.linalg.solve(singular, np.array([1.0, 2.0]))
except np.linalg.LinAlgError as e:
    print(e)        # Singular matrix

Never invert a matrix to solve a system. np.linalg.inv(A) @ b gives the same answer in theory and is slower and less accurate in practice, because it computes more than you asked for and loses precision doing it. Reach for solve, or for lstsq when the system is overdetermined.

Least squares: when there is no exact answer#

With more equations than unknowns, such as fitting a line through a hundred noisy points, there is usually no exact solution. Least squares finds the one that minimises the total squared error. This is linear regression.

rng = np.random.default_rng(0)
x = rng.uniform(0, 10, 100)
y = 3 * x + 2 + rng.normal(0, 1, 100)       # true slope 3, intercept 2, plus noise

X = np.c_[x, np.ones_like(x)]               # a column of ones gives the intercept
coef, *_ = np.linalg.lstsq(X, y, rcond=None)
print(coef)                                  # [2.9869275  1.99254002]

The column of ones is worth understanding rather than copying. Without it the model is y equals slope times x, a line forced through the origin. The ones column gives every row the same constant input of 1, so its weight is free to be whatever the intercept needs to be. Every framework does this internally under the name bias, and it is the same trick.

There is a closed-form solution, the normal equations, and it is worth computing once to see that lstsq is not doing anything mysterious. It is also worth not using afterwards, for the reason above: it inverts a matrix.

normal = np.linalg.inv(X.T @ X) @ X.T @ y
print(np.abs(coef - normal).max())      # about 1e-14: the same answer

residual = y - X @ coef
print(X.T @ residual)      # about [4e-13, 4e-14]: zero, to floating point

What to skip#

Determinants beyond knowing that zero means not invertible, computing an inverse by hand, row-reducing by hand, and every identity involving cofactors: none of it will appear in your work or in an engineering interview. What is worth having at your fingertips is this section's five calls — the dot product, the norm, the matrix product, solve and lstsq — and the ability to say what each one means geometrically when asked.