Hierarchical Balanced Clustering (HBC)

Overview

Hierarchical Balanced Clustering (HBC) is a tree-based vector indexing strategy designed for approximate nearest neighbor (ANN) search over high-dimensional embedding spaces. Unlike graph-based approaches such as HNSW, which maintain navigable small-world graphs with per-node neighbor lists, HBC organizes vectors into a balanced hierarchical tree where each internal node represents a cluster centroid and each leaf holds a compact set of data vectors. This structure enables efficient top-down traversal during queries while keeping build times predictable through recursive partitioning.

The core insight behind HBC is that balanced partitioning of vector space — where every subtree contains roughly the same number of vectors — yields trees with logarithmic depth and uniform workload distribution across nodes. By choosing among different splitting strategies (even, left-biased, and right-biased splits), the algorithm adapts to skewed data distributions that would otherwise degrade search quality. This flexibility makes HBC particularly well-suited for datasets where cluster densities vary significantly across regions of the embedding space.

In practice, HBC serves as the foundational vector index in production search systems, where it provides the backbone for similarity search over document embeddings, image feature vectors, and other high-dimensional data. The bulk build algorithm constructs the tree in a single pass with O(n log n) average-case complexity, while beam search traversal at query time achieves sub-linear lookup with tunable recall–latency trade-offs.

Core Concepts

Hierarchical Tree Structures for ANN

At its foundation, HBC builds a rooted tree over a corpus of n vectors in d dimensions. Internal nodes store a centroid vector — the mean of all vectors in their subtree — along with metadata such as the subtree’s bounding radius. Leaf nodes contain small buckets of actual data vectors (typically 32–256 vectors per leaf). During a query, the search algorithm uses the centroid hierarchy to prune entire subtrees whose centroids are far from the query vector, dramatically reducing the number of distance computations required.

Splitting Strategies: Even, Left, and Right

When a node’s vector set exceeds the maximum leaf capacity, it must be split into child nodes. HBC supports three splitting strategies that control how vectors are partitioned across children:

• Even Split: Divides the vector set into two (or more) roughly equal-sized partitions. This is the default strategy and produces the most balanced trees. Vectors are sorted along the principal axis of variance, then divided at the median. This yields subtrees of near-equal size regardless of the data distribution.
• Left Split: Assigns a larger fraction of vectors to the left child (e.g., 60–70%) and a smaller fraction to the right. This is useful when the data distribution has a dense core region — the left child absorbs the dense cluster while the right child handles outliers, allowing tighter bounding radii on the right subtree for faster pruning.
• Right Split: The mirror of left split, assigning more vectors to the right child. This is applied when the optimizer detects that right-biased partitioning produces tighter centroid bounds, which can occur with certain skewed distributions where the densest region lies on the upper end of the principal axis.

Centroid Computation

Each internal node’s centroid is the arithmetic mean of all vectors in its subtree. During the bulk build, centroids are computed bottom-up: leaf centroids are calculated directly from their vector buckets, and internal centroids are derived as the weighted mean of their children’s centroids (weighted by child subtree size). This approach avoids redundant passes over the data and ensures O(n · d) total work for centroid computation across the entire tree.

In addition to the centroid, each node stores a bounding radius — the maximum distance from the centroid to any vector in the subtree. This radius is critical for the triangle-inequality pruning used during beam search: if the distance from the query to a node’s centroid minus the bounding radius exceeds the current best distance, the entire subtree can be safely skipped.

Balanced Partitioning

True balance in HBC means that the ratio of the largest leaf subtree to the smallest stays within a configurable bound (commonly 2:1 or 3:1). The partitioning algorithm enforces this by checking the split ratio after each candidate split and falling back to an even split if the imbalance exceeds the threshold. This guarantee bounds the tree depth to O(log n) and ensures that beam search explores a predictable number of nodes regardless of which region of vector space the query falls into.

How It Works

Bulk Build via Recursive Partitioning

The HBC bulk build algorithm constructs the index in a single top-down pass. Given an input array of n vectors, it proceeds recursively:

1. Base case: If the vector count is at or below the leaf capacity L, create a leaf node containing all vectors and compute its centroid and bounding radius.
2. Principal axis selection: Compute the variance along each dimension and select the dimension with the highest variance as the splitting axis. This maximizes the separation between the resulting partitions.
3. Candidate split evaluation: For each of the three split strategies (even, left, right), compute the resulting partition sizes, child centroids, and bounding radii. Select the strategy that minimizes the maximum child bounding radius.
4. Partition: Sort or partition the vectors along the chosen axis according to the selected split strategy. This is done in-place using a nth_element-style selection algorithm for O(n) partitioning.
5. Recurse: Build child nodes recursively on each partition. Compute the internal node’s centroid as the weighted mean of child centroids.

The overall build complexity is O(n · d · log n) in the average case, dominated by the O(log n) levels of recursion with O(n · d) work at each level for variance computation and partitioning.

Beam Search Traversal

At query time, HBC uses beam search — a bounded-width priority search — to find the approximate nearest neighbors. The algorithm maintains a priority queue (the “beam”) of candidate nodes, ordered by distance from the query vector to the node’s centroid:

1. Initialize the beam with the root node. Set the beam width W (a tunable parameter that controls the recall–latency trade-off).
2. Pop the closest node from the beam. If it is a leaf, compute exact distances to all vectors in the leaf bucket and update the result set.
3. If it is an internal node, compute distances from the query to each child’s centroid. Apply pruning: skip any child whose centroid distance minus bounding radius exceeds the current k-th nearest distance (if k results have been found).
4. Add non-pruned children to the beam. If the beam exceeds width W, discard the most distant entries.
5. Repeat until the beam is empty or a maximum number of distance computations has been reached.

Pruning with Distance Bounds

The triangle inequality provides the theoretical basis for pruning. For a query vector q, a node centroid c, and bounding radius r, the minimum possible distance from q to any vector in the subtree is dist(q, c) - r. If this lower bound exceeds the distance to the current k-th nearest neighbor, no vector in the subtree can improve the result set. This pruning is what makes HBC competitive with graph-based methods for moderate-dimensional data (up to ~128 dimensions), where the bounding radii remain tight enough to prune effectively.

Practical Applications

Vector Indexing in a Production Search System

A production vector search system implementing HBC uses it as the primary vector index structure. When a collection of vectors is ingested — whether from document embeddings generated by a language model, image feature vectors from a vision encoder, or any other high-dimensional representation — the system builds an HBC tree to enable fast approximate nearest neighbor queries. A Zig implementation, for example, takes advantage of the language’s explicit memory management and SIMD-friendly data layouts to achieve high throughput during both bulk builds and queries.

The HBC index supports incremental updates through a combination of leaf-level insertions and periodic re-balancing. When a new vector arrives, it traverses the tree to find the most appropriate leaf, inserts into the leaf bucket, and triggers a split if the bucket overflows. Periodic background tasks re-evaluate split strategies for subtrees that have become imbalanced due to insertions, ensuring long-term query performance.

Hybrid Retrieval Systems

HBC trees integrate naturally with sparse retrieval methods (such as BM25 or SPLADE) in hybrid search pipelines. The tree structure supports efficient batch queries where multiple query vectors share beam state, reducing redundant centroid distance computations. This makes HBC a practical choice for re-ranking stages where a set of candidate documents must be scored against a query embedding.

Streaming & Bulk Ingestion

The bulk build algorithm is particularly valuable during initial data loading or migration scenarios. Rather than inserting vectors one at a time (which can degrade tree balance), a well-designed system buffers incoming vectors and triggers a full rebuild when the buffer reaches a configurable threshold. This amortized approach keeps the tree well-balanced and avoids the pathological cases that plague incremental tree construction.

Implementation Patterns

Memory Layout for Cache Efficiency

A performant HBC implementation stores the tree in a flat array using a breadth-first or van Emde Boas layout. Internal nodes are packed contiguously so that traversing from parent to child results in sequential memory access patterns. Centroid vectors for internal nodes at the same depth are stored adjacently, enabling SIMD-parallel distance computations when evaluating multiple children of a node simultaneously.

SIMD-Accelerated Distance Computation

Distance computation (typically L2 or cosine) is the bottleneck in both build and query phases. High-performance implementations use SIMD intrinsics to compute distances over 4, 8, or 16 dimensions per instruction. In Zig, the @Vector built-in type maps directly to hardware SIMD registers, enabling the compiler to auto-vectorize distance loops without manual intrinsic calls.

Configurable Branching Factor

While binary trees (branching factor 2) are the simplest HBC variant, higher branching factors (4, 8, or 16) reduce tree depth at the cost of more centroid comparisons per node. The optimal branching factor depends on the dimensionality d: for low-dimensional data (d < 32), higher branching factors reduce tree traversal overhead; for high-dimensional data (d > 128), binary splits with careful axis selection tend to produce tighter clusters.

Thread-Safe Concurrent Queries

Since the tree structure is immutable after the build phase (until a re-balance is triggered), multiple query threads can traverse the tree concurrently without locks. Each query thread maintains its own beam and result set on the stack, sharing only read access to the tree’s node array. This shared-nothing query model scales linearly with the number of CPU cores.

Code Examples

Zig: HBC Node Structure & Split Logic

The following Zig snippet illustrates the core node structure and the even-split partitioning logic used during the bulk build phase:

const std = @import("std");
const Allocator = std.mem.Allocator;

pub const SplitStrategy = enum {
    even,
    left,
    right,
};

pub const HbcNode = struct {
    centroid: []f32,
    bounding_radius: f32,
    left_child: ?*HbcNode,
    right_child: ?*HbcNode,
    vectors: ?[][]f32, // non-null only for leaf nodes

    pub fn isLeaf(self: *const HbcNode) bool {
        return self.vectors != null;
    }
};

/// Select the dimension with highest variance for splitting.
fn selectSplitAxis(vectors: [][]f32, dims: usize) usize {
    var best_axis: usize = 0;
    var best_var: f32 = 0.0;

    for (0..dims) |axis| {
        var sum: f32 = 0.0;
        var sum_sq: f32 = 0.0;
        for (vectors) |vec| {
            sum += vec[axis];
            sum_sq += vec[axis] * vec[axis];
        }
        const n: f32 = @floatFromInt(vectors.len);
        const variance = (sum_sq / n) - (sum / n) * (sum / n);
        if (variance > best_var) {
            best_var = variance;
            best_axis = axis;
        }
    }
    return best_axis;
}

/// Partition vectors using the specified split strategy.
/// Returns the split index into the sorted array.
fn computeSplitIndex(
    n: usize,
    strategy: SplitStrategy,
) usize {
    return switch (strategy) {
        .even => n / 2,
        .left => (n * 2) / 3,   // 66% left, 34% right
        .right => n / 3,         // 34% left, 66% right
    };
}

/// Build the HBC tree recursively.
pub fn buildTree(
    allocator: Allocator,
    vectors: [][]f32,
    dims: usize,
    leaf_capacity: usize,
) !*HbcNode {
    var node = try allocator.create(HbcNode);

    if (vectors.len <= leaf_capacity) {
        // Leaf node: store vectors directly
        node.* = .{
            .centroid = try computeCentroid(allocator, vectors, dims),
            .bounding_radius = computeBoundingRadius(vectors, node.centroid),
            .left_child = null,
            .right_child = null,
            .vectors = vectors,
        };
        return node;
    }

    const axis = selectSplitAxis(vectors, dims);
    // Sort along the chosen axis
    std.sort.pdq([]f32, vectors, axis, struct {
        fn lessThan(ctx: usize, a: []f32, b: []f32) bool {
            return a[ctx] < b[ctx];
        }
    }.lessThan);

    // Evaluate each split strategy and pick the best
    const strategy = chooseBestSplit(vectors, dims, axis);
    const split_idx = computeSplitIndex(vectors.len, strategy);

    node.left_child = try buildTree(allocator, vectors[0..split_idx], dims, leaf_capacity);
    node.right_child = try buildTree(allocator, vectors[split_idx..], dims, leaf_capacity);

    node.centroid = try computeWeightedCentroid(
        allocator, node.left_child.?, node.right_child.?, dims,
    );
    node.bounding_radius = computeBoundingRadius(vectors, node.centroid);
    node.vectors = null;

    return node;
}

Python: Beam Search Query

This Python example demonstrates the beam search traversal algorithm used at query time. It shows how the priority queue maintains candidates and how triangle-inequality pruning eliminates subtrees:

import heapq
import numpy as np
from dataclasses import dataclass, field
from typing import Optional


@dataclass
class HbcNode:
    centroid: np.ndarray
    bounding_radius: float
    left: Optional["HbcNode"] = None
    right: Optional["HbcNode"] = None
    vectors: Optional[np.ndarray] = None  # leaf only

    @property
    def is_leaf(self) -> bool:
        return self.vectors is not None


def beam_search(
    root: HbcNode,
    query: np.ndarray,
    k: int = 10,
    beam_width: int = 64,
) -> list[tuple[float, int]]:
    """
    Perform beam search over an HBC tree.

    Returns the k approximate nearest neighbors as
    (distance, vector_index) tuples.
    """
    # Priority queue: (distance_to_centroid, node)
    beam: list[tuple[float, HbcNode]] = []
    root_dist = float(np.linalg.norm(query - root.centroid))
    heapq.heappush(beam, (root_dist, root))

    # Result set: max-heap of (neg_distance, index)
    results: list[tuple[float, int]] = []
    kth_dist = float("inf")

    while beam:
        dist_to_centroid, node = heapq.heappop(beam)

        # Prune: if lower bound exceeds current kth distance
        lower_bound = dist_to_centroid - node.bounding_radius
        if lower_bound > kth_dist and len(results) >= k:
            continue

        if node.is_leaf:
            # Exhaustive scan of the leaf bucket
            for idx, vec in enumerate(node.vectors):
                d = float(np.linalg.norm(query - vec))
                if len(results) < k:
                    heapq.heappush(results, (-d, idx))
                    if len(results) == k:
                        kth_dist = -results[0][0]
                elif d < kth_dist:
                    heapq.heapreplace(results, (-d, idx))
                    kth_dist = -results[0][0]
        else:
            # Expand children into the beam
            for child in (node.left, node.right):
                if child is None:
                    continue
                child_dist = float(np.linalg.norm(query - child.centroid))
                child_lower = child_dist - child.bounding_radius

                # Prune children that cannot improve results
                if child_lower > kth_dist and len(results) >= k:
                    continue

                heapq.heappush(beam, (child_dist, child))

        # Trim beam to width W
        if len(beam) > beam_width:
            beam = heapq.nsmallest(beam_width, beam)
            heapq.heapify(beam)

    # Convert from max-heap to sorted results
    return sorted([(-d, i) for d, i in results])

Common Pitfalls

1. Choosing a leaf capacity that is too small. Setting the leaf bucket size below 32 vectors creates excessively deep trees where centroid traversal overhead dominates. The cost of visiting many internal nodes outweighs the savings from smaller leaf scans. Start with a leaf capacity of 64–128 and tune based on your dimensionality.
2. Using a fixed split strategy instead of adaptive selection. Applying even splits universally fails on real-world data with skewed distributions. Dense clusters end up sharing subtrees with sparse regions, producing loose bounding radii that defeat pruning. Always evaluate all three split strategies per node and select the one that minimizes bounding radius.
3. Neglecting to recompute centroids after incremental insertions. When vectors are inserted into leaf nodes without updating ancestor centroids, the centroid hierarchy drifts from the actual data distribution. This causes the pruning bounds to become stale, silently degrading recall. Implement either lazy centroid propagation or periodic full recomputation.
4. Setting beam width too low for high-dimensional data. In dimensions above 64, the curse of dimensionality loosens bounding radii, requiring wider beams to maintain acceptable recall. A beam width of 32 might achieve 95% recall at d=16 but only 70% at d=128. Benchmark recall at your target dimensionality and increase beam width accordingly.
5. Ignoring memory alignment for SIMD distance computations. Centroid and vector arrays that are not aligned to 16-byte (SSE) or 32-byte (AVX) boundaries force the CPU to perform unaligned loads, which can halve throughput on distance-intensive workloads. In Zig, use @alignCast and allocate vectors with explicit alignment to match your target SIMD width.

Key Takeaways