Hilbert Space-Filling Curves

Overview

A Hilbert curve is a continuous fractal curve that visits every point in a multi-dimensional space exactly once, mapping it to a single one-dimensional index. In the context of databases and vector search engines, Hilbert curves provide a powerful mechanism for converting high-dimensional data—such as embedding vectors or geographic coordinates—into a linear ordering that preserves spatial locality. Points that are close together in the original multi-dimensional space tend to have nearby Hilbert indices, which means range scans over a one-dimensional B-tree index can efficiently retrieve spatially clustered results.

The locality-preserving property is what sets Hilbert curves apart from simpler alternatives like Z-order (Morton) curves. Z-order curves interleave the bits of each coordinate directly, which produces a mapping that preserves locality reasonably well but suffers from significant discontinuities at certain boundaries. When traversing a Z-order curve, the path can make large jumps across the space at power-of-two boundaries, resulting in clusters of nearby points being split across distant index ranges. Hilbert curves eliminate these jumps by rotating and reflecting sub-quadrants at each level of recursion, producing a mapping with provably better worst-case locality.

In a vector database, Hilbert curves can be used to linearize multi-dimensional vector data so that standard one-dimensional indexing structures—B-trees, sorted runs in LSM trees—can serve as spatial indexes without requiring dedicated multi-dimensional data structures like R-trees or KD-trees. This approach simplifies the storage engine while still enabling efficient approximate nearest-neighbor queries and spatial range scans.

Core Concepts

Space-Filling Curves

A space-filling curve is a continuous curve whose range contains the entirety of a multi-dimensional unit square (or hypercube). Discovered by Giuseppe Peano in 1890 and refined by David Hilbert in 1891, these curves establish a bijection between one-dimensional and multi-dimensional spaces. The key insight is that a single scalar index can encode a position in n-dimensional space without losing spatial relationships.

Fractal Recursion

Hilbert curves are constructed recursively. At order 1, the curve visits the four quadrants of a 2D grid in a U-shaped pattern. At order 2, each quadrant is subdivided into four sub-quadrants, and the curve is applied recursively with specific rotations and reflections to maintain continuity. At order k, the curve covers a 2k × 2k grid, visiting all 4k cells. Each increase in order doubles the resolution along each axis.

Dimensionality Reduction

By mapping n-dimensional coordinates to a single Hilbert index, multi-dimensional search becomes a series of one-dimensional range queries. This is a form of dimensionality reduction that preserves enough spatial structure for clustering and approximate nearest neighbor search, while allowing the use of conventional sorted indexes.

How It Works

Recursive Construction

The 2D Hilbert curve is built by dividing the space into four quadrants and assigning each a sub-curve orientation. The base pattern at order 1 traces a U-shape through quadrants numbered 0 through 3. At each higher order, the algorithm:

1. Divides the current grid into four equal quadrants.
2. Applies the Hilbert curve recursively to each quadrant.
3. Rotates the sub-curve in quadrant 0 by 90° clockwise (swapping x and y).
4. Rotates the sub-curve in quadrant 3 by 90° counter-clockwise (swapping and flipping x and y).
5. Connects the four sub-curves into a single continuous path.

These rotations ensure that the exit point of one quadrant’s sub-curve aligns with the entry point of the next, maintaining continuity and locality.

Coordinate-to-Index (Encoding)

To encode a point (x, y) into a Hilbert index at order k, the algorithm processes the coordinates from the most significant bit to the least significant bit. At each level, it determines which quadrant the point falls into, appends the corresponding 2 bits to the index, and then applies the appropriate rotation to the coordinates before descending to the next level. The total index is a 2k-bit integer.

Index-to-Coordinate (Decoding)

Decoding reverses the process: starting from the most significant pair of bits in the index, the algorithm determines the quadrant, sets the corresponding coordinate bits, and applies the inverse rotation. After processing all k levels, the full (x, y) coordinates are recovered.

The Encoding State Machine

The 2D Hilbert encoding can be viewed as a finite state transducer with 4 states. Each state encodes a specific orientation of the U-shaped base curve. Given the current state and the 2-bit quadrant selector derived from the coordinates, the transducer emits 2 index bits and transitions to a new state. The entire encode/decode process requires exactly k iterations for an order-k curve, making it O(k) per point—or equivalently O(log N) where N is the grid resolution.

Practical Applications

Vector Database Clustering

In a typical LSM-tree-based vector store, Hilbert curves are used to sort high-dimensional embedding vectors into a linear order before inserting them. By computing the Hilbert index of each vector (after quantization to a fixed grid resolution), spatially similar vectors end up in the same or adjacent SSTable blocks. This dramatically improves cache hit rates during nearest-neighbor scans because a sequential read of the sorted index retrieves vectors that are also close in the embedding space.

R-Tree Alternatives

Traditional spatial databases use R-trees to index multi-dimensional data, but R-trees have significant drawbacks in LSM-tree-based storage engines: they require in-place updates, have complex split and merge logic, and are difficult to serialize into sorted runs. Hilbert curve indexing sidesteps these problems entirely. A Zig-based storage engine can use Hilbert indices as B-tree keys, turning spatial queries into simple range scans. Bulk loading is also simplified—sorting by Hilbert index and writing sequentially produces a naturally clustered index with no post-processing needed.

Geographic Indexing

For geospatial data, Hilbert curves provide an alternative to geohash-based indexing. Geohashes (which are essentially Z-order curves applied to latitude/longitude) suffer from the boundary discontinuity problem: points that are geographically adjacent may fall into different geohash prefixes if they straddle a power-of-two boundary. Hilbert curve indices mitigate this by ensuring smoother traversal across boundaries, reducing the number of index ranges needed to cover a circular query region.

Batch Ingestion Optimization

Storage engines using LSM trees commonly apply Hilbert curve sorting during bulk data ingestion. When a large batch of vectors arrives, the pipeline computes the Hilbert index for each vector, sorts the batch by index, and writes the sorted data as a single SSTable or B-tree leaf sequence. This ensures optimal spatial clustering without requiring expensive post-insertion reorganization.

Implementation Patterns

Lookup Tables for Small Dimensions

For 2D and 3D Hilbert curves, the state machine transitions can be precomputed into small lookup tables. A 2D curve needs a 4×4 table (4 states × 4 quadrants), encoding both the output bits and the next state in a single byte. This eliminates all branching from the inner loop, making encoding/decoding extremely fast—typically under 10 nanoseconds per point for order-16 curves.

Generalization to Higher Dimensions

The Hilbert curve generalizes to n dimensions, though the number of states grows as n! × 2n. For typical vector database workloads with 128 or 768 dimensions, direct Hilbert encoding is impractical. Instead, the common approach is:

• Apply dimensionality reduction (PCA or random projection) to reduce vectors to 8–16 dimensions.
• Quantize each reduced dimension to a fixed number of bits (e.g., 4–8 bits).
• Compute the Hilbert index over the reduced, quantized representation.

This trades some locality accuracy for computational feasibility, but experiments show that the combination of PCA + Hilbert ordering retains strong spatial clustering for approximate nearest-neighbor workloads.

Integration with B-Tree Indexes

Hilbert indices are stored as fixed-width unsigned integers (typically 64 or 128 bits) and used directly as B-tree keys. Range queries in the original space are translated into one or more Hilbert index ranges using a covering algorithm: the query region is recursively subdivided along the Hilbert curve until each sub-region maps to a contiguous index range. The resulting set of ranges is then executed as a union of B-tree range scans.

Batch Sorting

For bulk loading, the standard pattern is to compute Hilbert indices for all points, sort the array of (index, point) pairs using radix sort (since indices are fixed-width integers), and then write the sorted data sequentially. Radix sort runs in O(n) time for fixed-width keys, making it significantly faster than comparison-based sorts for large batches.

Code Examples

Zig: Hilbert Curve Coordinate Encoding/Decoding

This example implements 2D Hilbert curve encoding and decoding using a lookup-table-based state machine. The order parameter controls the grid resolution (2order × 2order).

const std = @import("std");

/// Hilbert curve state transition table for 2D encoding.
/// Each entry: { output_bits, next_state }
/// States represent the 4 orientations of the base U-curve.
const encode_table = [4][4]struct { u2, u2 }{
    // State 0 (standard U)
    .{ .{ 0, 1 }, .{ 1, 0 }, .{ 3, 3 }, .{ 2, 0 } },
    // State 1 (90deg CW rotation)
    .{ .{ 0, 0 }, .{ 3, 1 }, .{ 1, 1 }, .{ 2, 2 } },
    // State 2 (180deg rotation)
    .{ .{ 2, 2 }, .{ 1, 2 }, .{ 3, 0 }, .{ 0, 3 } },
    // State 3 (90deg CCW rotation)
    .{ .{ 2, 3 }, .{ 3, 3 }, .{ 1, 0 }, .{ 0, 1 } },
};

/// Encode (x, y) coordinates into a Hilbert index.
pub fn hilbertEncode(x: u32, y: u32, order: u5) u64 {
    var index: u64 = 0;
    var state: u2 = 0;

    var i: u5 = order;
    while (i > 0) {
        i -= 1;
        const bx: u2 = @truncate((x >> @intCast(i)) & 1);
        const by: u2 = @truncate((y >> @intCast(i)) & 1);
        const quadrant: u2 = (by << 1) | bx;

        const entry = encode_table[state][quadrant];
        index = (index << 2) | @as(u64, entry[0]);
        state = entry[1];
    }

    return index;
}

/// Decode a Hilbert index back into (x, y) coordinates.
pub fn hilbertDecode(index: u64, order: u5) struct { u32, u32 } {
    var x: u32 = 0;
    var y: u32 = 0;
    var state: u2 = 0;

    // Reverse lookup: for each state, map output bits -> quadrant
    const decode_table = comptime blk: {
        var table: [4][4]struct { u2, u2 } = undefined;
        for (0..4) |s| {
            for (0..4) |q| {
                const out = encode_table[s][q][0];
                const ns = encode_table[s][q][1];
                table[s][out] = .{ @intCast(q), ns };
            }
        }
        break :blk table;
    };

    var i: u5 = order;
    while (i > 0) {
        i -= 1;
        const bits: u2 = @truncate((index >> @intCast(i * 2)) & 3);
        const entry = decode_table[state][bits];
        const quadrant = entry[0];
        state = entry[1];

        x = (x << 1) | @as(u32, quadrant & 1);
        y = (y << 1) | @as(u32, (quadrant >> 1) & 1);
    }

    return .{ x, y };
}

test "hilbert round-trip" {
    const order: u5 = 8;
    const size: u32 = 1 << order;
    for (0..size) |xi| {
        for (0..size) |yi| {
            const x: u32 = @intCast(xi);
            const y: u32 = @intCast(yi);
            const idx = hilbertEncode(x, y, order);
            const decoded = hilbertDecode(idx, order);
            try std.testing.expectEqual(x, decoded[0]);
            try std.testing.expectEqual(y, decoded[1]);
        }
    }
}

Go: Using Hilbert Indices for Spatial Clustering

This example demonstrates how to sort a set of 2D points by their Hilbert index to achieve spatial clustering, then partition them into fixed-size clusters suitable for storage blocks.

package hilbert

import (
	"math/bits"
	"sort"
)

// HilbertIndex computes the Hilbert curve index for a 2D point
// at the given order (grid is 2^order x 2^order).
func HilbertIndex(x, y uint32, order int) uint64 {
	var index uint64
	for i := order - 1; i >= 0; i-- {
		bx := (x >> uint(i)) & 1
		by := (y >> uint(i)) & 1

		// Determine quadrant
		quadrant := by<<1 | bx

		// Rotate the coordinates for the next level
		switch quadrant {
		case 0: // bottom-left: swap x and y
			x, y = y, x
		case 1: // top-left: no rotation
			// nothing
		case 2: // top-right: no rotation
			// nothing
		case 3: // bottom-right: swap and flip
			n := uint32(1<<uint(i)) - 1
			x, y = n-y, n-x
		}

		// Map quadrant to Hilbert index bits
		hilbertBits := [4]uint64{0, 1, 3, 2}
		index = (index << 2) | hilbertBits[quadrant]
	}
	return index
}

// Point represents a 2D coordinate with associated data.
type Point struct {
	X, Y    uint32
	Payload []byte
	hilbert uint64 // cached Hilbert index
}

// ClusterPoints sorts points by Hilbert index and partitions
// them into clusters of the specified size. The final cluster
// may be smaller.
func ClusterPoints(points []Point, order int, clusterSize int) [][]Point {
	if len(points) == 0 {
		return nil
	}

	// Compute Hilbert indices
	for i := range points {
		points[i].hilbert = HilbertIndex(
			points[i].X, points[i].Y, order,
		)
	}

	// Sort by Hilbert index for spatial clustering
	sort.Slice(points, func(i, j int) bool {
		return points[i].hilbert < points[j].hilbert
	})

	// Partition into fixed-size clusters
	var clusters [][]Point
	for i := 0; i < len(points); i += clusterSize {
		end := i + clusterSize
		if end > len(points) {
			end = len(points)
		}
		clusters = append(clusters, points[i:end])
	}

	return clusters
}

// CoveringRanges computes the set of Hilbert index ranges that
// cover a rectangular query region. This allows translating a
// 2D range query into a set of 1D B-tree range scans.
func CoveringRanges(
	minX, minY, maxX, maxY uint32,
	order int,
) [][2]uint64 {
	var ranges [][2]uint64

	// Compute Hilbert indices for all grid cells in the region
	var indices []uint64
	for y := minY; y <= maxY; y++ {
		for x := minX; x <= maxX; x++ {
			indices = append(indices, HilbertIndex(x, y, order))
		}
	}

	// Sort indices
	sort.Slice(indices, func(i, j int) bool {
		return indices[i] < indices[j]
	})

	// Merge consecutive indices into ranges
	if len(indices) == 0 {
		return nil
	}
	start := indices[0]
	prev := indices[0]
	for _, idx := range indices[1:] {
		if idx == prev+1 {
			prev = idx
			continue
		}
		ranges = append(ranges, [2]uint64{start, prev})
		start = idx
		prev = idx
	}
	ranges = append(ranges, [2]uint64{start, prev})

	return ranges
}

// Ensure bits import is used (HilbertIndex uses bit operations)
var _ = bits.Len

Common Pitfalls

1. Choosing insufficient curve order. If the curve order is too low, multiple distinct points map to the same Hilbert index, destroying spatial resolution. The order must be at least ceil(log2(max_coordinate_value)) for each dimension.
2. Ignoring dimensionality limits. The Hilbert state machine grows factorially with dimensions. Applying Hilbert encoding directly to 128-dimensional vectors is computationally infeasible. Always reduce dimensionality first (e.g., via PCA or random projection to 8–16 dimensions).
3. Using Z-order when Hilbert is needed. Z-order curves are simpler and faster to compute, which makes them tempting. However, for workloads involving range queries or nearest-neighbor search, the boundary discontinuities in Z-order curves produce 20–30% more I/O than Hilbert curves. Profile your workload before choosing.
4. Forgetting to quantize floating-point coordinates. Hilbert curves operate on integer grids. Floating-point coordinates must be quantized to a fixed resolution before encoding. Inconsistent quantization (e.g., different scales per dimension) will distort the spatial mapping and degrade clustering quality.
5. Generating too many covering ranges for large queries. A spatial query that spans a large portion of the grid may decompose into hundreds or thousands of Hilbert index ranges. Each range becomes a separate B-tree scan, and the overhead of initiating many small scans can outweigh the benefit. Set a maximum range count and fall back to a full scan when the threshold is exceeded.

Key Takeaways