Linear Algebra Primer

An address, an RGB colour, a pair of GPS coordinates — all of them are lists where the order carries the meaning. (255, 0, 0) is red and (0, 255, 0) is green; (37.78, −122.42) is San Francisco and the swap is open ocean off Antarctica. A vector is that idea with one extra move: read the numbers as coordinates and the list becomes a place. Drag the tip of the arrow and watch its two numbers follow it:

Notice that nothing is stored twice. The arrow is the pair, drawn; its tail sits at the origin by convention, so the tip and the readout are one fact seen two ways. Slot one is the horizontal step and slot two the vertical one — swap them and you are somewhere else entirely.

The picture hands us length for free. Pythagoras is already drawn into the grid: the two components are the legs of a right triangle and the arrow is its hypotenuse. Raise the height and watch ‖v‖ climb:

Because the legs are squared before they are added, the longest component dominates: at the opening height the vertical leg contributes 16 of the total 25, and the arrow measures exactly 5. Flatten it to height 0 and only the horizontal leg is left. That is ‖v‖ = √(Σ vᵢ²) — the Euclidean or L2 norm — and it generalises to any number of components without a symbol changing.

Length and direction come apart cleanly. Divide a vector by its own length and what is left is the unit vector on the same ray: the direction with the size thrown away. Swing the angle and watch v and v̂ stay on one line through the origin:

Watch the readout as it turns: v̂'s two components are exactly the cosine and sine of the angle, and neither ever leaves [−1, 1]. Normalising is the first line of almost every retrieval pipeline, because it is what makes two documents comparable when one of them is ten times longer.

Multiplying by a single number — a scalar — slides the arrow along that ray without ever leaving it. Push the scalar below zero and the arrow flips through the origin instead:

Every multiple of one non-zero vector lands on the same line, so a single vector spans a one-dimensional space. At s = 0 the arrow vanishes: the zero vector has a length but no direction, which is why every formula that divides by a length needs an answer for the case where the length is zero.

Addition is walking. Put v's tail at w's tip and where you end up is the sum. Drag v and watch the parallelogram close:

The parallelogram is the proof that order does not matter: walking w then v ends where the other order ends, because both routes meet at the same far corner. Slot by slot it is just [3, 1] + [−1, 2] = [2, 3] — and the picture is why that stays obvious in 768 dimensions.

Which is where these actually live. A GPT-2 token is a 768-component vector; Llama 3 70B uses 8,192. High dimensions are not simply more of the same: two directions drawn at random there are almost always nearly perpendicular. Raise the dimension and watch the whole cloud of pairwise angles collapse onto 90°:

At d = 2 the angles are spread across the full range; by d = 1024 they sit within a couple of degrees ofperpendicular, and the spread shrinks like 1/√d. That is concentration of measure, and it is the reason cosine similarity is a usable signal at all: in a space whose default relationship is "unrelated, at 90°", a pair that measures 30° apart is saying something.

A matrix is a function that moves the whole plane, and it moves it in the one way that keeps grid lines straight, parallel and evenly spaced. That constraint is severe enough that you only have to say where two arrows go: î, which points at (1, 0), and ĵ, which points at (0, 1). Drag î and watch the lattice follow it:

Notice what you are actually dragging: the matrix's first column. The four numbers in the bracket are the two landing points written side by side — column one is where î goes, column two is where ĵ goes, and nothing else is stored anywhere.

Once the columns are fixed, so is every point. Any vector is v₁î + v₂ĵ, so its image has to be v₁(new î) + v₂(new ĵ) — which, multiplied out, is exactly the row-by-column arithmetic, arrived at from the picture rather than memorised. Pull the slider and watch v ride the bending grid:

Watch that v never leaves its own cell of the lattice. It started one step along each old basis arrow and it finishes one step along each new one: Mv = [2·1 + 1·1, 0.5·1 + 1.5·1] = [3, 2], which is the pair the readout shows once the transform is all the way on. The arithmetic and the picture are the same claim.

"Straight, parallel, evenly spaced" has a name: linearity, or M(u + w) = Mu + Mw and M(su) = s(Mu). Most functions break it — squaring does, thresholding does. Apply the transform and watch the parallelogram survive it:

Because addition survives, you may transform then add, or add then transform, and land on the same point. That is what lets one weight matrix be applied to a whole batch of token vectors in a single operation. It is also why a bias has to be added outside the matrix: a constant shift is precisely the thing a linear map cannot do, since the lattice would have to leave the origin, and a linear map is exactly a map that never does.

Two of the four numbers are enough to say how much area the transform creates or destroys. The unit square becomes a parallelogram, and its signed area is the determinant. Push the slider right until the area reaches zero:

At k = 2 the determinant is zero and the plane is squashed flat onto a line: the second column has become exactly twice the first, so the two arrows no longer span anything, and a whole direction of input now maps to a single point. Past that the determinant goes negative — the transform has flipped the plane over, and ĵ now sits clockwise of î instead of counter-clockwise. A negative determinant is an orientation flip, not an error.

One family changes direction and nothing else. A rotation puts its first column and its second on the unit circle, a quarter turn apart: [cos θ, sin θ] and [−sin θ, cos θ]. Turn it and watch the determinant refuse to move off 1:

Because both columns stay unit-length and stay perpendicular, every length and every angle survives the transform exactly. That is what an orthogonal matrix is, and QᵀQ = I is the same sentence in symbols. Orthogonal matrices are the ones you can apply a thousand times without anything growing or shrinking, which is why they turn up wherever numerical stability is the point.

Transforms compose by multiplication, and the order is part of the answer. AB means apply B first — the notation reads right to left and it catches everyone. Walk both steps, then swap the order and watch the landing point pull away from where the other order puts it:

Matrix multiplication does not commute, and the picture says why: shearing a rotated square is not the same shape as rotating a sheared one. This is the most common silent bug in model code — x @ W and W @ x both typecheck when the matrix is square, both produce finite numbers, and only the loss curve ever complains — the two answers differ exactly the way both steps in either order differ above.

A transform that has squashed nothing can be undone. M⁻¹ is the transform that puts every point back, so M⁻¹M = I. Step through applying M, then undoing it:

The grid returns exactly, which is only possible because det(M) ≠ 0 — nothing was lost on the way out. When the determinant is zero there is no inverse at all: two different inputs already landed on the same output, and no function can send one point back to two. "A singular matrix has no inverse" and "the plane got flattened" are the same sentence.

The definition is arithmetic you could do in your head: multiply matching slots and add them up, v · w = v₁w₁ + v₂w₂ + …. The definition hides the only thing about it that matters. Drag v around w and watch the number:

Notice where the sign turns over. The number is positive while v leans the same way as w, zero at exactly a right angle, and negative once it leans away. It is a measure of agreement, scaled by both lengths.

The geometric reading is a shadow. Drop a perpendicular from v's tip onto w's line; the foot lands ‖v‖ cos θ along it, and the dot product is that length times ‖w‖. Close the angle and watch the shadow grow:

Because v · w = ‖v‖ ‖w‖ cos θ, the definition you compute with and the one you picture are the same number — and the shadow is what makes them one. The slot-by-slot form is what the hardware runs and the cosine form is what you reason with. Neither is an approximation of the other.

One angle matters more than all the rest. Sweep v past a right angle and watch the number cross zero on the way through:

At exactly 90° the dot product is zero whatever the two lengths are. That is the definition of orthogonal, and it is the test the next two sections are built on: two directions carry independent information exactly when their dot product vanishes.

Dot a vector with itself and the angle is zero, so the cosine is 1 and all that is left is the length squared. Drag v and watch the square built on it:

So ‖v‖ = √(v · v), and the norm was never a separate idea: it is the dot product with both arguments the same — the square you just dragged. That is why L2 is the default norm across machine learning — it is the one the inner product hands you for free, and the one whose derivative is a clean 2v.

Divide both lengths out and what is left is pure agreement: cos θ = v·w / (‖v‖‖w‖), which lives in [−1, 1] however big the vectors are. Drag the short arrow and watch the cosine ignore how long you make it:

Watch the two dots on the unit circle. They are where the arrows land once their lengths are divided out, and the cosine knows about nothing else. This is why retrieval systems rank with cosine rather than raw dot product: without the division, a long document out-scores a relevant one for being long.

The cost is why this operation is everywhere: d multiplies and d − 1 adds, one pass, no branches, perfectly vectorisable. Raise the width and read the count off the curve:

At GPT-2's width of 768 one attention score costs 1,535 FLOPs — nothing on its own. But attention computes one for every ordered pair of tokens, so at a 1,024-token context that is 1,048,576 dot products per head per layer. The quadratic everyone complains about is this figure, run once per pair.

A matrix product introduces nothing new: cell (i, j) of AB is the dot product of row i of A with column j of B. Which is exactly why the shapes have to meet in the middle. Pull the inner dimension away from 4 and watch the product stop existing:

Notice that the failure is a shape failure, not a numeric one: at rows of B = 3 there is no answer for the product block to hold, so every framework raises here rather than producing garbage. It is the friendliest error in the whole stack, and the reason (3×4)(4×2) → 3×2 is worth learning as "the inner pair cancels".

With the shapes agreeing, the product fills in one cell at a time. Step through the six cells and watch which row of A and which column of B each one eats:

Each cell costs 4 multiplies and 3 adds, and there are 6 cells: 42 operations for a 3×4 by 4×2. In general an M×K by K×N product is exactly MN(2K − 1), which everyone rounds to 2MNK because the missing MN disappears next to it at real widths. Nothing in that count depends on the values — matrix multiply has no fast path and no early exit, which is precisely why hardware can be built for it.

Real workloads never push one row at a time. Stack N token vectors into a matrix and one product does all of them against the same weights. Raise the batch and watch only the input and output blocks grow:

Because the weights are read once and used N times, the work done per byte of weight fetched goes up with N — which is the whole reason batching makes an accelerator fast. At N = 1 you move 768² weights (1.2 MB in bf16) to perform 1.2 million multiply–adds: one FLOP per byte, and the arithmetic units sit idle waiting on memory. At N = 32 the same bytes do 32× the work.

The other axis is worse, because widening the model grows both sides of the weight block at once. The blocks below are drawn to scale — widen the model and watch the weights grow as the square while the token count stays put:

From 768 to 8,192 is 10.7× in width and 114× in arithmetic: that one projection over 2,048 tokens goes from 2.4 GFLOP to 275 GFLOP, or 0.28 ms on an H100's 989 TFLOP/s of dense BF16 (NVIDIA datasheet, 2023). And there are four such projections per attention block, per layer, before the feed-forward network doubles the bill again.

One more move appears in every attention implementation: the transpose. Aᵀ holds the same numbers with the indices swapped — but the bytes underneath do not move. Flip it and watch what happens to the read order:

Because the array is stored row after row, reading down a column steps a whole row width at a time instead of one element. On a real matrix with 768 floats to a row that stride is 3,072 bytes, so every element lands on its own 64-byte cache line and 15 of the 16 floats fetched are thrown away. The arithmetic is identical and the wall-clock time is not, which is why libraries fuse the transpose into the multiply instead of materialising it.

Almost everything a model does with a vector amounts to: keep the part of it that lies along some direction, throw the rest away. Drag v and watch it split into the part along the line and the part perpendicular to it:

Notice that the two pieces always add back to v, and that the perpendicular one always meets the line at a right angle. In symbols p = (v·a / a·a) a — the dot product picks the amount and the direction supplies the rest. At the opening position p = [1.84, 0.92] and r = [−0.84, 1.68].

The foot is not merely a point on the line, it is the closest one — and that is a claim you can test rather than accept. Slide the point along the line and watch the distance to v:

Watch the right-hand number as you pass the bottom. (v − q) · a is zero exactly where the distance is least and nowhere else: minimising a squared distance and making a residual orthogonal are the same condition. That equivalence is what turns calculus problems into linear algebra ones.

Least squares is this picture with more rows. Seven points, one line, seven residuals — and the best line is the one whose residual vector is orthogonal to everything the line can reach. Tilt the slope and watch the total:

The flat line leaves Σr² = 7.046; the least-squares slope is 0.575, and the nearest the slider gets to it, 0.58, leaves 0.837. The argument is the one you just drove: AᵀA x̂ = Aᵀb says precisely "the residual is orthogonal to every column of A", which is why the normal equations stop looking like a trick pulled out of a hat.

The error is quadratic in each parameter, so the surface it makes has one bottom and no others. Raise the line and watch the parabola pass through its minimum:

The bottom sits at intercept 0.08, and because the curve is a parabola, gradient descent on it cannot get stuck and a direct solve exists at all. That is a property of the squared error specifically: swap it for absolute error and the bottom becomes a crease with no derivative there, which is why L1 regression needs a different algorithm rather than a different constant.

The one thing a projection cannot do is remember. Every vector on a line perpendicular to the direction lands on the same foot, so the map is not reversible. Slide the vector along that line and watch the left-hand number refuse to move:

That set of vectors is the null space, and its existence is what a rank deficiency feels like from the inside. When a design matrix has two columns that are multiples of each other, the fit has a whole line of equally good answers, AᵀA is singular, and a solver either raises — or, worse, quietly returns whichever of them floating-point noise happened to favour.

Take the transform from the last two sections and put one vector through it. Most of the time Av points somewhere v does not. Swing the direction and hunt for the angles where the two arrows lie on top of each other:

There are exactly two, at 26.6° and 135°, and at those the readout hands you λ: 2.50 and 1.00. An eigenvector is a direction the matrix only scales, and its eigenvalue is the scale factor — Av = λv, a whole matrix collapsed into one number along one line.

Seen from the grid, those directions are the two lines the bending cannot move off. Apply the transform and watch the eigen lines hold their angles while everything between them swings:

Which is what A = PDP⁻¹ says: rewrite any vector in the eigenvector basis and the matrix becomes multiplication by two independent numbers. That is the entire payoff, and it is why Aᵏ costs one exponentiation per eigenvalue instead of k matrix multiplies.

It also predicts what repeated application does. Multiply any starting vector by A over and over and the component along the largest eigenvalue outgrows the rest exponentially. Step through eight applications and watch the direction settle:

From 80.5° off the top eigenvector, eight steps bring it to 0.075°: the error shrinks by λ₂/λ₁ = 0.4 every time. This is power iteration — how PageRank was actually computed — and it is the same arithmetic that makes a recurrent network with λ₁ > 1 explode and one with λ₁ < 1 forget.

Not every matrix has such a direction. A rotation moves every vector by the same angle, so none of them survives. Sweep the whole circle and watch the number refuse to reach zero:

The angle sits at 40.0° everywhere, because a real eigenvector would have to be a direction the rotation fixes and there is none. Its eigenvalues are complex, and those two statements are the same statement. So the answer to "does every matrix diagonalise" is no over the reals — and a defective matrix like [[1,1],[0,1]] does not diagonalise over the complexes either.

One family always behaves. When A = Aᵀ the eigenvectors are guaranteed real and mutually perpendicular, whatever the entries are. Push the off-diagonal anywhere and watch the right angle survive:

That right angle is the spectral theorem, and it is why covariance matrices, Gram matrices and AᵀA — symmetric by construction, every one of them — are the matrices every algorithm reaches for. It is also the bridge to the next section, because AᵀA is symmetric even when A is not square.

Eigenvectors need a square matrix and are not guaranteed to exist. The SVD needs neither, and its picture is one sentence: every matrix sends the unit circle to an ellipse. Apply the transform and watch the circle open out:

The two semi-axes are the singular values, σ₁ = 2.558 and σ₂ = 0.977: the most and the least this matrix can stretch anything. ‖Mv‖ ≤ σ₁‖v‖ for every v, with equality only along the long axis — which makes σ₁ the operator norm rather than a metaphor for one.

Read backwards, the ellipse gives the decomposition: the circle got here in three moves, a rotation of the input frame, a scaling along the axes, and a rotation of the result. Step through them:

So M = UΣVᵀ with U and V orthogonal and Σ diagonal and non-negative — and every matrix has one, rectangular or singular or rank-deficient. Notice the two arrows: they started perpendicular on the circle and are perpendicular at every step, because only Σ changes lengths and it changes them along its own axes.

Because Σ is sorted, the small end can be thrown away. Keep only the first k singular directions and what you get is the best rank-k approximation there is. Raise k and watch the shape come back:

At rank 1 the ellipse has collapsed to a segment and 35.7% of the matrix is missing; at rank 2 it is exact. That "best" is the Eckart–Young theorem and it is not a heuristic: no other rank-1 matrix is closer in either the Frobenius or the operator norm. It is the compression argument behind PCA, latent semantic analysis, and every low-rank adapter.

The ratio of the two axes has a name and a bill. Squeeze the short axis toward zero and watch the condition number run away:

κ = σ₁/σ₂ is the factor by which a relative error in the input can be multiplied on its way to the output. It is the number that decides whether a solve is trustworthy — and det ≠ 0 tells you nothing about it, because a matrix can have determinant 1 and a condition number of 10⁸. Drive σ₂ down to 0.016 and κ reaches 160.

The cost is measured in digits. Every factor of ten in κ eats one significant decimal digit of the answer, and a format only has so many. Raise the condition number and watch the digits go:

float32 carries about 7.2 significant decimal digits, so at κ = 10⁶ it has 1.2 left and the answer is noise with a decimal point in it. float64's 16.0 digits survive the same solve with 10.0 to spare. This is why linear solvers report a condition estimate, and why "it ran without an exception" is not evidence that the answer means anything.

The same truncation argument, run at model scale, is why low-rank fine-tuning works. A 4,096-square weight matrix is 16.8M parameters; a rank-r update is two thin blocks. Raise the rank and read the share:

At r = 8 the update is 65,536 parameters — 0.39% of the layer — and it can express any change whose singular values past the eighth are small. That is a bet about the shape of the update rather than its size, and it is exactly the bet LoRA makes.

One last reading of the same theorem. The singular directions of a centred data matrix are the directions the data actually spreads along, which is PCA. Swing the line and watch the spread of the projected points:

The spread peaks at 2.089 at 32°, against a total of 2.280 — one direction carries 91.6% of the variance, so this cloud is very nearly one-dimensional. Keeping the top axes of the covariance keeps the spread; that is PCA in a sentence, and it is the rank truncation above wearing different clothes.

Order. x @ W and W @ x both typecheck on a square weight and both return finite numbers; only the loss curve complains. Conditioning. A determinant far from zero says nothing about κ, so a solve can be numerically worthless without raising. Normalisation. Ranking on raw dot products silently prefers long vectors. Memory order. A transpose is free in index arithmetic and expensive in cache lines.