3.3 Embeddings and Vector Search

1. From Words to Vectors: The Distributional Hypothesis

A researcher types "compounds targeting the SARS-CoV-2 main protease" into a search engine. Within milliseconds, it returns papers about molecules she has never heard of, phrased in terminology she never used. Somewhere in a 384-dimensional space, her query and those papers point in nearly the same direction. The reason is the distributional hypothesis (Harris, 1954; Firth, 1957): words that occur in similar contexts tend to have similar meanings. "You shall know a word by the company it keeps," as Firth put it. This insight underpins every word embedding method. The lineage runs from early Latent Semantic Analysis (LSA) through Word2Vec (Mikolov et al., 2013) to modern transformer-based encoders like Bidirectional Encoder Representations from Transformers (BERT) and Sentence-BERT.

A word embedding maps each word \(w\) in a vocabulary \(V\) to a dense vector \(\mathbf{w} \in \mathbb{R}^d\), where \(d\) is typically 100 to 768 dimensions. The mapping is learned so that semantically related words are geometrically close. The celebrated vector arithmetic property of Word2Vec demonstrates this: \(\mathbf{v}(\text{king}) - \mathbf{v}(\text{man}) + \mathbf{v}(\text{woman}) \approx \mathbf{v}(\text{queen})\).

What an Embedding Captures

A word embedding is a learned lookup table. It assigns every word in a vocabulary a fixed-length vector, positioning it in a continuous geometric space where proximity encodes semantic relatedness. This matters because it converts discrete symbols into points. Computers cannot compare symbols for meaning, but they can compute dot products, distances, and averages on vectors, directly measuring and manipulating meaning. A shallow neural network learns the mapping by predicting a word from its surrounding context (or vice versa). Words appearing in similar neighborhoods converge to nearby vectors. Use embeddings whenever you need to compute semantic similarity, cluster concepts, or feed text into a downstream model. Use symbolic representations (ontologies, knowledge graphs) when you need explicit, human-auditable logical relationships.

For scientific discovery, we need embeddings that capture domain-specific semantics. "Cell" should be close to "neuron" in a biology context but close to "battery" in a materials science context. Contextual embeddings (from BERT and its descendants) handle this naturally: the same word gets different vectors depending on its surrounding text. Sentence embeddings extend this to entire passages: a research abstract becomes a single vector that captures its semantic content. In short: embeddings translate meaning into geometry, so that finding related concepts becomes finding nearby points.

from sentence_transformers import SentenceTransformer
import numpy as np

# Load a pre-trained sentence embedding model
model = SentenceTransformer("all-MiniLM-L6-v2")  # 384 dimensions

# Encode scientific abstracts as dense vectors
abstracts = [
    "CRISPR-Cas9 enables precise genome editing in human cells.",
    "Gene editing with CRISPR technology allows targeted DNA modifications.",
    "Convolutional neural networks achieve state-of-the-art image classification.",
    "Deep learning models for computer vision use hierarchical feature extraction.",
    "Graphene exhibits exceptional electrical conductivity and tensile strength.",
]

embeddings = model.encode(abstracts)
print(f"Embedding shape: {embeddings.shape}")  # (5, 384)

# Compute pairwise cosine similarities
def cosine_sim(a, b):
    return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b))

print("\nPairwise cosine similarities:")
labels = ["CRISPR-1", "CRISPR-2", "CNN", "DeepLearn", "Graphene"]
for i in range(len(abstracts)):
    for j in range(i + 1, len(abstracts)):
        sim = cosine_sim(embeddings[i], embeddings[j])
        print(f"  {labels[i]:10s} <-> {labels[j]:10s}: {sim:.3f}")

# Output:
# Embedding shape: (5, 384)
#
# Pairwise cosine similarities:
#   CRISPR-1   <-> CRISPR-2  : 0.874
#   CRISPR-1   <-> CNN       : 0.128
#   CRISPR-1   <-> DeepLearn : 0.103
#   CRISPR-1   <-> Graphene  : 0.091
#   CRISPR-2   <-> CNN       : 0.115
#   CRISPR-2   <-> DeepLearn : 0.098
#   CRISPR-2   <-> Graphene  : 0.107
#   CNN        <-> DeepLearn : 0.812
#   CNN        <-> Graphene  : 0.142
#   DeepLearn  <-> Graphene  : 0.119

2. Similarity Metrics for Vector Spaces

Cosine similarity is the most common metric for text embeddings, but not the only one. The choice of metric depends on how the embeddings were trained and what notion of "closeness" matters for your application:

• Cosine similarity / angular distance: measures directional alignment, invariant to vector magnitude. Preferred for text embeddings where the direction encodes meaning and the norm is an artifact of encoding.
• Euclidean distance (\(L_2\)): \(d(\mathbf{a}, \mathbf{b}) = \|\mathbf{a} - \mathbf{b}\|_2\). Sensitive to magnitude. Used by TransE (Section 3.2) and when absolute position in the space carries meaning.
• Inner product (dot product): \(\mathbf{a} \cdot \mathbf{b} = \sum_i a_i b_i\). Equivalent to cosine similarity for unit-normalized vectors. Used in maximum inner product search (MIPS), where MIPS is the task of finding the vector with the highest dot product against a query, which arises in recommendation systems and attention mechanisms.
• Manhattan distance (\(L_1\)): \(d(\mathbf{a}, \mathbf{b}) = \sum_i |a_i - b_i|\). More robust to outlier dimensions; sometimes used in sparse or binary embeddings.

For the hybrid search system we build in Section 3.4, we use cosine similarity for the dense component and BM25 (Best Matching 25, a statistical score, not a vector metric) for the sparse component.

3. The Nearest-Neighbor Problem

Choosing a similarity metric tells us how to score any pair of vectors, but it says nothing about how to find the best-scoring partner quickly when the candidate pool contains millions of entries.

Given a dataset of \(n\) vectors in \(d\) dimensions and a query vector \(\mathbf{q}\), the nearest-neighbor problem asks: find the vector \(\mathbf{x}^*\) in the dataset that maximizes \(\text{sim}(\mathbf{q}, \mathbf{x}^*)\). The brute-force approach computes all \(n\) similarities and takes \(O(nd)\) time. For a million 768-dimensional vectors, that is 768 million floating-point operations per query, roughly 2 milliseconds on a modern CPU. Acceptable for a single query, but not for serving thousands of concurrent users or indexing billions of vectors.

Three families of approximate nearest-neighbor (ANN) algorithms trade a small accuracy loss for dramatic speed gains:

• Tree-based: random projection trees (Annoy), KD-trees. Split the space with hyperplanes. Fast for low dimensions (\(d < 20\)), but performance degrades in high dimensions due to the curse of dimensionality, where the curse of dimensionality refers to the phenomenon that in high-dimensional spaces, distances between points become nearly uniform, so tree splits lose their discriminative power.
• Hash-based: locality-sensitive hashing (LSH). Map similar vectors to the same hash bucket. Theoretically elegant with provable guarantees, but often outperformed in practice by graph-based methods.
• Graph-based: navigable small-world graphs (NSW), HNSW. Build a proximity graph where each node connects to its approximate nearest neighbors, then search by greedy traversal. HNSW is the dominant algorithm in modern vector databases.

4. HNSW: Hierarchical Navigable Small World Graphs

When Semantic Scholar serves a query against 200 million papers, brute-force comparison is out of the question: scanning every vector would take minutes, not milliseconds. The algorithm that makes sub-millisecond retrieval possible at that scale is HNSW.

HNSW (Malkov and Yashunin, 2020) builds a multi-layer graph where each layer is a navigable small-world graph with exponentially fewer nodes. Figure 3.5 illustrates this hierarchical structure. The construction proceeds as follows:

1. Layer assignment: each new element is assigned a maximum layer \(\ell = \lfloor -\ln(\text{uniform}(0,1)) \cdot m_L \rfloor\), where \(m_L = 1/\ln(M)\) and \(M\) is the maximum number of connections per node. This gives an exponential distribution: most elements appear only in layer 0 (the bottom), and progressively fewer elements appear in higher layers.
2. Insertion: starting from an entry point at the top layer, greedily descend to layer \(\ell + 1\) by following the nearest neighbor at each layer. Then, at each layer from \(\ell\) down to 0, find the \(\text{efConstruction}\) nearest neighbors (efConstruction is a build-time parameter that controls how many candidates the algorithm evaluates when wiring a new node; higher values produce a better-connected graph at the cost of slower construction) and connect the new element to the closest \(M\) of them (bidirectionally).
3. Search: start at the entry point, greedily descend to layer 0, then perform a beam search at layer 0, where beam search is a bounded breadth-first traversal that maintains the \(\text{ef}\) best candidates seen so far, to find the \(k\) nearest neighbors.

The hierarchical structure provides logarithmic search complexity: each layer roughly halves the search space (like a skip list, a data structure with multiple linked-list levels where higher levels skip over elements for faster traversal), and the bottom layer provides fine-grained resolution. The result is \(O(\log n)\) search time with high recall, compared to \(O(n)\) for brute force. Figure 3.3.1 illustrates HNSW multi-layer graph search traversal.

Mental Model

Think of HNSW like navigating a city using progressively more detailed maps. The top layer is a highway map showing only major cities: you jump quickly across the country in a few hops, but you cannot find a specific street. The middle layer is a regional map with towns and neighborhoods. The bottom layer is a full street map with every address. When searching for a destination, you first pick the nearest major city on the highway map, then zoom in to the right neighborhood, then walk street by street to the exact address. Each "zoom level" has fewer, longer-range connections (highways) or more, shorter-range connections (local streets). This is why HNSW is fast: instead of checking every street in the country, you narrow from continent to city to block to building, visiting only a logarithmic fraction of the total addresses.

The following code builds a simplified HNSW index from scratch:

import numpy as np
import heapq
from collections import defaultdict

class SimpleHNSW:
    """A simplified HNSW index for educational purposes."""

    def __init__(self, dim, M=16, ef_construction=200, m_L=None):
        self.dim = dim
        self.M = M                    # max connections per node per layer
        self.M0 = 2 * M               # max connections at layer 0
        self.ef_construction = ef_construction
        self.m_L = m_L or 1.0 / np.log(M)
        self.vectors = []             # list of vectors
        self.graphs = defaultdict(     # layer -> {node_id: set of neighbor ids}
            lambda: defaultdict(set)
        )
        self.max_layer = -1
        self.entry_point = None

    def _distance(self, a, b):
        """Euclidean distance between two vectors."""
        return np.linalg.norm(self.vectors[a] - self.vectors[b])

    def _distance_to_query(self, query, b):
        """Euclidean distance from query vector to stored vector b."""
        return np.linalg.norm(query - self.vectors[b])

    def _random_level(self):
        """Assign a random layer using exponential distribution."""
        return int(-np.log(np.random.uniform()) * self.m_L)

    def _search_layer(self, query, entry_points, ef, layer):
        """Greedy search within a single layer, returning ef nearest neighbors."""
        visited = set(entry_points)
        # Min-heap of candidates (distance, id): closest first
        candidates = []
        # Max-heap of results (neg_distance, id): farthest-of-best first
        results = []

        for ep in entry_points:
            dist = self._distance_to_query(query, ep)
            heapq.heappush(candidates, (dist, ep))
            heapq.heappush(results, (-dist, ep))

        while candidates:
            d_c, c = heapq.heappop(candidates)
            # Farthest element in results
            d_f = -results[0][0]
            if d_c > d_f:
                break  # all remaining candidates are farther than worst result

            for neighbor in self.graphs[layer][c]:
                if neighbor not in visited:
                    visited.add(neighbor)
                    d_n = self._distance_to_query(query, neighbor)
                    if len(results) < ef or d_n < -results[0][0]:
                        heapq.heappush(candidates, (d_n, neighbor))
                        heapq.heappush(results, (-d_n, neighbor))
                        if len(results) > ef:
                            heapq.heappop(results)

        return [(node_id, -neg_dist) for neg_dist, node_id in sorted(results)]

    def insert(self, vector):
        """Insert a vector into the index."""
        idx = len(self.vectors)
        self.vectors.append(vector)
        level = self._random_level()

        if self.entry_point is None:
            self.entry_point = idx
            self.max_layer = level
            return idx

        # Phase 1: greedily descend from top to level+1
        curr_ep = [self.entry_point]
        for lc in range(self.max_layer, level, -1):
            nearest = self._search_layer(vector, curr_ep, ef=1, layer=lc)
            curr_ep = [nearest[0][0]]

        # Phase 2: insert at each layer from min(level, max_layer) down to 0
        for lc in range(min(level, self.max_layer), -1, -1):
            neighbors = self._search_layer(
                vector, curr_ep, ef=self.ef_construction, layer=lc
            )
            max_conn = self.M0 if lc == 0 else self.M
            selected = [n_id for n_id, _ in neighbors[:max_conn]]

            for n_id in selected:
                self.graphs[lc][idx].add(n_id)
                self.graphs[lc][n_id].add(idx)
                # Prune if too many connections
                if len(self.graphs[lc][n_id]) > max_conn:
                    # Keep only the closest max_conn neighbors
                    scored = [(self._distance(n_id, nb), nb)
                              for nb in self.graphs[lc][n_id]]
                    scored.sort()
                    self.graphs[lc][n_id] = set(nb for _, nb in scored[:max_conn])

            curr_ep = [n_id for n_id, _ in neighbors]

        if level > self.max_layer:
            self.max_layer = level
            self.entry_point = idx

        return idx

    def search(self, query, k=5, ef=50):
        """Find k approximate nearest neighbors of query."""
        if self.entry_point is None:
            return []
        curr_ep = [self.entry_point]
        for lc in range(self.max_layer, 0, -1):
            nearest = self._search_layer(query, curr_ep, ef=1, layer=lc)
            curr_ep = [nearest[0][0]]
        # Search layer 0 with larger beam width
        results = self._search_layer(query, curr_ep, ef=max(ef, k), layer=0)
        return results[:k]


# Build index with 1000 random vectors
np.random.seed(42)
dim = 128
n_vectors = 1000
data = np.random.randn(n_vectors, dim).astype(np.float32)
# Normalize to unit vectors (so Euclidean distance ~ angular distance)
data /= np.linalg.norm(data, axis=1, keepdims=True)

hnsw = SimpleHNSW(dim=dim, M=16, ef_construction=100)
for vec in data:
    hnsw.insert(vec)

# Query
query = np.random.randn(dim).astype(np.float32)
query /= np.linalg.norm(query)

# HNSW search
hnsw_results = hnsw.search(query, k=5, ef=50)

# Brute-force search for comparison
dists = np.linalg.norm(data - query, axis=1)
bf_top5 = np.argsort(dists)[:5]

print("HNSW top-5 (id, distance):")
for node_id, dist in hnsw_results:
    print(f"  id={node_id:4d}, dist={dist:.4f}")

print("\nBrute-force top-5 (id, distance):")
for idx in bf_top5:
    print(f"  id={idx:4d}, dist={dists[idx]:.4f}")

# Check recall
hnsw_ids = set(nid for nid, _ in hnsw_results)
bf_ids = set(bf_top5.tolist())
recall = len(hnsw_ids & bf_ids) / len(bf_ids)
print(f"\nRecall@5: {recall:.1%}")
# Output (typical): Recall@5: 100.0%

print(f"Layers used: {hnsw.max_layer + 1}")
print(f"Layer 0 avg connections: "
      f"{np.mean([len(v) for v in hnsw.graphs[0].values()]):.1f}")
# Output (typical):
# Layers used: 3
# Layer 0 avg connections: 12.4

Fun Note: Why "Small World"?

The "small world" in HNSW refers to the small-world network property discovered by Watts and Strogatz (1998): most nodes can be reached from any other node in a small number of hops, even in very large networks. Stanley Milgram's famous "six degrees of separation" experiment (1967) suggested this pattern in social networks. HNSW exploits the same property: by adding a few long-range connections (in higher layers), the graph becomes navigable in \(O(\log n)\) hops. The analogy is precise: the upper layers of HNSW act like the "weak ties" in social networks, the acquaintances who connect distant communities and let messages traverse the globe in six steps.

import faiss

dim, n = 128, 1_000_000
data = np.random.randn(n, dim).astype("float32")
faiss.normalize_L2(data)  # unit normalize

index = faiss.IndexHNSWFlat(dim, 32)  # 32 connections per node
index.hnsw.efConstruction = 200
index.add(data)

index.hnsw.efSearch = 64
distances, indices = index.search(query.reshape(1, -1), k=10)
# Searches 1M vectors in < 1ms

5. Vector Databases: Beyond In-Memory Indices

A vector database wraps an ANN index with the features needed for production systems: persistence, replication, filtering, multi-tenancy, and CRUD operations (create, read, update, delete). The key players as of 2025:

• Qdrant: Rust-based, supports payload filtering (metadata predicates applied during vector search), quantization, and distributed deployment. Its filtering is particularly useful for scientific search: "find papers similar to this query and published after 2020 and tagged with 'proteomics'."
• pgvector: a PostgreSQL extension that adds vector columns, cosine/Euclidean/inner-product operators, and IVFFlat (Inverted File with Flat quantization, a partitioning index that clusters vectors into Voronoi cells for faster search) or HNSW indices. Ideal when your vectors live alongside relational data (authors, institutions, citation counts) in a PostgreSQL database you already operate.
• Pinecone: fully managed cloud service with serverless scaling. No infrastructure to manage, but less control over index parameters.
• Weaviate: combines vector search with a built-in object model and module system for embedding generation.
• Milvus / Zilliz: open-source, designed for billion-scale vector search with GPU acceleration.

For the Discovery Workbench, we use Qdrant for its payload-filtering capability (essential for faceted scientific search) and pgvector for scenarios where vectors must join with relational metadata in SQL queries. Both appear in the following examples:

from qdrant_client import QdrantClient
from qdrant_client.models import Distance, VectorParams, PointStruct, Filter, FieldCondition, Range

# Connect to Qdrant (in-memory mode for demo)
client = QdrantClient(":memory:")

# Create a collection for scientific papers
client.create_collection(
    collection_name="papers",
    vectors_config=VectorParams(size=384, distance=Distance.COSINE),
)

# Insert paper embeddings with metadata payloads
papers = [
    {"title": "CRISPR-Cas9 genome editing in human cells",
     "year": 2023, "field": "biology", "citations": 142},
    {"title": "Transformer architectures for protein folding",
     "year": 2024, "field": "biology", "citations": 89},
    {"title": "Graph neural networks for molecular property prediction",
     "year": 2023, "field": "chemistry", "citations": 67},
    {"title": "Large language models for scientific reasoning",
     "year": 2024, "field": "AI", "citations": 203},
]

# Generate dummy embeddings (in practice, use Sentence-BERT)
np.random.seed(42)
dummy_embeddings = np.random.randn(len(papers), 384).astype(np.float32)
for i, emb in enumerate(dummy_embeddings):
    dummy_embeddings[i] = emb / np.linalg.norm(emb)

points = [
    PointStruct(id=i, vector=emb.tolist(), payload=meta)
    for i, (emb, meta) in enumerate(zip(dummy_embeddings, papers))
]
client.upsert(collection_name="papers", points=points)

# Search: find similar papers, filtered to biology and year >= 2024
query_vector = dummy_embeddings[0].tolist()  # use first paper as query
results = client.search(
    collection_name="papers",
    query_vector=query_vector,
    query_filter=Filter(
        must=[
            FieldCondition(key="field", match={"value": "biology"}),
            FieldCondition(key="year", range=Range(gte=2024)),
        ]
    ),
    limit=3,
)

print("Filtered search results (biology, year >= 2024):")
for result in results:
    print(f"  score={result.score:.3f} | {result.payload['title']} "
          f"({result.payload['year']})")
# Output:
# Filtered search results (biology, year >= 2024):
#   score=0.312 | Transformer architectures for protein folding (2024)

6. Quantization: Compressing Vectors for Scale

Vector databases solve the operational challenges of persistence, filtering, and multi-tenancy, but they still store every dimension of every vector at full precision, and that cost compounds fast as collections grow from millions to billions of entries.

At scale, memory becomes a bottleneck. One million 768-dimensional float32 vectors consume roughly 3 GB of RAM. A billion vectors would require 3 TB, well beyond typical server memory. (For perspective, that single index would consume more RAM than most database clusters use for all their tables combined.) Vector quantization compresses vectors while preserving approximate distances:

• Scalar quantization: reduce each float32 component to int8, cutting memory by 4x with minimal recall loss. Each dimension is independently mapped from its observed range to [0, 255].
• Product quantization (PQ): split each vector into \(m\) subvectors, quantize each subvector independently using a codebook (a fixed set of representative vectors, called centroids, learned by clustering) of \(k\) centroids. A 768-dimensional vector with \(m = 96\) subgroups and \(k = 256\) centroids compresses from 3,072 bytes to 96 bytes (32x compression). PQ enables asymmetric distance computation: the query remains uncompressed, and distances are computed via lookup tables.
• Binary quantization: reduce each dimension to a single bit (sign of the component). A 768-dimensional vector becomes 96 bytes. Recall drops significantly, but Hamming distance search is extremely fast, making binary quantization useful as a first-pass filter before rescoring with full-precision vectors.

In practice, most vector databases apply quantization as an optional layer: Qdrant supports scalar and product quantization; FAISS provides IVF-PQ and OPQ (optimized product quantization); pgvector supports half-precision vectors. The choice depends on your recall requirements and memory budget.