High-Dimensional Geometry

Chapter 1: The Thin Shell

Hold an orange. The peel is maybe three millimetres on a forty-millimetre fruit, so the peel is 7.5 percent of the radius, leaving an inner ball at 92.5 percent. Volume grows like the cube of the radius, so the peel is

— about 21 percent of the volume. Four fifths of an orange is flesh. Now imagine a 768-dimensional orange with the same proportions. Essentially all of it is peel. Not most. All, to seventeen decimal places. This chapter derives that, and then shows why it is the reason ticket two happened.

Step one: volume scales like the d-th power of radius

You do not need a formula for the volume of a d-dimensional ball. You need one fact: if you scale a shape by a factor s in every one of its d directions, its volume multiplies by s^d.

Check it where you can see it. Double a line segment (d = 1) and you get 2× the length. Double a square (d = 2) and you get 4× the area. Double a cube (d = 3) and you get 8× the volume. The pattern is not a coincidence: volume in d dimensions is a d-fold integral, and pulling a factor s out of each of the d integration variables gives s^d. So for a ball of radius R,

That is the entire mathematical content of this section. Everything shocking follows from raising a number slightly less than 1 to a large power.

Step two: how much is in the outer 5 percent?

Define the shell as the region between radius 0.95R and radius R — the outer 5 percent of the radius. The fraction of the ball’s volume inside the shell’s inner boundary is 0.95^d, so the shell holds 1 − 0.95^d. Let us do d = 100 by hand, because the arithmetic is short and the answer is worth feeling.

So in 100 dimensions, 0.59 percent of the volume is in the inner 95 percent of the radius, and 99.41 percent is in the outermost 5 percent. Push to d = 768 and the same computation gives 768 × (−0.051293) = −39.39, and e^−39.39 = 7.8 × 10⁻¹⁸. Eight parts in a hundred million billion. The inside of a 768-dimensional ball is, for every practical purpose, empty.

Step three: the version that applies to your data

Embeddings are not uniform in a ball, so let us redo this for the distribution that actually shows up in practice and in every simulation in this lesson: the isotropic Gaussian, written x ∼ N(0, I_d). That means each of the d coordinates is an independent draw from a standard bell curve with mean 0 and variance 1. “Isotropic” means it looks the same in every direction — no preferred axis.

We want the distribution of the length ‖x‖. Start with the squared length, which is easier because it is a sum:

Its mean. Expectation is linear, so E[‖x‖²] = ∑ E[x_i²]. Each x_i has mean 0 and variance 1, and variance is E[x²] − (E[x])² = E[x²] − 0. So E[x_i²] = 1, and

Its variance. The terms are independent, so variances add: Var[‖x‖²] = ∑ Var[x_i²]. For a standard normal, E[x⁴] = 3 (the fourth moment of a bell curve — a standard fact you can look up or derive by integrating by parts), so

From the square to the length. Here is the step that produces the surprise. If a random quantity S has mean m and a small spread around it, then √S has spread roughly sd(S) / (2√m) — this is the delta method, and it is just the statement that the function √· locally multiplies deviations by its derivative 1/(2√m). Plug in:

Look at what vanished: d. The absolute spread of the length of a Gaussian vector is about 0.707 no matter what the dimension is. The mean grows like √d; the spread stands still. That is concentration in one line.

20,000 samples per row. The sd column is pinned at 0.70 exactly as derived; only the mean moves.

Step four: from the shell to ticket two

Now the payoff. Your query is a vector q. Your corpus is a cloud of vectors y drawn from that same isotropic Gaussian. What does the distribution of distances from q to the corpus look like?

Expand the squared distance, which is the only algebra in this chapter:

Treat q as fixed — it is one specific query — and let y vary over the corpus. Take the three terms one at a time.

‖q‖² is a constant, and by the previous section it is about d. ‖y‖² has mean d and variance 2d. The cross term q · y is a weighted sum of independent standard normals with weights q_i, so it has mean 0 and variance ∑ q_i² = ‖q‖² ≈ d; the −2 multiplier turns that into variance 4d. The two varying terms are uncorrelated (one is odd in y, the other even), so their variances add:

Apply the delta method one more time to get from squared distance to distance:

Again the dimension cancels out of the spread. Every distance from your query to every document in the corpus is drawn from a distribution with mean √(2d) and standard deviation 0.87. The relative spread — the only thing that determines whether ranking is informative — is

For d = 768 that is 0.612/27.71 = 2.2%, and simulation with 10,000 corpus points gives 2.22%. For d = 3,072 it is 1.1%. Every document in your corpus sits at essentially the same distance from every query.

What that does to nearest-neighbor search

The quantity that decides whether “nearest neighbor” is a meaningful idea is the relative contrast: how much farther the worst match is than the best match, as a fraction of the best match.

Here it is, measured directly: one query, a corpus of 10,000 points, isotropic Gaussian, all dimensions.

This is the formal result of Beyer, Goldstein, Ramakrishnan and Shaft (1999): if the ratio of the variance of the distance to the square of its mean tends to zero as d grows, then the ratio d_max/d_min tends to 1, and they titled the paper with the conclusion — “When is nearest neighbor meaningful?”

Why your search still works — the honest caveat

Everything above assumed the coordinates are independent. Real embeddings are nothing like that. A sentence encoder does not produce 768 independent numbers; it produces 768 numbers that lie close to a curved, clustered surface whose intrinsic dimension — the number of directions the data actually varies in — is typically in the tens, not the hundreds.

And concentration depends on the intrinsic dimension, not the storage width. A corpus that lives on a 12-dimensional surface embedded in R⁷⁶⁸ concentrates like d = 12, which is to say barely at all. Add genuine cluster structure — billing documents genuinely near other billing documents — and the distance distribution becomes multimodal, which is exactly the contrast your ranker needs.

So the practical statement is not “high-dimensional search is hopeless.” It is:

Concept → realization: measuring it on your own corpus

python
import numpy as np

# X: (n, d) float32 array of your real embeddings, already L2-normalized or not
def concentration_report(X, n_probe=200, seed=0):
    rng = np.random.default_rng(seed)
    n, d = X.shape
    probes = rng.choice(n, n_probe, replace=False)

    rel, contrast = [], []
    for i in probes:
        # distances from one probe to everything else: shape (n-1,)
        dist = np.linalg.norm(X - X[i], axis=1)
        dist = np.delete(dist, i)
        rel.append(dist.std() / dist.mean())          # relative spread
        contrast.append((dist.max() - dist.min()) / dist.min())

    return {
        "relative_spread": float(np.median(rel)),
        "relative_contrast": float(np.median(contrast)),
        # the i.i.d. prediction for the SAME width, for comparison
        "iid_prediction": float(0.612 / np.sqrt(d)),
    }

Run it. If your measured relative spread is far larger than 0.612/√d, congratulations: your data has real structure and a low intrinsic dimension, and the curse is mostly theoretical for you. If it matches the i.i.d. prediction, your embeddings are behaving like noise and no amount of index tuning will save the ranking — the problem is upstream, in the encoder or the text you fed it.

Chapter 2: Near-Orthogonality

Start with a puzzle that should worry you. A transformer’s residual stream is 768 numbers wide in the small models and 4,096 in the large ones. Linear algebra says a 768-dimensional space has exactly 768 mutually perpendicular directions — that is what a basis is. And yet these models clearly represent far more than 768 distinguishable things. English has hundreds of thousands of words. A sentence encoder must separate millions of topics.

So either the models are wildly overloaded and confused, or perpendicular is the wrong requirement. It is the wrong requirement. Nothing downstream needs directions to be exactly perpendicular. It needs them to be perpendicular enough that reading one does not accidentally read another. And the number of “perpendicular enough” directions in d dimensions is not d. It is astronomically larger, and it grows exponentially in d.

This chapter derives that number. It is the single most important positive fact about high-dimensional space, and it is the reason embeddings work at all.

The distribution of a random cosine, from nothing

Take two vectors u and v drawn independently and uniformly from the unit sphere in d dimensions — the set of all vectors of length exactly 1. Their cosine similarity is just their dot product, since both have length 1:

The derivation takes three steps and no calculus.

Step 1 — freeze u. The uniform distribution on the sphere is rotationally invariant: if you rotate the whole sphere, the distribution is unchanged. So rotate until u points along the first coordinate axis, u = (1, 0, 0, …, 0). Now the dot product collapses to a single term:

The cosine between two random directions is just the first coordinate of one random unit vector. That reframing is the whole trick.

Step 2 — the mean is zero. If v is a valid draw then so is −v, with equal probability. Their first coordinates are v₁ and −v₁. Pair them up and the average is 0. So E[cos θ] = 0. Two random directions are perpendicular on average in every dimension — even in the plane. That was never the interesting part.

Step 3 — the spread is 1/√d. This is the interesting part, and it needs one line. Because v is on the unit sphere,

Take the expectation of both sides. The right side is the constant 1. The left side is a sum of d terms, and by symmetry every coordinate has the same expected square. So d × E[v₁²] = 1, giving

That is the result. The typical size of a random cosine shrinks like one over the square root of the dimension, and it does so for a reason you can state in a sentence: the total squared length is fixed at 1 and it has to be shared among d coordinates, so each coordinate gets a 1/d share.

The exact shape, and a surprise at d = 3

The variance tells you the width but not the shape. The full density of the cosine t between two random directions in d dimensions is

Read off the three regimes, because they explain why your intuition is calibrated so badly.

The d = 3 case is worth pausing on. It is the geometry you grew up in, and in it every cosine is equally likely. There is no notion of “typically perpendicular.” That is the intuition you carry to 768 dimensions, where the density has an exponent of 382.5 and looks like a needle standing on zero.

The typical absolute cosine follows from the standard deviation. For large d the density near 0 is essentially a bell curve with sd 1/√d, and the mean absolute value of a bell curve is √(2/π) times its sd:

Measured on 20,000 random pairs at d = 768: 0.0289. Two decimal places from a formula with no fitted constants.

How many nearly-orthogonal directions fit?

Now the capacity question. Call two unit vectors ε-nearly-orthogonal if |u · v| ≤ ε. How many vectors can you find in d dimensions that are pairwise ε-nearly-orthogonal?

The construction is almost insultingly simple: pick them at random. Then bound the probability that this fails. Two ingredients.

Ingredient 1 — the tail bound. For two independent uniform unit vectors on the sphere in d dimensions, a standard concentration inequality gives

The exponent is where the magic lives: the failure probability for a single pair dies exponentially in dε².

Ingredient 2 — the union bound. If you have N vectors there are N(N−1)/2 < N²/2 pairs, and the probability that any pair is bad is at most the sum of the individual probabilities. Demand that this total be at most δ = 0.01, so we succeed 99 times out of 100:

That is the capacity of high-dimensional space in one formula. Let us evaluate it by hand for the case you care about — d = 768, ε = 0.3, δ = 0.01.

Three point two million directions in 768 dimensions, every pair of them at a cosine of at most 0.3, and the construction is “draw them from a random number generator.” That bound is deliberately loose — using the sharper Gaussian tail instead of the inequality raises it to about 15 million — but even the pessimistic version is four orders of magnitude past 768.

Guaranteed counts from N ≤ √δ e^dε2/4 at δ = 0.01. A dash means the bound falls below 1 (it guarantees nothing) or the number stopped being interesting. Reading a column left to right: doubling d roughly squares the capacity.

Concept → realization: superposition, and the cost of interference

Capacity alone is not the whole story. If a model stores many features in nearly-orthogonal directions — a scheme the interpretability literature calls superposition — then reading one feature back picks up a little of every other feature that happens to be active. Let us compute exactly how much, because the answer tells you the design constraint.

Say the model has feature directions f₁, f₂, …, each a unit vector, pairwise near-orthogonal. On a given input, k of them are active with magnitude 1, so the residual stream holds

To read feature i, the next layer dots with f_i. Split the sum into the term you want and everything else:

The interference is a sum of k − 1 random cosines, each with mean 0 and standard deviation 1/√d. Independent zero-mean terms add in variance, not in magnitude — this is why a random walk of n steps travels √n and not n — so

Now put in numbers for d = 768. The signal is always 1.

Read the last column as a design rule. Superposition works when k « d — when the number of features active at once is far below the width, even though the number of features in the vocabulary can be in the millions. That condition has a name: sparsity. It is not a convenient assumption someone made up to get the math to work; it is the precise thing that has to be true for a d-wide vector to carry more than d things.

What this buys you as a practitioner

Three consequences you can act on today.

1. Linear probes work for a geometric reason. If concepts occupy nearly-orthogonal directions, then a single linear readout w can pick out one concept while being nearly blind to the others, and you can find w with logistic regression on a few hundred labelled examples. The reason it works is the capacity result on this page, not luck.

2. Adding a feature to a frozen embedding by concatenating a random direction is not crazy. If you append a small multiple of a fresh random unit vector to encode metadata, it will be nearly orthogonal to everything already there, so it perturbs existing similarities by O(1/√d). Quantify it before you ship it, but the operation is sound.

3. A mean random-pair cosine that is not near zero is a red flag. Everything above assumed isotropy. If your encoder produces vectors that all share a large common component — if the mean random-pair cosine is 0.6 instead of 0.03 — then your effective capacity is far below the table, because all your vectors are crammed into a narrow cone. That pathology is called anisotropy, and Chapter 8 shows how to measure and remove it.

Chapter 3: The Johnson–Lindenstrauss Lemma

Your index holds one million vectors of 3,072 float32 numbers. That is

which does not fit in the memory of the machine you wanted to run it on, and every similarity computation costs 3,072 multiply-adds. You want to make the vectors shorter. The question is how much shorter you can go before the answers change.

The obvious move is principal component analysis — find the directions in which your data varies most and keep those. It works, and Chapter 4 compares it fairly. But it requires reading all your data, computing a covariance matrix, and redoing the whole thing whenever the data drifts.

In 1984 William Johnson and Joram Lindenstrauss proved something that sounds like it must be wrong: you can multiply every vector by a random matrix — a matrix that has never seen your data, that you can generate from a seed before your first document arrives — and pairwise distances survive almost intact. And the number of dimensions you need does not depend on how many dimensions you started with.

The statement, in plain words first

And now in symbols, with the constant that comes from the standard modern proof (Dasgupta and Gupta, 2003). For any set of n points in R^D and any ε ∈ (0, 1), if

then there is a linear map f : R^D → R^k with, for every pair x, y in the set,

Three things in that statement are strange, and it is worth naming them before proving anything.

D is not in the formula. Not anywhere. Going from 3,072 dimensions and going from one million dimensions cost exactly the same k. The reason is that the proof never looks at the ambient space — it only ever looks at the n(n−1)/2 difference vectors between your points, and there are that many of them no matter how wide the container is.

n enters only through its logarithm. Multiplying your corpus by a thousand adds 4 ln(1000) / (…) to the budget — a fixed additive amount, not a multiplicative one. Corpus size is nearly free.

ε enters as 1/ε². Halving the distortion you will accept costs you four times the dimensions. Accuracy is expensive; scale is cheap. Larsen and Nelson proved in 2017 that this ε⁻² ln n is optimal — nobody will ever do better.

The construction, and why it preserves lengths

Here is the map. Build a k × D matrix R whose every entry is an independent draw from a normal distribution with mean 0 and variance 1/k. Then

Because f is linear, f(x) − f(y) = R(x − y), so preserving distances is the same problem as preserving lengths. Let us take one vector x and show that ‖Rx‖² has the right mean.

One row. Let r be a single row: D independent entries, each N(0, 1/k). The dot product r · x is a weighted sum of independent normals, so it is itself normal with mean 0 and variance

And since it has mean 0, E[(r · x)²] equals that variance: ‖x‖²/k.

All k rows. ‖Rx‖² is the sum of the k squared row-products, and expectation is linear, so

That is why the variance of the entries is 1/k and not 1: it is exactly the normalization that makes the projection length-preserving on average. Nothing was tuned.

How tight is it? The variance, by hand

On average is not enough — we need to know the spread. Write r_j · x = (‖x‖/√k) Z_j where Z_j is a standard normal. Then

We met that sum in Chapter 1: k independent squared normals, mean k and variance 2k. So

One clean number, and it is the most useful formula in this chapter. The typical relative error a random projection makes on a squared distance is √(2/k), independent of everything else.

Let us check it against an actual experiment. Take 300 random points in D = 1,000, project with a fresh Gaussian R, and measure the relative error on all 44,850 pairwise squared distances. The mean absolute error of a bell curve is 0.798 times its standard deviation, so the prediction for the typical error is 0.798 × √(2/k).

The theory nails the typical case. Look at the last column, though: at k = 768 the typical pair is off by 4% but some pair is off by 21%. That gap is the entire reason the JL bound has a logarithm in it. The worst of M samples from a bell curve sits about √(2 ln M) standard deviations out, and with M = 44,850 that is √(2 × 10.71) = 4.63 standard deviations — so 4.63 × 5.1% = 24%, which is what we measured. The lemma is not about the typical pair. It is about the worst one.

The worst-case step, and your actual budget

The concentration inequality for the sum of squared normals gives, for a single vector,

Apply the union bound over all n(n−1)/2 < n²/2 difference vectors and demand the total failure probability be below 1:

That is the lemma, derived. Now run your own numbers: a corpus of one million documents, and you want squared distances to hold to ±10%.

Read that carefully, because it is the opposite of the lesson people usually take from JL. To get a guarantee at ε = 0.1 on a million points you need 11,842 dimensions — nearly four times more than the 3,072 you started with. The lemma would have you expand.

Required k from the Dasgupta–Gupta bound. Notice how flat the columns are: a hundred-thousand-fold increase in corpus size costs less than 3× the dimensions. Notice how steep the rows are: that is the ε⁻².

Translating ε into something you can feel

The lemma is stated on squared distances, and people routinely forget it, which makes them think JL is far weaker than it is. Take the square root of both sides:

So the “ε = 0.5, k = 664” row of the table is not as bad as it looks: it guarantees that no distance in a million-point corpus is ever wrong by more than about 30%, using 664 numbers instead of 3,072. Whether that is good enough depends on the contrast in your data — which is exactly the quantity Chapter 1 taught you to measure.

The gap between the guarantee and the practice

Here is the part nobody tells you when they first show you the lemma. The bound protects the worst pair among half a trillion of them. Almost no application needs that.

If you are building a first-stage retriever that will hand its top 200 candidates to an exact reranker, you do not care that one obscure pair of documents had its distance mangled. You care that the typical distance is roughly right and that true neighbors usually survive. For that, the relevant number is √(2/k), and k = 128 — giving 12.5% typical distortion — is frequently enough. Practitioners routinely project to 64 or 128 dimensions with random matrices and measure recall directly, and the recall is far better than the worst-case bound would suggest.

Concept → realization: the projection you would actually ship

A dense Gaussian R costs D × k multiply-adds per vector and D × k floats to store. Dimitris Achlioptas showed in 2003 that you can replace the Gaussian entries with something far cheaper without weakening the theorem: draw each entry as

Two thirds of the matrix is zero, so you do a third of the arithmetic, and the nonzero entries are ±1, so the inner loop is additions and subtractions instead of multiplications. Li, Hastie and Church pushed this further in 2006 with very sparse random projections at density 1/√D, which for D = 3,072 means about 1.8% of entries are nonzero.

python
import numpy as np

def jl_matrix(D, k, seed=0, sparse=True):
    """Data-oblivious: no corpus needed, reproducible from a seed alone."""
    rng = np.random.default_rng(seed)
    if not sparse:
        return rng.normal(0.0, 1.0 / np.sqrt(k), size=(k, D)).astype(np.float32)
    # Achlioptas: +1 / 0 / -1 with probs 1/6, 2/3, 1/6, scaled by sqrt(3/k)
    u = rng.random((k, D))
    R = np.zeros((k, D), dtype=np.float32)
    R[u < 1/6]  = 1.0
    R[u > 5/6]  = -1.0
    return R * np.sqrt(3.0 / k)

# X: (n, D). One matmul. No fitting, no data pass, no state to version.
R = jl_matrix(D=3072, k=256, seed=42)
Z = X @ R.T                                   # (n, 256) float32

# Verify the promise on YOUR data instead of trusting the bound.
i, j = rng.integers(0, len(X), (2, 20000))
d_hi = ((X[i] - X[j])**2).sum(1)
d_lo = ((Z[i] - Z[j])**2).sum(1)
err  = np.abs(d_lo / d_hi - 1.0)
print("typical", np.median(err), "| p99", np.quantile(err, .99),
      "| theory", np.sqrt(2/256))

Two properties of that code are worth calling out because they are the real reason random projection survives in production. It takes no data pass — you can project the first document you ever receive. And the matrix is fully described by one integer, the seed, so there is no artifact to store, version, or keep in sync between your indexer and your query server.

Chapter 4: Projection in Production

Chapter 3 gave you a theorem. This chapter is about the three things people actually do with it, when each one is right, and the one number that decides between them.

There are exactly three ways to make an embedding shorter, and they differ in when they look at your data.

Why PCA usually wins on real data — and why that is not a contradiction

Principal component analysis finds the orthogonal directions along which your data varies most, ranks them by variance, and keeps the top k. It is optimal in a precise sense (the Eckart–Young theorem): no rank-k linear map has lower expected squared reconstruction error on that dataset.

Random projection has no such property, and it cannot: its guarantee is distribution-free, which means it must hold for the most hostile data you could hand it. The most hostile data is data with no structure at all — a cloud that varies equally in every direction. Real embeddings are the opposite: their variance decays fast.

Here is the comparison with numbers. Take a toy 8-dimensional cloud whose variance along its eight principal directions is

Keep three dimensions.

PCA wins overwhelmingly here, and it wins because the spectrum is steep. Flatten the spectrum — make all eight λ equal to 1 — and PCA keeps 3/8 = 37.5% of the variance, no better than any other choice of three directions, and the gap vanishes. The advantage of PCA over random projection is exactly the amount of structure in your data.

Matryoshka: moving the decision into the loss

Matryoshka representation learning (Kusupati et al., 2022) asks a different question: instead of compressing after the fact, why not train the encoder so that its own prefixes are good embeddings? The recipe is a single change to the objective. Pick a nesting list — say 768, 384, 192, 96, 48, 24, 12 — and train with

where x[:k] is the first k coordinates, renormalized. Gradient from the k = 12 term pushes the most important information into the first twelve slots; the k = 768 term keeps the full vector sharp. What comes out is a vector whose coordinates are ordered by importance by construction.

The operational payoff is that truncation is now a slice, not a matrix multiply, and you can pick k per request. That enables adaptive retrieval: shortlist with a cheap prefix, then rescore the shortlist at full width.

The number that actually gets you promoted: the memory bill

One million vectors. Every row below is the same corpus.

The two axes multiply, and they are independent techniques: dimensionality reduction removes coordinates, quantization shrinks each coordinate. Going from 3,072 float32 to 256 int8 is a 48× reduction — 12.29 GB down to 256 MB, which is the difference between “a dedicated machine” and “a process on your API server.”

Concept → realization: the two-stage funnel, costed

Nobody who has read Chapter 3 carefully uses a projected vector for the final ranking. Here is why, and it is a computation you should be able to do at a whiteboard.

A random projection also perturbs inner products, and by the same algebra as before you can get the size of the perturbation exactly. For unit vectors x and y with true cosine c,

So at k = 256, a cosine near zero picks up noise with standard deviation

Compare that with Chapter 2: the entire spread of random-pair cosines in an isotropic 768-dimensional space is 0.036. The projection noise is nearly twice the whole background distribution. That sounds fatal until you notice what it is being compared against: a genuine neighbor has a cosine of maybe 0.5, which is eight noise units above the background. True neighbors survive. What does not survive is the ordering within the top of the list — rank 3 versus rank 11 is decided by differences far smaller than 0.0625.

Now cost it. One million documents, 3,072-dimensional originals, shortlist of 500.

And the memory footprint of stage 1 is 64 MB at int8, which fits comfortably in a modest container while the full vectors live on disk and are read only for the 500 survivors.

python
# Two-stage retrieval. Stage 1 is cheap and approximate; stage 2 is exact and tiny.
def search(q_full, X_small, X_full, R, shortlist=500, topk=10):
    # q_full: (3072,)  X_small: (n, 64) int8  X_full: memmapped (n, 3072) float32
    q_small = (q_full @ R.T).astype(np.float32)         # project the QUERY the same way
    scores  = X_small.astype(np.float32) @ q_small       # (n,)  -- 6.4e7 MACs
    cand    = np.argpartition(-scores, shortlist)[:shortlist]

    # exact cosine on the survivors only -- 1.5e6 MACs, reads 500 rows from disk
    V = X_full[cand]
    exact = (V @ q_full) / (np.linalg.norm(V, axis=1) * np.linalg.norm(q_full))
    return cand[np.argsort(-exact)[:topk]]

The single most common bug in this pattern is forgetting that the query must go through the identical transformation as the documents. Same R, same seed, same normalization, same order of operations. A random projection is only distance-preserving between vectors that went through the same random projection; project the corpus with seed 42 and the query with seed 43 and you have built a very fast random-number generator.

Where the shortlist quietly fails

Two failure modes are worth knowing before you ship this.

Recall is not uniform across the corpus. Projection noise is the same size for everyone, but the margin between a document and its competitors is not. Documents in dense regions — where hundreds of near-ties sit within 0.05 cosine of each other — get shuffled by projection noise and drop out of the shortlist. Documents in sparse regions survive easily. So compressing the index does not lose recall evenly; it loses it preferentially in your most crowded topics, which are usually the ones that matter most.

The shortlist inherits hubness. Whatever universal neighbors exist in the full space are still there in the projected space, and they consume slots in every shortlist. If a hub occupies one of your 500 slots for every query, that is one fewer real candidate, forever. Chapters 5 and 6 are about that.

Chapter 5: Hubness

Back to ticket three. One generic article is in the top ten for a huge share of queries; a quarter of the corpus is in nobody’s top ten at all. This chapter shows that this is forced by the geometry, identifies exactly which documents become hubs and why, and gives you the three numbers to put on a dashboard.

The measurement: k-occurrence

Define, for every point x in your corpus, its k-occurrence N_k(x): the number of other points that have x in their k nearest neighbors. If x is nobody’s neighbor, N_k(x) = 0. If x is in every list, N_k(x) = n − 1.

Before anything else, notice a conservation law. Each of the n points hands out exactly k votes — one to each of its k nearest neighbors — so the votes must add up:

The mean is pinned at k in every dimension, for every dataset, forever. So hubness is never about the average. It is purely a question of shape: whether those n·k votes are spread evenly or collect on a few winners. Every extra vote a hub receives is a vote taken from somebody else, which is why hubs and orphans always arrive together.

A five-point toy you can do on paper

Put five points in a plane: A at the origin, and B, C, D, E at distance 3 along the four compass directions.

Take k = 1 and find each point’s nearest neighbor by hand.

A is a hub with four times the average, and C, D, E are orphans — and A did nothing special. It is simply the point nearest the centroid. Let us put a number on the skew, since we will use this statistic constantly. The skewness of a list is its average cubed deviation divided by the cube of its standard deviation.

Positive skewness means a long right tail: a few points far above the mean, many bunched below it. That single number is the standard hubness measure, written S_Nk, and a value above about 1 is the conventional threshold for “this dataset has real hubness.”

Why low dimensions physically cannot have hubs

The toy above needed a contrived layout to produce a hub of size 4. That is not an accident. In low dimensions there is a hard combinatorial ceiling on how many points can share a nearest neighbor, and it is one of the prettiest arguments in this lesson.

Suppose two points y and z both have x as their nearest neighbor. Claim: the angle at x between them must be at least 60°.

Proof. Suppose the angle ∠yxz were less than 60°. In any triangle the angles sum to 180°, so the other two angles sum to more than 120°, and therefore at least one of them — say the angle at z — exceeds 60°, which makes it strictly larger than the angle at x. In a triangle, a larger angle faces a longer side. The side facing the angle at z is xy, and the side facing the angle at x is yz. So ‖x − y‖ > ‖y − z‖. But x being y’s nearest neighbor says ‖x − y‖ ≤ ‖y − z‖. Contradiction. □

So all the points that choose x as their nearest neighbor stick out from x in directions that are pairwise at least 60° apart. The maximum number of directions in d dimensions with pairwise angles of at least 60° is a famous quantity: the kissing number τ(d), the number of unit spheres that can touch a central unit sphere without overlapping.

τ(d) grows exponentially — between roughly 2^0.21d and 2^0.40d — so by a few dozen dimensions the geometric ceiling has vanished. This is the reason hubness has an onset. Below about d = 5 the space forbids it; above about d = 20 the space permits it freely. What happens in between is the sharp rise you saw on the dimension dial in Chapter 0.

Who becomes a hub: the centroid channel

Take a query q and a candidate y_i, and let μ be the mean of the data cloud. Expand the squared distance around the mean rather than around the origin:

Three terms with three completely different characters, and the ranking only cares about how they vary across candidates.

The grey term is the same for every candidate. It shifts all scores equally and cannot change the ranking at all — ignore it.

The teal term is the honest signal. It depends on both q and y_i, and it is large exactly when the candidate lies in the direction the query points. This is what you want the ranking to use.

The orange term is the problem. Write ρ_i = ‖y_i − μ‖ for the candidate’s distance from the centroid. This term is ρ_i² — it depends on the candidate and not at all on the query. Every query, forever, applies the same penalty to the same documents. A document that sits slightly nearer the middle of the cloud than average gets a permanent head start in every single competition.

How big is the head start relative to the honest signal? Both are random over candidates, so compare variances. With y_i ∼ N(μ, σ²I) and using the moments from Chapter 1:

A third of everything that decides your search results is a fixed per-document number that has nothing to do with the query. And Chapter 1 already told you why that becomes fatal rather than merely annoying: distances concentrate, so the honest signal produces near-ties by the thousand, and a permanent one-third-of-the-variance tiebreaker wins essentially all of them.

The numbers, measured

1,000 points drawn from an isotropic Gaussian, k = 10, averaged over repeated runs. Remember the mean N₁₀ is exactly 10 in every row.

The skewness and the maximum are extreme-value statistics and jitter run to run — do not read the small wobbles as signal. The orphan share is the stable, monotone measurement, and it is the one to alert on.

Read the d = 300 row as an operations engineer. One document out of a thousand sits in 317 of the 1,000 top-10 lists — it is in a third of all results. Nearly a quarter of the corpus is in zero lists: those documents cannot be retrieved by any query in the corpus, and no error is logged anywhere when that happens.

Why it hurts, in four ways

1. Hubs consume result slots. If a hub is in the top ten for 30% of queries, it occupies a tenth of those result pages for no reason. Precision falls, and it falls on every query at once.

2. Orphans are silently unreachable. This is the worse failure, because nothing surfaces it. A document with N_k = 0 will never appear for any query drawn from the same distribution as your corpus. Your recall metric, computed over queries that happen to have labelled answers, will not notice. Your users experience it as “search never finds that page,” which is indistinguishable from “that page does not exist.”

3. Hubs are plausible. The documents that end up near the centroid of an embedding cloud are the generic ones: short, topically neutral, full of common words — a boilerplate page, a table of contents, an empty template. They look like reasonable results. Nobody files a bug for a plausible result, so hubness survives review in a way that a crash never would.

4. Bad hubs poison k-NN classification. Radovanović and colleagues introduced the term bad hub: a hub whose label disagrees with most of the points that select it. Because a bad hub appears in hundreds of neighbor lists, one mislabelled generic point can flip the predicted label of a large slice of the dataset. In a nearest-neighbor classifier or a retrieval-augmented pipeline that votes over retrieved documents, this is a single point of failure with an unusually large blast radius.

Concept → realization: the three numbers on your dashboard

python
import numpy as np
from scipy.stats import skew

def hubness_report(X, k=10, metric="cosine"):
    """X: (n, d). Returns the three numbers worth alerting on."""
    n = X.shape[0]
    if metric == "cosine":
        Xn = X / np.linalg.norm(X, axis=1, keepdims=True)
        S = Xn @ Xn.T                       # (n, n) similarity; higher is closer
    else:
        S = -((X**2).sum(1)[:, None] - 2*X@X.T + (X**2).sum(1)[None, :])
    np.fill_diagonal(S, -np.inf)      # a point is not its own neighbor

    nbrs = np.argpartition(-S, k, axis=1)[:, :k]        # (n, k) neighbor ids
    Nk = np.bincount(nbrs.ravel(), minlength=n)          # k-occurrence

    assert Nk.sum() == n * k                             # the conservation law
    return {
        "skew":        float(skew(Nk)),                    # > 1 means real hubness
        "orphan_frac": float((Nk == 0).mean()),           # silently unreachable docs
        "top_share":   float(np.sort(Nk)[-10:].sum() / (n*k)),  # slot share of the top 10
        "worst_hub":   int(Nk.max()),
        "hub_ids":     np.argsort(-Nk)[:10].tolist(),      # go read these documents
    }

The last line is the one people skip and should not. Print the ten worst hubs and read them. Nine times out of ten they are not a subtle geometric artefact at all: they are your empty-template page, your “page not found” boilerplate, a document that was truncated to its header, or a chunk that consists entirely of a navigation menu. Deleting or re-chunking those documents fixes more hubness than any algorithm in the next chapter — and unlike the algorithms, it costs nothing at query time.

Chapter 6: CSLS and Mutual Proximity

Chapter 5 identified the disease precisely: a query-independent term, ρ_i², contributes a third of the ranking variance and gives the same documents a permanent head start. The cure follows directly from that diagnosis — estimate each document’s head start and subtract it.

There are three standard ways to do that, and they differ only in what “head start” means. Local mean subtraction gives you CSLS. Local scale division gives you local scaling. Turning distances into probabilities gives you mutual proximity. We will derive the first, work the second and third numerically, and then cost all of them.

First, the free fix: center your data

Before any of the clever methods, do the cheap thing. If your embeddings share a large common component — if every vector has a big projection onto one particular direction — then every cosine is inflated by that shared part and the geometry is squeezed into a narrow cone. Subtracting the corpus mean removes it.

Measured effect, on 800 unit vectors in 100 dimensions with k = 10. The only difference between the rows is how much shared mean component the cloud has.

A shared mean component is manufacturing hubness. Centering costs one pass and one stored vector, it is exactly reversible, and it removes a large share of the problem before you touch anything sophisticated.

CSLS, derived from the diagnosis

Cross-domain Similarity Local Scaling was introduced by Conneau, Lample, Ranzato, Denoyer and Jégou (2018) for unsupervised word translation, where hubness was wrecking the mapping between two language embedding spaces. The construction is short.

For a point y, define its local density baseline: the mean similarity between y and its own K nearest neighbors.

Read that as “how similar is y to things in general.” A hub sitting in the crowded middle has a large r; a point out on the fringe has a small r. Now score a query x against a candidate y not by raw similarity but by how much that similarity exceeds both parties’ baselines:

The factor of 2 is bookkeeping: pull it out and the expression reads

— the similarity minus the average of the two local baselines. It is a z-score with the division left out. Hubs are penalized in proportion to how much of their similarity is explained by simply being in a crowded neighborhood.

A worked example you can check in your head

One query q, four candidates. A is the hub — a generic document in a dense region. D is out on the fringe. Take K = 10 for the baselines.

with r(q) = 0.50 throughout. The hub drops from first to joint second; B, whose 0.70 was earned rather than inherited from a crowded neighborhood, takes the top slot; and D, which was dead last on raw cosine, climbs to joint second because 0.50 is a lot for a document that is normally similar to nothing.

Does it actually work?

800 points, 100 dimensions, isotropic Gaussian, Euclidean neighbors, k = 10, and K = 10 for the baselines.

The busiest point falls from 234 of the 800 lists to 61, and the share of documents that are unreachable by any query falls from nearly a fifth to under two percent. On the strong-cone cosine cloud from earlier in this chapter the effect is similar: skewness 3.34 → 1.62, busiest 134 → 54, orphans 7.9% → 0.4%.

Mutual proximity: turn the distance into a probability

Schnitzer, Flexer, Schedl and Widmer (2012) attack the same asymmetry from a different angle. Instead of subtracting a local mean, ask a probabilistic question:

Estimate it by modelling the distances from each point as approximately normal with its own mean and spread — parameters you compute once per document — and assuming the two events are independent:

where Φ is the standard normal cumulative distribution and μ_x, σ_x are the mean and standard deviation of all distances from x. Work an example. A hub h is close to everything, so its mean distance is small; an ordinary point p has a typical mean distance.

On raw distance h wins: 12.5 beats 12.8. Now compute the mutual proximities.

The order reverses, and the reasoning is legible in the arithmetic. From q’s side both candidates look close (the first bracket is 0.93 versus 0.88 — almost a tie). From their own side the story splits completely: 12.5 is above the hub’s average distance of 12.0, so the hub does not consider q close at all (0.31); while 12.8 is well below p’s average of 14.0, so p does (0.88). The hub is close to q the way it is close to everything — which is to say, not specially.

The family, and what each one costs

Look at the third column, because it is the reason these techniques are practical and not merely elegant. Hubness correction costs four bytes per document. A hundred million documents is 400 MB of floats — noise next to the index itself — and the query-time cost is one addition per candidate, which disappears into the memory traffic you were already paying.

The precomputation is the real expense: to get r(y) you need the K nearest neighbors of every document, which is an all-pairs problem. Run it approximately with the index you already have, offline, once per rebuild. Approximate neighbors are entirely sufficient here — you are estimating a mean over K values, and getting eight of ten right barely moves it.

python
# ---- offline, once per index build ----------------------------------
def fit_csls(index, X, K=10):
    """r[i] = mean similarity of doc i to its own K nearest neighbors."""
    sims, _ = index.search(X, K + 1)      # (n, K+1); col 0 is the doc itself
    return sims[:, 1:].mean(axis=1).astype(np.float32)   # 4 bytes per doc

# ---- at query time ---------------------------------------------------
def search_csls(index, r, q, topk=10, pool=200, K=10):
    sims, ids = index.search(q[None], pool)          # over-retrieve, then rescore
    sims, ids = sims[0], ids[0]
    r_q = sims[:K].mean()                              # the query's own baseline: free
    csls = 2*sims - r[ids] - r_q
    order = np.argsort(-csls)[:topk]
    return ids[order], csls[order]

Two engineering notes hide in that code. First, you must over-retrieve: CSLS can promote a document from rank 150 into the top 10, and it can only do that if rank 150 was in the pool. A pool of 10× to 20× your topk is the usual setting, and it is the difference between CSLS working and CSLS appearing to do nothing. Second, r(q) comes free from the similarities you already fetched — no extra search.

What these corrections cost you

Your thresholds are gone. CSLS scores are not cosines. They can be negative, they do not live in [−1, 1], and any “reject below 0.7” rule you had calibrated is now meaningless. Recalibrate against a held-out set, and store the calibration next to r.

Cold start. A document indexed five seconds ago has no r(y) yet. Until the next offline pass it will be scored with a default — typically the corpus median r — which biases it slightly. For a corpus with steady ingestion this is invisible; for a live feed where most queries hit recent documents, compute r on insertion from the document’s own K nearest neighbors at insert time.

You cannot push it into the index. Approximate nearest-neighbor structures like HNSW are built on a metric. CSLS is not a metric — it is not even symmetric in its effect on ranking — so it must be applied as a rescoring pass over an over-retrieved pool, not as the index’s distance function. That is exactly the two-stage funnel from Chapter 4, doing double duty.

Chapter 7: t-SNE and UMAP Lie

Someone puts a slide on the screen. It shows your embeddings reduced to two dimensions: four beautifully separated blobs, one of them large and diffuse over on the right, the other three tight and close together on the left. The caption reads “the model has learned to separate the four product categories, and categories 1–3 are closely related while category 4 is distinct.”

Almost every claim in that caption is unsupported by the picture. Not “probably wrong” — unsupported, in the sense that the algorithm which produced the picture makes no attempt to preserve the properties being read off it. This chapter shows exactly which properties survive the trip to two dimensions, which are destroyed, and what you should put on the slide instead.

Start with what is impossible

Before criticising any particular algorithm, establish what no algorithm can do. Here is a small theorem with a two-line proof.

Proof. Suppose p₀, p₁, …, p_k are all at distance 1 from one another. Set v_i = p_i − p₀ for i = 1…k. Each has length 1, and for i ≠ j,

So the matrix of inner products of the v’s has 1 on the diagonal and ½ everywhere else. That matrix is ½I + ½J, whose eigenvalues are (1 + k)/2 once and ½ the rest of the time — all strictly positive. A positive-definite Gram matrix means the v_i are linearly independent, so k ≤ m, so at most m + 1 points. □

Now put numbers on it. A plane holds 3 pairwise equidistant points — an equilateral triangle, and you cannot add a fourth. Three-dimensional space holds 4, the vertices of a tetrahedron. And your embedding space of 768 dimensions holds 769 exactly equidistant points, and by Chapter 2 it holds millions that are equidistant to within a few percent.

What t-SNE actually optimizes

t-SNE (van der Maaten and Hinton, 2008) is a neighborhood embedding: it converts distances into probability distributions over neighbors and matches those distributions. Three steps.

Step 1 — neighbor probabilities in the original space. For each point i,

Everything you need to know about the honesty of t-SNE plots is in the subscript on σ_i. Each point gets its own bandwidth, tuned by binary search so that the entropy of its neighbor distribution hits a target you set: the perplexity, which is roughly “how many neighbors each point should feel it has.”

A point in a dense cluster gets a small σ_i. A point in a sparse cluster gets a large one. After this normalization every point in the dataset has the same effective number of neighbors, and therefore the same apparent density. Density and cluster size are destroyed in step one, before any optimization has happened, by design.

Step 2 — neighbor probabilities in the map. In the 2D map, t-SNE uses a Student-t distribution with one degree of freedom:

The heavy tail is a deliberate mismatch. It exists because of the crowding problem, which is Chapter 1 wearing a different hat: the volume available at moderate distance in 2D is tiny compared with the volume available in 768D, so a faithful map would have to pile all the moderately-distant points on top of each other. The heavy tail lets them spread out instead. It is a fudge factor whose entire purpose is to not preserve moderate distances.

Step 3 — match them. Minimize the Kullback–Leibler divergence by gradient descent on the map coordinates:

The asymmetry, worked

KL is not symmetric, and that asymmetry is the behaviour of the algorithm. Compare the cost of the two possible mistakes with actual numbers.

Exactly a hundred to one, and the ratio is just p_A/p_B: the cost of a mismatch is weighted by the true neighbor probability, so pairs that are not really neighbors contribute almost nothing to the loss no matter where the map puts them. Tearing a neighborhood apart is punished hard. Fabricating proximity between strangers is nearly free.

The ledger: what you may and may not read off the plot

Is UMAP better?

UMAP (McInnes, Healy and Melville, 2018) builds a weighted neighbor graph with a locally-adaptive radius, then optimizes a cross-entropy with negative sampling instead of a KL divergence. It is much faster and it has attractive repulsive terms that keep clusters compact. But it is in the same family: it is a neighborhood embedding, it adapts its bandwidth per point, and its loss weights local structure overwhelmingly.

The one difference that genuinely matters turns out not to be the loss at all. UMAP initializes from a spectral embedding of the neighbor graph; classic t-SNE initializes randomly. Kobak and Linderman showed in 2021 that most of UMAP’s reputation for preserving global structure comes from that initialization, and that t-SNE initialized from the first two principal components performs comparably. If you take one practical change from this chapter, make it this one:

Concept → realization: what to report instead

The picture is a navigation tool. The evidence lives in the original space. Three measurements, all cheap.

1. Neighborhood preservation. For each point, take its k nearest neighbors in the original space and in the 2D map and measure the overlap:

Compute it and put the number on the slide. It is routinely well below 1 — often in the 0.2 to 0.5 range at k = 10 for plots that look immaculate. A reader who sees “preservation@10 = 0.31” under the figure understands immediately that two thirds of the true neighborhoods are not in the picture.

2. Trustworthiness and continuity (Venna and Kaski). Trustworthiness penalizes points that the map brings close but which were far in the original space — the fabricated proximities. Continuity penalizes the opposite — the torn neighborhoods. Reporting both tells the reader which kind of lie this particular picture is telling.

3. Run the claim, not the picture. If the claim is “the classes are separable,” fit a logistic regression in the original space and report its accuracy. If the claim is “these two topics are related,” report the mean cosine between the two groups together with the random-pair baseline from Chapter 1. Both take less code than the plot did.

python
from sklearn.manifold import TSNE, trustworthiness
from sklearn.neighbors import NearestNeighbors

# PCA init, not random -- this is the single highest-value flag
Y = TSNE(n_components=2, init="pca", perplexity=30,
         random_state=0).fit_transform(X)

def preservation_at_k(X, Y, k=10):
    nn_hi = NearestNeighbors(n_neighbors=k+1).fit(X).kneighbors(return_distance=False)
    nn_lo = NearestNeighbors(n_neighbors=k+1).fit(Y).kneighbors(return_distance=False)
    hits = [len(set(a[1:]) & set(b[1:])) for a, b in zip(nn_hi, nn_lo)]
    return np.mean(hits) / k

print("preservation@10 =", preservation_at_k(X, Y))
print("trustworthiness  =", trustworthiness(X, Y, n_neighbors=10))

# The stability check that should accompany every published t-SNE figure:
# three perplexities x three seeds. If the story changes, it was never a fact.
for perp in (5, 30, 100):
    for seed in (0, 1, 2):
        ...

Chapter 8: Diagnosing Your Embedding Space

Everything so far has been a mechanism. This chapter is the instrument panel: six measurements that together tell you what kind of space your encoder actually produced, each one computable in a few lines on a sample of a few thousand vectors, and each one connected to a specific fix.

Run all six the day you adopt a new embedding model, and run them again on every corpus you point it at. They take minutes, and they catch problems that no retrieval metric will surface until a user complains.

Diagnostic 1 — the random-pair baseline, and why you must report z-scores

Sample ten thousand random pairs and compute the mean μ₀ and standard deviation σ₀ of their cosine similarity. This is the noise floor of your space: what “no relationship at all” scores. Every similarity you ever report should be expressed relative to it.

Two models, two very different worlds. Model X is well-centered and close to isotropic; Model Y produces a narrow cone.

The “weaker-looking” 0.31 is ten times the evidence of the impressive-looking 0.71. Now do the thing that actually matters — set a cutoff. For a false-positive rate around one in a thousand you need z ≈ 3.1, so

A team that ports the folk-wisdom threshold “accept above 0.8” from Model X to Model Y has just built a system that accepts nearly everything, and it will look like a relevance problem, a chunking problem, or a prompt problem for weeks. Measure the baseline. Publish it next to the model name.

Diagnostic 2 — anisotropy

That μ₀ = 0.62 is itself the diagnosis. In an isotropic space — one that looks the same in every direction — the mean random-pair cosine is 0 by the symmetry argument of Chapter 2. Anything substantially above 0 means all your vectors share a large common component: they are crammed into a narrow cone rather than spread over the sphere.

This is not a rare pathology. Ethayarajh (2019) measured it across contextual encoders and found the effect is severe in the upper layers of language models — in GPT-2’s final layer, two uniformly random words have an average cosine close to 1. Contrastively trained sentence encoders are far better, because the training objective explicitly pushes negatives apart, but many still sit at μ₀ of 0.5 or more.

Anisotropy costs you real resolution. Your usable similarity range is not [−1, 1]; it is the few standard deviations around μ₀. And Chapter 6 showed it manufactures hubness: the strong-cone cloud had 7.9% orphans where the centered one had none.

Diagnostic 3 — the spectrum and the effective rank

Center your vectors, take the covariance matrix, and look at its eigenvalues λ₁ ≥ λ₂ ≥ … ≥ λ_d. You do not need to read 768 numbers; one summary captures what matters, the participation ratio, also called the effective rank:

Read it as “how many directions is this data really using.” Two sanity checks, both by hand.

Perfectly flat spectrum, d = 8, every λ = 1.

A steep spectrum, d = 8.

Eight slots, three directions of real use. Now the number that will surprise you. A common shape for real embedding spectra is a power law, λ_j ∝ 1/j. Take d = 768:

A 768-dimensional embedding with a 1/j spectrum has an effective rank of about 32. You are paying to store, move and multiply 768 numbers to carry roughly 32 directions of information. That single computation is usually the strongest argument for the compression of Chapter 4 — and it is also the reassuring explanation for why your retrieval works despite everything in Chapter 1: your data never lived in 768 dimensions in the first place.

Diagnostic 4 — intrinsic dimension, in three lines of derivation

The participation ratio is a linear notion. If your data lies on a curved surface it will overstate the true number of degrees of freedom. The two-NN estimator (Facco et al., 2017) measures the nonlinear version, and its derivation is short enough to do here.

Consider a point and its two nearest neighbors at distances r₁ < r₂. If the data is locally uniform on a d-dimensional surface, then the number of points within radius r grows like r^d, and a standard property of such a process makes the ratio

follow a Pareto distribution with density f(μ) = d μ^−d−1 for μ ≥ 1 — and crucially, the local density cancels out, so you do not have to know or model it. Write the log-likelihood for n points and differentiate:

One division. Here it is validated on synthetic data with a known answer — a d-dimensional Gaussian rotated into a 768-dimensional ambient space, 2,000 points:

Exact at low dimension, and note the honest failure in the last row: the estimator is biased downward when the intrinsic dimension is large relative to the sample size, because you need roughly exponentially many points to fill a d-dimensional neighborhood. Treat a two-NN estimate above about 20 as “at least 20” rather than as a number, and cross-check it by rerunning on half the data: if the estimate rises with sample size, you are in the biased regime.

Diagnostic 5 — hubness

Chapter 5 gave you the code. The three numbers to alert on:

And the ordering of the fixes, cheapest first: read the top ten hubs and fix the data; center; then CSLS on an over-retrieved pool.

Diagnostic 6 — norms, and the metric mismatch that silently halves your recall

Compute the coefficient of variation of your vector norms, sd(‖x‖) / mean(‖x‖). Then check which metric your index uses and which one your scorer uses, because a mismatch here is one of the most common production bugs in vector search and it produces no error at all.

For normalized vectors the two metrics are equivalent, and the proof is one line:

Euclidean distance is a strictly decreasing function of cosine, so ranking by one is identical to ranking by the other. But drop the normalization and they diverge immediately: two documents that point in exactly the same direction with norms 3 and 12 have cosine 1 and Euclidean distance 9. An index built on L2 will consider them far apart; a scorer using cosine will consider them identical. Your ANN index will then faithfully return the wrong candidates and report excellent internal recall.

The triage table

Concept → realization: the whole panel

python
import numpy as np
from scipy.stats import skew

def embedding_panel(X, k=10, n_pairs=20000, seed=0):
    rng = np.random.default_rng(seed)
    n, d = X.shape
    Xn = X / np.linalg.norm(X, axis=1, keepdims=True)

    # 1 + 2: noise floor and anisotropy
    i, j = rng.integers(0, n, (2, n_pairs))
    c = (Xn[i] * Xn[j]).sum(1)[i != j]
    mu0, sd0 = float(c.mean()), float(c.std())

    # 3: spectrum + participation ratio (effective rank)
    Xc = X - X.mean(0)
    lam = np.linalg.svd(Xc, compute_uv=False)**2 / n
    pr = float(lam.sum()**2 / (lam**2).sum())

    # 4: two-NN intrinsic dimension
    D = np.linalg.norm(X[:, None] - X[None], axis=2)   # sample first if n is large
    np.fill_diagonal(D, np.inf)
    r = np.sort(D, axis=1)[:, :2]
    id2nn = float(n / np.log(r[:, 1] / r[:, 0]).sum())

    # 5: hubness
    S = Xn @ Xn.T; np.fill_diagonal(S, -np.inf)
    Nk = np.bincount(np.argpartition(-S, k, 1)[:, :k].ravel(), minlength=n)

    return {
        "mu0": mu0, "sd0": sd0,
        "cutoff_p001": mu0 + 3.1*sd0,          # use THIS, not 0.8
        "effective_rank": pr, "width": d,
        "intrinsic_dim_2nn": id2nn,
        "hub_skew": float(skew(Nk)),
        "orphan_frac": float((Nk == 0).mean()),
        "norm_cv": float(np.linalg.norm(X,axis=1).std() / np.linalg.norm(X,axis=1).mean()),
    }

Run it on 5,000 sampled vectors — the all-pairs steps are quadratic, and a sample is statistically ample for every statistic here. Store the output next to the index version. When retrieval quality moves and nobody can explain why, the diff between two of these panels usually says which knob turned.

Chapter 9: Connections

Nine chapters ago you had three tickets that looked like bugs. Here is the whole lesson on one card.

What this lesson does not cover, and where it stops being true

Everything here modelled data as isotropic and independent. That is the worst case, and it is the right place to start because it isolates the geometry from the data. Your corpus is not that. It has clusters, a curved low-dimensional support, correlated coordinates, and long tails. All of that helps — it is why retrieval works — and it is also why you must measure rather than assume. Chapter 8 exists for exactly this reason.

Hubness has causes beyond dimensionality. Class imbalance produces hubs in the majority class. Density variation produces hubs in dense regions independently of d. Encoders trained with a contrastive objective on a skewed corpus concentrate probability mass on frequent patterns. The corrections in Chapter 6 address the symptom regardless of the cause, which is convenient, but if you want the cause you have to look at the data.

Approximate indexes add their own geometry. HNSW, IVF and product quantization each impose structure on top of the metric, and each has its own recall profile as a function of dimension, of hubness, and of clustering. The measurements in this lesson describe the space; auditing an ANN index against exact brute-force search on a sample is a separate and equally mandatory exercise.

The bounds are worst-case. The JL budget protects the worst pair among n²/2; the capacity bound protects the worst pair among N²/2. Real systems almost always care about typical behaviour, which is better by a factor of √(2 ln M). Knowing which regime you are in is most of the skill.

Keep exploring

← Vector Embeddings — what these vectors are and where they come from
← Similarity Metrics — cosine, dot, Euclidean, and when each is the right question
← PCA — the spectrum and the effective rank, derived properly
← Embedding Layers — the lookup table where a model’s directions are stored
→ Vector Databases — HNSW, IVF, product quantization: the index on top of this geometry
→ RAG — where hubness and orphans turn into wrong answers
→ Embedding Benchmarks — how models are compared, and what the numbers hide
→ Contrastive Learning — the objective that pushes negatives apart and fights anisotropy
→ Interpretability — superposition and sparse features, the other side of Chapter 2
→ Matryoshka Representation Learning — nesting the compression into the loss
→ SimCSE — alignment and uniformity, which are Chapter 8’s diagnostics as a training objective
→ word2vec — where the idea of a direction meaning something began

References

Concentration and nearest neighbors
K. Beyer, J. Goldstein, R. Ramakrishnan, U. Shaft. “When Is ‘Nearest Neighbor’ Meaningful?” ICDT, 1999.
K. Ball. “An Elementary Introduction to Modern Convex Geometry.” Flavors of Geometry, MSRI, 1997. — the standard source for the sphere concentration bound used in Chapter 2.

Johnson–Lindenstrauss and random projection
W. B. Johnson, J. Lindenstrauss. “Extensions of Lipschitz mappings into a Hilbert space.” Contemporary Mathematics 26:189–206, 1984.
S. Dasgupta, A. Gupta. “An Elementary Proof of a Theorem of Johnson and Lindenstrauss.” Random Structures & Algorithms 22(1):60–65, 2003. — the constant used in Chapter 3.
D. Achlioptas. “Database-friendly random projections: Johnson–Lindenstrauss with binary coins.” Journal of Computer and System Sciences 66(4):671–687, 2003.
P. Li, T. Hastie, K. Church. “Very Sparse Random Projections.” KDD, 2006.
K. G. Larsen, J. Nelson. “Optimality of the Johnson–Lindenstrauss Lemma.” FOCS, 2017. arXiv:1609.02094
A. Kusupati et al. “Matryoshka Representation Learning.” NeurIPS, 2022. arXiv:2205.13147

Hubness
M. Radovanović, A. Nanopoulos, M. Ivanović. “Hubs in Space: Popular Nearest Neighbors in High-Dimensional Data.” JMLR 11:2487–2531, 2010. jmlr.org
D. Schnitzer, A. Flexer, M. Schedl, G. Widmer. “Local and Global Scaling Reduce Hubs in Space.” JMLR 13:2871–2902, 2012. jmlr.org
A. Conneau, G. Lample, M. Ranzato, L. Denoyer, H. Jégou. “Word Translation Without Parallel Data.” ICLR, 2018. arXiv:1710.04087 — CSLS.

Neighborhood embeddings
L. van der Maaten, G. Hinton. “Visualizing Data using t-SNE.” JMLR 9:2579–2605, 2008. jmlr.org
L. McInnes, J. Healy, J. Melville. “UMAP: Uniform Manifold Approximation and Projection.” 2018. arXiv:1802.03426
M. Wattenberg, F. Viégas, I. Johnson. “How to Use t-SNE Effectively.” Distill, 2016. distill.pub
D. Kobak, G. C. Linderman. “Initialization is critical for preserving global data structure in both t-SNE and UMAP.” Nature Biotechnology 39:156–157, 2021. arXiv:1905.05879
T. Chari, L. Pachter. “The specious art of single-cell genomics.” PLoS Computational Biology 19(8), 2023.
J. Venna, S. Kaski. “Local multidimensional scaling.” Neural Networks 19(6):889–899, 2006. — trustworthiness and continuity.

Anisotropy, intrinsic dimension, superposition
K. Ethayarajh. “How Contextual are Contextualized Word Representations?” EMNLP, 2019. arXiv:1909.00512
J. Mu, P. Viswanath. “All-but-the-Top: Simple and Effective Postprocessing for Word Representations.” ICLR, 2018. arXiv:1702.01417
E. Facco, M. d’Errico, A. Rodriguez, A. Laio. “Estimating the intrinsic dimension of datasets by a minimal neighborhood information.” Scientific Reports 7:12140, 2017. — the two-NN estimator.
N. Elhage et al. “Toy Models of Superposition.” Transformer Circuits Thread, 2022. transformer-circuits.pub