Continuum Website Continuum Applications Continuum Knowledge Axolotl Platform

Billion-scale similarity search with GPUs

PreviousGenerative Engine Optimisation (GEO)NextFOLLOWIR: Evaluating and Teaching Information Retrieval Models to Follow Instructions

Last updated 11 months ago

Was this helpful?

This famous 2017 paper focuses on improving the performance of similarity search on large-scale datasets using GPUs.

Similarity search is a fundamental operation in many applications, particularly those involving high-dimensional data like images and videos.

The paper proposes several optimisations to better utilise GPU resources and achieve significant speedups compared to prior state-of-the-art methods.

Key Points

Similarity Search

The goal is to find the k-nearest neighbours $(k-NN)$ of a query vector $x$ within a large collection of vectors $[yi]$ based on the Euclidean distance.
Exact search involves computing the full pairwise distance matrix between the query vectors and the database vectors, which is computationally expensive for large datasets.

Compressed-Domain Search

To improve search efficiency, approximate nearest-neighbor search methods are used, such as the IVFADC (Inverted File with Asymmetric Distance Computation) indexing structure.
IVFADC relies on two levels of quantization: a coarse quantizer $(q1)$ and a fine quantizer $(q2)$ that encodes the residual vector after the first level.
The search is performed in two steps: first, the coarse quantizer is used to select a subset of inverted lists (LIVF), and then the fine quantizer is used to compute distances only for the vectors in the selected inverted lists (LIVFADC).

Product Quantizer

A product quantizer is used for the fine quantizer $(q2)$ to provide a large number of reproduction values without increasing the processing cost.
The vector is divided into sub-vectors, each quantized with its own sub-quantizer, typically with 256 reproduction values to fit in one byte.
The product quantizer generates b-byte codes, where b is the number of sub-vectors, resulting in a large codebook size $(|C2| = 256^b)$ while keeping the quantization process computationally cheap.

GPU Optimisation

The paper proposes a novel GPU k-selection algorithm that operates in fast register memory and can be fused with other kernels, achieving up to 55% of theoretical peak performance.
The authors present an optimised algorithmic layout for exact and approximate k-NN search on GPUs, outperforming previous state-of-the-art methods by a large margin.

Experiments and Results

The proposed approach is evaluated on various similarity search scenarios, including brute-force search, approximate search, and compressed-domain search using product quantization.
The implementation achieves significant speedups compared to prior GPU-based methods, enabling the construction of a high-accuracy k-NN graph on 95 million images from the Yfcc100M dataset in 35 minutes and a graph connecting 1 billion vectors in less than 12 hours on 4 Maxwell Titan X GPUs.

GPU Architecture

GPUs (Graphics Processing Units) are specialised processors designed for parallel processing and are particularly well-suited for tasks that involve large amounts of data and computations.
GPUs consist of thousands of cores (processing units) that can execute instructions simultaneously, enabling massive parallelism.
The basic execution unit on a GPU is called a warp, which is a group of 32 threads that execute the same instruction in lockstep.
Threads within a warp are referred to as lanes, and each lane has its own registers for storing data.
Warps are organised into blocks (or cooperative thread arrays), which can contain up to 32 warps and have access to a shared memory space.
Blocks are further organised into a grid, which represents the entire workload of a GPU kernel (a function executed on the GPU).

Memory Hierarchy

GPUs have a memory hierarchy similar to CPUs but with some differences.
Global memory is the largest and slowest memory on the GPU, accessible by all threads across all blocks. It is used for communication between the CPU and GPU and for storing large datasets.
Shared memory is a fast on-chip memory that is shared among threads within a block. It is used for efficient data sharing and communication between threads.
Registers are the fastest memory on the GPU and are private to each thread. They are used for storing frequently accessed variables and temporary data.
Texture memory and constant memory are read-only caches optimized for specific access patterns.

Programming GPUs

CUDA (Compute Unified Device Architecture) is a popular programming model developed by NVIDIA for programming their GPUs.
CUDA allows developers to write code in languages like C, C++, or Fortran and provides extensions to define GPU kernels and manage GPU memory.
OpenCL (Open Computing Language) is an open standard for parallel programming across different devices, including GPUs, CPUs, and FPGAs.
When programming GPUs, you typically define kernels, which are functions that are executed in parallel on the GPU threads.
Kernels are launched with a specific grid and block configuration, determining the number of threads and their organisation.
Inside a kernel, you can access thread and block IDs to identify the specific thread and its position within the grid.
You can transfer data between the CPU (host) and GPU (device) using memory management functions provided by CUDA or OpenCL.
Synchronisation primitives, such as barriers and atomic operations, are used to coordinate the execution and data sharing between threads.

Example: Vector Addition on GPU

__global__ void vectorAdd(int *a, int *b, int *c, int size) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    if (tid < size) {
        c[tid] = a[tid] + b[tid];
    }
}

int main() {
    int size = 1024;
    int *a, *b, *c;
    cudaMallocManaged(&a, size * sizeof(int));
    cudaMallocManaged(&b, size * sizeof(int));
    cudaMallocManaged(&c, size * sizeof(int));

    // Initialize input vectors
    for (int i = 0; i < size; i++) {
        a[i] = i;
        b[i] = i * 2;
    }

    int blockSize = 256;
    int numBlocks = (size + blockSize - 1) / blockSize;
    vectorAdd<<<numBlocks, blockSize>>>(a, b, c, size);

    cudaDeviceSynchronize();

    // Print the result
    for (int i = 0; i < size; i++) {
        printf("%d ", c[i]);
    }

    cudaFree(a);
    cudaFree(b);
    cudaFree(c);
    return 0;
}

In this example, we define a GPU kernel called vectorAdd that performs element-wise addition of two vectors a and b and stores the result in vector c.

Each thread is responsible for adding one element based on its thread ID.

In the main function, we allocate memory on the GPU using cudaMallocManaged, initialise the input vectors, and launch the kernel with a specific grid and block configuration.

Finally, we synchronise the device to ensure all computations are complete and print the result.

This is just a simple example, but GPUs can be used for a wide range of applications, including image and video processing, machine learning, scientific simulations, and more.

By leveraging the massive parallelism and high memory bandwidth of GPUs, you can accelerate computationally intensive tasks and achieve significant performance gains compared to CPUs.

Computation layout for implementing the IVFADC indexing method

IVFADC is an approximate nearest neighbor search method that combines the inverted file index with asymmetric distance computation (ADC) based on product quantization.

The key aspects of the implementation are the distance computations and their integration with the k-selection process.

Exact search

Exhaustive search, also known as exact brute-force search, is used for small datasets or as a component of other indexing methods.
In the context of IVFADC, exact search is used for the coarse quantizer $(q1).$

Distance computation

The distance computation for exact search boils down to a matrix multiplication.
Optimised GEMM (General Matrix Multiplication) routines from the cuBLAS library are used to calculate the $-2<x_j, y_i>$ term for L2 distance.
The GEMM operation results in a partial distance matrix $D$ .

Fused k-selection kernel

To complete the distance calculation, a fused k-selection kernel is employed.
The fused kernel adds the $||y_i||^2$ term to each entry of the partial distance matrix $D$ and immediately submits the value to k-selection in registers.
The $||x_j||^2$ term is not required before k-selection and can be omitted.
Kernel fusion allows for only 2 passes over D (GEMM write and k-select read), reducing memory access compared to implementations that may require 3 or more passes.

Tiling and memory management

For realistic problem sizes, the distance matrix $D$ may not fit in GPU memory.
To handle this, the problem is tiled over the batch of queries, with $t_q ≤ n_q$ queries being run in a single tile.
Each tile represents an independent problem, but two tiles are run in parallel on different streams to better utilise the GPU.
The effective memory requirement for $D$ becomes $O(2 × t_q)$ .
The computation can also be tiled over the database size $(n)$ if necessary.

Data transfer and overlap

For very large input coming from the CPU, pinned memory buffering is used to overlap CPU-to-GPU data transfer with GPU computation.
This helps to hide the latency of data transfer and improve overall performance.

The efficient implementation of IVFADC on GPUs relies on several key optimisations:

Using optimised GEMM routines from cuBLAS for the distance computation matrix multiplication.
Fusing the k-selection kernel with the final distance calculation to reduce memory access and passes over the distance matrix.
Tiling the problem over the batch of queries and database size to handle large datasets that exceed GPU memory capacity.
Running multiple tiles in parallel on different streams to maximize GPU utilisation.
Overlapping CPU-to-GPU data transfer with computation using pinned memory buffering.

By carefully designing the computation layout and leveraging GPU-specific optimizations, the IVFADC method can be implemented efficiently on GPUs, outperforming other GPU-compliant approximate nearest neighbor search strategies.

The combination of optimized distance computations and tight integration with the k-selection process enables high-performance similarity search on large-scale datasets.

Non-exhaustive search with PQ codes

The key idea behind using product quantization (PQ) codes for non-exhaustive search is to compute the distance between a query vector and a set of reproduction values using lookup tables.

This approach significantly reduces the computational cost compared to explicit distance calculations.

Key points

PQ lookup tables

When the query vector x and the coarse quantizer reproduction value $q1(y)$ are known, all distances to the PQ reproduction values can be precomputed and stored in lookup tables $T1, ..., Tb$ .
Each lookup table is of size 256, corresponding to the number of reproduction values for each sub-quantizer.
Computing the distance between $x$ and $q(y)$ consists of $b$ lookups and additions

IVFADC lookup tables

When scanning the elements of an inverted list $IL$ , where $q1(y)$ is constant, the lookup table method can be applied.
The distance computation is further optimised by decomposing it into three terms
Term 1 is independent of the query and can be precomputed and stored in a table of size $|C1| × 256 × b.$
Term 2 is the distance to $q1's$ reproduction value and is a by-product of the first-level quantizer $q1$ .
Term 3 can be computed independently of the inverted list, with a cost of $d × 256$ multiply-adds.

GPU implementation

The inverted lists are stored as two separate arrays for $PQ$ codes and associated IDs.
A kernel is responsible for scanning the τ closest inverted lists for each query and calculating the per-vector pair distances using the lookup tables $Ti$ .
The lookup tables are stored in shared memory to enable fast access.
Multi-pass kernels are used to process the $nq × τ$ pairs of queries and inverted lists independently, with tiling to reduce memory consumption.
A two-pass k-selection is introduced to reduce the number of partial results written back to global memory.

The use of product quantization codes and lookup tables enables efficient non-exhaustive search on GPUs.

By precomputing distances and storing them in lookup tables, the computational cost of distance calculations is greatly reduced.

The GPU implementation leverages the massive parallelism and fast shared memory to scan inverted lists and perform k-selection efficiently.

Multi-GPU parallelism

To further scale the similarity search system, the authors propose two strategies for multi-GPU parallelism:

Replication

If an index instance fits in the memory of a single GPU, it can be replicated across $R$ different GPUs.
Each replica handles a fraction $nq/R$ of the queries, and the results are joined back together on a single GPU or in CPU memory.
Replication provides near-linear speedup, except for potential efficiency loss with small nq.

Sharding

If an index instance does not fit in the memory of a single GPU, it can be sharded across $S$ different GPUs.
For adding vectors, each shard receives $n/S$ of the vectors.
For querying, each shard handles the full query set $nq$ , and the partial results are joined on a single GPU or in CPU memory, requiring an additional round of k-selection.
Sharding provides a speedup compared to replication, but is usually less than pure replication due to fixed overhead and the cost of subsequent k-selection.

Replication and sharding can be used together, with $S$ shards and each shard having $R$ replicas, for a total of $S × R$ GPUs. These strategies can also be applied to distribute an index across multiple machines.

In summary, the paper presents a comprehensive approach to efficient similarity search on GPUs, covering both exact and approximate search methods.

The key contributions include the WarpSelect algorithm for fast k-selection, the optimised computation layout for IVFADC indexing, and the multi-GPU parallelism strategies.

The experimental results demonstrate significant speedups compared to state-of-the-art GPU implementations and showcase the effectiveness of the proposed methods on large-scale datasets.

The open-source implementation provided by the authors enables researchers and practitioners to benefit from these advances in similarity search on GPUs.

The goal is to find the k-nearest neighbours $(k-NN)$ of a query vector $x$ within a large collection of vectors $[yi]$ based on the Euclidean distance.

IVFADC relies on two levels of quantization: a coarse quantizer $(q1)$ and a fine quantizer $(q2)$ that encodes the residual vector after the first level.

A product quantizer is used for the fine quantizer $(q2)$ to provide a large number of reproduction values without increasing the processing cost.

The product quantizer generates b-byte codes, where b is the number of sub-vectors, resulting in a large codebook size $(|C2| = 256^b)$ while keeping the quantization process computationally cheap.

In the context of IVFADC, exact search is used for the coarse quantizer $(q1).$

Optimised GEMM (General Matrix Multiplication) routines from the cuBLAS library are used to calculate the $-2<x_j, y_i>$ term for L2 distance.

The GEMM operation results in a partial distance matrix $D$ .

The fused kernel adds the $||y_i||^2$ term to each entry of the partial distance matrix $D$ and immediately submits the value to k-selection in registers.

The $||x_j||^2$ term is not required before k-selection and can be omitted.

For realistic problem sizes, the distance matrix $D$ may not fit in GPU memory.

To handle this, the problem is tiled over the batch of queries, with $t_q ≤ n_q$ queries being run in a single tile.

The effective memory requirement for $D$ becomes $O(2 × t_q)$ .

The computation can also be tiled over the database size $(n)$ if necessary.

When the query vector x and the coarse quantizer reproduction value $q1(y)$ are known, all distances to the PQ reproduction values can be precomputed and stored in lookup tables $T1, ..., Tb$ .

Computing the distance between $x$ and $q(y)$ consists of $b$ lookups and additions

When scanning the elements of an inverted list $IL$ , where $q1(y)$ is constant, the lookup table method can be applied.

Term 1 is independent of the query and can be precomputed and stored in a table of size $|C1| × 256 × b.$

Term 2 is the distance to $q1's$ reproduction value and is a by-product of the first-level quantizer $q1$ .

Term 3 can be computed independently of the inverted list, with a cost of $d × 256$ multiply-adds.

The inverted lists are stored as two separate arrays for $PQ$ codes and associated IDs.

A kernel is responsible for scanning the τ closest inverted lists for each query and calculating the per-vector pair distances using the lookup tables $Ti$ .

Multi-pass kernels are used to process the $nq × τ$ pairs of queries and inverted lists independently, with tiling to reduce memory consumption.

If an index instance fits in the memory of a single GPU, it can be replicated across $R$ different GPUs.

Each replica handles a fraction $nq/R$ of the queries, and the results are joined back together on a single GPU or in CPU memory.

If an index instance does not fit in the memory of a single GPU, it can be sharded across $S$ different GPUs.

For adding vectors, each shard receives $n/S$ of the vectors.

For querying, each shard handles the full query set $nq$ , and the partial results are joined on a single GPU or in CPU memory, requiring an additional round of k-selection.