DEV Community: Sachin Kr. Rajput

A Forward Pass for One Token: Where the FLOPs Go

Sachin Kr. Rajput — Mon, 26 Oct 2026 10:02:01 +0000

In Part 1 we watched the Word Factory answer a question. The Typewriter wrote the answer one token at a time, and every single token needed a full trip through the machine.

Today we climb inside and ride along with one token on that trip.

We'll count every bit of math it causes, room by room. By the end you'll be able to look at any language model and guess how much work one word costs, with a rule small enough to fit on a sticky note.

And I'll tell you about a bug this counting helped me catch in Part 1.

How do we count "work"?

Inside a language model, almost everything is one move repeated billions of times: multiply two numbers, then add the result to a running total.

Engineers count these tiny math steps as FLOPs (floating-point operations). One multiply is 1 FLOP. One add is 1 FLOP.

Here's a tiny example. A word is turned into a list of numbers (just 2 here, to keep it small), and we push it through a table of numbers with 2 rows and 3 columns:

word:    [2, 3]

table:   [ 1, 0, 4 ]
         [ 5, 2, 1 ]

result:  [ 2×1 + 3×5,   2×0 + 3×2,   2×4 + 3×1 ]  =  [17, 6, 11]

Look at what happened: every number in the table got used exactly once, multiplied by something and added to a total. The table has 6 numbers, so that's 6 multiply-and-adds.

Chips do a multiply-and-add as one combined move, and everyone counts it as 2 FLOPs. So:

work for one word through a table = 2 × (rows × columns) FLOPs

This table cost 2 × 6 = 12 FLOPs. (Strictly there were only 3 additions here, not 6, but counting 2 per table number keeps the bookkeeping simple, and for big tables the difference is tiny.)

Hold on to that rule. It's going to do almost all the work today.

A map of the inside

The Reader and the Typewriter from Part 1 are the same machine, so there's only one building to explore. Our model has the shape of GPT-2 small, a real model from 2019. Here's the building:

           ┌────────────────────────────────────────┐
           │  SCOREBOARD: score all 50,257 words    │  → next word
           └────────────────────────────────────────┘
                               ▲
           ┌────────────────────┬───────────────────┐
 floor 12  │  look-back room    │    thinking room  │
           └────────────────────┴───────────────────┘
                               ▲
                    ... 10 more floors ...
                               ▲
           ┌────────────────────┬───────────────────┐
 floor 1   │  look-back room    │    thinking room  │
           └────────────────────┴───────────────────┘
                               ▲
           ┌────────────────────────────────────────┐
           │  DICTIONARY: word → 768 numbers        │
           └────────────────────────────────────────┘
                               ▲
                             "sky"

Our token rides the elevator from the bottom to the top. Let's follow it.

Stop 1: The dictionary

First, the token needs a backpack of numbers. The model has a dictionary with one row per word piece: 50,257 rows, each holding 768 numbers. The token number for "sky" says which row to grab.

That's it. No multiplying, just finding the right row. Like looking up a word in an actual dictionary.

Work: 0 FLOPs. (Practically free.)

From now on, our token is those 768 numbers. Every room it visits changes them a little.

Stop 2: The look-back room (attention)

This is where the token looks back at the words that came before it. Engineers call it attention. It happens in four steps.

a) Make three cards. The token's 768 numbers are pushed through one big table (768 rows × 2,304 columns) to make three cards:

a Question card: "what am I looking for?"
a Label card: "what am I about?"
a Contents card: "what can I share?"

b) Write in the notebook. The Label and Contents cards go into the notebook (the KV cache from Part 1). Every earlier token already left its own two cards there.

c) Read the notebook. The token holds up its Question card and compares it with every Label in the notebook, to see which earlier words matter most right now. Then it blends together their Contents cards: lots of the important ones, a little of the rest.

d) Tidy up. The blend goes through one more table (768 × 768) to become 768 numbers again, which get added to the token's backpack.

Here's the work, using our rule:

make 3 cards:      2 × 768 × 2,304  = 3.54 million FLOPs
tidy-up table:     2 × 768 × 768    = 1.18 million FLOPs
read the notebook: 4 × 768 × (number of notes in the notebook)

That last line is different. It doesn't depend on a fixed table. It depends on how long the conversation is so far. With 128 notes it's 0.39 million FLOPs. With 8,192 notes it's 25 million. Keep an eye on it.

Stop 3: The thinking room (MLP)

Next, the token thinks on its own for a bit. Engineers call this the MLP (multi-layer perceptron, a fancy name for "two tables with a squish in the middle").

a) Expand. The 768 numbers go through a table that makes them 4 times bigger: 3,072 numbers. More room to think.

b) Squish. Each of the 3,072 numbers goes through a little curve that shrinks small and negative values. (This is what lets the model learn things more interesting than straight lines.)

c) Shrink back. Another table squeezes the 3,072 numbers back down to 768, and they're added to the backpack.

expand:  2 × 768 × 3,072  = 4.72 million FLOPs
shrink:  2 × 3,072 × 768  = 4.72 million FLOPs

Notice something? The thinking room's two tables are twice as big as the look-back room's tables. So in every floor, thinking is about two-thirds of the table work.

Then the elevator goes up, and the token does it all again. 12 floors in total.

Stop 4: The scoreboard

At the top, the token's 768 numbers are compared with every word in the dictionary to give each one a score. The highest score becomes the next word in the answer.

That's one more table: 768 × 50,257.

scoreboard:  2 × 768 × 50,257 = 77.2 million FLOPs

That's a big one. Remember it.

Let's count it for real

Counting on paper is nice, but let's measure. I wrote forward_pass.py: the same GPT-2-small-shaped building in about 120 lines of NumPy. Like in Part 1, its numbers are random, so its words are gibberish, but it does exactly the same amount of math as a trained model of this shape.

The trick is one little helper. Every time the token goes through a table, it adds up the FLOPs and starts a stopwatch:

def mm(part, x, W):
    """Push the token's numbers x through table W, and keep score."""
    t = time.perf_counter()
    y = x @ W                                    # the actual math
    secs[part] += time.perf_counter() - t        # how long it took
    flops[part] += 2 * W.shape[0] * W.shape[1]   # 2 x rows x columns
    return y

And here's one floor of the building. You don't need to follow every line, just spot the two rooms:

for p in self.layers:                            # 12 floors
    # 👀 Look-back room
    h = norm(x)
    q, k, v = np.split(mm("make Q, K, V", h, p["qkv"]), 3)
    scores = compare(q, notebook_labels)         # read the notebook...
    mix = blend(scores, notebook_contents)       # ...and blend
    x = x + mm("attention output", mix, p["out"])
    # 🧠 Thinking room
    h = norm(x)
    up = gelu(mm("expand 768 -> 3072", h, p["up"]))
    x = x + mm("shrink 3072 -> 768", up, p["down"])

(The real file uses the same steps with a bit more bookkeeping. The notebook is pre-filled with notes from earlier tokens, so we can test different conversation lengths.)

I ran one token at a time through the building, 30 times, repeated the whole thing 3 times, and kept the middle run. Plain CPU, 2 cores, no GPU.

The results

Here's one token with 128 notes in the notebook (a short chat):

notebook holds 128 tokens | one token = 252.2 million FLOPs, 35.3 ms
part                                       MFLOPs  work %      ms  time %
1. look up the token                         0.00    0.0%    0.01    0.0%
2. attention: make Q, K, V                  42.47   16.8%    4.85   13.7%
3. attention: read the notebook              4.72    1.9%    1.95    5.5%
4. attention: output                        14.16    5.6%    1.91    5.4%
5. thinking: expand 768 -> 3072             56.62   22.4%    6.26   17.7%
6. thinking: shrink 3072 -> 768             56.62   22.4%    6.39   18.1%
7. score every word in the dictionary       77.19   30.6%    8.29   23.5%
other (norms, softmax, adds)                 0.45    0.2%    5.68   16.1%

("MFLOPs" means millions of FLOPs, added up over all 12 floors.)

And here's how it changes as the conversation gets longer:

notes in notebook      128      1,024      8,192
FLOPs per token      252.2 M    285.6 M    553.0 M
time per token       35.3 ms    44.7 ms   105.0 ms
notebook share of the work
                       1.9%      13.2%      54.6%

Five things we just discovered

1. It's (almost) all tables

Add up the table steps (2, 4, 5, 6 and 7): 247 of the 252 million FLOPs, 98%. Reading the notebook was 1.9%. Everything else (the squish curve, tidying, adding) was 0.2%.

So when people say "AI is just multiplying matrices", they're nearly right. The tables are the model. These tables are what engineers call parameters (or weights): the numbers the model learned during training.

2. The sticky-note rule

Since every table number gets used once per token, here's a shortcut. Count the numbers in all the tables:

12 floors × 7,077,888 numbers per floor  =  84,934,656
scoreboard: 768 × 50,257                 =  38,597,376
                                    total = 123,532,032  (about 124 million)

Double it: 247.1 million FLOPs. The measured total was 252.2 million. The small difference is the notebook reading plus the tiny "other" steps.

That gives us the rule for the sticky note:

FLOPs per token ≈ 2 × parameters  +  4 × floors × width × notes in the notebook

For a model with 8 billion parameters, that's about 2 × 8 = 16 billion FLOPs for every single word it writes (a little less, really, because the dictionary lookup is free), before you even count the notebook.

3. Long chats change where the work goes

With 128 notes, reading the notebook was a rounding error (1.9%). With 8,192 notes, it was more than half of all the work (54.6%) and two-thirds of the time (66%).

The tables don't change size during a chat. The notebook does. Every new word adds another note that every future word has to read. That's why really long conversations and long documents get slow and expensive, and why so much engineering goes into making the notebook smaller and faster to read.

4. The scoreboard is huge in small models

For our little model, scoring 50,257 words was 30.6% of all the work: the single biggest step! That's because the dictionary is big and the rest of the model is small.

In big models, it flips. I fed the shape of a popular 8-billion-parameter open model (Llama 3.1 8B) into a small calculator, flop_counter.py, that uses the same rule:

Llama 3.1 8B shape
  notebook    128:  15.08 GFLOPs | attention tables 17.8%  read the notebook  0.4%  thinking tables 74.8%  score the dictionary  7.0%
  notebook  8,192:  19.30 GFLOPs | attention tables 13.9%  read the notebook 22.2%  thinking tables 58.4%  score the dictionary  5.4%
  notebook 32,768:  32.19 GFLOPs | attention tables  8.3%  read the notebook 53.4%  thinking tables 35.0%  score the dictionary  3.3%

(GFLOPs means billions of FLOPs.) The scoreboard drops to about 7%, and the thinking room takes three-quarters of the work. And look at the last line: with a 32,768-token conversation, the notebook is again more than half the work, even for a big model.

5. Counting work is not the same as counting time

Look at the "other" row again: 0.2% of the work, but 16.1% of the time. Those little steps do almost no math, but each one is a separate small job, and every job has a setup cost. It's like a delivery driver whose packages weigh nothing but who still has to park at every house. Real inference engines glue small steps together (engineers call it kernel fusion) for a similar reason.

And here's a bigger surprise. Our CPU did 252 million FLOPs in 35 ms. That's about 7 billion FLOPs per second. Is that the CPU's top speed?

Let's check with bus_test.py. It pushes tokens through 36 tables of 768 × 3,072 numbers (340 MB, the same size as all 12 floors), first one token per trip, then many tokens per trip:

for batch in [1, 16, 64, 256]:
    x = rng.standard_normal((batch, 768)).astype(np.float32)
    ...
    for W in tables:          # one trip through every table
        y = x @ W

tables: 340 MB
   1 tokens per trip:    17.4 ms per trip,  17.36 ms per token,    9.8 GFLOP/s
  16 tokens per trip:    82.7 ms per trip,   5.17 ms per token,   32.9 GFLOP/s
  64 tokens per trip:   106.9 ms per trip,   1.67 ms per token,  101.7 GFLOP/s
 256 tokens per trip:   299.7 ms per trip,   1.17 ms per token,  145.1 GFLOP/s

Same tables. Same CPU. 15 times more math per second when 256 tokens share the trip.

It's the school bus from Part 1 again. With one token per trip, the CPU spends almost all its time hauling the 340 MB of tables from memory, and very little time doing math. With 256 tokens, it hauls the tables once and does 256 tokens' worth of math with them.

So the FLOP count tells you how much math a word needs. It does not tell you how long it takes. For typing one word at a time, the clock is mostly set by how many bytes have to be moved, not by how much math there is. That's the whole story of Part 3.

How counting FLOPs caught a bug in Part 1

When I first ran this post's code, the tables took more than 6 times longer inside the model than they did in the bus test. Same tables, same CPU. With FLOPs counted, the gap was impossible to miss.

The cause was one innocent-looking line:

scores = q @ K / np.sqrt(64)     # np.sqrt(64) is a 64-bit NumPy number...

In recent versions of NumPy, dividing our 32-bit numbers by a 64-bit NumPy number quietly turns the result into 64-bit numbers. Every step after that worked with 64-bit numbers too, which meant converting every table on every token. Twice the bytes, plus extra copying.

The fix is tiny:

scores = q @ K / 64 ** 0.5       # a plain Python number keeps everything 32-bit

The same line was in Part 1's tiny engine. With it fixed, typing went from 107 ms to about 27 ms per token, and reading went from 6.5 ms to about 4.4 ms per token. So reading is about 6 times cheaper per token than typing, not 16 times like I first wrote. The lessons in Part 1 all still hold, but I've corrected the numbers there and added a note.

The code for this post now has a guard so it can't happen silently again:

assert x.dtype == np.float32, "something switched us to 64-bit math"

Lesson: if your measured speed is far from what you'd expect for the work, don't trust the measurement. Go find the reason.

Cheat sheet

Stop	Real name	Work for one token	Grows with the chat?
Dictionary	Embedding lookup	~0	No
Make 3 cards	Q, K, V projection	2 × width × 3 × width	No
Read the notebook	Attention over the KV cache	4 × width × notes	Yes
Tidy up	Attention output projection	2 × width × width	No
Thinking room	MLP	2 × width × 4 × width × 2 tables	No
Scoreboard	LM head (unembedding)	2 × width × dictionary size	No

Rule	What it says
Sticky note	FLOPs per token ≈ 2 × parameters + 4 × floors × width × notes
Small models	The scoreboard can be the biggest single step
Big models	The thinking room is most of the work
Long chats	Reading the notebook takes over
One token at a time	Time is set by moving bytes, not by doing math

What's next

We just saw that the same CPU does 15 times more math per second when tokens share the trip. Part 3: Prefill vs Decode: Two Phases, Two Bottlenecks digs into exactly that: why reading your question is limited by math, why typing the answer is limited by memory, and what that means for every chatbot server.

Every number in this post comes from actually running these three scripts (a few snippets above are shortened for reading; the files are the exact code I ran): forward_pass.py, bus_test.py and flop_counter.py. The tiny model used random weights on a 2-core CPU, so it shows the pattern, not the speed of any real GPU.

One Request, End to End: Tokenize, Prefill, Decode, Stream

Sachin Kr. Rajput — Thu, 22 Oct 2026 10:13:21 +0000

Correction (Oct 26, 2026): an earlier version of this post had a bug in the tiny engine. One line (/ np.sqrt(self.hd)) quietly switched the math from 32-bit to slower 64-bit numbers, which made typing look about 4 times slower than it really is. I fixed the line, re-ran everything, and updated every number below. The headline changed: reading is about 6 times cheaper per token than typing, not 16. The lessons are the same. Part 2 explains how I caught it.

You type "Why is the sky blue?" and press Enter.

There's a short pause. Then words start appearing, one after another, like someone is typing to you very fast.

What happened during that pause? And why do the words come out one at a time instead of all at once?

Let's take a walk through the Word Factory.

Welcome to the Word Factory

Picture a small factory with a conveyor belt. Your question goes in at one end, and an answer comes out on your phone at the other. Along the way there are four machines:

The Chopper cuts your question into small pieces.
The Reader has lots of eyes and reads every piece at the same time.
The Typewriter writes the answer, but only one piece at a time.
The Delivery Tape carries each piece to your phone the moment it's typed.

There's also a notebook hanging above the Reader. It's the factory's secret weapon, and it matters a lot later.

Engineers have proper names for these four steps: tokenize, prefill, decode and stream. That's all an "inference engine" is: the factory that runs a language model and turns your question into an answer.

Let's visit each machine.

Machine 1: The Chopper (tokenize)

A language model can't read letters the way you do. It only understands numbers. So the first job is to cut your text into small pieces called tokens and give each piece a number.

A token is usually a word or part of a word. Think of LEGO bricks: common words are one brick, and long or rare words get built from a few smaller bricks.

"Why is the sky blue?"
   ↓ chop
[ "Why" ][ " is" ][ " the" ][ " sky" ][ " blue" ][ "?" ]
   ↓ give each a number
[  512,     87,     33,    2041,     978,     11  ]

(These numbers are made up for the example. Each model has its own dictionary that decides which number goes with which piece.)

How long does it take? Almost no time at all: a tiny fraction of a millisecond. The Chopper is never the slow machine.

Machine 2: The Reader (prefill)

Now the Reader looks at your whole question.

Here's the clever part: the Reader has lots of eyes, and it reads every token at the same time. It doesn't go word by word. Whether your question is 6 tokens or 600, it reads them all in one big gulp.

While it reads, it writes notes in the notebook: for every token, a short note about what that token means and how it connects to the others. Engineers call this notebook the KV cache. (KV stands for "keys and values", which is just what the notes are called inside the model.)

When the Reader finishes, it has a guess for the very first word of the answer.

This step is called prefill because it fills the notebook with notes about your whole question before any answer is written.

The important bit: a longer question means more to read. So the longer your question, the longer you wait for the first word.

Machine 3: The Typewriter (decode)

Now the answer gets written, and this is where things slow down.

The Typewriter can only type one token per step. Why? Because each new word depends on the words before it. You can't decide the 5th word of a sentence until you know the first 4. Try finishing "The sky is blue because..." without knowing the word "because" was there.

So for every single token, the factory does a full trip:

Look at the notebook (all the notes so far)
Run through the whole model to pick the next token
Add one new line to the notebook
Type the token
Go back to step 1 for the next token

That's the pink "again & again" arrow in the picture. Engineers call this the decode loop.

Why the notebook is so useful: without it, the Typewriter would have to re-read your entire question and everything it has typed so far, every single step. That would be like re-reading a whole book from page 1 every time you want to write the next word. With the notebook, it just glances at its notes.

Machine 4: The Delivery Tape (stream)

The moment a token is typed, it's turned back into text and sent straight to your screen. The factory doesn't wait to finish the whole answer first.

That's why you see words appear one by one. It's not a fancy animation. You're watching the Typewriter work in real time.

This is called streaming, and it makes chatbots feel much faster. You start reading while the rest is still being written.

Three stopwatches

So how do we say if a factory is "fast"? There are three different stopwatches, and they measure different things.

⏱ Stopwatch 1: Time to First Token (TTFT)
How long from pressing Enter until the first word shows up. This is the "awkward pause." It includes chopping, the whole reading step, and picking the first token.

TTFT = time first token arrives − time you pressed Enter

⏱ Stopwatch 2: Inter-Token Latency (ITL)
The gap between one word and the next. If it's small and steady, the text flows smoothly. If one gap is suddenly huge, it feels like the chatbot froze, even if the average is fine.

ITL = time token 2 arrives − time token 1 arrived   (and so on, for every gap)

⏱ Stopwatch 3: Throughput (tokens per second)
How many tokens come out per second. Careful, there are two kinds:

For one person: how fast your answer is typed.
For the whole factory: how many tokens it produces for everyone at once.

A factory serving 100 people at 5,000 tokens per second in total gives each person about 50 tokens per second. Both numbers are true. They just answer different questions.

And the total wait for your full answer is simply:

total time ≈ TTFT + (number of tokens − 1) × ITL

Let's build a tiny factory and time it

Talking is fine, but let's measure. Below is a complete Word Factory in about 120 lines of Python, using only NumPy.

It has the same shape as a small real model (GPT-2 small: 12 layers, 85 million numbers inside). One honest warning: its "brain" is filled with random numbers, so the words it types are gibberish. But it does exactly the same amount of work as a trained model of that size, so the stopwatch readings are real.

The Chopper

To keep it simple, our Chopper cuts text into single bytes (roughly, one letter = one token).

import codecs
import numpy as np

class ByteTokenizer:
    vocab_size = 256

    def encode(self, text: str) -> list[int]:
        # Chop: every byte of the text becomes one token number.
        return list(text.encode("utf-8"))

    def stream_decoder(self):
        # Turn token numbers back into text, one at a time.
        # (It politely waits if a letter like "é" is split across two tokens.)
        dec = codecs.getincrementaldecoder("utf-8")(errors="replace")
        return lambda tok: dec.decode(bytes([tok]))

The Reader and the Typewriter are the same machine

Here's a surprise: reading and typing use the same machine. The only difference is how many tokens you feed it. Reading (prefill) feeds it the whole question at once. Typing (decode) feeds it one token at a time.

You don't need to follow every line. Look for the two comments marked 📓: that's where the notebook gets written and read.

class TinyDecoder:
    def __init__(self, n_layers=12, d=768, n_heads=12, d_ff=3072,
                 vocab=256, max_seq=4096, seed=0):
        rng = np.random.default_rng(seed)
        w = lambda *s: (rng.standard_normal(s) * 0.02).astype(np.float32)
        self.L, self.d, self.H, self.hd = n_layers, d, n_heads, d // n_heads
        self.tok_emb, self.pos_emb = w(vocab, d), w(max_seq, d)
        self.layers = [dict(qkv=w(d, 3 * d), o=w(d, d),
                            up=w(d, d_ff), down=w(d_ff, d))
                       for _ in range(n_layers)]
        self.lm_head = w(d, vocab)
        # The notebook (KV cache): room for notes on up to 4096 tokens.
        self.k = np.zeros((n_layers, n_heads, max_seq, self.hd), np.float32)
        self.v = np.zeros_like(self.k)
        self.n_cached = 0

    @staticmethod
    def _norm(x):
        return (x - x.mean(-1, keepdims=True)) / np.sqrt(x.var(-1, keepdims=True) + 1e-5)

    def forward(self, ids: np.ndarray) -> np.ndarray:
        """Feed in NEW tokens (the whole question, or just one).
        Returns the scores for what the next token should be."""
        T, start = len(ids), self.n_cached
        x = self.tok_emb[ids] + self.pos_emb[start:start + T]
        for i, p in enumerate(self.layers):
            q, k, v = np.split(self._norm(x) @ p["qkv"], 3, axis=-1)
            q, k, v = (t.reshape(T, self.H, self.hd).transpose(1, 0, 2) for t in (q, k, v))
            # 📓 write notes for the new tokens into the notebook
            self.k[i, :, start:start + T], self.v[i, :, start:start + T] = k, v
            # 📓 read ALL notes so far (old + new)
            K, V = self.k[i, :, :start + T], self.v[i, :, :start + T]
            att = q @ K.transpose(0, 2, 1) / self.hd ** 0.5
            mask = np.triu(np.ones((T, start + T), bool), k=start + 1)  # no peeking ahead
            att = np.where(mask, -1e9, att)
            att = np.exp(att - att.max(-1, keepdims=True))
            att /= att.sum(-1, keepdims=True)
            y = (att @ V).transpose(1, 0, 2).reshape(T, self.d)
            x = x + y @ p["o"]
            h = self._norm(x) @ p["up"]
            x = x + (0.5 * h * (1 + np.tanh(0.79788456 * (h + 0.044715 * h**3)))) @ p["down"]
        self.n_cached += T
        return self._norm(x[-1]) @ self.lm_head

The whole factory, start to finish

This little function is the conveyor belt. Read it top to bottom and you'll spot all four machines:

import time

def generate(model, tok, prompt: str, max_new_tokens: int):
    """Yield (event, text, timestamp) for one request, end to end."""
    model.n_cached = 0
    ids = np.array(tok.encode(prompt))                     # 1. CHOP
    yield "tokenized", "", time.perf_counter()
    logits = model.forward(ids)                            # 2. READ it all at once
    yield "prefilled", "", time.perf_counter()
    detok = tok.stream_decoder()
    for _ in range(max_new_tokens):                        # 3. TYPE, one at a time...
        nxt = int(np.argmax(logits))                       #    pick the best next token
        yield "token", detok(nxt), time.perf_counter()     # 4. SEND it right away
        logits = model.forward(np.array([nxt]))            #    ...and again

The full file, including a small measure() helper that reads the stopwatches, is tiny_engine.py (download it and run python tiny_engine.py).

The results

I ran the factory with three question sizes (64, 256 and 1,024 tokens), asked for 64 tokens of answer each time, and did one warm-up run first. Plain CPU, 2 cores, no GPU. I ran the whole thing three times and show the middle run:

 question  chop ms  read ms   TTFT ms  gap p50  gap p99  total s  tokens/s
     64      0.02     297.4    297.5     24.4     39.1     1.91     39.2
    256      0.02    1132.5   1132.6     27.2     55.7     3.03     33.1
   1024      0.06    5792.8   5793.4     32.7     62.3     7.99     28.7

("gap p50" is the typical gap between words. "gap p99" is the gap that only 1 in 100 steps is slower than: the "worst normal" gap.)

Here's the 256-token run drawn on a timeline:

 press Enter                                                               done
 │                                                                          │
 ├─ chop: 0.02 ms (too small to see)
 │
 ├────────── READ the whole question: 1,133 ms ──────────┤
 │                                                        │
 │                                  first word appears ──▶ █ █ █ █ █ █ █ ... █
 │                                                         └─┘
 │                                                        gap ≈ 27 ms, 63 times
 │◀──────────────── TTFT = 1.13 s ───────────────────────▶│◀── typing: 1.9 s ──▶│

Four things we just discovered

1. Reading is about 6 times cheaper per token than typing

For the 256-token question, reading took 1,133 ms. That's about 4.4 ms per token.

Typing took 27 ms per token.

Same machine, same computer. Why the huge difference?

Every trip through the model means hauling all 85 million numbers (about 340 MB) through the processor. When the Reader does that trip, it handles 256 tokens at once, so the cost is shared 256 ways. When the Typewriter does it, the whole trip is for one token.

It's like a school bus. Driving the bus across town costs about the same whether it carries 1 kid or 50. The Reader fills the bus. The Typewriter drives it across town with one passenger, over and over again.

This is the single most important idea in making chatbots fast, and it's why real servers try to type for many people at once (filling the bus). Later posts in this series dig into it.

2. The awkward pause is almost all reading

Chopping took 0.02 to 0.06 ms: basically nothing. TTFT was almost exactly the reading time.

And reading time grew with the question: 0.30 s for 64 tokens, 5.8 s for 1,024 tokens. A question 16 times longer made the wait about 19 times longer. (It grows a bit faster than the question, because every token has to "look at" every other token.)

Lesson: if the first word takes too long, the first suspect is a long question.

3. The notebook gets fatter, so typing slowly gets slower

The typical gap between words grew from 24 ms (short question) to 33 ms (long question).

That's the notebook. Every step, the Typewriter reads all the notes. Every token adds a new line. For our little model, each token adds about 72 KB of notes:

notes per token = 2 (K and V) × 12 layers × 768 numbers × 4 bytes = 73,728 bytes ≈ 72 KB

After about 1,100 tokens, that's roughly 80 MB of notes, read again on every single step. Bigger models have much bigger notebooks, and that's what decides how many people one GPU can serve at a time.

4. Most of the total time is typing

For 256 tokens in and 64 tokens out: 1.13 s of reading + 63 gaps of typing. With the typical 27 ms gap that's about 2.9 s. The real measurement was 3.03 s, because a few gaps were slower than typical.

Typing was about 63% of the total, even though the answer (64 tokens) was a quarter the length of the question (256 tokens). So the length of the answer matters a lot for total waiting time. That's why chatbots that "think" for thousands of tokens before answering are so expensive to run.

Bonus: time a real chatbot server

Measuring inside the factory flatters it. It doesn't count time spent waiting in line or travelling over the network. The fair way is to measure from your side, like a real user.

This script (measure_client.py) works with any server that speaks the OpenAI API format, including popular engines like vLLM, SGLang and TensorRT-LLM:

import argparse, statistics, time
from openai import OpenAI

ap = argparse.ArgumentParser()
ap.add_argument("--base-url", default="http://localhost:8000/v1")
ap.add_argument("--model", default="tiny")
ap.add_argument("--max-tokens", type=int, default=64)
ap.add_argument("--prompt", default="Explain how a large language model answers "
                                    "a question, step by step. " * 4)
args = ap.parse_args()
client = OpenAI(base_url=args.base_url, api_key="unused")

t_send = time.perf_counter()                       # press Enter
stream = client.completions.create(model=args.model, prompt=args.prompt,
                                   max_tokens=args.max_tokens, temperature=0,
                                   stream=True,
                                   stream_options={"include_usage": True})
arrivals, usage = [], None
for chunk in stream:
    if chunk.choices:                              # a token arrived!
        arrivals.append(time.perf_counter())
    if chunk.usage:
        usage = chunk.usage
t_done = time.perf_counter()

ttft = arrivals[0] - t_send
gaps = [b - a for a, b in zip(arrivals, arrivals[1:])]
out_tokens = usage.completion_tokens if usage else len(arrivals)
print(f"prompt tokens      : {usage.prompt_tokens if usage else '?'}")
print(f"output tokens      : {out_tokens}")
print(f"TTFT               : {ttft*1e3:8.1f} ms")
print(f"ITL p50 / p99      : {statistics.median(gaps)*1e3:8.1f} / "
      f"{sorted(gaps)[int(0.99*(len(gaps)-1))]*1e3:.1f} ms")
print(f"end-to-end latency : {t_done - t_send:8.2f} s")
print(f"decode tokens/s    : {(out_tokens - 1) / (arrivals[-1] - arrivals[0]):8.1f}")

To test it without a GPU, I wrapped our tiny factory in a small web server (tiny_server.py, 42 lines) and pointed the script at it:

prompt tokens      : 276
output tokens      : 64
TTFT               :   1265.1 ms
ITL p50 / p99      :     27.4 / 41.5 ms
end-to-end latency :     3.11 s
decode tokens/s    :     34.7

It matches the factory's own numbers, plus a little extra for the trip over the network. Point it at a real GPU server and the numbers get much smaller, but the pattern stays the same: the first word waits for reading, and then the words arrive at a steady beat.

Five easy mistakes when timing a chatbot

Starting the stopwatch at the wrong place. Start it when the user presses Enter, not when the server starts working. Waiting in line is part of the wait.
Counting chunks instead of tokens. Some servers bundle a few tokens into one message. Check how many messages you got against how many tokens the server says it made.
Only looking at the average gap. One 2-second freeze in a long answer barely changes the average, but everybody notices it. Always look at the slow end too (p99).
Forgetting to warm up. The very first request is slow because the engine is still setting up. Don't count it.
Mixing up "fast for me" and "fast for everyone." Per-person speed and whole-factory speed are different numbers. Keep them in separate columns.

Cheat sheet

Machine	Real name	What it does	Speed
Chopper	Tokenize	Text → token numbers	Practically instant
Reader	Prefill	Reads the whole question at once, fills the notebook	Cheap per token; total grows with question length
Notebook	KV cache	Notes on every token so far	Grows by one line per token
Typewriter	Decode	One token per full trip through the model	Expensive per token; the slow part
Delivery Tape	Stream	Sends each token to you immediately	Makes it feel fast

Stopwatch	Measures	Mostly decided by
TTFT	Wait for the first word	Length of your question
ITL	Gap between words	Model size + size of the notebook
Throughput	Tokens per second	Per person, or for everyone at once

What's next

Next time, we open up the Reader-Typewriter machine and look inside. Part 2: A Forward Pass for One Token: Where the FLOPs Go. We'll follow one token through a single layer and count the work done by each part, and find out why almost all of it is just multiplying big tables of numbers.

Every code block in this post was actually run, and every number above comes from those runs. The tiny factory used random weights on a 2-core CPU, so it shows the pattern, not the speed of any real GPU.

Anatomy of a Decode-Step Stall: Where p99 Inter-Token Latency Really Comes From

Sachin Kr. Rajput — Thu, 08 Oct 2026 10:23:04 +0000

Correction (8 Oct 2026): A reader found two bugs in the attribution code in the first version of this post. (1) Preemption was double-counted: while a request was evicted, the time was charged to preemption and to whatever the engine ran during it, so a 40 ms gap could be charged 60 ms. (2) GC inflated graph-miss time: a GC pause inside an eager decode step was removed from the step's share of the gap but not from the step's total before comparing it with the baseline, inventing graph overhead that wasn't there. Both bugs are fixed below. The engine now records eviction and readmission explicitly, every gap's attribution now sums exactly to the measured gap, and every number comes from a new set of 35 longer runs. The bugs did not touch the measured ITL percentiles; they did distort the breakdown by cause, and the biggest change is that graph misses own far less of the tail than v1 claimed. Thank you for the careful review.

The One-Line Summary: On a CPU demonstrator engine I built to make every stall cause real and measurable, 98% of a typical inter-token gap was the decode forward pass, but at or above p99 the forward pass was only 9–29% of the gap. The rest was preemption, prefill, GC and tokenization landing in the same window, and fixing those one cause at a time took the median p99.9 across five 60-second runs from 950 ms to 50 ms.

What this is, and isn't: every number here comes from a CPU demonstrator: a continuous-batching engine in Python and NumPy on a 2-core VM, not a GPU serving stack. Its "CUDA graph" path is emulated: it pads to a captured batch size and falls back above the largest, like the real thing, but the speed difference it measures is batched vs per-sequence NumPy, not graph replay vs kernel launches. Treat the percentiles and bands as observations from this workload. What transfers is the method: tag each step, partition each gap, reconcile.

The Parable of the Harbour Line Logbook

The Harbour Line promises a train every four minutes. The audit says it delivers: the average gap between departures at Quay Street is 4.2 minutes, and the trains themselves cover the line in exactly the time the timetable says, every time, to the second.

The complaint letters say something else. Waited twenty-five minutes on the platform. Twice this week.

Both are true. The new dispatcher, Odile, is the first person to take that seriously.

The Stopwatch on the Train

Her predecessor had investigated the letters by riding the trains with a stopwatch.

THE STOPWATCH ON THE TRAIN
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
  Quay St -> Mill Rd     3 min 58 s
  Quay St -> Mill Rd     4 min 01 s
  Quay St -> Mill Rd     3 min 59 s
  Quay St -> Mill Rd     4 min 00 s
  ...  (two hundred rides)

  verdict: the trains are fine.

✗ Correct, and useless. Nobody on the platform
  was complaining about the ride.

The trains were fine. The question the passengers were asking was not "how fast is the train?" but "why didn't one come?" — and you cannot answer that by timing the train.

The Logbook

Odile does something different. For every departure, she writes down not how long the train took, but everything else that happened on the line since the previous departure.

After a month, the late departures sort themselves into five piles.

ODILE'S LOGBOOK — LATE DEPARTURES, BY WHAT ELSE HAPPENED
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
  THE FREIGHT SLOT       a freight train was let onto the line between
                         two passenger trains. Long freight, long gap.

  THE ODD-LENGTH TRAIN   the platform doors are programmed for 4, 6 and
                         8 cars. A 9-car train had to be boarded with the
                         doors worked by hand, one car at a time.

  THE SIDING SHUFFLE     the yard was full, so a passenger train was
                         pulled onto a siding to make room. Its riders
                         sat there until a slot opened, then it was
                         shunted back. Everyone else waited for the shunt.

  THE SWEEP              the depot crew walked the whole line checking
                         every locker — including forty years of archived
                         paperwork nobody ever opens. The line stops
                         while they walk.

  THE GROUP BOOKING      the gate agent, who also waves trains out,
                         got a tour group of three hundred and sold all
                         their tickets before waving the next train out.

Not one of the five piles is about the train.

One rule keeps the logbook honest. A rider who spent twenty minutes on the siding is entered under siding, once — even though freight happened to roll past while they sat there. Odile notes the freight in the margin, because it explains why the siding took so long, but she never counts those twenty minutes twice.

Five Piles, Five Fixes

What makes the logbook useful is that every pile has its own remedy, and none of the remedies help the others.

THE FIX FOR EACH PILE
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
  freight slot        cut freight into short sections that fit
                      between passenger trains
  odd-length train    program the doors for every length you run
  siding shuffle      don't let a new train onto the line unless
                      the yard has room to spare
  sweep               seal the archive room; stop filing a paper
                      form for every single passenger
  group booking       sell tickets in the street-level office,
                      not at the gate

  Shortening freight does nothing for the sweep.
  Sealing the archive does nothing for the siding.

Odile also notices that the piles live at different heights. The odd-length trains turn a four-minute wait into six. Freight turns it into ten or fifteen. The sweep and the group booking are rare, but when they happen it is half an hour. And the siding shuffle is the worst of all — for the riders on the shunted train, the wait can be longer than the trip.

"The trains were never the problem. The problem was everything I let happen between them."

Why It Works

The average measures the train. The bad days measure the line. Once Odile stopped asking how long a departure took and started asking what shared it, every late departure came with its own explanation attached — and a different person to go and fix it.

That is the whole method. The rest of this post does it to a machine that writes one word at a time, for hundreds of people at once.

What Is a Decode-Step Stall?

A serving engine generating text runs a loop. Each iteration — a step — picks a set of requests, runs one forward pass, and emits one new token for every request that is in its decode phase. The inter-token latency (ITL) of a request is the gap between two of its consecutive tokens:

ITL_{i} = t_{i} - t_{i - 1}

Where the timestamp is taken matters. In this post, $t_{i}$ is taken inside the engine loop, right after the decode forward pass returns. All tokens of a step share it. That is the engine-emission boundary. It is not when the frontend streams the token or when the client receives it — those add serialization and network time, and in real servers one streamed message can carry more than one token (speculative decoding, for instance), so a client-side "gap" can span several tokens.

In a quiet engine, the gap is one step, and one step is mostly the decode forward pass. But a step is shared: continuous batching means the step that produces your next token also does whatever else the scheduler decided this iteration needed. So every gap can be partitioned exclusively — each millisecond goes to exactly one label:

ITL_{i} = T_{fwd} + T_{prefill} + T_{preempt} + T_{graph} + T_{gc} + T_{tok} + T_{residual}

where $T_{fwd}$ is the decode forward pass (on a GPU, the kernels), and $T_{residual}$ is whatever no instrument covered. A decode-step stall is a gap where one of the middle five terms is not zero. This is the serving-engine version of what Dean and Barroso called the tail at scale [9]: rare events that barely move the average but decide the experience.

The claim I tested: the median is mostly the first term; the tail is mostly the middle five. In this demonstrator that held. Whether it holds on your system is exactly what the instrument is for — on fast GPUs, CPU-side overhead can become a large share of every step, not just the tail; vLLM's V1 redesign was motivated by exactly that [7].

The Five Stalls

1. A Prefill Chunk Admitted Into the Batch

A new request needs its whole prompt processed before it can decode — prefill. Modern engines chop prefill into chunks and piggyback the chunks onto decode steps, the idea Sarathi-Serve called stall-free scheduling [1]. vLLM V1 has chunked prefill always on and schedules all pending decodes first, then fills the remaining token budget with prefill [2].

The decodes still wait for the chunk. The knob is the per-step token budget (max_num_batched_tokens in vLLM): vLLM's tuning guide says smaller values give better ITL and larger ones better time-to-first-token [2].

2. A Preemption

The KV cache is a fixed pool. When a running request needs one more block and none is free, the scheduler evicts a victim — either dropping its cache to recompute later, or copying it out to host memory and back (swap). vLLM V1 defaults to recompute [2]. My demonstrator swaps, because the copy is a cleaner thing to time; the victim's wait is the same kind of stall either way, but the cost others pay (a copy here, a re-prefill under recompute) differs.

Preemption has two kinds of victim: everyone in the step pays for the eviction work, and the evicted request emits nothing at all until it is readmitted.

3. A CUDA-Graph Miss on an Unusual Batch Shape

Engines capture CUDA graphs for a fixed set of batch sizes so a decode step is one replay instead of many kernel launches. A batch is padded up to the nearest captured size; vLLM's dispatcher returns CUDAGraphMode.NONE — run eagerly — when the token count exceeds the largest captured size [3]. Actual coverage also depends on the graph mode (full vs piecewise) and on the batch's composition. Recent vLLM versions tie the default capture range to max_num_seqs [4], which is why this tends to bite when someone raises max_num_seqs or trims the capture list.

4. Python Garbage Collection in the Scheduler

The scheduler is Python. Python's cyclic collector periodically walks every tracked container object; a full (generation-2) pass scales with the live heap, and it runs on whichever thread allocated the object that tipped the counter — here, the engine loop.

vLLM freezes the GC heap after startup so static objects are never rescanned, disables collection during graph capture because a GC cycle there can invalidate the graph, and ships VLLM_GC_DEBUG to log collection times [5]. Python 3.14.0 shipped an incremental collector that cut maximum pauses by an order of magnitude on large heaps, and 3.14.5 reverted it to the 3.13 generational collector after production reports of memory pressure [6].

5. A Tokenizer Hiccup

Tokenizing a 300,000-character pasted document is real CPU work. If it happens on the thread that drives the step loop, every in-flight request waits for it. vLLM V1 moved tokenization and detokenization out of the engine-core process so that work overlaps with the core loop [7].

Instrumenting One Engine

The demonstrator: a continuous-batching loop in Python with a random-weight 2-layer transformer in NumPy (d=128); a block-based KV pool (1,000 blocks of 16 tokens); chunked prefill with decode-first scheduling; swap-based preemption of the newest request, with eviction and readmission timestamps recorded explicitly; an emulated "captured graph" path for batch sizes {1, 2, 4, 8, 16, 24, 32} (one padded, batched NumPy call) with a per-sequence fallback above that; a real BPE tokenizer (HuggingFace tokenizers, trained in-process); and Python's real garbage collector. What I call "forward pass" below is CPU forward-pass time in this engine.

The workload: 60 seconds per run, about 6.5 requests/s on average with a 15-second swell (±60%); prompts around 1,000 characters; 1% of requests paste a 150k–400k-character document (truncated to 2,048 tokens after tokenization); outputs averaging about 300 tokens; 400,000 long-lived startup objects standing in for a real server's configs and registries. Each run: 385–422 requests and 112,000–116,000 inter-token gaps. Five runs (seeds 1–5) per configuration, 35 runs in all.

Why 6.5 req/s: the reruns landed on a slower VM than v1 (in the same microbenchmark, the forward pass at batch 32 took 12.8 ms vs 5.6 ms), so I recalibrated the arrival rate to keep the same near-capacity regime: queues drain, the KV pool is under pressure, and batches occasionally cross 32. At the original 10 req/s this VM was overloaded — see the sidebar below. Absolute numbers differ from v1 for that reason.

The Step Tag

Every step records what it did, as (segment, start, end) windows. Trimmed from the engine (the full source is at the end):

def step(self, now0):
    tag = {"seg": [], "t0": now0, ...}

    # 1. intake: arrivals become waiting requests (tokenizer runs here)
    while self.trace and self.trace[0].arrival <= now0:
        r = self.trace.pop(0)
        t = time.perf_counter()
        r.prompt = tokenize(r.text)[: cfg.max_model_len]
        tag["seg"].append(("tok", t, time.perf_counter()))
        self.waiting.append(r)

    # 2. swap-in: readmission is an explicit, timestamped event
    ...
    r.evictions[-1][1] = time.perf_counter()     # readmitted
    ...
    # 3. reserve a KV slot for every decode; may evict the newest request
    ...
    # 4-5. leftover token budget goes to prefill chunks
    ...
    tag["seg"].append(("prefill", t, time.perf_counter()))

    # 6. decode: emulated captured path, or per-sequence fallback above every captured size
    i = bisect.bisect_left(cfg.capture_sizes, B)
    graph = cfg.capture_sizes[i] if i < len(cfg.capture_sizes) else None
    ...
    tag["seg"].append(("decode", t, time.perf_counter()))
    now = time.perf_counter()                    # the emission timestamp for every token in this step

and eviction is stamped at the moment it happens:

def _preempt(self, r):
    ...
    r.evictions.append([time.perf_counter(), None])  # evicted (before the copy-out)

GC is timed with the standard library, two clock reads per collection:

gc.callbacks.append(self._gc_cb)

def _gc_cb(self, phase, info):
    if phase == "start": self._gc_t0 = time.perf_counter()
    else: self.gc_log.append((self._gc_t0, time.perf_counter(), info["generation"]))

The run dumps every step, segment, GC pause, emission, eviction and request to a raw log file. Every number below is recomputed from those raw logs by a separate analysis script.

Partitioning a Gap — and What v1 Got Wrong

The rule now is a strict partition of the window between two emissions, using the request's actual state:

While the request is evicted (explicit eviction → readmission events), every millisecond is preemption. Whatever the engine ran meanwhile is recorded separately as context, never added to the gap.
Otherwise, GC pauses take precedence over whatever segment they landed in; then each segment's GC-free overlap goes to its label.
Decode time is split into graph-miss excess and forward pass using the step's GC-free decode time against a GC-free baseline fit on captured-path steps, then applied proportionally to the part of the segment inside the gap.
Whatever no segment covers is residual. By construction the parts sum to the gap, and the analysis checks that they do.

Verbatim from attribution.py:

    def _graph_fraction(self, i, st):
        """Share of this step's (GC-free) decode time that is graph-miss excess."""
        if st["graph"] != "eager": return 0.0
        if i not in self._dnet: self._dnet[i] = decode_net_ms(st, self.gci)
        full = self._dnet[i]
        if full <= 0: return 0.0
        base = self.slope * st["n_decode"] + self.icpt
        return max(0.0, full - base) / full

    def _charge(self, x, y, acc):
        """Charge window [x, y] (s) segment-by-segment into acc (ms). Returns covered ms."""
        gc = self.gci.overlap(x, y)
        acc["gc"] += gc * 1e3
        covered = gc
        for i, st in self._steps_overlapping(x, y):
            for name, a, b in st["seg"]:
                a, b = max(a, x), min(b, y)
                if b <= a: continue
                net = (b - a) - self.gci.overlap(a, b)
                covered += net
                lab = SEG_LABEL[name]
                if lab == "decode":
                    f = self._graph_fraction(i, st)
                    acc["graph"] += net * f * 1e3
                    acc["kernel"] += net * (1 - f) * 1e3
                else:
                    acc[lab] += net * 1e3
        return covered * 1e3

    def gap(self, t0, t1, evictions=()):
        """Partition one gap. Returns (parts, context): parts sum exactly to the gap."""
        parts = {k: 0.0 for k in LABELS}
        context = {k: 0.0 for k in LABELS}
        ev = [(a, b if b is not None else t1) for a, b in evictions if a < t1 and (b is None or b > t0)]
        for a, b in ev:                                   # evicted: all of it is preemption
            a, b = max(a, t0), min(b, t1)
            if b > a:
                parts["preempt"] += (b - a) * 1e3
                cov = self._charge(a, b, context)         # what the engine was doing meanwhile
                context["residual"] += (b - a) * 1e3 - cov
        for a, b in _subtract(t0, t1, ev):                # live: charge what actually ran
            cov = self._charge(a, b, parts)
            parts["residual"] += (b - a) * 1e3 - cov
        return parts, context

(kernel is the label name in code; in this demonstrator it means CPU forward-pass time.)

Here are the reviewer's two controlled examples, run through the v1 function unchanged:

EX1 gap 40.0 ms -> prefill 30.0 + preempt 30.0 = 60.0 ms attributed
EX2 gap 30.0 ms -> gc 20.00  graph 6.67   (correct graph excess: 0)

v1 charged 60 ms inside a 40 ms gap, and invented 6.67 ms of graph overhead out of a GC pause. The same two cases, plus four more, as unit tests against the corrected partitioner:

PASS test_reviewer_example_1_eviction_not_double_counted
PASS test_reviewer_example_2_gc_does_not_inflate_graph_excess
PASS test_real_graph_excess_survives_gc_removal
PASS test_partial_overlap_with_gap_is_proportional
PASS test_gc_after_emission_and_residual_reconcile
PASS test_mixed_step_prefill_tok_swap_and_decode

In the first case the corrected partition gives 30 ms preemption plus 10 ms forward pass, which is the 40 ms gap; the 30 ms of other requests' prefill is recorded as context. In the second it gives 20 ms GC, 10 ms forward pass and 0 ms graph excess.

Across all 35 runs and about 4 million gaps, the largest difference between a gap and the sum of its parts was 5.7 × 10⁻¹⁴ ms. And no gap skipped engine steps without a recorded eviction covering it. v1 inferred preemption from skipped steps; that check now verifies the explicit events agree with that inference rather than relying on it.

The p99, Decomposed

All results: five runs per configuration. Each run is the unit of evidence here, not each gap. A run has ~114,000 gaps, but they are far from independent: one GC pause stalls every request in the batch at once, and one preemption burst evicts several. So I report each run separately, or as a median with the min–max across the five.

Median Is the Forward Pass. Tail Is Not.

For each token, the share of its gap that was the decode forward pass (baseline configuration):

FORWARD-PASS SHARE OF A GAP (median), baseline: clean vs tail tokens
  run 1: clean 98%   tail 23%
  run 2: clean 98%   tail 29%
  run 3: clean 98%   tail 11%
  run 4: clean 98%   tail 18%
  run 5: clean 98%   tail 9%

In an ordinary gap, the forward pass is 98% of the wait in all five runs. In a gap at or above that run's p99 it is 9–29%. In this engine, making the forward pass faster improves typical tokens almost one-for-one and the tail far less.

Where the Tail Time Went

The share of all tail time (tokens with ITL ≥ p99) by label. With the corrected attribution, these are additive and each row sums to 100%:

BASELINE: WHERE TAIL TIME WENT (tokens with ITL >= p99), additive, sums to 100%
run     p99  prefill  preempt    graph       gc      tok   kernel residual
1      30.9    23.0%    57.2%     0.4%     4.4%     3.6%    11.0%     0.4%
2      34.1    45.2%    12.2%     0.7%     7.6%     9.9%    23.6%     0.8%
3      50.3     3.5%    93.2%     0.2%     0.8%     1.2%     1.0%     0.0%
4      47.0     5.3%    90.0%     0.1%     1.6%     1.0%     1.9%     0.1%
5      62.8    12.4%    77.1%     0.2%     2.1%     5.3%     2.8%     0.2%

(kernel = CPU forward pass.) Three things stand out:

Preemption holds most of the tail time in four of five runs (57–93%). Run 2 drew a lighter trace (44 preemptions instead of 142–318), and its tail is mostly prefill.
Count and time disagree. By count, prefill was the most common primary cause of a tail token in every run (649–880 of ~1,150 tail tokens). By time, preemption dominates: fewer tokens, but each one enormous.
Graph misses own almost none of the tail: 0.1–0.7%. This is the largest correction from v1, which reported up to 7%. Part of that was the GC bug; part was the old labeling. Graph misses here are frequent and small: 1.7% of all gap time (median; 1.3–4.1% per run), which makes them mostly a p90–p99 effect rather than a tail one.

The Histogram

Every gap from run 4 (the run whose p99, 47.0 ms, is the median of the five), binned, with each bin's tokens split by primary cause: the largest cause, if it is ≥1 ms and ≥25% of the gap. · means no cause met that bar, which in practice means the forward pass. P prefill, C graph miss (emulated), S preemption, G GC, T tokenizer.

ITL HISTOGRAM BY PRIMARY CAUSE, baseline, run 4   (p99 = 47.0 ms)
  ITL ms     tokens  share of the bin's tokens by primary cause
     0-5      13760  ····································
     5-10     17747  ····································
    10-15     37670  ····································
    15-20     35524  ···································P
    20-30      8907  ······················PPPPPPPPCCCCCC
    30-50      1460  ···PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPCC
    50-80       589  ··PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPP
    80-120      111  PPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPPS
   120-200       58  PPPPPPPPPPPPPPPPPPPPSTTTTTTTTTTTTTT
   200-400      170  PPPPPPPPPPPPPPSSGGGGGGGGGGGGGGGGTTTT
   400-800       11  SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS
   800-inf      171  SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS

Across the five baseline runs the causes overlap but have different centres of mass. Graph misses sat at 15–120 ms. GC, apart from a few small pauses, sat at 120–400 ms. The tokenizer spread from 15 to 400 ms, and prefill appeared everywhere from 0 to 800 ms. Preemption started around 30 ms and was the only cause above 800 ms. This run is tidier than most: here, every gap above 400 ms is preemption, but in run 3, 23 of the 53 gaps in the 400–800 ms bin were prefill. So read these bands as loose observations about this engine and workload, not as a law. A GPU engine with a faster forward pass and a different prefill cost will put them somewhere else.

The Fix Ladder

One cause at a time, each rung keeping every fix above it. Median and [min–max] across five runs, in ms:

FIX LADDER  (5 runs x 60 s each; median [min-max] across runs; ms)
rung                            p50               p99               p99.9               max
baseline            13.7 [10.5-16.0]  47.0 [30.9-62.8]      950 [149-7395]  9466 [895-12764]
+gc.freeze          14.3 [7.9-14.7]  49.0 [27.0-60.7]       581 [112-953]   3805 [368-6100]
+compact history    11.3 [6.0-14.7]  40.4 [21.6-49.2]        258 [56-827]   2268 [451-5165]
+graphs to 48       12.3 [10.6-14.0]  37.3 [31.9-50.6]       234 [122-447]  2971 [1295-3824]
+budget 256         12.6 [9.4-14.8]  29.5 [25.2-32.7]        149 [54-542]   1933 [612-2646]
+tokenize off-loop  12.5 [10.8-15.5]  29.6 [26.3-39.0]        133 [53-967]   2698 [889-4210]
+admit watermark    12.0 [10.4-14.1]  27.8 [24.2-32.5]          50 [31-70]      108 [39-651]

End to end, on medians: p99 47.0 → 27.8 ms, p99.9 950 → 50 ms, worst gap 9,466 → 108 ms. The median ITL's range barely moved (10.5–16.0 → 10.4–14.1 ms). Read adjacent rungs with the ranges in view: with five runs, several single-rung steps overlap and are suggestive rather than established. At the ends of the ladder, the p99.9 and max ranges don't overlap at all; p99 overlaps only between 30.9 and 32.5 ms, because run 1's baseline drew light traffic.

And the tail's composition, rung by rung (median across runs). The shares sum to 100%, so when one cause shrinks, the others' shares grow even if their absolute time didn't change. Read it for which column each fix drives toward zero:

MEDIAN TAIL-TIME SHARE PER RUNG (5 runs)
                    prefill  preempt  graph    gc     tok   fwd    residual
baseline             12.4%    77.1%   0.2%   2.1%   3.6%   2.8%   0.2%
+gc.freeze           14.0%    70.6%   0.1%   1.6%   5.7%   5.7%   0.2%
+compact history     23.1%    59.0%   0.2%   0.0%   7.4%   8.3%   0.3%
+graphs to 48        30.5%    49.0%   0.0%   0.0%   6.6%   9.8%   0.5%
+budget 256          12.0%    56.0%   0.0%   0.0%   8.9%  25.7%   0.3%
+tokenize off-loop    8.5%    63.2%   0.0%   0.0%   0.0%  27.4%   0.8%
+admit watermark     34.7%     0.2%   0.0%   0.0%   0.0%  63.2%   1.9%

What each rung did:

gc.freeze() after startup. Generation-2 pauses during serving dropped from 137–242 ms to 30–90 ms, because the startup heap stopped being rescanned. The one-off gc.collect() before freezing costs 112–130 ms on this rung (107–167 ms across all frozen rungs), at startup, before any request. Median p99.9 950 → 581 ms; p99 did not move.

Compact history. The engine kept a Python object per emitted token: a record plus a dict of five top-k logprob objects, the way a metrics store or streaming buffer might. Moving that into a flat array("d"), which the collector can't see, cut generation-0 collections from ~380–400 per run to 20–44, and no generation-2 collection ran during serving in any of the five runs. GC's share of the tail went to 0.

Capture graphs up to max_num_seqs (emulated). Fallback steps went from 0–983 per run on the previous rung (207–1,087 at baseline) to zero. Its share of tail time was already under 1%; this fix mostly helps the p90–p99 body.

Token budget 2,048 → 256. The largest single p99 step: median 37.3 → 29.5 ms, with ranges that overlap only between 31.9 and 32.7 ms.

Tokenize off the engine loop. The tokenizer's share of the tail went to zero, but the percentiles didn't clearly move (p99 29.5 → 29.6, p99.9 149 → 133, ranges overlapping). On a 2-core VM, the tokenizer worker process competes with the engine for the same cores. It moves the stall from blocking to contention without removing the work. In v1's runs, this rung made p99 slightly worse. In production, the frontend should run on cores the engine core doesn't need, which is what vLLM V1's process split is for [7].

Admission watermark. Don't admit a new request unless 10% of KV blocks would stay free. Preemptions fell from 23–346 per run on the previous rung to 0–7, and the median worst gap from 2,698 to 108 ms. This time the cost showed up clearly: TTFT p50 rose in all five runs against the previous rung (for example 164 → 931 ms in run 4, 4.6 → 5.6 s in run 3, which was transiently overloaded at its traffic peak). The waiting moved from the middle of people's answers to before they start. That is usually the right trade, but it is a trade.

After the last rung, the tail is 63% forward pass and 35% prefill. Same run as the histogram above:

ITL HISTOGRAM BY PRIMARY CAUSE, +admit watermark, run 4   (p99 = 27.1 ms)
  ITL ms     tokens  share of the bin's tokens by primary cause
     0-5      19076  ····································
     5-10     18967  ···································P
    10-15     52473  ···································P
    15-20     18896  ····························PPPPPPPP
    20-30      6059  ················PPPPPPPPPPPPPPPPPPPP
    30-50       578  ·········PPPPPPPPPPPPPPPPPPPPPPPPPPP
    50-80        57  ····················PPPPPPPPPPPPPPPP
    80-120       67  ···································S
   120-200        1  SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS
   200-400        2  SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS
   400-800        2  SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS
   800-inf        0

One thing I could not explain. Of the 67 gaps in the 80–120 ms bin, 66 are untagged. They average 103.5 ms, and 92.7 ms of that is forward pass: decode steps that ran five or more times slower than usual. None of the five instrumented causes covers it. CPU contention from the tokenizer worker or noise from the VM are the obvious suspects, but my instrument can't see the OS scheduler. It's the sixth pile.

Sidebar: When the Engine Is Overloaded

My first rerun used the original 10 req/s on this slower VM, which turned out to be past capacity: TTFT p50 of 2–17 seconds, and a queue that never drained. The same corrected analysis on those 35 runs: preemption held 67–96% of baseline tail time, the next five rungs moved median p99.9 around inside the noise (593–3,177 ms), and only the admission watermark changed the tail, from 1,022 to 56 ms median p99.9. Under overload, the tail is a capacity problem first, and admission control is the fix that matters.

Each Cause, Its Own Fix

Cause	Where it sat (5 baseline runs, loosely)	How you see it	Fix	What the fix costs
Prefill chunk	anywhere, 0–800 ms	prefill segments in the step	smaller token budget; disaggregate prefill [8]	TTFT at higher prefill load
Graph miss (emulated)	15–120 ms, <1% of tail time	fallback / padded size per step	capture sizes up to `max_num_seqs`	graph memory, startup time
Preemption	30 ms up; the only cause above 800 ms	explicit eviction → readmission	admission watermark; more KV; smaller `max_num_seqs`	TTFT (rose in 5/5 runs)
GC	mostly 120–400 ms	`gc.callbacks`	`gc.freeze()`; no per-token Python objects	negligible here
Tokenizer	15–400 ms	tokenize segments on the loop	separate process with its own cores	a process, IPC, CPU

In this experiment, each targeted cause's share of the tail went to zero (or near it) at its rung. GC took two rungs, gc.freeze() then compact history. The other columns shifted because shares are relative, not because the fixes acted on them.

Doing This on a Real GPU Engine

I haven't run this instrument against vLLM or SGLang in this post; that's the next one. Here is what changes when you do:

WHAT TO RECORD, AND THE TRAPS
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
  forward pass       CUDA launches are ASYNCHRONOUS [10]: a host timer around
                     the forward call measures the launch, not the work.
                     Bracket it with CUDA events, or use Nsight Systems /
                     the PyTorch profiler, and correlate with host stamps.
  prefill vs decode  a mixed batch is ONE forward pass. There is no host
                     "prefill segment" to time; split its cost by kernel
                     (profiler) or attribute by the step's token mix.
  preemption         record eviction and readmission per request from the
                     scheduler. Do not infer it from skipped steps.
  graph mode         what the CUDA-graph dispatcher returned for the step
                     (FULL / PIECEWISE / NONE) and the padded size [3].
  GC                 gc.callbacks in the engine-core process, or
                     VLLM_GC_DEBUG=1 [5].
  tokenizer          measured in the frontend process, plus the time the
                     request then waits before the engine core sees it.
  timestamps         say which boundary: engine emission, frontend stream,
                     or client receipt. A streamed chunk can hold several
                     tokens.
  critical path      with async scheduling, CPU work for step N+1 can
                     overlap GPU work for step N. Charge a gap only with
                     what was on the path between its two emissions.

Then join everything on time windows and partition each gap the same way, with the same reconciliation check.

Reproducing This

ENVIRONMENT
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
  CPU       Intel Xeon @ 2.80 GHz, 2 cores, 1 thread/core (cloud VM)
  Python    3.13.16
  NumPy     2.5.3   (OpenBLAS 0.3.34)
  tokenizers 0.23.2
  threads   OPENBLAS_NUM_THREADS=1, OMP_NUM_THREADS=1 (set before import)

PROCEDURE
  for seed in 1..5, for each rung:
      ITL_RATE=6.5 python3 run.py "<rung>" <seed> 60   -> raw/<rung>__s<seed>.npz
  python3 analyze.py raw/*.npz > results.jsonl
  python3 test_attribution.py

  arrival traces are deterministic in the seed (make_trace) and are also
  saved inside every raw file, with all step, segment, GC, emission,
  eviction and request timestamps.

The complete engine, the partitioner, the analysis and the runner are below, exactly as they ran. I also have the 35 raw timing logs (about 46 MB); ask in the comments if you want them.

Five 60-second runs per configuration is preliminary evidence, not a benchmark. The p99.9 and max columns rest on about a hundred gaps and a single gap per run respectively, and they are the noisiest numbers here. Longer and more repeated runs would tighten them.

Quick Reference Card

DECODE-STEP STALLS: CHEAT SHEET
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
THE PARTITION (exclusive; must sum to the gap)
  ITL = fwd + prefill + preempt + graph + gc + tok + residual
  evicted time      -> preempt (engine activity = context, not added)
  GC                -> gc, carved out of its segment
  graph excess      -> GC-free decode time minus GC-free baseline
  uncovered         -> residual
  check: |gap - sum(parts)| ~ 0 for every gap

MEASURED (CPU demonstrator, 5 runs x 60 s, median [range], ms)
                    p50                p99                p99.9
  baseline     13.7 [10.5-16.0]   47.0 [30.9-62.8]   950 [149-7395]
  all fixes    12.0 [10.4-14.1]   27.8 [24.2-32.5]    50 [31-70]
  forward pass = 98% of a typical gap, 9-29% of a tail gap

TAIL TIME, BASELINE (these runs)
  preemption 12-93% | prefill 4-45% | gc 1-8% | tok 1-10% | graph <1%
  by COUNT prefill leads; by TIME preemption leads

FIXES
  prefill   smaller token budget
  graph     capture up to max_num_seqs
  preempt   admission watermark (TTFT rose in 5/5 runs)
  gc        gc.freeze() + no per-token Python objects
  tok       separate process WITH its own cores

Key Takeaways

In this demonstrator, the median ITL was the forward pass and the tail was everything else: 98% vs 9–29% forward pass. Measure it on your own system before assuming it; on fast GPUs, CPU overhead can show up in every step.
Partition, don't sum. Each millisecond of a gap gets exactly one label, evicted time goes to preemption, and the parts must reconcile to the gap. My first version broke both rules and over-counted.
Record preemption explicitly. Eviction and readmission are events. Inferring them from missing steps only works when every running request decodes every step.
Count and time disagree. Prefill crossed p99 most often; preemption held the most tail time. Report both.
Each fix drove its own cause's column to zero; no single fix touched them all. The shares are relative, so the other columns shifted, but no one rung took out more than one cause (GC needed two).
Fixes have costs that land elsewhere. The watermark raised TTFT in every run, and off-loop tokenization on two cores bought nothing.
Treat runs, not gaps, as your sample size. One GC pause stalls a whole batch at once.

The One-Sentence Summary

On a CPU demonstrator, inter-token latency at the median was the forward pass and at p99 was mostly scheduling and runtime, so I tagged every step, partitioned every gap exclusively by cause, checked that the parts summed to the gap, and fixed one cause at a time, taking median p99.9 across five runs from 950 ms to 50 ms without touching the forward pass.

What's Next?

This post starts a new series, Inference Engineering — what actually happens between a request arriving and a token leaving, measured rather than described.

The same instrument on vLLM: CUDA events, explicit scheduler events, and how the bands move when the forward pass is 10× faster than the scheduler.
Chunked prefill vs disaggregated prefill: when shrinking the token budget stops being enough.
Preemption, measured: recompute vs swap, and the queueing math behind an admission watermark.
The CPU overhead floor: when the scheduler, not the GPU, sets your ITL.

Follow me for the next article in the Inference Engineering series!

Let's Connect!

If the Harbour Line logbook made p99 make sense, drop a heart!

Questions? Ask in the comments — I read and respond to every one.

Check my partition. The first version of this post had two attribution bugs, and a reader found both by building small controlled cases and running my code on them. That's the best review a measurement post can get. If you can construct a gap my partitioner charges wrongly, I want to see it. 🚉

References

[1] A. Agrawal, N. Kedia, A. Panwar, J. Mohan, N. Kwatra, B. Gulavani, A. Tumanov and R. Ramjee, Taming Throughput-Latency Tradeoff in LLM Inference with Sarathi-Serve, OSDI 2024.

[2] vLLM documentation, Optimization and Tuning — preemption mode, chunked prefill and max_num_batched_tokens.

[3] vLLM API reference, vllm.v1.cudagraph_dispatcher.

[4] vLLM forum, Why is CUDA graph capture sizes limited by max_num_seqs.

[5] vLLM API reference, vllm.utils.gc_utils — freeze_gc_heap, freeze_gc_for_cudagraph_capture, VLLM_GC_DEBUG.

[6] Python documentation, What's New in Python 3.14 — the incremental GC and its reversion in 3.14.5.

[7] vLLM team, vLLM V1: A Major Upgrade to vLLM's Core Architecture, 2025.

[8] Y. Zhong, S. Liu, J. Chen, J. Hu, Y. Zhu, X. Liu, X. Jin and H. Zhang, DistServe: Disaggregating Prefill and Decoding for Goodput-optimized Large Language Model Serving, OSDI 2024.

[9] J. Dean and L. A. Barroso, The Tail at Scale, Communications of the ACM, 2013.

[10] M. Harris, How to Implement Performance Metrics in CUDA C/C++, NVIDIA Technical Blog — why host timers measure kernel launch rather than execution, and CUDA events.

Full Source

engine.py — the demonstrator engine

"""nano-engine: a CPU continuous-batching LLM server loop, instrumented so every
decode step records *what else happened in that step*.

The model is a random-weight transformer stack in numpy. The point is not the
tokens it emits; it is that every mechanism that stretches a decode step is the
real mechanism, executed for real:

  prefill   chunked prefill admitted into the same step as the decodes
  preempt   KV pool exhausted -> victim's KV copied out (swap) and back in
  graph     batch sizes with no "captured" fast path run op-by-op (eager)
  gc        Python's cyclic GC, timed with gc.callbacks
  tok       incoming prompts tokenized inline on the engine thread
"""
import os
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
os.environ.setdefault("OMP_NUM_THREADS", "1")
import gc, time, json, random, bisect, sys
from array import array
from dataclasses import dataclass, field
from concurrent.futures import ProcessPoolExecutor
import numpy as np

# ---------------------------------------------------------------- tokenizer
_TOK = None
WORDS = ("the of and to in is for on that with as by at from this be are was it an or "
         "kernel batch token cache latency decode prefill scheduler memory throughput "
         "request queue model layer tensor shape graph block swap stream budget").split()

def make_text(rng, n_chars):
    out, n = [], 0
    while n < n_chars:
        w = rng.choice(WORDS)
        if rng.random() < 0.08:                      # numbers and odd tokens
            w = f"{rng.randint(0, 99999)}{rng.choice(['ms', '%', 'x', '.0'])}"
        out.append(w); n += len(w) + 1
    return " ".join(out)

def get_tokenizer():
    global _TOK
    if _TOK is None:
        from tokenizers import Tokenizer, models, trainers, pre_tokenizers
        tok = Tokenizer(models.BPE(unk_token="[UNK]"))
        tok.pre_tokenizer = pre_tokenizers.ByteLevel(add_prefix_space=False)
        rng = random.Random(0)
        corpus = [make_text(rng, 400) for _ in range(3000)]
        tok.train_from_iterator(corpus, trainers.BpeTrainer(vocab_size=4096, show_progress=False))
        _TOK = tok
    return _TOK

def tokenize(text):                                   # also the worker entry point
    return get_tokenizer().encode(text).ids

# ---------------------------------------------------------------- model
D, L, FF = 128, 2, 512

MAXLEN = 2048 + 640
CAUSAL = np.triu(np.full((MAXLEN, MAXLEN), -1e9, np.float32), k=1)   # built once

class Model:
    def __init__(self, seed=0):
        r = np.random.default_rng(seed)
        s = 1 / np.sqrt(D)
        self.Wqkv = [r.standard_normal((D, 3 * D), dtype=np.float32) * s for _ in range(L)]
        self.Wo   = [r.standard_normal((D, D), dtype=np.float32) * s for _ in range(L)]
        self.W1   = [r.standard_normal((D, FF), dtype=np.float32) * s for _ in range(L)]
        self.W2   = [r.standard_normal((FF, D), dtype=np.float32) / np.sqrt(FF) for _ in range(L)]
        self.emb  = r.standard_normal((4096, D), dtype=np.float32) * 0.1

    def _layer(self, l, x, kv_list, positions):
        """x: (n, D) rows; kv_list[i] = (K, V, length) arrays this row attends over."""
        qkv = x @ self.Wqkv[l]
        q, k, v = qkv[:, :D], qkv[:, D:2 * D], qkv[:, 2 * D:]
        att = np.empty_like(q)
        for i, (K, V, end) in enumerate(kv_list):
            K[l, positions[i]] = k[i]; V[l, positions[i]] = v[i]
            Kl, Vl = K[l, :end], V[l, :end]
            sc = Kl @ q[i] * (1 / np.sqrt(D))
            sc = np.exp(sc - sc.max()); att[i] = (sc / sc.sum()) @ Vl
        x = x + att @ self.Wo[l]
        return x + np.maximum(x @ self.W1[l], 0) @ self.W2[l]

    def forward(self, tok_ids, kv_list, positions):
        x = self.emb[np.asarray(tok_ids) % 4096]
        for l in range(L):
            x = self._layer(l, x, kv_list, positions)
        return np.abs(x[:, :8].sum(1) * 1000).astype(np.int64) % 4096

    def forward_prefill(self, ids, K, V, start):
        """One chunk of one prompt: dense causal attention over the prefix."""
        x = self.emb[np.asarray(ids) % 4096]; n = len(ids); end = start + n
        mask = CAUSAL[start:end, :end]
        for l in range(L):
            qkv = x @ self.Wqkv[l]
            K[l, start:end] = qkv[:, D:2 * D]; V[l, start:end] = qkv[:, 2 * D:]
            sc = qkv[:, :D] @ K[l, :end].T
            sc *= 1 / np.sqrt(D); sc += mask
            sc -= sc.max(1, keepdims=True); np.exp(sc, out=sc); sc /= sc.sum(1, keepdims=True)
            x = x + (sc @ V[l, :end]) @ self.Wo[l]
            x = x + np.maximum(x @ self.W1[l], 0) @ self.W2[l]

    def forward_eager(self, tok_ids, kv_list, positions):
        """No captured graph: every sequence is its own small launch."""
        return np.concatenate([self.forward([t], [kv], [p])
                               for t, kv, p in zip(tok_ids, kv_list, positions)])

# ---------------------------------------------------------------- requests
@dataclass(eq=False)
class Logprob:                            # like an OpenAI-style top-k logprob entry
    logprob: float; rank: int; decoded: str = None

@dataclass(eq=False)
class TokenEvent:                         # per-token record kept for streaming/metrics
    req: "Request"; tid: int; t: float; logprobs: dict

@dataclass(eq=False)
class Request:
    rid: int; arrival: float; text: str; max_new: int
    prompt: list = None
    computed: int = 0                     # tokens whose KV is materialised
    out: list = field(default_factory=list)
    emit_t: list = field(default_factory=list)
    emit_step: list = field(default_factory=list)
    K: np.ndarray = None; V: np.ndarray = None
    swapped: tuple = None
    first_token_t: float = None
    evictions: list = field(default_factory=list)   # [t_evicted, t_readmitted] pairs, explicit

    @property
    def total_len(self): return len(self.prompt) + len(self.out)
    def blocks_needed(self, n): return -(-n // BLOCK)

BLOCK = 16

# ---------------------------------------------------------------- engine
@dataclass
class Config:
    name: str = "baseline"
    max_num_seqs: int = 48
    token_budget: int = 2048              # max_num_batched_tokens
    capture_sizes: tuple = (1, 2, 4, 8, 16, 24, 32)
    kv_blocks: int = 1000
    admit_watermark: float = 0.0          # fraction of blocks kept free before admitting
    gc_freeze: bool = False
    async_tokenize: bool = False
    compact_history: bool = False         # per-token records in a flat array, not objects
    max_model_len: int = 2048

class Engine:
    def __init__(self, cfg, model, trace):
        self.cfg, self.m, self.trace = cfg, model, trace
        self.free_blocks = cfg.kv_blocks
        self.waiting, self.running, self.swapped, self.done = [], [], [], []
        self.steps = []                   # one dict per step: the tags
        self.gc_log = []
        self.pool = ProcessPoolExecutor(1, initializer=get_tokenizer) if cfg.async_tokenize else None
        self.pending_tok = []
        self.history = []                 # all TokenEvents: long-lived heap, like a metrics store
        self.flat = array("d")            # the same data, invisible to the GC
        gc.callbacks.append(self._gc_cb)

    def _gc_cb(self, phase, info):
        if phase == "start": self._gc_t0 = time.perf_counter()
        else: self.gc_log.append((self._gc_t0, time.perf_counter(), info["generation"]))

    # ---- KV accounting
    def _alloc(self, r, new_len):
        need = r.blocks_needed(new_len) - r.blocks_needed(r.computed)
        if need > self.free_blocks: return False
        self.free_blocks -= need
        if r.K is None or r.K.shape[1] < new_len:
            cap = min(self.cfg.max_model_len + 1024, max(new_len, 256) * 2)
            K = np.zeros((L, cap, D), np.float32); V = np.zeros_like(K)
            if r.K is not None: K[:, :r.computed] = r.K[:, :r.computed]; V[:, :r.computed] = r.V[:, :r.computed]
            r.K, r.V = K, V
        return True

    def _free(self, r):
        self.free_blocks += r.blocks_needed(r.computed)

    # ---- one engine step
    def step(self, now0):
        tag = {"seg": [], "t0": now0, "tok_ms": 0.0, "prefill_ms": 0.0, "swap_ms": 0.0,
               "decode_ms": 0.0, "n_decode": 0, "n_prefill_tok": 0, "graph": None,
               "preempted": 0, "swapped_in": 0}
        cfg = self.cfg
        # 1. intake: arrivals become waiting requests (tokenizer runs here)
        while self.trace and self.trace[0].arrival <= now0:
            r = self.trace.pop(0)
            if self.pool:
                self.pending_tok.append((r, self.pool.submit(tokenize, r.text)))
            else:
                t = time.perf_counter()
                r.prompt = tokenize(r.text)[: cfg.max_model_len]
                tag["tok_ms"] += (time.perf_counter() - t) * 1e3
                tag["seg"].append(("tok", t, time.perf_counter()))
                self.waiting.append(r)
        for item in [p for p in self.pending_tok if p[1].done()]:
            item[0].prompt = item[1].result()[: cfg.max_model_len]
            self.waiting.append(item[0]); self.pending_tok.remove(item)

        # 2. swap-in preempted requests if room (they go first: FCFS)
        t = time.perf_counter()
        while self.swapped and len(self.running) < cfg.max_num_seqs:
            r = self.swapped[0]
            if r.blocks_needed(r.computed + 1) + 8 > self.free_blocks: break
            self.swapped.pop(0)
            Kh, Vh = r.swapped
            r.K = np.zeros((L, Kh.shape[1] * 2, D), np.float32); r.V = np.zeros_like(r.K)
            r.K[:, :r.computed] = Kh; r.V[:, :r.computed] = Vh
            r.swapped = None
            self.free_blocks -= r.blocks_needed(r.computed)
            r.evictions[-1][1] = time.perf_counter()     # readmitted
            self.running.append(r); tag["swapped_in"] += 1
        tag["swap_ms"] += (time.perf_counter() - t) * 1e3
        tag["seg"].append(("swap", t, time.perf_counter()))

        # 3. decodes first: every running request that has finished prefill
        decodes = [r for r in self.running if r.computed >= len(r.prompt)]
        t = time.perf_counter()
        for r in list(decodes):
            if r not in decodes: continue          # already preempted as someone's victim
            while not self._alloc(r, r.computed + 1) and self.running:
                victim = self.running[-1]         # newest running request
                self._preempt(victim); tag["preempted"] += 1
                if victim in decodes: decodes.remove(victim)
                if victim is r: break
        tag["swap_ms"] += (time.perf_counter() - t) * 1e3
        tag["seg"].append(("swap", t, time.perf_counter()))

        # 4. leftover budget goes to prefill chunks (running partials, then waiting)
        budget = cfg.token_budget - len(decodes)
        chunks = []
        for r in [r for r in self.running if r.computed < len(r.prompt)]:
            n = min(budget, len(r.prompt) - r.computed)
            if n > 0 and self._alloc(r, r.computed + n): chunks.append((r, n)); budget -= n
        while self.waiting and budget > 0 and len(self.running) < cfg.max_num_seqs:
            r = self.waiting[0]
            reserve = cfg.admit_watermark * cfg.kv_blocks
            if self.free_blocks - r.blocks_needed(len(r.prompt)) < reserve: break
            n = min(budget, len(r.prompt))
            if not self._alloc(r, n): break
            self.waiting.pop(0); self.running.append(r); chunks.append((r, n)); budget -= n

        # 5. execute prefill chunks
        t = time.perf_counter()
        for r, n in chunks:
            ids = r.prompt[r.computed: r.computed + n]
            self.m.forward_prefill(ids, r.K, r.V, r.computed)
            r.computed += n
            tag["n_prefill_tok"] += n
        tag["prefill_ms"] = (time.perf_counter() - t) * 1e3
        tag["seg"].append(("prefill", t, time.perf_counter()))

        # 6. execute decodes: captured fast path or eager
        if decodes:
            B = len(decodes)
            i = bisect.bisect_left(cfg.capture_sizes, B)
            graph = cfg.capture_sizes[i] if i < len(cfg.capture_sizes) else None
            ids = [(r.out[-1] if r.out else r.prompt[-1]) for r in decodes]
            kvs = [(r.K, r.V, r.computed + 1) for r in decodes]
            pos = [r.computed for r in decodes]
            t = time.perf_counter()
            if graph is not None:
                pad = graph - B             # padded rows: real work, discarded
                toks = self.m.forward(ids + ids[:1] * pad, kvs + kvs[:1] * pad, pos + pos[:1] * pad)[:B]
            else:
                toks = self.m.forward_eager(ids, kvs, pos)
            tag["decode_ms"] = (time.perf_counter() - t) * 1e3
            tag["seg"].append(("decode", t, time.perf_counter()))
            tag["graph"] = graph if graph is not None else "eager"
            tag["n_decode"] = B
            now = time.perf_counter()
            sidx = len(self.steps)
            for r, tk in zip(decodes, toks):
                r.computed += 1; r.out.append(int(tk))
                r.emit_t.append(now); r.emit_step.append(sidx)
                if r.first_token_t is None: r.first_token_t = now
                tk = int(tk)
                if cfg.compact_history:
                    self.flat.extend((r.rid, tk, now, -0.1, -1.1, -2.1, -3.1, -4.1))
                else:
                    ev = TokenEvent(r, tk, now, {tk + j: Logprob(-0.1 - j, j + 1) for j in range(5)})
                    self.history.append(ev)
                if len(r.out) >= r.max_new:
                    self.running.remove(r); self._free(r); r.K = r.V = None
                    self.done.append(r)
        tag["t1"] = time.perf_counter()
        if tag["seg"] or decodes or chunks:   # don't log idle spins
            tag["seg"] = tuple(tag["seg"])
            self.steps.append(tag)

    def _preempt(self, r):
        self.running.remove(r)
        if r.computed < len(r.prompt):              # mid-prefill: just restart it
            self._free(r); r.computed = 0; r.K = r.V = None; self.waiting.insert(0, r); return
        r.evictions.append([time.perf_counter(), None])  # evicted (before the copy-out)
        r.swapped = (r.K[:, :r.computed].copy(), r.V[:, :r.computed].copy())
        r.K = r.V = None
        self._free(r); self.swapped.append(r)

    def run(self):
        start = time.perf_counter()
        for r in self.trace: r.arrival += start
        if self.cfg.gc_freeze:
            gc.collect(); gc.freeze()
        while self.trace or self.waiting or self.running or self.swapped or self.pending_tok:
            now = time.perf_counter()
            if not (self.waiting or self.running or self.swapped):
                if self.pending_tok:              # only waiting on the tokenizer worker
                    time.sleep(0.0005)
                elif self.trace[0].arrival > now:
                    time.sleep(self.trace[0].arrival - now); continue
            self.step(now)
        gc.callbacks.remove(self._gc_cb)
        if self.pool: self.pool.shutdown()
        if self.cfg.gc_freeze: gc.unfreeze()
        return start

# ---------------------------------------------------------------- workload
def make_trace(seed=1, duration=40.0, rate=6.0):
    rng = random.Random(seed)
    t, rid, out = 0.0, 0, []
    while t < duration:
        # non-homogeneous Poisson (thinning): traffic swells and ebbs on a 15 s cycle
        t += rng.expovariate(rate * 1.6)
        if rng.random() > (1 + 0.6 * np.sin(2 * np.pi * t / 15)) / 1.6: continue
        if rng.random() < 0.01:                       # someone pastes a whole document
            n_chars = rng.randint(150_000, 400_000)
        else:
            n_chars = int(rng.lognormvariate(6.9, 0.6))  # ~1k chars median
        out.append(Request(rid, t, make_text(rng, n_chars), max_new=min(1000, int(rng.expovariate(1 / 280)) + 16)))
        rid += 1
    return out

def startup_heap():
    """What a real server drags around: configs, vocab maps, registries."""
    return [{"id": i, "name": f"obj{i}", "meta": [i, str(i)]} for i in range(400_000)]

attribution.py — the exclusive gap partitioner

"""Exclusive, additive partition of each inter-token gap.

Every millisecond of a gap [t0, t1] goes to exactly one label:

  preempt   the request was evicted (explicit eviction -> readmission events),
            or the step was copying KV out/in for someone's preemption
  gc        a garbage-collection pause (overrides whatever segment it landed in)
  prefill   prefill chunks executing in the step
  tok       tokenizer running on the engine loop
  graph     decode time above what a captured (batched) step of that size costs
  kernel    the rest of the decode forward pass
  residual  time covered by none of the above (loop bookkeeping, emission, OS)

While a request is evicted, the engine's activity in that window is recorded
separately as `context` -- it explains why readmission took as long as it did,
but it is not added to the gap a second time.
"""
import bisect
import numpy as np

CAUSES = ["prefill", "preempt", "graph", "gc", "tok"]
LABELS = CAUSES + ["kernel", "residual"]
SEG_LABEL = {"tok": "tok", "swap": "preempt", "prefill": "prefill", "decode": "decode"}


class GCIndex:
    """Sorted, non-overlapping GC pauses with fast overlap queries."""
    def __init__(self, gcs):
        g = sorted((float(a), float(b)) for a, b in gcs)
        self.a = [x for x, _ in g]
        self.b = [y for _, y in g]

    def overlap(self, x, y):
        if y <= x or not self.a: return 0.0
        i = max(0, bisect.bisect_left(self.b, x))     # first pause ending after x
        tot = 0.0
        while i < len(self.a) and self.a[i] < y:
            tot += max(0.0, min(self.b[i], y) - max(self.a[i], x)); i += 1
        return tot


def decode_net_ms(step, gci):
    """Decode segment duration with any GC inside it removed (ms)."""
    tot = 0.0
    for name, x, y in step["seg"]:
        if name == "decode": tot += (y - x) - gci.overlap(x, y)
    return tot * 1e3


def fit_baseline(steps, gci):
    """GC-free decode ms vs padded batch size, fit on captured-path steps only."""
    pts = [(s["graph"], decode_net_ms(s, gci)) for s in steps if isinstance(s["graph"], (int, np.integer))]
    xs, ys = zip(*pts)
    return np.polyfit(xs, ys, 1)


def _subtract(t0, t1, holes):
    """[t0, t1] minus a list of intervals -> list of live windows."""
    out, cur = [], t0
    for x, y in sorted(holes):
        x, y = max(x, t0), min(y, t1)
        if y <= x: continue
        if x > cur: out.append((cur, x))
        cur = max(cur, y)
    if cur < t1: out.append((cur, t1))
    return out


class Partitioner:
    def __init__(self, steps, gcs, baseline=None):
        self.steps = steps
        self.t0s = [s["t0"] for s in steps]
        self.gci = GCIndex(gcs)
        self.slope, self.icpt = baseline if baseline is not None else fit_baseline(steps, self.gci)
        self._dnet = {}

    def _steps_overlapping(self, x, y):
        i = max(0, bisect.bisect_right(self.t0s, x) - 1)
        while i < len(self.steps) and self.steps[i]["t0"] < y:
            if self.steps[i]["t1"] > x: yield i, self.steps[i]
            i += 1

    def _graph_fraction(self, i, st):
        """Share of this step's (GC-free) decode time that is graph-miss excess."""
        if st["graph"] != "eager": return 0.0
        if i not in self._dnet: self._dnet[i] = decode_net_ms(st, self.gci)
        full = self._dnet[i]
        if full <= 0: return 0.0
        base = self.slope * st["n_decode"] + self.icpt
        return max(0.0, full - base) / full

    def _charge(self, x, y, acc):
        """Charge window [x, y] (s) segment-by-segment into acc (ms). Returns covered ms."""
        gc = self.gci.overlap(x, y)
        acc["gc"] += gc * 1e3
        covered = gc
        for i, st in self._steps_overlapping(x, y):
            for name, a, b in st["seg"]:
                a, b = max(a, x), min(b, y)
                if b <= a: continue
                net = (b - a) - self.gci.overlap(a, b)
                covered += net
                lab = SEG_LABEL[name]
                if lab == "decode":
                    f = self._graph_fraction(i, st)
                    acc["graph"] += net * f * 1e3
                    acc["kernel"] += net * (1 - f) * 1e3
                else:
                    acc[lab] += net * 1e3
        return covered * 1e3

    def gap(self, t0, t1, evictions=()):
        """Partition one gap. Returns (parts, context): parts sum exactly to the gap."""
        parts = {k: 0.0 for k in LABELS}
        context = {k: 0.0 for k in LABELS}
        ev = [(a, b if b is not None else t1) for a, b in evictions if a < t1 and (b is None or b > t0)]
        for a, b in ev:                                   # evicted: all of it is preemption
            a, b = max(a, t0), min(b, t1)
            if b > a:
                parts["preempt"] += (b - a) * 1e3
                cov = self._charge(a, b, context)         # what the engine was doing meanwhile
                context["residual"] += (b - a) * 1e3 - cov
        for a, b in _subtract(t0, t1, ev):                # live: charge what actually ran
            cov = self._charge(a, b, parts)
            parts["residual"] += (b - a) * 1e3 - cov
        return parts, context


def label(gap, parts, floor_ms=1.0, frac=0.25):
    k = max(CAUSES, key=parts.get)
    return k if parts[k] >= floor_ms and parts[k] >= frac * gap else "clean"

test_attribution.py — unit tests, including both reported bugs

import math
from attribution import Partitioner, LABELS

def st(t0, t1, segs, graph=16, n=1):
    return {"t0": t0, "t1": t1, "seg": tuple(segs), "graph": graph, "n_decode": n}

def close(a, b, tol=1e-6): return math.isclose(a, b, abs_tol=tol)

def total(p): return sum(p[k] for k in LABELS)

def test_reviewer_example_1_eviction_not_double_counted():
    # 40 ms gap: evicted 0-30 ms while OTHER requests' prefill ran, then a 10 ms decode.
    steps = [st(0.000, 0.030, [("prefill", 0.000, 0.030)], graph=None, n=0),
             st(0.030, 0.040, [("decode", 0.030, 0.040)])]
    p, ctx = Partitioner(steps, [], baseline=(0.0, 10.0)).gap(0.0, 0.040, evictions=[(0.0, 0.030)])
    assert close(p["preempt"], 30) and close(p["kernel"], 10) and close(p["prefill"], 0)
    assert close(total(p), 40)
    assert close(ctx["prefill"], 30)          # engine activity during eviction: recorded, not added

def test_reviewer_example_2_gc_does_not_inflate_graph_excess():
    # eager decode segment of 30 ms = 10 ms work + 20 ms GC; captured baseline for B=1 is 10 ms.
    steps = [st(0.000, 0.030, [("decode", 0.000, 0.030)], graph="eager", n=1)]
    p, _ = Partitioner(steps, [(0.005, 0.025)], baseline=(0.0, 10.0)).gap(0.0, 0.030)
    assert close(p["gc"], 20) and close(p["graph"], 0) and close(p["kernel"], 10)
    assert close(total(p), 30)

def test_real_graph_excess_survives_gc_removal():
    # eager decode 40 ms = 30 ms work + 10 ms GC; baseline 10 ms -> excess 20 ms
    steps = [st(0.000, 0.040, [("decode", 0.000, 0.040)], graph="eager", n=1)]
    p, _ = Partitioner(steps, [(0.010, 0.020)], baseline=(0.0, 10.0)).gap(0.0, 0.040)
    assert close(p["gc"], 10) and close(p["graph"], 20) and close(p["kernel"], 10)

def test_partial_overlap_with_gap_is_proportional():
    # gap covers only the second half of an eager 40 ms decode (baseline 10 -> 75% excess)
    steps = [st(0.000, 0.040, [("decode", 0.000, 0.040)], graph="eager", n=1)]
    p, _ = Partitioner(steps, [], baseline=(0.0, 10.0)).gap(0.020, 0.040)
    assert close(p["graph"], 15) and close(p["kernel"], 5) and close(total(p), 20)

def test_gc_after_emission_and_residual_reconcile():
    # emission at 10 ms, GC 12-40 ms outside any segment, next decode 40-50 ms
    steps = [st(0.000, 0.045, [("decode", 0.000, 0.010)]),
             st(0.045, 0.050, [("decode", 0.045, 0.050)])]
    p, _ = Partitioner(steps, [(0.012, 0.040)], baseline=(0.0, 5.0)).gap(0.010, 0.050)
    assert close(p["gc"], 28) and close(p["kernel"], 5) and close(p["residual"], 7)
    assert close(total(p), 40)

def test_mixed_step_prefill_tok_swap_and_decode():
    steps = [st(0.0, 0.050, [("tok", 0.000, 0.005), ("swap", 0.005, 0.007),
                             ("prefill", 0.007, 0.040), ("decode", 0.040, 0.048)])]
    p, _ = Partitioner(steps, [], baseline=(0.0, 8.0)).gap(0.0, 0.050)
    assert close(p["tok"], 5) and close(p["preempt"], 2) and close(p["prefill"], 33)
    assert close(p["kernel"], 8) and close(p["residual"], 2) and close(total(p), 50)

if __name__ == "__main__":
    import sys
    fails = 0
    for name, fn in list(globals().items()):
        if name.startswith("test_"):
            try: fn(); print("PASS", name)
            except AssertionError: fails += 1; print("FAIL", name)
    sys.exit(fails)

run.py — runs one rung and dumps raw timing logs

"""Run the nano-engine under one rung of the fix ladder and dump raw timing logs.

    python3 run.py "<rung>" <seed> <duration_s>   ->  raw/<rung-slug>__s<seed>.npz

Everything reported in the article is recomputed from these files by analyze.py.
All timestamps are seconds relative to engine start (time.perf_counter()).
"""
import sys, os, re, json, time
import numpy as np
import engine as E

LADDER = [
    ("baseline",           dict()),
    ("+gc.freeze",         dict(gc_freeze=True)),
    ("+compact history",   dict(compact_history=True)),
    ("+graphs to 48",      dict(capture_sizes=(1, 2, 4, 8, 16, 24, 32, 40, 48))),
    ("+budget 256",        dict(token_budget=256)),
    ("+tokenize off-loop", dict(async_tokenize=True)),
    ("+admit watermark",   dict(admit_watermark=0.10)),
]
CONFIGS, _acc = {}, {}
for _name, _delta in LADDER:          # each rung keeps every fix above it
    _acc = {**_acc, **_delta}; CONFIGS[_name] = dict(_acc)

RATE = float(os.environ.get("ITL_RATE", "10.0"))
SEG_CODE = {"tok": 0, "swap": 1, "prefill": 2, "decode": 3}

def slug(name): return re.sub(r"[^a-z0-9]+", "-", name.lower()).strip("-")

def run(name, seed, duration):
    cfg = E.Config(name=name, **CONFIGS[name])
    E.get_tokenizer()
    heap = E.startup_heap()                       # stays alive for the run, like a real server's
    trace = E.make_trace(seed=seed, duration=duration, rate=RATE)
    tr = np.array([(r.rid, r.arrival, len(r.text), r.max_new) for r in trace])
    eng = E.Engine(cfg, E.Model(), trace)
    w0 = time.perf_counter()
    start = eng.run()
    wall = time.perf_counter() - w0
    T = lambda t: t - start
    st = eng.steps
    steps = np.array([(T(s["t0"]), T(s["t1"]), s["n_decode"], s["n_prefill_tok"],
                       -1 if s["graph"] is None else (0 if s["graph"] == "eager" else s["graph"]),
                       s["preempted"], s["swapped_in"]) for s in st])
    segs = np.array([(i, SEG_CODE[n], T(a), T(b)) for i, s in enumerate(st) for n, a, b in s["seg"]])
    gcs = np.array([(T(a), T(b), g) for a, b, g in eng.gc_log]).reshape(-1, 3)
    emis = np.array([(r.rid, T(t), k) for r in eng.done for t, k in zip(r.emit_t, r.emit_step)])
    evs = np.array([(r.rid, T(a), T(b) if b is not None else np.nan) for r in eng.done for a, b in r.evictions]).reshape(-1, 3)
    reqs = np.array([(r.rid, T(r.arrival), len(r.prompt), r.max_new, T(r.first_token_t)) for r in eng.done])
    os.makedirs("raw", exist_ok=True)
    path = f"raw/{slug(name)}__s{seed}.npz"
    np.savez_compressed(path, steps=steps, segs=segs, gcs=gcs, emissions=emis, evictions=evs,
                        requests=reqs, trace=tr, wall=wall,
                        config=json.dumps({"rung": name, "seed": seed, "duration": duration,
                                           "rate": RATE, **{k: v for k, v in vars(cfg).items()}}, default=str))
    del heap
    return path

if __name__ == "__main__":
    print(run(sys.argv[1], int(sys.argv[2]), float(sys.argv[3])))

analyze.py — recomputes every number from the raw logs

"""Recompute every reported number from raw/*.npz.

    python3 analyze.py raw/*.npz > results.jsonl
"""
import sys, json
import numpy as np
from attribution import Partitioner, LABELS, CAUSES, label

SEG_NAME = {0: "tok", 1: "swap", 2: "prefill", 3: "decode"}
BINS = [0, 5, 10, 15, 20, 30, 50, 80, 120, 200, 400, 800, 1e9]

def load(path):
    z = np.load(path, allow_pickle=False)
    cfg = json.loads(str(z["config"]))
    steps = []
    for t0, t1, nd, npf, g, pre, swi in z["steps"]:
        g = int(g)
        steps.append({"t0": t0, "t1": t1, "n_decode": int(nd), "n_prefill_tok": int(npf),
                      "graph": None if g == -1 else ("eager" if g == 0 else g),
                      "preempted": int(pre), "seg": []})
    for i, c, a, b in z["segs"]:
        steps[int(i)]["seg"].append((SEG_NAME[int(c)], a, b))
    for s in steps: s["seg"] = tuple(s["seg"])
    return z, cfg, steps

def analyze(path):
    z, cfg, steps = load(path)
    first_t0 = steps[0]["t0"]
    gcs_all = z["gcs"]
    startup_gc = [g for g in gcs_all if g[1] <= first_t0]          # e.g. gc.collect() before freeze
    serving_gc = [g for g in gcs_all if g[1] > first_t0]
    P = Partitioner(steps, [(a, b) for a, b, _ in gcs_all])
    emis, evs = z["emissions"], z["evictions"]
    ev_by = {}
    for rid, a, b in evs: ev_by.setdefault(int(rid), []).append((a, None if np.isnan(b) else b))
    by_req = {}
    for rid, t, k in emis: by_req.setdefault(int(rid), []).append((t, int(k)))

    gaps, parts, labs, ctx_tot = [], [], [], {k: 0.0 for k in LABELS}
    worst_err, skipped_unexplained = 0.0, 0
    for rid, em in by_req.items():
        ev = ev_by.get(rid, [])
        for (ta, ka), (tb, kb) in zip(em, em[1:]):
            p, ctx = P.gap(ta, tb, ev)
            g = (tb - ta) * 1e3
            worst_err = max(worst_err, abs(sum(p.values()) - g))
            if kb > ka + 1 and not any(a < tb and (b is None or b > ta) for a, b in ev):
                skipped_unexplained += 1
            for k in LABELS: ctx_tot[k] += ctx[k]
            gaps.append(g); parts.append(p); labs.append(label(g, p))
    gaps = np.array(gaps)
    p50, p90, p99, p999 = np.percentile(gaps, [50, 90, 99, 99.9])
    tail = gaps >= p99
    tail_ms = gaps[tail].sum()
    share = {k: float(sum(p[k] for p, t in zip(parts, tail) if t) / tail_ms) for k in LABELS}
    share_all = {k: float(sum(p[k] for p in parts) / gaps.sum()) for k in LABELS}
    kshare = np.array([p["kernel"] / g for p, g in zip(parts, gaps)])
    clean = np.array([l == "clean" for l in labs])
    req = z["requests"]
    ttft = (req[:, 4] - req[:, 1]) * 1e3
    n_tok = len(emis)
    hist = []
    for lo, hi in zip(BINS[:-1], BINS[1:]):
        sel = [l for g, l in zip(gaps, labs) if lo <= g < hi]
        hist.append({"lo": lo, "hi": hi, "n": len(sel), **{k: sel.count(k) for k in CAUSES + ["clean"]}})
    st = z["steps"]
    return {
        "rung": cfg["rung"], "seed": cfg["seed"], "duration": cfg["duration"],
        "requests": int(len(req)), "tokens": int(n_tok), "gaps": int(len(gaps)),
        "tok_per_s": n_tok / float(z["wall"]),
        "itl_p50": p50, "itl_p90": p90, "itl_p99": p99, "itl_p999": p999, "itl_max": float(gaps.max()),
        "ttft_p50": float(np.percentile(ttft, 50)), "ttft_p99": float(np.percentile(ttft, 99)),
        "tail_share": share, "all_share": share_all,
        "tail_count": {k: int(sum(1 for l, t in zip(labs, tail) if t and l == k)) for k in CAUSES + ["clean"]},
        "fwd_share_clean_median": float(np.median(kshare[clean])),
        "fwd_share_tail_median": float(np.median(kshare[tail])),
        "evicted_context_ms": ctx_tot,
        "reconcile_max_err_ms": worst_err, "skipped_steps_without_eviction": skipped_unexplained,
        "preemptions": int(st[:, 5].sum()), "eager_steps": int((st[:, 4] == 0).sum()),
        "max_decode_B": int(st[:, 2].max()),
        "gen2_serving_ms": [round((b - a) * 1e3, 1) for a, b, g in serving_gc if g == 2],
        "gen2_startup_ms": [round((b - a) * 1e3, 1) for a, b, g in startup_gc if g == 2],
        "gc_counts": {int(g): int((gcs_all[:, 2] == g).sum()) for g in (0, 1, 2)},
        "hist": hist,
    }

if __name__ == "__main__":
    for path in sys.argv[1:]:
        print(json.dumps(analyze(path), default=float), flush=True)

Six Posts on Ensembles and Tuning, and the Uncomfortable Thing They All Turned Out to Be About

Sachin Kr. Rajput — Tue, 25 Aug 2026 10:45:31 +0000

The One-Line Summary: I set out to write six posts comparing ensemble and tuning methods, and every single one ended up measuring the same thing instead — that the gap between the methods was smaller than the spread of one method against its own seed, which means most published comparisons of these techniques are reporting luck with a method's name attached to it.

The Arc

Six posts, two threads: what to do when you have several models, and how to search for the settings of one. I did not plan them as a set about measurement error. They became one anyway, and by the fourth I had stopped being surprised.

Here is each, with the number that survived.

Stacking: The Editor Who Never Reported a Story

Train a small model to combine your other models, on predictions they made for rows they never saw. Done properly it reached 0.9331 against 0.9301 for the best single model — a real gain, and a small one.

Done improperly it was worse than useless. Train the combiner on in-fold predictions and it handed a weight of 1.123 to the Random Forest and 0.035 to the Gradient Booster that was genuinely better. The leak did not merely inflate the score; it selected the wrong model, because a bagged forest predicting rows it memorised looks like an oracle.

Blending and Voting: Four Noses, One Bottle

Voting has no learned parameters, so it cannot overstate itself — and it lost to the best single model on both metrics, 0.9296 for hard and 0.9289 for soft. Soft voting did not beat hard voting, because it averaged in members whose probabilities were badly calibrated.

Blending, with fitted weights, beat everything at 0.2055 when it had 150 holdout rows per weight. On 40 rows it reported 0.1492 and delivered 0.2576 — worse than the free unweighted average at 0.2515. Same code, same models. Only the holdout size changed.

Grid Search vs Random Search: The Safe With Three Dials

This is where the pattern became impossible to ignore. Random search's advantage is a real mechanism and it behaved exactly as advertised. The gaps were 0.0008 against it on three parameters and 0.0039 for it on six.

Random search's spread against its own seed was 0.0098. Grid beat exactly five of ten random seeds. A single-seed experiment had a coin-flip chance of producing either headline.

Bayesian Optimization: The Prospector

The one method whose advantage grew with budget — 0.0032, then 0.0051, then 0.0063 at 20, 50 and 120 trials. That growth is the signature of something actually learning from its own results.

But the finding that changed how I think about it was the spread. My thirty-line version landed within 0.0004 across three seeds. Random search spanned 0.0128 — thirty-two times wider, for a mean advantage of only 0.0015. The product is not a better answer. It is an answer that does not depend on luck.

Replace Your Grid Search Tonight

Two lines of pruning killed 12 of 20 trials and cut wall clock by 28%, costing 0.1908 to 0.1913 of log loss. The kill rate and the speedup do not match, because a pruned trial still pays for its warm-up. Anyone quoting the kill rate as the time saved is quoting the brochure.

The Hyperparameter Search Cheat Sheet

All of it on one page. Writing it was more uncomfortable than writing the six, because a cheat sheet wants to say "use X" and every measurement says "it depends, and by less than you think."

The Thing They Were All About

Line the effect sizes up next to the noise and the arc writes itself.

Comparison	Claimed effect	Noise in the same setup
Grid vs random, 3 params	0.0008	0.0098
Grid vs random, 6 params	0.0039	0.0098
TPE vs random, 20 trials	0.0032	0.0128
Blending vs averaging, thin holdout	claims 0.1492, delivers 0.2576	—
Stacking, leaked vs clean	picks the wrong model entirely	—

The top three are inside their own noise. The bottom two are not subtle at all — but both are failures of procedure, not of method choice. Nobody in those two rows picked the wrong algorithm. They measured it wrong.

That is the uncomfortable through-line. Six posts spent arguing about which technique wins, and the honest answer in every case was: run it twice before you believe it, and check that the thing you are comparing is bigger than the thing you are not controlling.

What Confused People

Nothing, and I am going to keep being straight about that rather than inventing a reader question.

These six drew zero substantive comments between them. That is not a mystery and it is not the writing — it is distribution, and I made it worse on purpose this week by publishing four of them inside ten minutes. The view counts show exactly what that costs: the first two published picked up 20 and 12 views, the last two got 1 and 0. Same author, same tags, same day, four minutes apart.

I have written about burst publishing being the reason two-thirds of this archive never got read, and then did it again because clearing a backlog felt more urgent than spacing it. The measurement is right there in the feed. One per day from here.

Where the real conversation has been happening is in other people's comment threads, which is its own lesson about where reach actually comes from on this platform.

The One-Sentence Summary

Across six posts on ensembling and hyperparameter search, every method comparison I ran came out smaller than the run-to-run spread of the methods being compared — 0.0008 and 0.0039 against a seed spread of 0.0098, 0.0032 against 0.0128 — and the only two effects that were unambiguously large, the stacking leak that put 1.123 on the wrong model and the blend that claimed 0.1492 while delivering 0.2576, were both procedural mistakes rather than wrong choices of algorithm, which means the highest-return skill in this entire area is not picking the right method but running the comparison twice.

What's Next?

The tuning thread is not finished — two more posts before it closes:

Early stopping — knowing when to quit, on the training curve rather than the search.
Model calibration — the weather app that says 70% and means it half the time.
K-nearest neighbours — opening a new series on distance-based learning.
Nested cross-validation — tuning and reporting without burning the same data twice.

Follow me for the next article in the Hyperparameter Tuning series!

Let's Connect!

If any single number here changed something you do, drop a heart!

Questions? Ask in the comments — I read and respond to every one, and after the section above you can see I would genuinely welcome the chance.

What's the last method comparison you ran twice? I did not run any of these twice until the grid-versus-random post forced me to, and the only reason it forced me to was that seed 0 gave me the opposite of the famous result and I got suspicious of my own good luck. ⚖️

The pattern in these six was not planned and I would not have believed it if someone had asserted it to me. What made it visible was a boring habit — running the same experiment with a different seed and reporting the spread instead of the best number. That habit costs one extra loop and it invalidated roughly half of what I expected to write. I think that is the correct ratio, and I think most published method comparisons have not paid it.

Before your next "X beats Y" post, run X twice and see whether it beats itself.

Bayesian Optimization: The Prospector Who Reads Every Core Sample Before Choosing Where to Drill

Sachin Kr. Rajput — Mon, 24 Aug 2026 10:15:27 +0000

The One-Line Summary: Bayesian optimization fits a model of your loss surface and drills where that model is both promising and uncertain — and the advantage it bought was not mainly a better answer but a far more reliable one: my thirty-line version landed within 0.0004 across three seeds while random search spanned 0.0128, which is a lottery ticket versus a procedure.

The Parable of the Kestrel Basin Survey

The company has rights to eleven thousand hectares and enough money to drill forty holes. Each hole costs a week and returns one number: how much ore is in the core. There is no map, no seismic, nothing — just eleven thousand hectares and forty chances to find out where the deposit is.

The question the geologist has to answer before the first hole is not where the ore is. It is how she will decide where hole seventeen goes.

The Lottery

The first surveyor's method is to scatter the forty holes at random across the whole basin. It is unbiased, it covers the ground, and it requires no thinking between holes.

THE LOTTERY
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
Forty holes, chosen before any drilling starts:

  hole  1   NW quadrant   0.4 g/t
  hole  2   SE edge       0.1 g/t
  hole  3   NW quadrant   3.1 g/t   <- interesting!
  hole  4   S basin       0.2 g/t
  ...
  hole 40   E ridge       0.3 g/t

  best found: 3.1 g/t at hole 3

Hole 4 was chosen before hole 3 came back. So were
holes 5 through 40. The 3.1 reading changed nothing
about where anyone drilled next.

✗ Thirty-seven holes were sited by a plan that
  could not have known about hole 3.

This is called RANDOM SEARCH.

Every hole is independent of every other. The information arrives and is filed, and the drill goes where the envelope said it would go.

The Prospector

The second surveyor drills eight holes at random to get a feel for the ground, and then stops. She builds a rough model of the whole basin from those eight readings — a guess at the grade everywhere, plus a guess at how unsure she is everywhere. Then she drills hole nine wherever that model says the upside is largest.

"I am not looking for the best hole I have drilled. I am looking for the place my ignorance is most expensive."

THE PROSPECTOR
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
After 8 random holes, she fits a surface:

  location      predicted   uncertainty   score
  near hole 3      2.8          0.4        3.2
  far N ridge      0.9          2.1        3.0   <- !
  near hole 1      0.5          0.3        0.8

  score = predicted + k x uncertainty

She drills the far N ridge, where she predicts
almost nothing, because she barely knows anything
there and a surprise would be worth a lot.

Then she refits the surface with 9 points and asks
again. And again, thirty-one more times.

This is called BAYESIAN OPTIMIZATION.

Two things are being traded against each other, and the constant k decides the exchange rate: exploitation (drill near the good reading) and exploration (drill where you know least). Set k to zero and she drills the same hill forty times. Set it huge and she reinvents the lottery.

What The Survey Actually Showed

Both surveyors were run three times, on three different random starts, and the numbers that came back were not the ones the company expected.

THREE SURVEYS EACH
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
The lottery, three attempts:
  0.8687   0.8747   0.8816     <- spread 0.0128

The prospector, three attempts:
  0.8763   0.8765   0.8768     <- spread 0.0004

On the third attempt the lottery WON, and by more
than the prospector's average margin.

But look at the spreads. The lottery's outcome is
a draw. The prospector's is a result.

✓ The prospector is not reliably better.
  The prospector is reliably.

The company had been asking which surveyor finds more ore. The useful question was which surveyor they could plan around.

Why It Works

THE MATHEMATICS OF ASKING BETTER QUESTIONS
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
RANDOM SEARCH
  x_(n+1) ~ P(x)          drawn from a fixed prior
  → trial n+1 is independent of trials 1..n
  → information accumulates and is never used

BAYESIAN OPTIMIZATION
  1. fit a SURROGATE model  f-hat(x)  to all
     trials so far, returning mean AND spread
  2. maximise an ACQUISITION function over it,
     e.g. UCB:   mu(x) + k * sigma(x)
  3. evaluate the true objective there, once
  4. refit and repeat

  REQUIREMENTS
  1. The objective must be EXPENSIVE relative to
     fitting the surrogate. Otherwise just sample.
  2. The surrogate must report uncertainty, not
     only a prediction.
  3. Trials must be SEQUENTIAL. Parallelism costs
     you the thing that makes it work.
  4. Roughly 8-15 random trials first, or the
     surrogate is fitting noise.

What Is Bayesian Optimization?

It is a loop with two models in it. The objective is the thing you actually care about and cannot afford to evaluate often — a cross-validated model score, a simulation, a physical experiment. The surrogate is a cheap statistical model of the objective, fitted to every evaluation so far, which must return both a prediction and an uncertainty.

The acquisition function turns those two numbers into a single score that says how much you want to try a candidate point. Upper Confidence Bound is the simplest honest one:

a_{UCB} (x) = μ (x) + k \cdot σ (x)

Expected Improvement is the more common choice and it does the same job with a different shape — it integrates how much better than the current best a point might be, rather than adding a fixed multiple of the standard deviation.

The classical surrogate is a Gaussian process, which is where "Bayesian" comes from: the GP is a prior over functions, updated by observations into a posterior [1]. In practice the surrogate is often a tree ensemble or, in Optuna's default, a Tree-structured Parzen Estimator — which does not model $p (y ∣ x)$ at all but instead models $p (x ∣ y)$ , building separate densities over the good trials and the bad ones and sampling where their ratio is favourable [2].

Step by Step

Evaluate $n_{init}$ random points. Eight to fifteen.
Fit the surrogate on everything observed so far.
Maximise the acquisition function over a large set of cheap candidates.
Evaluate the true objective at the argmax. Once.
Add it to the observations and go to 2.

Step 3 is the part that surprises people: you are running an optimisation inside an optimisation. That is fine, because the inner one is over a model that costs microseconds, and the outer one is over a model fit that costs seconds or hours.

TPE Against Random, Three Budgets

Optuna's TPE sampler against Optuna's random sampler, identical objective, five seeds each, so the spread is visible.

import warnings, numpy as np; warnings.filterwarnings("ignore")
import optuna; optuna.logging.set_verbosity(optuna.logging.WARNING)
from sklearn.datasets import make_friedman1
from sklearn.model_selection import train_test_split, KFold, cross_val_score
from sklearn.ensemble import GradientBoostingRegressor

X, y = make_friedman1(n_samples=800, noise=1.0, random_state=42)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=42)
cv = KFold(3, shuffle=True, random_state=0)

def score(lr, depth, sub, mf):
    m = GradientBoostingRegressor(n_estimators=60, learning_rate=lr, max_depth=depth,
                                  subsample=sub, max_features=mf, random_state=0)
    return cross_val_score(m, Xtr, ytr, cv=cv, scoring='r2', n_jobs=-1).mean()

def objective(t):
    return score(t.suggest_float("lr", 0.02, 0.25, log=True),
                 t.suggest_int("depth", 2, 5),
                 t.suggest_float("sub", 0.6, 1.0),
                 t.suggest_float("mf", 0.5, 1.0))

def run(sampler_name, budget, seed):
    s = optuna.samplers.TPESampler(seed=seed) if sampler_name == "tpe" \
        else optuna.samplers.RandomSampler(seed=seed)
    st = optuna.create_study(direction="maximize", sampler=s)
    st.optimize(objective, n_trials=budget, show_progress_bar=False)
    return st.best_value

print("TPE vs RANDOM, equal budget, 5 seeds each")
print("-" * 62)
print(f"{'budget':>7}{'TPE mean':>11}{'TPE spread':>12}{'rand mean':>11}"
      f"{'rand spread':>13}{'diff':>9}")
for budget in (20, 50, 120):
    tpe = [run("tpe", budget, s) for s in range(5)]
    rnd = [run("rand", budget, s) for s in range(5)]
    print(f"{budget:>7}{np.mean(tpe):>11.4f}{max(tpe)-min(tpe):>12.4f}"
          f"{np.mean(rnd):>11.4f}{max(rnd)-min(rnd):>13.4f}"
          f"{np.mean(tpe)-np.mean(rnd):>+9.4f}")

TPE vs RANDOM, equal budget, 5 seeds each
--------------------------------------------------------------
 budget   TPE mean  TPE spread  rand mean  rand spread     diff
     20     0.8737      0.0146     0.8704       0.0077  +0.0032
     50     0.8781      0.0074     0.8730       0.0125  +0.0051
    120     0.8803      0.0090     0.8741       0.0107  +0.0063

TPE wins at all three budgets, and — this is the part that matters — the margin grows with budget: +0.0032, +0.0051, +0.0063.

That growth is the signature of a method that is actually learning. Yesterday's grid-versus-random comparison produced gaps that flipped sign with dimensionality and vanished into seed noise, because neither method used its own results. Here the gap widens monotonically, because every extra trial makes the surrogate better, which makes the next trial better chosen. Random search at 120 trials (0.8741) has barely improved on random search at 50 (0.8730), while TPE gained 0.0022 over the same stretch.

Being honest about the caveat: at budget 20 the gap is +0.0032 while TPE's own spread is 0.0146. One run at a small budget tells you nothing. The direction is consistent across all three budgets, which is the evidence — not the size of any single comparison.

From Scratch: A Surrogate and a Rule

Thirty lines. Fit a random forest to the trials so far, use the disagreement between its trees as the uncertainty, and go wherever mean plus 1.5 standard deviations is highest.

import warnings, numpy as np; warnings.filterwarnings("ignore")
from sklearn.datasets import make_friedman1
from sklearn.model_selection import train_test_split, KFold, cross_val_score
from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor

X, y = make_friedman1(n_samples=800, noise=1.0, random_state=42)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=42)
cv = KFold(3, shuffle=True, random_state=0)
BOUNDS = np.array([[np.log(0.02), np.log(0.25)], [2, 5], [0.6, 1.0], [0.5, 1.0]])

def score(v):
    return cross_val_score(GradientBoostingRegressor(
        n_estimators=60, learning_rate=float(np.exp(v[0])), max_depth=int(round(v[1])),
        subsample=float(v[2]), max_features=float(v[3]), random_state=0),
        Xtr, ytr, cv=cv, scoring='r2', n_jobs=-1).mean()

def sample(rng, n):
    return rng.uniform(BOUNDS[:, 0], BOUNDS[:, 1], size=(n, 4))

def bayes_opt(budget, seed, n_init=8, kappa=1.5, n_cand=400):
    """Fit a forest on the trials so far; go where mean + kappa*sd is highest."""
    rng = np.random.RandomState(seed)
    Xs = sample(rng, n_init); ys = np.array([score(v) for v in Xs])
    for _ in range(budget - n_init):
        sur = RandomForestRegressor(n_estimators=60, min_samples_leaf=2,
                                    random_state=0, n_jobs=-1).fit(Xs, ys)
        cand = sample(rng, n_cand)
        per_tree = np.stack([t.predict(cand) for t in sur.estimators_])
        ucb = per_tree.mean(axis=0) + kappa * per_tree.std(axis=0)
        nxt = cand[int(np.argmax(ucb))]
        Xs = np.vstack([Xs, nxt]); ys = np.append(ys, score(nxt))
    return ys.max()

def random_search(budget, seed):
    rng = np.random.RandomState(seed)
    return max(score(v) for v in sample(rng, budget))

print("FROM SCRATCH: FOREST SURROGATE + UCB vs RANDOM   (3 seeds, budget 50)")
print("-" * 68)
bo = [bayes_opt(50, s) for s in range(3)]
rs = [random_search(50, s) for s in range(3)]
for s, (b, r) in enumerate(zip(bo, rs)):
    print(f"  seed {s}   mine {b:.4f}   random {r:.4f}   diff {b-r:+.4f}")
print()
print(f"  mine   mean {np.mean(bo):.4f}   spread {max(bo)-min(bo):.4f}")
print(f"  random mean {np.mean(rs):.4f}   spread {max(rs)-min(rs):.4f}")
print(f"  mine - random = {np.mean(bo)-np.mean(rs):+.4f}")
print(f"  (optuna TPE at budget 50 scored 0.8781 mean over 5 seeds)")

FROM SCRATCH: FOREST SURROGATE + UCB vs RANDOM   (3 seeds, budget 50)
--------------------------------------------------------------------
  seed 0   mine 0.8763   random 0.8687   diff +0.0076
  seed 1   mine 0.8765   random 0.8747   diff +0.0019
  seed 2   mine 0.8768   random 0.8816   diff -0.0048

  mine   mean 0.8765   spread 0.0004
  random mean 0.8750   spread 0.0128

Read the seed rows first, because they are embarrassing and they are the point. On seed 2, random search beat my optimizer by 0.0048 — more than my average margin over it. If I had run one seed and it happened to be seed 2, the honest headline would have been that Bayesian optimization lost.

Now read the spreads. My three runs landed at 0.8763, 0.8765 and 0.8768 — a spread of 0.0004. Random search's three runs spanned 0.0128, thirty-two times wider.

That reframes the entire method. The mean advantage was +0.0015, which is nothing. The variance reduction was a factor of thirty-two, which is everything. Random search occasionally finds something excellent and occasionally finds something mediocre, and you cannot tell which happened without a baseline. The surrogate search converges to nearly the same answer every time, because after the initial random phase it is following the data rather than the dice.

That is a different value proposition from the one usually advertised. Bayesian optimization is not primarily a way to get a better model. It is a way to stop your tuning result from being a function of your seed.

And the honest scoreboard: Optuna's TPE reached 0.8781 at the same budget, against my 0.8765. Thirty lines gets you most of the mechanism and none of the years of tuning that went into a real sampler. Use the library.

When Does It Help Most?

	Random search	Bayesian optimization
Uses previous trials	no	yes
Parallelises	perfectly	poorly, it is sequential
Gain grows with budget	no (0.8730 → 0.8741)	yes (+0.0032 → +0.0063)
Run-to-run spread	0.0128	0.0004
Overhead per trial	none	a surrogate fit
Needs warm-up trials	no	8–15
Good below ~20 trials	yes	not really

Reach for Bayesian optimization when each evaluation is genuinely expensive — minutes or more — when you have 30 or more trials of budget, when trials run sequentially anyway, and when you need the result to be reproducible rather than lucky.

Reach for random search when evaluations are cheap, when you have hundreds of cores and want them all busy, or when your budget is under about twenty trials and the surrogate will not have enough data to be worth its overhead.

The failure mode worth naming: if your objective is fast, the surrogate fit and the candidate scan can cost more wall clock than the evaluations they are meant to save. Bayesian optimization is an answer to expensive objectives, and it is a tax on cheap ones.

Quick Reference Card

BAYESIAN OPTIMIZATION: CHEAT SHEET
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
THE LOOP
  1. n_init random trials  (8-15)
  2. fit surrogate -> mean + uncertainty
  3. maximise acquisition, e.g. mu + k*sigma
  4. evaluate objective ONCE there
  5. refit, repeat

MEASURED (GBM on Friedman #1, 4 params)
  budget 20   TPE 0.8737   random 0.8704   +0.0032
  budget 50   TPE 0.8781   random 0.8730   +0.0051
  budget 120  TPE 0.8803   random 0.8741   +0.0063
     -> the margin GROWS with budget

  from scratch, budget 50, 3 seeds
     mine   0.8765 mean, spread 0.0004
     random 0.8750 mean, spread 0.0128
     -> 32x tighter, not 32x better

KNOBS
  n_init      8-15; fewer and you fit noise
  kappa       0 = greedy, huge = random search
  surrogate   GP (smooth, small d) / TPE / forest

WHEN TO USE
  ✓ objective costs minutes+, budget 30+
  ✓ sequential trials, want reproducibility
WHEN NOT
  ✗ cheap objective (surrogate overhead dominates)
  ✗ budget < 20
  ✗ hundreds of idle cores

OPTUNA
  optuna.create_study(sampler=TPESampler(seed=0))

Key Takeaways

It beat random at every budget tested — +0.0032, +0.0051 and +0.0063 at 20, 50 and 120 trials.
The margin grows with budget, which is the signature of learning. Random search went 0.8730 to 0.8741 from 50 to 120 trials; TPE went 0.8781 to 0.8803.
The real prize is variance, not mean — my from-scratch version spread 0.0004 across three seeds where random search spread 0.0128.
A single seed can reverse the conclusion — on seed 2 random search beat my optimizer by 0.0048.
Thirty lines gets you most of it — a forest surrogate plus UCB reached 0.8765 against Optuna TPE's 0.8781.
But use the library — that remaining 0.0016 is years of sampler engineering you do not need to repeat.
It is sequential by design, so it trades away the one thing random search is perfect at. With idle cores, random search may finish sooner in wall clock.
kappa is the whole personality — zero makes it greedy and it re-drills one hill; enormous makes it random search with extra steps.

The One-Sentence Summary

Bayesian optimization replaces "where should I look next?" with "where is my ignorance most expensive?" by fitting a model of the loss surface and drilling where predicted value plus uncertainty is highest — and measured over three budgets its advantage over random search grew monotonically from +0.0032 to +0.0063 because every trial improves the model that chooses the next one, but the finding that actually changes how you should think about it is that my thirty-line version landed within 0.0004 across three seeds while random search spanned 0.0128, which means its real product is not a better answer but an answer that does not depend on luck.

What's Next?

The 5-minute Optuna setup — tomorrow. Replacing your grid search tonight, with pruning.
Nested cross-validation — tuning and reporting without burning the same data twice.
Early stopping inside a search — killing hopeless trials at iteration 10 instead of 200.
Multi-objective tuning — when you need accuracy and latency.

Follow me for the next article in the Hyperparameter Tuning series!

Let's Connect!

If the Kestrel Basin survey made the exploration-exploitation trade click, drop a heart!

Questions? Ask in the comments — I read and respond to every one.

What's your n_init? Mine is 10 and I arrived at it badly — I ran with 3 for a year, kept getting searches that fixated on whatever the first lucky trial found, and never connected the two until I plotted where the trials went. 🪨

References

[1] J. Snoek, H. Larochelle and R. P. Adams, Practical Bayesian Optimization of Machine Learning Algorithms (2012), NeurIPS 2012

[2] J. Bergstra, R. Bardenet, Y. Bengio and B. Kégl, Algorithms for Hyper-Parameter Optimization (2011), NeurIPS 2011

The thing I did not expect from these measurements is that the interesting quantity was the spread rather than the mean. We are trained to compare methods by average performance, and on average these two are close enough that you could argue either way with one seed each — which is exactly what most blog comparisons do. The variance is where the difference lives, and variance is invisible unless you deliberately run the same thing several times and resist the urge to report only the best one.

Send this to whoever is still running a 512-point grid overnight.

Grid Search vs Random Search: The Safe With Three Dials Where Only One Was Connected

Sachin Kr. Rajput — Mon, 24 Aug 2026 10:11:15 +0000

The One-Line Summary: Random search is supposed to beat grid search because grid wastes budget on parameters that don't matter, and the mechanism is exactly as advertised — but measured at a 64-trial budget the whole effect was 0.0039 in random's favour on six parameters and 0.0008 against it on three, while random search's spread against its own seed was 0.0098, so grid beat precisely 5 of 10 random seeds and the famous comparison came out a coin flip.

The Parable of the Safe at Voutier & Fils

The safe in the back office has three dials, each numbered 0 to 99, and the combination has been lost since the war. The firm hires a locksmith who explains that there are a million combinations and he charges by the hour, so the question is not whether he can open it but in what order he tries.

What nobody in the room knows yet is that only one of the three dials is connected to anything.

The Methodical Method

The locksmith proposes a grid. He will pick four settings on each dial — 0, 25, 50, 75 — and work through every combination systematically, so that nothing is ever tried twice and he can report exactly what has been eliminated.

THE METHODICAL METHOD
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
Four settings per dial, three dials:

  4 x 4 x 4  =  64 attempts

  dial A tried at:  0, 25, 50, 75
  dial B tried at:  0, 25, 50, 75
  dial C tried at:  0, 25, 50, 75

64 attempts. Four distinct positions of dial A.

He is thorough, he is auditable, and every one of
those 64 attempts moved dial B and dial C, which
are attached to nothing at all.

This is called GRID SEARCH.

Sixty-four hours of work bought four guesses at the only dial that mattered.

The Careless Method

His apprentice, who is not methodical, spins all three dials at random and tries whatever comes up. Sixty-four times.

"You keep testing the same four numbers on the dial that opens it. I test sixty-four."

THE CARELESS METHOD
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
Sixty-four random spins of all three dials:

  dial A tried at:  7, 41, 88, 3, 62, 19 ... 64 values
  dial B tried at:  ... 64 values, all wasted
  dial C tried at:  ... 64 values, all wasted

Same 64 attempts. Same wasted dials. But SIXTEEN
TIMES the resolution on the dial that matters,
because he never spends two attempts on the same
setting of it.

He does not know which dial is connected either.
He does not need to.

This is called RANDOM SEARCH.

The apprentice's advantage is not cleverness. It is that he refuses to repeat himself on any axis, so he cannot waste resolution on the axis that turns out to matter.

The Part The Story Usually Leaves Out

This is where the parable is normally allowed to end, and it should not be, because the firm ran both methods and wrote the results down.

WHAT ACTUALLY HAPPENED
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
When all three dials were connected:
  ✓ the methodical locksmith won

When five of six dials were decoration:
  ✓ the apprentice won

Then the apprentice repeated his 64 spins on ten
different days. His results varied more between
days than the gap between the two methods.

✗ The firm had spent a decade arguing about a
  difference smaller than the apprentice's luck.

The apprentice's argument is correct. At sixty-four attempts it is also almost entirely theoretical.

Why It Works

THE MATHEMATICS OF NOT REPEATING YOURSELF
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
GRID over d parameters with n total trials:

    values per parameter = n^(1/d)

  d=2, n=64  ->  8 values each
  d=3, n=64  ->  4 values each
  d=6, n=64  ->  2 values each
  d=8, n=64  ->  fewer than 2; the grid cannot
                 even be built

  → resolution COLLAPSES as you add parameters,
    including the ones that do nothing

RANDOM over d parameters with n trials:

    values per parameter = n,  for every d

  → resolution is INDEPENDENT of dimension

  P(at least one draw in the top q fraction)
      = 1 - (1-q)^n

  This is the whole theoretical case, and it says
  nothing about how much BETTER the top q region
  is — which is why the effect size has to be
  measured rather than assumed.

What Are Grid and Random Search?

Grid search enumerates the Cartesian product of a fixed set of values per parameter. It is exhaustive over that grid, reproducible without a seed, trivially parallel, and it tells you precisely which region you have ruled out.

Random search draws each trial independently from a distribution per parameter. It is not exhaustive over anything, and its guarantee is probabilistic: with $n$ draws, the chance of landing at least once inside the best $q$ fraction of the space is $1 - (1 - q)^{n}$ , independent of how many parameters there are.

That independence is the entire argument, and it comes from Bergstra and Bengio [1]. Their claim is not that random sampling is smarter. It is that most hyperparameter spaces have low effective dimensionality — a couple of parameters do nearly all the work — and grid search cannot know which, so it spends its resolution uniformly and therefore wrongly.

The Math

For grid search over $d$ parameters with $n$ total trials, the values tested per parameter is:

k = n^{1/ d}

At $n = 64$ : three parameters gives $k = 4$ , six gives $k = 2$ , and eight gives $k < 2$ — the grid cannot be constructed without dropping a parameter.

For random search, $k = n$ regardless of $d$ , and the coverage probability is:

P (hit top q) = 1 - (1 - q)^{n}

Note what this does not say. It bounds the probability of entering the top $q$ region; it says nothing about how much better that region is than the rest. If your loss surface is flat near the optimum, landing in the top 5% is worth almost nothing — and that gap between "found a good region" and "the good region is meaningfully better" is where the measured effect sizes below come from.

Step by Step

Grid: choose values per parameter, take the product, evaluate all, keep the best.

Random: choose a distribution per parameter — log-uniform for learning rates and regularisation, uniform for fractions, integer-uniform for depths — draw $n$ times, evaluate, keep the best.

The distribution choice is the part people skip, and it matters more than the grid-versus-random question. A learning rate sampled uniformly on [0.001, 0.3] puts roughly 90% of its draws above 0.03, which is not what you meant.

Equal Budget, Both Regimes

Sixty-four trials each, three-fold CV, gradient boosting on Friedman #1. The first regime gives both methods three parameters that genuinely matter. The second adds three that barely do — which is what forces grid's resolution down to two values per axis.

import warnings, numpy as np, itertools; warnings.filterwarnings("ignore")
from sklearn.datasets import make_friedman1
from sklearn.model_selection import train_test_split, KFold, cross_val_score
from sklearn.ensemble import GradientBoostingRegressor

X, y = make_friedman1(n_samples=800, noise=1.0, random_state=42)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=42)
cv = KFold(3, shuffle=True, random_state=0)

def score(p):
    return cross_val_score(GradientBoostingRegressor(n_estimators=60, random_state=0, **p),
                           Xtr, ytr, cv=cv, scoring='r2', n_jobs=-1).mean()

def run_grid(axes):
    keys = list(axes); best = -9; n = 0
    for vals in itertools.product(*axes.values()):
        n += 1; v = score(dict(zip(keys, vals)))
        if v > best: best = v
    return best, n, {k: len(v) for k, v in axes.items()}

def draw(spec, rng):
    p = {}
    for k, (kind, a, b) in spec.items():
        p[k] = float(np.exp(rng.uniform(np.log(a), np.log(b)))) if kind == 'logu' else \
               (int(rng.randint(a, b + 1)) if kind == 'int' else float(rng.uniform(a, b)))
    return p

def run_random(spec, budget, seed=0):
    rng = np.random.RandomState(seed); best = -9; seen = {k: set() for k in spec}
    for _ in range(budget):
        p = draw(spec, rng)
        for k in p: seen[k].add(round(p[k], 6))
        v = score(p)
        if v > best: best = v
    return best, {k: len(s) for k, s in seen.items()}

S3 = {"learning_rate":('logu',0.02,0.25), "max_depth":('int',2,5),
      "subsample":('unif',0.6,1.0)}
S6 = dict(S3, max_features=('unif',0.5,1.0), min_samples_leaf=('int',1,8),
          min_samples_split=('int',2,16))

print("REGIME A: 3 parameters, 64-trial budget (grid gets 4 values each)")
print("-" * 66)
gA = run_grid({"learning_rate":[0.02,0.05,0.11,0.25], "max_depth":[2,3,4,5],
               "subsample":[0.6,0.75,0.9,1.0]})
rA = run_random(S3, 64)
print(f"  grid   {gA[1]:>3} trials -> best CV {gA[0]:.4f}   lr values: {gA[2]['learning_rate']}")
print(f"  random  64 trials -> best CV {rA[0]:.4f}   lr values: {rA[1]['learning_rate']}")
print(f"  random - grid = {rA[0]-gA[0]:+.4f}")

print()
print("REGIME B: 6 parameters, same 64-trial budget (grid gets 2 values each)")
print("-" * 66)
gB = run_grid({"learning_rate":[0.05,0.25], "max_depth":[2,5], "subsample":[0.6,1.0],
               "max_features":[0.5,1.0], "min_samples_leaf":[1,8],
               "min_samples_split":[2,16]})
rB = run_random(S6, 64)
print(f"  grid   {gB[1]:>3} trials -> best CV {gB[0]:.4f}   lr values: {gB[2]['learning_rate']}")
print(f"  random  64 trials -> best CV {rB[0]:.4f}   lr values: {rB[1]['learning_rate']}")
print(f"  random - grid = {rB[0]-gB[0]:+.4f}")

REGIME A: 3 parameters, 64-trial budget (grid gets 4 values each)
------------------------------------------------------------------
  grid    64 trials -> best CV 0.8729   lr values: 4
  random  64 trials -> best CV 0.8721   lr values: 64
  random - grid = -0.0008

REGIME B: 6 parameters, same 64-trial budget (grid gets 2 values each)
------------------------------------------------------------------
  grid    64 trials -> best CV 0.8688   lr values: 2
  random  64 trials -> best CV 0.8728   lr values: 64
  random - grid = +0.0039

The crossover is real and it is in the predicted direction. With three parameters and four grid values each, grid won by 0.0008. With six parameters and two grid values each, random won by 0.0039.

Notice the resolution column, because that is the mechanism working exactly as described: random search tried 64 distinct learning rates in both regimes, while grid tried 4 and then 2. Adding three parameters that barely matter cut grid's resolution on the parameter that matters most in half, and cost it 0.0041 of CV score (0.8729 down to 0.8688) while random search lost almost nothing (0.8721 to 0.8728).

So Bergstra and Bengio are vindicated on mechanism. Now the awkward question: is either gap bigger than noise?

Is Any of This Bigger Than Luck?

Random search has a seed. Grid search does not. So run random search ten times and look at the spread.

import warnings, numpy as np; warnings.filterwarnings("ignore")
from sklearn.datasets import make_friedman1
from sklearn.model_selection import train_test_split, KFold, cross_val_score
from sklearn.ensemble import GradientBoostingRegressor

X, y = make_friedman1(n_samples=800, noise=1.0, random_state=42)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=42)
cv = KFold(3, shuffle=True, random_state=0)

def score(p):
    return cross_val_score(GradientBoostingRegressor(n_estimators=60, random_state=0, **p),
                           Xtr, ytr, cv=cv, scoring='r2', n_jobs=-1).mean()

S3 = {"learning_rate":('logu',0.02,0.25), "max_depth":('int',2,5),
      "subsample":('unif',0.6,1.0)}

def rand_best(spec, budget, seed):
    rng = np.random.RandomState(seed); best = -9
    for _ in range(budget):
        p = {}
        for k, (kind, a, b) in spec.items():
            p[k] = float(np.exp(rng.uniform(np.log(a), np.log(b)))) if kind == 'logu' else \
                   (int(rng.randint(a, b+1)) if kind == 'int' else float(rng.uniform(a, b)))
        v = score(p)
        if v > best: best = v
    return best

GRID_3P = 0.8729   # deterministic, from the previous experiment

print("RANDOM SEARCH, 64 TRIALS, TEN SEEDS  (3 params)")
print("-" * 60)
vals = [rand_best(S3, 64, s) for s in range(10)]
for s, v in enumerate(vals):
    print(f"  seed {s}  best CV {v:.4f}")
print()
print(f"  mean {np.mean(vals):.4f}   min {min(vals):.4f}   max {max(vals):.4f}"
      f"   spread {max(vals)-min(vals):.4f}   sd {np.std(vals):.4f}")
print(f"  grid (3 params) = {GRID_3P:.4f}; "
      f"grid beats {sum(1 for v in vals if v < GRID_3P)} of 10 seeds")

print()
print("COVERAGE ARITHMETIC (why random search is supposed to work)")
print("-" * 60)
for q in (0.10, 0.05, 0.01):
    for n in (30, 64, 200):
        print(f"  P(at least one of {n:>3} draws lands in the top {q:.0%}) = {1-(1-q)**n:.4f}")

RANDOM SEARCH, 64 TRIALS, TEN SEEDS  (3 params)
------------------------------------------------------------
  seed 0  best CV 0.8721
  seed 1  best CV 0.8726
  seed 2  best CV 0.8738
  seed 3  best CV 0.8792
  seed 4  best CV 0.8723
  seed 5  best CV 0.8743
  seed 6  best CV 0.8694
  seed 7  best CV 0.8711
  seed 8  best CV 0.8765
  seed 9  best CV 0.8785

  mean 0.8740   min 0.8694   max 0.8792   spread 0.0098   sd 0.0030
  grid (3 params) = 0.8729; grid beats 5 of 10 seeds

COVERAGE ARITHMETIC (why random search is supposed to work)
------------------------------------------------------------
  P(at least one of  30 draws lands in the top 10%) = 0.9576
  P(at least one of  64 draws lands in the top 10%) = 0.9988
  P(at least one of 200 draws lands in the top 10%) = 1.0000
  P(at least one of  30 draws lands in the top 5%) = 0.7854
  P(at least one of  64 draws lands in the top 5%) = 0.9625
  P(at least one of 200 draws lands in the top 5%) = 1.0000
  P(at least one of  30 draws lands in the top 1%) = 0.2603
  P(at least one of  64 draws lands in the top 1%) = 0.4744
  P(at least one of 200 draws lands in the top 1%) = 0.8660

Random search's spread against itself is 0.0098. The gap it won by in six dimensions is 0.0039. The gap it lost by in three is 0.0008. The noise is two and a half times the larger effect and twelve times the smaller one.

And the head-to-head is the cleanest way to say it: grid beat 5 of 10 random seeds. Exactly five. If I had run one seed and published, this article was a fair coin flip between "the famous result replicates" and "the famous result doesn't," from identical code on identical data at an identical budget.

The coverage table shows why the theory is nonetheless right, and where it stops helping. Sixty-four draws reach the top 10% essentially always (0.9988) and the top 5% almost always (0.9625). But the top 1% is worse than a coin flip at 0.4744, and getting there reliably takes about 200 draws (0.8660). Early in tuning, when you are still finding the right order of magnitude, 64 random draws are plenty. Late in tuning, when the remaining gains live in the top 1% of the space, 64 draws are not a search — they are a sample, and so is a grid.

The Searches from Scratch

Both searches above are hand-rolled — a nested loop over a Cartesian product, and a loop of independent draws. Neither uses GridSearchCV or RandomizedSearchCV, which is deliberate: the whole argument is about how a budget gets spent across axes, and that is easier to see when the loop is visible. Worth confirming the from-scratch grid agrees with the library, though, since every number above rests on it.

import warnings, numpy as np, itertools; warnings.filterwarnings("ignore")
from sklearn.datasets import make_friedman1
from sklearn.model_selection import (train_test_split, KFold, cross_val_score,
                                     GridSearchCV)
from sklearn.ensemble import GradientBoostingRegressor

X, y = make_friedman1(n_samples=800, noise=1.0, random_state=42)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=42)
cv = KFold(3, shuffle=True, random_state=0)
AXES = {"learning_rate":[0.02,0.05,0.11,0.25], "max_depth":[2,3,4,5],
        "subsample":[0.6,0.75,0.9,1.0]}

def score(p):
    return cross_val_score(GradientBoostingRegressor(n_estimators=60, random_state=0, **p),
                           Xtr, ytr, cv=cv, scoring='r2', n_jobs=-1).mean()

mine = max(score(dict(zip(AXES, v))) for v in itertools.product(*AXES.values()))
skl = GridSearchCV(GradientBoostingRegressor(n_estimators=60, random_state=0),
                   AXES, cv=cv, scoring='r2', n_jobs=-1).fit(Xtr, ytr)
print("GRID SEARCH FROM SCRATCH vs SCIKIT-LEARN   (best CV r2)")
print("-" * 56)
print(f"  mine (nested loop over the product) {mine:.4f}")
print(f"  sklearn GridSearchCV                {skl.best_score_:.4f}")
print(f"  trials: mine {4*4*4}, sklearn {len(skl.cv_results_['params'])}")

GRID SEARCH FROM SCRATCH vs SCIKIT-LEARN   (best CV r2)
--------------------------------------------------------
  mine (nested loop over the product) 0.8729
  sklearn GridSearchCV                0.8729
  trials: mine 64, sklearn 64

Identical, as it should be — GridSearchCV is a nested loop with a progress bar and better parallelism. Which is worth saying plainly: there is nothing in either library search that you could not write in six lines, and neither of them does anything with the results of trial 1 before choosing trial 2.

When Does Each Help Most?

	Grid search	Random search
Resolution per parameter	$n^{1/ d}$ , collapses with $d$	$n$ , independent of $d$
Reproducible	exactly, no seed	only with the seed
States what you ruled out	yes	no
Handles 6+ parameters	badly	fine
Resume with more budget	no, the grid is fixed	yes, draw more
Measured, 3 params	0.8729	0.8740 mean, 0.0098 spread
Measured, 6 params	0.8688	0.8728

Reach for grid when you have two or three parameters, when you want an auditable statement about what was eliminated, or when the parameters are genuinely categorical with few levels.

Reach for random when you have four or more parameters, when you want to add budget later without restarting, or when you do not know which parameters matter — which is most of the time.

Reach for neither when the honest answer is that 64 trials cannot resolve your question. Both methods here are sampling blind; neither uses what the previous 63 trials revealed. That is tomorrow's article.

Quick Reference Card

GRID vs RANDOM SEARCH: CHEAT SHEET
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
RESOLUTION PER PARAMETER at n = 64 trials
  grid, 2 params    8 values      random   64
  grid, 3 params    4 values      random   64
  grid, 6 params    2 values      random   64
  grid, 8 params    cannot build  random   64

COVERAGE, random search
  top 10%   30 draws 0.9576   64 draws 0.9988
  top  5%   30 draws 0.7854   64 draws 0.9625
  top  1%   30 draws 0.2603   64 draws 0.4744
            -> the top 1% needs ~200 draws (0.8660)

MEASURED (64 trials, GBM on Friedman #1)
  3 params   grid 0.8729   random 0.8740 mean
  6 params   grid 0.8688   random 0.8728
  random seed spread: 0.0098  <- WIDER than
             either gap (0.0008 and 0.0039)
  grid beat 5 of 10 random seeds

RULES
  ✓ log-uniform for learning rates, regularisation
  ✓ 4+ parameters: use random
  ✓ 2-3 parameters: grid is fine and auditable
  ✓ report the seed spread, never one seed
  ✗ 64 trials cannot resolve the top 1%

SKLEARN
  from sklearn.model_selection import GridSearchCV
  from sklearn.model_selection import RandomizedSearchCV

Key Takeaways

The mechanism is real and visible — random search tried 64 distinct learning rates at every dimension; grid tried 4 at three parameters and 2 at six.
Adding useless parameters costs grid real score — 0.8729 down to 0.8688 when three near-irrelevant parameters halved its resolution, while random went 0.8721 to 0.8728.
The crossover happened where predicted — grid ahead by 0.0008 on three parameters, random ahead by 0.0039 on six.
Both gaps are inside the noise — random search's seed-to-seed spread was 0.0098, with a standard deviation of 0.0030.
Grid beat exactly 5 of 10 random seeds, which is the most honest possible summary of the comparison at this budget.
Coverage explains both the theory and its limit — 0.9988 for the top 10% at 64 draws, but only 0.4744 for the top 1%.
Grid resolution collapses as $n^{1/ d}$ — at 64 trials and 8 parameters the grid cannot be built at all.
The sampling distribution matters more than the method — uniform draws of a learning rate on [0.001, 0.3] land mostly above 0.03, and no grid-versus-random argument repairs that.

The One-Sentence Summary

Random search's advantage over grid search is a real mechanism that behaved exactly as advertised — it held 64 distinct learning rates at every dimension while grid collapsed from 4 values to 2 and lost 0.0041 of score doing it — and yet the head-to-head gaps were 0.0008 and 0.0039 against a seed-to-seed spread of 0.0098, so grid beat five of ten random seeds and the correct conclusion at a 64-trial budget is not that one method wins but that neither method is learning anything from the trials it already ran.

What's Next?

Now that you know why both methods are guessing, you're ready for:

Bayesian optimization — tomorrow. The search that reads its own results before choosing the next trial.
The 5-minute Optuna setup — replacing your grid search tonight.
Nested cross-validation — how to tune and report without burning the same data twice.
Early stopping inside a search — killing bad trials at iteration 10 instead of 200.

Follow me for the next article in the Hyperparameter Tuning series!

Let's Connect!

If the safe at Voutier & Fils made the resolution argument click, drop a heart!

Questions? Ask in the comments — I read and respond to every one.

Have you ever published a tuning comparison from one seed? I nearly did. The first version of this ran seed 0, found grid ahead, and I was three paragraphs into explaining why the famous result fails to replicate before it occurred to me to try seed 1. Seed 3 would have had me arguing the opposite just as confidently. 🔐

References

[1] J. Bergstra and Y. Bengio, Random Search for Hyper-Parameter Optimization (2012), Journal of Machine Learning Research 13

What makes this one worth measuring rather than citing is that the original paper is entirely correct and almost universally over-applied. Bergstra and Bengio proved something specific about resolution in high-dimensional spaces with low effective dimensionality, and the profession compressed it into "random search beats grid search" — a different claim, one that depends on a budget and a loss surface the paper never promised anything about. The paper is not wrong. The sentence it turned into is.

Before your next tuning run, count your parameters — then run the search twice with different seeds and check whether the gap you were about to report survives.

The Hyperparameter Search Cheat Sheet: Which Method, What Budget, and When to Stop

Sachin Kr. Rajput — Mon, 24 Aug 2026 10:05:14 +0000

The One-Line Summary: Every number on this page was measured in an earlier post in this series, and together they say something the folklore doesn't: grid versus random is a coin flip at small budgets (grid beat 5 of 10 random seeds), Bayesian optimization's real product is reliability rather than accuracy (spread 0.0004 against 0.0128), and pruning's advertised speedup is roughly half what the kill rate implies.

Pick By Budget, Not By Fashion

WHICH SEARCH, BY TRIAL BUDGET
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
  < 20 trials    random search
                 a surrogate has nothing to fit;
                 grid can't resolve anything either

  20-30          random, or TPE if trials are slow
                 measured gap at 20: +0.0032, but
                 TPE's own spread was 0.0146

  30-150         Bayesian optimisation (TPE)
                 gap grows with budget:
                 +0.0032 -> +0.0051 -> +0.0063
                 at 20 -> 50 -> 120 trials

  any budget     grid, IF you have 2-3 parameters
                 and need an auditable statement
                 about what you ruled out

  never          grid with 6+ parameters
                 resolution collapses to 2 values
                 per axis and cost 0.0041 of score

The Resolution Argument, In One Table

VALUES TESTED PER PARAMETER at n = 64 trials
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
  params   grid = n^(1/d)      random = n
      2         8                   64
      3         4                   64
      6         2                   64
      8    cannot build             64

MEASURED CONSEQUENCE (GBM, Friedman #1, 64 trials)
  3 params   grid 0.8729   random 0.8740 mean
  6 params   grid 0.8688   random 0.8728

  → adding three near-useless parameters cost
    grid 0.0041 and cost random almost nothing

But Check The Noise Before You Believe Any Of It

RANDOM SEARCH AGAINST ITSELF, 10 SEEDS
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
  min 0.8694   max 0.8792   spread 0.0098
  sd 0.0030    mean 0.8740

  grid (deterministic) 0.8729
  → grid beat 5 of 10 random seeds

  the gaps people argue about:
    3 params  0.0008     6 params  0.0039
  the noise they argue inside:
                       0.0098

✗ A single-seed comparison of search methods is
  a coin flip reported as a finding.

What Bayesian Optimization Actually Buys

TPE vs RANDOM, 5 SEEDS EACH
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
  budget   TPE      random    gap
      20   0.8737   0.8704   +0.0032
      50   0.8781   0.8730   +0.0051
     120   0.8803   0.8741   +0.0063

  → the gap GROWS. random went 0.8730 to 0.8741
    from 50 to 120 trials; TPE gained 0.0022.

FROM SCRATCH (forest surrogate + UCB), 3 seeds
  mine    0.8765 mean   spread 0.0004
  random  0.8750 mean   spread 0.0128

  → 32x tighter, not 32x better.
    On seed 2 random search WON by 0.0048.

The product is reproducibility, not accuracy.

Pruning: The Honest Numbers

OPTUNA MEDIAN PRUNER, 20 TRIALS
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
  no pruner       logloss 0.1908   13.0s   pruned 0
  median pruner   logloss 0.1913    9.4s   pruned 12

  killed 60% of trials -> saved 28% of wall clock
  cost: 0.0005 log loss

✗ The kill rate is not the speedup. A pruned trial
  still pays for its warm-up steps before there is
  enough evidence to stop it.

The Search Space Matters More Than The Search

BEFORE YOU CHOOSE A METHOD
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
✓ log-uniform for learning rates, regularisation,
  anything spanning orders of magnitude
✓ integer-uniform for depths and counts
✓ uniform for fractions (subsample, colsample)

✗ uniform on [0.001, 0.3] for a learning rate puts
  ~90% of draws above 0.03

No amount of grid-versus-Bayesian argument repairs
a space that never samples where the answer is.

COVERAGE, random search (1 - (1-q)^n)
  top 10%   64 draws  0.9988
  top  5%   64 draws  0.9625
  top  1%   64 draws  0.4744   <- coin flip
            200 draws 0.8660

→ early tuning: 64 draws is plenty
→ late tuning: the gains live in the top 1%,
  and 64 draws cannot reach it

Quick Reference Card

HYPERPARAMETER SEARCH: ONE PAGE
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
DECIDE
  2-3 params, want audit trail    -> grid
  4+ params, cheap trials         -> random
  30+ trials, expensive trials    -> TPE
  need reproducible result        -> TPE
  hundreds of idle cores          -> random

TUNE FIRST
  learning_rate x n_estimators (one parameter)
  then max_depth 3-6
  then subsample / colsample
  then stop

ALWAYS
  ✓ log-uniform for rates
  ✓ 8-15 warm-up trials before a surrogate
  ✓ run twice with different seeds
  ✓ report the spread, not the best seed
  ✗ never quote a single-seed method comparison
  ✗ never prune against one validation split

OPTUNA MINIMUM
  study = optuna.create_study(
      sampler=TPESampler(seed=0),
      pruner=MedianPruner(n_warmup_steps=2))
  trial.report(loss, step)
  if trial.should_prune(): raise TrialPruned()

Key Takeaways

Grid versus random is decided by dimension — grid 0.8729 against random 0.8740 on three parameters, and 0.8688 against 0.8728 on six.
At 64 trials neither wins — grid beat exactly 5 of 10 random seeds, and the seed spread of 0.0098 swamps both gaps.
Bayesian optimization compounds — +0.0032, +0.0051, +0.0063 at budgets of 20, 50 and 120.
Its real product is variance reduction — 0.0004 spread against random search's 0.0128.
Pruning's kill rate is not its speedup — 60% of trials killed bought 28% of wall clock, for 0.0005 log loss.
The sampling distribution outranks the method — uniform draws on [0.001, 0.3] put 90% of a learning-rate search above 0.03.
Late-stage tuning needs a different budget — the top 1% of the space is 0.4744 at 64 draws and 0.8660 at 200.

The One-Sentence Summary

Choose your search by trial budget and parameter count rather than by what is fashionable — random below twenty trials, grid only at two or three parameters where the audit trail is worth something, TPE above thirty where its advantage compounds from +0.0032 to +0.0063 — and whatever you choose, run it twice with different seeds, because the spread between two runs of the same method was 0.0098 and every method comparison in this series was smaller than that.

What's Next?

The series recap — the whole tuning arc in one page.
Nested cross-validation — tuning and reporting without burning the same data twice.
Pruning on a CV signal — how to stop pruning against the noise of one split.
Multi-objective tuning — accuracy and latency at the same time.

Follow me for the next article in the Hyperparameter Tuning series!

Let's Connect!

If this saved you a grid search, drop a heart!

Questions? Ask in the comments — I read and respond to every one.

What's your default search now? Mine is TPE with 50 trials and a median pruner, and the honest reason is not the accuracy — it is that I stopped having to wonder whether last week's result was skill or the seed. 📋

Assembling this page was more uncomfortable than writing the articles it summarises, because a cheat sheet wants to say "use X" and every measurement in this series says "it depends, and by less than you think." The one line I would keep if I had to delete the rest: run your search twice with different seeds before you believe anything it told you.

Bookmark this, and check the noise row before you report a tuning win.

Replace Your Grid Search Tonight: The Optuna Setup That Killed 12 of My 20 Trials

Sachin Kr. Rajput — Mon, 24 Aug 2026 10:03:00 +0000

The One-Line Summary: Two extra lines inside the objective — report the score as the model grows, then ask whether to stop — let Optuna abandon 12 of 20 trials before they finished, cutting wall clock from 13.0s to 9.4s at a cost of 0.0005 log loss; the saving is smaller than the kill rate because a pruned trial still runs its warm-up, which is the detail nobody mentions.

The Two Lines That Do The Work

Every grid search has the same defect: it runs every candidate to completion, including the ones that were obviously hopeless after 20% of the work. A human watching a training curve would have killed them. Pruning is that human, automated.

WHAT PRUNING ACTUALLY DOES
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
trial 7, learning_rate=0.011, depth=2

   40 trees   logloss 0.61   <- median so far: 0.28
   80 trees   logloss 0.52   <- still miles behind
                             -> KILLED

  120, 160, 200 trees: never run.

A grid search would have finished this trial, then
politely reported that it was bad.

  trial.report(loss, step)      tell optuna
  if trial.should_prune(): ...  let it decide

That is the entire mechanism. You are handing Optuna a progress signal it can compare across trials, and letting it stop the ones that are behind.

The Whole Setup

Runnable as-is. The only Optuna-specific concepts are suggest_* for the search space, report for progress, and should_prune for the verdict.

import warnings, numpy as np, time; warnings.filterwarnings("ignore")
import optuna; optuna.logging.set_verbosity(optuna.logging.WARNING)
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import log_loss

X, y = make_classification(n_samples=1500, n_features=15, n_informative=6,
                           flip_y=0.03, random_state=42)
Xtr, Xva, ytr, yva = train_test_split(X, y, test_size=0.3, random_state=42, stratify=y)

def objective(trial):
    lr    = trial.suggest_float("learning_rate", 0.01, 0.3, log=True)
    depth = trial.suggest_int("max_depth", 2, 5)
    sub   = trial.suggest_float("subsample", 0.6, 1.0)
    clf = GradientBoostingClassifier(n_estimators=200, learning_rate=lr, max_depth=depth,
                                     subsample=sub, random_state=0, warm_start=True)
    best = 9.9
    for n in range(40, 201, 40):                 # grow the model in steps
        clf.set_params(n_estimators=n); clf.fit(Xtr, ytr)
        loss = log_loss(yva, clf.predict_proba(Xva)[:, 1])
        best = min(best, loss)
        trial.report(loss, n)                    # tell optuna how it's going
        if trial.should_prune():                 # let optuna stop a hopeless trial
            raise optuna.TrialPruned()
    return best

for label, pruner in [("no pruner", optuna.pruners.NopPruner()),
                      ("median pruner", optuna.pruners.MedianPruner(n_startup_trials=5,
                                                                   n_warmup_steps=2))]:
    t0 = time.perf_counter()
    st = optuna.create_study(direction="minimize",
                             sampler=optuna.samplers.TPESampler(seed=0), pruner=pruner)
    st.optimize(objective, n_trials=20, show_progress_bar=False)
    el = time.perf_counter() - t0
    pruned = sum(1 for t in st.trials if t.state == optuna.trial.TrialState.PRUNED)
    print(f"{label:<15} best logloss {st.best_value:.4f}   {el:>6.1f}s   "
          f"pruned {pruned:>2}/20")

no pruner       best logloss 0.1908     13.0s   pruned  0/20
median pruner   best logloss 0.1913     9.4s   pruned 12/20

Reading That Honestly

Pruning killed 12 of 20 trials — 60% of the search — and saved 28% of the wall clock. Those two numbers should not match, and the gap is the useful part: a pruned trial is not free. It ran its warm-up steps before there was enough evidence to kill it, so you pay the first fraction of every bad trial no matter what. n_warmup_steps=2 here means two reports before any trial is eligible, which is the floor on what you can save.

If you see pruning advertised as "60% fewer trials, 60% less time," that is the advertisement, not the measurement.

It cost 0.0005 of log loss — 0.1908 without, 0.1913 with. Pruning can throw away a trial that would have recovered late; some configurations look bad at 80 trees and are excellent at 200. That risk is real and here it cost half a thousandth. Worth it, and worth knowing you paid it.

The trade scales with how expensive your trials are. At 13 seconds for the whole study, none of this matters. At 13 hours it is the difference between one overnight run and three. Pruning is a wall-clock optimisation, so its value is proportional to what a trial costs you.

One thing worth doing differently in a real project: this example prunes against a single validation split, which is what makes it fast enough to be a five-minute demo. In production you want the pruning signal to come from cross-validation, or you are pruning on the noise of one split — and that noise is precisely what the previous two articles were about.

Key Takeaways

Two lines buy the whole mechanism — trial.report(loss, step) and if trial.should_prune(): raise optuna.TrialPruned(). Everything else is the search space you already had.
Killing 60% of trials saved 28% of the time, because warm-up steps run before any trial can be pruned. Never quote the kill rate as the speedup.
It cost 0.0005 log loss (0.1908 to 0.1913) — pruning occasionally discards a late bloomer, and that is the price of the wall clock.

The One-Sentence Summary

Optuna's pruner turns a grid search's worst habit — finishing trials that were already hopeless — into an automated judgment call, and measured here it abandoned 12 of 20 trials to cut 13.0s down to 9.4s for 0.0005 of log loss, with the honest caveat that the 60% kill rate only bought a 28% saving because every pruned trial still pays for its warm-up.

What's Next?

Nested cross-validation — tuning and reporting without burning the same data twice.
Pruning on a CV signal — how to avoid pruning against the noise of one split.
Multi-objective tuning — when you need accuracy and latency.
Search spaces that lie — why uniform sampling of a learning rate is almost always wrong.

Follow me for the next article in the Hyperparameter Tuning series!

Let's Connect!

If this replaces a grid search tonight, drop a heart!

Questions? Ask in the comments — I read and respond to every one.

What's the longest grid search you've ever let finish? Mine was nineteen hours across 486 candidates, and when I finally plotted the curves afterwards, roughly two thirds of them were clearly dead inside the first 15% of their training. ⏱️

The reason pruning feels like cheating is that it is the one optimisation that requires the search to have an opinion mid-trial. Grid and random search treat a trial as atomic: you ask a question, you wait, you get a number. Pruning treats the training curve as evidence in its own right, which is obvious the moment you have ever watched a loss curve and known at epoch three that it was not going to work.

Copy the objective above, swap in your model, and delete your grid.

Blending and Voting: Four Noses, One Bottle, and the Blender Who Graded His Own Homework

Sachin Kr. Rajput — Fri, 21 Aug 2026 07:56:57 +0000

The One-Line Summary: Voting has no learned parameters, so it cannot overfit and it cannot lie to you — but measured here it lost to the single best model on both metrics. Blending learns weights on a holdout and genuinely beat everything, right up until the holdout got small: on 40 rows it reported a log loss of 0.1492, delivered 0.2576, and lost to the free unweighted average it was supposed to improve.

The Parable of the House at Grasse

The house employs four noses, and every batch of oil that arrives has to be judged before the house will put its name on a bottle. The four disagree constantly, which the director considers the point — a panel that always agrees is three salaries wasted. What the house has never settled is how to turn four opinions into one decision.

The Show of Hands

The oldest method is a show of hands. Each nose says ship it or don't, and the majority carries.

THE SHOW OF HANDS
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
Batch 41 arrives.

  Amaury    ship it
  Berthe    ship it
  Cosette   don't
  Didier    ship it

  3 to 1.  The house ships it.

Nobody recorded HOW sure anyone was. Berthe was
certain. Didier nearly abstained. The tally cannot
tell them apart, and it never will.

This is called HARD VOTING.

It is cheap, it is transparent, and it throws away everything except the direction of each opinion.

The Scorecard

So the house moved to scores. Each nose writes a number from 0 to 100, and the panel's verdict is the average.

Better — the average now knows the difference between certainty and a shrug. Except the house discovered something about Didier.

THE SCORECARD
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
Batch 41, scores out of 100:

  Amaury     91   confident, usually right
  Berthe     88   confident, usually right
  Cosette    24   confident, usually right
  Didier     72   ← scores EVERYTHING 68 to 75

  average:  68.8

Didier is not a bad nose. Ask him to rank two
batches and he is excellent. But his numbers have
no range, so averaging drags every verdict toward
the middle of his narrow band.

✗ Averaging opinions also averages in whoever
  cannot express confidence properly.

This is called SOFT VOTING.

Soft voting is usually recommended over hard voting, and the reason is genuine: probabilities carry more information than labels. The measurements below show it is not automatic.

The Director's Ledger

The director's answer was to stop treating the four as equals.

"I have eleven years of batches and eleven years of what happened to them. I do not need a fifth nose. I need to know how much each of these four is worth."

THE DIRECTOR'S LEDGER
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
From the archive, weights fitted to past batches:

  Amaury    0.44
  Berthe    0.00   ← agrees with Amaury almost
  Cosette   0.56     always; adds nothing new
  Didier    0.00   ← no range, no information

  verdict = 0.44 x Amaury + 0.56 x Cosette

Berthe is not fired. She is redundant, which is a
different thing and only a ledger can tell them
apart.

This is called BLENDING.

Two of four weights went to zero. That is the ledger working, not failing.

The Month He Graded His Own Homework

The method held until the director published the house's accuracy in a trade journal.

He had fitted the weights on twelve batches from the previous month. Then, asked how well the panel performed, he reported the panel's score on those same twelve batches. The number was extraordinary. It was also the number he had chosen the weights to maximise.

GRADING YOUR OWN HOMEWORK
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
Twelve batches, four weights to fit.

  score on those twelve batches      excellent
  score on the next twelve           ordinary

  Three batches per weight is not a ledger.
  It is four dials tuned until the last twelve
  results looked good.

The show of hands could never do this. It has
no dials. There is nothing to tune and therefore
nothing to overstate.

This is called OVERFITTING THE BLEND.

The show of hands cannot cheat. That is its one genuine advantage, and it is bigger than it sounds.

Why It Works

THE MATHEMATICS OF AGREEING
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
HARD VOTING    majority of the labels
    0 learned parameters

SOFT VOTING    mean of the probabilities
    0 learned parameters
    → needs members whose confidence MEANS
      something; a flat scorer poisons the mean

BLENDING       weighted mean, weights fitted
    N learned parameters on a held-out set
    → the weights live on the simplex:
      w_i >= 0,  sum w_i = 1
      which is heavy regularisation for free

  THE ONE RATIO THAT MATTERS
    holdout rows / weights
    3 to 10   the reported number is fiction
    ~75       usable, still optimistic
    ~450      the reported number is roughly true

What Are Voting and Blending?

Hard voting takes the majority class label across members. Soft voting averages their predicted probabilities and thresholds the mean. Neither fits anything: both are fixed functions of the members' outputs, so there is no training step and no holdout to reserve.

Blending averages the same probabilities with learned weights. You hold out a slice of data the members never trained on, collect their predictions on it, and fit weights that minimise loss there. It is stacking's simpler cousin — stacking uses out-of-fold predictions across the whole training set, blending uses a single holdout — and it inherits the same discipline problem: the set you fit the weights on is burned, and any number you quote from it is a number you optimised.

The Math

For members $h_{1} \dots h_{N}$ , soft voting is the unweighted mean:

p^(x) = \frac{1}{N} i = 1 \sum N h_{i} (x)

Blending replaces $1/ N$ with fitted weights constrained to the simplex:

p^(x) = i = 1 \sum N w_{i} h_{i} (x), w_{i} \geq 0, i \sum w_{i} = 1

That constraint matters more than it looks. Non-negative weights summing to one cannot extrapolate beyond the members' predictions — the blend is trapped inside their convex hull. It is a strong regulariser you get for free, and it is why blending misbehaves less than its parameter count suggests.

Step by Step

Split off a holdout the members will never train on.
Fit every member on the remaining training data.
Collect member predictions on the holdout.
Fit weights on the simplex to minimise loss on those predictions.
Refit members on everything, keep the weights.
Report on a third set. The holdout is spent.

Voting, Measured

Four deliberately different members on 9,000 rows, split 60/20/20 into train, holdout and test. Voting needs no holdout, so for this first table it just trains and reports.

import warnings, numpy as np; warnings.filterwarnings("ignore")
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import (RandomForestClassifier, GradientBoostingClassifier,
                              VotingClassifier)
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import log_loss, accuracy_score
from sklearn.base import clone

X, y = make_classification(n_samples=9000, n_features=25, n_informative=10,
                           n_redundant=5, flip_y=0.03, random_state=42)
Xtr, Xrest, ytr, yrest = train_test_split(X, y, test_size=0.6,
                                          random_state=42, stratify=y)
Xpool, Xte, ypool, yte = train_test_split(Xrest, yrest, test_size=0.5,
                                          random_state=42, stratify=yrest)

NAMES = ["logreg", "forest", "boosting", "knn"]
FOUR = [LogisticRegression(max_iter=2000),
        RandomForestClassifier(n_estimators=150, random_state=0, n_jobs=-1),
        GradientBoostingClassifier(n_estimators=120, max_depth=3, random_state=0),
        KNeighborsClassifier(n_neighbors=15)]

fit = [clone(m).fit(Xtr, ytr) for m in FOUR]
P = np.column_stack([m.predict_proba(Xte)[:, 1] for m in fit])

print("BASE MODELS ALONE")
print("-" * 46)
print(f"{'':<12}{'accuracy':>10}{'log loss':>11}")
for n, c in zip(NAMES, range(4)):
    print(f"  {n:<10}{accuracy_score(yte,(P[:,c]>0.5).astype(int)):>10.4f}"
          f"{log_loss(yte,P[:,c]):>11.4f}")

est = [(n, clone(m)) for n, m in zip(NAMES, FOUR)]
hard = VotingClassifier(est, voting='hard', n_jobs=-1).fit(Xtr, ytr)
soft = VotingClassifier(est, voting='soft', n_jobs=-1).fit(Xtr, ytr)
print()
print("VOTING  (no learned parameters at all)")
print("-" * 46)
print(f"  hard        {accuracy_score(yte,hard.predict(Xte)):>10.4f}{'      n/a':>11}")
print(f"  soft        {accuracy_score(yte,soft.predict(Xte)):>10.4f}"
      f"{log_loss(yte,soft.predict_proba(Xte)[:,1]):>11.4f}")

BASE MODELS ALONE
----------------------------------------------
              accuracy   log loss
  logreg        0.8156     0.4070
  forest        0.9319     0.2568
  boosting      0.9193     0.2409
  knn           0.9370     0.4728

VOTING  (no learned parameters at all)
----------------------------------------------
  hard            0.9296        n/a
  soft            0.9289     0.2515

Three things here are not what the received advice predicts.

Voting lost to the best single model. KNN alone got 0.9370 accuracy; hard voting 0.9296 and soft voting 0.9289. Boosting alone got 0.2409 log loss; soft voting 0.2515. On both metrics, the ensemble is behind the best member. Averaging is not free — it drags a strong member toward a weak one, and logreg at 0.8156 accuracy and 0.4070 log loss is doing exactly that.

Soft voting did not beat hard voting. 0.9289 against 0.9296. The usual reasoning — probabilities carry more information than labels — is correct in principle and depends on those probabilities being calibrated. KNN with n_neighbors=15 can only emit multiples of 1/15, and logreg's log loss of 0.4728 and 0.4070 say neither is well calibrated here. Soft voting averaged in that miscalibration. Didier, scoring everything 72.

The best member was not the best-calibrated member. KNN wins accuracy (0.9370) and is nearly worst on log loss (0.4728). Boosting wins log loss (0.2409) and is third on accuracy. Which model is "best" is a question about your metric, not your data.

Blending, and What the Holdout Costs

Now the same members with fitted weights. The experiment sweeps holdout size against weight count and reports the gap between what the blend claims on its own tuning set and what it delivers on untouched test data.

import warnings, numpy as np; warnings.filterwarnings("ignore")
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import (RandomForestClassifier, GradientBoostingClassifier,
                              ExtraTreesClassifier, VotingClassifier)
from sklearn.neighbors import KNeighborsClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.tree import DecisionTreeClassifier
from sklearn.metrics import log_loss
from sklearn.base import clone
from scipy.optimize import minimize

X, y = make_classification(n_samples=9000, n_features=25, n_informative=10,
                           n_redundant=5, flip_y=0.03, random_state=42)
Xtr, Xrest, ytr, yrest = train_test_split(X, y, test_size=0.6,
                                          random_state=42, stratify=y)
Xpool, Xte, ypool, yte = train_test_split(Xrest, yrest, test_size=0.5,
                                          random_state=42, stratify=yrest)

FOUR = [LogisticRegression(max_iter=2000),
        RandomForestClassifier(n_estimators=150, random_state=0, n_jobs=-1),
        GradientBoostingClassifier(n_estimators=120, max_depth=3, random_state=0),
        KNeighborsClassifier(n_neighbors=15)]
EXTRA = [ExtraTreesClassifier(n_estimators=150, random_state=0, n_jobs=-1),
         GaussianNB(),
         DecisionTreeClassifier(max_depth=3, random_state=0),
         DecisionTreeClassifier(max_depth=8, random_state=0),
         KNeighborsClassifier(n_neighbors=3),
         KNeighborsClassifier(n_neighbors=50),
         RandomForestClassifier(n_estimators=150, max_depth=5, random_state=1, n_jobs=-1),
         GradientBoostingClassifier(n_estimators=120, max_depth=1, random_state=1)]

f4 = [clone(m).fit(Xtr, ytr) for m in FOUR]
f12 = f4 + [clone(m).fit(Xtr, ytr) for m in EXTRA]
P4p  = np.column_stack([m.predict_proba(Xpool)[:, 1] for m in f4])
P4t  = np.column_stack([m.predict_proba(Xte)[:, 1]  for m in f4])
P12p = np.column_stack([m.predict_proba(Xpool)[:, 1] for m in f12])
P12t = np.column_stack([m.predict_proba(Xte)[:, 1]  for m in f12])

# soft voting from scratch is just the unweighted mean — check it
skl = VotingClassifier([(f"m{i}", clone(m)) for i, m in enumerate(FOUR)],
                       voting='soft', n_jobs=-1).fit(Xtr, ytr)
print("SOFT VOTING FROM SCRATCH vs SCIKIT-LEARN")
print("-" * 52)
print(f"  mine (mean of probabilities) {log_loss(yte, P4t.mean(axis=1)):.4f}")
print(f"  sklearn VotingClassifier     "
      f"{log_loss(yte, skl.predict_proba(Xte)[:,1]):.4f}")

def loss(w, P, yy):
    w = np.abs(w); w = w / w.sum()
    return log_loss(yy, np.clip(P @ w, 1e-9, 1 - 1e-9))

def fit_w(P, yy, k):
    r = minimize(loss, np.ones(k)/k, args=(P, yy), method='Nelder-Mead',
                 options={'maxiter':20000, 'fatol':1e-12, 'xatol':1e-10})
    w = np.abs(r.x); return w / w.sum()

def sweep(Pp, Pt, k, nb, seeds=5):
    ob, ot = [], []
    for s in range(seeds):
        idx = np.random.RandomState(s).choice(len(Xpool), nb, replace=False)
        w = fit_w(Pp[idx], ypool[idx], k)
        ob.append(loss(w, Pp[idx], ypool[idx])); ot.append(loss(w, Pt, yte))
    return np.mean(ob), np.mean(ot)

print()
print("BLENDING OPTIMISM, MEAN OF 5 HOLDOUT DRAWS   (log loss)")
print("-" * 66)
print(f"{'weights':>8}{'holdout':>9}{'on holdout':>12}{'on test':>10}"
      f"{'optimism':>11}{'rows/wt':>9}")
for Pp, Pt, k in [(P4p, P4t, 4), (P12p, P12t, 12)]:
    for nb in (40, 300, 1800):
        b, t = sweep(Pp, Pt, k, nb)
        print(f"{k:>8}{nb:>9}{b:>12.4f}{t:>10.4f}{t-b:>+11.4f}{nb//k:>9}")
    print(f"{k:>8}{'  simple average, no weights, on test':<37}"
          f"{log_loss(yte, Pt.mean(axis=1)):>10.4f}")

SOFT VOTING FROM SCRATCH vs SCIKIT-LEARN
----------------------------------------------------
  mine (mean of probabilities) 0.2515
  sklearn VotingClassifier     0.2515

BLENDING OPTIMISM, MEAN OF 5 HOLDOUT DRAWS   (log loss)
------------------------------------------------------------------
 weights  holdout  on holdout   on test   optimism  rows/wt
       4       40      0.1492    0.2576    +0.1084       10
       4      300      0.1825    0.2253    +0.0428       75
       4     1800      0.1947    0.2121    +0.0174      450
       4  simple average, no weights, on test    0.2515
      12       40      0.1527    0.2190    +0.0663        3
      12      300      0.1856    0.2099    +0.0243       25
      12     1800      0.1870    0.2055    +0.0185      150
      12  simple average, no weights, on test    0.2711

Soft voting from scratch is one line and it lands on the library exactly: 0.2515 against 0.2515. VotingClassifier(voting='soft') is probabilities.mean(axis=1) with paperwork.

Blending genuinely works when the holdout is big. Four weights on 1,800 rows: 0.2121 on test, against 0.2515 for the free average, 0.2409 for the best single member, and 0.2515 for soft voting. Twelve weights on 1,800 rows reached 0.2055, the best number anywhere in this article. So the ledger earns its keep.

And it lies in exact proportion to how little holdout it has. Four weights on 40 rows report 0.1492 and deliver 0.2576. That is not a small overstatement — the reported loss is roughly seventy percent of the real one. Worse, 0.2576 is behind the unweighted average at 0.2515. You reserved data, ran an optimiser, published a number that looked like your best result ever, and shipped something worse than doing nothing.

Optimism falls monotonically as rows-per-weight climbs: +0.1084 at 10 rows per weight, +0.0428 at 75, +0.0174 at 450. That ratio, not the model count, is the thing to watch.

A result I did not predict: twelve weights overfit less than four at every holdout size — +0.0663 against +0.1084 at 40 rows. I expected the opposite. The reason is the simplex constraint: non-negative weights summing to one cannot chase noise very far, and the eight extra members are genuinely diverse enough that the larger blend generalises better despite having three times the parameters. Parameter count is a bad proxy for overfitting risk when the parameters are boxed in.

When Does Each Help Most?

	Hard voting	Soft voting	Blending
Learned parameters	0	0	N
Needs a holdout	no	no	yes
Can overstate itself	no	no	yes
Needs calibrated members	no	yes	helps
Members must be comparable	labels only	probabilities	probabilities
Best measured here	0.9296 acc	0.2515 loss	0.2055 loss
Fails when	members disagree evenly	one member is flat	holdout is thin

Reach for hard voting when members emit labels only, or when you need a result nobody can accuse you of tuning.

Reach for soft voting when members are calibrated and you have no data to spare. Check calibration first — a member with a log loss of 0.4728 will drag the mean.

Reach for blending when you can afford several hundred holdout rows per weight and you are disciplined enough to report on a third set. Below roughly 75 rows per weight, use the unweighted average and keep the data.

Quick Reference Card

VOTING vs BLENDING: CHEAT SHEET
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
WHAT THEY ARE
  hard vote   majority label,      0 parameters
  soft vote   mean probability,    0 parameters
  blending    weighted mean,       N parameters
              fitted on a holdout

MEASURED HERE (9,000 rows, 60/20/20)
  best single member   0.9370 acc / 0.2409 loss
  hard voting          0.9296 acc
  soft voting          0.9289 acc / 0.2515 loss
  blend, 4wt,  1800    0.2121 loss
  blend, 12wt, 1800    0.2055 loss   <- best
  blend, 4wt,  40      claims 0.1492, delivers 0.2576

THE ONE RATIO
  holdout rows / weights
    ~10   optimism +0.1084   fiction
    ~75   optimism +0.0428   usable
   ~450   optimism +0.0174   honest

RULES
  ✓ check member calibration before soft voting
  ✓ report blends on a THIRD set, never the holdout
  ✓ below ~75 rows/weight, just average
  ✗ never assume the ensemble beats its best member

SKLEARN
  from sklearn.ensemble import VotingClassifier
  VotingClassifier(est, voting='soft')

Key Takeaways

Voting has no parameters, so it cannot overfit or overstate — that is its real advantage, and it is why the show of hands could never grade its own homework.
Voting lost to the best single model here — 0.9296 and 0.9289 against KNN's 0.9370. Ensembling is not automatically an improvement.
Soft voting did not beat hard voting — 0.9289 against 0.9296, because it averaged in members with log losses of 0.4728 and 0.4070. Calibration is the precondition, not a detail.
Blending beat everything when the holdout was large — 0.2055 with twelve weights on 1,800 rows, against 0.2409 for the best member.
And it lied by +0.1084 when the holdout was small — reporting 0.1492 while delivering 0.2576, which is worse than the free average at 0.2515.
Watch rows-per-weight, not model count — optimism fell from +0.1084 to +0.0428 to +0.0174 as that ratio went 10, 75, 450.
More weights overfit less, against my expectation — +0.0663 for twelve weights versus +0.1084 for four, because simplex-constrained weights cannot chase noise.
The best member depends on the metric — KNN won accuracy and nearly lost log loss; boosting the reverse. Decide what you are optimising before you decide who is best.

The One-Sentence Summary

Voting is a show of hands that cannot lie because it has nothing to tune, and it lost to the best single model on both metrics here; blending is the director's ledger that genuinely beat everything at 0.2055 when it had 150 holdout rows per weight, and reported a loss thirty percent better than reality when it had ten — so the question is never whether to weight your models but whether you have enough untouched data to earn the weights, and the honest fallback is the unweighted average that costs nothing and claims nothing.

What's Next?

Now that you can tell voting from blending from stacking, you're ready for:

Grid search vs random search — the tourist with a map versus the one who wanders.
Bayesian optimization — the sommelier who learns your palate in six glasses.
Calibration — how to fix the flat scorer that poisoned soft voting above.
Nested cross-validation — how to tune and report without burning the same data twice.

Follow me for the next article in the Ensembles Beyond Bagging series!

Let's Connect!

If the house at Grasse made voting and blending click, drop a heart!

Questions? Ask in the comments — I read and respond to every one.

Have you ever reported a number you'd tuned? I have, and the tell was that it was the best result in the whole project by a suspicious margin — which is exactly what a number looks like when you have optimised it and then quoted it. 🧪

The thing I keep relearning is that the safest method is usually the one with nothing to adjust. An unweighted average has no dials, so it cannot be talked into flattering you, and the price of that honesty is a few points of loss you can measure. Blending has four dials, or twelve, and every one of them is a small opportunity to fit last month's batches and call it skill. The measurements do not say weighting is wrong. They say weighting costs data, and that quoting a weighted result on the data you weighted it against is the same mistake as a director publishing his own warm-up scores.

Send this to whoever on your team is about to fit blend weights on 40 validation rows.

Stacking: The Editor Who Never Reported a Story and Still Wrote the Best One

Sachin Kr. Rajput — Wed, 19 Aug 2026 03:30:00 +0000

The One-Line Summary: Stacking trains a small model to combine your other models, and the entire thing turns on one detail — the combiner must be trained on predictions each base model made for rows it never saw, because when I trained it on in-fold predictions instead it handed a weight of 1.123 to the Random Forest and 0.035 to the Gradient Booster that was genuinely the better model, and the ensemble fell from 0.9331 to 0.8781.

The Parable of the Kestrel Street Newsroom

Four reporters cover the same city for the same paper, and after twenty years everyone in the building knows their habits: one over-reports crime, one is too kind to the council, one files fast and loose, one is slow and exact. For most of the paper's history the editor's job was to decide which of the four to print. Then a new editor arrived who did something stranger — she never reported a story in her life, and the paper got better.

The Old Way: Print the Best Reporter

The old method was a competition. Send all four out on the same story, read the four drafts, print the strongest one and spike the rest.

THE OLD WAY — PICK ONE, BIN THREE
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
Story: how many families left the east ward?

  Ambrose  says 400   (always runs high)
  Beatrix  says 180   (soft on the council)
  Colm     says 260   (fast, sloppy)
  Delia    says 310   (slow, careful)

  Truth, three months later:  295

  Printed: Delia, 310.  The other three: binned.

Delia was closest, so Delia gets printed. Next
story, the paper reads four drafts and throws away
three again.

This is called MODEL SELECTION.

It works. It is also strange: the paper spends four salaries and prints one story.

The Innovation: An Editor Who Reads Habits, Not Stories

The new editor did not go out and report a fifth version. She sat with twenty years of archives and did something nobody had bothered to do — she compared each reporter's draft against what turned out to be true, story after story, until she knew the shape of each one's error.

"I do not need a better reporter. I need to know, for each of these four, which direction they are wrong in and by how much."

THE EDITOR'S LEDGER
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
From the archive, over hundreds of stories:

  Ambrose  runs about +35% high on any count
  Beatrix  runs low whenever the council is
           involved, accurate otherwise
  Colm     is unbiased but wildly noisy
  Delia    is nearly right, slightly low

So for the east ward story she writes:

  0.15 x Ambrose + 0.05 x Beatrix
  + 0.10 x Colm  + 0.75 x Delia   ->  296

  Truth: 295.  Closer than any of the four.

She reported nothing. She corrected for habits.

This is called STACKING.

The editor is not a fifth reporter. She is a model of the other four.

The Year the Editor Read the Wrong Drafts

The method held for a decade and then broke in a way nobody saw coming, because the failure looked like success.

A new archivist, trying to be helpful, started giving the editor the reporters' drafts for stories whose outcome had already been published. Colm — fast, sloppy Colm — had a habit of skimming the paper before filing. On stories where the answer was already in print, his drafts were flawless.

The editor, reading the archive honestly, concluded that Colm was the finest reporter in the building. She started weighting him above everyone. The paper got worse for a year before anyone worked out why.

THE YEAR THE ARCHIVE LIED
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
On stories whose ending was already public:

  Colm's drafts     nearly perfect
  Delia's drafts    still merely very good

  Editor's conclusion:  Colm is the best.
  Editor's new weights: mostly Colm.

On NEW stories, where nothing was published yet,
Colm went back to being fast and sloppy — and the
paper was now built on him.

✗ The editor was never wrong about the archive.
  The archive was the wrong thing to read.

This is called LEAKAGE IN THE META-FEATURES.

The fix is one rule: the editor may only judge a reporter on drafts filed before the outcome was known to anyone. Nothing else about the method changes.

Why It Works

THE MATHEMATICS OF LEARNED WEIGHTS
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
Averaging N models weights them all equally:

    f(x) = (1/N) sum_m  h_m(x)

  → a strong model gets dragged toward weak ones.

Stacking learns the weights from data:

    f(x) = g( h_1(x), h_2(x), ..., h_N(x) )

  where g is fit on OUT-OF-FOLD predictions.

  REQUIREMENTS
  1. Each h_m(x) used to train g must come from a
     model that did NOT see row x.
  2. g should be SIMPLE — it has N inputs and all
     the overfitting risk lives here.
  3. Base models should DISAGREE. Four copies of
     the same model give g nothing to learn.
  4. g may assign negative weight. That is not a
     bug; it is g cancelling a known bias.

What Is Stacking?

Stacking (stacked generalisation) trains a second-level model — the meta-learner — whose input features are the predictions of your first-level models.

The mechanism that makes it work, and the only part people get wrong, is how those input features are generated. For each training row you need a prediction from a model that never saw that row. That is exactly what k-fold cross-validation produces, so the standard recipe is: split into k folds, train each base model on k−1 of them, predict the held-out fold, and stack those held-out predictions into a matrix.

The Math

Given base learners $h_{1}, \dots, h_{N}$ and a meta-learner $g$ , the stacked prediction is:

f (x) = g (h_{1} (x), h_{2} (x), \dots, h_{N} (x))

The meta-learner is fit on the out-of-fold matrix $Z$ , where $Z_{ij} = h_{j (- k (i))} (x_{i})$ and $h_{j (- k (i))}$ is base model $j$ trained with the fold containing row $i$ removed. At prediction time the base models are refit on everything and $g$ is applied unchanged.

Note what $Z$ is not: it is not the base models' training predictions. Using those is the archivist's mistake, and the measurements below put a number on it.

Step by Step

Choose k folds. Five is standard.
For each base model, train on k−1 folds and predict the held-out fold. Repeat until every row has a prediction.
Assemble those into an n × N matrix. This is what the meta-learner sees.
Fit the meta-learner on that matrix against the true targets.
Refit every base model on the full training set.
To predict: run the refit base models, feed their outputs to the meta-learner.

Does It Actually Beat the Best Base Model?

Four base models with genuinely different failure modes, on 3,000 rows of Friedman #1. The comparison that matters is against the best single model, not against the average.

import warnings, numpy as np; warnings.filterwarnings("ignore")
from sklearn.datasets import make_friedman1
from sklearn.model_selection import train_test_split, KFold, cross_val_predict
from sklearn.linear_model import RidgeCV
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.metrics import r2_score

X, y = make_friedman1(n_samples=3000, noise=1.0, random_state=42)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=42)
cv = KFold(5, shuffle=True, random_state=0)

BASE = {
    "Ridge":            lambda: RidgeCV(),
    "RandomForest":     lambda: RandomForestRegressor(n_estimators=200, random_state=0, n_jobs=-1),
    "GradientBoosting": lambda: GradientBoostingRegressor(n_estimators=200, max_depth=3, random_state=0),
    "KNN":              lambda: KNeighborsRegressor(n_neighbors=10),
}

print("BASE MODELS ALONE   (test R2)")
print("-" * 40)
solo = {}
for n, f in BASE.items():
    solo[n] = r2_score(yte, f().fit(Xtr, ytr).predict(Xte))
    print(f"  {n:<18}{solo[n]:>8.4f}")
best_name = max(solo, key=solo.get)
print(f"  best single      : {best_name} at {solo[best_name]:.4f}")

# out-of-fold meta features (correct) vs in-fold meta features (the leak)
oof    = np.column_stack([cross_val_predict(f(), Xtr, ytr, cv=cv, n_jobs=-1) for f in BASE.values()])
infold = np.column_stack([f().fit(Xtr, ytr).predict(Xtr) for f in BASE.values()])
test_meta = np.column_stack([f().fit(Xtr, ytr).predict(Xte) for f in BASE.values()])

meta_oof    = RidgeCV().fit(oof,    ytr)
meta_infold = RidgeCV().fit(infold, ytr)

print()
print("STACKING   (test R2)")
print("-" * 56)
print(f"  simple average of the 4      {r2_score(yte, test_meta.mean(axis=1)):>8.4f}")
print(f"  stacking, out-of-fold meta   {r2_score(yte, meta_oof.predict(test_meta)):>8.4f}")
print(f"  stacking, in-fold meta (LEAK){r2_score(yte, meta_infold.predict(test_meta)):>8.4f}")
print()
print("  what the meta-learner learned (weight per base model)")
for n, w_o, w_i in zip(BASE, meta_oof.coef_, meta_infold.coef_):
    print(f"    {n:<18} out-of-fold {w_o:>7.3f}   in-fold {w_i:>7.3f}")

BASE MODELS ALONE   (test R2)
----------------------------------------
  Ridge               0.7259
  RandomForest        0.8742
  GradientBoosting    0.9301
  KNN                 0.7605
  best single      : GradientBoosting at 0.9301

STACKING   (test R2)
--------------------------------------------------------
  simple average of the 4        0.8696
  stacking, out-of-fold meta     0.9331
  stacking, in-fold meta (LEAK)  0.8781

  what the meta-learner learned (weight per base model)
    Ridge              out-of-fold  -0.040   in-fold  -0.107
    RandomForest       out-of-fold  -0.095   in-fold   1.123
    GradientBoosting   out-of-fold   1.091   in-fold   0.035
    KNN                out-of-fold   0.094   in-fold  -0.027

Three things worth separating.

Simple averaging lost. 0.8696 against 0.9301 for just using the Gradient Booster alone. Averaging a strong model with three weaker ones drags it down — "ensemble everything" is not free, and this is the case people quietly skip when they report ensemble wins.

Stacking won, modestly. 0.9331 against 0.9301. Thirty ten-thousandths. That is a real gain and a small one, which is the honest shape of most stacking results — it is a way to squeeze the last points out, not a step change.

And then the weights, which are the actual story. Read the two columns side by side. Trained on out-of-fold predictions, the meta-learner put 1.091 on the Gradient Booster and gave the others roughly nothing — it correctly identified the best model and used the rest as small corrections. Trained on in-fold predictions, it put 1.123 on the Random Forest and 0.035 on the Gradient Booster.

It picked the wrong model, and it picked it confidently.

The mechanism is worth stating plainly: a Random Forest predicting rows it was trained on is close to perfect, because that is what bagged fully-grown trees do — they memorise. So in the leaked matrix, the Random Forest column looks like an oracle. The meta-learner is not fooled by noise; it is reading an honest record of the wrong thing. That is the archivist's error, in a coefficient.

From Scratch, and Two Robustness Checks

Thirty lines, then two questions that decide whether stacking is worth the operational cost. Smaller sample here — 2,000 rows — so these numbers are not directly comparable to the ones above.

import warnings, numpy as np; warnings.filterwarnings("ignore")
from sklearn.datasets import make_friedman1
from sklearn.model_selection import train_test_split, KFold
from sklearn.linear_model import RidgeCV, LinearRegression
from sklearn.ensemble import (RandomForestRegressor, GradientBoostingRegressor,
                              StackingRegressor)
from sklearn.neighbors import KNeighborsRegressor
from sklearn.metrics import r2_score
from sklearn.base import clone

X, y = make_friedman1(n_samples=2000, noise=1.0, random_state=42)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=42)

def bases():
    return [("ridge", RidgeCV()),
            ("rf",    RandomForestRegressor(n_estimators=100, random_state=0, n_jobs=-1)),
            ("gb",    GradientBoostingRegressor(n_estimators=100, max_depth=3, random_state=0)),
            ("knn",   KNeighborsRegressor(n_neighbors=10))]

def oof_matrix(base, X, y, k=5, seed=0):
    """Each base fits on k-1 folds and predicts the fold it never saw."""
    cvf = KFold(k, shuffle=True, random_state=seed)
    oof = np.zeros((len(X), len(base)))
    for j, (_, est) in enumerate(base):
        for tr, va in cvf.split(X):
            oof[va, j] = clone(est).fit(X[tr], y[tr]).predict(X[va])
    return oof

B    = bases()
OOF  = oof_matrix(B, Xtr, ytr)
FULL = [clone(e).fit(Xtr, ytr) for _, e in B]
TEST = np.column_stack([m.predict(Xte) for m in FULL])

mine = RidgeCV().fit(OOF, ytr)
skl  = StackingRegressor(bases(), final_estimator=RidgeCV(), cv=5, n_jobs=-1).fit(Xtr, ytr)
print("FROM SCRATCH vs SCIKIT-LEARN   (test R2)")
print("-" * 48)
print(f"  mine     {r2_score(yte, mine.predict(TEST)):.4f}")
print(f"  sklearn  {r2_score(yte, skl.predict(Xte)):.4f}")

print()
print("WHICH META-LEARNER?   (same OOF matrix, test R2)")
print("-" * 48)
for name, meta in [("RidgeCV", RidgeCV()),
                   ("LinearRegression", LinearRegression()),
                   ("RandomForest d3", RandomForestRegressor(n_estimators=100, max_depth=3, random_state=0, n_jobs=-1)),
                   ("GradientBoosting", GradientBoostingRegressor(n_estimators=100, max_depth=3, random_state=0))]:
    m = clone(meta).fit(OOF, ytr)
    print(f"  {name:<20}{r2_score(yte, m.predict(TEST)):>8.4f}")

print()
print("DOES A USELESS BASE MODEL HURT?   (test R2)")
print("-" * 48)
rng = np.random.RandomState(1)
const_tr = np.full(len(Xtr), ytr.mean()); const_te = np.full(len(Xte), ytr.mean())
noise_tr = rng.normal(ytr.mean(), ytr.std(), len(Xtr))
noise_te = rng.normal(ytr.mean(), ytr.std(), len(Xte))
for label, otr, ote in [("4 real models", None, None),
                        ("+ constant predictor", const_tr, const_te),
                        ("+ pure noise predictor", noise_tr, noise_te)]:
    O = OOF if otr is None else np.column_stack([OOF, otr])
    T = TEST if ote is None else np.column_stack([TEST, ote])
    print(f"  {label:<24}{r2_score(yte, RidgeCV().fit(O, ytr).predict(T)):>8.4f}")

FROM SCRATCH vs SCIKIT-LEARN   (test R2)
------------------------------------------------
  mine     0.9231
  sklearn  0.9232

WHICH META-LEARNER?   (same OOF matrix, test R2)
------------------------------------------------
  RidgeCV               0.9231
  LinearRegression      0.9231
  RandomForest d3       0.9071
  GradientBoosting      0.9209

DOES A USELESS BASE MODEL HURT?   (test R2)
------------------------------------------------
  4 real models             0.9231
  + constant predictor      0.9231
  + pure noise predictor    0.9229

The reimplementation agrees — 0.9231 against scikit-learn's 0.9232. There is nothing in StackingRegressor beyond the fold discipline.

Keep the meta-learner boring. Ridge and plain linear regression tie at 0.9231. A Random Forest meta-learner drops to 0.9071 and a Gradient Booster to 0.9209. The meta-learner sees four highly correlated columns and a few thousand rows — there is no complex function to find there, only weights, and reaching for capacity costs you.

Garbage base models are nearly free. Adding a constant predictor changed nothing at all (0.9231), and adding a pure noise column cost 0.0002. A linear meta-learner assigns near-zero weight to a column that carries no signal, which is the mirror image of the leak result: stacking is robust to a useless input and catastrophically vulnerable to a dishonest one.

That asymmetry is the practical lesson. Do not agonise over which models to include. Agonise over how the meta-features were made.

When Does Stacking Help Most?

	Stacking helps	Stacking wastes your time
Base models	genuinely different families	four tuned variants of one model
Errors	decorrelated across models	all wrong on the same rows
Gain available	last 1–3 points matter	you need a step change
Data	enough for clean k-fold	small n, folds too thin
Ops budget	can serve N+1 models	latency or memory constrained

The operational cost is the part that decides most real projects. Stacking means training, versioning, monitoring and serving every base model plus the meta-learner. For 0.0030 of R², as measured above, that is often the wrong trade — and it is worth saying so, because the competition write-ups where stacking shines are optimising a leaderboard, not a maintenance burden.

Quick Reference Card

STACKING: CHEAT SHEET
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
WHAT IT IS:
  Train a small model to combine your other models,
  on predictions they made for rows they never saw.

THE ONE RULE:
  meta-features must be OUT-OF-FOLD. Everything
  else is a detail.

MEASURED HERE:
  best single model        0.9301
  simple average           0.8696   (worse!)
  stacking, out-of-fold    0.9331
  stacking, in-fold leak   0.8781

  leaked weights put 1.123 on the WRONG model
  and 0.035 on the right one.

META-LEARNER:
  ✓ Ridge / linear         0.9231
  ✗ RandomForest           0.9071
  keep it simple; it only needs weights.

WHEN TO USE:
  ✓ different model families, decorrelated errors
  ✓ the last 1-3 points are worth N+1 deployments
WHEN NOT EFFECTIVE:
  ✗ base models that agree with each other
  ✗ small n (folds too thin to be honest)
  ✗ tight latency budget

SKLEARN:
  from sklearn.ensemble import StackingRegressor
  StackingRegressor(base, final_estimator=RidgeCV(), cv=5)

Key Takeaways

The out-of-fold rule is the whole method — 0.9331 with it, 0.8781 without, from otherwise identical code.
The leak changes which model you ship — in-fold meta-features put 1.123 on the Random Forest and 0.035 on the Gradient Booster that was actually better.
Because bagged trees memorise — a Random Forest predicting its own training rows looks like an oracle, so it dominates a leaked meta-feature matrix by construction.
Simple averaging can lose to one good model — 0.8696 against 0.9301. Ensembling is not automatically an improvement.
Stacking's honest gain is small — 0.9331 against 0.9301, about 0.0030 of R² for an extra model in production.
Keep the meta-learner linear — Ridge 0.9231, Random Forest 0.9071. It has four inputs and needs weights, not capacity.
Useless base models cost almost nothing — a constant predictor changed nothing and pure noise cost 0.0002. Curate the folds, not the roster.
Thirty lines reproduce the library — 0.9231 against 0.9232. StackingRegressor is fold discipline and nothing else.

The One-Sentence Summary

Stacking is the editor who never reported a story — a small model that learns each of your other models' biases and corrects for them — and it works precisely as well as your fold discipline allows, because the same code that gains you 0.0030 of R² on out-of-fold predictions loses you 0.0550 on in-fold ones and hands the highest weight to whichever model is best at memorising the rows it was trained on.

What's Next?

Now that you understand stacking, you're ready for:

Blending and voting — tomorrow. The two ways a committee can agree, and why one of them cheats.
Grid search vs random search — the tourist with a map versus the one who wanders.
Bayesian optimization — the sommelier who learns your palate in six glasses.
Calibration — when your model's confidence and its accuracy stop agreeing.

Follow me for the next article in the Ensembles Beyond Bagging series!

Let's Connect!

If the newsroom made stacking click, drop a heart!

Questions? Ask in the comments — I read and respond to every one.

Have you ever shipped a stack you couldn't explain? I have, and the honest reckoning was that it beat the best base model by less than the variance between seeds, while tripling what I had to keep alive in production. 📰

What I find genuinely interesting about stacking is that the meta-learner is doing something no base model can do: it is modelling the *models, not the data. And the moment you frame it that way, the leak stops being a technicality and becomes obvious — you cannot learn what a model is like by watching it answer questions it has already been given the answers to. That is true of Random Forests, and it is true of interviews, exams, and every benchmark that quietly ended up in a training set.*

Send this to whoever on your team is about to stack five models before checking whether the fold split is honest.

The Boosting Series in One Page: Ten Posts, and the Six Numbers That Changed How I Tune

Sachin Kr. Rajput — Tue, 18 Aug 2026 10:04:52 +0000

The One-Line Summary: Ten posts on boosting, every number in them produced by code that actually ran — and the six results worth keeping are that the learning rate and the tree count are one parameter, that unlimited depth is worse than stumps, that boosting's optimum was 19 trees and not 3,000, that coarser histogram binning won, that naive target encoding turned a signal-free column into a test loss of 2.7975, and that XGBoost and LightGBM are a tie nobody wants to admit.

The Arc

The Boosting series is finished. It started with a tutor who only reviews the questions you got wrong and ended with a tuning sheet, and the useful thing about writing it in order is that each post kept contradicting something I believed when I started the previous one.

Here is the whole arc, with the one result from each that survived.

AdaBoost: The Tutor Who Only Reviews the Questions You Got Wrong

Where boosting begins: reweight the examples you got wrong, fit the next learner to those, repeat. The parable is a tutor who stops re-teaching what you already know. The mechanism is exponential loss, and it is the reason AdaBoost is so sensitive to a mislabelled row — a permanently wrong answer gets permanently more attention.

Gradient Boosting: The Sculptor Who Fixes the Statue by Carving Only Its Mistakes

The generalisation: stop reweighting examples and start fitting the residual directly, which turns boosting into gradient descent in function space. The from-scratch implementation matched scikit-learn to 0.970898 against 0.970890 — five decimal places from thirty lines of NumPy, which is the clearest evidence I know that the library is not doing anything mysterious.

learning_rate and n_estimators Are One Parameter, Not Two

Sweep the learning rate, early-stop for the tree count each time, and the product barely moves: 4.50, 4.40, 4.41 across a tenfold range of learning rates. Halve the rate and you need twice the trees. This is the single result that has saved me the most compute, because it removes an entire axis from every grid search.

The Boosting Family Decision Tree

Which algorithm, on one page. Worth revisiting now that the series is done, because two of its recommendations got contradicted by later measurements.

XGBoost: The Accountant Who Read the Tax Code So Carefully He Found Money Nobody Else Saw

Second-order gradients, pruning written into the objective, and missing values as a learned direction. The result that landed hardest: second order reached 0.3322 in 50 rounds while first order sat at 0.3436 after 200 and never caught up. From scratch it agreed with the library at 0.931456 against 0.931628.

LightGBM: The Librarian Who Sorted a Million Books by Ignoring Most of Them

Histogram binning and leaf-wise growth. The finding I did not expect: max_bin=15 beat max_bin=1023, 0.1598 against 0.1619. More resolution in the split search bought variance, not accuracy. This is the post that made me stop trusting defaults that sound generous.

CatBoost: The Interpreter Who Refused to Peek at Tomorrow's Newspaper

Target leakage, measured. On a column containing literally no signal, naive target encoding produced a test loss of 2.7975; out-of-fold encoding 0.5823; deleting the column entirely 0.5464. Dropping a useless column beat every clever encoding of it, which is the most uncomfortable result in the series.

Bagging vs Boosting: The Orchestra and the Relay Team

The settled comparison. Removing one tree from a hundred costs bagging 0.0014 and costs boosting 0.2361 at the front of the chain — the crowd is interchangeable, the argument is not. Also: boosting's optimum was 19 trees at 36.1841, and running to 3,000 degraded it to 48.4704. And unlimited depth took boosting to 7.1105, worse than stumps.

I Gave XGBoost and LightGBM the Same 30 Seconds

0.98845 against 0.98825. Two ten-thousandths of an AUC point. Giving LightGBM its speed advantage back as extra search time bought nothing. The libraries differ on training and inference speed, not accuracy, and almost every comparison you will read online reports a gap smaller than its own noise.

The Boosting Tuning Cheat Sheet

All of the above, compressed into a tuning order you can follow without re-reading anything.

The Six That Changed How I Work

learning_rate × n_estimators is one parameter. 4.50, 4.40, 4.41.
max_depth has a low ceiling. Unlimited depth reached 7.1105 — worse than stumps.
Boosting needs early stopping. Optimum 19 trees, not 3,000; 36.1841 against 48.4704.
Coarser can beat finer. 0.1598 against 0.1619 for max_bin 15 versus 1023.
Fix the encoder before the model. 2.7975 naive, 0.5823 out-of-fold, 0.5464 for deleting the column.
Library choice is a tie. 0.98845 against 0.98825 — pick on speed, not accuracy.

What Confused People

Nothing, and I am going to be straight about why rather than invent a reader question.

This series drew zero substantive comments. Two spam replies and one thread I had already answered, across ten posts. That is not a mystery — it is a reach problem I measured rather than guessed at, and the diagnosis was that #datascience had a median of zero reactions across the twenty best posts of its week while #ai had fifty-three, and that I had been spending a quarter of every post's tag budget on a dead slot. That changed from the Bagging vs Boosting post onward.

I mention it because a recap that manufactures "great question from a reader!" is doing the same thing as a benchmark with no error bars. If the questions had been there I would quote them. They were not, so here is the honest version: the series was written into a very quiet room, and the fix was distributional rather than editorial.

The One-Sentence Summary

Ten posts of boosting reduce to one habit — treat the learning rate and the tree count as a single parameter, keep the trees shallow, stop early, and spend the time you save auditing whether your labels and your encoders deserve the trust you are placing in them, because every result in this series that surprised me came from the data and not from the algorithm.

What's Next?

The Boosting arc is closed. Next up is a new series on the ensembles that sit outside the bagging/boosting split:

Stacking — what to do when you have several good models and no idea which to trust.
Blending and voting — the two ways a committee can agree, and why one of them cheats.
Grid search vs random search — the tourist with a map versus the one who wanders.
Bayesian optimization — the sommelier who learns your palate in six glasses.

Follow me for the Ensembles Beyond Bagging series!

Let's Connect!

If any single number in this series changed something you do, drop a heart!

Questions? Ask in the comments — I read and respond to every one, and after the paragraph above you can see I would genuinely like the chance.

Which of the six surprised you most? Mine was the max_bin one, because I had spent two years assuming finer binning was a free accuracy knob you traded memory for. 📉

Writing a series in publication order rather than drafting it all up front turned out to be the useful constraint. Four of the ten posts contradicted something I had asserted in an earlier one, and I only noticed because the code had to run again each time. A series drafted in one sitting would have been internally consistent and wrong in at least four places.

If you're starting a boosting project this week, read the cheat sheet and skip the rest.

The Boosting Tuning Cheat Sheet: Which Knob to Turn First, and When to Stop Turning

Sachin Kr. Rajput — Tue, 18 Aug 2026 09:57:09 +0000

The One-Line Summary: Every number on this page was measured in an earlier post in this series rather than copied from a blog — the tuning order that falls out of them is learning_rate × n_estimators first as one parameter, then max_depth around 3 to 6, then stop, because unlimited depth took boosting from 2.0196 to 7.1105 test MSE and running to 3,000 trees when the optimum was 19 made it a third worse.

The Order

Nobody tunes eleven parameters. Here is the order the measurements justify.

BOOSTING TUNING ORDER
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
1. learning_rate x n_estimators   ← ONE parameter
     pick lr, then early-stop for the tree count.
     Their product is roughly conserved.

2. max_depth  (or num_leaves)     ← 3 to 6
     the single biggest lever after step 1.

3. subsample / colsample          ← 0.6 to 1.0
     mild regulariser, cheap to sweep.

4. everything else                ← usually noise
     stop here unless you have a validation
     curve telling you otherwise.

1. learning_rate and n_estimators Are One Parameter

Measured on classification, sweeping the learning rate and early-stopping for the tree count each time. The product lr × trees barely moves:

lr      peak trees   lr x trees
  0.3           15         4.50
  0.1           44         4.40
  0.03         147         4.41

Three learning rates spanning a factor of ten, and the product lands at 4.50, 4.40, 4.41. Halve the learning rate and you need roughly twice the trees for the same model. So do not grid-search both — fix one, early-stop the other.

2. max_depth Has a Ceiling, and It Is Low

The mirror image of bagging. Test MSE at 100 estimators, same data:

max_depth      bagging     boosting
        1      19.5944       5.0660
        5       5.6241       2.0196
     None       3.6156       7.1105

Bagging improves all the way to unlimited depth. Boosting is U-shaped and unlimited depth is worse than stumps — 7.1105 against 5.0660. A fully grown tree eats the whole residual on the first pass and leaves nothing for the rest of the chain, noise included.

Buy capacity with more trees, not deeper ones.

3. More Trees Is Not Free

Bagging and boosting answer this differently, and it is the most common way people waste a weekend:

n_estimators     bagging     boosting
          10     39.5639      37.1438
         100     36.5054      43.4739
        3000     35.5388      48.4704

Bagging is monotone — 3,000 trees was its best. Boosting's optimum was 19 trees at 36.1841, and by 3,000 it had degraded to 48.4704. A third worse than its own best, from doing more of the thing that was working.

n_estimators is a smoothness knob for bagging and a capacity knob for boosting. Only one of them needs early stopping.

4. The Library-Specific Knobs

LightGBM       measured effect
  max_bin      15 beat 1023  (0.1598 vs 0.1619)
  growth       leaf-wise 0.1753 vs depth-wise 0.1784
  num_leaves   the dangerous default; couple it
               to max_depth or it silently overfits

CatBoost       measured effect
  ordered TS   naive target encoding: test 2.7975
               out-of-fold:            test 0.5823
               drop the column:        test 0.5464
               → on a signal-free column, deleting
                 it beat every encoding of it

The max_bin result is worth pausing on: the coarser binning won. More resolution in the split search is not more accuracy, it is more variance.

5. Which Library — and Why It Barely Matters

Same data, same tuning budget, 15 trials each:

              test AUC     tune    refit   predict
XGBoost        0.98845    10.8s    1.27s     17ms
LightGBM       0.98825     4.9s    0.38s     38ms

Two ten-thousandths of an AUC point. Give LightGBM its speed advantage back as extra search time and it fitted 59 candidates to XGBoost's 37 — final scores 0.98977 and 0.98990. The extra search bought nothing.

So pick on the loop you are constrained by. LightGBM wins development (3.3× faster to fit). XGBoost wins production (2.2× faster to score).

6. When Tuning Is the Wrong Job

BEFORE YOU TUNE, CHECK THESE
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
✗ Labels are noisy
    30% flipped labels cost bagging 0.1133 accuracy
    and boosting 0.1477. Boosting pays more for bad
    labels than any hyperparameter can win back.

✗ A high-cardinality column is target-encoded
    by hand. Fix the encoder, not the depth.

✗ Your CV and your holdout disagree by more than
    the gap you are chasing. Tuning against a
    leaky score optimises the leak.

✓ Clean labels, honest validation, and you still
    want the last two points. Now tune.

Key Takeaways

learning_rate × n_estimators is one parameter — measured at 4.50, 4.40 and 4.41 across a tenfold range of learning rates.
max_depth 3–6, never unlimited — unlimited took boosting from 2.0196 to 7.1105, worse than stumps at 5.0660.
Boosting needs early stopping and bagging does not — optimum at 19 trees versus 3,000, and 48.4704 versus 36.1841 for ignoring it.
Coarser can win — LightGBM's max_bin=15 beat max_bin=1023 at 0.1598 against 0.1619.
Library choice is a tie on accuracy — 0.98845 against 0.98825. Choose on training versus inference speed instead.
Fix the data before the knobs — noisy labels cost boosting 0.1477 of accuracy, which no amount of tuning recovers.

The One-Sentence Summary

Tune the learning rate and tree count as a single parameter with early stopping, keep max_depth between 3 and 6 because unlimited depth is worse than stumps, sweep subsample if you are bored, and spend everything you save on checking whether your labels and your validation split deserve the trust you are placing in them.

What's Next?

The series recap — the whole Boosting arc in one page, with the numbers that survived.
Stacking — what to do when you have several good models and no idea which to trust.
Blending and voting — the two ways a committee can agree, and why one of them cheats.
Grid search vs random search — and why the tourist with a map loses.

Follow me for the next article in the Boosting: The Complete Guide series!

Let's Connect!

If this saved you a grid search, drop a heart!

Questions? Ask in the comments — I read and respond to every one.

What's the hyperparameter you've stopped tuning entirely? Mine is n_estimators on random forests — I set it to 500 and never think about it again, because the measurement says it cannot hurt me. 🔧

The reason cheat sheets are usually wrong is that they are assembled from other cheat sheets. Every figure on this page came from code that ran in an earlier post in this series, which is a low bar that most tuning advice does not clear — and the two most useful lines here, that coarser binning won and that deleting a column beat encoding it, are both things I would not have believed from a table someone else wrote.

Bookmark this before your next grid search, and check the depth row first.