anydist

Approximate nearest-neighbor search on arbitrary data with arbitrary dissimilarity functions, using Hierarchical Navigable Small World graphs (HNSW; Malkov & Yashunin).

Unlike most HNSW libraries, which index numeric vectors under a fixed set of metrics, this one accepts any Python objects and any dissimilarity function you write — including non-metric ones. It also has the usual fast path: a set of built-in metrics over numeric vectors, computed in SIMD C kernels with no Python on the hot path, with parallel bulk insertion and parallel batch queries. Between the two there is a third door: hand it a C function pointer — from numba, ctypes, Cython or an extension module of your own — and your dissimilarity runs on that same GIL-free hot path (see C metrics).

The index stays live. You can add to it and delete from it after it is built, interleaved with queries, without a rebuild — which not every library offers: nmslib and annoy are batch indexes, where a point added after the build costs a full rebuild to become findable (see Living indexes).

Two implementations:

• HNSWFast — the one you want: a Cython implementation over typed, array-backed data structures. HNSW is an alias for it, and a normal install builds it.
• HNSWPure — the pure-Python, dict-based reference implementation, no dependencies beyond the standard library. It exists to document the algorithm readably and to check the fast path against, not as a deployment option: with a built-in metric it is ~30x slower to build and ~40x slower to query (measured on 5k 20-dimensional vectors), and even with a Python callback — where the callback itself dominates — it is ~1.5x slower.

HNSW never falls back to HNSWPure. Without the compiled extension there is no HNSW, and asking for it raises ImportError saying so — an order of magnitude is not something to lose silently to a build that did not happen. If you want the pure implementation, name it:

from anydist import HNSWPure as HNSW      # deliberately, penalty and all

In a hurry? Quickstart is two code blocks, and the API reference documents every constructor argument and method in one place — go straight there if you would rather read signatures than prose.

This code has been developed together with fishdbc, an incremental density-based clustering algorithm that leverages this library for scalability and allowing arbitrary dissimilarity functions.

Installation

pip install anydist

or, from a checkout of this repository:

pip install .

Windows is not supported: MSVC, the compiler Windows Pythons are built with, lacks the atomic builtins the parallel code needs. Under WSL, the Linux wheels work as usual.

PyPI has prebuilt wheels for the common platforms. Anywhere else, and from a checkout, pip compiles the Cython extension, so you need a C compiler with OpenMP; on macOS, brew install libomp first (or point LIBOMP_PREFIX at another libomp). That extension is what gives you HNSWFast. A failed build shows up the first time you ask for HNSW, since there is no fallback to hide it, but you can check directly:

import anydist
assert anydist.HNSWFast is not None

On x86-64 Linux, the portable build already selects AVX2 kernels at load time on CPUs that have them, so it is not leaving the obvious speedup on the table. For a further boost tied to the machine you build on (at the cost of a non-portable binary):

HNSW_NATIVE=1 pip install .

Quickstart

Arbitrary objects, your own dissimilarity function:

from anydist import HNSW

def jaccard(a, b):
    return 1 - len(a & b) / len(a | b)

index = HNSW(jaccard)                       # the dissimilarity is the 1st argument
index.update(my_sets)                       # insert them all
neighbors = index.search(query_set, k=10)   # [(id, distance), ...], nearest first

Ids are insertion order: the id in a result is the position of that element in the sequence you inserted. The index is a read-only sequence of those elements, so index[id] gives one back, len(index) counts them, and iteration and slicing work.

Numeric vectors and a built-in metric — no Python callback on the hot path:

index = HNSW('euclidean')            # vector length inferred from the first row
index.update(X)                      # X: (n, dim) array, inserted in parallel
I, D = index.search_batch(Q, k=10)   # (n_queries, k) ids and distances, in parallel

With a built-in metric, update inserts across all cores; with a callback metric it stays serial by default, since only you can say whether your function is safe to run on several threads — pass num_threads=N to opt in. add(elem) is still there for genuinely incremental use, where you interleave inserts and queries.

index[id] works in both modes, but they store elements differently: a callback metric keeps your objects, while a built-in metric keeps rows of a typed matrix and hands you back a copy of one — rounded to float32, since that is the default storage there (see float32).

Multithreading

update and search_batch both take num_threads: None (the default) picks for you, a positive value is a thread count, and -1 (or any value ≤ 0) means one thread per CPU.

The rule is that anything we wrote defaults to every core, anything you wrote defaults to one, because only you can say whether your dissimilarity function d — the callable you passed as metric — is thread-safe. What you get once you opt in depends on the metric, because the GIL does:

• Built-in metric — ours, and stateless, so the kernels are nogil and the threads are real on any interpreter build. Measured on 12 CPUs: 5.5x build / 5.0x query at 50k x 32, and 4.3x / 3.7x at Fashion-MNIST's 60k x 784. Don't read those as an efficiency: the machine's cores are not all the same speed, so what a given box can offer is its own question.
• C metric (a function pointer you supplied) — the call is nogil too, so opting in is the only way a user-written distance gets real threads on a stock interpreter: measured 4.3x at one-per-core (n=4k, dim=16, a deliberately slow jit_metric). It still defaults to serial, because nogil buys the ability to run your d in parallel, not the safety of doing so — mutable state reached through user_data, a memo table, a counter, all race, with no GIL here to accidentally serialise them.
• Callback metric — the callback runs under a thread state re-attached per distance. On a free-threaded interpreter that is real parallelism; on a stock GIL build it pays off only when each distance holds the GIL briefly and releases it for a long time — measured 3.0x at 2 threads for a 40k-dimensional NumPy distance, but slower than serial at 128 dimensions, where contention and the thread-state attach cost more than the op. Measure before believing.

These figures are indicative, not a spec — see About these numbers for the method and its resolution limits.

index.update(X)                              # built-in metric: all cores
index.update(objs)                           # your metric: one thread
index.update(objs, num_threads=-1)           # opt in: "my d is thread-safe"
index.update(objs, num_threads=4)            # opt in, four threads
I, D = index.search_batch(Q, num_threads=1)  # force a single-threaded query

search_batch is worth calling for a callback metric even serially: it is the same graph walk as a search loop, but it hands back the whole batch as two (n_queries, k) arrays — neighbor ids and their distances — instead of a list of lists.

chunksize, the constructor's or a call's own, interacts with this: workers take chunks of consecutive ids, which is what keeps parallel insertion correct when the input is order-local (sorted, time-series, grouped records) — two near neighbors inserted at the same instant never see each other, and the edge between them is lost. The default adapts to the batch; raise it if your input is strongly ordered.

Writing a callback metric that scales

On a free-threaded interpreter there is one rule, and it is not about how you write the function — it is about what the function reads:

A parallel metric must not read a shared mutable object on every call. Reading self.something or a closed-over variable costs 3-4x at eight threads. Read parameters from module globals or default arguments instead.

CPython gives some objects immortal or deferred reference counts, so touching them from many threads is free: None, small ints, literals, module-level functions, and all classes (including ones built at runtime by type(), namedtuple, or a class statement inside a function). Everything else is ordinary, and an object that every thread increfs on every distance call turns one cache line into a contention point. Measured on the same laptop (n=10 000 tuples of 8 floats, ef=25, 8 threads, free-threaded 3.13):

metric written as	speedup at 8 threads
module-level def, reads only its arguments	3.53x
closure that reads its captured variable	0.82x — slower than serial
closure that ignores it	3.08x
bound method that reads self	1.19x
bound method that ignores self	2.94x

Note the pairs: same callable form, same arithmetic, only the read differs. The form is irrelevant — functools.partial and lambdas are fine. So a metric with parameters should carry them like this:

# slow in parallel: `self.p` and the closure cell are one shared object each
class Metric:
    def __call__(self, a, b):
        return sum((x - y) ** self.p for x, y in zip(a, b))

# all of these scale (3.3-3.5x measured):
P = 2
def metric(a, b):                    # module-level global
    return sum((x - y) ** P for x, y in zip(a, b))

def make_metric(p):
    def metric(a, b, p=p):           # default-argument binding
        return sum((x - y) ** p for x, y in zip(a, b))
    return metric

Per-element objects are not a problem: your data points are thousands of distinct objects, so the traffic spreads. It is the single shared one that hurts.

If threading makes things slower, suspect numpy. There is a bug in numpy (numpy#32298) where an array that has been through pickle, copy.deepcopy, or a joblib cache ends up sharing one small internal object that every thread has to synchronise on. The array is otherwise identical, and single-threaded work is unaffected — but on a free-threaded build the contention can turn a speedup into a slowdown.

Usually you need not care: update() detects and repairs it whenever you hand it an array. The case it cannot reach is a list of rows you extracted yourself, because then the library never sees the array. Repair it at the source:

from anydist import canonicalize_dtype
rows = list(canonicalize_dtype(X))       # X came from pickle / joblib / a cache
index.update(rows, num_threads=-1)

So: parallel inserts slower than serial, on data that came out of a cache, is worth one try of the two lines above — see BUGS.md, which also lists the other issue worth knowing about.

Built-in metrics

Passing a name instead of a callable computes the distance in a SIMD C kernel with no Python on the hot path. Available names (aliases in parentheses):

Name	Distance	Notes
euclidean (l2)	sqrt(sum (x-y)^2)
l2sqr (sqeuclidean)	sum (x-y)^2	ranks identically to l2, skips the sqrt
cosine (cosinesimil)	1 - cos_sim	normalizes for you; the value scipy and sklearn report
angular (angulardist)	arccos(cos_sim) / pi	normalizes for you; ranks as cosine does, and obeys the triangle inequality
dot	1 - dot(x, y)	normalizes nothing — cosine distance only if you pre-normalized
negdotprod	-dot(x, y)	the same space as dot, shifted by 1
l1 (manhattan)	sum |x-y|
linf (chebyshev)	max |x-y|

The names follow nmslib's dense-vector spaces. dot and negdotprod differ by an additive constant, and every decision the search makes is a comparison, so they build identical graphs and return identical orders — choose between them by the number you want back. dot is non-negative on unit-length rows, which a downstream consumer of the weights may need; negdotprod is what code ported from nmslib expects.

Cosine distance

Three names get you cosine-like behavior, and they differ in whether they normalize for you and in what number comes back. All three return the same neighbor lists on unit-length rows, so this is a choice about values and cost, not about quality.

Just want cosine distance? Use cosine. It divides by the length of each of the two vectors itself, so your rows do not have to be unit-length, and it returns 1 - cos_sim — the same value scipy, scikit-learn, nmslib (cosinesimil) and hnswlib report.

index = HNSW('cosine')
index.update(X)                          # any lengths; normalizing is on us
index.search(q, k=10)                    # distances are 1 - cos(q, x), in [0, 2]

Need the triangle inequality? Use angular. Same normalization, but it returns the angle between the vectors scaled to [0, 1] instead of 1 - cos_sim. Since a smaller angle is exactly a larger cosine, it ranks candidates identically to cosine; unlike cosine, it is a true metric on directions. The extra arccos is not free — it is the only libm call in any built-in kernel, and it costs roughly a third of query time at 50 dimensions — so reach for angular when something downstream needs metric behavior, not by default.

Already normalized, and want the cheapest option? Use dot. It normalizes nothing, so on unit rows 1 - dot(x, y) is cosine distance, computed with a single accumulator instead of three.

import numpy as np
from anydist import HNSW

X = np.asarray(X, dtype=np.float64)
X /= np.linalg.norm(X, axis=1, keepdims=True)   # unit rows (drop zero rows first)

index = HNSW('dot')
index.update(X)

q = q / np.linalg.norm(q)                # the query too: dot normalizes nothing
index.search(q, k=10)

A zero row has no direction, so cosine and angular both report distance 0 for it rather than a NaN.

C metrics: numba, ctypes, an extension module

A built-in metric is fast because it never enters Python; a callback is flexible because it is Python. This section is for when you want both — and it asks more of you than the two options above, so skip it unless a Python callback has turned out to be your bottleneck.

Give HNSW a C function pointer instead of a Python callable, and your distance runs with no GIL and no Python on the hot path, exactly where a built-in metric runs. Parallel insertion and parallel queries then become possible on a stock interpreter — but you have to ask for them:

index.update(X, num_threads=-1)              # -1 = one thread per core
I, D = index.search_batch(Q, k=10, num_threads=-1)

A C metric still defaults to one thread, like a Python callback and unlike a built-in metric. nogil buys the ability to run your distance in parallel, not the safety of doing so: state reached through user_data, a memo table or a counter all race, and there is no GIL here to serialise them by accident. Only you can say whether your d is thread-safe, so the opt-in is yours to make.

The signature is fixed, and it is the one scipy's LowLevelCallable users will recognise:

double f(const double *a, const double *b, Py_ssize_t n, void *user_data);

Writing one in Python: jit_metric

Nobody wants to spell that signature in numba by hand, so the decorator generates it from the two-argument form:

from anydist import HNSW, jit_metric

@jit_metric                          # needs numba installed
def manhattan(a, b):                 # a, b: float64 arrays of length dim
    s = 0.0
    for i in range(a.shape[0]):
        s += abs(a[i] - b[i])
    return s

index = HNSW(manhattan)
index.update(X, num_threads=-1)      # opt in: "my d is thread-safe"
I, D = index.search_batch(Q, k=10, num_threads=-1)

It returns a numba.cfunc, which HNSW recognises by its .address. The body must compile in nopython mode — loops over the two arrays are the idiom, arbitrary Python is not; if you need arbitrary Python, pass the plain function and pay the GIL per distance. @jit_metric(cache=True) persists the compiled kernel across processes, which matters for a benchmark harness that runs one cell per process; it requires the function to live in an importable file.

What the pointer buys, measured on the same 12 CPUs (stock-GIL CPython 3.14, n=5000, dim=32, m=16, ef=100, k=10) with the same distance written both ways — the callback being the obvious numpy one-liner:

metric	threads	build	query
Python callback	1	5.34 s	1.63 ms/q
C metric (jit_metric)	1	0.20 s	0.032 ms/q
C metric	12	0.044 s	0.008 ms/q

The 12-thread row is opt-in (num_threads=-1) — a C metric is serial by default, see Multithreading — and it is safe here because this metric is a pure function of its two arguments.

26x on build and 50x on query before threading, 4.6x and 4.3x more from it. The serial ratio is the part that moves with your metric: a callback whose body does real work narrows it, one written in pure Python widens it. The threading factor is the part a callback cannot have on a stock build at all.

Pointing it at a library you already have

The pointer is duck-typed, so nothing here is imported unless you use it. Any of these shapes is accepted:

• scipy LowLevelCallable — .function is a PyCapsule, .user_data carries parameters. scipy already normalises numba, ctypes and Cython into it, so this one shape covers all three.
• numba cfunc — has .address.
• ctypes function pointer — how an existing .so reaches Python.
• anything else exposing .address or a .function capsule (plus an optional .user_data), including a PyCapsule from your own extension module.

import ctypes

dbl = ctypes.POINTER(ctypes.c_double)
lib = ctypes.CDLL("./libdtw.so")
lib.dtw_metric.restype = ctypes.c_double
lib.dtw_metric.argtypes = [dbl, dbl, ctypes.c_ssize_t, ctypes.c_void_p]

index = HNSW(lib.dtw_metric)         # DTW, in C, across all cores

Parameters go through user_data, which is read from every one of those shapes: keep them in a struct (or a single ctypes cell) that outlives the index, and expose its address as .user_data next to the pointer. jit_metric does not expose user_data — a closure over module-level constants covers the parameterised case with less ceremony, and a hand-built cfunc or a LowLevelCallable is there when you want the pointer itself.

One thing to know before you take this door: getting the signature wrong crashes the process instead of raising an exception. How much protection you get depends on what you hand over:

All three forms declare their signature somewhere, and all three are checked against double (const double *, const double *, Py_ssize_t, void *), raising a CSignatureWarning if they disagree:

• A capsule carries it in its name — a scipy.LowLevelCallable, or anything else whose .function is a named PyCapsule. Spelling is not the test: const, spacing, parameter names and Py_ssize_t/ssize_t/ npy_intp all pass; the types have to line up.
• A ctypes function pointer carries it as _restype_/_argtypes_, and a numba cfunc carries the same on its .ctypes view. Either pointer spelling is accepted for the vectors — POINTER(c_double) or a bare c_void_p. jit_metric builds the signature for you, so that path cannot disagree in the first place.

What is genuinely unchecked is a bare address: an unnamed capsule, or a pointer whose declaration you never wrote down. There the crash contract stands in full.

The check warns rather than refuses, because it recognises a list of spellings rather than parsing C, and an unfamiliar typedef must not be able to block a pointer that is in fact correct. If that is your situation:

warnings.simplefilter("ignore", anydist.CSignatureWarning)

and please open an issue with the spelling, so the list can grow.

One limit is worth knowing before you write the metric rather than after: your vectors must all be the same length. n is a single dimension, inferred from the first row and passed to every call. Variable-length data stays on the callback path, which imposes nothing on what your objects are.

Otherwise a C metric follows the built-in-metric path rather than the callback one, which the API reference covers argument by argument. One of those defaults is worth a second look, because it is the only place that path assumes something about your function rather than ours: select is chosen from the vector length, on the reasoning that a short vector is cheap to re-measure. That holds for the built-in kernels, whose work is linear in the length. If yours is expensive at any length — an alignment, a search, a simulation — pass select='edges' and it will read stored weights instead of calling you O(m²) times per row.

API

Constructor

HNSW(metric, m=5, ef=200, m0=None, level_mult=None,
     level_strategy='random', keep_pruned=False, select=None, cache=None,
     float32=None, keep_workspace=False, chunksize=None)

• metric — the dissimilarity, in one of three forms (it need not satisfy the metric axioms in any of them — the name follows the scikit-learn/hdbscan convention, not a metricity requirement): a callable d(x, y) over arbitrary data; a built-in name (see above) over numeric vectors, whose length is inferred from the first one inserted; or a C function pointer over such vectors (numba cfunc, scipy.LowLevelCallable, ctypes pointer, capsule), which runs GIL-free like a built-in one — see C metrics. A pointer is recognised before callability, so an object that is both takes the fast door.

• m, ef, m0, level_mult — the paper's parameters; see Malkov & Yashunin.

• level_strategy — how insertion levels are assigned. 'random' is the paper's strategy and the best choice for approximate nearest-neighbor search; 'balanced' is a deterministic variant (an element is promoted when its neighborhood is full and no neighbor sits above it) that yields a more balanced structure.

• keep_pruned — the paper's optional keepPrunedConnections (Alg. 4, lines 15–17). Off by default.

• select — the neighbor-selection rule. 'edges' answers the paper's pruning heuristic (Alg. 4) from edges the graph already has, spending zero extra distance evaluations; 'distances' is the heuristic as written, recomputing the candidate-to-result distances; 'naive' just takes the m nearest (Alg. 3). HNSWFast defaults to 'edges' for a callback metric, and for a built-in one decides at the first vector: 'distances' while dim <= 4*m (recomputing a short vector beats fetching a stored weight), 'edges' above that — a C metric follows that same first-vector rule. HNSWPure always defaults to 'edges', since every distance there is a Python call.

  That first-vector rule takes vector length as a stand-in for what a distance costs, which is true of the built-in kernels. A C metric that is expensive out of proportion to its vector length should be given select='edges' explicitly — otherwise short rows get 'distances', which evaluates your function O(m²) times per row instead of reading stored weights.

• cache — reuse, within a single insert or query, distances computed at a previous level of the same descent. Defaults on for a callback metric (13% faster builds, measured), off for a built-in one — and off for a C metric, which is cheap for the same reason — where the check costs about what it saves.

• float32 — store built-in-metric vectors as float32. Defaults on for a built-in metric, as hnswlib and faiss do: it halves both the vector store and the bandwidth of a query's dominant read. It does round the data, so pass False when the exact float64 distances are the point. True is rejected for the other two metric forms rather than quietly overridden: a callback metric keeps Python objects and has no typed store to narrow, and a C metric reads its rows as const double *, so a float32 store would be reinterpreted rather than converted. Both leave the vectors as float64, and index[i] hands back a copy of the stored row either way.

• keep_workspace — hold the threads' working memory between calls instead of allocating it fresh each time. Only affects update() and search_batch() with more than one thread.

  Each worker needs a set of buffers sized by the number of nodes in the index, so every threaded call pays an allocation proportional to the index — not to the work you asked for. One big update() never notices. A loop of small ones pays it every time, and it grows as the index does. At 100k nodes and four threads:

  	fresh each call	kept
  update() in chunks of 200	4.01 ms per chunk	2.46 ms
  search_batch() of 64 rows	2.11 ms	0.37 ms

  Off by default because it is a trade, not a free win: the buffers stay allocated — threads × nodes × 32 bytes, or 44 with cache=True, so 370 to 500 MB at twelve threads and a million nodes — until you call release_workspace(). Peak memory is the same either way; what changes is the floor between calls. Turn it on if you build or query in a stream of small batches, leave it off otherwise, and call release_workspace() when that phase is over.

  Results are identical either way, and the setting survives pickling while the buffers themselves do not.

• chunksize — how many consecutive ids each worker takes in a threaded update() that does not pass its own; see update below for what the chunk guards. None, the default, adapts it to each batch. HNSWPure accepts and keeps it, with no parallel insertion to schedule.

Methods

• update(X, num_threads=None, ef=None, chunksize=None) — insert many elements. X is a numeric array for a built-in or C metric — anything array-like will do, and it is converted to C-contiguous float64, so integer rows, float32 rows and lists of lists are all accepted and all stored as float64 — or any sequence of objects for a callback metric. Parallel by default for a built-in metric, serial by default for the two you supply (see Multithreading). chunksize sets how many consecutive ids each worker takes: the chunk is what keeps parallel insertion correct on order-local input (sorted, time-series, grouped), where two near neighbors inserted concurrently would never see each other. None takes the constructor's chunksize, and if that is None too the chunk adapts to the batch; set it, on the index or per call, if you know your input is strongly order-local.

• add(elem, ef=None) — insert one element, returning nothing; its id is the previous len(index).

• search(q, k=None, ef=None) — the (approximate) k closest elements to q as (id, distance) pairs, sorted by distance. k=None returns all ef candidates found. Soft-deleted nodes are excluded (see below).

• search_batch(Q, k=None, ef=None, num_threads=None) — answer many queries at once, optionally in parallel. Q is (n_queries, dim) for a built-in or C metric, or a sequence of objects for a callback one; returns (indices, distances) of shape (n_queries, k) — k defaulting to ef — padded with -1/inf where fewer than k neighbors were found. Both classes have it, with the same rows, widths and padding, so I, D = index.search_batch(Q, k) reads the same either way. HNSWFast returns numpy arrays and can use threads; HNSWPure returns lists of lists, having no numpy, and loops over search, accepting num_threads for signature parity and ignoring it as its update does.

• index[i], len(index), iteration, slicing — the index is a read-only sequence of the elements you inserted, in id order (negative ids count from the end; a built-in or C metric returns copies of the stored rows). Note this is not what index[i] meant in this package's predecessor, where it yielded node i's neighbors; graphs is where those live now.

• mark_deleted(i) / unmark_deleted(i) — soft deletion (see below); index.deleted_count reports how many are currently marked.

• add_recording(elem, ef=None, harvest_cap=0) — insert one element and return every distance computed during the insertion as three parallel sequences (i, j, dist), where i is always the new id — numpy arrays from HNSWFast, and from HNSWPure the standard library's array.array (typecodes 'i', 'i', 'd'), as compact at 16 bytes a triple. A plain add records nothing and pays no overhead.

• update_recording(X, ..., harvest_cap=0) — the batch form, returning the same three sequences for the whole insertion. It is add_recording in a loop, and given the same insertion levels the two harvest the same triples. Note that HNSWFast's comes back in per-worker order rather than insertion order, so sort by i if you need it.

  harvest_cap bounds what comes back, identically on both: 0 is off, and a positive N keeps only each inserted node's N shortest triples. It does not change the index — the same distances are computed and the same graph is built either way, so it trades recorded detail for the downstream cost of consuming it. Selection is by distance alone, which can keep the wrong bridge when a nearer cluster fills the quota first; leave it off for streaming workloads with concept drift, where insertion order cannot shuffle that away. Any cap >= 1 keeps every node attached, since a node is only ever matched against already-inserted ones.

• freeze() — give up inserting in exchange for the memory only insertion needed (the stored edge weights: n*m0*8 bytes at the base layer, plus the smaller per-node blocks above it). There is no way back, as with a frozenset. Under select='distances', which stores no weights to begin with, it only installs the guard. Returns None: it freezes the index you call it on rather than handing back a frozen copy.

• release_workspace() — free the threads' working memory that keep_workspace=True holds on to. Safe to call at any time, including when nothing is held; the next threaded call simply allocates again. Present on both classes, so code holding either can call it without asking which.

• graphs — the adjacency as a list of {node: {neighbor: dist}}, one per level, for inspection.

Living indexes

An index you can keep changing is not universal among ANN libraries. Adding 2,000 vectors one at a time to a 20k-node index, with a k=10 query every 100 inserts (dim=32, m=16, serial):

library	insert after build	μs/insert	μs/delete
this library	yes	61	mark_deleted, 0.1
hnswlib	yes	53	yes, 0.1
usearch	yes	73	yes, 0.6
faiss (IndexHNSWFlat)	yes	129	none

nmslib and annoy are not in the table because they are batch indexes and say so — nmslib: "only static data sets are supported"; annoy: "after calling build, no more items can be added". A point added later is kept, but the structure that answers queries was built once and does not know about it, so getting it in means rebuilding the whole index. That is a different operation from an incremental insert, not a slower one, and there is nothing here to compare it against.

Deletion

mark_deleted(i) drops node i from results; unmark_deleted(i) puts it back, and deleted_count says how many are out. Deletion is soft, as in hnswlib: the node stops appearing in results but stays a routing waypoint, so the graph never gets disconnected and recall for everything else is untouched. The cost is that deleted nodes keep their memory and still cost distance computations while routing, so an index that is mostly deleted wants a rebuild.

index.mark_deleted(7)
ids, dists = index.search_batch(Q, k=10)   # 7 will not appear

The exclusion happens inside the search, not over its output: deleted nodes are still traversed, so what comes back is the true nearest live neighbors — not the unfiltered answer with rows struck out (measured >0.99 recall against brute force over the live set, with half the index deleted). You get min(k, ef, live) results: the search fills its pool with ef live candidates before it is willing to stop, so a short list means there really are that few live nodes, not that the search gave up. (k > ef truncates to ef, deletions or not — that is the pool size, so keep ef comfortably above k as usual.)

Persistence

anydist.save and anydist.load take a path or a binary file object:

import anydist

anydist.save(index, "index.pkl")
index = anydist.load("index.pkl")          # same graph, same results

There is no saves()/loads() pair returning bytes — pass an io.BytesIO(), as with numpy.save. Plain pickle (or dill) still works on an index too; save/load exist for the one thing pickle cannot do on its own.

The metric is the only part that may not travel. A built-in name or a module-level function pickles; a lambda, a closure or a C function pointer does not, and an address would be meaningless in the loading process anyway. Everything expensive — the vectors, the adjacency, the deletions — always serialises, so rebuilding would be a poor answer. Instead, save without the metric and give it back on load:

anydist.save(index, "index.pkl", include_metric=False)
index = anydist.load("index.pkl", metric=my_cfunc)

Including the metric is the default and raises if it cannot be done, so save never quietly writes a file that load(f) cannot open. You get a self-contained file or an exception naming the way out — never a surprise in the process that loads it.

Speed. Saving and loading are close to memory-bandwidth bound, so an index round-trips far faster than it builds. Measured on the machine described under About these numbers, writing to a real file in a temporary directory — as every library in the comparison does, so the figures are like-for-like:

index	on disk	save	load
Fashion-MNIST, 60k x 784	213 MB	0.18 s	0.17 s
Gaussian, 100k x 128	92 MB	0.09 s	0.06 s
Gaussian, 50k x 32	14 MB	0.02 s	0.01 s

Roughly 1.2 GB/s each way — the same order as a memory copy, and two orders of magnitude cheaper than rebuilding the index, which is the comparison that matters. (Nothing is fsynced, so this is the cost of handing the bytes to the OS rather than of a durable write.)

hnswlib is about twice as quick at both, and its file is ~8% smaller (197 MB against 213 MB on Fashion-MNIST): it writes a fixed-size record per element straight out of its own arena, where this library pickles, and it stores no edge weights. Ours keeps them whenever select does — under select='distances', which stores none, the sizes match. Both are far below the cost of building, so the difference is unlikely to be what you optimise.

Only load files you trust. save/load are built on pickle, and unpickling any file runs code that the file chooses. That is pickle's design, not something this library adds or can screen for: load sees the constructor arguments only after the interpreter has already executed whatever the stream asked for. Treat an index file exactly as you would a .py you were about to run — the format is not a safe interchange format for data from elsewhere.

load is strict in both directions: it refuses a metric= for a file that already has one (that would replace a metric known to be right with one it cannot check), and refuses a metric of the wrong kind — the three forms store their elements differently, so a callback cannot be attached to a graph built for a built-in name. What it cannot check is that you passed the same distance; supply a different one and the graph answers for a metric it was not built for. That much is yours to get right.

Under the hood this is a plain pickle using persistent_id, the stdlib's own mechanism for an object the loader must re-supply — there is no bespoke format here to version. The pickle carries the vectors (or the elements), the adjacency, the construction parameters and any soft deletions, so it is self-contained but not small: expect roughly the in-memory size of the index. There is no separate index-only/mmap format.

The two classes

They take the same constructor arguments and answer the same calls, and both accept all three metric forms. What only HNSWFast has is parallelism (update across cores, search_batch across queries), numpy arrays out of the batch API, and the float32 vector store.

HNSWPure accepts num_threads and inserts serially, and rejects float32=True having no typed store to narrow. It accepts keep_workspace=True and ignores it, rather than rejecting it: float32 changes the numbers, so ignoring a request for it would mislead about results, while keep_workspace only decides when memory is allocated and a class with no threads runs the same program either way. It does call C metrics, through ctypes, marshalling both operands on every distance. It is there so the reference implementation can run the metric you are actually indexing with — but that marshalling is a fixed per-call cost, so whether it beats writing the distance in Python depends on how much work the distance itself does. Against the plain Python loop a jit_metric wraps, the pointer measured ~2.8x slower at dim=8 and ~5x faster at dim=64. If you just want a working index at low dimension, that original two-argument function is still reachable on .py_func_2arg.

The two produce equivalent results only distributionally (ties are broken differently), so compare them via k-NN recall rather than exact neighbor lists.

Performance

On the dense-vector case this is in the same class as the established C++ libraries — which is the point, since that case is the one they optimize.

What the generality is and is not. It is not that nobody else indexes non-vector data: nmslib ships jaccard_sparse, leven and others, and is a real competitor on anything in its catalogue. The difference is that a catalogue is fixed at compile time and a Python callable is not. If your dissimilarity is in nmslib's list, use nmslib and expect it to win; if it is a business rule, a weighted ensemble, or anything you would otherwise have to write C++ for, this library will run it — and, with jit_metric, run it at a speed the section below quantifies.

Against the other libraries

50k Gaussian vectors in 32 dimensions, m=16, efConstruction=200, k=10, 12 threads. Each figure is a ratio to this library, formed within a repetition and shown with its range over three of them. A ratio whose range spans 1.0 is not a difference; the resolution limit here is 1.06x on build and 1.07x on query, measured by running the identical build twice as a control arm.

Both columns are times, relative to this library, so in both of them lower is faster and 1.00x is us:

library	build time	query time	recall@ef160
this library	1.00x	1.00x	.976
hnswlib 0.8.0	1.03x [0.93-1.10]	0.95x [0.94-1.06]	.980
faiss IndexHNSWFlat	1.06x [1.00-1.07]	0.89x [0.88-0.93]	.979
nmslib 2.1.2 (AVX2)	1.44x [1.42-1.47]	1.09x [1.06-1.12]	.980
usearch	1.36x [1.26-1.41]	1.54x [1.49-1.54]	.978
annoy (trees, not HNSW)	0.17x [0.16-0.18]	100x	1.000

hnswlib is indistinguishable from this library on both, which is the result that matters: the reference implementation of the algorithm, and we are inside the noise of it. faiss queries ~11% faster. We build ~1.4x faster than both nmslib and usearch, and query ~1.5x faster than usearch; nmslib's query time is within a whisker of ours. annoy is a tree index included for scale, and is the only arm whose parameters were not swept the way the others' ef was.

Recall is identical to three decimals across every HNSW arm at matched ef, which is the check that the comparison is like-for-like: the same graph quality, differing only in how fast it is reached.

Two caveats. nmslib's PyPI wheel is built without SIMD, and warns so on import; the row above is from a source build (pip install --no-binary :all: --no-build-isolation nmslib — plain --no-binary fails on Python 3.14). And the larger datasets we measured, at 100k x 128 and Fashion-MNIST, had control spreads of 2.0x to 5.6x on this machine, so they are omitted rather than reported: nothing there was resolvable.

An arbitrary metric, measured

kosarak from ann-benchmarks: 74,962 sets of ints, mean size 55.6, Jaccard, recall@10 against the published ground truth, M=16 efC=100. The comparison is nmslib's jaccard_sparse, built from source so it has AVX2 — its purpose-built C++ space against a Python function of ours.

Three ways to give this library the same distance. Every ratio below is formed inside one run, from the minimum of seven alternating repetitions, with a control arm — the same index timed a second time — so the resolution is measured rather than hoped for. The two columns of ratios span ef=100 and ef=400, that is recall 0.90 and 0.95:

how	against @jit_metric	against nmslib	resolution
plain Python callable	7.5–8.1x slower	—	1.02x
@jit_metric, 1 thread (pinned)	—	1.18–1.20x slower	1.002x
@jit_metric, 12 threads	—	1.24–1.26x slower	1.02–1.04x

So compiling the callable is worth about 8x, and once compiled the gap to a purpose-built C++ space is about 1.2x, whether serial or threaded. That is what writing your distance in Python instead of C++ costs.

Each arm queries through whichever API is faster for it: search_batch for the compiled metric, a search loop for the Python callable, which serially on this GIL build beat batching by ~6%. Matching them would widen the 8x slightly, not narrow it.

For scale rather than for comparison, the serial pinned run measured 0.61 ms/query for @jit_metric against nmslib's 0.50 at recall 0.90, and 0.095 against 0.075 at twelve threads. Absolute figures move between runs on this machine — the same compiled index measured 0.61 in one run and 0.73 in another, after a long build had heated the box — which is exactly why the ratios above are formed within a run and the absolutes are not compared across them.

The serial rows are pinned to one core with taskset: this machine's cores are not all the same speed, so an unpinned "1 thread" figure silently depends on which one it landed on — measured 2.6x apart for the same code. See About these numbers.

Variable-length objects reach a C metric through a fixed double* signature by storing [offset, length] in the row and closing over one flat array, so the vector store is 1.2 MB rather than the ~1.5 GB a padded matrix would need, and the distance stays exact:

FLAT = np.concatenate(sorted_sets)      # every set, back to back

@jit_metric
def jaccard(a, b):                      # each row is [offset, length]
    ao, an = int(a[0]), int(a[1])
    bo, bn = int(b[0]), int(b[1])
    i = j = inter = 0
    while i < an and j < bn:            # sorted-merge intersection
        x, y = FLAT[ao + i], FLAT[bo + j]
        if x == y:   inter += 1; i += 1; j += 1
        elif x < y:  i += 1
        else:        j += 1
    return 1.0 - inter / (an + bn - inter)

When an index is worth building at all. On an expensive metric the answer depends on n, because brute force is linear and the graph is not. Jaccard over random sets, ef=400, serial:

n	brute force	this library	speedup
5,000	5.04 ms/q	4.68 ms/q	1.1x
20,000	22.25 ms/q	9.50 ms/q	2.3x
80,000	87.72 ms/q	12.81 ms/q	6.8x

Below ~10k elements an index barely pays for itself; by ~100k it is most of the cost. (Recall drifts 0.955 -> 0.825 across those rows at fixed ef, so part of the last speedup is bought with quality.)

About these numbers

Every figure in this README comes from one throttling laptop, so the method matters more than usual. What was done, and what it means for how much weight a number can carry:

• A control arm. Every comparison includes the same build measured twice. The spread between those two runs is the resolution limit, and no effect smaller than it is claimed. This is not a formality: it caught a 10x "difference" that was two benchmark processes accidentally running at once, and a 3.8x one that was a metric interpolating across the steep part of a recall curve.
• Interleaved, paired repetitions. Arms alternate within a repetition and ratios are formed inside it, so drift in the machine is common-mode and cancels. Summarising each arm separately and dividing would not.
• One process per measurement, with the thread count pinned before any import and the observed value recorded, plus load average, CPU-seconds and timestamps on every row.
• Pinning, where it decides anything. The cores here are not all the same speed, so an unpinned single-threaded figure depends on which one it landed on — measured 2.6x apart for identical code. Serial comparisons use taskset; thread-count speedups are reported as measured, without being turned into an efficiency against a core count that would not mean much.

What this buys is roughly 3-13% resolution single-threaded and 6-7% on the 50k x 32 dataset at twelve threads. On the larger datasets at twelve threads it buys nothing: identical builds differed by 2.0x to 5.6x, so those cells are reported as unresolvable rather than given numbers. Redoing all of this on a quiet, homogeneous machine would tighten every bound here, and is on the list.

Reproducing any of this. Everything above is produced by benchmarks/ against the anydist you have installed — python run.py then python analyze.py, with every comparison library optional. The harness prints your machine's resolution before it prints any comparison, so you can see which of its own numbers to believe.

Author

Matteo Dell'Amico — della@linux.it

License

BSD 3-clause; see the LICENSE file.