00 · The jargon

Key concepts before we start

These terms come up over and over. Here’s what each one means, visually.

The starter toolbox

These tiny functions appear in every section from here on. See them once now, then watch them reappear.

import numpy as np

# Softmax: turn any raw scores into percentages that sum to 100%
def softmax(x):
    ex = np.exp(x - np.max(x))
    return ex / ex.sum()

# LayerNorm: center and scale values to keep them stable
def layer_norm(x, eps=1e-5):
    return (x - x.mean()) / (x.std() + eps)

# ReLU: simplest activation — zero out negatives, keep positives
def relu(x):
    return np.maximum(0, x)

# Logit → probability: the raw-to-percentage pipeline
logits = np.array([3.2, 1.1, 0.5])
probs = softmax(logits)
print(f"logits: {logits} → probs: {probs.round(3)}")
# Output: logits: [3.2 1.1 0.5] → probs: [0.761 0.254 0.146]

Notice how softmax and layer_norm are each just 3 lines. They are the only “magic” you need to understand everything that follows.

01 · The foundation

Every neural network, at its core: multiply, add, repeat

An LLM has billions of numbers inside it. But those billions are just the same tiny machine — the neuron — repeated over and over. Understand the neuron, and you understand the foundation of everything.

One neuron = one weighted vote

A single neuron does one simple calculation:

Think of it like voting:

• Each input is a person casting a vote
• Each weight (w) is how much that person’s vote counts
• The bias is a fixed advantage one candidate gets regardless of the vote

So: the neuron multiplies each input by its weight, adds them all up, adds the bias, and spits out one number. That’s it. This operation is called a linear layer — but that’s just a fancy name for “multiply each input by its importance, add them up.”

# One neuron = multiply inputs by weights, add bias
inputs = [2.0, 4.0, 1.5]
weights = [0.8, -0.3, 1.2]
bias = 0.5

output = 0
for i in range(len(inputs)):
    output += inputs[i] * weights[i]
output += bias

print(f"output = {output:.2f}")  # output = 3.10

Now imagine many neurons side by side, all looking at the same inputs but each with its own weights and bias — that’s a layer. Writing it as a group of equations is the same as a matrix multiplication.

Why one layer isn’t enough

Here’s the trap: if you stack linear layers on top of each other, $W_2\cdot (W_1\cdot x+b_1)+b_2$ just simplifies to $W^{\prime }\cdot x+b^{\prime }$ — one big linear equation. No matter how many layers you stack, you’re still just drawing straight lines. The depth is useless.

The fix: after each linear layer, apply an activation function — a simple “bend” that breaks the straight line. The most common ones:

• ReLU: if the output is negative, make it zero. Positive stays as-is. Like a light switch — either on or off.
• GELU: a smooth version of ReLU (used in GPT-2, BERT)
• SwiGLU: has two paths — one that creates information and another that acts like a gate, deciding what passes through (used in LLaMA, most modern LLMs)

Without these bends, a 100-layer network is no smarter than a 1-layer network. With them, each layer can build on the previous one in increasingly abstract ways — that’s where the depth earns its power.

import math

def relu(x):      return max(0, x)
def gelu(x):      return 0.5 * x * (1 + math.erf(x / math.sqrt(2)))
def swish(x):     return x / (1 + math.exp(-x))

x = [0.5, -1.2, 0.8]
print(f"ReLU:  {[relu(v) for v in x]}")
print(f"GELU:  {[round(gelu(v), 2) for v in x]}")
print(f"Swish: {[round(swish(v), 2) for v in x]}")

Now here’s the key insight: if you stack multiple neurons side by side — all looking at the same inputs but each with its own weights — that’s a layer. And writing that layer as a group of equations is the same as a matrix multiplication.

Where $w_{ij}$ means “weight from input j to neuron i.” The row index = target neuron, column index = source input.

# One layer = matrix multiply (many neurons at once)
import numpy as np

X = np.array([2.0, 4.0])                          # inputs
W = np.array([[0.8, -0.3],                        # weights: 3 neurons × 2 inputs
              [1.2,  0.5],
              [-0.4, 0.9]])
b = np.array([0.5, -0.2, 0.1])                    # biases

output = X @ W.T + b                               # @ = matrix multiply
print(f"layer output: {output}")                   # [ 2.1 ,  2.7 ,  0.9 ]

The @ is matrix multiplication. The .T transposes W (flips rows ↔ columns) so the dimensions line up. This one line replaces writing out every neuron’s equation by hand.

Every core computation in an LLM — embedding lookup, attention projection, MLP layers, output prediction — is built from $W\cdot x+b$ operations (with activation functions breaking linearity at key points). Keep that in your head as we go.

02 · The data problem

How to eat the internet

Computers can’t read words — they need numbers. So how do we turn “the cat sat on the mat” into numbers without losing meaning?

Simple idea: give every word a number. “the” = 1, “cat” = 2, “sat” = 3. This fails because:

• “run” and “running” get different numbers even though they’re related
• There are too many words (and new ones created daily)
• Misspellings and names would need entirely new numbers

The real solution is BPE (Byte Pair Encoding): start with individual letters, then merge the most frequent letter pairs into new tokens. Repeat until you have a useful vocabulary.

# BPE tokenizer: one merge step
corpus = "the cat and the dog sat on the mat and watched the cat run"

# Step 1: start with individual characters
tokens = list(corpus)
print(f"Characters: {''.join(tokens)}")

# Step 2: find the most frequent adjacent pair
def most_frequent_pair(tokens):
    pairs = {}
    for i in range(len(tokens) - 1):
        pair = (tokens[i], tokens[i+1])
        pairs[pair] = pairs.get(pair, 0) + 1
    return max(pairs, key=pairs.get)

# Step 3: merge that pair everywhere
pair = most_frequent_pair(tokens)
new_tokens, i = [], 0
while i < len(tokens):
    if i < len(tokens) - 1 and (tokens[i], tokens[i+1]) == pair:
        new_tokens.append(tokens[i] + tokens[i+1])
        i += 2
    else:
        new_tokens.append(tokens[i])
        i += 1

print(f"After merging '{pair[0]}' + '{pair[1]}': {''.join(new_tokens)}")
print(f"Vocabulary: {sorted(set(new_tokens))}")

The merge rules are learned once, before training ever starts. After that, every input text is tokenized the same way.

Once BPE is done, every piece of text becomes a sequence of token IDs. “the cat sat” becomes something like [5, 182, 93, 12, 5, 44]. The model’s entire world is these integers.

03 · Embeddings

Turning token IDs into meaningful vectors

We have token IDs: [5, 182, 93, …]. But the number 182 means nothing — it’s just an arbitrary label. “cat” = 182 tells the model nothing about what a cat is.

We need each token to become a vector — a list of numbers — where those numbers actually encode meaning. Similar tokens should have similar vectors. This vector is called an embedding.

Think of it like a city map:

• Each token gets coordinates (its vector)
• Tokens used in similar ways end up in similar neighborhoods
• “cat” and “dog” live near each other; “math” and “poetry” are in different districts

The embeddings are stored in a giant table. Looking up “cat” gives you its vector. This is just a special case of the $W\cdot x+b$ we learned about: you’re multiplying a one-hot vector (all zeros except one 1) by the embedding matrix. The matrix row you pick is the embedding.

A neat trick many LLMs use: the same table is used both to turn token IDs into vectors (at the input) AND to turn vectors back into token predictions (at the output). The input mapping and output prediction share the same “meaning” space — this is called tied embeddings and cuts the parameter count roughly in half.

import numpy as np

vocab_size = 6  # tiny: ["the", "cat", "dog", "sat", "mat", "run"]
embed_dim = 4   # 4 numbers per token (real LLMs use 768–4096)

# Embedding table: one row per token, one column per dimension
E = np.random.randn(vocab_size, embed_dim) * 0.1

# Look up "cat" (index 1) — just grab its row
cat_vector = E[1]
print(f"cat's embedding vector: {cat_vector}")

# Same lookup as one-hot × matrix multiplication
one_hot = np.array([0, 1, 0, 0, 0, 0])
cat_vector2 = one_hot @ E  # same result
print(f"same via one-hot: {cat_vector2}")

Before training, all vectors are random. During training, the model adjusts them: words used in similar ways drift toward similar vectors. “cat” and “dog” end up nearby because they appear before “sat” and after “the.” The model never sees a dictionary — meaning emerges from usage patterns alone.

04 · Positional encoding

How the model knows which word came first

Think about the sentence “dog bit man” vs “man bit dog.” Same three words. Completely different meaning. The model needs to know which word came first — order matters.

But the Transformer reads all tokens at once (not one by one). It has zero built-in sense of order. So we must manually inject position info.

The simplest idea: just add 1 to token 1’s vector, add 2 to token 2’s vector, etc. Why doesn’t this work? Read on.

What actually happens, step by step

We directly modify the token’s vector. There is no separate “position input.” Here’s exactly what happens:

1. Token “dog” at position 3 has an embedding vector, like [0.2, -0.5, 0.8, ...] (a list of 4096 numbers)
2. We compute a position vector for position 3 using sine/cosine waves, giving us another vector of the same length, like [0.14, -0.99, -0.76, ...]
3. We add the two vectors together: new_vector = embedding_vector + position_vector
4. This summed vector is what the model works with from that point on

The position info is mixed into the token’s numbers. The model doesn’t know “here are 50 numbers for meaning and 50 numbers for position” — they’re all blended together, and the model learns to use whichever parts it needs.

The solution: sine + cosine waves (together they make a fingerprint)

Imagine you’re writing a long exam and every student needs a unique ID. If you give them just one number (1, 2, 3…), you run out of space or the numbers get huge. Instead, imagine giving each student two numbers: one from a fast-ticking clock and one from a slow-ticking clock. Together they form a unique pair.

Similarly, for each position in the sequence, we generate a unique fingerprint using two complementary waves:

• One wave uses sine (sin) — think of it as the “vertical coordinate”
• The other uses cosine (cos) — think of it as the “horizontal coordinate”

Why both? Because one alone isn’t enough. Just like a GPS needs both latitude AND longitude to pinpoint a location, sin gives you one value and cos gives you another. Together they form a unique pair for every position. If you used only sin, many different positions would have the same value (sin repeats every 360°). Using sin+cos as a pair creates a unique signature.

Each dimension of the position vector uses a wave at a different frequency:

• Dimension 0: very fast wave (changes every 1-2 positions) — captures local relationships like “the” before “cat”
• Dimension 1: slightly slower wave — captures slightly longer patterns
• Dimension 10: slow wave (changes every 50+ positions) — captures long-range patterns like a pronoun referring to a noun 30 words back
• Dimension 100+: extremely slow wave — barely changes across the whole sequence

All these waves oscillate between -1 and +1 forever. Position 1,000,000 gets a valid fingerprint just as well as position 1.

The formula, broken down simply

Don’t be scared. Here’s what each symbol means:

• p = position number (0, 1, 2, 3… which word in the sequence)
• i = which dimension of the vector we’re computing (0, 1, 2, 3… up to 4096)
• d = total number of dimensions in the vector (e.g., 4096)
• 2i and 2i+1 = even and odd dimension numbers (sin for even, cos for odd)

Now the inside: $\frac{p}{10000^{2i/d}}$

Think of it like this:

• 10000 is just a big reference number. It’s like saying “the slowest wave takes 10,000 positions to complete one full cycle”
• 2i/d is the fraction of “how far through the dimensions we are.” For dimension 0, it’s 0/4096 = 0. For dimension 256, it’s 512/4096 ≈ 0.125. For the last dimension, it’s 4096/4096 = 1.
• $10000^{2i/d}$ for dimension 0 gives $10000^0=1$, so the wave speed is p/1 = p — this is the fastest wave (changes with every position)
• For the last dimension, $10000^1=10000$, so the wave speed is p/10000 — this is the slowest wave (barely changes)

Different dimensions get different denominators: 1, 10000¹ᐟ⁴⁰⁹⁶, 10000²ᐟ⁴⁰⁹⁶, … up to 10000. This creates a smooth range from “fast wave” to “slow wave.”

If the formula still feels fuzzy, just remember: each position gets a unique set of numbers from -1 to +1, and that set is added directly to the token’s embedding vector. The model learns to use those numbers to tell where each word sits in the sentence.

import math

def positional_encoding(seq_len, d_model):
    PE = []
    for pos in range(seq_len):
        pos_vec = []
        for i in range(d_model):
            denom = 10000 ** ((2 * (i // 2)) / d_model)
            if i % 2 == 0:
                pos_vec.append(math.sin(pos / denom))
            else:
                pos_vec.append(math.cos(pos / denom))
        PE.append(pos_vec)
    return PE

seq_len, d_model = 50, 6
PE = positional_encoding(seq_len, d_model)

for pos in range(8):
    print(f"pos {pos}: {[round(v, 2) for v in PE[pos]]}")

# Adding position to embedding (this is what the model actually does):
embedding_cat = [0.2, -0.5, 0.8, 0.1, -0.3, 0.6]
final_vector = [embedding_cat[i] + PE[3][i] for i in range(d_model)]
print(f"\nembedding + pos=3: {[round(v, 2) for v in final_vector]}")

Each position gets a unique fingerprint. Position 0 is different from position 5, and position 1000 works just as well as position 1. The vector is added directly to the token’s embedding — they’re blended together.

05 · Attention

The core: letting tokens look at each other

Now every token has meaning + position. But it still doesn’t know which other tokens are relevant to it.

Attention is the mechanism that answers: “Which other words in this sentence should I look at to understand myself better?”

Think of it like a classroom discussion:

• Each student (token) can hear all other students
• But not everyone is equally relevant to what you’re saying
• You need to figure out who to listen to and how much

Why attention replaced the old way (RNNs)

Before Transformers, models called RNNs processed words one at a time — word 1 finishes before word 2 starts. Like reading a book one letter at a time through a tiny straw. Worse, information from word 1 had to travel through every intermediate step to reach word 100 — it degraded along the way (like the telephone game).

Attention changed everything: it compares all words at once, directly. Word 1 connects to word 100 in a single operation, with no degradation. This is why Transformers train faster and handle longer text.

How it works: Q, K, V (like a library)

Every word produces three special vectors from its embedding (using our familiar $W\cdot x+b$):

Think of a library:

• Query (Q): what I am searching for. Like typing in the search bar.
• Key (K): what I contain. Like the title and tags on a book.
• Value (V): the actual information I carry. Like the book’s contents.

Each word asks a question (Query), and every other word’s Key answers whether they’re relevant. If there’s a match, the Value (the actual info) is passed along.

import numpy as np

x = np.array([0.2, 0.8, 0.1, 0.3])  # embedding of "dog"

# Three different weight matrices (learned during training)
W_q = np.random.randn(4, 3) * 0.1
W_k = np.random.randn(4, 3) * 0.1
W_v = np.random.randn(4, 3) * 0.1

Q = x @ W_q
K = x @ W_k
V = x @ W_v

print(f"Query:  {Q.round(2)}")
print(f"Key:    {K.round(2)}")
print(f"Value:  {V.round(2)}")

Same input x, three different projections. Each word does this. Then attention computes how well each word’s Query matches every other word’s Key.

Why we scale: keeping the distribution fair

The raw Query·Key scores get large when vectors are long. Imagine multiplying big numbers together — the results can be huge. When we then apply softmax (which converts raw scores into percentages that add up to 100%), huge scores make softmax too confident — it picks one word and ignores everything else.

The fix: divide all scores by $\sqrt{d_k}$ (the square root of the dimension size). This keeps scores in a reasonable range so softmax produces a smooth distribution.

def softmax(x):
    ex = np.exp(x - np.max(x))       # subtract max for numerical stability
    return ex / ex.sum()

def attention(Q, K, V):
    scores = Q @ K.T                  # dot product: how relevant is each key?
    scores = scores / np.sqrt(K.shape[1])  # scale so softmax isn't too extreme
    weights = softmax(scores)         # convert to percentages
    return weights @ V                # weighted blend of values

# Toy sentence: "the dog chased" → 3 tokens, 4-dim embeddings
X = np.array([[0.2, 0.8, 0.1, 0.3],  # "the"
              [0.9, 0.2, 0.7, 0.1],  # "dog"
              [0.3, 0.6, 0.2, 0.9]]) # "chased"
W_q = np.random.randn(4, 4) * 0.1
W_k = np.random.randn(4, 4) * 0.1
W_v = np.random.randn(4, 4) * 0.1

Q, K, V = X @ W_q, X @ W_k, X @ W_v
output = attention(Q, K, V)

print(f"Attention output:\n{output.round(3)}")
# Each token's new vector = blend of all values, weighted by relevance

The scaling by $\sqrt{d_k}$ prevents one token from hogging all the attention. Without it, softmax becomes a near one-hot vote — the model only “hears” one token and misses everything else.

Multiple heads: many specialists, not one generalist

One type of relationship is not enough. The word “bank” might connect to “river” through one relationship and to “money” through another. So the Transformer runs attention multiple times in parallel — these are called heads. Each head has its own Q, K, V weights, so each can learn a different type of relationship:

• Head 1: “which word is the subject of this verb?”
• Head 2: “what does this pronoun refer to?”
• Head 3: “which adjective describes this noun?”
• etc.

Modern models use 32+ heads per layer. Each generates its own attention distribution, and their outputs are mashed together into one final representation.

# Multi-head attention: multiple specialists in parallel
n_heads = 4
d_k = 3  # per-head dimension

# 4 heads, each with its own Q, K, V projections
heads = []
for h in range(n_heads):
    W_q_h = np.random.randn(4, d_k) * 0.1
    W_k_h = np.random.randn(4, d_k) * 0.1
    W_v_h = np.random.randn(4, d_k) * 0.1
    heads.append((W_q_h, W_k_h, W_v_h))

# Each head computes attention independently
all_head_outputs = []
for h, (W_q, W_k, W_v) in enumerate(heads):
    Q_h, K_h, V_h = X @ W_q, X @ W_k, X @ W_v
    out = attention(Q_h, K_h, V_h)   # same attention() from above
    all_head_outputs.append(out)
    print(f"Head {h} output: {out.round(2)}")

# Concatenate all heads and project back
combined = np.concatenate(all_head_outputs, axis=1)
W_o = np.random.randn(n_heads * d_k, 4) * 0.1
final = combined @ W_o
print(f"\nFinal (all heads combined):\n{final.round(3)}")

Each head runs the same attention mechanism but with its own learned weights. One might learn syntax, another pronouns, another copy behavior — specialization emerges from the training data, not from human labels.

The causal mask: when generating text, the model must predict the next token using only tokens that came before it. If position 5 could look at position 6, it would “cheat” by seeing the answer. A mask sets the attention scores for future positions to $-\infty$, so after softmax they become 0%.

def causal_attention(Q, K, V, mask=None):
    scores = Q @ K.T / np.sqrt(K.shape[1])
    if mask is not None:
        scores = scores + mask  # adding -inf zeros out future positions
    return softmax(scores) @ V

# Mask: upper triangle = -inf (can't see future), lower = 0 (can see past)
seq_len = 4
mask = np.triu(np.full((seq_len, seq_len), -np.inf), k=1)
print(f"Causal mask:\n{mask}")

# With mask: token 3 can see tokens 0,1,2 but NOT token 3's own future
Q = np.random.randn(seq_len, 4)
K = np.random.randn(seq_len, 4)
V = np.random.randn(seq_len, 4)
out = causal_attention(Q, K, V, mask)
print(f"\nWith causal mask:\n{out.round(3)}")
# Without mask: every token sees every other token (for comparison)
out2 = causal_attention(Q, K, V)
print(f"\nWithout mask:\n{out2.round(3)}")

How attention weights are trained

The three projection matrices $W_q,W_k,W_v$ are learned through backpropagation through softmax — the same gradient descent as every other parameter, but with a unique challenge.

During the forward pass, attention computes:

When the loss gradient arrives at the attention output, it must flow backward through three distinct paths:

1. Through $W_v$ (simplest): the gradient tells $W_v$ “the value you produced for token 3 was multiplied by a large attention weight from token 5 — that weight made your contribution to the error large.” $W_v$ adjusts its rows so tokens that should not be attended to produce smaller or less relevant value vectors.

2. Through the attention matrix A (trickier): the softmax couples all tokens — increasing one weight forces others to decrease (since they sum to 1). This means $W_q$ and $W_k$ learn a push-pull dynamic: if the model repeatedly predicts “sat” after “cat” but gets it wrong, the gradient flows back through the softmax to nudge $W_q$ of “cat” and $W_k$ of “sat” to produce a higher dot product, while pushing down the keys of irrelevant tokens.

3. Competition between heads: each head’s $W_q$ and $W_k$ receive gradients that discourage redundancy — if two heads attend to the same pattern, their outputs are near-identical and the later projection $W_o$ can’t use both. The training dynamics push heads into complementary specializations automatically.

The full formula involves three chained matrix derivatives, but the intuition: weights that produce attention scores closer to the correct blend of context get rewarded; weights that amplify irrelevant context get penalized.

How attention behaves during inference

During generation, attention works differently than during training:

• Training: the model sees the full sequence at once. Attention scores are computed for every pair of tokens simultaneously (masked so tokens can’t peek ahead).
• Inference: the model generates one token at a time. At step 5, it has only 5 tokens. The newest token computes Q against all previous K, but those K don’t change — they were computed at earlier steps.

This is where the KV cache comes in. Instead of recomputing K and V for every token at every step, the model stores them:

1. Token 1: compute $K_1,V_1$, cache them. Generate token 2.
2. Token 2: compute $Q_2$, retrieve cached $K_1$. Compute score, retrieve cached $V_1$ for blending. Generate token 3.
3. Token 3: compute $Q_3$, retrieve cached $K_1,K_2$. Compute scores, retrieve $V_1,V_2$. Generate token 4.

This reduces the per-step computation from $O(N^2)$ to O(N) — linear in sequence length instead of quadratic. For a 128K-token context, this is the difference between milliseconds and seconds per token.

FlashAttention is an algorithmic optimization that makes the attention computation itself faster by being smarter about GPU memory access. Standard attention materializes the full N×N attention matrix in GPU memory, then reads it again for softmax and the weighted sum. FlashAttention instead processes attention in blocks that fit in fast GPU SRAM, never writing the full matrix to slow HBM. This makes attention 2–4× faster and uses much less memory — which is why it enabled the jump from 4K to 128K+ context lengths in production models. Almost every modern LLM training and inference stack uses FlashAttention or a variant of it.

The trained attention heads behave differently at inference time too: heads that learned syntactic patterns (e.g., “attend to the subject of this verb”) apply those patterns to whatever prompt you give, even if they’ve never seen that sentence before. Heads that learned copying behavior can attend to any token in the prompt and copy its value, which is how in-context learning works — the prompt’s examples become values that attention heads can retrieve.

06 · MLP & residual stream

Knowledge storage and the highway through the model

Attention is group work — tokens exchange information. The MLP is individual work — each token processes what it learned on its own.

Think of it like studying: attention is the group discussion where you hear others’ ideas; the MLP is the quiet desk time where you make sense of what you heard and connect it to what you already know.

The MLP has two linear layers with a bend in between:

First, the token’s vector is expanded into a wider space (4×–8× wider) — like taking notes and spreading them out on a bigger desk to see patterns. Then the activation function filters out noise (negative values become zero). Then it’s compressed back to the original size.

Modern LLMs use SwiGLU, which adds a third “gate” projection that learns what to let through and what to block — like having a smart filter that decides which notes are worth keeping.

def mlp_forward(x, W1, b1, W2, b2):
    # Expand: project into wider space (4x)
    hidden = x @ W1 + b1
    hidden = np.maximum(0, hidden)   # ReLU: kill negatives
    # Contract: project back to original size
    return hidden @ W2 + b2

def swiglu_forward(x, W1, W3, W2):
    # Gate: one path computes, another controls flow
    hidden = x @ W1
    gate = x @ W3
    gate = 1 / (1 + np.exp(-gate))   # sigmoid: values between 0-1
    gated = hidden * gate             # element-wise multiply
    return gated @ W2

# Toy example: 4-dim input → 8-dim hidden → 4-dim output
d_model, d_ff = 4, 8
x = np.array([0.5, -0.2, 0.8, 0.1])

W1 = np.random.randn(d_model, d_ff) * 0.1
b1 = np.zeros(d_ff)
W2 = np.random.randn(d_ff, d_model) * 0.1
b2 = np.zeros(d_model)

result = mlp_forward(x, W1, b1, W2, b2)
print(f"MLP input:  {x}")
print(f"MLP output: {result.round(3)}")

The MLP stores most of the model’s knowledge — facts, grammar patterns, style. Attention decides what context to gather; the MLP decides what to do with it.

The residual stream: after the MLP (and after attention), the original input is added back:

This is the residual connection. It acts like a highway: information flows through the model even if individual layers have weak signals. Every token starts as its embedding at layer 0, and each layer adds a refinement. The final representation is the sum of contributions from every layer.

def layer_norm(x, eps=1e-5):
    """Stabilize values: center (subtract mean), scale (divide by std)."""
    mean = x.mean(axis=-1, keepdims=True)
    std = x.std(axis=-1, keepdims=True)
    return (x - mean) / (std + eps)

def transformer_block(x, attn_weights, mlp_weights):
    """One complete transformer block = attention + MLP + residuals."""
    # Pre-norm: normalize before each sublayer
    xn = layer_norm(x)

    # Attention (gather context from other tokens)
    attn_out = attention(xn, xn, xn)  # self-attention: Q, K, V all = xn
    x = x + attn_out                   # residual: add back original

    # MLP (process individually)
    xn = layer_norm(x)
    mlp_out = mlp_forward(xn, *mlp_weights)
    x = x + mlp_out                    # residual: add back again

    return x

# One token passes through one block
x = np.array([0.5, -0.2, 0.8, 0.1])
dummy_mlp = (np.random.randn(4, 8) * 0.1, np.zeros(8),
             np.random.randn(8, 4) * 0.1, np.zeros(4))
result = transformer_block(x, None, dummy_mlp)
print(f"Input:  {x}")
print(f"Output: {result.round(3)}")
# Notice: output ≈ input + small refinements from each sublayer

Note: modern LLMs apply LayerNorm before each sublayer (attention and MLP), not after. This “pre-norm” configuration stabilizes training because the normalization happens on the clean input rather than on the potentially noisy output.

Output projection: from residual to vocabulary

After the last transformer block, the residual stream holds a vector for each token. To predict the next token, feed the last token’s residual through the same embedding table (transposed) to get logits, then softmax to get probabilities.

This is called tied embeddings — the input lookup table doubles as the output projection. Every LLM uses this or a variant.

def unembed(residual, embed_table):
    """Turn the final token's residual into vocabulary probabilities."""
    logits = residual @ embed_table.T  # dot product against every word vector
    return softmax(logits)             # probabilities over vocabulary

Why the transpose? The embedding table has shape (vocab_size, d_model). We need (d_model,) @ (vocab_size, d_model).T = (vocab_size,) — one score per word in the vocabulary.

The complete forward pass: everything wired together

Now we have every piece. Let’s stitch them into one function that takes token IDs and returns next-token probabilities. This is a tiny, working transformer with 2 blocks, 2 attention heads, and 32-dimensional vectors — small enough to train on a laptop.

Every function below was introduced in earlier sections.

import numpy as np
import math

# ─── From Section 0 (starter toolbox) ───
def softmax(x):
    ex = np.exp(x - np.max(x, axis=-1, keepdims=True))
    return ex / ex.sum(axis=-1, keepdims=True)

def layer_norm(x, eps=1e-5):
    return (x - x.mean(axis=-1, keepdims=True)) / (x.std(axis=-1, keepdims=True) + eps)

def relu(x):
    return np.maximum(0, x)

# ─── From Section 4 (positional encoding) ───
def positional_encoding(seq_len, d_model):
    PE = np.zeros((seq_len, d_model))
    for pos in range(seq_len):
        for i in range(d_model):
            denom = 10000 ** ((2 * (i // 2)) / d_model)
            if i % 2 == 0:
                PE[pos, i] = math.sin(pos / denom)
            else:
                PE[pos, i] = math.cos(pos / denom)
    return PE

# ─── From Section 5 (attention) ───
def attention(Q, K, V):
    d_k = K.shape[-1]
    scores = Q @ K.T / np.sqrt(d_k)
    weights = softmax(scores)
    return weights @ V

# ─── From Section 6 (MLP + block) ───
def mlp_forward(x, W1, b1, W2, b2):
    hidden = relu(x @ W1 + b1)
    return hidden @ W2 + b2

def transformer_block(x, W_q, W_k, W_v, W_o, W1, b1, W2, b2):
    # Pre-norm → attention → residual
    xn = layer_norm(x)
    Q = xn @ W_q
    K = xn @ W_k
    V = xn @ W_v
    attn_out = attention(Q, K, V)
    attn_out = attn_out @ W_o     # project concatenated heads back
    x = x + attn_out

    # Pre-norm → MLP → residual
    xn = layer_norm(x)
    mlp_out = mlp_forward(xn, W1, b1, W2, b2)
    x = x + mlp_out
    return x

# ─── NEW: unembedding (output projection) ───
def unembed(residual, embed_table):
    logits = residual @ embed_table.T
    return softmax(logits)

# ─── Complete forward pass ───
def forward(token_ids, embed_table, pos_enc, params, n_layers=2):
    seq_len = len(token_ids)
    d_model = embed_table.shape[1]

    # Step 1: Embed + position
    x = embed_table[token_ids]       # (seq_len, d_model)
    x = x + pos_enc[:seq_len]        # add position fingerprint

    # Step 2: N transformer blocks
    for layer in range(n_layers):
        p = params[layer]
        x = transformer_block(x, p['W_q'], p['W_k'], p['W_v'],
                              p['W_o'], p['W1'], p['b1'], p['W2'], p['b2'])

    # Step 3: Predict next token (use last position's residual)
    probs = unembed(x[-1], embed_table)
    return probs

# ─── Cross-entropy loss (Section 7) ───
def cross_entropy(probs, correct_idx):
    return -math.log(max(probs[correct_idx], 1e-10))

# ─── Instantiate a tiny model ───
VOCAB = ["the", "cat", "dog", "sat", "mat", "run", "."]
vocab_size = len(VOCAB)
d_model = 32
n_layers = 2
d_k = 16
d_ff = 64

np.random.seed(42)
E = np.random.randn(vocab_size, d_model) * 0.1  # tied embeddings
PE = positional_encoding(20, d_model)            # supports up to 20 tokens

params = []
for _ in range(n_layers):
    params.append({
        'W_q': np.random.randn(d_model, d_k) * 0.1,
        'W_k': np.random.randn(d_model, d_k) * 0.1,
        'W_v': np.random.randn(d_model, d_k) * 0.1,
        'W_o': np.random.randn(d_k, d_model) * 0.1,
        'W1':  np.random.randn(d_model, d_ff) * 0.1,
        'b1':  np.zeros(d_ff),
        'W2':  np.random.randn(d_ff, d_model) * 0.1,
        'b2':  np.zeros(d_model),
    })

# Test it: predict next token after "the cat"
ids = [0, 1]  # "the cat"
probs = forward(ids, E, PE, params)
print(f"After 'the cat', next-word probabilities:")
for i, p in enumerate(probs):
    print(f"  {VOCAB[i]}: {p:.1%}")

The next-token probabilities are random (weights are untrained). That’s the problem Section 07 solves.

07 · Training

Learning by predicting the next token

Every weight in the model starts as random noise. The model knows nothing — it’s as smart as a baby randomly babbling. How does it learn?

The training goal is embarrassingly simple: predict the next word. Given “the cat sat on the ___”, predict “mat”. That’s the entire objective.

For every position in a sentence, the model guesses — it assigns a probability to every word in its vocabulary. If the correct word gets high probability, good. If not, the model gets “penalized” and adjusts.

The penalty is cross-entropy loss:

Don’t fear the formula. It just means:

• If the model is 90% sure about the right word → loss ≈ 0.1 (small penalty)
• If the model is 90% sure about the wrong word → loss ≈ 2.3 (big penalty)

The model’s goal: minimize this loss across billions of sentences. That’s it.

import math

def cross_entropy(probs, correct_idx):
    """Penalize the model: low loss if correct token is likely, high if not."""
    p_correct = probs[correct_idx]          # probability assigned to correct word
    return -math.log(max(p_correct, 1e-10)) # log(1)=0, log(0.1)=2.3, log(0.01)=4.6

# Example: model predicts next word after "the cat sat on the"
# Correct answer is "mat", model gives it 80% probability
probs = np.array([0.05, 0.80, 0.03, 0.02, 0.10])  # [the, mat, sat, dog, run]
correct_idx = 1  # "mat"
loss = cross_entropy(probs, correct_idx)
print(f"Correct = mat, P(correct) = {probs[correct_idx]:.0%}, loss = {loss:.3f}")

# If model was only 10% confident about "mat":
probs_bad = np.array([0.70, 0.10, 0.05, 0.05, 0.10])
loss_bad = cross_entropy(probs_bad, correct_idx)
print(f"Correct = mat, P(correct) = {probs_bad[correct_idx]:.0%}, loss = {loss_bad:.3f}")

Backpropagation: once we know the loss, we need to adjust every weight to reduce it. The key insight: for each weight, we calculate “if I nudge this weight up by a tiny bit, does the loss go up or down?”

Weights that caused more error get bigger adjustments. The calculation flows backward through the entire network — from output → MLP → attention → embeddings — assigning blame to each weight along the way. It’s like tracing a mistake backward through a chain of reasoning: “the output was wrong because the MLP gave bad info, which happened because the attention gathered wrong context, which happened because the embedding was unclear…”

# Training the real forward() with numerical gradients

# Numerical gradient: (f(x+h) - f(x-h)) / (2h)
# This estimates how each weight affects the loss by nudging it
# up and down and measuring the difference. Slow but exact.

def numerical_grad_params(params, E, PE, token_ids, correct_idx, h=1e-4):
    """Return gradient of loss w.r.t. every parameter, one at a time."""
    base_loss = cross_entropy(forward(token_ids, E, PE, params), correct_idx)
    grads = []
    for layer_params in params:
        layer_grads = {}
        for key in layer_params:
            original = layer_params[key].copy()
            grad = np.zeros_like(original)
            # Perturb each element and measure loss change
            flat = original.ravel()
            grad_f = grad.ravel()
            for j in range(len(flat)):
                flat[j] += h
                loss_p = cross_entropy(forward(token_ids, E, PE, params), correct_idx)
                flat[j] -= 2 * h
                loss_m = cross_entropy(forward(token_ids, E, PE, params), correct_idx)
                flat[j] += h  # restore
                grad_f[j] = (loss_p - loss_m) / (2 * h)
            layer_grads[key] = grad
        grads.append(layer_grads)
    return grads

# ⚠ WARNING: this loop is SLOW — it calls forward() twice per parameter
# per step. We train only the embedding table E for demonstration.
# Real training uses automatic differentiation (PyTorch, JAX) which
# computes all gradients in one backward pass, not one param at a time.

learning_rate = 0.5
sentence_ids = [0, 1, 3, 4, 6]  # "the cat sat mat ."

for step in range(20):
    total_loss = 0
    for i in range(len(sentence_ids) - 1):
        context = sentence_ids[:i+1]
        correct = sentence_ids[i+1]

        # Forward
        probs = forward(context, E, PE, params)
        loss = cross_entropy(probs, correct)
        total_loss += loss

        # Numerical gradient
        grads = numerical_grad_params(params, E, PE, context, correct, h=1e-4)

        # Update weights
        for layer in range(len(params)):
            for key in grads[layer]:
                params[layer][key] -= learning_rate * grads[layer][key]

    if step % 5 == 0:
        test_probs = forward([0, 1], E, PE, params)  # "the cat" → ?
        print(f"step {step}: loss = {total_loss:.3f},  "
              f"P('sat'|'the cat') = {test_probs[3]:.1%}")

Notice: each training step calls forward() twice for every single weight — thousands of forward passes. This is why it’s called “numerical” and not used at scale. Real training uses backpropagation, which computes all gradients in a single backward pass by walking the chain rule backward through the computation graph. But the idea is identical: measure how nudging each weight changes the loss, then nudge in the direction that reduces it.

Training at scale: real LLM training involves:

• Batches: processing thousands of sequences at once, not one at a time
• Mixed precision: using 16-bit numbers instead of 32-bit for speed
• Learning rate schedule: starting with large weight updates and gradually shrinking them — like coarse-tuning then fine-tuning. A warmup phase slowly increases the rate for the first ~2000 steps (when gradients are chaotic), then cosine decay smoothly decreases it to near zero.
• Checkpointing: saving the model every few hours (training can and will crash)
• Data parallelism: splitting work across hundreds of GPUs, each computing gradients, then averaging them

08 · Using the trained model

From next-token predictor to assistant

After pre-training, the model is an expert at continuing text. Ask it a question and it will ramble — it was trained to predict the next word, not to answer helpfully.

To turn it into an assistant, we do fine-tuning: continued training on curated question-answer pairs.

Only the answer part gets penalized — the question is just context. The learning rate is tiny (0.1× or less) because we want to gently polish the model’s behavior, not retrain it from scratch.

# Generation loop: the model writes one word at a time
def generate(model, prompt_ids, n_tokens=5, T=1.0):
    generated = list(prompt_ids)
    for _ in range(n_tokens):
        probs = model(generated)           # predict next word probabilities
        probs = softmax_with_temp(np.log(probs + 1e-10), T)  # apply temperature
        next_id = np.random.choice(len(probs), p=probs)      # sample
        generated.append(next_id)
        if next_id == 0:                    # <end> token
            break
    return generated

# Toy example: a model that predicts next word based on last word
toy_vocab = ["<end>", "the", "cat", "sat", "mat"]

def toy_model(context):
    # In a real model, this would be the full transformer forward pass
    # Here we just return fixed probabilities for demo
    return np.array([0.01, 0.40, 0.30, 0.20, 0.09])

prompt = [1]  # start with "the"
output = generate(toy_model, prompt, n_tokens=4)
print(f"Generated: {' '.join(toy_vocab[i] for i in output)}")

This is autoregressive generation: each new word becomes part of the input for the next prediction. The model never goes back to rewrite — it’s like writing a sentence one word at a time, never erasing.

How to control what the model says (inference techniques)

The model assigns probabilities to every possible next word. How you pick from those probabilities changes the output:

• Temperature (T): controls randomness. Low T (0.1) = always picks the most likely word (boring but safe). High T (1.5) = sometimes picks less likely words (creative but risky). T=1 is the default.
• Top-k: only consider the k most likely words, ignore the rest.
• Top-p: consider words until their combined probability reaches p (e.g., 90%), cut off the long tail of unlikely words.

# Temperature: adjust how "sharp" the probability distribution is
def apply_temperature(logits, T):
    return logits / T                       # divide raw scores by T

def softmax_with_temp(logits, T=1.0):
    logits = logits / T                     # apply temperature
    ex = np.exp(logits - np.max(logits))
    return ex / ex.sum()

# Example: raw scores (logits) for next word
logits = np.array([3.2, 1.1, 0.5, -0.2])   # raw: higher = more likely

for T in [0.1, 1.0, 2.0]:
    probs = softmax_with_temp(logits, T)
    print(f"T={T:.1f}: {[f'{p:.0%}' for p in probs]}")
# T=0.1: [~100%, ~0%, ~0%, ~0%]  — almost always picks top word
# T=1.0: [76%, 14%, 8%, 2%]       — balanced
# T=2.0: [48%, 22%, 16%, 14%]     — more random

# Top-k: keep only the k most probable tokens, zero out the rest
def top_k_sampling(probs, k=2):
    top_indices = np.argsort(probs)[-k:]    # indices of k highest
    new_probs = np.zeros_like(probs)
    new_probs[top_indices] = probs[top_indices]
    return new_probs / new_probs.sum()

# Top-p (nucleus): keep tokens until cumulative probability reaches p
def top_p_sampling(probs, p=0.9):
    sorted_idx = np.argsort(probs)[::-1]
    cumsum = np.cumsum(probs[sorted_idx])
    cutoff = np.searchsorted(cumsum, p) + 1
    new_probs = np.zeros_like(probs)
    new_probs[sorted_idx[:cutoff]] = probs[sorted_idx[:cutoff]]
    return new_probs / new_probs.sum()

probs = np.array([0.50, 0.20, 0.15, 0.10, 0.05])
print(f"\nTop-2: {top_k_sampling(probs, 2).round(2)}")   # keeps only top 2
print(f"Top-p (0.9): {top_p_sampling(probs, 0.9).round(2)}")  # keeps ~90% prob mass

Three ways to steer the model (without retraining)

1. Fine-tuning: change the weights by training on curated examples
2. Prompting: don’t change weights — just write a good instruction. The model figures out what you want from the context alone. This is called in-context learning — giving examples in the prompt works because attention can use those examples.
3. Chain-of-thought (CoT): ask the model to “think step by step.” The intermediate reasoning steps give attention more material to work with, dramatically improving math and logic problems.
4. RAG (Retrieval-Augmented Generation): before answering, look up relevant documents and inject them into the prompt. The model then answers based on real, up-to-date info instead of relying on its imperfect memory. This is how production chatbots stay factual.

09 · Interpretability

Looking inside the black box

After all this, the model is just a bunch of numbers — weights and activations. But patterns emerge that we can visualize and interpret.

Inspecting our own model

We built a tiny transformer in Section 06. Let’s modify forward() to also return the attention weights — the “who-looks-at-whom” scores — and print them.

def forward_with_attn(token_ids, embed_table, pos_enc, params, n_layers=2):
    """Same as forward(), but also returns attention weights per layer."""
    seq_len = len(token_ids)
    d_model = embed_table.shape[1]

    x = embed_table[token_ids]
    x = x + pos_enc[:seq_len]

    all_attn = []
    for layer in range(n_layers):
        p = params[layer]
        xn = layer_norm(x)
        Q = xn @ p['W_q']
        K = xn @ p['W_k']
        V = xn @ p['W_v']

        d_k = K.shape[-1]
        scores = Q @ K.T / np.sqrt(d_k)
        attn_weights = softmax(scores)     # (seq_len, seq_len)
        all_attn.append(attn_weights)

        attn_out = attn_weights @ V
        attn_out = attn_out @ p['W_o']
        x = x + attn_out

        xn = layer_norm(x)
        x = x + mlp_forward(xn, p['W1'], p['b1'], p['W2'], p['b2'])

    probs = softmax(x[-1] @ embed_table.T)
    return probs, all_attn

# Inspect: "the cat sat mat ."
ids = [0, 1, 3, 4, 6]
probs, attn_layers = forward_with_attn(ids, E, PE, params)

VOCAB = ["the", "cat", "dog", "sat", "mat", "run", "."]
words = [VOCAB[i] for i in ids]

for layer_idx, attn in enumerate(attn_layers):
    print(f"\nLayer {layer_idx} attention matrix:")
    print(f"         " + "  ".join(f"{w:>6s}" for w in words))
    for i, word in enumerate(words):
        row = "  ".join(f"{attn[i,j]:6.2f}" for j in range(len(ids)))
        print(f"  {word:6s}  {row}")