Approximate Nearest Neighbour Search

A course of study, leading to the construction of a vector search index in C++

Chapter 0. How to use this document

This document teaches a subject and then supervises a project. The subject is approximate nearest neighbour search, which is the technique underlying the retrieval component of most modern artificial intelligence systems. The project is a working implementation of a particular algorithm for it, written in C++, without reference to any existing implementation.

Nothing specific to this field is assumed. Every term is defined at the point where it is first needed, and every claim of the form "this method is too slow" is supported by arithmetic that you can check yourself. What is assumed is ordinary programming competence in C++ and a working familiarity with elementary data structures such as arrays, graphs and priority queues. Where a structure is used in an unfamiliar way, that use is explained.

The document is arranged so that the theory precedes the work. Chapters 1 through 8 establish why the algorithm has the form it has. Chapters 9 through 13 direct the implementation. The temptation to begin coding at Chapter 10 should be resisted, because the implementation stages are written on the assumption that the reader understands what each component is for, and a component built without that understanding cannot be debugged when it misbehaves.

A remark on method. Each implementation stage is accompanied by a verification procedure. These are not suggestions. The characteristic way this project fails is that a defect is introduced early, is not detected because nothing tested for it, and is then diagnosed several stages later in the wrong component. A day spent confirming that stage three is correct is cheaper than a week spent searching stage six for a fault that is not there.

The reference paper for the algorithm is Y. Malkov and D. Yashunin, "Efficient and robust approximate nearest neighbour search using Hierarchical Navigable Small World graphs" (arXiv:1603.09320). It should be read after Chapter 7 of this document and not before. Read in that order, it will appear as a compact restatement of ideas already understood. Read first, it is difficult.

Chapter 1. The problem, developed from nothing

1.1 Points and coordinates

A list of numbers can be understood as the coordinates of a point. The pair of numbers (3, 4) describes a position on a plane: three units along the horizontal axis, four units along the vertical. The triple (3, 4, 2) describes a position in ordinary space, adding a height of two.

There is no reason to stop at three. A list of 128 numbers describes a point in what is called a 128-dimensional space. This cannot be visualised, and no attempt should be made to visualise it. It is nevertheless entirely well defined, because everything that will be done with such points is arithmetic on their coordinates, and arithmetic does not care how many coordinates there are.

Throughout this document, the number of coordinates is written D, and a point with D coordinates is called a vector of dimension D. For this project D is 128.

1.2 The distance between two points

On a plane, the distance between two points is obtained by the theorem of Pythagoras. Take the points x = (1, 2) and y = (4, 6). The difference along the first coordinate is 4 − 1 = 3, and along the second it is 6 − 2 = 4. The distance is the square root of the sum of the squares of these differences: the square root of 9 + 16, which is the square root of 25, which is 5.

In three dimensions the rule is identical with one more term. In D dimensions it is identical with D terms. To compute the distance between two vectors, subtract them coordinate by coordinate, square each of the D differences, add the squares together, and take the square root of the total.

This formula is the only piece of geometry required for the entire project. It will be evaluated several billion times before the project is finished, and consequently it deserves more attention than its simplicity suggests.

1.3 Why the square root is discarded

Consider two candidate points, one at squared-sum 25 from the query and another at squared-sum 36. Taking square roots gives distances of 5 and 6. Whether the square roots are taken or not, the first point is the nearer one.

This holds in general. The square root is a strictly increasing function on non-negative numbers, meaning that if a is less than b then the square root of a is less than the square root of b. Therefore taking square roots never changes which of two points is nearer. Since every use of distance in this system is a comparison, and no part of it requires the true geometric distance as a number, the square root can be omitted entirely.

The quantity actually computed, the sum of squared coordinate differences without the root, is called the squared Euclidean distance. From this point onward the word "distance" means this quantity.

The saving is not decorative. The square root is among the slower arithmetic operations available, and removing it from the innermost loop of a system that performs billions of distance evaluations is a measurable improvement obtained at no cost in correctness.

1.4 Embeddings: where the vectors come from

So far the vectors are abstract. Their origin is what makes the problem worth solving.

A neural network can be trained to accept an object of some kind, such as a paragraph of text or a photograph, and to emit a fixed-length list of numbers. Such a list is called an embedding of that object. The training procedure is arranged so that objects which resemble one another in meaning are assigned lists of numbers that are close together under the distance of Section 1.2, while objects unrelated in meaning are assigned lists that are far apart.

Consider two questions submitted to a support system: "how do I reset my password" and "I have forgotten my login details". These share almost no words. A system matching on words would not connect them. An embedding model, however, assigns them nearby positions, because it has been trained on the relationship between language and meaning rather than on the surface form of the words. A sentence about cellular biology would be placed somewhere else entirely.

The importance of this is that it converts a hard problem into an easier one. Judging whether two passages mean similar things is difficult. Measuring whether two points are close together is arithmetic. The embedding model performs the difficult conversion once, and thereafter the question "what do I have that resembles this" becomes the question "which stored points lie nearest to this point". That second question involves no machine learning whatsoever. It is a problem in data structures, and it is the subject of this course.

This explains where such systems appear. When a language model answers a question using a company's internal documents, some component first selected the relevant documents, and that component almost certainly performed a nearest neighbour search over embeddings. The same mechanism underlies semantic search, recommendation systems, detection of near-duplicate documents in training corpora, and the retrieval stage of nearly every retrieval-augmented generation system in production. In such systems the neural network typically runs on specialised hardware, while the search component runs on ordinary processors, is required to sustain thousands of queries per second over collections of millions to billions of points, and is usually written in C++.

1.5 Statement of the problem

A collection S contains N points, each of dimension D. A query point q is supplied, which need not belong to S. An integer k is given. The task is to report the k members of S whose distance to q is smallest.

For this project, N is one million, D is 128 and k is 10. These values correspond to a public benchmark dataset, which permits the finished system to be compared directly against published results for established software.

1.6 A worked example

The following six points lie on a plane, so that the situation can be followed concretely.

A = (1, 1)    B = (2, 3)    C = (8, 8)
D = (9, 7)    E = (1, 4)    F = (7, 9)

Let the query be q = (2, 2), and let k be 2.

The squared distance from q to A is (1 − 2)² + (1 − 2)² = 1 + 1 = 2. To B it is (2 − 2)² + (3 − 2)² = 0 + 1 = 1. To C it is 36 + 36 = 72. To D it is 49 + 25 = 74. To E it is 1 + 4 = 5. To F it is 25 + 49 = 74.

The two nearest points are therefore B, at distance 1, and A, at distance 2.

This is the entire problem. Everything that follows exists because performing this calculation for every one of a million points, ten thousand times a second, is not possible by direct means.

Chapter 2. Why the straightforward solutions fail

2.1 Exhaustive search and its exact cost

The method of Section 1.6, applied to the full collection, is called exhaustive search or linear scan. Compute the distance from the query to all N points and retain the k smallest.

It is correct, it requires perhaps fifteen lines of code, and it is the first component this project will construct, because it provides the reference answers against which everything else is checked. It is nevertheless unusable as a production method, and the arithmetic showing why is worth performing.

One distance evaluation between vectors of dimension 128 requires 128 subtractions, 128 multiplications and 127 additions, which is approximately 384 arithmetic operations. One query against one million points therefore requires approximately 3.84 × 10^8 arithmetic operations.

The memory traffic is more restrictive still. One million vectors of 128 coordinates, each coordinate occupying four bytes, is 512 million bytes. Every query must read all of it. Contemporary processors sustain memory bandwidth in the region of tens of gigabytes per second, which places a floor on the time for a single query in the region of ten to twenty milliseconds even if the arithmetic were free.

A service handling ten thousand queries per second would need to read five terabytes per second. The method is short of the requirement by roughly three orders of magnitude, and the shortfall grows in direct proportion to the size of the collection.

2.2 Why sorting does not help

In one dimension the problem is easy. Sort the numbers once, then locate the query by binary search and inspect its neighbours in the sorted order. The cost per query falls from N to the logarithm of N.

This approach does not extend beyond one dimension, and the reason is worth stating explicitly because the instinct to sort is strong. Sorting requires a total order, meaning an arrangement in which every element has a definite position relative to every other. Points on a plane admit no such arrangement consistent with distance. Sorting by the first coordinate places (0, 100) next to (0, 0) although they are far apart. Sorting by the second coordinate has the symmetric defect. No single ordering can preserve proximity in more than one dimension, because proximity in two dimensions is not a one-dimensional property.

2.3 Dividing the space: grids and trees

The next natural approach is to divide the space into regions. Overlay a grid on the plane and record which points fall in which cell. To answer a query, examine the cell containing the query and the cells immediately surrounding it. If the nearest point found is closer than the distance to the edge of the examined region, that point is the true nearest neighbour and the remaining cells can be ignored.

The tree-based refinements of this idea, of which the k-d tree is the best known, divide the space recursively rather than uniformly, but the principle is identical: partition the space, descend to the relevant region, and exclude the remaining regions without examining them. The efficiency of the whole family depends entirely on that exclusion. Regions must be ruled out without being looked at.

In two or three dimensions this works extremely well. In 128 dimensions it fails completely, and the failure has two independent causes.

2.4 First cause: the number of neighbouring regions

Consider the grid method. In two dimensions, a cell has eight neighbouring cells, so a query examines nine cells in total. The general formula is that a cell in D dimensions has 3^D − 1 neighbours, since along each of the D axes a neighbouring cell may be one step below, level, or one step above.

In two dimensions this is 3² − 1 = 8, which is manageable. In three dimensions it is 26, still manageable. In 128 dimensions it is 3^128 − 1, which is approximately 1.2 × 10^61.

This number exceeds the number of atoms in the observable universe by many orders of magnitude. There is no sense in which a query can examine the neighbouring cells. The method does not become slow; it ceases to be defined.

The same objection applies to trees. Each level of a k-d tree constrains one coordinate. To narrow a query to a small region of a 128-dimensional space, a substantial fraction of the 128 coordinates must be constrained, so the tree must be deep, and the number of leaf regions grows as a power of two in that depth. For the exclusion rule to eliminate most regions, the number of points N would have to greatly exceed 2^128. No dataset of that size can exist.

2.5 Second cause: distances stop distinguishing

The second cause is subtler and more damaging, because it undermines the exclusion rule itself rather than the bookkeeping around it.

Here is an illustration that can be checked by arithmetic. Take the unit cube in D dimensions, that is, the set of points all of whose coordinates lie between 0 and 1, and consider the shell consisting of points within 0.1 of the surface. The interior that remains, excluding the shell, is a smaller cube whose side is 0.8, so its volume is 0.8^D.

For D = 2, this interior volume is 0.64, so 36 per cent of the area lies in the shell. For D = 3, it is 0.512, so about 49 per cent lies in the shell. For D = 128, it is 0.8^128, which is approximately 4 × 10^-13.

That is to say, in 128 dimensions, essentially every point of the cube lies near its surface. This is characteristic of high-dimensional spaces: the volume is concentrated at the boundary, and intuitions formed in two or three dimensions are not merely imprecise but systematically wrong.

A related and more directly relevant effect is distance concentration. When points are distributed through a high-dimensional space, the distances from a query to all of them become increasingly similar as D grows. The nearest point is only slightly nearer than a typical point.

This is fatal for space partitioning. The rule that makes such methods fast has the form: "the closest possible point within this unexamined region is further away than the best candidate already found, therefore skip the region." When all distances are nearly equal, this condition is almost never satisfied, so almost no region is skipped. The search examines nearly everything, having added the overhead of tree traversal to the cost of exhaustive search. It is slower than the method it was intended to improve upon.

2.6 Conclusion

These observations are known collectively as the curse of dimensionality. Their joint consequence is that for high-dimensional data there is no known exact method substantially faster than examining everything, and there are strong reasons to believe none exists.

The field therefore abandoned exactness. That decision, and what is obtained in exchange, is the subject of the next chapter.

Chapter 3. Approximation, and how it is measured

3.1 The exchange

An approximate nearest neighbour method does not promise to return the true k nearest points. It returns k points which are usually the true ones, or very nearly. In exchange it can be faster than exhaustive search by a factor of a thousand or more.

Stated baldly this sounds like an evasion of the problem. It is defensible, and the reason lies in where the vectors came from. The embeddings are produced by a trained model whose notion of similarity is itself imperfect. The system is already approximate before the search begins. In a document retrieval application, the practical difference between returning the true tenth-nearest document and the thirteenth-nearest is generally invisible to the user, whereas the difference between thirty queries per second and thirty thousand determines whether the system can be built at all.

The guarantee being surrendered was never worth very much. The speed obtained in exchange is worth a great deal.

3.2 Recall

The accuracy of an approximate method is measured by a quantity called recall.

The measurement requires knowing the correct answers. These are computed once, in advance and without time pressure, by exhaustive search, and stored. This stored set is called the ground truth. For a benchmark dataset it is supplied along with the data.

For a single query, let T be the set of true k nearest neighbours and let A be the set of k results the method actually returned. The recall for that query is the number of points appearing in both sets, divided by k.

As an example, suppose k is 10 and the true ten nearest neighbours are the points numbered

17, 42, 103, 256, 301, 477, 512, 688, 701, 894

and the method returns

17, 42, 103, 256, 301, 477, 512, 688, 701, 905

Nine of the ten returned points are correct. The recall for this query is 9 / 10 = 0.9. The final figure, written recall@10, is this value averaged across every query in the benchmark, commonly ten thousand of them.

Systems in production typically operate between recall 0.90 and 0.99.

3.3 Why a single number is never adequate

Every approximate method has at least one parameter controlling how hard it searches. This has an important consequence for how results must be reported.

Suppose a method is described as achieving twenty thousand queries per second. This tells you nothing, because that rate could have been obtained by searching carelessly and returning poor answers. Suppose instead a method is described as achieving recall 0.99. This also tells you nothing, because that accuracy could have been obtained by exhaustive search.

Neither number is meaningful alone. They are meaningful only in combination, and since the controlling parameter can be set to any value, the honest description of a method is not a pair of numbers but a curve.

The curve is produced by varying the search-effort parameter across its useful range and recording, for each setting, both the recall achieved and the query rate sustained. Plotting rate against recall traces out the method's behaviour across its whole operating range. One method is superior to another when its curve lies above and to the right, meaning it achieves higher throughput at every level of accuracy.

3.4 The standard benchmark

This methodology is used by a public benchmark suite published at ann-benchmarks.com, which evaluates the established libraries on common datasets and publishes their curves. Producing such a curve for your own implementation, on the same dataset, is the final deliverable of this project. It permits your work to be placed directly alongside published results, which is a considerably stronger claim than a demonstration that the program runs.

Chapter 4. Searching by walking: the graph approach

4.1 The change of viewpoint

Space partitioning failed because it tried to reason about regions of a space whose geometry defeats such reasoning. The successful approach abandons the space and reasons about the points themselves.

Construct a graph. Each point of the collection becomes a vertex. Each vertex is joined by edges to a small number of other vertices which lie near it. A structure of this kind is called a proximity graph. Nothing is stored about the space between the points; only the points and the connections among them.

A query is then answered by walking. Begin at some vertex and move along edges, at each step choosing a vertex that lies nearer to the query than the current one, until no further improvement is available.

The significance of this reformulation lies in what determines the cost. At each step the algorithm computes the distance from the query to each neighbour of the current vertex. The number of such computations depends on how many edges a vertex has, which is a design choice, and not on how many points the collection contains. The total cost of a query is therefore the number of edges per vertex multiplied by the number of steps taken. The design problem becomes that of building a graph in which few steps are needed, and this is a problem about which something can be done.

4.2 The greedy walk

The procedure described informally above is called greedy routing, and it is stated precisely as follows.

Maintain a current vertex, initially some fixed starting vertex called the entry point. Compute the distance from the query to every neighbour of the current vertex. If the nearest of these neighbours is closer to the query than the current vertex is, move to that neighbour and repeat. Otherwise stop and report the current vertex.

A useful way to picture this is to assign to each vertex a height equal to its distance from the query. The walk then proceeds downhill, and stops at the bottom of a valley.

4.3 Where the greedy walk fails

The walk stops at the bottom of a valley, but not necessarily the deepest valley. A vertex may be nearer to the query than all of its own neighbours while still being far from the true nearest point, simply because it has no edge pointing in the direction that would improve matters.

The following arrangement illustrates the failure. Suppose the graph contains a vertex P whose four neighbours all lie further from the query than P does, and suppose that the true nearest neighbour is a vertex R which is not connected to P and is reachable only by first moving to one of P's neighbours, which is a step in the wrong direction. Greedy routing refuses to take that step, because the rule requires every step to be an improvement. It stops at P and reports it.

This situation is not rare. In a graph of a million vertices where each vertex has perhaps sixteen edges, there are many such local valleys, and greedy routing alone therefore produces poor recall. The method must be modified so that it can pass through a temporarily worse vertex in order to reach a better one.

4.4 Beam search: exploring several possibilities at once

The remedy is to stop committing to a single path, and instead to keep track of a number of promising vertices simultaneously, expanding them in order of promise. This procedure is the central component of the entire system, and the remainder of the algorithm is machinery surrounding it.

Two collections are maintained.

The first is the candidate set. It holds vertices that have been discovered but whose neighbours have not yet been examined. It must support extracting the vertex nearest to the query, which is the ordinary function of a priority queue arranged so that the smallest element is extracted first.

The second is the result set. It holds the best vertices found so far, and its size is capped at a fixed limit. It must support extracting the vertex furthest from the query, so that when the cap is exceeded the worst member can be discarded. This is a priority queue arranged in the opposite direction, so that the largest element is extracted first.

The procedure is then as follows. Extract the nearest vertex from the candidate set and examine each of its neighbours. For any neighbour not previously visited, mark it visited, compute its distance to the query, and if it is nearer than the furthest member of the result set, or if the result set is not yet full, add it to both the candidate set and the result set. Whenever the result set exceeds its cap, discard its furthest member. Repeat.

The visited marking matters. Without it the procedure would revisit vertices endlessly, since the graph contains cycles.

4.5 When to stop, and why that rule is correct

The procedure needs a termination condition, and the correct one deserves explanation because it is what separates this from an exhaustive traversal of the graph.

Stop when the nearest remaining candidate is further from the query than the furthest member of a full result set.

The justification is as follows. Candidates are extracted in increasing order of distance from the query. If the nearest remaining candidate is already worse than the worst accepted result, then every remaining candidate is also worse, since they are all further away still. Continuing would mean expanding vertices which are themselves not good enough to be kept, in the hope that their neighbours might be better. The algorithm gives up that hope, and this is precisely the approximation being made: it is possible that a good vertex lay just beyond a bad one and was missed. This is the source of the lost recall, and increasing the cap on the result set is what reduces the loss.

4.6 The control knob

The cap on the result set is written ef. It is the search-effort parameter promised in Section 3.3.

Setting ef to 1 reduces the procedure exactly to greedy routing, since only one result is ever kept and the search stops as soon as no improvement is found. This is fast and inaccurate. Increasing ef causes more alternatives to be retained, allows the search to escape more local valleys, raises recall and costs more time. Sweeping ef across a range and recording the resulting pairs of recall and throughput is what produces the curve described in Section 3.3.

A warning that will be repeated in Chapter 12. The two collections are ordered in opposite directions: the candidate set gives up its smallest element, the result set its largest. If one of them is built the wrong way round, the program will compile, run, terminate normally, and return results that are merely somewhat worse than they should be. Nothing will crash. There will be no obvious place to begin looking. Take care with this.

Chapter 5. Making the walk short

5.1 The difficulty with uniformly short edges

The cost of a query is the number of edges per vertex multiplied by the number of steps. Chapter 4 addressed the quality of the search. This chapter addresses the number of steps.

Suppose every edge of the proximity graph joins two genuinely close points. Each step of the walk then advances a short distance. If the walk begins at a vertex far from the query, it must cross the entire collection in short steps, and the number of steps is large. The search is accurate but slow.

Suppose instead that the graph also contains a small number of long edges joining distant regions. The walk can then use a few long edges to arrive in approximately the right region, and short edges to refine its position once there.

5.2 Small-world graphs

A graph possessing both properties, namely that most connections are local and a few are long-range, is called a small-world graph. The characteristic result about such graphs is that adding a small number of long-range edges to a locally connected graph reduces the typical path length between two vertices from proportional to N to proportional to the logarithm of N. The same effect accounts for the short chains of acquaintance observed between arbitrary pairs of people.

The predecessor of the algorithm treated here, called the Navigable Small World graph, obtains this structure without explicitly arranging for it. Vertices are inserted one at a time, and each new vertex is connected to the nearest vertices present at the moment of its insertion. A vertex inserted early, when the graph holds only a handful of points, is connected to whatever exists, which in absolute terms is far away. These early connections become the long-range edges. Vertices inserted later, into a populated graph, receive short edges. The required mixture arises from the construction order.

5.3 The remaining inefficiency

The Navigable Small World graph works, but it wastes effort. Long edges and short edges are mixed together in a single graph, so at every step of the walk the algorithm computes distances to long-range neighbours. This is useful at the beginning, when the search is crossing the space. It is useless at the end, when the search has already arrived in the correct region and is attempting to refine its answer, yet the cost continues to be paid at every step.

5.4 Separating edges into layers

The Hierarchical Navigable Small World structure separates edges by length, placing them in different layers.

Each vertex is assigned an integer called its height. Every vertex appears in layer zero. A vertex appears in layer l if and only if its height is at least l. Heights are assigned randomly from a distribution which decays geometrically, so that layer zero contains all N vertices, layer one contains roughly N/M of them, layer two roughly N/M², and so on, until the topmost layer contains a single vertex.

Within each layer, a proximity graph is built over only the vertices present in that layer. Because the upper layers contain few vertices scattered across the whole space, their edges are necessarily long. Because layer zero contains every vertex, its edges are short.

A search begins at the single vertex of the topmost layer. Within each layer above the bottom it performs greedy routing with ef set to 1, which is cheap, until no improvement is available. The vertex where it stops becomes the starting point for the layer below. On reaching layer zero, it performs the full beam search of Chapter 4 with the user's chosen ef.

The result is a search that is coarse at first and fine at the end. The upper layers place the search in the correct region using very few distance computations. The bottom layer performs the careful local search that determines the accuracy of the final answer.

5.5 The correspondence with skip lists

This structure is not new in its essentials. It is a skip list applied to a graph.

A skip list stores a sorted sequence as a linked list, together with additional lists above it. Each element appears in the bottom list; each element appears in the list above with some fixed probability, and so on upward. A search travels as far as it can along a sparse upper list, taking long strides, then drops to a denser list to refine, and repeats. The expected search cost is logarithmic, and the heights of elements are chosen by repeated random trials.

The hierarchical proximity graph applies exactly this mechanism, with ordering by geometric proximity replacing ordering by key. The random height assignment, the descent from sparse to dense, and the resulting logarithmic cost all correspond directly.

5.6 How heights are assigned, with the arithmetic

The height of a vertex is obtained as follows. Draw a random number u uniformly from the open interval between 0 and 1. Compute the negative of its natural logarithm, multiply by a constant written mL, and take the floor of the result. That is the height.

The constant mL is set to the reciprocal of the natural logarithm of M. It is worth confirming that this produces the intended layer occupancies, since the derivation is short and makes the choice of constant intelligible rather than arbitrary.

A vertex reaches height 1 or more when the computed value is at least 1, that is, when the negative logarithm of u multiplied by mL is at least 1. Rearranging, this requires that the negative logarithm of u be at least 1/mL, which with mL equal to 1/ln(M) means at least ln(M). Equivalently, u must be at most e raised to the power of minus ln(M), which is exactly 1/M.

Since u is uniform on the interval from 0 to 1, the probability that u is at most 1/M is exactly 1/M. Therefore precisely a fraction 1/M of vertices reach layer one, and by the same argument a fraction 1/M of those reach layer two, and so forth.

With M set to 16, one vertex in sixteen appears in layer one, one in two hundred and fifty six in layer two, and so on. A collection of one million vertices produces roughly five layers. This is why the search reaches the bottom quickly.

Chapter 6. Choosing which edges a vertex should have

6.1 The obvious rule, and why it is wrong

One question remains: given that a vertex is to receive M edges, which M vertices should they join it to?

The obvious answer is the M nearest. This answer is defective, it is the single most commonly omitted refinement in implementations of this algorithm, and understanding its failure is worth more than any other single insight in this document.

Consider a group of points packed tightly together and well separated from everything else. Under the nearest-M rule, every vertex inside the group spends all M of its edges on other members of the group, because those are the nearest vertices available. The group therefore has almost no edges connecting it to the rest of the graph.

Two consequences follow, and both are damaging. A search that has not entered the group cannot find its way in, because scarcely any edge from outside points into it. A search that has entered cannot leave, because scarcely any edge from inside points out. Portions of the collection become unreachable from the entry point, and no increase in ef can repair this, because ef governs how widely the search explores along existing edges and cannot create an edge that was never built.

The defect is structural. It is present in the graph, not in the search.

6.2 The diversity rule

The correct rule prefers variety of direction over nearness.

Consider the candidates in increasing order of distance from the vertex being connected. Accept a candidate only if it is nearer to that vertex than it is to any candidate already accepted. If a candidate is nearer to some already-accepted neighbour than it is to the vertex itself, reject it.

The reasoning behind the rejection is this. If candidate c lies nearer to an accepted neighbour r than to the vertex v, then c occupies a direction that r already covers. A search arriving at v can reach c by stepping to r and then onward, so a direct edge from v to c would be largely redundant. The edge budget is better spent on a direction not yet represented.

The effect is that a vertex's M edges are distributed among distinct directions, rather than being concentrated wherever the points happen to be densest. Dense groups therefore retain edges leading out of them, and remain reachable from outside.

6.3 A small illustration

Suppose vertex v is to receive two edges, and there are three candidates: a and b lying close together to the north of v, and c lying to the south, slightly further away than either.

Under the nearest-M rule, v connects to a and b. Both edges point north. Nothing to the south of v is reachable from it, and the search cannot travel southward through v at all.

Under the diversity rule, a is accepted first, being nearest. Then b is considered: b is nearer to a than it is to v, since a and b lie close together and both are some distance from v, so b is rejected as covering a direction already served. Then c is considered: c is further from a than it is from v, since they are on opposite sides, so c is accepted. Vertex v ends with one edge north and one edge south, and the graph remains navigable in both directions.

This rule appears as Algorithm 4 in the reference paper. Two optional refinements are described there: extending the candidate pool with the neighbours of the candidates before applying the rule, which helps when the data is strongly clustered, and retaining some rejected candidates when the rule proves so restrictive that a vertex would be left with too few edges.

Chapter 7. The parameters

Four quantities control the system. Each represents a definite trade, and being able to state each trade precisely is a fair test of whether the preceding chapters have been understood.

M is the number of edges each vertex receives in the layers above zero. A larger M produces a better connected graph in which higher recall is attainable, at the cost of more memory, since each vertex stores M neighbour identifiers per layer, and slower steps, since each step computes the distance to every neighbour. Values between 8 and 48 are usual.

Mmax0 is the same quantity for layer zero, conventionally set to twice M. Layer zero is where the final and most accuracy-determining phase of every search takes place, and where a failure of connectivity does the most damage, so it is given a larger budget.

efConstruction is the value of ef used during insertion, when the algorithm is searching for good neighbours for a new vertex. A larger value produces a better graph and a slower build. It is paid once.

ef is the value used during a query. It is supplied per query rather than fixed when the index is built, which means one index can serve both a fast approximate query and a slow accurate one. This distinction between build-time and query-time parameters has direct operational consequences and is a common subject of examination.

mL is the height-distribution constant of Section 5.6, set to the reciprocal of the natural logarithm of M.

Chapter 8. What must be understood about the hardware

The algorithm as described is complete. An implementation following it faithfully will nevertheless be several times slower than it should be unless a few facts about memory are taken into account. These facts are not usually covered in an algorithms course, so they are set out here.

8.1 Memory is not uniformly fast

A processor does not read individual bytes from main memory. It reads fixed-size blocks called cache lines, which on current hardware are 64 bytes, and holds recently used lines in small fast memories called caches.

Reading a value already held in cache takes on the order of one to twenty processor cycles. Reading a value that must be fetched from main memory takes on the order of two to three hundred cycles. The difference is roughly two orders of magnitude, and it is the dominant consideration in the performance of this system.

Sixty-four bytes holds sixteen four-byte floating-point numbers. A vector of dimension 128 therefore occupies eight cache lines.

8.2 The consequence for layout

Two arrangements of the same one million vectors behave very differently.

In the first, all vectors are stored in one continuous block of memory, with vector number i beginning at offset i × 128. Reading a vector fetches eight consecutive cache lines. Consecutive access is a pattern the hardware recognises, and it will fetch subsequent lines before they are requested, a mechanism called prefetching.

In the second, each vector is allocated separately, as occurs when a container of containers is used. The vectors are scattered across memory at unrelated addresses. Reading one requires first reading a pointer, then following it to an arbitrary location, which the hardware cannot predict and therefore cannot prefetch.

The second arrangement is straightforward to write and is the natural thing to reach for. It costs somewhere between five and ten times the performance of the first. The cost does not appear in any one place that profiling will identify as the culprit; it is distributed across every operation that touches a vector. Correcting it once the system is complete requires revising nearly everything.

The same reasoning applies to the adjacency lists of the graph, which are read constantly during search and should likewise occupy short continuous runs.

This is the reason Stage 0 of the implementation, which appears to be trivial file-reading work, is treated as a design decision.

8.3 Performing several operations at once

Processors provide instructions operating on several values simultaneously, a facility called SIMD, for single instruction, multiple data. The AVX2 instruction set available on current x86 processors operates on eight four-byte floating-point values per instruction.

Since the distance computation performs the same operation on each of 128 coordinates independently, it is well suited to this treatment, and the compiler will often apply it automatically when the loop is written simply and the data is laid out continuously.

One point is not obvious. A direct implementation accumulates the running total into a single variable, which means each addition must wait for the previous addition to finish before it can begin. The processor is then limited by the time an addition takes to complete rather than by how many it can start. Keeping two to four separate running totals, and adding them together only at the end, allows these chains to proceed in parallel and yields a further substantial improvement.

Chapter 9. The data

The standard benchmark is the SIFT collection published by the TEXMEX research group. Two subsets are required, and both should be downloaded before implementation begins.

The smaller subset, siftsmall, contains 10,000 vectors, 100 queries and the correct 100 nearest neighbours of each query. All development should be carried out on this subset. Building an index over it completes almost instantly, which permits a rate of experimentation that is impossible with the full collection. Attempting to develop against a million vectors means waiting minutes to discover faults that a ten-thousand-vector run would have revealed in a second.

The larger subset, ANN_SIFT1M, contains 1,000,000 vectors of dimension 128, 10,000 queries and ground truth. It occupies approximately 516 megabytes. All reported results should come from this collection.

The file format is called .fvecs and is deliberately simple. Each vector is stored as a four-byte integer giving the dimension, followed by that many four-byte floating-point values, and the vectors are placed one after another with no header. The ground truth files use the extension .ivecs and have exactly the same layout with four-byte integers as the payload. Reading these files requires no library of any kind, which is why this format is preferred here over the HDF5 packaging used by the public benchmark suite.

Chapter 10. The implementation, stage by stage

Each stage states what is to be built, the decisions it involves, and how to establish that it is correct before proceeding.

Stage 0. Reading the data, and the layout decision

Objective. Read the vector and ground-truth files into memory.

Decision. Store all vectors in a single continuous allocation, with the vector of index i beginning at offset i × D, aligned to a 64-byte boundary. The reasoning is given in Section 8.2 and should be re-read if it is not fresh.

Verification. On siftsmall, confirm that exactly 10,000 vectors, 100 queries and 100 ground-truth records are read, and that every vector has dimension 128.

Stage 1. Exhaustive search, which becomes your reference

Objective. Implement exact k-nearest-neighbour search by scanning the whole collection.

This component is not scaffolding to be discarded. It produces the correct answers against which every approximate result will be measured, and it is the only component whose correctness can be established without reference to anything else.

Verification, which must not be skipped. Your exhaustive search must reproduce the supplied ground truth for siftsmall exactly, allowing for arbitrary ordering among points at equal distance. If it does not, the fault lies in the file reader or in the distance function. Both must be eliminated before a single line of graph code is written. If you proceed without this check, you will later observe poor recall and begin searching for the fault in your graph construction, where it is not.

Record. Queries per second on the full million-vector collection. This is the baseline against which the final result will be expressed, and you should expect to improve on it by roughly a factor of a thousand.

Stage 2. The distance function

Objective. Implement squared Euclidean distance, first plainly and then using the vector instructions described in Section 8.3.

Verification. The plain and vectorised versions must agree to within floating-point tolerance across at least a million randomly generated pairs. They will not agree exactly, because floating point addition is not associative and the two versions add the terms in different orders. A small tolerance is correct here; exact equality is not the right test.

Record. Repeat the Stage 1 measurement. The improvement obtained from layout and vectorisation alone is worth documenting, since it is a real result and demonstrates the material of Chapter 8.

Stage 3. The beam search, built and tested alone

Objective. Implement the procedure of Sections 4.4 and 4.5, operating on a single graph layer, with no graph construction code present.

This one procedure accounts for most of the algorithm's behaviour and most of its running time. It is to be built and proven correct before anything is built on top of it.

A necessary technique. The procedure must record which vertices it has already examined. The direct approach, allocating a boolean array of size N per query or clearing one between queries, will dominate the running time completely, since a single search examines only a few hundred vertices out of a million and would be paying for a million-element operation each time.

The standard solution is to keep one array of 32-bit integers of length N, together with a counter that is incremented once per query. A vertex counts as visited when its stored value equals the current counter value. Marking a vertex means storing the counter. Nothing is ever cleared, because incrementing the counter invalidates all previous marks at a stroke. Both marking and testing are single operations.

Verification. Build a temporary graph in which every vertex is joined to its exact k nearest neighbours, computed using your Stage 1 implementation. Such a graph is about as easy to navigate as a graph can be. Running your beam search on it with a generous ef must produce recall@10 of at least 0.99. If it does not, the fault is in the procedure itself, and finding it now is far easier than finding it after construction code exists to obscure the diagnosis.

Stage 4. Building a single-layer graph

Objective. Construct a Navigable Small World graph by inserting vertices one at a time, with no hierarchy yet.

Procedure. For each new vertex, run the Stage 3 search from a fixed entry point with ef set to efConstruction. This yields a set of candidate neighbours. Select M of them and connect. Add each edge in both directions, since a neighbour that cannot reach back to the new vertex provides no help in navigating toward it. Adding the reverse edge may push an existing vertex over its edge budget, in which case that vertex's neighbour set must be reduced.

Decision. Store the adjacency as a flat array giving a fixed number of neighbour slots per vertex, together with a count of how many are used. The reasoning parallels Stage 0: a vertex's neighbour list is read at every step of every search and should arrive in one or two cache lines.

A frequent error. When a vertex exceeds its budget, the tempting remedy is to discard its single furthest edge. This is not the same as re-applying the neighbour selection procedure to its complete set of neighbours, and the difference degrades recall in a way that looks like a parameter tuning problem and is not one. Apply the selection procedure properly.

Verification. Plot recall@10 against ef on siftsmall. The curve should rise and then level off below the accuracy you want. That ceiling is the deficiency that Stages 5 and 6 exist to remove, and seeing it for yourself is the point of this stage.

Stage 5. Adding the hierarchy

Objective. Extend Stage 4 to multiple layers as described in Section 5.4.

Procedure. Assign each vertex a height by the method of Section 5.6. The vertex of greatest height is the global entry point. When inserting a vertex of height l, first descend from the top layer down to layer l + 1 using greedy routing with ef set to 1, which serves only to find a good starting position. Then, from layer l down to layer zero, run the full search with ef set to efConstruction and establish edges. Layer zero uses Mmax0 rather than M.

Verification. Produce a histogram of vertex heights and confirm that the population falls by a factor of about M at each successive layer, as derived in Section 5.6. If a large fraction of your vertices reach layer one, either mL is wrong or your random number generator is being reset, most commonly by constructing it inside the insertion function so that every vertex receives the same random value. This check takes a few minutes and catches a class of fault that is otherwise very hard to notice, because the index still works and is merely slow.

Record. At equal recall, both the query rate and the average number of distance computations per query should improve relative to Stage 4.

Stage 6. The diversity rule

Objective. Replace nearest-M selection with the procedure of Section 6.2.

Verification. Run both selection rules and plot the two recall-against-throughput curves on the same axes. The diversity rule should show a clear advantage in the high-recall region, which is the region in which such systems are actually deployed.

Keep this comparison. It is the clearest evidence you will produce that you understood the algorithm rather than transcribed it, and it belongs at the top of whatever you write about the project.

Stage 7. The parameter study and the final result

Objective. Characterise the finished index across its operating range on the full collection.

Procedure. Vary M over 8, 16, 32 and 48; efConstruction over 100, 200 and 500; and ef over a range from roughly 10 to 800. For each configuration record recall@10 and sustained throughput, and plot the curve.

Comparison. Install the hnswlib library and run it on the same machine, the same data and the same thread count, and plot its curve alongside yours. Reaching throughput within roughly a factor of two of it at recall 0.95 is a strong result for an independent implementation.

Measuring honestly. Accumulate your search results into a value that is printed at the end. Without this, an optimising compiler may observe that the results are never used and remove the search entirely, and you will report an extraordinary figure that measures nothing. Report the minimum across several runs rather than the average, since you are trying to measure the computation rather than the interference from the rest of the operating system. State the thread count and whether the caches were warm.

Stage 8. Deletion

Objective. Support removing vertices from a built index.

This stage is included precisely because the algorithm handles deletion badly. Being able to explain a genuine weakness of a structure is worth more in conversation than another feature that works without incident.

Implement and compare two approaches. The first marks a vertex as deleted and filters it out of results while leaving it in the graph. This remains correct, but deleted vertices continue to consume steps during every search, and as deletions accumulate an increasing proportion of each vertex's edges lead to vertices that can never be returned, so navigability degrades. The second repairs the graph when a vertex is removed by connecting its former neighbours to one another.

Verification. Delete 10 per cent, then 30, then 50, and plot recall over the surviving points for both strategies.

Stage 9. Saving the index

Objective. Write the index to a file arranged so that it can be memory-mapped, that is, made available to the program by the operating system without being copied or parsed.

Decision. The file should contain a header recording the dimension, the vertex count, M, the entry point and the per-layer offsets, followed by the adjacency data and then the vectors. It must contain offsets rather than pointers, since the operating system may place the file at any address and stored pointers would be meaningless there.

Verification. Save, load and confirm that queries return results identical to those from the in-memory index.

Record. Time to map the file compared with time to rebuild the index from scratch.

Stage 10. Optional extensions, if you want more

Product quantisation reduces the memory the vectors occupy. Each 128-dimensional vector is cut into, say, eight pieces of 16 dimensions. For each piece position, the k-means clustering algorithm is run across the whole collection to find 256 representative values, and each piece is then replaced by the single byte identifying which representative it is closest to. A vector of 512 bytes becomes 8 bytes. Distances are computed approximately from small precomputed tables. This is the technique that allows indexes of a billion vectors to fit in memory, and combining it with a graph index is standard practice in production systems.

Filtered search answers queries of the form "the nearest neighbours among those items satisfying some condition", such as being in stock or belonging to a particular user. This appears to be a small addition and is not. The graph's navigability depends on edges that the filter may forbid the search to follow, so a restrictive condition can leave the search stranded in a region from which no permitted edge leads anywhere useful. This remains a subject of active research rather than a solved problem, which makes it interesting ground if you want to work on something unsettled.

Parallel construction, using one lock per vertex, should be attempted only after correctness is fully established.

Chapter 11. Suggested order of work

Your current machine is available until 5 October. Both datasets should be downloaded immediately. Stages 0 through 2 should be finished and Stage 3 verified before that date.

The reason for choosing that particular stopping point is that Stage 3 is the structural core of the project. Once the beam search has been proven correct against an exact nearest-neighbour graph, every later stage adds to it without modifying it, and the work can be resumed on any machine without needing to re-establish context.

Stages 4 through 6 contain the real algorithmic content and should be expected to take the longest. Stage 7 produces the result you show to other people.

Chapter 12. Errors that will cost you days

These are listed roughly in order of how much time they consume when committed.

Storing vectors or adjacency lists in a container of containers. The cost is large, spread across the whole program, and awkward to correct once everything is built. See Section 8.2.

Allocating or clearing the visited set on every query, instead of using the counter technique described in Stage 3.

Ordering one of the two priority queues the wrong way round. The search still runs, still terminates, and still returns plausible results that are quietly worse than they should be. Nothing fails, so there is no natural point at which to start investigating. The only symptom is a recall figure lower than you expected, which is easily mistaken for a tuning problem.

Leaving the square root in the distance function. See Section 1.3.

Looking for the fault in the wrong component. Investigating graph construction when the defect is in the file reader. This is exactly what the Stage 1 verification exists to prevent.

Discarding the furthest edge when pruning rather than re-applying the selection procedure.

Measuring recall against your own earlier output rather than against ground truth. This establishes only that your program agrees with itself.

Benchmarking a loop the compiler removed because nothing consumed its results.

Developing against the full collection. Every experiment then costs minutes that ought to cost seconds. The small subset exists for this purpose.

Chapter 13. Questions to test yourself

You should be able to answer each of these without notes once the project is finished. They are close to what someone familiar with the field would ask.

Why do grid and tree methods fail in 128 dimensions? Give both arguments: the one about the number of neighbouring regions, and the one about distances ceasing to distinguish.

What ordinary data structure is the layer hierarchy an instance of, and in what respects does the correspondence hold?

What exactly goes wrong when neighbours are chosen purely by nearness? Describe the failure in terms of the connectivity of the graph, not in terms of a measured recall figure.

Why does layer zero receive twice the edge budget of the other layers?

Which parameters are fixed when the index is built, and which are chosen per query? What practical flexibility follows from that distinction?

How many bytes does each vector cost in your index, counting the vector itself and its edges separately? If you had to halve that, which would you attack first and by what method?

Why is deletion difficult in this structure, and what are the two available strategies?

Which part of a query is limited by memory access rather than by arithmetic? Explain how you established that by measurement rather than by reasoning about it.