Product Quantization

Product quantization (PQ) compresses a vector into a short code by cutting it into $m$ pieces and replacing each piece by the index of its nearest centroid in a small codebook of its own. A 128-dimensional SIFT descriptor, a VLAD or GIST vector, or a network embedding becomes, for example, eight bytes, and the distance between a query and millions of such codes can be estimated with a few table lookups per code, without decompressing them. Jégou, Douze and Schmid (2011) proposed it, as Ge et al. (2013) and the approximate nearest neighbor search article describe, and it is the compression stage of the inverted-file indexes described in approximate nearest neighbor search, which covers how PQ codes are searched and how their recall is measured. This article covers the quantizer itself: how it is trained, how much it distorts the data and why, how the choice of subspaces and a rotation before quantization (optimized product quantization) change that distortion, how biased the distances estimated from codes are, and the variants that relax its constraints: locally optimized, additive and anisotropic quantization. Other names for the same codes are PQ codes and compact codes; FAISS calls an index of them IndexPQ.

Problem

Large image collections are searched by comparing vectors: global descriptors and embeddings in image retrieval and image–text retrieval, aggregated descriptors such as VLAD, local descriptors in matching. A billion 128-dimensional float vectors take 512 GB, more than the memory of one machine, and comparing a query with all of them reads every byte. The problem is to store each vector in a few bytes such that distances to an uncompressed query can still be estimated accurately and quickly.

A quantizer replaces a vector $x$ by the nearest of a finite set of centroids, so storing $x$ costs only the centroid’s index. With $b$ bits per vector there are $2^b$ centroids. k-means learns good centroids, but for a 64-bit code it would need $2^{64}$ of them, which can be neither trained nor stored nor searched. Product quantization obtains a codebook of that size as a Cartesian product of small codebooks that can.

Inputs and Outputs

Training: a sample of $n$ vectors of dimension $D$ from the distribution to be compressed; the number of subquantizers $m$ (dividing $D$ ) and the number of centroids $k^*$ per subquantizer, usually $2^8 = 256$ . Output: $m$ codebooks of $k^*$ centroids of dimension $D/m$ , and optionally an orthogonal $D \times D$ rotation.
Encoding: a vector; output: a code of $m$ indices, $m \log_2 k^*$ bits (8 bytes for $m = 8$ , $k^* = 256$ ).
Decoding: a code; output: the reconstruction $q(x)$ , the concatenation of the $m$ selected centroids.
Distance estimation: a query and a set of codes; output: an approximate squared distance to each.

Intuition

Cutting a vector into $m$ subvectors and quantizing each with 256 centroids lets every combination of the $m$ choices be a reconstruction: $256^m$ of them from only $256\,m$ stored centroids. The price is a constraint. The reconstructions form a grid, the product of the subspace codebooks, so a combination of subvectors that never occurs in the data still has a centroid, and a dependency between subvectors cannot be represented. Locally optimized product quantization’s authors illustrate this with a 2D mixture, where many grid points of PQ fall where there are no data (Kalantidis & Avrithis, 2014, Fig. 1).

How the dimensions are split into subspaces therefore matters. If two strongly correlated coordinates fall into different subvectors, each is quantized as if the other did not exist; if one subvector receives most of the variance, its 256 centroids must cover a large spread while other subquantizers waste theirs on almost nothing. Optimized product quantization rotates the data before splitting it, so that subspaces are as independent as possible and carry balanced variance.

Algorithm

Training.

Optionally learn an orthogonal matrix $R$ (identity for plain PQ; see Optimized Product Quantization below) and replace every training vector $x$ by $Rx$ .
Split each training vector into $m$ consecutive subvectors of $D/m$ dimensions.
For each $j = 1, \dots, m$ , run k-means with $k^*$ centroids on the $j$ -th subvectors of all training vectors. FAISS’s ProductQuantizer trains each subquantizer with its k-means, whose defaults (Clustering.h, v1.15.1) are 25 iterations and at most 256 training points per centroid, warning below 39.

Encoding. Rotate the vector, split it, and record for each subvector the index of its nearest centroid. Because the codebook is a product, the nearest product centroid is the concatenation of the nearest subcentroids (Ge et al., 2013, §2.2), so encoding costs $k^* D$ multiply–adds rather than a search over $k^{*m}$ codewords.

Distance estimation. For a query $x$ , compute the $m \times k^*$ table of squared distances between the query’s subvectors and every subcentroid. The asymmetric distance to a code is then the sum of $m$ table entries; the symmetric distance quantizes the query too and reads a precomputed table of distances between subcentroids. The approximate nearest neighbor search article describes both as search methods, with the inverted file (IVFADC) that encodes residuals from a coarse quantizer, and, as that article reports, Jégou et al. found the asymmetric distance clearly more accurate; Ge et al. (2013) likewise note that it is more accurate at the same cost.

Mathematical Formulation

The quantizer and its objective

Write $x \in \mathbb{R}^D$ as $m$ subvectors $u_1(x), \dots, u_m(x)$ of dimension $D/m$ , and let $\mathcal{C}_j$ be the $j$ -th codebook, with $q_j$ mapping a subvector to its nearest centroid in $\mathcal{C}_j$ . The product quantizer is

q(x) = \big(q_1(u_1(x)), \dots, q_m(u_m(x))\big), \qquad q(x) \in \mathcal{C} = \mathcal{C}_1 \times \cdots \times \mathcal{C}_m ,

and $|\mathcal{C}| = k^{*m}$ . Ge, He, Ke and Sun (2013, Eq. 2) pointed out that PQ minimizes the same objective as k-means, the quantization distortion

E = \frac{1}{n} \sum_{x} \lVert x - q(x) \rVert^2 ,

with the codebook constrained to a Cartesian product; they note that the original paper did not state this objective explicitly. Because $\lVert x - q(x) \rVert^2 = \sum_j \lVert u_j(x) - q_j(u_j(x)) \rVert^2$ , the problem splits into $m$ independent k-means problems. Ge et al. also show that the hashing method ITQ is a vector quantizer of the same kind, with the codebook constrained to the vertices of a rotated hypercube, so that one distortion compares k-means, PQ and binary codes (§2.2–2.3). Across these methods, at 16 bits on SIFT1M and GIST1M, they observed that the distortion was tightly correlated with the mean average precision of nearest-neighbour search (their Fig. 1), and they used it as the measure of a quantizer’s quality.

Distance estimates and their bias

The asymmetric (ADC) and symmetric (SDC) estimates of the squared distance between a query $x$ and a database vector $y$ are

\tilde d_{\mathrm{ADC}}(x, y)^2 = \sum_{j=1}^{m} \lVert u_j(x) - q_j(u_j(y)) \rVert^2, \qquad \tilde d_{\mathrm{SDC}}(x, y)^2 = \sum_{j=1}^{m} \lVert q_j(u_j(x)) - q_j(u_j(y)) \rVert^2 .

Both are biased low on average, by a short argument (ours). If a centroid $c$ is the mean of the vectors $y$ in its cell, which is what a k-means update makes it on the training set, then averaging over those vectors gives

\overline{\lVert x - y \rVert^2} = \lVert x - c \rVert^2 + \overline{\lVert y - c \rVert^2},

since the cross term $-2 (x - c)^\top \overline{(y - c)}$ vanishes. The asymmetric estimate therefore falls short of the true squared distance by the database vector’s quantization error on average, and the symmetric estimate, applying the same argument to the query’s cell as well, by the errors of both vectors. Adding the mean squared error of each cell, stored once per centroid, removes this average bias on the training set. On new data a centroid is not exactly the mean of its cell, and the example below measures what remains. For ranking the database by distance to one query, a shift that is the same for all database vectors changes nothing; what matters is the part of the error that varies from vector to vector.

A Gaussian bound on the distortion

Ge et al. (2013, §3.2) analyse PQ for data drawn from $\mathcal{N}(0, \Sigma)$ . From rate–distortion theory, a quantizer with $k$ codewords of equal length cannot have a distortion below

E \ge k^{-2/D}\, D\, |\Sigma|^{1/D} \qquad \text{(their Eq. 9)},

an expression they derive under the assumption that every variance is sufficiently large, and point to Cover and Thomas for the general form. They consider it an approximation to k-means distortion: with $10^5$ samples, $k = 256$ and variances drawn in $[0.5, 1]$ , k-means came within about 5% of the bound for $D = 32$ , 64 and 128 (k-means 17.1, 39.9 and 88.5 against bounds of 16.2, 38.8 and 86.7). The general form is the distortion-rate function obtained by reverse water-filling over the eigenvalues of $\Sigma$ , which the example below computes as the reference that no quantizer of Gaussian data at that code length can beat. Applying Eq. 9 to each subspace, with $\Sigma_{jj}$ the covariance of the $j$ -th subvector, gives the bound for a product quantizer with $k^*$ centroids per subspace:

E_{\mathrm{PQ}} = k^{*\,-2m/D}\, \frac{D}{m} \sum_{j=1}^{m} |\Sigma_{jj}|^{m/D} \qquad \text{(their Eq. 11)} .

At a fixed code length $b = m \log_2 k^*$ the factor $k^{*\,-2m/D} = 2^{-2b/D}$ does not depend on $m$ ; the split enters only through the determinants of the diagonal blocks of $\Sigma$ .

Optimized Product Quantization

The problem. Jégou et al. had found that PQ’s accuracy depends on the order of the dimensions, as Ge et al. report: SIFT works best with its natural order, where each subvector holds neighbouring histograms, and GIST with a “structural” order. For unknown data, the dimensions were ordered at random, or the data were aligned by PCA and then randomly rotated, which Ge et al. credit to Jégou, Douze, Schmid and Pérez (2010) (Ge et al., 2013, §1 and §4). Ge et al. made the space decomposition a parameter: an orthogonal matrix $R$ applied before splitting, chosen with the codebooks to minimize distortion,

\min_{R,\, \mathcal{C}_1, \dots, \mathcal{C}_m} \sum_x \lVert R x - \hat c(R x) \rVert^2 \quad \text{s.t.} \quad \hat c \in \mathcal{C}_1 \times \cdots \times \mathcal{C}_m,\ R^\top R = I ,

which they called optimized product quantization (OPQ; their Eq. 4). A permutation of the dimensions is one such $R$ . Norouzi and Fleet (2013) developed a closely related generalization of PQ, Cartesian k-means, at the same conference, and each paper notes the other as concurrent work.

Non-parametric solution (OPQ $_\mathrm{NP}$ ). Alternate between the codebooks and the rotation (Ge et al., Algorithm 1). With $R$ fixed, the problem is PQ on $Rx$ . With the codebooks and codes fixed, let $y$ be each vector’s reconstruction; minimizing $\sum \lVert R x - y \rVert^2$ over orthogonal $R$ is the orthogonal Procrustes problem, solved by the singular value decomposition $X Y^\top = U S V^\top$ and $R = V U^\top$ (the problem the Kabsch algorithm solves when it aligns point sets, there restricted to proper rotations). The authors found that one k-means iteration per alternation sufficed for convergence, that 100 iterations were good enough for search, and that the result is a local optimum that depends on the initialization (§3.1). FAISS implements this method as OPQMatrix; in v1.15.1 it starts from a random rotation and runs 50 outer iterations, with 40 k-means iterations in the first and 4 in each later one.

Parametric solution (OPQ $_\mathrm{P}$ ). For Gaussian data the rotation can be chosen from the covariance alone, by minimizing Eq. 11 over $R$ , that is, $\sum_j |\hat\Sigma_{jj}|^{m/D}$ with $\hat\Sigma = R \Sigma R^\top$ (their Eq. 12). By the inequality of arithmetic and geometric means and by Fischer’s inequality $\prod_j |\hat\Sigma_{jj}| \ge |\hat\Sigma| = |\Sigma|$ ,

\sum_{j=1}^{m} |\hat\Sigma_{jj}|^{m/D} \;\ge\; m \prod_{j=1}^{m} |\hat\Sigma_{jj}|^{1/D} \;\ge\; m\, |\Sigma|^{1/D} \qquad \text{(their Eqs. 13–15)},

with equality when two criteria hold: independence, the subspaces are uncorrelated, which PCA alignment achieves; and balance, every subspace has the same product of eigenvalues. Ge et al. called the resulting procedure eigenvalue allocation: align the data by PCA without reducing the dimension, sort the eigenvalues in decreasing order, and assign each in turn to the bucket, of $m$ buckets with room for $D/m$ eigenvalues each, whose product of eigenvalues is smallest; each bucket’s principal directions form one subspace. On SIFT and GIST covariances with $m = 8$ and $k^* = 256$ they reported that this greedy rule reached the theoretical minimum to four or five digits (for SIFT an objective of $2.9287 \times 10^3$ against the minimum $2.9286 \times 10^3$ ; for GIST $1.9870 \times 10^{-3}$ for both). Kalantidis and Avrithis (2014, Eq. 12) state the same rule. Both criteria explain earlier practice in Ge et al.’s reading: PCA projections give independence, and random rotations or adaptive bit allocation give balance, but a random rotation after PCA gives up independence.

How the greedy rule depends on scale (our reasoning). The rule as written in the text compares products of eigenvalues, and multiplying by an eigenvalue smaller than 1 makes a product smaller. Once a bucket’s product falls below the others, it remains the smallest after it receives the next eigenvalue if that eigenvalue is below 1, so it keeps receiving them until it is full, and the buckets end up unbalanced rather than balanced. Rescaling the data, for instance measuring intensities in $[0, 255]$ rather than $[0, 1]$ , multiplies a bucket’s product by a power of the scale that depends on how many eigenvalues it holds, so the allocation can change. The sources leave such details open: the papers state the rule without discussing eigenvalues below 1, and FAISS’s PCAMatrix option balanced_bins uses the same greedy on sums of eigenvalues instead. The example therefore prints three allocation rules: the product rule as written, the sum rule, and a search that swaps eigenvalues between buckets as long as Eq. 12 decreases.

Results reported by Ge et al. On 128-dimensional Gaussian data with eigenvalues $e^{-0.1 d}$ , OPQ $_\mathrm{NP}$ and OPQ $_\mathrm{P}$ performed almost the same, and OPQ $_\mathrm{P}$ reached the theoretical minimum of Eq. 15; random ordering and PCA with a random rotation were substantially worse. On SIFT1M, GIST1M and MNIST, with $k^* = 256$ and $m = b/8$ (their Figs. 4–6), both OPQ variants substantially outperformed random ordering and random rotation for symmetric and asymmetric distances; the non-parametric solution improved further on SIFT1M and MNIST, which the authors attribute to clusters in those data, and the two were comparable on GIST1M. Without any prior knowledge of the descriptor layout, OPQ $_\mathrm{NP}$ also outperformed PQ with the recommended natural or structural order at 64 bits with symmetric distances (their Fig. 7).

Parameters

Code length $b = m \log_2 k^*$ bits. It fixes memory per vector and, at a fixed $b$ , the theoretical distortion factor $2^{-2b/D}$ . FAISS stores a code in $\lceil m \cdot \mathrm{nbits} / 8 \rceil$ bytes and refuses more than 24 bits per subquantizer.
Centroids per subquantizer $k^*$ . With $k^* = 256$ each index is one byte and the lookup tables stay small; the experiments of OPQ and LOPQ use 256. Encoding and table costs grow in proportion to $k^*$ (see Complexity), and each centroid needs enough training vectors.
Subquantizers $m$ at a fixed $b$ . Fewer subquantizers with more centroids each keep more of the dependencies between dimensions; the example measures this trade-off.
Rotation. Identity, a known good order, or a learned rotation (OPQ). For data with a single dominant cluster structure the parametric solution is cheap; the non-parametric one costs a few PQ trainings.
Training set. It should come from the distribution being encoded. Distortion on held-out vectors is higher than on the training vectors (E2 and E3 in the example); FAISS subsamples to at most 256 training points per centroid.

Complexity

Memory: $b$ bits per vector, plus $k^* D$ floats for the codebooks and $D^2$ for a rotation.
Training: each k-means iteration costs $n k^* D$ multiply–adds over all subquantizers, the same as one k-means with $k^*$ centroids in $D$ dimensions. OPQ $_\mathrm{NP}$ adds a $D \times D$ SVD and a rotation of the training set per iteration.
Encoding: $k^* D$ multiply–adds per vector, plus $D^2$ for a rotation.
Asymmetric search: $k^* D$ operations per query for the table, then $m$ lookups and additions per code, the cost that makes exhaustive search over a million codes take milliseconds (Ge et al. quote about 20 ms per million 64-bit codes on their 2013 CPU).

Implementation

The script below measures three things about the quantizer on data whose statistics are known: how distortion depends on the split at a fixed code length, what a rotation before quantization gains, and how far distance estimates from codes are from the true distances. It measures distortion and distance error, not recall, which the approximate nearest neighbor search example measures for 8-byte PQ codes of SIFT-like descriptors.

Data A are $8 \times 8$ patches of synthetic images with an exactly known correlation structure: white Gaussian noise convolved with a unit-norm anisotropic Gaussian kernel (standard deviation 3 pixels along a direction at 30°, 1 pixel across it), plus independent pixel noise of standard deviation 0.1. Each patch is flattened row by row into a 64-dimensional vector; its covariance is the kernel’s autocorrelation at the pixel offsets plus $0.01 I$ , with trace 64.64, and the vectors are Gaussian with mean zero. Training patches come from four $256 \times 256$ images (20,000 patches), test patches from two others (10,000). (The principal component analysis article whitens similar patches of a smoothed-noise texture.) Data B are $8 \times 8$ patches of scikit-image’s bundled camera photograph, with intensities in $[0, 1]$ , centred by the training mean: 20,000 training patches from the even-numbered 64-row bands of the image and 10,000 test patches from the odd-numbered bands, so the two sets share no pixels. These data are not Gaussian.

The k-means is Lloyd’s algorithm started from $k^*$ random training vectors, run for 25 iterations (FAISS’s default), and checked against scikit-learn’s. Distortion is the mean squared error per vector summed over the 64 dimensions, on the test patches unless marked as training distortion. Seeds 0, 1 and 2 change the k-means initialization, the random permutation and rotation, and the swap search; tables give the mean and the range over the three seeds.

E1 quantizes data A in natural pixel order at 32 and 64 bits with every split of the form $k^* = 2^{b/m}$ , from 256 centroids per subspace down to 2, and prints Eq. 11 for that split and the Gaussian distortion-rate function $D(R)$ at the code length.
E2 compares rotations at $m = 8$ , $k^* = 256$ (64 bits) on both data sets: the natural order; a random permutation of the dimensions; PCA with the eigenvalues sorted, so that consecutive principal components share a subspace; PCA followed by a random rotation; the three eigenvalue-allocation rules described above; and OPQ $_\mathrm{NP}$ started from the swap search, with 50 outer iterations, each a Procrustes step and one k-means step per subspace. Rotations are computed from the training covariance; Eq. 11 is evaluated with the true covariance for data A and with the training covariance for data B.
E3 encodes data A in natural order ( $m = 8$ , $k^* = 256$ , seed 0), uses 1,000 test patches as queries and the other 9,000 as the database, and measures the error of the four distance estimates: ADC, SDC, and each with the mean training error of the cells added.

The checks printed first confirm the set-up: the sample covariances of data A against the true one, the total variance of data B’s two sets, agreement with scikit-learn’s k-means, the swap search reaching the Eq. 15 minimum, every rotation being orthogonal, OPQ $_\mathrm{NP}$ ’s objective never rising, and the averaging identity of the previous section on the training set. The script uses NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2 and scikit-learn 1.7.2, runs in about four minutes on one thread, and two runs printed identical output.

import os
os.environ["OMP_NUM_THREADS"] = os.environ["OPENBLAS_NUM_THREADS"] = "1"   # before NumPy loads
import numpy as np
from scipy.signal import fftconvolve
from skimage import data as skdata
from sklearn.cluster import KMeans

P, D = 8, 64                       # 8x8 patches, flattened row by row into 64-D vectors
SEEDS = (0, 1, 2)                  # seeds for k-means initialization and random rotations

# ---------- data A: patches of a Gaussian random field whose covariance is known exactly ----------
def kernel(K=25, s_major=3.0, s_minor=1.0, theta=np.deg2rad(30)):
    y, x = np.mgrid[-(K // 2):K // 2 + 1, -(K // 2):K // 2 + 1].astype(float)
    u, v = x * np.cos(theta) + y * np.sin(theta), -x * np.sin(theta) + y * np.cos(theta)
    h = np.exp(-(u ** 2 / (2 * s_major ** 2) + v ** 2 / (2 * s_minor ** 2)))
    return h / np.linalg.norm(h)                    # unit-norm kernel: field variance 1

SIGMA_N = 0.1                                       # white pixel noise on top of the field

def true_cov(h):
    """Exact covariance of a patch: the autocorrelation of h at every pixel offset, plus noise."""
    ac = fftconvolve(h, h[::-1, ::-1], mode="full")
    c = ac.shape[0] // 2
    r, q = np.divmod(np.arange(D), P)
    return ac[c + r[:, None] - r[None, :], c + q[:, None] - q[None, :]] + SIGMA_N ** 2 * np.eye(D)

def data_A(seed=2026, n_train=20000, n_test=10000):
    rng, h = np.random.default_rng(seed), kernel()
    def patches(n, size=256):                       # one new 256x256 image, n random patches
        img = fftconvolve(rng.standard_normal((size + 24, size + 24)), h, mode="valid")
        img += SIGMA_N * rng.standard_normal(img.shape)
        r, c = rng.integers(0, size - P + 1, (2, n))
        return np.stack([img[i:i + P, j:j + P].ravel() for i, j in zip(r, c)])
    train = np.concatenate([patches(n_train // 4) for _ in range(4)])
    test = np.concatenate([patches(n_test // 2) for _ in range(2)])
    return train, test, true_cov(h)

# ---------- data B: patches of the bundled camera image (train and test from alternate bands) ----------
def data_B(seed=2027, n_train=20000, n_test=10000):
    rng, img = np.random.default_rng(seed), skdata.camera() / 255.0
    def patches(bands, n):                          # patches lie inside 64-row bands
        pos = [(64 * b + i, j) for b in bands for i in range(64 - P + 1) for j in range(512 - P + 1)]
        return np.stack([img[pos[k][0]:pos[k][0] + P, pos[k][1]:pos[k][1] + P].ravel()
                         for k in rng.choice(len(pos), n, replace=False)])
    train, test = patches([0, 2, 4, 6], n_train), patches([1, 3, 5, 7], n_test)
    mu = train.mean(0)
    return train - mu, test - mu

# ---------- k-means and product quantization (Jegou, Douze & Schmid 2011) ----------
def sqdist(X, C):
    return (X ** 2).sum(1)[:, None] - 2 * X @ C.T + (C ** 2).sum(1)[None, :]

def update(X, a, C):
    """Move each centroid to the mean of its points (empty clusters keep their centroid)."""
    n = np.bincount(a, minlength=len(C))
    S = np.stack([np.bincount(a, X[:, i], len(C)) for i in range(X.shape[1])], 1)
    C[n > 0] = S[n > 0] / n[n > 0, None]
    return C

def kmeans(X, k, rng, iters=25, C=None):
    """Lloyd's algorithm from k random training points, 25 iterations (FAISS's default)."""
    C = X[rng.choice(len(X), k, replace=False)].copy() if C is None else C
    for _ in range(iters):
        a = np.argmin(sqdist(X, C), 1)
        C = update(X, a, C)
    return C, a                                      # C holds the means of assignment a

def pq_train(X, m, k, rng):
    d = X.shape[1] // m
    out = [kmeans(X[:, j * d:(j + 1) * d], k, rng) for j in range(m)]
    return [C for C, _ in out], np.stack([a for _, a in out], 1)

def encode(X, books):
    d = books[0].shape[1]
    return np.stack([np.argmin(sqdist(X[:, j * d:(j + 1) * d], B), 1) for j, B in enumerate(books)], 1)

def decode(codes, books):
    return np.concatenate([B[codes[:, j]] for j, B in enumerate(books)], 1)

def distortion(X, books):
    """Mean squared error per vector, summed over the 64 dimensions."""
    return ((X - decode(encode(X, books), books)) ** 2).sum(1).mean()

# ---------- Gaussian reference values ----------
def bound_pq(S, m, k):
    """Ge et al. (2013), Eq. (11): k^(-2m/D) (D/m) sum_j |S_jj|^(m/D) for the blocks of S."""
    d = D // m
    logdet = np.array([np.linalg.slogdet(S[j * d:(j + 1) * d, j * d:(j + 1) * d])[1] for j in range(m)])
    return k ** (-2 * m / D) * (D / m) * np.exp(logdet * m / D).sum()

def shannon(S, bits):
    """Distortion-rate function of N(0, S) by reverse water-filling: no quantizer with
    2^bits codewords can have a lower expected distortion."""
    lam, lo, hi = np.linalg.eigvalsh(S), 1e-12, np.linalg.eigvalsh(S).max()
    for _ in range(200):
        th = np.sqrt(lo * hi)
        lo, hi = (th, hi) if 0.5 * np.log2(np.maximum(lam / th, 1)).sum() > bits else (lo, th)
    return np.minimum(th, lam).sum()

# ---------- rotations R (the encoder quantizes R x) ----------
def allocate(lam, m, rule, rng=None):
    """Group the principal axes into m subspaces of D/m axes. 'product': Ge et al.'s greedy as
    written (next eigenvalue to the non-full bucket with the smallest product); 'sum': the same
    greedy on sums (FAISS PCAMatrix.balanced_bins); 'search': pairwise swaps that minimize Eq. (12)."""
    d = D // m
    if rule == "search":
        order = allocate(lam, m, "sum")
        obj = lambda o: np.exp(np.log(lam[o]).reshape(m, d).sum(1) * m / D).sum()
        best = obj(order)
        for _ in range(20000):
            i, j = rng.integers(0, D, 2)
            if i // d != j // d:
                o = order.copy(); o[i], o[j] = o[j], o[i]
                if obj(o) < best:
                    order, best = o, obj(o)
        return order
    buckets = [[] for _ in range(m)]
    for i in range(D):
        score = [np.inf if len(b) == d else (np.prod(lam[b]) if rule == "product" else lam[b].sum())
                 for b in buckets]
        buckets[int(np.argmin(score))].append(i)
    return [i for b in buckets for i in b]

def rotations(S, m, rng):
    lam, V = np.linalg.eigh(S)
    lam, V = lam[::-1], V[:, ::-1]                    # principal axes, largest variance first
    Q, T = np.linalg.qr(rng.standard_normal((D, D)))
    Q *= np.sign(np.diag(T))                          # a random orthogonal matrix
    return {"natural order": np.eye(D),
            "random permutation": np.eye(D)[rng.permutation(D)],
            "PCA, sorted": V.T,
            "PCA + random rotation": Q @ V.T,
            "EA greedy, product": V[:, allocate(lam, m, "product")].T,
            "EA greedy, sum": V[:, allocate(lam, m, "sum")].T,
            "EA swap search": V[:, allocate(lam, m, "search", rng)].T}

def opq_np(X, R, m, k, rng, outer=50):
    """Ge et al., Algorithm 1, initialized with R: PQ on RX, then each iteration solves the
    orthogonal Procrustes problem for R and runs one k-means step in every subspace."""
    d = D // m
    books, codes = pq_train(X @ R.T, m, k, rng)
    trace = []
    for _ in range(outer):
        U, _, Vt = np.linalg.svd(X.T @ decode(codes, books))
        R = Vt.T @ U.T                               # argmin_R sum ||R x - y||^2, R^T R = I
        Xr = X @ R.T
        for j in range(m):
            Z = Xr[:, j * d:(j + 1) * d]
            books[j] = update(Z, codes[:, j], books[j])
            codes[:, j] = np.argmin(sqdist(Z, books[j]), 1)
        trace.append(((Xr - decode(codes, books)) ** 2).sum(1).mean())
    return R, books, np.array(trace)

def stats(v):
    v = np.array(v)
    return f"{v.mean():8.4f} [{v.min():.4f}-{v.max():.4f}]"

# ---------- checks ----------
trA, teA, S_true = data_A()
trB, teB = data_B()
S_A, S_B = trA.T @ trA / len(trA), trB.T @ trB / len(trB)
rel = lambda S: np.linalg.norm(S - S_true) / np.linalg.norm(S_true)
print(f"data A: |S - S_true|/|S_true| = {rel(S_A):.4f} (train), {rel(teA.T @ teA / len(teA)):.4f} (test); "
      f"trace S_true = {np.trace(S_true):.2f}")
print(f"data B: total variance {np.trace(S_B):.4f} (train), {(teB ** 2).sum(1).mean():.4f} (test)")
X0 = trA[:, :8]
C0 = X0[np.random.default_rng(5).choice(len(X0), 256, replace=False)]
C, a = kmeans(X0, 256, None, C=C0.copy())
km = KMeans(256, init=C0, n_init=1, max_iter=25, tol=0, algorithm="lloyd").fit(X0)
print(f"k-means vs scikit-learn (8-D, k=256, 25 iterations): max centroid difference "
      f"{np.abs(C - km.cluster_centers_).max():.1e}")

# ---------- E1: subquantizers m against centroids k at a fixed code length (data A, natural order) ----------
print(f"\nE1 data A, natural pixel order: test distortion, mean [min-max] over seeds {SEEDS[0]}-{SEEDS[-1]}")
print(f"{'bits':>4} {'m':>3} {'k':>4} {'dims/sub':>8} {'distortion':>26} {'Eq. 11':>8} {'Shannon':>8}")
for bits, ms in ((32, (4, 8, 16, 32)), (64, (8, 16, 32, 64))):
    for m in ms:
        k = 2 ** (bits // m)
        e = [distortion(teA, pq_train(trA, m, k, np.random.default_rng(s))[0]) for s in SEEDS]
        print(f"{bits:4d} {m:3d} {k:4d} {D // m:8d} {stats(e):>26} {bound_pq(S_true, m, k):8.3f} "
              f"{shannon(S_true, bits):8.3f}")

# ---------- E2: the rotation before quantization (m = 8, k = 256, 64 bits) ----------
worst_orth, worst_rise = 0.0, -np.inf
for name, tr, te, S_ref in (("A (Gaussian field)", trA, teA, S_true), ("B (camera)", trB, teB, S_B)):
    S_tr = tr.T @ tr / len(tr)
    R = rotations(S_tr, 8, np.random.default_rng(0))["EA swap search"]
    eq15 = D * 256 ** (-16 / D) * np.exp(np.linalg.slogdet(S_tr)[1] / D)
    print(f"\ndata {name}: EA swap search on the training covariance, Eq. 11 / Eq. 15 - 1 = "
          f"{bound_pq(R @ S_tr @ R.T, 8, 256) / eq15 - 1:.1e}")
    rows = {}
    for s in SEEDS:
        rng = np.random.default_rng(s)
        Rs = rotations(S_tr, 8, rng)
        for r, R in Rs.items():
            worst_orth = max(worst_orth, np.abs(R @ R.T - np.eye(D)).max())
            books = pq_train(tr @ R.T, 8, 256, rng)[0]
            rows.setdefault(r, []).append((distortion(te @ R.T, books), distortion(tr @ R.T, books),
                                           bound_pq(R @ S_ref @ R.T, 8, 256)))
        R, books, trace = opq_np(tr, Rs["EA swap search"], 8, 256, rng)
        worst_rise = max(worst_rise, np.max(np.diff(trace) / trace[:-1]))
        worst_orth = max(worst_orth, np.abs(R @ R.T - np.eye(D)).max())
        rows.setdefault("OPQ_NP (from EA swap search)", []).append(
            (distortion(te @ R.T, books), distortion(tr @ R.T, books), bound_pq(R @ S_ref @ R.T, 8, 256)))
    print(f"E2 data {name}, m = 8, k = 256: distortion, mean [min-max] over seeds {SEEDS[0]}-{SEEDS[-1]}")
    print(f"{'rotation':32s} {'test distortion':>26} {'train':>9} {'Eq. 11':>9}")
    for r, v in rows.items():
        v = np.array(v)
        print(f"{r:32s} {stats(v[:, 0]):>26} {v[:, 1].mean():9.4f} {v[:, 2].mean():9.4f}")
    print(f"{'Eq. 15 (minimum over R)':32s} {'':36s} {D * 256 ** (-16 / D) * np.exp(np.linalg.slogdet(S_ref)[1] / D):9.4f}")
    print(f"{'Shannon D(R), 64 bits':32s} {'':36s} {shannon(S_ref, 64):9.4f}")
print(f"\nlargest |R R^T - I| = {worst_orth:.1e}; largest relative rise of the OPQ_NP objective = {worst_rise:.1e}")

# ---------- E3: asymmetric against symmetric distances (data A, natural order, m=8, k=256, seed 0) ----------
books, codes = pq_train(trA, 8, 256, np.random.default_rng(0))
E_train = ((trA - decode(codes, books)) ** 2).sum(1).mean()
cell_err = [np.bincount(codes[:, j], ((trA[:, 8 * j:8 * j + 8] - B[codes[:, j]]) ** 2).sum(1), 256)
            / np.maximum(np.bincount(codes[:, j], minlength=256), 1) for j, B in enumerate(books)]
Q, X = teA[:1000], teA[1000:]
true_tr = np.maximum(sqdist(Q, trA), 0)
adc_tr = sqdist(Q, decode(codes, books))
print(f"\nE3 training set as database: E_train = {E_train:.4f}; "
      f"max over queries |mean(d^2 - ADC) - E_train| / E_train = "
      f"{np.abs((true_tr - adc_tr).mean(1) - E_train).max() / E_train:.1e}")
cq, cx = encode(Q, books), encode(X, books)
yq, yx = decode(cq, books), decode(cx, books)
corr = lambda c: sum(e[c[:, j]] for j, e in enumerate(cell_err))
true = np.maximum(sqdist(Q, X), 0)
E_q, E_x = ((Q - yq) ** 2).sum(1).mean(), ((X - yx) ** 2).sum(1).mean()
est = {"ADC": sqdist(Q, yx), "SDC": sqdist(yq, yx)}
est["ADC + cell error"] = est["ADC"] + corr(cx)[None, :]
est["SDC + cell errors"] = est["SDC"] + corr(cq)[:, None] + corr(cx)[None, :]
print(f"1,000 queries x 9,000 test vectors; mean true squared distance {true.mean():.3f}; "
      f"E_query = {E_q:.4f}, E_database = {E_x:.4f}")
print(f"{'estimate':18s} {'bias':>8} {'bias/E_db':>9} {'error sd':>9}")
for n, v in est.items():
    err = v - true
    print(f"{n:18s} {err.mean():8.4f} {err.mean() / E_x:9.3f} {err.std():9.4f}")

data A: |S - S_true|/|S_true| = 0.0250 (train), 0.0228 (test); trace S_true = 64.64
data B: total variance 5.4156 (train), 5.4250 (test)
k-means vs scikit-learn (8-D, k=256, 25 iterations): max centroid difference 2.2e-15

E1 data A, natural pixel order: test distortion, mean [min-max] over seeds 0-2
bits   m    k dims/sub                 distortion   Eq. 11  Shannon
  32   4  256       16     6.6629 [6.6494-6.6763]    2.559    2.317
  32   8   16        8  13.1051 [13.0828-13.1301]    3.906    2.317
  32  16    4        4  17.6683 [17.5928-17.7576]    6.039    2.317
  32  32    2        2  25.4962 [25.4962-25.4963]   13.334    2.317
  64   8  256        8     3.9113 [3.9052-3.9169]    1.953    0.714
  64  16   16        4     6.0638 [6.0544-6.0686]    3.020    0.714
  64  32    4        2  10.2567 [10.2547-10.2593]    6.667    0.714
  64  64    2        1  23.6334 [23.6334-23.6334]   16.160    0.714

data A (Gaussian field): EA swap search on the training covariance, Eq. 11 / Eq. 15 - 1 = 1.3e-04
E2 data A (Gaussian field), m = 8, k = 256: distortion, mean [min-max] over seeds 0-2
rotation                                    test distortion     train    Eq. 11
natural order                        3.9179 [3.9082-3.9279]    3.4964    1.9529
random permutation                10.7068 [10.0653-11.0447]    9.5614    7.2188
PCA, sorted                       15.3702 [15.3300-15.4191]   13.5812   11.5528
PCA + random rotation             10.4500 [10.1268-10.7313]    9.3664    7.5656
EA greedy, product                   2.0554 [2.0538-2.0571]    1.8677    1.5704
EA greedy, sum                       1.4509 [1.4475-1.4526]    1.3198    1.0409
EA swap search                       0.9927 [0.9680-1.0266]    0.9048    0.7120
OPQ_NP (from EA swap search)         0.9929 [0.9679-1.0191]    0.8961    0.7130
Eq. 15 (minimum over R)                                                  0.7108
Shannon D(R), 64 bits                                                    0.7139

data B (camera): EA swap search on the training covariance, Eq. 11 / Eq. 15 - 1 = 4.2e-07
E2 data B (camera), m = 8, k = 256: distortion, mean [min-max] over seeds 0-2
rotation                                    test distortion     train    Eq. 11
natural order                        0.1099 [0.1093-0.1107]    0.0547    0.0702
random permutation                   0.1463 [0.1444-0.1480]    0.0770    0.1055
PCA, sorted                          0.1597 [0.1590-0.1603]    0.0834    0.1190
PCA + random rotation                0.1535 [0.1517-0.1562]    0.0837    0.1045
EA greedy, product                   0.1064 [0.1059-0.1070]    0.0606    0.0616
EA greedy, sum                       0.0591 [0.0590-0.0592]    0.0319    0.0266
EA swap search                       0.0569 [0.0566-0.0571]    0.0315    0.0240
OPQ_NP (from EA swap search)         0.0547 [0.0541-0.0555]    0.0287    0.0242
Eq. 15 (minimum over R)                                                  0.0240
Shannon D(R), 64 bits                                                    0.0241

largest |R R^T - I| = 4.8e-15; largest relative rise of the OPQ_NP objective = -3.1e-05

E3 training set as database: E_train = 3.4959; max over queries |mean(d^2 - ADC) - E_train| / E_train = 8.9e-16
1,000 queries x 9,000 test vectors; mean true squared distance 131.014; E_query = 3.9476, E_database = 3.9005
estimate               bias bias/E_db  error sd
ADC                 -4.6158    -1.183   10.6020
SDC                 -9.2274    -2.366   15.0329
ADC + cell error    -1.1264    -0.289   10.2948
SDC + cell errors   -2.2275    -0.571   14.6273

The checks. The sample covariances of data A are within 2.5% of the exact one, data B’s training and test sets have nearly equal total variance, and the hand-written k-means agrees with scikit-learn’s to rounding error. The swap search reaches the Eq. 15 minimum to within $1.3 \times 10^{-4}$ (data A) and $4.2 \times 10^{-7}$ (data B), every rotation is orthogonal, and the OPQ $_\mathrm{NP}$ objective decreased at every iteration, as each of its steps minimizes the same objective over one block of parameters. With the training set as database, the mean gap between true and asymmetric squared distances equals the training distortion for every query to $8.9 \times 10^{-16}$ , the identity derived above.

E1: fewer subquantizers with more centroids. At both code lengths the distortion grows with $m$ : from 6.66 with four 16-dimensional subquantizers of 256 centroids to 25.50 with 32 subquantizers of two centroids at 32 bits, and from 3.91 to 23.63 at 64 bits, with seed ranges that do not overlap. The hypothesis written before the run, that the distortion increases with $m$ at a fixed code length, held at both lengths. Eq. 11 for the natural split rises with $m$ too (in our reading, because each further split ignores more correlation between pixels), but the measured distortions stay far above it: at 64 bits with $m = 8$ , 3.91 against 1.95. Every measured value is above the distortion-rate function (2.317 at 32 bits, 0.714 at 64 bits), as it must be, which the plan recorded as a check. In our reading, the gap between Eq. 11 and the measurements has two sources. Eq. 11 inherits Eq. 9’s assumption that every variance is sufficiently large, which we read as large compared with the distortion per dimension, but at 64 bits the distortion-rate function allows about $0.714 / 64 \approx 0.011$ per dimension, about the pixel-noise variance of 0.01 that, by construction, every direction of data A has at least; and a small quantizer, such as one with two centroids per pixel, is far from the bound’s ideal.

E2: the rotation decides the distortion. On data A the test distortion ranges from 15.37 with sorted PCA components, where the first subquantizer receives the eight largest eigenvalues and the last the eight smallest, to 0.99 after the swap search: the same code length and the same k-means, with a fifteen-fold difference set by $R$ alone. The natural pixel order (3.92) is better than a random permutation (10.71) or PCA plus a random rotation (10.45); in our reading, because each row of an $8 \times 8$ patch keeps neighbouring, strongly correlated pixels together. The natural-order row differs slightly from E1’s row with the same $m$ and $k^*$ (3.9179 against 3.9113) because in E2 the random generator first draws the permutation, rotation and swap search, so k-means starts from different points.

The plan’s hypothesis that the measured order of the seven fixed rotations on data A would follow their order under Eq. 11 failed: the order held for five of the six neighbouring pairs, with seed ranges that do not overlap, but the random permutation (mean 10.71, range 10.07–11.04) was worse than PCA plus a random rotation (10.45, range 10.13–10.73), while their mean Eq. 11 values are in the opposite order (7.22 against 7.57); the seed ranges of these two overlap. Among the eigenvalue allocations, the product rule as written gives 2.06 on data A, the sum rule 1.45 and the swap search 0.99, in the same order as their Eq. 11 values (1.57, 1.04 and 0.712, against the minimum 0.711); the corresponding hypothesis held on both data sets, where data B gives 0.1064, 0.0591 and 0.0569. The product rule’s Eq. 11 of 1.57 against the attainable 0.711 on data A is the imbalance described above.

On the Gaussian data A, OPQ $_\mathrm{NP}$ started from the swap search ends at almost the same mean test distortion, 0.02% higher (0.9929 against 0.9927, with overlapping seed ranges 0.9679–1.0191 and 0.9680–1.0266), while lowering its training distortion (0.8961 against 0.9048), so the hypothesis that its test distortion would be within 2% of the parametric solution’s held; in our reading, for Gaussian data the parametric rotation is already as good as what alternation finds, as Ge et al. reported for their Gaussian data. On data B it lowers the mean test distortion by 3.9%, from 0.0569 to 0.0547, with seed ranges that do not overlap, so the hypothesis that the non-parametric method would reduce the test distortion by a larger fraction on the non-Gaussian camera patches than on data A, where it reduced nothing, held, consistent with Ge et al.’s finding on SIFT1M and MNIST. On data B, Eq. 11 and $D(R)$ are only Gaussian reference values, not bounds: several training distortions fall below Eq. 11 for their rotation (0.0547 against 0.0702 in natural order). Test distortions on data B are 1.8 to 2.0 times the training distortions in every row; in our reading, the gap reflects both centroids fitted to the training patches and differences in content between the bands of the photograph, which the code does not separate.

E3: distances estimated from codes are biased low. Relative to a mean true squared distance of 131, ADC underestimates by 4.62 on average and SDC by 9.23, about twice as much, and SDC’s error has the larger spread (standard deviation 15.03 against 10.60), as the plan expected; that hypothesis held. The hypotheses that the ADC bias would be between 0.9 and 1.1 times the database distortion and the SDC bias between 1.8 and 2.2 times it both failed: the printed ratios are $-1.183$ and $-2.366$ . On the training set the identity held to rounding, so the derivation is right where its assumption holds. In our reading, the larger held-out bias reflects centroids fitted to the training vectors, whose distortion (3.4959) is lower than that of the test vectors (3.9005 for the database): a centroid is not the mean of the new vectors in its cell. E3 uses seed 0 only. Adding the mean training error of each cell removes most of the bias, leaving $-0.289$ and $-0.571$ times the database distortion; the hypothesis that the correction would leave less than 25% of the uncorrected bias held, narrowly, at 24.4% for ADC and 24.1% for SDC. The correction also reduced the error’s standard deviation slightly in both cases.

Properties and Behavior

The objective is distortion. PQ, OPQ and the variants below minimize reconstruction error, which Ge et al. found tightly correlated with search accuracy; the anisotropic variant changes the objective for inner-product search.
Training is k-means and inherits its behaviour: local optima, dependence on the initialization and empty clusters. In E1, where only the k-means initialization changes with the seed, the distortions varied by about 1% or less across seeds.
Distance estimates are biased low, by about one quantization error for ADC and two for SDC, and more on new data (example, E3). A constant offset does not change a ranking, but the part that varies between database vectors does.
The split matters more than the codebooks. At a fixed code length the distortion changed fifteen-fold with the rotation in the example; Ge et al. found that PQ-based methods strongly depend on the space decomposition on real descriptors too.
Fixed rate. Every vector gets the same number of bits, whatever its distance from the data’s typical values; Ge et al. note that variable-length codes could lower the average rate.

Limitations

Independence between subspaces. PQ cannot represent dependencies between subvectors; Babenko and Lempitsky (2014) note that its performance decreases when those dependencies are strong. A rotation can decorrelate the subspaces, but decorrelation is not independence for non-Gaussian data.
One global quantizer. Rotations and codebooks are fitted to the whole distribution; on a synthetic mixture of 1,000 anisotropic Gaussians, a global OPQ rotation did not improve over IVFADC (Kalantidis & Avrithis, 2014, Fig. 2).
Reconstruction error is not the search objective. For maximum inner-product search, errors on the vectors most similar to the query matter most (Guo et al., 2020).
Retraining. New data from another distribution, or embeddings from a new model, need new codebooks and new codes.
Accuracy ceiling. Distances from codes rank near neighbours imperfectly; systems re-rank a shortlist with exact vectors when they are available, as the approximate nearest neighbor search article describes.

Variants

Optimized product quantization and Cartesian k-means (Ge et al., 2013; Norouzi & Fleet, 2013) learn the rotation with the codebooks, described above.
Locally optimized product quantization (LOPQ) (Kalantidis & Avrithis, 2014) builds on the inverted file: a coarse k-means quantizer assigns each vector to a cell, and the residuals of each cell, which are much more unimodal than the data, get their own OPQ $_\mathrm{P}$ rotation and their own product quantizer. The space overhead is $O(K(D^2 + D k^*))$ for $K$ cells and the query overhead $O(w D^2)$ for $w$ probed cells. On MNIST with 64 cells and 64-bit codes the average distortion per point fell from 70.1 for IVFADC and 12.6 for a global OPQ $_\mathrm{P}$ on the same inverted file to 8.13 (their Fig. 3); a variant for the inverted multi-index optimizes per subspace cell instead of per cell.
Additive quantization (AQ) (Babenko & Lempitsky, 2014) represents a vector as the sum of $m$ codewords, each of full dimension $D$ and from its own codebook, rather than a concatenation. It drops the independence assumption and, on SIFT-1M, gave lower coding error than PQ and OPQ at 4, 8 and 16 bytes (at 8 and 16 bytes in a product variant, APQ), with the improvement particularly favourable for short codes of four or eight bytes, at the cost of encoding by beam search and an extra term $\lVert x \rVert^2$ in distance computation, which the authors can store in one extra byte per vector to avoid its time overhead. FAISS implements additive quantizers such as residual and local-search quantization alongside PQ.
Anisotropic vector quantization (Guo et al., 2020), released as open source in the scann directory of Google Research’s repository, changes the loss for maximum inner-product search. For queries uniformly distributed on the unit sphere, they show that the expected weighted error decomposes into the squared residual parallel to the data point and the squared residual orthogonal to it, and that for a weight on the inner product that is zero for negative values and non-decreasing, the parallel part weighs at least as much, equally only for a constant weight (their Theorems 3.2 and 3.3); the codebooks are trained by a Lloyd-like algorithm with a closed-form update under this loss, also for product quantization. On Glove1.2M they reported higher recall than with the reconstruction loss, and in a speed–recall comparison with 11 baselines in the ANN-Benchmarks set-up their implementation was fastest in the high-recall regime (their Fig. 4b).
Inverted files and multi-indexes. IVFADC encodes residuals from a coarse quantizer; the inverted multi-index applies the product idea to the coarse quantizer itself, with two subspace quantizers whose cells form a grid (as LOPQ’s authors describe it). The approximate nearest neighbor search article covers IVFADC.
Binary codes. Hashing methods such as ITQ produce bit strings compared by the Hamming distance, the codes of binary descriptors and locality-sensitive hashing; Ge et al. show ITQ as a constrained vector quantizer and found it less accurate than OPQ at equal length. FAISS’s IndexPQ also offers Hamming-based and “polysemous” filtering on PQ codes, besides asymmetric (default) and symmetric search.

Use in Computer Vision

Large-scale image search. Global descriptors, aggregated local descriptors such as VLAD, and network embeddings are compressed to a few bytes for image retrieval and near-duplicate detection over billions of images, for example inside an inverted file.
Cross-modal search. Galleries of image embeddings for image–text retrieval are searched with PQ-based indexes; the embeddings of LAION-5B were served from a product-quantized (PQ128) index, as the LAION-5B paper article describes.
Descriptor databases for local feature matching and localization, where the descriptors of a large map must fit in memory.
Compressed features for learning. Babenko and Lempitsky (2014) trained and tested image classifiers on compressed vectors and found that better coding accuracy gave higher classification accuracy.