Recommendation, Ranking, and Personalization Questions

Approach

Maximal marginal relevance (MMR) builds the list one item at a time. The first pick is the most relevant candidate. Every later pick maximises a blend of "how relevant is this item" and "how different is it from what I already picked":

• $S$ is the list built so far and $i$ ranges over candidates not yet in it.
• $\tilde r(i)$ is the relevance score of item $i$ rescaled to [0, 1] within the request, so that $\lambda$ means the same thing for every request.
• $\mathbf{e}_i$ is the embedding (a vector of numbers describing the item); $\cos$ is cosine similarity, the dot product of the two vectors after each is divided by its length.
• $\lambda = 1$ ignores diversity and returns the candidates in descending score order, which is the original ranking when the input is already ranked best first (as in the example below). $\lambda = 0$ ignores relevance after the first pick.

The key efficiency step: instead of recomputing the maximum over $S$ for every candidate at every step, keep one array max_sim holding each candidate's best similarity to anything chosen so far. When a new item nxt is added, update it with np.maximum(max_sim, sim[nxt]), which costs $O(n)$.

One step by hand

Three candidates: A (score 0.90), B (score 0.85, nearly the same item as A, cosine similarity 0.95) and C (score 0.60, unrelated to A, similarity 0). A is picked first because it is the most relevant. Rescaled to [0, 1] the relevances are A = 1, B = 0.833, C = 0. For the second pick the program evaluates the formula for B and C at two settings of $\lambda$:

import numpy as np

scores = {"A": 0.90, "B": 0.85, "C": 0.60}
lo, hi = min(scores.values()), max(scores.values())
rel = {k: (v - lo) / (hi - lo) for k, v in scores.items()}   # rescale to [0, 1]
sim_to_A = {"B": 0.95, "C": 0.0}                               # the only chosen item is A
for lam in (0.7, 0.4):
    vals = {k: lam * rel[k] - (1 - lam) * sim_to_A[k] for k in ("B", "C")}
    print(f"lam={lam}:", {k: round(v, 3) for k, v in vals.items()}, "-> second pick", max(vals, key=vals.get))

Output:

lam=0.7: {'B': 0.298, 'C': 0.0} -> second pick B
lam=0.4: {'B': -0.237, 'C': 0.0} -> second pick C

At $\lambda = 0.7$, B still wins ($0.7 \times 0.833 - 0.3 \times 0.95 = 0.298$ against 0 for C): it is relevant enough to pay its redundancy penalty. At $\lambda = 0.4$ the penalty weighs more ($0.4 \times 0.833 - 0.6 \times 0.95 = -0.237$) and C wins. That flip is what the full program shows on 12 items.

Code

If only ranks and no scores are available, pass -rank or 1/rank as the score. Run with Python 3.12, numpy 2.1.3, the program below defines the function, a second slow version that recomputes everything from the definition, and the checks, and prints:

import numpy as np

def mmr_rerank(scores, embeddings, k, lam=0.7):
    """Maximal marginal relevance. Returns the indices of the k chosen items, in order.

    scores:     (n,) relevance from the ranker, higher is better
    embeddings: (n, d) item vectors; cosine similarity is used
    lam:        1.0 = pure relevance, 0.0 = pure diversity
    Complexity: O(n^2 d) for the similarity matrix, then O(k n) for the greedy loop.
    """
    scores = np.asarray(scores, dtype=float)
    X = np.asarray(embeddings, dtype=float)
    n = len(scores)
    if n == 0 or k <= 0:
        return []
    k = min(k, n)
    norms = np.linalg.norm(X, axis=1, keepdims=True)
    X = X / np.where(norms == 0, 1.0, norms)       # zero vectors stay zero: similarity 0
    sim = X @ X.T                                  # (n, n) cosine similarity
    # put relevance on [0, 1] so lam has the same meaning on every request
    span = scores.max() - scores.min()
    rel = (scores - scores.min()) / span if span > 0 else np.zeros(n)

    chosen = [int(np.argmax(rel))]                 # argmax returns the lowest index on ties
    max_sim = sim[chosen[0]].copy()                # best similarity to anything chosen so far
    available = np.ones(n, dtype=bool)
    available[chosen[0]] = False
    while len(chosen) < k:
        mmr = lam * rel - (1 - lam) * max_sim
        mmr[~available] = -np.inf
        nxt = int(np.argmax(mmr))
        chosen.append(nxt)
        available[nxt] = False
        max_sim = np.maximum(max_sim, sim[nxt])    # incremental update, O(n) per step
    return chosen

def mmr_reference(scores, embeddings, k, lam, agg=max):
    """Plain-Python version that recomputes everything from the definition at every step."""
    n = len(scores)
    lo, hi = min(scores), max(scores)
    rel = [(s - lo) / (hi - lo) if hi > lo else 0.0 for s in scores]
    def cos(a, b):
        na, nb = np.linalg.norm(a), np.linalg.norm(b)
        return 0.0 if na == 0 or nb == 0 else float(a @ b / (na * nb))
    chosen = [max(range(n), key=lambda i: (rel[i], -i))]
    while len(chosen) < min(k, n):
        cands = [i for i in range(n) if i not in chosen]
        def val(i):
            return lam * rel[i] - (1 - lam) * agg([cos(embeddings[i], embeddings[j]) for j in chosen])
        chosen.append(max(cands, key=lambda i: (val(i), -i)))
    return chosen

# ---- driver ----
rng = np.random.default_rng(0)
n, d = 12, 4
scores = np.sort(rng.uniform(0.2, 1.0, n))[::-1]          # already ranked, best first
topic = np.array([0, 0, 0, 0, 1, 1, 1, 2, 2, 3, 3, 3])    # the top 4 are the same topic
centres = rng.normal(size=(4, d))
emb = centres[topic] + 0.15 * rng.normal(size=(n, d))

print("topic of each item, in ranked order:", topic.tolist())
for lam in [1.0, 0.7, 0.4, 0.0]:
    out = mmr_rerank(scores, emb, 5, lam)
    print(f"lam={lam}: items {out}  topics {topic[out].tolist()}")

# 1) lam=1 must be the original ranking
assert mmr_rerank(scores, emb, 5, 1.0) == list(range(5))
# 2) vectorised version equals the from-the-definition version on 200 random problems
agree = 0
for t in range(200):
    r = np.random.default_rng(t)
    s, e = r.uniform(size=15), r.normal(size=(15, 3))
    lam = float(r.choice([0.2, 0.5, 0.8]))
    agree += mmr_rerank(s, e, 6, lam) == mmr_reference(list(s), e, 6, lam)
print("agree with reference on 200 random problems:", agree)
# 3) the check can fail: a mean-similarity penalty is a different algorithm
diff = sum(mmr_reference(list(s), e, 6, lam, agg=lambda v: sum(v) / len(v)) != mmr_rerank(s, e, 6, lam)
           for s, e, lam in [(np.random.default_rng(t).uniform(size=15), np.random.default_rng(t).normal(size=(15, 3)), 0.5) for t in range(50)])
print("mean-similarity variant differs from MMR on", diff, "of 50 problems")
# 4) edge cases
print("k > n:", mmr_rerank(scores[:3], emb[:3], 10, 0.5))
print("empty:", mmr_rerank([], np.zeros((0, 2)), 5))
print("all-equal scores:", mmr_rerank(np.ones(4), [[1, 0], [1, 0], [0, 1], [1, 1]], 3, 0.5))
print("zero vector item, k=2:", mmr_rerank([0.9, 0.8, 0.7], [[1, 0], [0, 0], [1, 0.1]], 2, 0.5))

Output:

topic of each item, in ranked order: [0, 0, 0, 0, 1, 1, 1, 2, 2, 3, 3, 3]
lam=1.0: items [0, 1, 2, 3, 4]  topics [0, 0, 0, 0, 1]
lam=0.7: items [0, 4, 1, 7, 2]  topics [0, 1, 0, 2, 0]
lam=0.4: items [0, 5, 9, 7, 1]  topics [0, 1, 3, 2, 0]
lam=0.0: items [0, 9, 5, 8, 6]  topics [0, 3, 1, 2, 1]
agree with reference on 200 random problems: 200
mean-similarity variant differs from MMR on 46 of 50 problems
k > n: [0, 1, 2]
empty: []
all-equal scores: [0, 2, 3]
zero vector item, k=2: [0, 1]

Reading the code, line by line where it is not obvious:

• np.argmax(rel) returns the index of the largest value (the lowest index if several tie); it picks the first item. np.linalg.norm(X, axis=1) is each row's length, and dividing a vector by its length gives a unit vector, so the dot product of two unit vectors is exactly their cosine similarity. np.where(norms == 0, 1.0, norms) swaps a zero length for 1 before dividing, so a zero vector stays all zeros instead of becoming NaN (0/0), which makes its similarity to everything 0.
• sim = X @ X.T multiplies the unit-vector table by its transpose, giving an $n \times n$ table whose entry $(i, j)$ is the cosine similarity of items $i$ and $j$.
• Min-max rescaling, $(s - s_{\min}) / (s_{\max} - s_{\min})$, maps the lowest score to 0 and the highest to 1; if all scores are equal the span is 0, so the code uses zeros.
• max_sim = sim[chosen[0]].copy() starts as the first pick's row because, with one item chosen, the best similarity to the chosen set is just the similarity to that item. After each new pick np.maximum(max_sim, sim[nxt]) keeps, for every candidate, the larger of its old best and its similarity to the new item (the incremental update: one $O(n)$ pass instead of a fresh maximum over the whole chosen set).
• mmr[~available] = -np.inf gives already chosen items the lowest possible value so argmax can never select them again.

What the output shows:

• The top 4 items all share topic 0. At $\lambda = 1$ the list is items 0 to 4, four items from topic 0 and one from topic 1. At $\lambda = 0.7$ the second pick jumps to item 4 (topic 1) and the list holds three topics; at $\lambda = 0.4$ it holds four topics.
• The first check (an assert) proves the $\lambda = 1$ limit is the original ranking.
• The vectorised function agrees with the from-the-definition reference on all 200 random problems. The reference uses a different code path (Python loops, recomputing the maximum over chosen items at each step), so agreement is evidence the incremental update is right.
• The line beginning "mean-similarity variant" shows the comparison can fail: replacing the maximum similarity by the mean similarity (a different algorithm) disagrees with MMR on 46 of 50 problems, so the "agree with reference" line is not trivially true.
• Edge cases: $k > n$ returns all $n$ items; an empty candidate list returns []; with equal scores relevance drops out and diversity decides: the printed result is [0, 2, 3], so item 1, an exact duplicate of item 0, is left out, and any remaining tie goes to the lowest index; a zero vector (an item with no embedding) is treated as having similarity 0 to everything instead of producing NaN: in the zero-vector line, item 1 (the zero vector) is chosen second, ahead of item 2, which is almost a copy of item 0 and is penalised for it.

Complexity

• Similarity matrix: $O(n^2 d)$ time and $O(n^2)$ memory for $n$ candidates and embedding size $d$.
• Greedy loop: $k$ steps, each $O(n)$, so $O(kn)$.
• If $n$ is large, only compute similarities between candidates and the chosen items: $O(k n d)$ time and $O(nd)$ memory (the embeddings themselves plus one length-$n$ array), because each step needs only the row of the newly chosen item.

Edge cases and pitfalls

• Scale. Relevance and similarity must be on comparable scales, which is why the code min-max normalises the scores. If the ranker outputs a calibrated probability, you may skip the normalisation, but then $\lambda$ depends on the score distribution.
• Constant scores. If every score is equal, the code sets relevance to zero for all items and the pick order is driven by diversity alone.
• Duplicates. Two items with identical embeddings have similarity 1, so after one is chosen the other is penalised by $(1-\lambda)$ in full.
• Negative cosine. Items pointing in opposite directions have similarity near -1, which gives them a bonus; clip similarities at 0 if "opposite" should not count as "very different".
• Greedy is a heuristic. It can miss the best set; do not present the output as the optimum.
• Mutating input. The function allocates new arrays and does not modify its arguments.

Running the code

Save the listings above to files and run each with python3. Python 3.12 and numpy 2.1.3 are all that is needed.

Direct answer

Ask each source for its own ranked list, convert every list to a rank-based score (so sources with different score scales become comparable), add up the weighted scores per track, de-duplicate on a canonical track id (so a remaster or live version of a track already proposed or heard counts as the same track), apply diversity limits, and return the best $N$. The weights say how much you trust each source; the sum rewards tracks that several sources agree on. The function never pads with weak filler when sources run dry, and a source that fails contributes nothing instead of failing the request.

Approach

1. Sources. The artist index returns tracks by artists the user follows, most popular first. The co-listen index returns tracks that other listeners played alongside this user's tracks. Embedding similarity returns tracks whose embeddings (vectors, meaning lists of numbers, learned as a description of each track) point in a similar direction to the user's vector, measured by cosine similarity (the cosine of the angle between two vectors).

2. Rank-based scoring with reciprocal rank fusion (RRF). Each source contributes $w_s / (k + \text{rank}_s)$ for a track at position $\text{rank}_s$ in source $s$'s list, where $w_s$ is the source weight, ranks start at 1 for the first item in a source's list, and $k$ is a smoothing constant (10 here). The weight multiplies the whole term, so tripling a source's weight triples every contribution it makes. A track's score is the sum over sources. Raw scores (co-listen counts, cosine values, popularity) are on incompatible scales; ranks are not.

   The constant $k$ controls how steeply the first positions outrank the next ones. With $k = 10$ the rank-1 term $1/11$ is only about 1.09 times the rank-2 term $1/12$, so a track that three sources all place mid-list can beat a track that one source puts first. With $k = 1$ the ratio is $(1/2)/(1/3) = 1.5$, and a single first place counts for much more. The program below prints both ratios.

3. De-duplicate in layers.

   • Exact id: the same track from two sources merges into one entry whose score is the sum, and the provenance (the record of which sources proposed it) is kept for logging.
   • Canonical id: a lookup maps remasters, live versions and re-releases to one canonical track; dedupe and the "already heard" filter both use the canonical id.
   • Near-duplicates not in the lookup (the same song under a different id) can be caught by an embedding-similarity threshold, which needs tuning because it also merges legitimate covers.

4. Exclusions and diversity. Drop tracks the user already heard. Cap the number of tracks per artist (2 here) so one source or artist cannot fill the list.

5. Deterministic output. Sort by score descending, then canonical id, so ties never vary between calls.

Code

The catalogue is 40 synthetic tracks over 8 artists; tracks 30 to 33 are remasters of tracks 2, 5, 9 and 12, and tracks 0, 1 and 2 are already heard (the last demo lines also try a user who heard a remaster). The embedding source is an exact cosine search over random 8-dimensional vectors, a stand-in for an ANN (approximate nearest-neighbour) index, which finds close vectors quickly by not comparing against every item; the artist and co-listen lists are fixed stand-ins for real index lookups. The asserts check the stated properties; the cosine used is compared against scikit-learn's.

import numpy as np
from collections import defaultdict
from sklearn.metrics.pairwise import cosine_similarity

rng = np.random.default_rng(0)
N_TRACKS = 40
track_artist = {t: f"artist_{t % 8}" for t in range(N_TRACKS)}
# tracks 30..33 are remasters / live versions of tracks 2, 5, 9, 12
canonical = {t: t for t in range(N_TRACKS)}
canonical.update({30: 2, 31: 5, 32: 9, 33: 12})
emb = rng.normal(0, 1, (N_TRACKS, 8))
user_vec = rng.normal(0, 1, 8)
heard = {0, 1, 2}

# --- three sources, each returns an ordered list of track ids (best first)
def artist_source(user, limit):          # tracks by followed artists, most popular first
    return [t for t in range(N_TRACKS) if track_artist[t] in {"artist_1", "artist_2"}][:limit]

def co_listen_source(user, limit):       # tracks other listeners played with this user's tracks
    return [5, 31, 9, 14, 3, 17, 20, 33][:limit]

def embedding_source(user, limit):       # exact cosine search; a real system would use an ANN index
    sims = (emb @ user_vec) / (np.linalg.norm(emb, axis=1) * np.linalg.norm(user_vec))
    return [int(t) for t in np.argsort(-sims, kind="stable")[:limit]]

SOURCES = {"artist": artist_source, "co_listen": co_listen_source, "embedding": embedding_source}


def generate_candidates(user, n, weights, heard, per_artist_cap=2, rrf_k=10, per_source_limit=20):
    """Weighted reciprocal-rank fusion over sources, de-duplicated on canonical track id."""
    if n <= 0:
        return []
    heard_c = {canonical.get(h, h) for h in heard}   # a heard remaster blocks the original, and the reverse
    score = defaultdict(float)
    prov = defaultdict(set)              # which sources proposed each canonical track
    rep = {}                             # canonical id -> representative track id to serve
    for name, fn in SOURCES.items():
        w = weights.get(name, 0.0)
        if w <= 0:
            continue
        try:
            ranked = fn(user, per_source_limit)
        except Exception:                # a failing source must not fail the request
            ranked = []
        seen = set()                     # a source may list a track and its remaster: count the better rank once
        for rank, t in enumerate(ranked, start=1):
            c = canonical[t]
            if c in heard_c or c in seen:
                continue
            seen.add(c)
            score[c] += w / (rrf_k + rank)
            prov[c].add(name)
            rep.setdefault(c, c)         # always serve the canonical (original) version
    out, per_artist = [], defaultdict(int)
    for c in sorted(score, key=lambda c: (-score[c], c)):          # deterministic tie-break
        a = track_artist[rep[c]]
        if per_artist[a] >= per_artist_cap:
            continue
        per_artist[a] += 1
        out.append((rep[c], round(score[c], 4), sorted(prov[c])))
        if len(out) == n:
            break
    return out


res = generate_candidates("u1", 8, {"artist": 1.0, "co_listen": 1.0, "embedding": 1.0}, heard)
for t, s, p in res:
    print(f"track {t:2d} {track_artist[t]:9s} score {s:.4f} from {p}")

def trace(c):                            # one track's fused score, term by term (ranks start at 1)
    total = 0.0
    for name, fn in SOURCES.items():
        pos = [r for r, t in enumerate(fn("u1", 20), start=1) if canonical[t] == c]
        if pos:
            term = 1.0 / (10 + pos[0]); total += term
            print(f"  {name:9s} rank {pos[0]:2d}: 1/(10+{pos[0]}) = {term:.4f}")
    print(f"  track {c} total {total:.4f}")
trace(17); trace(9)
print("rank 1 vs rank 2 term ratio, k=10:", round((1 / 11) / (1 / 12), 3), " k=1:", round((1 / 2) / (1 / 3), 3))

ids = [t for t, _, _ in res]
assert len(ids) <= 8 and len(set(canonical[t] for t in ids)) == len(ids)       # no duplicates incl. remasters
assert not (set(ids) & heard)                                                 # nothing already heard
cnt = defaultdict(int)
for t in ids: cnt[track_artist[t]] += 1
assert max(cnt.values()) <= 2                                                 # per-artist cap
assert res == generate_candidates("u1", 8, {"artist": 1.0, "co_listen": 1.0, "embedding": 1.0}, heard)
print("checks passed: no duplicates, nothing heard, artist cap, deterministic")

# the cosine used above agrees with scikit-learn
print("cosine matches sklearn:", bool(np.allclose(
    (emb @ user_vec) / (np.linalg.norm(emb, axis=1) * np.linalg.norm(user_vec)),
    cosine_similarity(emb, user_vec[None]).ravel())))

# weights change the mix: count how many of the 8 slots each source proposed
def mix(weights):
    c = defaultdict(int)
    for _, _, p in generate_candidates("u1", 8, weights, heard):
        for s in p: c[s] += 1
    return dict(sorted(c.items()))
eq = [t for t, _, _ in generate_candidates("u1", 8, {"artist": 1, "co_listen": 1, "embedding": 1}, heard)]
x3 = [t for t, _, _ in generate_candidates("u1", 8, {"artist": 1, "co_listen": 3, "embedding": 1}, heard)]
print("equal weights ids  :", eq)
print("co_listen x3 ids   :", x3)
print("equal weights      :", mix({"artist": 1, "co_listen": 1, "embedding": 1}))
print("co_listen x3       :", mix({"artist": 1, "co_listen": 3, "embedding": 1}))
print("embedding only     :", mix({"embedding": 1}))

# robustness: a failing source, an empty source, fewer than n available
def broken(user, limit): raise TimeoutError
SOURCES["co_listen"] = broken
print("co_listen down, still returns:", len(generate_candidates("u1", 8, {"artist": 1, "co_listen": 1, "embedding": 1}, heard)))
SOURCES["co_listen"] = lambda user, limit: []
SOURCES["embedding"] = lambda user, limit: []
print("only artist source, n=8 asked:", len(generate_candidates("u1", 8, {"artist": 1, "co_listen": 1, "embedding": 1}, heard)))
SOURCES["artist"] = lambda user, limit: [30, 31]       # 30 remasters heard track 2, 31 remasters track 5
print("remasters served:", [t for t, _, _ in generate_candidates("u1", 8, {"artist": 1}, heard)])
SOURCES["artist"] = lambda user, limit: [2, 5]
print("user heard remaster 30, original 2 proposed:", [t for t, _, _ in generate_candidates("u1", 8, {"artist": 1}, {30})])
print("n = 0 returns:", generate_candidates("u1", 0, {"artist": 1}, heard))

Output:

track 17 artist_1  score 0.2006 from ['artist', 'co_listen', 'embedding']
track  9 artist_1  score 0.1938 from ['artist', 'co_listen', 'embedding']
track 12 artist_4  score 0.1749 from ['artist', 'co_listen', 'embedding']
track  5 artist_5  score 0.1534 from ['co_listen', 'embedding']
track 20 artist_4  score 0.1422 from ['co_listen', 'embedding']
track 18 artist_2  score 0.1181 from ['artist', 'embedding']
track 14 artist_6  score 0.1059 from ['co_listen', 'embedding']
track 34 artist_2  score 0.0955 from ['artist', 'embedding']
  artist    rank  5: 1/(10+5) = 0.0667
  co_listen rank  6: 1/(10+6) = 0.0625
  embedding rank  4: 1/(10+4) = 0.0714
  track 17 total 0.2006
  artist    rank  3: 1/(10+3) = 0.0769
  co_listen rank  3: 1/(10+3) = 0.0769
  embedding rank 15: 1/(10+15) = 0.0400
  track 9 total 0.1938
rank 1 vs rank 2 term ratio, k=10: 1.091  k=1: 1.5
checks passed: no duplicates, nothing heard, artist cap, deterministic
cosine matches sklearn: True
equal weights ids  : [17, 9, 12, 5, 20, 18, 14, 34]
co_listen x3 ids   : [9, 5, 17, 12, 20, 14, 3, 18]
equal weights      : {'artist': 5, 'co_listen': 6, 'embedding': 8}
co_listen x3       : {'artist': 4, 'co_listen': 7, 'embedding': 7}
embedding only     : {'embedding': 8}
co_listen down, still returns: 8
only artist source, n=8 asked: 5
remasters served: [5]
user heard remaster 30, original 2 proposed: [5]
n = 0 returns: []

Reading the output:

• One score decomposed. Track 17 has the top score, 0.2006. The trace lines show where it comes from: the artist source lists it at rank 5, giving $1/(10+5) = 0.0667$; co-listen at rank 6, $1/(10+6) = 0.0625$; the embedding source at rank 4, $1/(10+4) = 0.0714$. Their sum is 0.2006 (all weights are 1). Track 9 is different: two sources rank it 3rd (0.0769 each) but the embedding source has it only 15th (0.0400), so it totals 0.1938. A rank is the position in the source's own list, and the already-heard tracks 1 and 2 still occupy positions 1 and 2 of the artist list, so the first unheard artist track sits at rank 3.
• Tracks 17, 9 and 12 lead because all three sources propose them. The artist source and the co-listen source list track 12 only as its remaster, track 33, and the canonical lookup merges that into the original, so the remaster neither takes a second slot nor splits the score.
• Within one source, a track and its remaster count once, at the better rank: co-listen lists track 5 at rank 1 and its remaster 31 at rank 2, and only rank 1 is scored. Without that, one source could double-count one song.
• The remasters served line shows the heard filter working on canonical ids: asked for a remaster of heard track 2 (track 30) and a remaster of unheard track 5 (track 31), the function serves only track 5. The filter also works the other way round: the heard set is canonicalised too, so a user who heard the remaster (track 30) is not served the original (track 2); the next line shows only track 5 coming back. A request for n = 0 returns an empty list instead of the whole fused list.
• The artist cap holds: artists 1, 2 and 4 hold two tracks each, and no artist has more.
• Weights change the mix. The numbers count how many of the 8 returned tracks each source proposed; a track proposed by several sources counts for each, so the counts add to more than 8. With equal weights: artist 5, co-listen 6, embedding 8. With co-listen weighted 3: artist 4, co-listen 7, embedding 7. The id lists show why. Tripling the co-listen weight lets track 3, which only co-listen proposes (rank 5, so $3/(10+5) = 0.2$), jump into the list, and it pushes out track 34, which was proposed by the artist and embedding sources and scored only 0.0955. That one swap removes a track from both the artist and embedding counts and adds one to co-listen. The embedding count falls by one for the same reason, not because the embedding source got worse.
• If the co-listen index times out, 8 candidates still return. If only the artist source works, only 5 candidates return for the 8 requested: the cap and the heard filter leave fewer, and the function returns 5 instead of padding with arbitrary tracks. A caller that needs exactly $N$ should add a fallback list (popular in the user's region) and mark it as a fallback.

Choosing the weights

Start from equal weights, then tune on logged outcomes: for each source, how often are its tracks played to completion, and how often does it uniquely contribute a track that gets played? Give a higher weight to sources with high completion, and watch coverage (a source with a high weight can crowd out new-artist discovery). Tune offline with logged propensities, confirm with an online test, and keep the weights in configuration so they can change without a code release. Rank-based fusion is a good first version because it needs no score calibration; if you later want finer control, replace it with a ranker model that takes the source ranks and scores as features.

Complexity

Each source costs its own lookup, bounded by per_source_limit. Fusion visits at most $S \times L$ entries ($S$ sources, $L$ the per-source limit) and sorts the $M$ distinct tracks in $O(M \log M)$. The output is bounded by $N$.

Edge cases and pitfalls

• Same track in several sources is a feature (agreement raises the score), but only if ids are canonical first; otherwise the same song takes several slots.
• Source failure or timeout: isolate each call, give it a time budget, and log failures.
• Fewer than $N$ available: return fewer, or top up from a labelled fallback.
• Cold start (a new user, or a new track, with no history): a new user has no co-listen or follow data; the embedding source then needs a default vector or the weights shift to a popularity source.
• Artist cap interacts with the shortage: it can reduce the count below $N$.
• Position and source bias in logs: log which source proposed each shown track and its rank in that source, so later training can separate source quality from position.

Direct answer

Train one model with a separate output for each objective (a multi-task model with a click head and a value head; a head is the final branch of the network that produces one prediction), keep both outputs calibrated so they mean real probabilities and real money, and combine them at serving time. When the business statement is "maximise clicks, but average revenue per click must stay above a floor", use the constraint form: add a Lagrange multiplier $\mu$ (a price on violating the constraint, explained with numbers below) to the score and tune $\mu$ until the constraint holds. When there is no hard floor, use a weighted sum and pick the weight from the Pareto curve (the set of trade-offs where you cannot improve one objective without losing the other). The trade-off is chosen by a business budget ("how many clicks may we give up for how much revenue"), not by a model.

Training options

Keep the heads separate and combine them at serving: the trade-off is a business setting that will change, and changing a serving-time weight is a config change while changing a training-time weight is a retrain.

Serving: how to combine

Scalarising (turning two objectives into one number) comes in two forms.

Weighted sum. With $p$ the click probability and $v$ the revenue per click, the expected revenue of showing an item is $p \cdot v$, so

Constraint with a multiplier. The goal is to maximise total expected clicks $\sum p_i$ subject to revenue per click at least $R$, that is $\sum p_i (v_i - R) \ge 0$ over the shown items. The algebra: revenue per click over the shown items is $\sum p_i v_i / \sum p_i$; requiring it to be at least $R$ and multiplying both sides by $\sum p_i > 0$ gives $\sum p_i v_i \ge R \sum p_i$, and moving the right side over gives the form above. The quantity $\sum p_i (v_i - R)$ is the slack: positive when the floor is met with room to spare, negative when it is missed. Form the Lagrangian (the objective plus $\mu$ times the constraint slack, with $\mu \ge 0$; $\mu$ must not be negative because a negative price would reward missing the floor):

Items whose value per click is above $R$ get scored up, below $R$ scored down, and $\mu$ sets the strength. Nothing else changes: the ranker still orders by one number.

Two items by hand

Before the large simulation, here is the multiplier score on two items, so you can follow every number. Item X is a niche item (click probability 0.10, revenue per click 2.00) and item Y is a popular one (0.30 and 0.20); the floor is $R = 0.42$. The second half of the program shows why straight-line scores like these can miss a good trade-off: three possible list outcomes are described by (clicks, revenue) and scored with the weighted sum.

import numpy as np

# Two candidate items: p = click probability, v = revenue per click. Floor R on revenue per click.
items = {"item X": (0.10, 2.00), "item Y": (0.30, 0.20)}
R = 0.42
print(f"floor R = {R}")
for mu in (0.0, 0.3, 1.0):
    scores = {k: p * (1 + mu * (v - R)) for k, (p, v) in items.items()}
    order = sorted(scores, key=scores.get, reverse=True)
    detail = "  ".join(f"{k}: {p:.2f}*(1+{mu}*({v:.2f}-{R}))={scores[k]:.4f}" for k, (p, v) in items.items())
    print(f"mu={mu:<3}  {detail}  -> order {order}")
print(f"X and Y swap when mu = {(0.30 - 0.10) / (0.10 * (2.00 - R) - 0.30 * (0.20 - R)):.3f}")

# Why a straight-line (linear) score can miss a good trade-off: three possible list outcomes (clicks, revenue).
outcomes = {"all-clicks list": (10.0, 2.0), "balanced list": (4.0, 4.0), "all-revenue list": (2.0, 10.0)}
print("\nlinear score (1-w)*clicks + w*revenue, best outcome for each w:")
for w in (0.0, 0.2, 0.4, 0.5, 0.6, 0.8, 1.0):
    s = {k: (1 - w) * c + w * r for k, (c, r) in outcomes.items()}
    print(f"  w={w:.1f}  best = {max(s, key=s.get)}")
floor_rev = 4.0
ok = {k: c for k, (c, r) in outcomes.items() if r >= floor_rev}
print(f"most clicks with revenue >= {floor_rev}: {max(ok, key=ok.get)} ({ok[max(ok, key=ok.get)]} clicks)")

floor R = 0.42
mu=0.0  item X: 0.10*(1+0.0*(2.00-0.42))=0.1000  item Y: 0.30*(1+0.0*(0.20-0.42))=0.3000  -> order ['item Y', 'item X']
mu=0.3  item X: 0.10*(1+0.3*(2.00-0.42))=0.1474  item Y: 0.30*(1+0.3*(0.20-0.42))=0.2802  -> order ['item Y', 'item X']
mu=1.0  item X: 0.10*(1+1.0*(2.00-0.42))=0.2580  item Y: 0.30*(1+1.0*(0.20-0.42))=0.2340  -> order ['item X', 'item Y']
X and Y swap when mu = 0.893

linear score (1-w)*clicks + w*revenue, best outcome for each w:
  w=0.0  best = all-clicks list
  w=0.2  best = all-clicks list
  w=0.4  best = all-clicks list
  w=0.5  best = all-clicks list
  w=0.6  best = all-revenue list
  w=0.8  best = all-revenue list
  w=1.0  best = all-revenue list
most clicks with revenue >= 4.0: balanced list (4.0 clicks)

Reading it: with no price ($\mu = 0$) the popular item Y wins, 0.3000 against 0.1000. At $\mu = 0.3$ item X's score rises from 0.1000 to 0.1474 because its revenue per click (2.00) is 1.58 above the floor, while Y's falls from 0.3000 to 0.2802 because 0.20 is 0.22 below it; the order does not change yet. At $\mu = 1.0$ they swap (0.2580 against 0.2340), and the swap point is $\mu = 0.893$. So $\mu$ is the exchange rate: how many clicks you are willing to give up per unit of revenue above the floor. In the second half, the balanced list (4 clicks, 4 revenue) is the best list that meets a revenue floor of 4, yet no weight $w$ ever selects it: below $w = 0.5$ the all-clicks list wins, above it the all-revenue list wins. That is what "linear weights reach only the convex part of the frontier" means (the convex part is the set of outcomes that some straight-line score can pick; an outcome that sits below the line joining two others is a dent that no weight reaches).

Both forms need $p$ to be calibrated (a predicted 0.1 means about 10% of such items are clicked), because $p_i v_i$ is a product and an uncalibrated $p$ distorts it.

Worked example

Synthetic data (a stand-in for real logs): 2000 training requests with 30 candidates each, click probability drawn from Beta(2, 12), and value per click that falls as popularity rises, so the objectives genuinely conflict. A request shows 5 items; clicks are position-weighted by $1/\log_2(\text{rank}+1)$. A second set of 2000 requests is held out (never used to choose the trade-off, so it plays the role of live traffic). Run with Python 3.12, numpy 2.1.3:

import numpy as np

K, N_CAND = 5, 30
POS_W = 1 / np.log2(np.arange(2, K + 2))          # position weights 1/log2(rank+1), ranks 1..K

def make_requests(n_req, seed):
    """Synthetic stand-in data: p = click probability, v = revenue per click, negatively related to p."""
    rng = np.random.default_rng(seed)
    p = rng.beta(2, 12, size=(n_req, N_CAND))
    v = np.clip(1.5 - 5 * p + rng.normal(0, 0.6, size=p.shape), 0.1, None)   # cheap, popular items earn less
    return p, v

def evaluate(p, v, score):
    """Expected position-weighted clicks and revenue per click of the top-K by `score`."""
    idx = np.argsort(-score, axis=1)[:, :K]
    pc = np.take_along_axis(p, idx, 1) * POS_W      # click mass per slot
    rev = (pc * np.take_along_axis(v, idx, 1)).sum()
    clicks = pc.sum()
    return clicks / len(p), rev / clicks

def lagrangian_score(p, v, mu, target):
    # maximise clicks subject to revenue/click >= target  <=>  sum p*(v - target) >= 0
    return p * (1 + mu * (v - target))

def weighted_score(p, v, w):
    return (1 - w) * p + w * p * v

p_tr, v_tr = make_requests(2000, 0)      # used to choose the trade-off
p_te, v_te = make_requests(2000, 1)      # held-out requests

c0, r0 = evaluate(p_tr, v_tr, p_tr)
print(f"click-only ranking (train): clicks/request={c0:.4f}  revenue/click={r0:.4f}")
print("\nweighted-sum sweep, score = (1-w)*p + w*p*v   (train)")
for w in [0.0, 0.25, 0.5, 0.75, 1.0]:
    c, r = evaluate(p_tr, v_tr, weighted_score(p_tr, v_tr, w))
    print(f"  w={w:.2f}  clicks/request={c:.4f}  revenue/click={r:.4f}  clicks kept={c / c0:.3f}")

target = 1.5 * r0                           # business constraint: revenue per click at least 50% above click-only

def solve_mu(target_rpc):
    lo, hi = 0.0, 500.0
    for _ in range(50):                     # bisection on mu: revenue/click rises with mu
        mid = (lo + hi) / 2
        c, r = evaluate(p_tr, v_tr, lagrangian_score(p_tr, v_tr, mid, target))
        lo, hi = (mid, hi) if r < target_rpc else (lo, mid)
    return hi

print(f"\nconstraint: revenue/click >= {target:.4f}")
mus = {}
for label, margin in [("calibrated exactly", 1.00), ("with 3% safety margin", 1.03)]:
    mu = solve_mu(target * margin)
    mus[label] = mu
    print(f"{label}: mu={mu:.3f}")
    for name, (p, v) in {"train": (p_tr, v_tr), "held-out": (p_te, v_te)}.items():
        c, r = evaluate(p, v, lagrangian_score(p, v, mu, target))
        cb, rb = evaluate(p, v, p)
        print(f"  {name:9s} clicks kept={c / cb:.4f}  revenue/click={r:.4f}  meets target: {r >= target}")

# the same two multipliers on 20 more fresh request sets: the held-out result above is one draw
print("\nsame multipliers on 20 other held-out sets of 2000 requests")
for label, mu in mus.items():
    ratios = []
    for s in range(2, 22):
        p, v = make_requests(2000, s)
        ratios.append(evaluate(p, v, lagrangian_score(p, v, mu, target))[1] / target)
    print(f"  {label:24s} misses the floor on {sum(r < 1 for r in ratios)} of 20 sets, revenue/target from {min(ratios):.3f} to {max(ratios):.3f}, std {np.std(ratios):.3f}")

click-only ranking (train): clicks/request=0.9109  revenue/click=0.2828

weighted-sum sweep, score = (1-w)*p + w*p*v   (train)
  w=0.00  clicks/request=0.9109  revenue/click=0.2828  clicks kept=1.000
  w=0.25  clicks/request=0.8898  revenue/click=0.4303  clicks kept=0.977
  w=0.50  clicks/request=0.7702  revenue/click=0.7445  clicks kept=0.845
  w=0.75  clicks/request=0.5945  revenue/click=1.1522  clicks kept=0.653
  w=1.00  clicks/request=0.5304  revenue/click=1.3144  clicks kept=0.582

constraint: revenue/click >= 0.4241
calibrated exactly: mu=0.282
  train     clicks kept=0.9785  revenue/click=0.4242  meets target: True
  held-out  clicks kept=0.9786  revenue/click=0.4195  meets target: False
with 3% safety margin: mu=0.304
  train     clicks kept=0.9749  revenue/click=0.4369  meets target: True
  held-out  clicks kept=0.9743  revenue/click=0.4346  meets target: True

same multipliers on 20 other held-out sets of 2000 requests
  calibrated exactly       misses the floor on 12 of 20 sets, revenue/target from 0.965 to 1.021, std 0.012
  with 3% safety margin    misses the floor on 1 of 20 sets, revenue/target from 0.997 to 1.054, std 0.012

What it shows:

• Click-only ranking earns 0.2828 revenue per click.
• Weighted-sum sweep: a small weight ($w=0.25$) lifts revenue per click to 0.4303 and keeps 97.7% of clicks; at $w=0.5$ it is 0.7445 with 84.5% of clicks; at $w=1$ clicks fall to 58.2% of the baseline. Each row is a point on the frontier, and the best choice is a business decision on the exchange rate.
• Constraint form: with the floor set at 1.5 times the click-only value (0.4241), $\mu = 0.282$ meets the floor on the training requests with 97.85% of clicks kept, but on the held-out requests revenue per click is 0.4195 and misses the floor. With a 3% safety margin in the calibration ($\mu = 0.304$) both sets meet it (0.4369 and 0.4346) and 97.4% of clicks are kept. Lesson: the constraint holds only in expectation over a distribution, and the held-out result is a single random draw, so it needs a margin or a controller that keeps adjusting $\mu$ from live data. The extra block of the program, which re-uses both multipliers on 20 other held-out sets, shows the scatter: the exactly calibrated $\mu$ misses the floor on 12 of those 20 sets (revenue per click between 96.5% and 102.1% of the target, standard deviation about 1.2% of it), and the 3% margin version misses on 1 of 20 (at worst 0.3% below the floor). A margin lowers the chance of a miss without removing it.

How solve_mu works: revenue per click only goes up as $\mu$ goes up, so bisection (repeatedly halving an interval that must contain the answer) brackets $\mu$ between 0 and 500, tries the midpoint, moves the lower end up if revenue per click is still below the target and the upper end down otherwise. After 50 halvings the bracket is far narrower than any difference that matters.

How to pick the margin in practice: size it from the scatter of the live number across several fresh sets (other days, other splits), not from one shortfall. Across the 20 fresh sets above, revenue per click moves with a standard deviation of about 0.012 of the target for either multiplier, so a 3% margin is roughly two and a half standard deviations: enough to make a miss rare (1 of 20 here), not impossible. The single held-out miss earlier (0.0046 on 0.4241, about 1.1%) is one draw from that same scatter. Re-estimate the scatter from live data, and let the feedback controller (a loop that measures the live number and nudges $\mu$ up or down) absorb what the margin does not.

The weighted-sum sweep is a diagnostic for drawing the frontier, not the recommendation: its score adds $p$ (clicks) to $p v$ (money), which mix units, so use it to see the shape of the trade-off and ship the multiplier form when the business states a floor.

All values here are simulated, so the specific percentages say nothing about a real product.

Choosing the trade-off

1. Plot the frontier offline, as in the sweep, and mark the points that satisfy the hard business constraints.
2. Take the point with the best primary objective inside the constraint (here, the largest clicks kept), not the maximum of the sum. There is no universal good value of revenue per click or of clicks kept: the business sets the floor $R$ (for example the revenue per click needed to hit the plan) and you ship the smallest $\mu$ that meets it with margin, then report "clicks kept at the floor" (97.4% in the simulation, a synthetic number).
3. Confirm online with an A/B test of two or three arms, with long-run guardrails (retention, repeat visits), because pushing toward revenue per click can train users to click less.
4. Re-tune: the constraint multiplier drifts as the catalogue and prices change. Run it as a feedback controller (increase $\mu$ when measured revenue per click falls below the floor, decrease when it is comfortably above).

Trade-offs and failure modes

• Linear weights reach only the convex part of the frontier. The multiplier score above is also a linear combination, so if the true frontier has a dent, no setting of $w$ or $\mu$ lands inside it; serving a mix of two settings across traffic reaches points between neighbours, and constrained learning or an explicit per-list optimiser is needed to do better.
• Units. Combine in the same units: clicks and money have different scales, so a weight of 0.5 is not "half and half". Prefer expressing everything as expected money or as a constraint with an explicit floor.
• Aggregate versus per-request constraint. The example constrains the average over all requests. Requiring every list to meet the floor is much stricter and costs more clicks.
• Feedback. Boosting high-value items changes what users see and so changes future click data; keep some unboosted traffic as a holdout to keep estimating the true $p$.
• Position effects: the clicks you observe depend on position, so estimates of $p$ must be corrected for it, or the multiplication by $v$ amplifies the bias.

Common mistakes

• Training one model on a weighted label and then being unable to move the trade-off without a retrain.
• Using uncalibrated scores in a product $p \cdot v$.
• Reading a gain on training requests as a guarantee for live traffic.

The idea

A catalog can hold millions of items, but the model that orders items well is expensive per item because it reads many features. So the system runs two stages. Candidate generation (also called retrieval) uses a cheap model to cut millions of items down to a few hundred or thousand that are plausibly good for this user. Ranking then spends heavy compute on only those candidates and puts them in the final order. The two stages are optimized for different things: retrieval for not losing good items (recall), ranking for getting the order at the top right (precision and ordering quality).

What each stage optimizes and what suits it

The first three rows (goal, items scored, model shape) carry the core idea; the last three rows say how each stage is trained and measured. Terms in the table: gradient-boosted trees are many small decision trees added one after another, each correcting the errors of the ones before; a cross feature is a feature built from both user and item, such as how often this user viewed this seller; in-batch or sampled negatives are items the user did not engage with, taken from the other examples in the same training batch or from a random sample, and used as the wrong answers a positive must beat; learning-to-rank trains a model on the ordering of a list rather than on each item alone; AUC is the probability that a randomly chosen engaged item is scored above a randomly chosen unengaged one; log loss penalises predicted engagement probabilities more heavily the more confidently wrong they are.

The towers can only be cheap because the item vector does not depend on the user. A ranker may use a feature such as "this user viewed this seller three times this week", which cannot be baked into an index. That is why the cheap stage is weaker and the expensive stage cannot be run on everything.

Real systems usually run several retrieval sources in parallel (an embedding source, a "recently viewed similar" source, a popularity source, a source for new items), merge and deduplicate them, and hand the union to the ranker. A final re-ranking step (a last reorder of the ranker's output) may then apply business rules such as diversity.

Choosing the size of the candidate set

Two quantities pull in opposite directions. Recall rises with the number of candidates K but with diminishing returns. Ranker cost rises linearly with K, so the latency budget gives a hard ceiling:

$K_{\max} = \frac{T_{\text{rank}}}{c}$

where $T_{\text{rank}}$ is the time the product allows for the ranking stage and $c$ is the ranker's cost per scored item. The procedure is: take the ranker as the reference, find for a sample of requests which items it would put in its top results if it scored the whole catalog, then measure how many of those the retrieval stage keeps at each K. Pick the smallest K whose recall meets the target, and check it is below $K_{\max}$. For example, if the product allows $T_{\text{rank}} = 30$ ms for ranking and the ranker costs $c = 0.02$ ms per item, then $K_{\max} = 30 / 0.02 = 1500$ candidates; a target met at K = 1000 passes this check and a target that needs K = 2000 does not. If it does not pass, the fix is a better retriever or a lighter ranker, not a larger K.

Worked example

The program builds 20,000 items and 200 users with 32-dimensional vectors whose early dimensions carry more signal. The stand-in "ranker" scores with all 32 dimensions. The stand-in "retriever" keeps only the first 8 dimensions, so it is cheaper but weaker. The program measures how many of the ranker's top 10 items per user fall inside the retriever's top K. The per-item ranker cost of 0.02 ms and the 30 ms ranking allowance are assumed figures for illustration, not measurements. The retriever stands in for any cheap first stage that sees a lossy view of what the ranker sees: in the program the weakness is a truncated view of the same vectors, while in a real system it is that the retriever cannot use features that mix user and item. What carries over is the shape (recall rising with diminishing returns while ranker cost rises linearly) and the procedure of measuring recall against the ranker's own top results, not these particular numbers.

import numpy as np

rng = np.random.default_rng(0)
n_items, n_users, d = 20000, 200, 32
scale = 1.0 / np.sqrt(np.arange(1, d + 1))      # early dimensions carry more signal
items = rng.normal(size=(n_items, d)) * scale
users = rng.normal(size=(n_users, d))

# Stand-in for the heavy ranker: it sees all 32 dimensions of the affinity.
full_scores = users @ items.T                   # shape (200, 20000)
# Stand-in for the cheap retriever: it keeps only the first 8 dimensions.
cheap_scores = users[:, :8] @ items[:, :8].T

top_true = np.argsort(-full_scores, axis=1)[:, :10]
order_cheap = np.argsort(-cheap_scores, axis=1)

print("K     recall of ranker top-10 inside the top-K candidates")
for K in [50, 100, 500, 1000, 2000]:
    cand = order_cheap[:, :K]
    hit = [len(set(top_true[u]) & set(cand[u])) / 10 for u in range(n_users)]
    print(f"{K:<5d} {np.mean(hit):.3f}")

# Latency arithmetic under an ASSUMED per-item ranker cost of 0.02 ms.
cost_ms_per_item = 0.02
for K in [500, 1000, 10000, 1000000]:
    print(f"ranker over {K:>7d} items: {K * cost_ms_per_item:>9.1f} ms")

# The ranking-stage allowance (assumed) gives the ceiling on K.
T_RANK_MS = 30.0
print(f"K_max = {T_RANK_MS:.0f} ms / {cost_ms_per_item} ms per item = {T_RANK_MS / cost_ms_per_item:.0f} items")

# Memory for an in-memory index of 1,000,000 item vectors, 64 dimensions, float32 (4 bytes each).
print(f"index memory: {1_000_000 * 64 * 4 / 1e6:.0f} MB")

Output:

K     recall of ranker top-10 inside the top-K candidates
50    0.489
100   0.613
500   0.852
1000  0.917
2000  0.963
ranker over     500 items:      10.0 ms
ranker over    1000 items:      20.0 ms
ranker over   10000 items:     200.0 ms
ranker over 1000000 items:   20000.0 ms
K_max = 30 ms / 0.02 ms per item = 1500 items
index memory: 256 MB

Reading it: the ceiling is $K_{\max} = 30 / 0.02 = 1500$ items. With K=50 the retriever keeps 0.489 of the ranker's top 10, with K=500 it keeps 0.852, K=1000 gives 0.917 and K=2000 gives 0.963. Doubling from 1000 to 2000 buys 0.046 recall for 20 ms of extra ranker time under the assumed cost, so the curve has flattened. If the target is recall around 0.9, K=1000 is the choice here: it costs 20.0 ms of ranking, below the 30 ms allowance, so the check passes. Scoring all 1,000,000 items with the ranker would cost 20000.0 ms under the same assumption, which is why the two-stage split exists.

Latency budgets, evaluation and a 10M-user, 1M-item startup

• Per-stage latency. Fix the end-to-end budget first, subtract network, feature fetching and re-ranking, and give retrieval and ranking what is left. Retrieval time is dominated by the ANN lookup, which grows far more slowly with catalog size than a full scan does. Ranking time is $K \times c$, as above. Parallel retrieval sources run concurrently, so the slowest one sets the retrieval time.
• Offline evaluation per stage. Retrieval: recall@K against logged positives (items the user later engaged with) and against the ranker's top results, sliced by user segment and item age. Ranking: NDCG, AUC and log loss on logged candidate sets. Evaluate the ranker on candidates drawn from the same retriever that serves, because that is the distribution it will face.
• Online evaluation. The end-to-end A/B test on engagement is the arbiter, with guardrail metrics (numbers that must not get worse, such as latency and diversity). A diagnostic that isolates the bottleneck: score a larger K than production in a shadow run (the larger candidate set is scored in the background and compared, without being shown to users). If quality jumps, retrieval is the bottleneck; if it does not, the ranker is.
• Why a learned ranker beats heuristics. A hand-written score such as popularity times recency uses two or three signals with weights someone guessed. A learned ranker combines hundreds of signals, learns interactions between them (a price that is fine for one segment and a deal-breaker for another) from logged outcomes, and can be retrained as behaviour shifts. Heuristics still earn a place as retrieval sources and as guardrails.
• Startup example. With 1M items, a 64-dimensional float32 vector per item is 256 MB (printed above), small enough to hold an exact or approximate index in memory on one machine. The 10M users do not need per-user storage in the index: compute the user vector at request time from recent activity. Start with one embedding source plus a popularity source feeding a gradient-boosted ranker over about 1000 candidates, and add sources only when the shadow run shows retrieval is the bottleneck. If the catalog were 10,000 items, scoring everything costs 200.0 ms under the assumed cost, which may fit the budget and make the split unnecessary.

Trade-offs and common mistakes

• Training the retriever and the ranker on the same loss and features removes the reason for having two stages.
• Tuning K from the ranker's NDCG alone hides losses upstream. Measure retrieval recall directly.
• A single popularity-heavy source starves new and niche items. Items that never enter the candidate set can never be learned about.
• The ranker trained on shown items only sees a selected sample. Log the whole candidate set with positions so training and evaluation reflect what serving does.
• Cascading more than two stages (a chain of rankers, each lighter than the next, that cut the list step by step) is common when K must be large, at the price of more components that can drift apart.

Running the code

Use Python 3.12 with numpy 2.1.3. Save the listing to a file and run it with python; the seed is fixed inside the program.

Direct answer

Log what the production system chose and with what probability, then estimate a new policy's value by re-weighting the logged outcomes by how much more (or less) likely the new policy was to take the logged action. This is off-policy evaluation (OPE: estimating a policy from data collected by a different policy). The framework has four parts: a logging contract that records propensities, a set of estimators (inverse propensity scoring, self-normalised, direct method, doubly robust), diagnostics that tell you when a number cannot be trusted, and a gate that sends only the best one or two candidates to a live A/B test.

A policy is the rule that decides what to show: given a context it gives a probability for each eligible item. The logging policy (or production policy) is the one that generated the log; the candidate policy is the new one you want to score. Overlap means the logging policy gave non-zero probability to every item the candidate might pick. An estimator is consistent if it converges to the true value as the log grows. Bias is a systematic miss; variance is how much the estimate wobbles between different logs; RMSE (root mean squared error) combines both. The effective sample size is how many equally weighted rows would be as noisy as your weighted ones. A slate is an ordered list of several items shown together. Replay here means checking logged behaviour against logged probabilities; a paired bootstrap by user resamples whole users with replacement and compares candidates on the same resamples.

1. Logging contract (the part that cannot be repaired later)

For every impression, store:

• the context (user features, time, device) as the policy saw it;
• the full candidate set that was eligible, not only the shown items;
• the item(s) shown, their positions, and the propensity: the probability that the production policy chose that item or slate (a list), given the context;
• the outcome (click, watch time, purchase) and when it happened.

Propensities exist only if the production policy is stochastic. A deterministic ranker has propensity 1 for what it showed and 0 for everything else, so no other policy can be evaluated on that log. Keep a small randomised slice (for example a few percent of traffic where items are sampled with known probabilities, or a softmax over scores with a temperature) and log it. This is the cost of the framework: some exploration traffic.

2. Why implicit feedback needs care

A click is not a preference. A user can only click what was shown, and top positions get clicks from being at the top. Missing interactions are not negatives, they are unobserved. The framework therefore treats the logged action and its position as part of the treatment, includes position in the context or models it (a position-based click model: the chance to examine slot $k$ times the chance the item is attractive), and defines the outcome in terms of something the product values (for example a watch of at least 60% or a purchase) and not raw clicks.

3. Estimators

Let the logging policy be $\mu(a \mid x)$, the candidate policy $\pi(a \mid x)$, the logged action $a_i$, reward $r_i$ and importance weight $w_i = \pi(a_i \mid x_i) / \mu(a_i \mid x_i)$.

Why dividing by the propensity works. Suppose production shows item B 10 percent of the time and the candidate would show it 50 percent of the time. B's clicks are under-represented in the log by a factor of 5, so each logged B impression is counted 5 times: weight $0.5 / 0.1 = 5$. Exactly, the expected value of one term is $\sum_a \mu(a)\,\frac{\pi(a)}{\mu(a)}\,q(a) = \sum_a \pi(a)\,q(a)$, the candidate's true value, because $\mu$ cancels, which is why $\mu > 0$ wherever $\pi > 0$ is needed. The first program below traces this with numbers.

Which to learn first. Start with IPS (the idea), SNIPS (its practical form) and DR (the safe default when you have a reward model); DM is the baseline to know and clipping is a tuning knob.

• IPS (inverse propensity scoring) is unbiased when $\mu(a \mid x) > 0$ wherever $\pi(a \mid x) > 0$ (overlap) and the logged propensities are correct, but has high variance when weights are large.
• SNIPS (self-normalised IPS) divides by the sum of weights; slightly biased, much lower variance.
• DM (direct method) fits a reward model $\hat q(x,a)$ and plugs it in; low variance, but its bias is whatever the model gets wrong.
• DR (doubly robust) uses the model as a baseline and corrects it with weighted residuals (residual: click minus what the model predicted for the logged item); it is consistent if either the propensities or the reward model is right. In plain words: if the model is right, the residuals are near 0 and the correction adds nothing; if the model is wrong, the weighted residuals have the average error of the model on the candidate's traffic and subtract it out.

A hand-sized trace (one context, three items, production probabilities 0.5, 0.1, 0.4, candidate probabilities 0.2, 0.5, 0.3, true click rates 0.3, 0.8, 0.5, and a model that is wrong on item 1):

import numpy as np

# One context, three items. mu = production policy, pi = candidate policy, q = true click rate, qhat = a slightly wrong model.
mu = np.array([0.5, 0.1, 0.4])
pi = np.array([0.2, 0.5, 0.3])
q = np.array([0.3, 0.8, 0.5])
qhat = np.array([0.3, 0.6, 0.5])
w = pi / mu
print("importance weights pi/mu:", np.round(w, 3))
print("true value of pi:", round(float((pi * q).sum()), 3), "| expected IPS term E_mu[w r] =", round(float((mu * w * q).sum()), 3), "(identical: mu cancels)")

# three logged rows: (item shown, click)
log = [(1, 1), (0, 0), (2, 1)]
a = np.array([x[0] for x in log]); r = np.array([x[1] for x in log], float)
wi = w[a]
dm = float((pi * qhat).sum())
print("logged items", a, "clicks", r, "weights", np.round(wi, 3))
print("IPS   =", round(float(np.mean(wi * r)), 4), " (5.0 + 0.0 + 0.75) / 3")
print("SNIPS =", round(float((wi * r).sum() / wi.sum()), 4), " 5.75 / 6.15")
print("DM    =", round(dm, 4), " sum of pi x qhat")
terms = dm + wi * (r - qhat[a])
print("DR terms:", np.round(terms, 3), "-> DR =", round(float(terms.mean()), 4))

# the same estimators on a large log: DM keeps the model's error, DR does not
rng = np.random.default_rng(0)
n = 200000
A = rng.choice(3, size=n, p=mu); R = (rng.random(n) < q[A]).astype(float); W = w[A]
print(f"n={n}: IPS {np.mean(W * R):.4f}  SNIPS {(W * R).sum() / W.sum():.4f}  DM {dm:.4f}  DR {np.mean(dm + W * (R - qhat[A])):.4f}  (truth {float((pi * q).sum()):.4f})")
ess = W.sum() ** 2 / (W ** 2).sum()
print(f"effective sample size {ess:.0f} of {n} ({ess / n:.1%})")

Output:

importance weights pi/mu: [0.4  5.   0.75]
true value of pi: 0.61 | expected IPS term E_mu[w r] = 0.61 (identical: mu cancels)
logged items [1 0 2] clicks [1. 0. 1.] weights [5.   0.4  0.75]
IPS   = 1.9167  (5.0 + 0.0 + 0.75) / 3
SNIPS = 0.935  5.75 / 6.15
DM    = 0.51  sum of pi x qhat
DR terms: [2.51  0.39  0.885] -> DR = 1.2617
n=200000: IPS 0.6062  SNIPS 0.6081  DM 0.5100  DR 0.6082  (truth 0.6100)
effective sample size 71241 of 200000 (35.6%)

Reading it: the weights are 0.4, 5 and 0.75. With three logged rows the estimates are noisy (IPS 1.9167, SNIPS 0.935, DR 1.2617; the real value is 0.61), so do not read anything into three rows; the large log shows what the formulas do on average. IPS, SNIPS and DR land near 0.61, while DM stays at 0.51 because its model is wrong and nothing corrects it. This is the same pattern as the simulation below, in the case where every number can be checked by hand.

4. Worked example: do the estimators recover the truth?

A synthetic world with 8 items and 3 context features, where the true click probability is known, so the true value of each candidate policy can be computed. In the code, W fixes how each context feature raises or lowers each item's click chance, q(x) is the resulting true click probability, and V is the logging policy's imperfect belief about W. Each candidate is a softmax over the true click log-odds with a temperature: pi_A uses weight 2, pi_B weight 5 (sharper, nearly always the best item), and pi_C weight -1 (it prefers the worst items on purpose, as a bad candidate that offline evaluation must rank last). A stochastic logging policy (an imperfect softmax) generates 3000 impressions per log. Candidates: pi_A and pi_B favour items with a high true click probability (pi_B more sharply), and pi_C prefers low ones. Each estimator runs on 60 independently generated logs. The reward model in DM and DR is a logistic regression with the item one-hot and context-by-item features; "nocontext" uses a deliberately wrong model that ignores context. This is a single-item simulation, not slates; a stand-in for a real recommender log. Run with Python 3.12, numpy 2.1.3, scikit-learn 1.5.2:

import numpy as np
from sklearn.linear_model import LogisticRegression

A, D = 8, 3                                   # 8 candidate items per request, 3 context features
rng_w = np.random.default_rng(123)
W = rng_w.normal(0, 1.0, size=(D, A))          # fixed world: click logit of item a for context x is x @ W[:, a] - 1.5
V = W + rng_w.normal(0, 1.0, size=(D, A))      # the logging policy's imperfect idea of W

def q(x):                                      # true click probability, shape (n, A)
    return 1 / (1 + np.exp(-(x @ W - 1.5)))

def softmax(z):
    z = z - z.max(1, keepdims=True)
    e = np.exp(z)
    return e / e.sum(1, keepdims=True)

def logging_policy(x):                         # stochastic, so every item has a non-zero propensity
    return softmax(1.0 * (x @ V))

TARGETS = {                                     # candidate policies to be compared offline
    "pi_A": lambda x: softmax(2.0 * np.log(q(x) / (1 - q(x)))),
    "pi_B": lambda x: softmax(5.0 * np.log(q(x) / (1 - q(x)))),
    "pi_C": lambda x: softmax(-1.0 * np.log(q(x) / (1 - q(x)))),
}
xs = np.random.default_rng(7).normal(size=(200000, D))
true_value = {k: float((f(xs) * q(xs)).sum(1).mean()) for k, f in TARGETS.items()}
print("true value (expected click rate):", {k: round(v, 4) for k, v in true_value.items()},
      "| logging policy:", round(float((logging_policy(xs) * q(xs)).sum(1).mean()), 4))

def feats(x, a):
    onehot = np.eye(A)[a]
    return np.hstack([onehot, (x[:, :, None] * onehot[:, None, :]).reshape(len(x), -1)])

def estimate(seed, n=3000, clip=10.0):
    r = np.random.default_rng(1000 + seed)
    x = r.normal(size=(n, D))
    mu = logging_policy(x)
    a = np.array([r.choice(A, p=p) for p in mu])
    click = (r.random(n) < q(x)[np.arange(n), a]).astype(int)
    prop = mu[np.arange(n), a]                                  # logged propensity of the shown item
    model = LogisticRegression(C=10.0, max_iter=500).fit(feats(x, a), click)
    qhat = np.stack([model.predict_proba(feats(x, np.full(n, j)))[:, 1] for j in range(A)], 1)
    # deliberately misspecified reward model: ignores the context, one average click rate per item
    mean_by_item = np.array([click[a == j].mean() if (a == j).any() else click.mean() for j in range(A)])
    qbad = np.tile(mean_by_item, (n, 1))
    out = {}
    for k, f in TARGETS.items():
        pi = f(x); w = pi[np.arange(n), a] / prop
        dm = (pi * qhat).sum(1)
        out[k] = {"IPS": np.mean(w * click),
                  "IPS_clip10": np.mean(np.minimum(w, clip) * click),
                  "SNIPS": np.sum(w * click) / np.sum(w),
                  "DM": dm.mean(),
                  "DR": np.mean(dm + w * (click - qhat[np.arange(n), a])),
                  "DM_nocontext": (pi * qbad).sum(1).mean(),
                  "DR_nocontext": np.mean((pi * qbad).sum(1) + w * (click - qbad[np.arange(n), a])),
                  "ESS_frac": w.sum() ** 2 / (w ** 2).sum() / n}
    return out

SEEDS = 60
runs = [estimate(s) for s in range(SEEDS)]
print(f"\n{SEEDS} logs of 3000 impressions each; target = pi_B")
print(f"{'estimator':13s}{'mean':>8s}{'bias':>9s}{'RMSE':>8s}")
for e in ["IPS", "IPS_clip10", "SNIPS", "DM", "DR", "DM_nocontext", "DR_nocontext"]:
    v = np.array([r["pi_B"][e] for r in runs])
    print(f"{e:13s}{v.mean():8.4f}{v.mean() - true_value['pi_B']:+9.4f}{np.sqrt(((v - true_value['pi_B']) ** 2).mean()):8.4f}")
print(f"effective sample size as a fraction of n (pi_B): {np.mean([r['pi_B']['ESS_frac'] for r in runs]):.3f}")

order = sorted(true_value, key=true_value.get, reverse=True)
print("\ntrue order of the three candidates:", order)
for e in ["IPS", "SNIPS", "DM", "DR"]:
    ok = np.mean([sorted(TARGETS, key=lambda k: -r[k][e]) == order for r in runs])
    print(f"  {e:6s} recovers the true order in {ok:.2f} of the logs")

true value (expected click rate): {'pi_A': 0.6006, 'pi_B': 0.6354, 'pi_C': 0.0783} | logging policy: 0.5038

60 logs of 3000 impressions each; target = pi_B
estimator        mean     bias    RMSE
IPS            0.6370  +0.0016  0.0560
IPS_clip10     0.5807  -0.0548  0.0584
SNIPS          0.6330  -0.0024  0.0269
DM             0.6316  -0.0038  0.0142
DR             0.6365  +0.0011  0.0254
DM_nocontext   0.5596  -0.0758  0.0766
DR_nocontext   0.6344  -0.0011  0.0434
effective sample size as a fraction of n (pi_B): 0.155

true order of the three candidates: ['pi_B', 'pi_A', 'pi_C']
  IPS    recovers the true order in 0.98 of the logs
  SNIPS  recovers the true order in 1.00 of the logs
  DM     recovers the true order in 1.00 of the logs
  DR     recovers the true order in 1.00 of the logs

Reading it:

• IPS is close to unbiased (+0.0016) with the largest spread (RMSE, root mean squared error, 0.0560). The effective sample size, $(\sum w)^2 / \sum w^2$, is only 15.5% of n for pi_B: of the 3000 impressions, about 465 behave like independent samples, which is the signature of a candidate policy that departs a lot from the logging policy. There is no universal cutoff for "too low": 15.5 percent (a loss of about 6.5-fold in data) is workable here, while a candidate at 2 percent of $n$ would be flagged as unmeasurable under step 5 of the comparison procedure below. Pick a threshold before looking at the estimates, for example a few percent of $n$ or fewer than a few hundred effective rows, and report it.
• Clipping weights at 10 (IPS_clip10) removes the large-weight variance on paper but at a cost of bias: -0.0548, and its RMSE (0.0584) is about the same as IPS (0.0560): the bias it adds is about as large as the variance it removes. Clipping trades variance for bias and is a tuning choice. In practice a cap is chosen above most weights so that it trims only the rare huge ones; in this toy pi_B is sharp, so many weights are large and a cap of 10 cuts a lot of real signal. Do not read this as "clipping never helps"; read it as "a cap set too low costs more bias than it saves in noise".
• SNIPS cuts the RMSE from 0.0560 to 0.0269 in this run, at the price of a small bias (-0.0024).
• DM has the lowest RMSE (0.0142) here only because its model family matches how the data was generated, which a real recommender rarely achieves; this is why DR is preferred as the default despite DM's lower variance in this table. With the context-free model (DM_nocontext) the bias is -0.0758 and RMSE 0.0766, the worst of the table.
• DR with the same wrong model (DR_nocontext) has bias -0.0011 and RMSE 0.0434: the importance-weighted correction fixes the model's error, at higher variance than a correct DM. This is the doubly robust property in action.
• Ranking candidates: the true order is pi_B, pi_A, pi_C. SNIPS, DM and DR recover it in every one of the 60 logs, and IPS in 0.98 of them. The gap between pi_A and pi_B is only 0.035 in true value, so the ordering is the easier job than estimating each value precisely. The recovery fractions, and RMSE differences of a few thousandths (IPS against clipped IPS, SNIPS against DR), belong to this particular set of 60 logs and move when the logs are redrawn; the large gaps (DM_nocontext against every other row, the clipping bias of about -0.05) do not. This toy is deliberately friendly: the logging policy covers every item, the propensities are exact, and the world is stationary. Real OPE fails when those fail, for example when the candidate favours items production almost never showed (low overlap), when logged propensities are wrong, or when behaviour has drifted since the log was written.

5. Efficient comparison of many candidate models

Evaluating $M$ candidates does not need $M$ passes over raw logs or $M$ new experiments:

1. Store the log once, with $\mu(a_i \mid x_i)$ attached to each row.
2. For each candidate, compute only $\pi_m(a_i \mid x_i)$ for the logged action (one forward pass of the candidate on the logged contexts, a vectorised batch), giving the weights $w_{mi}$.
3. Evaluate all candidates' SNIPS and DR estimates from these weights in one pass. The reward model $\hat q$ is fitted once and shared.
4. Use paired comparisons: bootstrap by user and compute the difference of two candidates on the same resamples, which cancels the shared noise and gives a tighter interval than two separate intervals.
5. Report the effective sample size and the largest weight per candidate. A candidate whose effective sample size is a few percent of $n$ is flagged as unmeasurable, not ranked.
6. Pass the top one or two candidates to an online A/B test. Offline numbers choose which to test, they do not replace the test.

6. Diagnostics and failure modes

• Overlap violations. If the candidate often picks items the logging policy rarely or never showed, weights explode or estimates are silent about those items. Check the weight distribution and the share of the candidate's mass on items with $\mu$ near zero.
• Wrong propensities. Logged $\mu$ must be the true sampling probability at decision time, including filters and business rules applied later. Validate with a replay: do the logged shown-item frequencies match the logged propensities?
• Slates. For lists, the action space is huge. Use per-position or user-behaviour-model estimators (assumptions about independent positions) and say which assumptions are made.
• Non-stationarity. Logs from last month evaluate last month's world. Use recent logs and check for drift.
• Calibration against reality. Validate the whole framework by checking, on past launches, whether the offline estimate predicted the online result. Report that agreement.

Common mistakes

• Logging only what was shown and not the propensity or the eligible set.
• Treating a click-through in the logged data as the policy value for a new ranker (this is the policy value for the old one).
• Reporting a point estimate without an interval or the effective sample size.