Making np.searchsorted up to 25× Faster in NumPy 2.5

NumPy's `searchsorted` function, which implements the binary search algorithm, has been significantly optimized in NumPy 2.5, resulting in up to a 25x speedup in benchmarks. This operation is crucial in the Python scientific ecosystem, used for tasks such as histogram computation and interval lookups. Several techniques, including branch elimination, batching, and cache-friendly data layouts, have been explored to improve binary search performance.

In this report, we delve into how these optimizations can be applied using NumPy's vectorized primitives.

The `searchsorted` function locates the insertion position of a sequence of query keys in a static sorted array. Traditionally, a pure-Python implementation runs a binary search for each key, with the running time for each query growing logarithmically as the input size increases. To benchmark this, we generated two random arrays of uniformly distributed 32-bit integers.

The keys, which are the elements being searched, had a fixed length of 10,000, while the length of the sorted array varied up to $2^{30}$ (4 GiB). Both keys and values arrays were contiguous in memory. Each benchmark was repeated 50 times, and the minimum execution time was reported. The benchmarks were conducted on a MacBook Pro with an Apple M1 Pro and 32 GB of memory.

In the baseline implementation, a binary search is performed for each query sequentially in Python. To enhance performance, we adapted the algorithm to process multiple binary searches simultaneously in batches. By representing the state of all searches as arrays and updating them concurrently using vectorized operations, we removed Python overhead and enabled the CPU to efficiently handle large batches of independent work.

NumPy's compiled C++ loops execute operations on arrays, which eliminates Python overhead and allows for efficient processing of large batches of independent work. To vectorize the algorithm, we treated scalar variables as array states, maintaining the search intervals (lo and hi) for each query instead of scalar variables. This allows for simultaneous updates, improving overall performance.

However, different queries may require different numbers of iterations to converge, especially when searching for keys that are not present in the array. To address this, we introduced an active mask to identify queries whose search intervals have not yet collapsed to a single position. By updating only the active queries at each iteration, we can significantly reduce the number of iterations required.

The iteration process continues until the insertion position of each key lies within the range [lo, hi), with the interval shrinking to roughly half its size at each step. After approximately np.ceil(np.log2(n)) iterations, every interval collapses to a single position.

This vectorized implementation outperforms NumPy's native `searchsorted` function, which executes a sequential binary search for each key independently. The native implementation, while straightforward, suffers from dependency chains within each search, where each iteration depends on the result of the previous one. This creates a chain of dependencies that can lead to cache-unfriendly memory access patterns, particularly for large arrays.

Additionally, sequential binary searches tend to be cache-unfriendly due to the uneven distribution of memory accesses required for each step. In contrast, the vectorized implementation performs the same logical step across all queries simultaneously, allowing the CPU to issue several memory accesses in parallel. This improves cache utilization and reduces the impact of cache misses, resulting in faster execution times once the array size exceeds the L1 and L2 cache sizes.