Vectors, matrices, and why GPUs are shaped the way they are

A vector is mechanically just an ordered list of numbers. What makes it useful is the second layer: a vector represents something with both a magnitude and a direction, or equivalently, a point’s coordinates in space. In machine learning, a vector is almost always “a thing’s properties encoded as numbers”: a word’s meaning, an image’s pixels, a user’s preferences. Every embedding, every hidden state, every gradient in a neural network is a vector, and every one of those vectors is stored using numeric formats: binary place value, floating point, and whatever precision it got quantized to.

Dot products measure agreement, not just multiply-and-sum#

The operation you probably memorized as “multiply corresponding entries, sum them” is doing something more specific: it measures directional agreement. For vectors $u, v \in \mathbb{R}^n$:

$u \cdot v = \sum_{i=1}^{n} u_i v_i$

Geometrically this equals $\|u\|\|v\|\cos\theta$, where $\theta$ is the angle between them. Large and positive means the vectors point in similar directions. Near zero means roughly perpendicular, unrelated. Negative means opposite directions. This is the mechanism behind attention scores in transformers: “how much should token A attend to token B” is computed as a dot product between their vector representations.

The norm $\|v\| = \sqrt{\sum_i v_i^2} = \sqrt{v \cdot v}$ is the Pythagorean theorem generalized to $n$ dimensions: the straight-line distance from the origin to the point $v$. Dividing dot product by both norms strips out magnitude entirely, leaving pure direction agreement:

$\cos\theta = \frac{u \cdot v}{\|u\|\|v\|}$

This, not the raw dot product, is what embedding comparisons actually use. Two vectors can have a large dot product just because they’re both long, even if they point in barely-related directions; cosine similarity is what corrects for that.

Worked example. $u = [1, 2, 3]$, $v = [4, 0, -1]$.

dot product: 1(4) + 2(0) + 3(-1) = 4 + 0 - 3 = 1
‖u‖ = √14 ≈ 3.742
‖v‖ = √17 ≈ 4.123
cosine similarity: 1 / (3.742 × 4.123) ≈ 0.0648

The dot product looks like a nonzero, meaningful “1”, but the cosine similarity reveals these vectors are nearly perpendicular: barely related in direction at all. Raw dot product conflates size and direction; cosine similarity separates them.

import numpy as np

u = np.array([1, 2, 3], dtype=np.float32)
v = np.array([4, 0, -1], dtype=np.float32)

dot = np.dot(u, v)
cosine_sim = dot / (np.linalg.norm(u) * np.linalg.norm(v))
print(dot, cosine_sim)   # 1.0 0.0648...

np.dot calls into BLAS (Basic Linear Algebra Subprograms) under the hood rather than looping in Python. That gap, between “the mathematical definition” and “how it actually executes on real hardware,” is the thread the rest of this note pulls on.

A matrix is a transformation, not just a grid#

A matrix is a grid of numbers, but the useful way to think about it is: a matrix is a function that transforms vectors. Feed a vector in, a (possibly differently-sized) vector comes out, and the matrix defines exactly how. A neural network layer, $y = Wx$, is nothing more than “take vector $x$, transform it into vector $y$ using the rule encoded in $W$.”

The rule for multiplying matrices, row of $A$ dotted with column of $B$, isn’t an arbitrary convention. It’s the unique rule that makes applying transformation $B$ then transformation $A$ equal a single matrix $AB$ that does both at once. If matrix multiplication worked elementwise instead, it would stop corresponding to “do this transformation, then that one,” and the entire reason matrices are useful for stacking neural network layers would collapse. The dimension-compatibility rule, $A$‘s columns must equal $B$‘s rows, follows directly: $A$ transforms length-$k$ vectors into length-$m$ vectors, so $B$ must output length-$k$ vectors for $A$ to consume next.

$C_{ij} = \sum_{p=1}^{k} A_{ip} B_{pj}$

$A$ is $m \times k$, $B$ is $k \times n$, $C$ is $m \times n$. $C_{ij}$ is exactly the dot product from the previous section, just relabeled: row $i$ of $A$ as one vector, column $j$ of $B$ as the other, repeated for every $(i,j)$ pair.

Each of the $m \times n$ output entries needs $k$ multiplications and roughly $k$ additions, so total floating-point operations is approximately $2mnk$. This is the formula behind every “how many FLOPs did training this model take” number you’ve ever seen quoted for a frontier model. It isn’t a separate metric; it’s this sum, counted.

Worked example.

$A = \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}, \quad B = \begin{bmatrix} 5 & 6 \\ 7 & 8 \end{bmatrix}$

C₁₁ = (1)(5) + (2)(7) = 19
C₁₂ = (1)(6) + (2)(8) = 22
C₂₁ = (3)(5) + (4)(7) = 43
C₂₂ = (3)(6) + (4)(8) = 50

C = [19 22; 43 50]

def matmul_naive(A, B):
    m, k = A.shape
    k2, n = B.shape
    assert k == k2, "inner dimensions must match"
    C = [[0.0] * n for _ in range(m)]
    for i in range(m):
        for j in range(n):
            C[i][j] = sum(A[i, p] * B[p, j] for p in range(k))
    return C

The triple-nested structure (i, j, and the sum over p) isn’t a simplified stand-in for the formula above, it is the formula, executed one scalar at a time. That’s deliberately the slowest possible correct implementation, and that’s the point: it’s the baseline every hardware optimization in the rest of this note is measured against.

Four implementations of the exact same $C_{ij} = \sum_p A_{ip}B_{pj}$ formula exist in production, and they differ by orders of magnitude in speed for reasons that have nothing to do with the math:

Approach	Speed	Who writes it	Trade-off
Naive triple loop	Extremely slow	Anyone learning the math	Correct, but ignores the memory hierarchy entirely; never used in production
BLAS (OpenBLAS, MKL)	Fast	Numerical library engineers	Tuned for CPU cache hierarchy, no tensor-core acceleration
cuBLAS / CUTLASS (GPU)	Very fast	NVIDIA and kernel engineers	Requires reasoning about tiling, memory hierarchy, and instruction shapes
Hand-written kernel for one exact shape	Fastest possible for that shape	Frontier-lab performance engineers	Highest engineering cost, brittle the moment the shape changes

Every row computes the identical numbers. The gap between rows is the entire subject of the next two sections.

Why GPUs have tensor cores at all#

An ordinary GPU compute core (a CUDA core, in NVIDIA’s terms) does one fused multiply-add at a time: one element of $A$ times one element of $B$, accumulated, one instruction. A tensor core is a specialized unit that instead does a whole tile of the sum in one instruction: multiply a small block of $A$ by a small block of $B$ and accumulate directly into a block of $C$. That single distinction, scalar-at-a-time versus tile-at-a-time, is the entire reason tensor cores exist as separate silicon from general-purpose GPU cores, and it’s the hardware answer to why numeric formats like BF16 and FP8 matter so much: those formats exist largely to feed tensor cores data in the shape and width they’re built to consume.

The tiles themselves have grown fast enough to make the point vivid. A single tensor-core instruction handled a 16×16×16 tile on the Ampere generation; on Hopper that grew to as much as 64×256×16 computed per warp-group in one instruction; on Blackwell, two streaming multiprocessors can now cooperate on a single instruction covering a 256×256×16 tile. Three hardware generations, same underlying formula, three very different answers to “how much of $C_{ij} = \sum_p A_{ip}B_{pj}$ can one instruction finish before it has to ask memory for more data.”

The real bottleneck is usually memory, not arithmetic#

A GPU has a memory hierarchy: global memory (large, comparatively slow), shared memory (smaller, much faster, and explicitly managed by whoever writes the kernel), and registers (tiny, fastest, private to each thread). A naive matrix multiply that re-reads elements of $A$ and $B$ from global memory over and over will hit a memory bandwidth ceiling long before it hits a compute ceiling. That’s called being memory-bound rather than compute-bound, and telling the two apart is arguably the single most useful diagnostic in performance engineering.

Regime	Symptom	Fix
Memory-bound	GPU utilization is low, memory bandwidth is near saturated	Bigger tiles or more reuse, fuse operations to avoid redundant reads, or use a narrower numeric format
Compute-bound	Memory bandwidth has headroom, compute is near saturated	Lower-precision formats with hardware acceleration (FP8 tensor cores), better kernel or tile-shape selection

Tiling is the fix for the first row: load a block of $A$ and a block of $B$ into fast shared memory once, then reuse that block across many output entries before moving on, instead of re-fetching from slow memory for every single multiply. The quantity that captures this trade-off precisely is arithmetic intensity, the ratio of compute to memory traffic:

$I = \frac{\text{FLOPs}}{\text{bytes moved}}$

Low intensity means you’re memory-bound: the GPU spends most of its time waiting on data, not computing. High intensity means you’re compute-bound: data is reused enough that arithmetic, not memory bandwidth, is the limiting factor. Tiling’s entire job is pushing $I$ up by making each loaded byte do more work before it’s discarded. This is the same shape of trade-off as everything in the numeric-formats notes: there, you traded precision for memory footprint; here, you trade memory traffic for arithmetic reuse, and the skill in both cases is knowing which resource is actually your bottleneck before you spend effort optimizing the wrong one.

def arithmetic_intensity(m, n, k, bytes_per_element=4, reuse_factor=1.0):
    flops = 2 * m * n * k
    naive_bytes = (m * k + k * n) * bytes_per_element
    return flops / (naive_bytes / reuse_factor)

print(arithmetic_intensity(4096, 4096, 4096, reuse_factor=1.0))    # low: memory-bound
print(arithmetic_intensity(4096, 4096, 4096, reuse_factor=64.0))   # much higher: tiling's effect

The uncomfortable part: newer hardware doesn’t always win#

Here’s the finding that should make you distrust spec sheets. Going from NVIDIA’s A100 to its B200, raw compute throughput scaled by roughly 7.2x, but memory bandwidth scaled by only about 4.5x. Every generation, GPUs get relatively more compute-rich and relatively more memory-starved, which means arithmetic intensity and tiling strategy matter more with each hardware generation, not less. That alone is worth sitting with: the skill this note is teaching doesn’t get obsolete as hardware improves, it gets more load-bearing.

It gets stranger. Microbenchmarks across matrix sizes from 1024³ up to 8192³ found Hopper consistently outperforming Blackwell on wall-clock runtime for FP8 GEMM, in nearly all tested configurations, with the gap widening at larger sizes and Blackwell showing latency spikes attributed to immature kernel scheduling, this despite Blackwell having the higher theoretical peak on paper. A well-tuned kernel library can get remarkably close to that theoretical peak too: reported figures put one GEMM implementation at around 1,639 TOPS on Hopper (83% of INT8 peak) and 3,654 TOPS on Blackwell (81% of peak), when the kernels are mature enough to reach it.

Put those two facts together and you get the actual lesson: a brand-new chip’s peak-FLOPs number on a spec sheet is a ceiling, not a promise, and the kernels that let real code approach that ceiling typically mature well after the hardware ships. If you’re choosing hardware for a training run, the responsible move isn’t reading the spec sheet, it’s benchmarking your actual matrix shapes on actual current kernels.

Common mistakes#

Assuming “matrix multiply” is one operation with one fixed cost: the real cost depends on the numeric format and the hardware’s tile shape, not just $m$, $n$, $k$, which is exactly how a benchmark like the one above can surprise people who only looked at peak FLOPs.

Not checking dimension compatibility before assuming a bug is elsewhere: a large fraction of real shape-mismatch errors happen at matmul boundaries, and the compatibility rule from earlier, inner dimensions must match, makes them diagnosable instead of mysterious.

Confusing elementwise multiplication with matrix multiplication: * versus @ in NumPy or PyTorch are different operations with different shape rules, and mixing them up is a disproportionately common source of code that runs without error but silently computes the wrong thing.

Try it yourself#

Beginner. By hand, multiply a $2 \times 3$ matrix by a $3 \times 2$ matrix, and confirm the result is $2 \times 2$ using the compatibility rule above.

Intermediate. Rewrite matmul_naive so that a shape mismatch raises a clear error naming which dimensions disagreed, instead of failing on a bare assert.

Advanced. Using arithmetic_intensity above, sweep reuse_factor upward for a fixed $4096^3$ matmul and find roughly where the ratio crosses from clearly memory-bound into plausibly compute-bound. Then ask yourself which side of that line the naive triple loop from earlier actually sits on, and why.

---

The one-sentence version: a dot product measures how much two vectors agree in direction, matrix multiplication is that same operation applied row-by-column so that stacked transformations compose correctly, and how fast either one runs in practice comes down to a single ratio, arithmetic intensity, that tells you whether you’re waiting on memory or waiting on arithmetic, a ratio that gets more decisive, not less, with every new GPU generation. One more thing every linear layer’s matrix quietly is: its own Jacobian, which is the entire reason backpropagation through a transformer is tractable at all.