P08 hands-on — Recommender system, block by block

Three of four models lose to popularity. This page is about why.

Source: handson/h08_recsys.py --- run it with python3 handson/h08_recsys.py
Full project spec: P08 — End-to-End Recommendation System

This page is mostly a record of being wrong, which is why it is the most useful of the fifteen.

The plan was straightforward: build matrix factorisation, beat the popularity baseline, show the two-stage architecture. What the measurements said instead was that an under-regularised model scores below popularity, that copying word2vec's negative-sampling constant makes it worse still, that more training makes it worse rather than better, and that the whole question of how much a model can win is decided by a property of the data rather than by the model.

Every one of those findings is in the file, with the experiment that produced it and, where a mechanism was proposed, a prediction tested against a measurement. Block 6 is the clearest case: a theory about p(i|u)/q(i), a prediction it implies, and a verdict of partially confirmed with the residual explained.

Contents

How to read this page

Each block below is a self-contained lego piece: it builds one mechanism, proves it works on its own, and returns what the next block needs. The code is the real source, sliced out of the script. The output underneath it is the real output, captured by running that script --- not transcribed, not idealised. Where a measurement contradicted what I expected, the contradiction is in the output and the prose says so.

The assembly at the end wires every block into one working thing and measures it.

Block 1 — Synthesise a world whose answer you know

Teaches: you cannot debug a recommender on real data

The problem. You cannot debug a recommender on real data, because every bug looks like "the model is bad". Synthesise a world whose ground truth you know, and the same symptom becomes diagnosable.

@block(1, "Synthesise a world whose answer you know", "you cannot debug a recommender on real data")
def b1(s, show):
    NU, ev = generate()
    if show:
        n = sum(len(v) for v in ev.values()); c = counts(ev); top = np.sort(c)[::-1]
        print(f"  {n} interactions, {NU} users, {NI} items, true latent dim {D}")
        print(f"  head mass: top 1% of items = {100*top[:12].sum()/n:>4.1f}% of events;"
              f"  top 10% = {100*top[:120].sum()/n:.1f}%")
        print(f"  {int((c==0).sum())} items ({100*(c==0).mean():.0f}%) never touched -- "
              "the cold tail is the default state")
        print("  Ground truth U, V and the popularity mixture beta are KNOWN here. On")
        print("  MovieLens they are not, so every bug looks like 'the model is bad'.")
    return {"NU": NU, "ev": ev}

Reading the implementation

The generator has one parameter that matters: β, the fraction of choice mass that is pure popularity rather than taste. At β=0.5, half of every user's selections come from a global Zipf popularity distribution and half from their latent affinity. That single knob turns out to determine how much a model can possibly win (block 8), which is why it is a parameter rather than a constant.

Zipf with exponent ~0.9 is not arbitrary — it is the shape of essentially every observed catalogue interaction distribution, and it is why the head-mass numbers below look extreme but are conservative relative to real platforms.

Sampling without replacement per user is a deliberate simplification with a real consequence: as a user's history grows, the remaining candidates are drawn from a depleted head, so very active users have systematically harder held-out items. That artefact is visible if you sweep events-per-user, and it is the kind of thing that would be invisible on real data.

What the numbers say

Output:

  28961 interactions, 800 users, 1200 items, true latent dim 8
  head mass: top 1% of items = 13.8% of events;  top 10% = 39.9%
  0 items (0%) never touched -- the cold tail is the default state
  Ground truth U, V and the popularity mixture beta are KNOWN here. On
  MovieLens they are not, so every bug looks like 'the model is bad'.

The line that matters most is the cold tail: a large fraction of the catalogue is never touched. That is the default state of a catalogue, not an edge case, and any evaluation that only scores items with interaction history is measuring a subset of the problem.

Beyond the toy

Real datasets and what each hides: MovieLens is pre-filtered to users with ≥20 ratings, which removes the cold-start population that dominates production; Amazon reviews are extremely sparse and heavily biased toward extremes; Criteo is the only large public dataset with realistic feature sparsity. Every public benchmark is filtered in a way that makes the problem easier than production, which is why the offline/online gap is not a subtle effect.

Block 2 — The split decides the number

Teaches: a random split leaks the future and inflates everything

The problem. The single largest lever on a reported recommender metric is not the model. It is the split — and a random split leaks the future in a way that inflates every number.

@block(2, "The split decides the number", "a random split leaks the future and inflates everything")
def b2(s, show):
    def split(ev, mode="temporal", seed=3):
        rng = np.random.default_rng(seed)
        tr, te = defaultdict(set), {}
        for u, xs in ev.items():
            j = len(xs) - 1 if mode == "temporal" else int(rng.integers(0, len(xs)))
            te[u] = xs[j]; tr[u] = {x for k, x in enumerate(xs) if k != j}
        return tr, te
    tr, te = split(s["ev"], "temporal")
    rtr, rte = split(s["ev"], "random")
    if show:
        print("  Two splits of the SAME data, one held-out event per user:")
        print("    temporal — hold out each user's LAST event (honest)")
        print("    random   — hold out a uniformly chosen event (leaks the future)")
        print(f"  train sizes are identical ({sum(map(len,tr.values()))} vs "
              f"{sum(map(len,rtr.values()))} pairs), so any difference in score is")
        print("  purely the leak. Block 7 measures how big it is.")
    return {"split": split, "tr": tr, "te": te, "rtr": rtr, "rte": rte}

Reading the implementation

Two splits of identical data, holding out one event per user:

  • Temporal — hold out the user's last event. Honest: the model sees only the past, as it would in production.
  • Random — hold out a uniformly chosen event. The training set then contains events that happened after the one being scored.

The training set sizes are identical, so any difference in score is purely the leak. That control is what makes block 7's inflation measurement meaningful.

The leak is subtler than "the model sees the answer". It sees the user's later behaviour, which reveals their preferences more completely than the past does, and it sees items in a co-occurrence pattern that includes the future. Both make the held-out item easier to predict for reasons that will not exist at serving time.

What the numbers say

Output:

  Two splits of the SAME data, one held-out event per user:
    temporal — hold out each user's LAST event (honest)
    random   — hold out a uniformly chosen event (leaks the future)
  train sizes are identical (28161 vs 28161 pairs), so any difference in score is
  purely the leak. Block 7 measures how big it is.

Beyond the toy

Even temporal-per-user is generous. The strictest and most realistic protocol is a global time split: pick a timestamp, train on everything before, evaluate on everything after. It is harsher because it includes users with no history at all (true cold start) and items that did not exist during training, and it is the only protocol that matches how the model will actually be deployed.

Other leaks that survive a correct split, in rough order of frequency: features computed over the full dataset (a popularity count that includes the test period), negative sampling from the future, hyperparameter tuning on the test set, and early stopping on the test set. The general rule is that anything derived from the whole dataset before splitting is a leak, and normalisation statistics are the classic offender.

Block 3 — The baseline that embarrasses you

Teaches: popularity is not a strawman

The problem. Popularity is not a strawman. It costs one bincount, it requires no training, and it beats a substantial fraction of published models. Any system that does not compute this row does not know whether its model works.

@block(3, "The baseline that embarrasses you", "popularity is not a strawman")
def b3(s, show):
    def evaluate(score, tr, te, k=20):
        rec = ndcg = 0.0
        for u in te:
            order = [i for i in score(u) if i not in tr[u]][:k]
            if te[u] in order:
                rec += 1; ndcg += 1 / np.log2(order.index(te[u]) + 2)
        return rec / len(te), ndcg / len(te)
    pop_order = np.argsort(-counts(s["tr"]))
    r, n = evaluate(lambda u: pop_order, s["tr"], s["te"])
    if show:
        print(f"  popularity, temporal split:   recall@20={r:.4f}   ndcg@20={n:.4f}")
        print(f"  uniform random guessing:      recall@20={20/NI:.4f}")
        print(f"  Popularity is {r/(20/NI):.0f}x random and costs one bincount. Any model")
        print("  that does not clear this line has learned popularity and nothing")
        print("  else -- and you only find out by computing this row.")
    return {"evaluate": evaluate, "pop_order": pop_order}

Reading the implementation

Rank all items by interaction count, exclude what the user has seen, return the top k. That is the whole baseline.

The evaluation function is worth reading carefully because two details decide whether the numbers mean anything:

  • Excluding already-seen items (if i not in tr[u]). Without this, a recommender that returns the user's own history scores brilliantly and is useless.
  • The denominator is \(k\), so a method returning fewer than \(k\) results is correctly penalised — the failure mode of naive post-filtering (P03) and of over-aggressive business-rule filters.

NDCG is reported alongside recall because they answer different questions: recall asks "was it in the list", NDCG asks "how high". A change can improve one and worsen the other, and a system that only tracks one will ship it.

What the numbers say

Output:

  popularity, temporal split:   recall@20=0.1775   ndcg@20=0.0908
  uniform random guessing:      recall@20=0.0167
  Popularity is 11x random and costs one bincount. Any model
  that does not clear this line has learned popularity and nothing
  else -- and you only find out by computing this row.

Popularity is over 10× random. That ratio, not the absolute number, is the reference point for everything that follows — and in the assembly it turns out to be the largest single jump in the entire system.

Beyond the toy

Stronger non-personalised baselines worth beating before claiming a model works: recent-popularity (a trailing window, which handles trends), item-item collaborative filtering (no training, often within a few percent of deep models on public data), and most-popular-in-the-user's-last-category. Dacrema et al. found that most published neural recommenders failed to beat properly tuned versions of these — which is a claim about evaluation discipline, not about neural networks.

Block 4 — BPR: rank, don't predict

Teaches: the loss must match the task, or the metric punishes you

The problem. The loss function decides what the model optimises, and the metric decides how it is graded. When they disagree, the model does exactly what you asked and scores near random.

@block(4, "BPR: rank, don't predict", "the loss must match the task, or the metric punishes you")
def b4(s, show):
    NU, tr, te, ev = s["NU"], s["tr"], s["te"], s["evaluate"]
    P, Q = train_bpr(NU, tr); Pm, Qm = train_mse(NU, tr)
    if show:
        rp, np_ = ev(lambda u: s["pop_order"], tr, te)
        rm, nm = ev(lambda u: np.argsort(-(Pm[u] @ Qm.T)), tr, te)
        rb, nb = ev(lambda u: np.argsort(-(P[u] @ Q.T)), tr, te)
        print(f"  {'model':<28}{'recall@20':>11}{'ndcg@20':>10}{'vs popularity':>15}")
        for lbl, r, n in (("popularity", rp, np_), ("MF + squared error", rm, nm),
                          ("MF + BPR (pairwise)", rb, nb)):
            print(f"  {lbl:<28}{r:>11.4f}{n:>10.4f}{r/rp:>14.2f}x")
        print("  Same architecture, same d, same data, same epochs. Only the loss")
        print("  differs, and squared error lands near random. It asks 'what score?';")
        print("  BPR asks 'which of these two?' -- the question recall@20 grades.")
    return {"P": P, "Q": Q}

Reading the implementation

Two losses, same architecture, same dimensions, same data, same epochs:

  • Squared error treats every observed interaction as a target of 1.0. It asks "what score?", and with only positive examples the trivial optimum is to predict 1.0 for everything — which carries no ranking information at all.
  • BPR samples an unobserved item \(j\) and maximises \(\log \sigma(\hat{x}{ui} - \hat{x}{uj})\). It asks "which of these two?", which is precisely the question recall@k grades.

The gradient makes the difference concrete. BPR's update is \(\sigma(-x_{uij}) \cdot \partial(\hat{x}{ui}-\hat{x}{uj})/\partial\theta\), so correctly ordered pairs produce almost no gradient and the model spends its capacity on the pairs it currently gets wrong. Squared error weights every observation equally regardless of whether the ranking is already right.

Implementation note: the minibatch version computes all three gradients from the batch-start parameters (Pu, Qi, Qj are read before any add.at). That is standard Hogwild-style staleness and is fine at this batch size; it is worth knowing it is a deliberate approximation rather than an oversight.

What the numbers say

Output:

  model                         recall@20   ndcg@20  vs popularity
  popularity                       0.1775    0.0908          1.00x
  MF + squared error               0.0175    0.0137          0.10x
  MF + BPR (pairwise)              0.1900    0.0867          1.07x
  Same architecture, same d, same data, same epochs. Only the loss
  differs, and squared error lands near random. It asks 'what score?';
  BPR asks 'which of these two?' -- the question recall@20 grades.

Squared error lands near random — a 10× gap from the same architecture. This is the clearest demonstration in the track that the loss is not a detail.

Beyond the toy

The ranking-loss family, and when each applies: pointwise (predict a score — appropriate only for explicit ratings), pairwise (BPR, WARP, RankNet — the right default for implicit feedback), and listwise (ListNet, LambdaRank, softmax cross-entropy over the catalogue — optimises the whole list and is what sampled-softmax retrieval models use). WARP is worth knowing specifically: it samples until it finds a violating negative and weights the update by how many attempts that took, which approximates optimising precision@k directly.

Block 5 — Regularisation: the cliff I fell off

Teaches: an under-regularised MF scores BELOW popularity

The problem. This block exists because the first version of this file shipped with reg=0.01 and concluded that matrix factorisation cannot beat popularity. That was not a fact about matrix factorisation.

@block(5, "Regularisation: the cliff I fell off", "an under-regularised MF scores BELOW popularity")
def b5(s, show):
    NU, tr, te, ev = s["NU"], s["tr"], s["te"], s["evaluate"]
    if show:
        rp, _ = ev(lambda u: s["pop_order"], tr, te)
        print(f"  popularity baseline = {rp:.4f}. Recall@20 as training proceeds:")
        print(f"  {'epochs':>8}{'reg=0.01':>11}{'reg=0.05':>11}")
        for epochs in (10, 30, 60, 150):
            row = []
            for reg in (0.01, 0.05):
                P, Q = train_bpr(NU, tr, epochs=epochs, reg=reg)
                row.append(ev(lambda u: np.argsort(-(P[u] @ Q.T)), tr, te)[0])
            print(f"  {epochs:>8}{row[0]:>11.4f}{row[1]:>11.4f}"
                  f"{'   <-- below baseline' if row[0] < rp else ''}")
        print("  This block exists because the first version of this file shipped")
        print("  reg=0.01 and concluded 'MF cannot beat popularity'. It was not a")
        print("  fact about matrix factorisation; it was one hyperparameter. With 32")
        print("  free parameters per user fit from ~35 events, the model memorises")
        print("  the training set and pushes every unobserved item down -- including")
        print("  the held-out one. MORE training makes it WORSE, which is the")
        print("  signature of overfitting and not of a bad architecture.")
    return {}

Reading the implementation

32 free parameters per user, fit from ~35 observations. The model has roughly as many degrees of freedom as data points, so it memorises the training set — and because BPR's objective pushes every unobserved item down, including the held-out one, memorisation actively harms the metric.

The signature is unmistakable once you know it: more training makes it worse. Recall peaks early and declines. That is overfitting, not a bad model class, and the fix is a hyperparameter rather than an architecture.

The general lesson is about attribution. A model that underperforms has many possible causes — wrong loss (block 4), wrong regularisation (here), wrong negative distribution (block 6), or genuinely insufficient signal (block 8) — and they are distinguishable only by controlled experiment. Concluding "MF does not work here" after one configuration is the most common analytical error in applied ML, and it is expensive because it redirects months of effort.

What the numbers say

Output:

  popularity baseline = 0.1775. Recall@20 as training proceeds:
    epochs   reg=0.01   reg=0.05
        10     0.1713     0.1737   <-- below baseline
        30     0.1487     0.1900   <-- below baseline
        60     0.1200     0.1725   <-- below baseline
       150     0.1275     0.1675   <-- below baseline
  This block exists because the first version of this file shipped
  reg=0.01 and concluded 'MF cannot beat popularity'. It was not a
  fact about matrix factorisation; it was one hyperparameter. With 32
  free parameters per user fit from ~35 events, the model memorises
  the training set and pushes every unobserved item down -- including
  the held-out one. MORE training makes it WORSE, which is the
  signature of overfitting and not of a bad architecture.

Beyond the toy

Regularisation in recommenders is unusual in one respect: the right amount depends on per-user support, which varies by orders of magnitude across the user base. A single global \(\lambda\) is simultaneously too strong for heavy users and too weak for light ones. The principled fixes are weighted regularisation (\(\lambda \cdot n_u\), as in ALS-WR), hierarchical Bayesian priors, or simply fewer dimensions for low-support users — which is what mixed-dimension embeddings do at industrial scale for a memory reason and get the regularisation benefit for free.

Block 6 — Negative sampling: a prediction, then a test

Teaches: BPR's optimum ranks by p(i|u)/q(i)

The problem. A theory, a prediction it implies, and a test that could refute it. This block is the 14-step loop compressed into one page — and the verdict is partially confirmed, which is more instructive than a clean win.

@block(6, "Negative sampling: a prediction, then a test", "BPR's optimum ranks by p(i|u)/q(i)")
def b6(s, show):
    if show:
        print("  Theory: BPR with negatives drawn from q converges to a ranking by")
        print("  p(i|u)/q(i) -- the same importance-weighting that makes NCE work.")
        print("  PREDICTION: sampling negatives proportional to popularity DIVIDES OUT")
        print("  the popularity signal. If the truth is popularity-heavy that should")
        print("  be catastrophic; if the truth has no popularity component (beta=0)")
        print("  it should be harmless.\n")
        print(f"  {'beta (popularity mass)':<24}{'uniform q':>11}{'q ~ pop^0.75':>14}"
              f"{'damage':>9}")
        rows = {}
        for beta in (0.5, 0.2, 0.0):
            NU, ev_ = generate(beta=beta)
            tr, te = s["split"](ev_, "temporal")
            P1, Q1 = train_bpr(NU, tr)                 # uniform negatives
            P2, Q2 = train_bpr(NU, tr, alpha=0.75)     # q ~ popularity^0.75
            po = np.argsort(-counts(tr))
            a = s["evaluate"](lambda u: np.argsort(-(P1[u] @ Q1.T)), tr, te)[0]
            b = s["evaluate"](lambda u: np.argsort(-(P2[u] @ Q2.T)), tr, te)[0]
            pr = s["evaluate"](lambda u: po, tr, te)[0]
            rows[beta] = (pr, a, b)
            print(f"  {beta:<24.1f}{a:>11.4f}{b:>14.4f}{b/a:>8.2f}x")
        print("  VERDICT: partially confirmed. The damage shrinks monotonically as the")
        print("  popularity mass falls (0.42x -> 0.68x -> 0.75x), exactly as predicted,")
        print("  but it does not vanish at beta=0. A second mechanism is also present:")
        print("  under q ~ pop, tail items are almost never sampled as negatives, so")
        print("  their embeddings stay near random initialisation and rank spuriously.")
        print("  word2vec uses pop^0.75 because there discounting frequency is the")
        print("  GOAL. Copying the constant into a recommender inverts its purpose.")
        s["beta_rows"] = rows
    return {}

Reading the implementation

Theory. BPR with negatives drawn from distribution \(q\) converges to a ranking by \(p(i|u)/q(i)\) — the same importance weighting that makes noise- contrastive estimation work. The sampling distribution does not merely affect convergence speed; it changes what the optimum is.

Prediction. Sampling \(q \propto \text{pop}^{0.75}\) therefore divides out the popularity signal. If the truth is popularity-heavy (β=0.5) that should be catastrophic; if the truth has no popularity component (β=0) it should be harmless.

Test. Sweep β and measure the damage ratio.

Implementation detail worth noting: the sampler uses cumsum + searchsorted rather than rng.choice(p=...). The latter is \(O(N)\) per call and made the experiment take minutes; the former is \(O(\log N)\) and made it seconds. Same distribution, verified by the results matching exactly.

What the numbers say

Output:

  Theory: BPR with negatives drawn from q converges to a ranking by
  p(i|u)/q(i) -- the same importance-weighting that makes NCE work.
  PREDICTION: sampling negatives proportional to popularity DIVIDES OUT
  the popularity signal. If the truth is popularity-heavy that should
  be catastrophic; if the truth has no popularity component (beta=0)
  it should be harmless.

  beta (popularity mass)    uniform q  q ~ pop^0.75   damage
  0.5                          0.1900        0.0800    0.42x
  0.2                          0.2000        0.1350    0.68x
  0.0                          0.3362        0.2525    0.75x
  VERDICT: partially confirmed. The damage shrinks monotonically as the
  popularity mass falls (0.42x -> 0.68x -> 0.75x), exactly as predicted,
  but it does not vanish at beta=0. A second mechanism is also present:
  under q ~ pop, tail items are almost never sampled as negatives, so
  their embeddings stay near random initialisation and rank spuriously.
  word2vec uses pop^0.75 because there discounting frequency is the
  GOAL. Copying the constant into a recommender inverts its purpose.

Verdict: partially confirmed. The damage shrinks monotonically as popularity mass falls (0.42× → 0.68× → 0.75×), exactly as predicted — but it does not vanish at β=0. A second mechanism is present: under \(q \propto \text{pop}\), tail items are almost never sampled as negatives, so their embeddings stay near random initialisation and can rank spuriously high.

A prediction that lands directionally but incompletely, with the residual explained, is a working model. One that lands exactly is usually a coincidence you have not noticed yet.

Beyond the toy

word2vec uses \(0.75\) because in language modelling, discounting frequency is the goal — you want "the" to stop dominating the objective. Copying the constant into a recommender inverts its purpose, because there popularity is signal rather than nuisance. This is the general failure mode the whole track keeps hitting: importing a result without importing the conditions that made it true.

The correct production practice is mixed sampling (some uniform, some popularity-proportional) with the ratio tuned, or logQ correction — explicitly subtracting \(\log q(i)\) from the logits so the bias is removed analytically rather than by choosing \(q\) carefully. The latter is what large-scale sampled-softmax retrieval models do.

Block 7 — Two-stage: the ceiling you cannot re-rank past

Teaches: stage 2 can only reorder what stage 1 returned

The problem. Serving cannot score a million items per request. The two-stage architecture is forced by latency — and it introduces a hard ceiling that is the first thing to check when quality is bad.

@block(7, "Two-stage: the ceiling you cannot re-rank past", "stage 2 can only reorder what stage 1 returned")
def b7(s, show):
    NU, tr, te, P, Q = s["NU"], s["tr"], s["te"], s["P"], s["Q"]
    pop_order = s["pop_order"]
    def two_stage(u, C, k=20):
        cands = [int(i) for i in pop_order[:C] if i not in tr[u]]   # cheap, no user model
        sc = P[u] @ Q[cands].T                                       # expensive, per-user
        return [cands[j] for j in np.argsort(-sc)][:k]
    if show:
        full = s["evaluate"](lambda u: np.argsort(-(P[u] @ Q.T)), tr, te)[0]
        print(f"  stage 1 = popularity top-C (one bincount, shared by all users)")
        print(f"  stage 2 = the BPR model, scoring only those C items\n")
        print(f"  {'C':>6}{'stage-1 ceiling':>17}{'after re-rank':>15}"
              f"{'scores/user':>13}")
        for C in (20, 50, 200, 600, NI):
            ceil = sum(1 for u in te if te[u] in set(int(i) for i in pop_order[:C])) / len(te)
            hit = sum(1 for u in te if te[u] in two_stage(u, C)) / len(te)
            print(f"  {C:>6}{ceil:>17.4f}{hit:>15.4f}{C:>13}")
        print(f"  single-stage full scan: {full:.4f} using {NI} scores per user")
        print("  Re-ranking never exceeds the ceiling -- it is a hard cap, not a")
        print("  tendency. Before blaming the ranker for a miss, check whether the")
        print("  item was in the candidate set at all. This is also the join to P02:")
        print("  swap popularity top-C for an HNSW query and the ceiling becomes")
        print("  recall@C of the index, which is the number that project measured.")
    return {"two_stage": two_stage}

Reading the implementation

Stage 1 is deliberately made cheap and user-independent (popularity top-C, one bincount shared by all users) so the ceiling actually binds. Stage 2 is the BPR model scoring only those C items.

The measurement is set up to make one point unmissable: re-ranking recall tracks the stage-1 ceiling exactly and can never exceed it. It is a hard cap, not a tendency. Doubling the ranker's quality changes nothing if the item was not in the candidate set.

The practical diagnostic that follows: before investigating a ranking model for missed recommendations, check whether the item was in the candidates. In production this means logging the candidate set, which is expensive and almost always worth it — the alternative is optimising a stage that is not the bottleneck.

What the numbers say

Output:

  stage 1 = popularity top-C (one bincount, shared by all users)
  stage 2 = the BPR model, scoring only those C items

       C  stage-1 ceiling  after re-rank  scores/user
      20           0.1512         0.1512           20
      50           0.2425         0.1850           50
     200           0.4688         0.1875          200
     600           0.7975         0.1900          600
    1200           1.0000         0.1900         1200
  single-stage full scan: 0.1900 using 1200 scores per user
  Re-ranking never exceeds the ceiling -- it is a hard cap, not a
  tendency. Before blaming the ranker for a miss, check whether the
  item was in the candidate set at all. This is also the join to P02:
  swap popularity top-C for an HNSW query and the ceiling becomes
  recall@C of the index, which is the number that project measured.

Beyond the toy

  • Multi-source retrieval is the standard production answer: union several candidate generators (ANN on embeddings, recent-popularity, same-category, collaborative-filtering neighbours, editorial), each covering a different failure mode of the others. Recall is the union's, and no single generator has to be good at everything.
  • The connection to P02 is direct: swap popularity top-C for an HNSW query and the ceiling becomes recall@C of the index, which is exactly the number that project measures. The index's recall and the recommender's recall compose multiplicatively.
  • Latency budget decides C. At 50 ms for ranking and ~50 µs per item scored by a deep model, C ≈ 1000. That arithmetic, not model quality, is what sets the candidate count in most production systems.

Block 8 — How much headroom exists at all

Teaches: the data, not the model, sets the ceiling

The problem. How much can a model possibly win? This block answers it with the same model and the same hyperparameters on three different worlds — and the answer is that the data, not the model, sets the ceiling.

@block(8, "How much headroom exists at all", "the data, not the model, sets the ceiling")
def b8(s, show):
    if show:
        print(f"  {'beta (popularity mass)':<24}{'popularity':>12}{'MF+BPR':>9}"
              f"{'model uplift':>14}")
        for beta in (0.5, 0.2, 0.0):
            pr, mf, _ = s["beta_rows"][beta]
            print(f"  {beta:<24.1f}{pr:>12.4f}{mf:>9.4f}{mf/pr:>13.2f}x")
        print("  Identical model, identical hyperparameters, three worlds. When half")
        print("  the choices are pure popularity there is a 1.07x model to be won;")
        print("  when none are, there is a 6x one. Personalisation uplift is a")
        print("  property of the DOMAIN. Before a quarter of modelling work, estimate")
        print("  the head mass -- it tells you the size of the prize.")
    return {}

Reading the implementation

Identical model, identical hyperparameters, three values of β. The only thing that changes is how much of user behaviour is popularity-driven versus taste-driven.

What the numbers say

Output:

  beta (popularity mass)    popularity   MF+BPR  model uplift
  0.5                           0.1775   0.1900         1.07x
  0.2                           0.1325   0.2000         1.51x
  0.0                           0.0550   0.3362         6.11x
  Identical model, identical hyperparameters, three worlds. When half
  the choices are pure popularity there is a 1.07x model to be won;
  when none are, there is a 6x one. Personalisation uplift is a
  property of the DOMAIN. Before a quarter of modelling work, estimate
  the head mass -- it tells you the size of the prize.

1.07× when half the choice mass is popularity; 6.11× when none of it is. Personalisation uplift is a property of the domain, and it is measurable before any modelling work: estimate the head mass, and you have estimated the size of the prize.

That reframes the assembly's headline. The tuned model winning by 1.07× looks unimpressive until you see the same model win 6.11× on data with no popularity mass. The model was never the limiting factor.

Beyond the toy

Domains ordered roughly by available personalisation uplift: news and trending video (very low — recency and popularity dominate, and personalisation mostly helps with diversity), general e-commerce (moderate), music and long-tail retail (high — taste is idiosyncratic and catalogues are enormous), and dating or job matching (highest, and effectively unsolvable by popularity).

The practical recommendation: spend a day estimating head mass before a quarter of modelling work. It is one bincount, it is the cheapest analysis in this curriculum, and if the answer is that popularity explains most behaviour, the correct engineering decision may be to ship the bincount and work on something else.

The assembly

Every block above, wired together into one working system:

def assembly(s):
    print("\nEight blocks = a recommender. One table, temporal split, honest rows.\n")
    NU, tr, te, ev = s["NU"], s["tr"], s["te"], s["evaluate"]
    P, Q = s["P"], s["Q"]
    Pm, Qm = train_mse(NU, tr)
    Pu, Qu = train_bpr(NU, tr, reg=0.01, epochs=150)
    Pa, Qa = train_bpr(NU, tr, alpha=0.75)
    base = ev(lambda u: s["pop_order"], tr, te)
    rows = [("uniform random", (20/NI, 0.0)),
            ("popularity", base),
            ("MF, squared error", ev(lambda u: np.argsort(-(Pm[u]@Qm.T)), tr, te)),
            ("MF, BPR, reg=0.01, 150ep", ev(lambda u: np.argsort(-(Pu[u]@Qu.T)), tr, te)),
            ("MF, BPR, q ~ pop^0.75", ev(lambda u: np.argsort(-(Pa[u]@Qa.T)), tr, te)),
            ("MF, BPR, tuned", ev(lambda u: np.argsort(-(P[u]@Q.T)), tr, te))]
    hit = sum(1 for u in te if te[u] in s["two_stage"](u, 600)) / len(te)
    rows.append(("  served two-stage, C=600", (hit, float("nan"))))
    print(f"  {'system':<28}{'recall@20':>11}{'ndcg@20':>10}{'vs popularity':>15}")
    for lbl, (r, n) in rows:
        nn = "   -- " if n != n else f"{n:.4f}"
        print(f"  {lbl:<28}{r:>11.4f}{nn:>10}{r/base[0]:>14.2f}x")
    rt = ev(lambda u: np.argsort(-(P[u]@Q.T)), tr, te)[0]
    Pr, Qr = train_bpr(NU, s["rtr"]); rr = ev(lambda u: np.argsort(-(Pr[u]@Qr.T)),
                                              s["rtr"], s["rte"])[0]
    print(f"\n  the same tuned model, scored on the RANDOM split: {rr:.4f} "
          f"({rr/rt:.2f}x)")
    print("  Nothing changed but which event was hidden. A number reported without")
    print("  naming its split is not comparable to anything.")
    print("\n  Read the table top to bottom. The largest single jump is random ->")
    print("  popularity, and it required no model at all. Three of the four MF rows")
    print("  score BELOW that baseline -- one for the wrong loss, one for weak")
    print("  regularisation, one for copying word2vec's sampling constant. The tuned")
    print("  row wins by 1.07x -- an unimpressive number until block 8, where the")
    print("  IDENTICAL model and hyperparameters win 6.11x on data with no popularity")
    print("  mass. The model was never the limiting factor here; the domain was.")
    print("\n  Built: synthetic ground truth -> split discipline -> popularity ->")
    print("  loss choice -> regularisation -> negative-sampling distribution ->")
    print("  two-stage retrieval -> headroom analysis.")
    print("  Missing, on the project page: real MovieLens/Amazon ingest (m1), item")
    print("  and user features for cold start (m7), a served HNSW index in place of")
    print("  popularity top-C (m9, reusing P02), latency budgets under load (m11),")
    print("  and E9 -- the diversity/accuracy tradeoff where recall@20 goes DOWN and")
    print("  the system gets better.")

Output:

Eight blocks = a recommender. One table, temporal split, honest rows.

  system                        recall@20   ndcg@20  vs popularity
  uniform random                   0.0167    0.0000          0.09x
  popularity                       0.1775    0.0908          1.00x
  MF, squared error                0.0175    0.0137          0.10x
  MF, BPR, reg=0.01, 150ep         0.1275    0.0530          0.72x
  MF, BPR, q ~ pop^0.75            0.0800    0.0352          0.45x
  MF, BPR, tuned                   0.1900    0.0867          1.07x
    served two-stage, C=600        0.1900       --           1.07x

  the same tuned model, scored on the RANDOM split: 0.2300 (1.21x)
  Nothing changed but which event was hidden. A number reported without
  naming its split is not comparable to anything.

  Read the table top to bottom. The largest single jump is random ->
  popularity, and it required no model at all. Three of the four MF rows
  score BELOW that baseline -- one for the wrong loss, one for weak
  regularisation, one for copying word2vec's sampling constant. The tuned
  row wins by 1.07x -- an unimpressive number until block 8, where the
  IDENTICAL model and hyperparameters win 6.11x on data with no popularity
  mass. The model was never the limiting factor here; the domain was.

  Built: synthetic ground truth -> split discipline -> popularity ->
  loss choice -> regularisation -> negative-sampling distribution ->
  two-stage retrieval -> headroom analysis.
  Missing, on the project page: real MovieLens/Amazon ingest (m1), item
  and user features for cold start (m7), a served HNSW index in place of
  popularity top-C (m9, reusing P02), latency budgets under load (m11),
  and E9 -- the diversity/accuracy tradeoff where recall@20 goes DOWN and
  the system gets better.

The design space

Recommenders are a pipeline, and each stage has a different cost model. Confusing the stages is the most common architectural mistake.

StageCandidatesLatency budgetModel classWhat it optimises
Retrieval\(10^6\text{--}10^9 \to 10^2\text{--}10^3\)1--10 mstwo-tower, ANN (P02), popularityrecall@C
Filtering\(10^3 \to 10^3\)<1 msbusiness rules, dedup, seen-listcorrectness
Ranking\(10^3 \to 10^2\)10--50 msGBDT, DLRM, cross-attentionpointwise/pairwise accuracy
Re-ranking\(10^2 \to 10\)1--10 msdiversity (MMR, DPP), calibration, business objectivesslate value

The critical property, measured in block 7, is that stage \(n\) cannot exceed stage \(n-1\)'s ceiling. A ranker that is blamed for a miss usually never saw the item. This is the same hard-ceiling relationship as P03's planner and P02's recall@C, and it means the first diagnostic for any recommender quality problem is retrieval recall, not ranking metrics.

Model families, and what each is really for

FamilyExampleStrengthCost
Neighbourhooditem-item CFstrong baseline, interpretable, no training\(O(N^2)\) similarity, cold items excluded
Matrix factorisationALS, BPRdense, fast, good with implicit feedbackno features, cold start fails
Factorisation machinesFM, FFMfeature interactions with shared embeddingsquadratic in fields
Deep + embeddingsDLRM, DCN, Wide&Deeparbitrary features, cross termsembedding tables dominate memory
SequentialGRU4Rec, SASRec, BERT4Recmodels order and intent driftexpensive at serve time
GraphPinSage, LightGCNpropagates signal to cold nodesneighbourhood sampling is the bottleneck
Two-towerYouTube retrievaluser and item encode independently → ANN-ableno cross features, weaker than a ranker

The two-tower/ranker split is forced by latency: a cross-feature model must score each candidate against the user, which is \(O(C)\) forward passes; a two-tower model encodes the user once and the items offline, so retrieval is one ANN query. The architecture is a consequence of the latency budget, not of model quality.

What the blocks actually found

Three of the four models here lose to a bincount, and the reasons generalise:

  1. Loss mismatch. Squared error on implicit feedback scores near random, because it answers "what value?" while recall@k grades "which of these two?".
  2. Regularisation, not architecture. 32 free parameters per user fit from ~35 events memorises the training set and pushes every unobserved item down — including the held-out one. More training makes it worse, which is the signature of overfitting and is easy to mistake for a bad model class.
  3. The negative-sampling distribution changes what is learned. BPR with negatives from \(q\) converges to a ranking by \(p(i|u)/q(i)\), the same importance weighting that makes NCE work. Sampling \(q \propto \text{pop}^{0.75}\) therefore divides out the popularity signal. word2vec uses that exponent because discounting frequency is the goal there; copying the constant into a recommender inverts its purpose.

And the structural finding from block 8: uplift is a property of the domain. Identical model and hyperparameters win 1.07× when half the choice mass is popularity and 6.11× when none of it is. Estimating head mass costs one bincount and tells you the size of the prize before any modelling work.

Latency, memory and the embedding-table problem

Industrial recommenders are dominated by embedding tables, not by compute.

ComponentTypical sizeBound by
Embedding tables100 GB--10 TBmemory capacity and random-access bandwidth
MLP layers10--100 MBcompute
Serving p99 budget50--200 ms end to endthe whole pipeline

An embedding lookup is a random gather: 100 lookups × 128 dims × 4 B = 51 KB, but scattered — so ~100 DRAM round trips at 121 ns ≈ 12 µs of pure latency for one example's features. Batching turns those into parallel gathers, which is why inference batch size matters as much here as in P01, for the same memory-level-parallelism reason as P02.

At scale the tables exceed one machine, and the standard answers are: hashing tricks (mod the ID space into a fixed table, accepting collisions), mixed- dimension embeddings (frequent IDs get more dimensions), quantisation (int8/int4 rows), and sharding across hosts — which turns a lookup into a network round trip and makes the recommender a distributed system.

Hardware note: TPUs have dedicated embedding hardware (SparseCore) precisely because the gather pattern is hostile to a systolic array, and NVIDIA's Merlin/HugeCTR exists for the same reason on GPU. This is one of the few ML workloads where the memory system, not the matmul unit, defines the chip.

Advanced algorithms and evaluation

  • Sampled softmax with logQ correction — the retrieval-side analogue of the \(p/q\) argument above; without the correction the model learns popularity inverted.
  • MMR and DPPs for diversity: a determinantal point process scores a set by the volume its item vectors span, which is the principled version of "do not show ten near-duplicates". Expect recall@k to go down when it works.
  • Calibration: a ranker's scores must be probabilities if downstream business logic (bidding, thresholds) uses them. Isotonic regression or Platt scaling, checked with a reliability diagram, not with AUC.
  • Position bias correction in training labels — see P09; training on raw clicks teaches the model the previous ranker's layout.
  • Evaluation discipline: temporal split, not random. The blocks measure the inflation directly, and the same model scored on a random split looks better by a factor that has nothing to do with users.
  • Offline/online gap. Offline metrics are computed on logged data collected under the old policy, so they systematically favour models that agree with it. This is why P10 exists and why offline recall is a hypothesis rather than a result.

How this connects to the rest of the track

  • P02 is the retrieval index; its recall@C is this system's ceiling.
  • P09 simulates the feedback loop this system creates once deployed.
  • P10 is the only instrument that can tell you whether a model change helped users.
  • P07 computes the real-time features and counters this consumes.
  • P14 explains why embedding lookups are memory-bound and batching is the only lever.

Failure modes at scale

  • Feedback loops — the model trains on data it generated (P09 block 4 makes this visible: 599 items alive → 96 in sixty days).
  • Popularity collapse in the candidate generator, so the ranker only ever sees the head.
  • Training/serving skew: a feature computed one way in the batch pipeline and another way at serve time. The single most common production defect, and the reason feature stores exist.
  • Stale embeddings for new items — cold start is not an edge case; block 1 measures that a large fraction of the catalogue is never touched.
  • Metric gaming: optimising CTR produces clickbait, optimising watch time produces long boring content. The metric is a proxy and the system will find its gap.

Primary sources

  • Rendle et al., BPR: Bayesian Personalized Ranking from Implicit Feedback (UAI 2009) — the loss in block 4.
  • Hu, Koren & Volinsky, Collaborative Filtering for Implicit Feedback Datasets (ICDM 2008).
  • Covington, Adams & Sargin, Deep Neural Networks for YouTube Recommendations (RecSys 2016) — the two-stage architecture.
  • Naumov et al., DLRM (2019) — the embedding-table cost model.
  • Dacrema, Cremonesi & Jannach, Are We Really Making Much Progress? (RecSys 2019) — the paper that showed most reported gains vanish against tuned baselines, which is what blocks 3--5 reproduce in miniature.
  • Chen et al., Bias and Debias in Recommender System: A Survey (2020).

Running it

python3 handson/h08_recsys.py            # every block, then the assembly
python3 handson/h08_recsys.py --block 3  # just block 3 and its prerequisites
python3 handson/h08_recsys.py --quiet    # the assembly only

What to do with this

Run block 8 on your own data. Estimate the head mass --- what fraction of interactions go to the top 1% of items --- before writing any modelling code. It tells you the size of the prize, and it is the single cheapest analysis in this entire curriculum. If the answer is that popularity explains most of the behaviour, the correct engineering decision may be to ship the bincount.

Milestones, experiments, readings and exit criteria for this project: P08 — End-to-End Recommendation System.