Embeddings in NumPy Before You Add a Database

A vector database is a real answer to a real problem, and most projects adopt one about two orders of magnitude before they have that problem. Below roughly a hundred thousand vectors, a NumPy array and a matrix multiply are simpler, exact, and fast enough — and the arithmetic that shows why fits on one screen.

The arithmetic, first

Two quantities decide whether brute force is viable: how much memory the vectors occupy, and how much arithmetic one query costs. Both are exactly calculable. Take 100,000 vectors of 1,536 dimensions, which is a common embedding size.

MEMORY
  values          = 100,000 x 1,536              = 153,600,000 floats
  float32         = 153,600,000 x 4 bytes        = 614,400,000 B  ~ 586 MiB
  float16         = 153,600,000 x 2 bytes        = 307,200,000 B  ~ 293 MiB
  int8            = 153,600,000 x 1 byte         = 153,600,000 B  ~ 147 MiB

ARITHMETIC PER QUERY (one dot product against every vector)
  multiply-adds   = 100,000 x 1,536              = 153,600,000 MACs
  FLOPs           = 2 x 153,600,000              = 307,200,000    ~ 0.31 GFLOP

MEMORY TRAFFIC PER QUERY
  every stored value is read exactly once        = 586 MiB (float32)

The third block is the one that predicts the wall-clock time, because a single dot-product pass is memory-bound, not compute-bound: 0.31 GFLOP is a rounding error for any modern CPU, while reading 614 MB is not. If your machine sustains B gigabytes per second of memory bandwidth, one query costs roughly 0.614 / B seconds before overheads. That is the whole model, and it is why the honest thing to do next is measure B on your machine rather than trust a number from somebody else’s.

Three consequences fall out of the same arithmetic. Storing float16 halves the traffic and therefore roughly halves the query time. The cost is linear in the number of vectors, so 1,000,000 vectors is ten times the time and 5.7 GiB of RAM — which is where brute force stops being reasonable. And batching queries is nearly free: a hundred queries as one matrix multiply reads the corpus once instead of a hundred times, which is the single largest optimisation available here.

Storing the vectors

One .npy file for the matrix and one JSONL file for the metadata, with row index as the join key. This is unglamorous and it survives restarts, version control and being handed to somebody else.

# build_index.py
import json

import numpy as np


def build(records: list[dict], embed: callable, path: str) -> None:
    """records: [{"id": ..., "text": ...}]. embed: list[str] -> list[list[float]]"""
    vectors: list[list[float]] = []
    with open(f"{path}.jsonl", "w", encoding="utf-8") as meta:
        for start in range(0, len(records), 256):        # batch the API calls
            batch = records[start:start + 256]
            vectors.extend(embed([r["text"] for r in batch]))
            for record in batch:
                meta.write(json.dumps(record, ensure_ascii=False) + "\n")

    matrix = np.asarray(vectors, dtype=np.float32)
    matrix /= np.linalg.norm(matrix, axis=1, keepdims=True)   # normalise ONCE
    np.save(f"{path}.npy", matrix)
    print(f"{matrix.shape[0]} vectors x {matrix.shape[1]} dims,"
          f" {matrix.nbytes / 2**20:.1f} MiB")

Normalising at build time is the trick that makes the query cheap. Once every row has unit length, cosine similarity is the dot product — no division, no norms recomputed per query. Do it once for 100,000 rows instead of 100,000 times per query.

The embedding call itself is the same HTTP shape as everything else in this cluster:

def embed(texts: list[str]) -> list[list[float]]:
    response = client.post("/embeddings",
                           json={"model": EMBED_MODEL, "input": texts})
    response.raise_for_status()
    data = response.json()["data"]
    data.sort(key=lambda item: item["index"])   # do not assume input order
    return [item["embedding"] for item in data]

The sort is not paranoia — the response carries an index field precisely because the order is not guaranteed, and a misaligned batch produces an index where every answer is subtly wrong and nothing errors.

The search, in four lines

# search.py
import json

import numpy as np


class Index:
    def __init__(self, path: str):
        self.matrix = np.load(f"{path}.npy")             # (n, d), unit rows
        with open(f"{path}.jsonl", encoding="utf-8") as fh:
            self.meta = [json.loads(line) for line in fh]
        assert self.matrix.shape[0] == len(self.meta), "index and metadata differ"

    def search(self, query_vector: list[float], k: int = 5) -> list[dict]:
        q = np.asarray(query_vector, dtype=np.float32)
        q /= np.linalg.norm(q)
        scores = self.matrix @ q                          # (n,) cosine similarities
        top = np.argpartition(-scores, min(k, len(scores) - 1))[:k]
        top = top[np.argsort(-scores[top])]               # order the k, not the n
        return [{"score": float(scores[i]), **self.meta[i]} for i in top]

    def search_many(self, query_vectors: np.ndarray, k: int = 5) -> np.ndarray:
        """(m, d) queries -> (m, k) row indices. Reads the corpus once."""
        q = query_vectors / np.linalg.norm(query_vectors, axis=1, keepdims=True)
        scores = q @ self.matrix.T                        # (m, n)
        top = np.argpartition(-scores, k, axis=1)[:, :k]
        ordered = np.take_along_axis(
            top, np.argsort(-np.take_along_axis(scores, top, axis=1), axis=1), axis=1
        )
        return ordered

self.matrix @ q is the entire search. It is exact — there is no recall trade-off, unlike every approximate index — and it dispatches to your platform’s BLAS, which is heavily optimised C and not Python.

The assert on the second line of __init__ has saved more time than the rest of the class. The most common failure of a file-backed index is a matrix rebuilt without its metadata, which does not error — it returns confident results about the wrong documents.

Timing it on your own machine

The arithmetic above bounds the work; it cannot tell you the seconds, because that depends on your memory bandwidth, your BLAS build and whether the array fits in cache. Run this and you will have the number for your hardware, which is the only one worth having.

# bench.py
import time

import numpy as np

RUNS = 20


def bench(n: int, d: int, dtype=np.float32) -> None:
    rng = np.random.default_rng(0)
    matrix = rng.standard_normal((n, d)).astype(dtype)
    matrix /= np.linalg.norm(matrix, axis=1, keepdims=True)
    query = matrix[0].copy()

    matrix @ query                                   # warm up: page in, plan BLAS

    single = []
    for _ in range(RUNS):
        start = time.perf_counter()
        scores = matrix @ query
        np.argpartition(-scores, 5)[:5]
        single.append(time.perf_counter() - start)
    single.sort()

    batch = rng.standard_normal((100, d)).astype(dtype)
    batch /= np.linalg.norm(batch, axis=1, keepdims=True)
    start = time.perf_counter()
    batch @ matrix.T
    batched = time.perf_counter() - start

    gib = matrix.nbytes / 2**30
    median = single[len(single) // 2]
    print(f"n={n:,} d={d} dtype={np.dtype(dtype).name}")
    print(f"  resident      {matrix.nbytes / 2**20:8.1f} MiB")
    print(f"  1 query p50   {median * 1e3:8.2f} ms"
          f"   -> {gib / median:6.1f} GiB/s effective")
    print(f"  100 queries   {batched * 1e3:8.2f} ms"
          f"   ({batched / 100 * 1e3:.3f} ms each)")


if __name__ == "__main__":
    for n in (10_000, 100_000, 1_000_000):
        bench(n, 1536)

Three things to read out of the output. Whether single-query latency is acceptable for your interface at all. The effective GiB/s figure, which you can compare against your machine’s specification to see whether NumPy is reaching the hardware. And the ratio between the batched and single per-query times, which is how much you gain by batching — it is usually large, and it changes how you design the calling code.

Top-k without a full sort

np.argsort(-scores)[:k] sorts all 100,000 scores to keep five. That is O(n log n) work for an O(n) problem. np.argpartition does a partial selection instead: it puts the k best in the first k positions, in arbitrary order, in linear time. Sort only those k afterwards.

# wrong, and gets slower as the corpus grows
top = np.argsort(-scores)[:k]

# right: partition n, then sort k
top = np.argpartition(-scores, k)[:k]
top = top[np.argsort(-scores[top])]

Two edge cases the second form has: argpartition raises if k is not smaller than the array length, so clamp it, and the resulting indices are unordered until the second line runs, so do not use them before then.

A metadata filter is a boolean mask applied to the scores, which stays exact where an approximate index has to choose between filtering before or after the search:

mask = np.array([m["lang"] == "en" for m in self.meta])
scores = np.where(mask, scores, -np.inf)

Keeping the index current

Documents change, and re-embedding everything nightly is both slow and an unnecessary bill. Two operations cover almost every case, and both are ordinary array manipulation.

# update.py
import json

import numpy as np


def append(path: str, records: list[dict], embed) -> None:
    """Add new documents without touching the existing vectors."""
    matrix = np.load(f"{path}.npy")
    fresh = np.asarray(embed([r["text"] for r in records]), dtype=np.float32)
    fresh /= np.linalg.norm(fresh, axis=1, keepdims=True)
    np.save(f"{path}.npy", np.vstack([matrix, fresh]))
    with open(f"{path}.jsonl", "a", encoding="utf-8") as meta:
        for record in records:
            meta.write(json.dumps(record, ensure_ascii=False) + "\n")


def rewrite(path: str, keep: np.ndarray) -> None:
    """Drop rows. keep is a boolean mask over the current rows."""
    matrix = np.load(f"{path}.npy")
    with open(f"{path}.jsonl", encoding="utf-8") as fh:
        meta = [json.loads(line) for line in fh]
    np.save(f"{path}.npy", matrix[keep])
    with open(f"{path}.jsonl", "w", encoding="utf-8") as fh:
        for row, keeping in zip(meta, keep):
            if keeping:
                fh.write(json.dumps(row, ensure_ascii=False) + "\n")

Deletion has to rewrite both files together, and that is the operation to be careful with: an interrupted rewrite leaves the matrix and the metadata describing different sets of documents, which the assert in the loader catches only if the counts happen to differ. Write to temporary paths and rename both at the end, or accept that a crash means a full rebuild.

Store the embedding model id and the content hash alongside each record. The model id is what tells you which rows need re-embedding after a model change — vectors from two different models are not comparable and mixing them silently degrades every result, not just the new rows. The content hash is what lets an incremental job embed only the documents whose text actually changed, rather than everything with a new modification time.

When you do need the database

Brute force stops being the right answer at fairly specific boundaries, and it is worth knowing which one you have hit rather than migrating out of unease.

• The corpus no longer fits in RAM. This is the hard one. At 1,536 dimensions in float32 that is about 5.7 GiB per million vectors; once the array is swapping, every query pays disk latency and the arithmetic above stops applying.
• Latency at your size exceeds your budget. Measured with the harness above, not assumed. If a query must return in 10 ms and it takes 60, an approximate index such as HNSW trades a few per cent of recall for a large factor of speed.
• Writes are continuous. Rebuilding a .npy file is fine nightly and hopeless if documents arrive every second. Live insert and delete is what a real index gives you.
• Several processes must query it. One array per process means one copy of 586 MiB per process. That is the point at which a shared server pays for itself.

Until then, keeping the index as two files means it can be rebuilt from scratch in a script, committed to object storage, diffed and thrown away. Do you need a vector database argues the same case from the operational side.