Arrays, shapes, and the first linear model

Course map · NumPy foundations · Original session 001

A neural network does arithmetic on arrays, but an array is not a complete model of intelligence—or even a complete model of a program. Start with a smaller promise: represent two observations, apply the same rule to each, and know what every axis means.

One prediction, then a batch

Suppose each row contains hours studied and exercises completed. Our invented score rule is $z=2x_1-x_2+0.5$. With rows $(1,2)$ and $(3,4)$, the predictions are $0.5$ and $2.5$. These weights are chosen, not learned; their units depend on the features.

Write observations as $X\in\mathbb R^{B\times D}$, weights as $W\in\mathbb R^{D\times H}$ and biases as $b\in\mathbb R^H$. Then $Z=XW+b$ has shape $(B,H)$. Each entry is a sum over features:

The repeated index $k$ disappears; batch index $i$ and output index $j$ survive. This is why matrix multiplication is not elementwise multiplication. For two outputs, choose columns $(2,-1)$ and $(1,3)$. The first row produces $(0,7)$ before bias; the second produces $(2,15)$. A layer can therefore learn several feature combinations at once.

The bias arrow is sharing, not a second observation. Draw axis names beside dimensions before debugging values.

Broadcasting is a rule, not a guess

NumPy's broadcasting rules compare dimensions from the right: equal sizes or a size of one are compatible. Missing leading dimensions behave as ones. Thus $(2,3)+(3,)$ yields $(2,3)$, but $(2,3)+(2,)$ fails. Most dangerously, $(B,1)-(B,)$ yields $(B,B)$, not one error per observation. When $B=1$, the mistake can hide in tests.

A reduction removes an axis unless keepdims=True. To subtract each row's mean, reduce the feature axis and retain its singleton dimension. Do not memorize “axis zero means columns”; ask which axis is being collapsed.

import numpy as np
 
X = np.array([[1., 2.], [3., 4.]])
W = np.array([[2., 1.], [-1., 3.]])
b = np.array([.5, -1.])
Z = X @ W + b
np.testing.assert_allclose(Z, [[.5, 6.], [2.5, 14.]])
centered = X - X.mean(axis=1, keepdims=True)
np.testing.assert_allclose(centered.mean(axis=1), [0., 0.])
 
pred = Z[:, :1]             # (2, 1)
y = np.array([0., 2.])      # (2,)
assert (pred - y).shape == (2, 2)  # legal but wrong objective
error = pred - y[:, None]
assert error.shape == (2, 1)
print('predictions', Z.tolist(), 'squared error', float((error**2).mean()))

Expected squared error: 0.25. Shape assertions are part of the mathematical specification, not decorative defensive programming.

Arrays are not always contiguous

An ndarray couples storage with shape, dtype and strides. A transpose can describe the same data with different strides. A later reshape may need a copy. Basic slicing usually gives the array view you expect; advanced indexing returns a copy. NumPy's copies-and-views guide explicitly distinguishes these operations. Test overlap with shares_memory; b.base is a is not a reliable test when a is itself a view.

import numpy as np
 
base = np.arange(12, dtype=np.float64)
a = base.reshape(3, 4)
assert np.shares_memory(a, a.T)
assert np.shares_memory(a, a[:, 1:3])
assert not np.shares_memory(a, a[[0, 2]])
flat = a.T.reshape(-1)
assert not np.shares_memory(a, flat)  # this layout needs a copy
assert a.nbytes == 96
print('shape', a.shape, 'strides', a.strides, 'bytes', a.nbytes)

A broadcast need not copy the smaller operand, but the result still occupies memory. Subtracting an array of 50,000 vectors from another 50,000 vectors via a pairwise broadcast can allocate billions of entries. Blocking the outer calculation is then sensible; “never use Python loops” is not a performance law.

Geometry without mysticism

For $u=(3,4)$, $\|u\|_2=5$. For $v=(4,-3)$, $u^Tv=0$. Their nonzero lengths make the right-angle interpretation valid. The zero vector has zero dot product with everything but no defined direction; cosine similarity is undefined there. A dot product between embeddings is a learned score, not proof that two sentences mean the same thing.

Exercises and answer sketches

1. Predict the shape of (4,1,3) + (2,3). Answer: align to (4,1,3) and (1,2,3); output (4,2,3).
2. Why can subtracting labels of shape (B,) silently corrupt a binary classifier with output (B,1)? Answer: it compares every prediction with every label. Match shapes before reducing.
3. How many bytes does a dense (1024,1024) float64 result need, excluding overhead? Answer: $1024^2\times8=8{,}388{,}608$ bytes, or 8 MiB.
4. Replace the bias with one separate bias per training row. What changes? Answer: the rule can memorize row identity and cannot directly apply to an unseen batch. Parameters should express the intended sharing.

Exit test: annotate X @ W + b with all three shapes and compute its first row without running it. Next, DL002 gives these arrays geometric meaning; gradients come later in the sequence.

The following sections develop the operations used in the opening linear-model example; the sequence continues with vectors and matrices in DL002.

The story we're starting with

The small linear rule above is the destination of this first pass. Now develop the array operations underneath it, rather than treating their syntax as magic. An ndarray is a programming data structure; a matrix is a two-dimensional mathematical object; “tensor” has both mathematical and framework-specific meanings. They are related, not universally interchangeable. A complete language model also needs nonlinearities, normalization, data and an objective—not just large matrix products.

Prerequisites

Python indexing, loops and functions are sufficient. Predict each small example before executing it. The experiments here use NumPy only and no large allocations. The GPU and research sections are conceptual previews with primary sources, not instructions to rent hardware or claims of reproduced benchmarks.

1 · What is a NumPy array, really?

A numerical ndarray combines a data buffer with metadata: shape, dtype, strides and an address/offset interpretation. For an ordinary dense array of four int64 values, the element storage is 32 bytes. nbytes counts the logical element bytes, not the Python header, all retained parent storage or total process memory. For object arrays it counts references, not all objects reachable through them.

A Python list stores references to Python objects; a fixed-width numerical array stores values in an implementation suitable for bulk operations. That helps NumPy avoid per-element Python interpretation, but does not imply a universal 100× speedup or that every ndarray is contiguous. Views may step through a buffer with gaps, reversed strides or zero strides from broadcasting.

1.1 · The stride trick — how one array can be many shapes

For a = np.arange(12, dtype=np.int64).reshape(3, 4), strides (32, 8) mean one row step moves 32 bytes and one column step moves eight. The element (i,j) is reached by adding i*32 + j*8 to the first element's address. For a.T, shape becomes (4,3) and strides (8,32). The storage is shared even though a.T.base is a need not be true: a may itself be a view of an earlier owner. The protected test above uses np.shares_memory for exactly this reason.

A view is possible when the requested indexing order can be represented by shape, offset and strides over the same storage. Equal element count is necessary for reshape but does not guarantee a view. Contiguity is not simply “strides decrease”: singleton axes, negative strides and traversal order matter. Use flags.c_contiguous and flags.f_contiguous, and test sharing when ownership matters.

1.2 · Same story on the GPU (a preview of S014–S015)

PyTorch's view documentation says view operations share storage and move no data; transpose is an example that can produce a non-contiguous tensor. It also says reshape may return a view or a new tensor, while contiguous() returns the same tensor if already contiguous and otherwise copies. That is the useful contract. The previous explanation naming a particular CUDA copy routine, a fixed header size and a guaranteed multi-fold speed difference went beyond that contract and is not retained.

A metadata-only operation may affect the cost of the next operation through layout. To investigate, record shapes/strides and profile the actual backend with representative inputs. This lesson does not execute PyTorch, CUDA or a GPU profiler.

2 · Making arrays (five ways)

2.1 From a Python list

np.array([1,2,3], dtype=np.int64) has shape (3,), one axis and 24 logical data bytes. Specify a dtype when your arithmetic or memory contract depends on it; a platform/version default is not a portable guarantee. np.array(3) is a zero-dimensional array with shape (), unlike the length-one array np.array([3]).

2.2 A 2D array — nested lists become a matrix

np.array([[1,2,3],[4,5,6]], dtype=np.float32) has shape (2,3), six values and 24 data bytes. A ragged list does not automatically become a regular numerical matrix. Establish how to pad, pack or otherwise represent variable-length data instead of assuming all nested lists have the same shape.

A one-dimensional (3,) array has no separate row or column orientation. Add an axis to choose (1,3) or (3,1). Image and sequence axes also need semantic names: a convention such as (batch, channels, height, width) is chosen by an API, not enforced by the bytes.

2.3 Zero-filled, one-filled, uninitialized

Use np.zeros((2,3)) for initialized zeros, np.ones((2,3), dtype=np.float32) for initialized float32 ones, and np.empty((2,3)) only when you will overwrite every numerical entry before reading it. Empty numerical arrays contain unspecified values; a convenient-looking zero from one run is not evidence of initialization. Do not assume empty is faster in every workload or print its contents as meaningful data.

2.4 Ranges

np.arange(0,6,2) gives 0, 2, 4: the stop is excluded. np.linspace(0,1,5) gives five evenly spaced samples including both endpoints by default. Floating-point step arithmetic makes endpoint expectations for arange delicate; identical inputs are not randomly assigned different lengths. Use a specified count when that is the intended contract.

np.array is not a general generator-consumption API. For (i*i for i in range(4)), use np.fromiter(..., dtype=np.int64, count=4) to obtain a numerical array. Supplying a count also lets the implementation preallocate.

2.5 Random arrays (the workhorse)

Create a local rng = np.random.default_rng(7) and pass it to code needing randomness. rng.normal(size=(2,3)) samples a distribution; six samples need not have sample mean zero or standard deviation one. A seed makes a controlled sequence reproducible within the relevant environment, not every parallel training run or every NumPy release. NumPy documents no cross-version bit-stream compatibility guarantee for Generator.

The older np.random.seed interface is a legacy global-state interface, not an API this lesson can simply declare removed or deprecated. A local Generator avoids accidental coupling between unrelated examples. Record NumPy version, seed and draw order when exact reproduction matters.

import numpy as np
 
a = np.array([1, 2, 3], dtype=np.int64)
b = np.array([[1, 2, 3], [4, 5, 6]], dtype=np.float32)
assert (a.shape, a.nbytes) == ((3,), 24)
assert (b.shape, b.nbytes) == ((2, 3), 24)
assert np.array(3).shape == ()
assert np.array([3]).shape == (1,)
z = np.zeros((2, 3))
o = np.ones((2, 3), dtype=np.float32)
e = np.empty((2, 3)); e[:] = 5
assert np.all(z == 0) and np.all(o == 1) and np.all(e == 5)
assert np.array_equal(np.arange(0, 6, 2), [0, 2, 4])
assert np.allclose(np.linspace(0, 1, 5), [0, .25, .5, .75, 1])
g = np.fromiter((i*i for i in range(4)), dtype=np.int64, count=4)
assert np.array_equal(g, [0, 1, 4, 9])
r1, r2 = np.random.default_rng(7), np.random.default_rng(7)
assert np.array_equal(r1.normal(size=(2, 3)), r2.normal(size=(2, 3)))
print('creation, dtype, byte-count and seeded-draw checks passed')

3 · Indexing and slicing — where the shapes come from

3.1 The one-dimensional case

For a=np.arange(10), a[2:5] gives a view containing 2, 3, 4; a[::2] selects alternate values; a[::-1] reverses traversal with a negative stride. a[0] is a NumPy scalar, not an array slice you can mutate to change the parent. The original statement that every indexing example was a view was too broad.

3.2 The two-dimensional case — this is where beginners get lost

For shape (3,4), b[0] and b[0,:] have shape (4,); b[0:1,:] has (1,4); b[:,0] has (3,); b[1:3,1:3] has (2,2). Basic integer indices consume axes, while slices preserve them with their selected length. For a (5,3,4) array, a[1:4,:,::-1] has (3,3,4)—reversal changes order, not length.

3.3 The three-dot ellipsis ...

The ellipsis fills the unspecified axes with full slices. For (2,3,4), x[...,2] has shape (2,3), while x[0,...] has (3,4). Work out which original axes remain; do not count punctuation as dimensions.

3.4 Boolean masks — the killer feature

a[a>5] retrieves selected elements in a new array. But a[a>5]=0 assigns into a itself: advanced-index retrieval and indexed assignment are different operations. A mask used for data selection also changes the number of samples, so keep labels aligned with the selected rows.

3.5 Fancy indexing with integer arrays

a[[0,2,4]] gathers those positions and returns a copy. Repeated indices gather repeated values. Indexed += with repeated positions is not a general scatter-sum: use np.add.at when each repeated occurrence must contribute. This distinction will matter when several tokens share one embedding row.

import numpy as np
 
x = np.arange(24).reshape(2, 3, 4)
assert [x[0].shape, x[:, 1].shape, x[..., 2].shape,
        x[0, 1, :].shape, x[:, :, 1:3].shape] == [
            (3, 4), (2, 4), (2, 3), (4,), (2, 3, 2)]
a = np.arange(60).reshape(5, 3, 4)
assert a[1:4, :, ::-1].shape == (3, 3, 4)
a = np.arange(10)
v = a[2:5]; v[0] = 999
assert a[2] == 999 and np.shares_memory(a, v)
c = a[[2, 3, 4]]; c[0] = -1
assert a[2] == 999 and not np.shares_memory(a, c)
mask_copy = a[a > 5]
a[a > 5] = 0
assert 999 in mask_copy and 999 not in a
counts = np.zeros(3, dtype=int)
np.add.at(counts, [0, 0, 2], 1)
assert np.array_equal(counts, [2, 0, 1])
print('slice shapes, views, gathers, mask assignment and repeated updates passed')

4 · Axes — what does axis=0 actually do?

For a reduction such as sum, the named axes are the ones combined. This is not a rule for every operation taking an axis argument: sorting along an axis does not remove it, and softmax retains the input shape.

4.1 Worked numeric example

For rows [1,2,3] and [4,5,6], summing over axis 0 combines the two rows into [5,7,9], shape (3,). Summing over axis 1 combines each row's three columns into [6,15], shape (2,). Describe the retained dimensions before naming the result “row” or “column” sums.

4.2 Same rule works for higher dimensions

For shape (4,5,6), a sum over axis 0 yields (5,6), axis 1 yields (4,6), axis 2 yields (4,5), and axes (0,2) yield (5,). Summing all entries gives a scalar result with shape (). Axis -1 refers to the final axis, which is useful only if that axis really has the semantic role you intend.

4.3 keepdims — the flag that saves your neural network

Reducing (B,D) over the feature axis with keepdims=True gives (B,1). That makes subtraction or division by one number per row broadcast as intended. Without it, shape (B,) aligns with the last dimension: it may fail or, if sizes happen to align, compute the wrong operation.

For softmax, first subtract each row's maximum, exponentiate, then divide by that row's exponential sum, retaining singleton dimensions for both reductions. The complete softmax does not collapse the class axis. In exact arithmetic, finite input rows yield strictly positive normalized probabilities. In floating-point arithmetic, very negative shifted logits can underflow to zero: float64 softmax of [0., -1000.] gives [1., 0.]. Subtracting the maximum prevents exponential overflow but does not guarantee strict positivity or make log(softmax(x)) numerically safe; compute log probabilities with a stable log-softmax/log-sum-exp formulation instead. An entirely masked row of negative infinities needs a separate policy because subtracting its maximum is undefined.

5 · Reshaping — same data, different shape

np.arange(12).reshape(3,4) can share storage. .reshape(-1,4) infers the missing dimension from the total count; only one unknown dimension can be inferred. For general layouts, reshape may copy. ravel produces a contiguous one-dimensional view when possible and otherwise copies; flatten always copies.

5.1 Transpose

For a NumPy ndarray, .T reverses all axes, not merely the last two. Thus a (2,3,4) array becomes (4,3,2). Use swapaxes(-2,-1) to swap only the matrix axes in a batch, or give the complete permutation to transpose. A 1-D transpose is still 1-D.

5.2 Adding and removing axes

For v.shape==(3,), v[:,None] is (3,1), v[None,:] is (1,3), and v[None,:,None] is (1,3,1). Use squeeze(axis=...) when only a particular singleton axis should disappear. An unrestricted squeeze can remove the batch axis when the batch happens to contain one example.

import numpy as np
 
X = np.array([[1., 2., 3.], [4., 5., 6.]])
assert np.array_equal(X.sum(axis=0), [5, 7, 9])
assert np.array_equal(X.sum(axis=1), [6, 15])
T = np.arange(120).reshape(4, 5, 6)
assert T.sum(axis=(0, 2)).shape == (5,)
assert T.sum(axis=1, keepdims=True).shape == (4, 1, 6)
assert T.T.shape == (6, 5, 4)
assert T.swapaxes(-2, -1).shape == (4, 6, 5)
v = np.arange(3)
assert v[:, None].shape == (3, 1) and v[None, :].shape == (1, 3)
batch = np.ones((1, 4, 1))
assert batch.squeeze(axis=-1).shape == (1, 4)
assert batch.squeeze().shape == (4,)
a = np.arange(12).reshape(3, 4)
assert np.shares_memory(a, a.ravel())
assert not np.shares_memory(a, a.T.ravel())
assert not np.shares_memory(a, a.flatten())
# A contiguous slice can still retain its parent through ascontiguousarray.
small = a[0:1]
assert np.shares_memory(a, np.ascontiguousarray(small))
assert not np.shares_memory(a, small.copy())
logits = np.array([[1000., 1001., 1002.], [-1000., -1000., -1000.]])
exp = np.exp(logits - logits.max(axis=-1, keepdims=True))
probs = exp / exp.sum(axis=-1, keepdims=True)
assert probs.shape == logits.shape and np.allclose(probs.sum(axis=-1), 1)
print('axes, transpose, ownership and stable-softmax checks passed')

6 · Element-wise operations and vectorization

Vectorization moves repeated numerical work into library operations. It is not a ban on Python loops: blockwise loops can prevent a huge temporary array, and small computations may be dominated by call/allocation overhead. Compare equal computations on equal data before comparing timing.

6.1 The slow way

A Python column-sum implementation visits each element and updates a Python accumulator. That makes the semantics easy to inspect and provides a reference for the bulk operation.

6.2 The fast way

X.sum(axis=0) performs the same mathematical reduction through NumPy. Floating-point addition order may differ, so compare with tolerance. The bounded timing experiment below prints measurements from this run rather than preprinting an expected speedup. Repetition, dtype, shape and thread settings affect results; these small timings do not extrapolate to GPUs or full training workloads.

6.3 The vocabulary

a+b, a*b, a**2, np.exp(a), np.log(a) and np.maximum(a,0) act elementwise with applicable broadcasting. a@b contracts matrix dimensions, broadcasting leading batch dimensions. For 1-D operands, matmul has special vector rules: two vectors give a scalar dot product. np.dot agrees for ordinary 1-D/2-D examples but differs from matmul for higher-dimensional inputs. A log requires a suitable positive real input if a finite real output is intended; a shape-correct calculation can still be numerically invalid.

7 · Reductions and stats — the daily-driver ops

sum, mean, min, max, var, std, argmin, argmax and cumsum answer different questions. Argmax returns a position, not the value or a calibrated confidence. On (batch,classes) logits, argmax(axis=-1) gives one class index per row; the default flattened argmax does not.

NumPy variance uses divisor N by default (ddof=0); ddof=1 uses N−1 and requires enough observations. That distinction does not by itself solve dependence or make standard deviation unbiased. Float32 accumulation can lose small contributions; sum(dtype=np.float64) can help, with extra cost. Fixed-width integers can overflow. The earlier personal performance anecdote is replaced by observable arithmetic rather than an unverifiable claimed accuracy loss.

import numpy as np
import timeit
 
X = np.arange(8000, dtype=np.float64).reshape(400, 20)
def column_loop():
    out = np.zeros(X.shape[1])
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            out[j] += X[i, j]
    return out
 
def column_bulk():
    return X.sum(axis=0)
 
np.testing.assert_allclose(column_loop(), column_bulk())
column_loop(); column_bulk()  # warm up both
loop_s = min(timeit.repeat(column_loop, number=3, repeat=3)) / 3
bulk_s = min(timeit.repeat(column_bulk, number=3, repeat=3)) / 3
print('numpy', np.__version__, 'shape', X.shape, 'dtype', X.dtype)
print('seconds/call, loop and bulk:', loop_s, bulk_s)
print('observed local ratio:', loop_s / bulk_s)
a = np.array([1e8, 1., -1e8], dtype=np.float32)
assert a.sum(dtype=np.float64) == 1
assert np.array_equal(np.array([[1, 3], [4, 2]]).argmax(axis=-1), [1, 0])
assert np.isclose(np.array([1., 3.]).var(), 1)
assert np.isclose(np.array([1., 3.]).var(ddof=1), 2)

8 · A tiny worked example — normalize a batch of vectors

For row (3,4), the norm is √(9+16)=5, so the normalized row is (0.6,0.8). For five rows of width three, the norm array must have shape (5,1) to divide each row by its own norm. The zero row has no direction: declare whether it is rejected, masked or kept at zero. Keeping it at zero is a policy, not a unit normalization.

The following also solves the original cosine-similarity stretch exercise on a bounded 16-row batch. Center columns first, normalize the resulting rows, multiply by the transpose, and exclude the query itself before selecting its neighbours. A zero-row policy and a deterministic tie order make the output interpretable.

import numpy as np
 
V = np.array([[3., 4.], [0., 0.], [4., -3.]])
norms = np.linalg.norm(V, axis=1, keepdims=True)
U = np.divide(V, norms, out=np.zeros_like(V), where=norms != 0)
assert np.allclose(U[0], [.6, .8])
assert np.all(U[1] == 0)
assert np.allclose(np.linalg.norm(U[[0, 2]], axis=1), 1)
 
X = np.random.default_rng(1).normal(size=(16, 5))
centered = X - X.mean(axis=0, keepdims=True)
assert np.allclose(centered.mean(axis=0), 0)
norms = np.linalg.norm(centered, axis=1, keepdims=True)
assert np.all(norms > 0)  # this particular generated batch
unit = centered / norms
similarity = unit @ unit.T
assert similarity.shape == (16, 16)
assert np.allclose(similarity, similarity.T)
assert np.allclose(np.diag(similarity), 1)
scores = similarity[0].copy(); scores[0] = -np.inf
order = np.lexsort((np.arange(len(scores)), -scores))
top5 = order[:5]
assert 0 not in top5 and len(set(top5)) == 5
print('top-five neighbours of row zero:', top5.tolist())

For 1,000 float64 rows, a full (1000,1000) score matrix uses eight million data bytes, excluding inputs and intermediates. If only neighbours of row zero are needed, compute unit @ unit[0] instead. For much larger collections, blocking or a retrieval index avoids materializing every pair. Centering is not automatically beneficial for a pretrained embedding space; this exercise teaches operations, not a universal retrieval recipe.

8b · From array operations to model mechanisms

Shapes are necessary evidence that an implementation is well formed, not proof that you understand its objective, learned representation or numerical behaviour. The following previews repair the original field report without asserting that architecture research is merely a broadcasting trick.

8b.1 · The attention core, in NumPy shapes

For q, k and v of shape (B,H,T,D), q @ k.swapaxes(-2,-1) gives (B,H,T,T). After scaling by √D and applying a causal mask, normalize over keys (the last axis); multiplying weights by v returns (B,H,T,D). NumPy has no ndarray masked_fill method, and NumPy's transpose(-2,-1) is not PyTorch's two-axis swap API for a 4-D tensor.

import numpy as np
rng = np.random.default_rng(11)
B, H, T, D = 1, 2, 3, 4
q, k, v = (rng.normal(size=(B, H, T, D)) for _ in range(3))
scores = (q @ k.swapaxes(-2, -1)) / np.sqrt(D)
allowed = np.tril(np.ones((T, T), dtype=bool))
scores = np.where(allowed, scores, -np.inf)
weights = np.exp(scores - scores.max(axis=-1, keepdims=True))
weights /= weights.sum(axis=-1, keepdims=True)
out = weights @ v
assert out.shape == (B, H, T, D)
assert np.allclose(weights.sum(axis=-1), 1)
assert np.all(weights[..., ~allowed] == 0)
assert np.allclose(out[:, :, 0], v[:, :, 0])
print('causal attention shapes and first-position identity passed')

This tiny example explicitly materializes the score matrix. It does not implement a fused attention kernel or test its speed.

8b.2 · FlashAttention: logical matrix versus materialized matrix

The original FlashAttention and FlashAttention-2 already avoid storing the full quadratic attention matrix in high-bandwidth memory. It was incorrect to present that as a new difference between FA-2 and FA-3. FA-2's paper emphasizes work partitioning, reducing non-matmul operations and improving occupancy/communication. FA-3's paper emphasizes Hopper asynchrony, overlapping computation/data movement, interleaving blockwise matmul/softmax, and FP8-related techniques.

The array expression specifies a computation; a tiled implementation can combine partial sums and normalization statistics without allocating every logical intermediate. Floating-point implementations need not be bit-identical. Paper-reported utilization and speedups belong to their specified hardware and benchmark configurations, not to this CPU notebook. No H100 benchmark is reproduced or promised here.

8b.3 · Latent KV compression: learned maps, not just reshape

The DeepSeek-V3 technical report describes MLA as a learned low-rank joint compression of attention keys and values. A down-projection maps each token representation to a smaller latent; up-projections produce content keys and values. Reshape only organizes those projected coordinates into heads—it does not itself learn or perform the compression.

For a deliberately simplified comparison, conventional K and V storage contains $2BHTD$ values, while a latent-only store contains $BTd_c$. With B=1, T=8192, H=128, D=128 and $d_c=512$, their ratio is 64, not 14. At two bytes per value these are 512 MiB versus 8 MiB per layer. This is toy arithmetic, not the complete DeepSeek cache: the report explicitly caches both the compressed KV latent and the decoupled RoPE key. Actual head dimensions, positional terms, dtype and runtime buffers alter the total.

Additional projections or a changed contraction order also affect compute. Fewer stored values does not automatically mean proportionally less latency. The old multi-million-dollar headline is omitted rather than treated as total research cost; the report's own cost table excludes prior research and ablations.

8b.4 · Positional angle grids and long-context limits

For T positions and D/2 rotary frequencies, positions[:,None] * frequencies[None,:] creates a (T,D/2) angle grid. For positions0,1,2 and frequencies1,.1, rows are(0,0),(1,.1),(2,.2). Cosine and sine define rotations of coordinate pairs. Multiplying a frequency by1/s changes the phase at a fixed position to the original phase at position t/s; this is position interpolation, not a proof of useful long-context reasoning.

Model identity matters. Meta's official Llama4 card, released2025-04-05, lists Scout10M and Maverick1M context, not a10M limit for every Llama4. Google's Gemma3 card lists128K total context for4B/12B/27B and32K for1B (and the later270M variant in the accessed card). Input and output share the context budget; a headline length is not an independent output allowance. These are dated model-card specifications, not tested recall performance.

YaRN and LongRoPE retain genuine technical substance beyond the grid. YaRN is a RoPE extension method using interpolation/scaling choices and attention rescaling to improve the tested extension regime. LongRoPE searches nonuniform positional interpolation, extends progressively through an intermediate long-context model and readjusts for short contexts; its original paper reports2048K extension in specified LLaMA2/Mistral experiments. Neither paper alone proves that Meta or Google used that method to achieve its own announced limit. A single constant multiplication of every frequency is not a full implementation of either method.

The shape calculation exposes costs: a dense angle table uses T(D/2) entries; its sine and cosine tables add storage if materialized. Dense attention scores grow as T², while ordinary cached keys/values grow linearly with T. At one million positions and head width128, one float32 angle table costs256,000,000 bytes; sine and cosine together add512,000,000 if separately retained. That is just positional tables, not model memory. Algorithms may generate slices on demand. Position encoding, cache size, attention implementation and measured retrieval accuracy are separate engineering questions.

8b.5 · Quantization — dtype fluency pays off

For one million stored values, float64 uses8MB decimal and float32 uses4MB. Ideal packed four-bit codes use0.5MB before scales, zero points, padding, metadata and runtime buffers. Reducing weight storage need not halve total inference memory because activations, caches, workspaces and unquantized tensors remain.

A symmetric quantizer can choose scale s=max(abs(w))/7 for signed codes−7 through7, then q=clip(round(w/s),-7,7) and reconstructed weights s*q. For w=(−1,−.3,.2,1), s=1/7 gives q=(−7,−2,1,7), reconstructed(−1,−2/7,1/7,1). Away from clipping, nearest rounding bounds absolute error by s/2. A fixed scale can clip outliers and exceed that bound. Asymmetric formats add a zero point; per-channel/group scales change both the error and metadata cost. This toy stores codes in int8 for inspection—it does not pack nibbles or implement a fast integer matrix kernel.

The original named formats describe different mechanisms:

• FP8 E4M3 and E5M2: the primary format paper distinguishes four-exponent/three-mantissa and five-exponent/two-mantissa encodings, with a sign bit and different special-value/range conventions. They are floating-point formats, not INT8. Their tradeoff concerns range and precision; scaling and accumulation remain important. A weight payload going from two bytes to one halves that payload, not necessarily runtime or all memory.
• GPTQ: one-shot post-training weight quantization uses approximate second-order information to manage quantization error. It is not a NumPy dtype cast.
• AWQ: activation statistics identify salient channels and guide equivalent channel scaling before low-bit weight-only quantization. The source does not make every channel equally important or require preserving a small subset as a separate high-precision array.
• QuIP#: randomized Hadamard incoherence processing, lattice-based vector codebooks and fine-tuning are distinct ingredients. Calling it merely four-bit integers loses the algorithm.
• BitNet b1.58: Ma and colleagues'2024 paper, The Era of 1-bit LLMs: All Large Language Models are in1.58 Bits, studies ternary weights in a training recipe. It is not the older BitNet: Scaling1-bit Transformers paper, nor evidence that casting any pretrained model to−1,0,1 preserves quality. log₂3≈1.58496 is the alphabet information count; actual packed representation, scales and other tensors determine bytes. Reported quality equivalence is confined to the paper's matched model/token experiments.

Casting floats directly to np.int8 truncates fractional values and supplies no calibration scale. Distinguish storage dtype, reconstructed values, compute dtype and accumulation dtype before trusting a memory or accuracy claim.

import numpy as np
 
positions=np.arange(3);frequencies=np.array([1.,.1])
angles=positions[:,None]*frequencies[None,:]
np.testing.assert_allclose(angles,[[0,0],[1,.1],[2,.2]])
np.testing.assert_allclose(positions[:,None]*(frequencies/4)[None,:],(positions/4)[:,None]*frequencies)
assert 1_000_000*(128//2)*4==256_000_000
w=np.array([-1.,-.3,.2,1.]);scale=np.max(np.abs(w))/7
q=np.clip(np.rint(w/scale),-7,7).astype(np.int8);reconstructed=q*scale
np.testing.assert_array_equal(q,[-7,-2,1,7])
assert np.max(np.abs(w-reconstructed)) <= scale/2+1e-12
assert not np.allclose(w,w.astype(np.int8))
fixed_scale=.1;outlier=3.
clipped=np.clip(np.rint(outlier/fixed_scale),-7,7)*fixed_scale
assert abs(outlier-clipped)>fixed_scale/2
zero=np.zeros(4);zero_scale=np.max(np.abs(zero))/7
# A zero group requires an explicit convention rather than dividing by zero.
qzero=np.zeros_like(zero,dtype=np.int8) if zero_scale==0 else np.rint(zero/zero_scale).astype(np.int8)
assert np.all(qzero==0)
assert np.isclose(np.log2(3),1.584962500721156)
assert q.nbytes==4  # int8 teaching codes, not a two-byte packed INT4 array
print('rotary grid, interpolation, positional-table bytes and quantization error controls passed')

9 · Reproducible bugs and ownership tradeoffs

• All-pairs metric: predictions (B,) compared to labels (B,1) produce (B,B). Assert equal intended shapes; use reshape(-1) only when one scalar label per row is truly the contract.
• Mutated dataset: a basic crop is a view. For an image of shape (3,224,224), img[:,16:240,16:240] is (3,208,208), not (3,224,224): slice bounds clip to the actual extent. Copy before an augmentation that must not alter its source.
• Wrong multiplication: with a square X, X.T * X can produce a plausible-shaped elementwise product. A shape check alone will not detect that the intended Gram matrix was X.T @ X; use a numeric oracle too.
• Wrong softmax axis: normalize over classes, not observations. The final result retains all axes; only its intermediate maxima/sums reduce.
• Unexpected precision: test the actual storage and reduction dtype. Neither “always float64” nor “always lowest precision” is a substitute for an error budget.

Passing a writable view is useful for an explicitly shared output buffer. Copying isolates mutation and can release a small result from a large owner's lifetime. Setting view.flags.writeable=False blocks writes through that view but does not stop another alias mutating the same memory. And np.ascontiguousarray is not guaranteed to break ownership: if the input is already contiguous it can return it unchanged. Use .copy() when independence is required; the tested example above demonstrates the difference.

For (batch, sequence, embedding), each row is one token's feature vector and each batch item contains a sequence of such rows. Reducing the embedding axis combines features within a token; reducing the sequence axis combines token positions. A reshape never supplies those meanings by itself. Draw the axis roles, then separately ask whether mutation should propagate through a shared view.

11 · Retention scaffold

Recall — answer from memory (no scrolling)

1. What are the shapes and data bytes of a (2,3) float32 array and a (3,) int64 array?
2. For x.shape==(2,3,4), predict x[:,1], x[...,2] and x.T shapes.
3. Why does a (B,) row-norm array not generally divide a (B,D) matrix correctly?
4. Which creates a copy: a basic slice, an integer-array gather, or flatten?
5. Why can base is original and ascontiguousarray both mislead an ownership check?
6. What is the latent-only toy cache ratio above, and which real MLA cache component did it omit?

Original warm-up questions — worked answers

For a (3,4) X, sum(axis=0) combines the three entries in each column and returns four values. For (5,3) embeddings, retained norms are (5,1), so each norm is reused over exactly three features. An unreshaped (5,) divisor aligns with the feature axis and fails here; square examples can hide the same semantic error.

The two characteristic silent bugs are an unintended all-pairs metric and a mutation leaking through a view. Check the residual/comparison shape before reducing; for ownership, test sharing and define whether the callee may write. Copy when independence or release of a large retained parent matters. A read-only view guards only writes through that alias.

For (10000,50) column sums, both the Python loop and NumPy sum process500,000 values; the library avoids per-element Python dispatch and may use optimized iteration. The actual ratio depends on dtype/layout/hardware/threading, so the bounded timing program reports measurements instead of a guessed25×/100×. Comparing unequal computations or timing one cold call is not a defensible benchmark.

The original MLA exercise's latent-only ratio is64, with512MiB conventional KV versus8MiB latent under its two-byte toy assumptions. Extra positional cache terms change the real count; up-projection or an altered contraction can also add compute. Smaller storage does not entail a64-fold speedup.

Stretch — for one extra hour of practice

Rebuild the centering, normalization and top-five-neighbour exercise without looking. Test a zero vector, duplicate vectors, one-row input and a non-contiguous slice. For one row, there are no other neighbours: do not return the query as its own answer. For duplicates, specify a tie policy. Reduce the full similarity-matrix implementation to a single-query vector and verify identical rankings on non-tied inputs.

In your own words

Explain how dtype, shape and strides answer different questions. Then explain why reducing Python-loop overhead is a mechanism, not a universal speed multiplier.

Spaced-review pointer

Tomorrow, predict the five slice shapes before running the indexing test. A week later, recreate row normalization including its zero policy and the repeated-index add.at test. Use those checks to distinguish retention from merely recognizing code.

Next-session teaser

DL002 gives these arrays geometric meaning: dot products combine direction and magnitude, matrices transform basis vectors, and an all-pairs score table contracts a feature axis. Keep the zero-vector and shape exceptions from this lesson when turning those calculations into similarities.

What to bring back tomorrow — sticky note

Write the axis roles beside every array. Check the intended reduction axis and retain dimensions when broadcasting back. Ask whether mutation should propagate. Preserve a tiny numeric oracle and test a batch larger than one.

Stretch edge-case solutions

A zero centered row has no cosine direction; reject it or explicitly exclude it rather than claiming it became a unit vector. With a one-row batch, column centering makes that row zero and there are no non-self neighbors. Duplicate normalized rows can tie at1: use row ID as a deterministic secondary key and exclude self. A noncontiguous input should preserve numerical results but may alter cost.

For the original1000×50 exercise, C=X-X.mean(axis=0,keepdims=True), U=C/norm(C,axis=1,keepdims=True) for nonzero rows, and S=U@U.T solve centering and pairwise cosine. For only row0's neighbors, U@U[0] avoids the full matrix; sort by descending score then ascending ID, discard self, and return at most the number of eligible rows. The original16-row implementation demonstrates exactly this mechanism at a safe size.

Sources and tested boundary

• NumPy 2.2 copies and views: indexing, reshape, ravel and flatten contracts.
• NumPy 2.2 Generator and fromiter: local randomness, compatibility limits and generator consumption.
• PyTorch 2.9 tensor views: documentation-supported view/copy preview; no PyTorch execution.
• FlashAttention-2, v1 and FlashAttention-3, v2: distinct kernel contributions, not CPU-measured speed claims.
• DeepSeek-V3 technical report, v2, §2.1.1: learned latent projections and the additional cached positional key.

The small NumPy programs and the two previously protected programs are tested locally. Timing output is specific to this run and environment. No large language model, GPU kernel, PyTorch runtime, model quality or cloud service is reproduced by these examples.

Additional primary context: Llama4 official model card, Gemma3 model card, YaRN, LongRoPE, FP8 formats, GPTQ, AWQ, QuIP#, and BitNet b1.58. Context specifications and paper results are source-supported comparisons, not downloads or runtime reproductions of those models.

Next: Vectors, matrices and geometry