Operational Quantization Theory

Table of Contents

• 1. Scalar Quantization
  • 1.1 Encoder, Decoder, and Quantization Cells
  • 1.2 MSE of the Uniform Quantizer
  • 1.3 Lloyd-Max Optimality Conditions
• 2. High-Rate Quantization Asymptotics
  • 2.1 The Panter-Dite Formula
  • 2.2 Gaussian Special Case and the 2.72 dB Penalty
• 3. Vector Quantization
  • 3.1 Codebook Definition and Voronoi Cells
  • 3.2 The Lloyd Algorithm
  • 3.3 Asymptotic Optimality and the Curse of Dimensionality
• 4. Product Quantization and Residual VQ
  • 4.1 Subvector Decomposition and Distortion Additivity
  • 4.2 Multi-Stage Residual VQ
  • 4.3 Application to Transformer KV-Cache Compression
• 5. Random Rotation Preprocessing
  • 5.1 Equalizing Marginals via Orthogonal Transforms
  • 5.2 Coordinate Distribution After Rotation
  • 5.3 Why Rotation Lets Scalar Quantizers Behave Like VQ
• 6. The Johnson-Lindenstrauss Lemma
  • 6.1 Statement
  • 6.2 Proof Sketch via MGF and Concentration
  • 6.3 1-Bit JL and Inner Product Preservation
• 7. TurboQuant
  • 7.1 Distortion Bound and Information-Theoretic Lower Bound
  • 7.2 The Constant Gap and Connection to Panter-Dite
  • 7.3 Two-Stage Pipeline for Inner Product Preservation
• References

1. Scalar Quantization 📐

1.1 Encoder, Decoder, and Quantization Cells

Let \(X\) be a real-valued random variable with density \(f\). A scalar quantizer with \(N = 2^b\) reconstruction levels (equivalently, \(b\) bits) consists of two maps:

Definition 1.1 (Scalar quantizer). A \(b\)-bit scalar quantizer is a pair \((Q_e, Q_d)\) where: - The encoder \(Q_e : \mathbb{R} \to \{1, \ldots, 2^b\}\) partitions \(\mathbb{R}\) into \(N\) disjoint quantization cells (or bins) \(\mathcal{R}_1, \ldots, \mathcal{R}_N\) via \(Q_e(x) = k\) iff \(x \in \mathcal{R}_k\). - The decoder \(Q_d : \{1, \ldots, 2^b\} \to \mathbb{R}\) maps each index \(k\) to a reconstruction level (or codeword) \(y_k \in \mathbb{R}\).

The composite map \(Q = Q_d \circ Q_e : \mathbb{R} \to \{y_1, \ldots, y_N\}\) is the quantizer. For a contiguous partition of \(\mathbb{R}\), cell \(\mathcal{R}_k = [t_{k-1}, t_k)\) for breakpoints \(-\infty = t_0 < t_1 < \cdots < t_{N-1} < t_N = +\infty\). The mean squared error (MSE) distortion is

\[D = \mathbb{E}[(X - Q(X))^2] = \sum_{k=1}^{N} \int_{t_{k-1}}^{t_k} (x - y_k)^2 f(x)\, dx.\]

1.2 MSE of the Uniform Quantizer 🔑

The uniform quantizer sets all cell widths equal to \(\Delta\) (the step size): \(t_k = t_0 + k\Delta\) and \(y_k = t_{k-1} + \Delta/2\). Within each cell, \(x - y_k\) ranges uniformly over \([-\Delta/2, \Delta/2]\).

Proposition 1.2 (Granular MSE of uniform quantizer). For a uniform quantizer with step size \(\Delta\) applied to a density \(f\) supported on \(N\Delta\) contiguous cells, as \(\Delta \to 0\),

\[D \approx \sum_k \int_{t_{k-1}}^{t_k} (x - y_k)^2 f(x)\, dx \approx \frac{\Delta^2}{12}.\]

Proof sketch. On cell \(k\), \(f(x) \approx f(y_k)\) (constant to leading order in \(\Delta\)). The within-cell error \(x - y_k\) is approximately \(\mathrm{Uniform}(-\Delta/2, \Delta/2)\) with variance \(\Delta^2/12\). Summing over all \(k\) and using \(\sum_k f(y_k)\Delta \approx 1\):

\[D \approx \sum_k f(y_k) \Delta \cdot \frac{\Delta^2}{12} = \frac{\Delta^2}{12} \cdot 1 = \frac{\Delta^2}{12}. \quad \square\]

The MSE of a \(b\)-bit uniform quantizer on \([-A, A]\) satisfies \(D = \Delta^2/12 = A^2 / (3 \cdot 4^b)\).

1.3 Lloyd-Max Optimality Conditions 🔑

To minimize MSE over all choices of breakpoints \(\{t_k\}\) and reconstruction levels \(\{y_k\}\), differentiate \(D\) with respect to each parameter and set to zero.

Theorem 1.3 (Lloyd-Max conditions). The MSE-optimal scalar quantizer satisfies two necessary conditions simultaneously:

1. Centroid condition. Each reconstruction level is the conditional mean of \(X\) within its cell: \[y_k = \mathbb{E}[X \mid X \in \mathcal{R}_k] = \frac{\int_{t_{k-1}}^{t_k} x f(x)\, dx}{\int_{t_{k-1}}^{t_k} f(x)\, dx}.\]

2. Nearest-neighbor condition. Each breakpoint is the midpoint of adjacent reconstruction levels: \[t_k = \frac{y_k + y_{k+1}}{2}.\]

Derivation sketch. Differentiate \(D\) with respect to \(y_k\): setting \(\partial D / \partial y_k = 0\) gives \(y_k = \mathbb{E}[X \mid X \in \mathcal{R}_k]\). Differentiate with respect to \(t_k\) (a shared boundary between cells \(k\) and \(k+1\)): the integrand equals \((t_k - y_k)^2 f(t_k) - (t_k - y_{k+1})^2 f(t_k)\) at the boundary, and setting this to zero (assuming \(f(t_k) > 0\)) yields \(t_k = (y_k + y_{k+1})/2\).

The Lloyd-Max algorithm iterates these two conditions to convergence (analogous to \(k\)-means). Convergence to a global optimum is not guaranteed; only local optima are found in general.

2. High-Rate Quantization Asymptotics 📊

2.1 The Panter-Dite Formula

For a \(b\)-bit scalar quantizer, increasing \(b\) by 1 adds one bit of resolution. How does optimal distortion scale with \(b\)? The Panter-Dite formula (1951) gives the precise answer in the high-rate limit.

Theorem 2.1 (Panter-Dite). Let \(X\) have density \(f\) with finite variance. For the optimal \(b\)-bit scalar quantizer in the limit \(b \to \infty\),

\[D^*(b) \sim \frac{1}{12} \cdot 4^{-b} \cdot \left(\int_{-\infty}^{\infty} f(x)^{1/3}\, dx\right)^3.\]

Derivation sketch. In the high-rate limit, the quantizer spacing function \(\Delta(x)\) can be treated as locally constant. The optimal cell density \(\lambda(x)\) (cells per unit interval) must satisfy \(\lambda(x) \propto f(x)^{1/3}\) — the result of minimizing total distortion \(\int (1/12)\Delta(x)^2 f(x)\, dx = \int f(x)/(12\lambda(x)^2)\, dx\) subject to the total rate constraint \(\int \lambda(x)\, dx = N = 2^b\). Applying calculus of variations yields \(\lambda(x) = c \cdot f(x)^{1/3}\) with \(c = 2^b / \int f^{1/3}\). Substituting back gives the formula. \(\square\)

The exponent \(1/3\) in \(f(x)^{1/3}\) has a clean interpretation: the optimal quantizer allocates more resolution (smaller cells) where the density is higher, but with diminishing returns — doubling the density only warrants \(2^{1/3} \approx 1.26\times\) more cells.

2.2 Gaussian Special Case and the 2.72 dB Penalty

For \(X \sim \mathcal{N}(0, \sigma^2)\), the integral \(\int f(x)^{1/3}\, dx\) can be evaluated in closed form.

Proposition 2.2 (Gaussian Panter-Dite). For \(f = \mathcal{N}(0, \sigma^2)\):

\[\int_{-\infty}^{\infty} f(x)^{1/3}\, dx = \left(\frac{1}{2\pi\sigma^2}\right)^{1/6} \int_{-\infty}^{\infty} e^{-x^2/(6\sigma^2)}\, dx = \left(\frac{1}{2\pi\sigma^2}\right)^{1/6} \cdot \sqrt{6\pi\sigma^2}.\]

Cubing and substituting into the Panter-Dite formula:

\[D^*_{\text{Gaussian}}(b) \sim \frac{\sigma^2 \pi \sqrt{3}}{2} \cdot 4^{-b}.\]

The asymptotically optimal scalar quantizer for a Gaussian source achieves \(D \approx (\sigma^2 \pi\sqrt{3}/2) \cdot 4^{-b}\).

Now compare against Shannon’s rate-distortion theory bound. For a Gaussian source with \(R = b\) bits per sample, the rate-distortion function gives

\[D_{\text{Shannon}}(b) = \sigma^2 \cdot 4^{-b}.\]

The ratio is

\[\frac{D^*_{\text{Gaussian}}(b)}{D_{\text{Shannon}}(b)} = \frac{\pi\sqrt{3}}{2} \approx 2.7207,\]

which in decibels is \(10 \log_{10}(\pi\sqrt{3}/2) \approx 4.35\) dB. (Since we are comparing MSE values and each bit adds 6 dB, the gap corresponds to \(10\log_{10}(2.72) \approx 4.35\) dB, often called the “2.72 dB penalty” in the sense that the optimal scalar quantizer requires \(\log_4(\pi\sqrt{3}/2)/2 \approx 0.72\) extra bits per sample relative to Shannon’s limit.)

3. Vector Quantization 🗂️

3.1 Codebook Definition and Voronoi Cells

Definition 3.1 (Vector quantizer). A \(b\)-bit vector quantizer of dimension \(d\) is a pair \((Q_e^d, Q_d^d)\) where: - The codebook \(\mathcal{C} = \{y_1, \ldots, y_N\} \subset \mathbb{R}^d\) with \(N = 2^{bd}\) total bits (\(b\) bits per dimension). - The encoder \(Q_e^d(\mathbf{x}) = \arg\min_{k} \|\mathbf{x} - \mathbf{y}_k\|^2\) maps each input \(\mathbf{x} \in \mathbb{R}^d\) to its nearest codeword index. - The decoder \(Q_d^d(k) = \mathbf{y}_k\) maps indices back to codewords.

The Voronoi cell of codeword \(\mathbf{y}_k\) is \(\mathcal{V}_k = \{\mathbf{x} : \|\mathbf{x} - \mathbf{y}_k\| \leq \|\mathbf{x} - \mathbf{y}_j\|\ \forall j\}\). The per-dimension MSE is

\[D = \frac{1}{d}\,\mathbb{E}[\|\mathbf{X} - Q^d(\mathbf{X})\|^2].\]

3.2 The Lloyd Algorithm

The Lloyd algorithm alternates between the two necessary optimality conditions for VQ, directly analogous to Lloyd-Max for scalar quantizers.

import numpy as np

def lloyd_vq(X: np.ndarray, n_codes: int, max_iter: int = 300, tol: float = 1e-6) -> np.ndarray:
    """
    Lloyd's algorithm for training a VQ codebook.
    X: (n_samples, d) array of training vectors
    Returns codebook: (n_codes, d)
    """
    # Random initialisation from data
    rng = np.random.default_rng(0)
    idx = rng.choice(len(X), size=n_codes, replace=False)
    codebook = X[idx].copy()

    prev_distortion = np.inf
    for _ in range(max_iter):
        # Step 1: Nearest-neighbour assignment (Voronoi partition)
        dists = np.linalg.norm(X[:, None, :] - codebook[None, :, :], axis=-1)  # (n, k)
        assignments = np.argmin(dists, axis=1)

        # Step 2: Centroid update
        new_codebook = np.array([
            X[assignments == k].mean(axis=0) if (assignments == k).any() else codebook[k]
            for k in range(n_codes)
        ])

        distortion = np.mean(np.min(dists, axis=1) ** 2)
        if abs(prev_distortion - distortion) < tol:
            break
        codebook = new_codebook
        prev_distortion = distortion

    return codebook

Each iteration is guaranteed to decrease (or maintain) distortion, so the algorithm converges. It finds only a local minimum.

3.3 Asymptotic Optimality and the Curse of Dimensionality

Theorem 3.2 (VQ achieves Shannon R(D)). For an i.i.d. source \(X^d\) with a continuous density and MSE distortion, the optimal \(b\)-bit-per-dimension VQ distortion satisfies

\[\lim_{d \to \infty} D^*_{\text{VQ}}(b, d) = D_{\text{Shannon}}(b),\]

where \(D_{\text{Shannon}}(b)\) is the distortion at rate \(b\) on Shannon’s rate-distortion curve.

This is a consequence of the rate-distortion theorem: a good \((d, 2^{bnd})\) source code exists for block length \(d\) large enough, and VQ over \(\mathbb{R}^d\) is exactly such a code. The scalar penalty \(\pi\sqrt{3}/2\) vanishes as \(d \to \infty\).

💡 Curse of dimensionality. Training a codebook of size \(N = 2^{bd}\) in \(\mathbb{R}^d\) requires at least \(\Omega(N)\) training vectors and \(O(Nd)\) storage. For \(b = 4\) bits per dimension and \(d = 256\), \(N = 2^{1024}\) — entirely infeasible. This motivates structured alternatives (Section 4) and preprocessing tricks (Section 5).

4. Product Quantization and Residual VQ 🔧

4.1 Subvector Decomposition and Distortion Additivity

Product quantization (PQ) circumvents the exponential codebook by decomposing the ambient space into \(M\) subspaces of dimension \(d/M\).

Definition 4.1 (Product quantizer). Partition \(\mathbf{x} \in \mathbb{R}^d\) as \(\mathbf{x} = (\mathbf{x}^{(1)}, \ldots, \mathbf{x}^{(M)})\) with \(\mathbf{x}^{(m)} \in \mathbb{R}^{d/M}\). Given \(M\) independent codebooks \(\mathcal{C}_m\) of size \(K = 2^{b_m}\), the PQ encoder assigns each subvector independently:

\[Q_{\text{PQ}}(\mathbf{x}) = (Q^{(1)}(\mathbf{x}^{(1)}), \ldots, Q^{(M)}(\mathbf{x}^{(M)})).\]

The total codebook has \(K^M\) effective codewords using only \(MK \cdot (d/M)\) storage.

Proposition 4.2 (Distortion additivity). If the subvector quantizers are trained independently, the total MSE decomposes exactly:

\[D_{\text{PQ}} = \frac{1}{d}\sum_{m=1}^{M} \mathbb{E}[\|\mathbf{x}^{(m)} - Q^{(m)}(\mathbf{x}^{(m)})\|^2].\]

This holds because the per-block errors are orthogonal: they live in disjoint coordinate subspaces, so cross-terms in \(\|\mathbf{x} - Q_{\text{PQ}}(\mathbf{x})\|^2\) vanish.

4.2 Multi-Stage Residual VQ

Residual vector quantization (RVQ) uses a cascade of \(L\) codebooks, each quantizing the residual left by the previous stage.

Definition 4.3 (Residual VQ). Define the \(\ell\)-th stage residual recursively: \(\mathbf{r}_1 = \mathbf{x}\), and for \(\ell \geq 2\),

\[\mathbf{r}_\ell = \mathbf{r}_{\ell-1} - Q_{\ell-1}(\mathbf{r}_{\ell-1}).\]

The total reconstruction is \(\hat{\mathbf{x}} = \sum_{\ell=1}^{L} Q_\ell(\mathbf{r}_\ell)\).

With \(b\) bits per stage and \(L\) stages, RVQ uses \(Lb\) bits total. Distortion decreases geometrically in \(L\): after \(L\) stages the residual variance is roughly \(\sigma^2 \cdot (D^*/\sigma^2)^L\) for appropriately normalized \(D^*\). Each additional stage subtracts one codebook’s MSE from the remaining residual.

4.3 Application to Transformer KV-Cache Compression

In autoregressive transformers, the key-value (KV) cache stores attention keys \(\mathbf{k}_t \in \mathbb{R}^d\) and values \(\mathbf{v}_t \in \mathbb{R}^d\) for every past token \(t\) in every layer. Memory is linear in sequence length, making long-context inference expensive.

Compressing each key and value vector from FP16 to \(b\) bits per dimension reduces memory by a factor of \(16/b\). The distortion metric is the per-token inner product error \(|{\langle \mathbf{q}, \hat{\mathbf{k}}_t \rangle - \langle \mathbf{q}, \mathbf{k}_t \rangle}|^2\) rather than MSE, because attention scores depend on inner products, not squared norms.

def rvq_encode(x: np.ndarray, codebooks: list[np.ndarray]) -> list[int]:
    """
    x: (d,) vector
    codebooks: list of L arrays, each (K, d)
    Returns list of L indices
    """
    indices = []
    residual = x.copy()
    for C in codebooks:
        dists = np.linalg.norm(C - residual, axis=1)
        idx = np.argmin(dists)
        indices.append(idx)
        residual = residual - C[idx]
    return indices

5. Random Rotation Preprocessing 🌀

5.1 Equalizing Marginals via Orthogonal Transforms

The Panter-Dite formula reveals that the scalar quantizer’s performance depends on the shape of the coordinate marginals through \(\int f(x)^{1/3}\, dx\). A source with highly non-uniform, heavy-tailed, or skewed marginals will have large distortion relative to a source with uniform or Gaussian marginals. The key observation is that orthogonal transforms preserve \(\ell^2\) norm and hence MSE, but they can dramatically reshape the marginals.

Definition 5.1 (Random rotation). Let \(\mathbf{U} \in \mathbb{R}^{d \times d}\) be drawn uniformly from the Haar measure on the orthogonal group \(O(d)\) (equivalently, obtain \(\mathbf{U}\) via the QR decomposition of a matrix with i.i.d. \(\mathcal{N}(0,1)\) entries). The randomly rotated vector is \(\tilde{\mathbf{x}} = \mathbf{U}\mathbf{x}\).

Since \(\mathbf{U}\) is orthogonal, \(\|\tilde{\mathbf{x}}\| = \|\mathbf{x}\|\) exactly, so no information is lost and MSE is preserved under inversion \(\hat{\mathbf{x}} = \mathbf{U}^\top \hat{\tilde{\mathbf{x}}}\).

5.2 Coordinate Distribution After Rotation

Proposition 5.2 (Rotated Gaussian). If \(\mathbf{x} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}_d)\), then for any fixed orthogonal \(\mathbf{U}\), \(\tilde{\mathbf{x}} = \mathbf{U}\mathbf{x} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}_d)\). The Gaussian distribution is rotationally invariant; rotation changes nothing.

For non-Gaussian sources, however, the marginals of \(\tilde{\mathbf{x}}\) will differ from those of \(\mathbf{x}\). For a fixed unit vector \(\mathbf{x}/\|\mathbf{x}\|\), the marginals of \(\mathbf{U}\mathbf{x}\) are determined by the distribution of a single coordinate of a random point on the \(\ell^2\) sphere of radius \(\|\mathbf{x}\|\).

Proposition 5.3 (Coordinate distribution on the sphere). Let \(\mathbf{z}\) be uniformly distributed on \(S^{d-1}\) (the unit sphere in \(\mathbb{R}^d\)). Then each coordinate \(z_i\) has the marginal density

\[f_{z_i}(t) = \frac{\Gamma(d/2)}{\Gamma((d-1)/2)\,\sqrt{\pi}} \left(1 - t^2\right)^{(d-3)/2}, \quad t \in [-1, 1].\]

This is the Beta distribution \(\mathrm{Beta}((d-1)/2, (d-1)/2)\) in the transformed variable; it concentrates around zero. As \(d \to \infty\), by the central limit theorem applied to \(z_i = U_{i,:}\mathbf{x}/\|\mathbf{x}\|\) with \(U_{i,:}\) a random unit vector:

\[\sqrt{d}\, z_i \xrightarrow{d} \mathcal{N}(0, 1), \quad \text{so } z_i \approx \mathcal{N}(0, 1/d).\]

For a vector \(\mathbf{x}\) with \(\|\mathbf{x}\|^2 = \sigma^2 d\) (i.e., per-dimension variance \(\sigma^2\)), the \(i\)-th coordinate of \(\mathbf{U}\mathbf{x}\) is approximately \(\mathcal{N}(0, \sigma^2)\) for large \(d\).

5.3 Why Rotation Lets Scalar Quantizers Behave Like VQ

After random rotation, two properties hold approximately for large \(d\):

1. Equal marginal variances: \(\mathrm{Var}[\tilde{x}_i] \approx \sigma^2\) for all \(i\).
2. Near-independence: \(\tilde{x}_i\) and \(\tilde{x}_j\) are approximately uncorrelated (and approximately independent for Gaussian sources).

Under these conditions, applying independent per-coordinate scalar quantizers to \(\tilde{\mathbf{x}}\) is nearly equivalent to VQ on the original \(\mathbf{x}\). The per-coordinate Panter-Dite integral \(\int f^{1/3}\) is approximately the same for every coordinate (since all marginals are approximately \(\mathcal{N}(0, \sigma^2)\)), so a single uniform codebook calibrated to \(\mathcal{N}(0, \sigma^2)\) achieves near-optimal distortion on all coordinates simultaneously.

This is the core argument: random rotation + \(d\) independent scalar quantizers approximates optimal VQ, at the cost of one matrix-vector multiply \(O(d^2)\) per vector.

import numpy as np
import time

def haar_rotation(d: int) -> np.ndarray:
    G = np.random.randn(d, d)
    U, _ = np.linalg.qr(G)
    return U

def profile_rotation(d: int, n: int = 1000) -> float:
    U = haar_rotation(d)
    X = np.random.randn(n, d)
    t0 = time.perf_counter()
    X_rot = X @ U.T   # (n, d) @ (d, d) -> O(n * d^2)
    return time.perf_counter() - t0

# Structured alternative: Hadamard-Rademacher transform
# H_d * diag(s) where s are random ±1 signs, H_d is Walsh-Hadamard matrix
# Applied via fast Walsh-Hadamard transform: O(d log d) per vector

6. The Johnson-Lindenstrauss Lemma 📏

6.1 Statement

Theorem 6.1 (Johnson-Lindenstrauss, 1984). For any \(0 < \varepsilon < 1\) and any set of \(n\) points \(\{\mathbf{x}_1, \ldots, \mathbf{x}_n\} \subset \mathbb{R}^d\), there exists a map \(f : \mathbb{R}^d \to \mathbb{R}^k\) with

\[k = O\!\left(\frac{\log n}{\varepsilon^2}\right)\]

such that for all \(i, j \in [n]\):

\[(1 - \varepsilon)\|\mathbf{x}_i - \mathbf{x}_j\|^2 \leq \|f(\mathbf{x}_i) - f(\mathbf{x}_j)\|^2 \leq (1 + \varepsilon)\|\mathbf{x}_i - \mathbf{x}_j\|^2.\]

The map \(f(\mathbf{x}) = \frac{1}{\sqrt{k}}\mathbf{A}\mathbf{x}\) with \(\mathbf{A} \in \mathbb{R}^{k \times d}\) having i.i.d. \(\mathcal{N}(0,1)\) entries achieves this with high probability.

6.2 Proof Sketch via MGF and Concentration

Fix \(\mathbf{u} = \mathbf{x}_i - \mathbf{x}_j\), \(\|\mathbf{u}\| = 1\) (WLOG). The \(k\)-th projected norm is \(Y = \|\mathbf{A}\mathbf{u}\|^2/k = (1/k)\sum_{i=1}^k (\mathbf{a}_i^\top \mathbf{u})^2\) where \(\mathbf{a}_i^\top \mathbf{u} \sim \mathcal{N}(0,1)\) i.i.d. So \(Y = \chi^2_k / k\) with \(\mathbb{E}[Y] = 1\).

Step 1 (MGF bound). For \(t < 1/2\), the MGF of \(\mathbf{a}_i^\top \mathbf{u}\) squared satisfies \(\mathbb{E}[e^{t Z_i}] = (1 - 2t)^{-1/2}\) for \(Z_i = (\mathbf{a}_i^\top \mathbf{u})^2 \sim \chi^2_1\). Thus

\[\mathbb{E}[e^{t Y}] = (1-2t)^{-k/2}.\]

Step 2 (Chernoff bound). By Markov’s inequality applied to \(e^{tY}\):

\[P(Y \geq 1 + \varepsilon) \leq \inf_{t > 0} e^{-t(1+\varepsilon)} (1 - 2t)^{-k/2}.\]

Optimizing over \(t\) gives \(P(Y \geq 1+\varepsilon) \leq e^{-k[\varepsilon^2/4 - \varepsilon^3/6 + \cdots]}\). For \(k \gtrsim C\log(n)/\varepsilon^2\), this decays faster than \(1/n^2\), and a union bound over the \(\binom{n}{2}\) pairs finishes the proof.

6.3 1-Bit JL and Inner Product Preservation

Definition 6.2 (1-bit JL sketch). The 1-bit JL sketch of \(\mathbf{x} \in \mathbb{R}^d\) is

\[s(\mathbf{x}) = \mathrm{sign}(\mathbf{A}\mathbf{x}) \in \{-1, +1\}^k,\]

with \(\mathbf{A} \in \mathbb{R}^{k \times d}\) as above.

Proposition 6.3 (Angular distance preservation). For unit vectors \(\mathbf{x}, \mathbf{y} \in S^{d-1}\),

\[\mathbb{E}\!\left[\frac{1}{k}\langle s(\mathbf{x}), s(\mathbf{y})\rangle\right] = 1 - \frac{2\,\theta(\mathbf{x}, \mathbf{y})}{\pi},\]

where \(\theta(\mathbf{x}, \mathbf{y}) = \arccos(\langle \mathbf{x}, \mathbf{y}\rangle)\) is the angle between the vectors.

This is because \(\text{sign}(\mathbf{a}^\top \mathbf{x}) = \text{sign}(\mathbf{a}^\top \mathbf{y})\) iff \(\mathbf{a}\) does not cross the hyperplane separating \(\mathbf{x}\) and \(\mathbf{y}\), which happens with probability \(1 - \theta/\pi\).

Thus 1-bit JL sketches preserve angular (and hence cosine) similarity, enabling unbiased inner product estimation from sign bits alone.

7. TurboQuant 🚀

TurboQuant (2025) is an online, data-oblivious vector quantization algorithm that achieves near-optimal MSE distortion via random rotation followed by per-coordinate scalar quantization, and further provides an unbiased inner product quantizer via a two-stage residual scheme.

7.1 Distortion Bound and Information-Theoretic Lower Bound

TurboQuant operates on unit-sphere vectors \(\mathbf{x} \in S^{d-1}\) (equivalently, \(\|\mathbf{x}\|^2 = 1\), so per-dimension variance \(\sigma^2 = 1/d\)).

Theorem 7.1 (TurboQuant MSE upper bound). For any \(\mathbf{x} \in S^{d-1}\) and \(b \geq 1\) bits per coordinate, TurboQuant achieves

\[D_{\text{MSE}} \leq \frac{\sqrt{3}\pi}{2} \cdot 4^{-b}.\]

Proof sketch. After random rotation, each coordinate has marginal distribution \(\mathrm{Beta}((d-1)/2, (d-1)/2)\) scaled to \([-1, 1]\), which converges to \(\mathcal{N}(0, 1/d)\) in large \(d\). Applying the Lloyd-Max optimal quantizer to this marginal (whose Panter-Dite integral evaluates to \((\pi\sqrt{3}/2d) \cdot 4^{-b}\) per coordinate), summing over \(d\) coordinates and normalizing by \(d\):

\[D_{\text{MSE}} = \frac{1}{d}\sum_{i=1}^d D_i \leq \frac{1}{d} \cdot d \cdot \frac{\pi\sqrt{3}}{2d} \cdot 4^{-b} \cdot d = \frac{\pi\sqrt{3}}{2} \cdot 4^{-b}. \quad \square\]

Theorem 7.2 (Information-theoretic lower bound). For any (possibly randomized) \(b\)-bit quantizer \(Q\) acting on \(\mathbf{x} \in S^{d-1}\),

\[D_{\text{MSE}}(Q) \geq 4^{-b}.\]

Proof sketch. The rate-distortion function for \(d\) i.i.d. \(\mathcal{N}(0, 1/d)\) coordinates at rate \(b\) bits per coordinate gives minimum distortion \(4^{-b}/d\) per coordinate (using \(R(D) = (1/2)\log(\sigma^2/D)\) with \(\sigma^2 = 1/d\)). Multiplying by \(d\) to get total normalized MSE gives \(\sigma^2 \cdot 4^{-b} = (1/d) \cdot d \cdot 4^{-b} = 4^{-b}\). No quantizer can do better. \(\square\)

The TurboQuant distortion bound \(D \leq (\sqrt{3}\pi/2) \cdot 4^{-b}\) with the lower bound \(D \geq 4^{-b}\) implies the algorithm is within a constant factor \(\pi\sqrt{3}/2 \approx 2.72\) of optimal.

7.2 The Constant Gap and Connection to Panter-Dite

The multiplicative constant \(\pi\sqrt{3}/2\) in the TurboQuant upper bound is exactly the Gaussian Panter-Dite constant from Section 2.2. This is not a coincidence: after random rotation, the coordinate marginals converge to \(\mathcal{N}(0, 1/d)\), and the Panter-Dite formula for a Gaussian source yields \(D \sim (\sigma^2 \pi\sqrt{3}/2) \cdot 4^{-b} = (1/d)(\pi\sqrt{3}/2) \cdot 4^{-b}\) per coordinate. Summing over \(d\) coordinates and normalizing recovers \((\pi\sqrt{3}/2) \cdot 4^{-b}\) for the per-dimension distortion.

In decibels, \(10\log_{10}(\pi\sqrt{3}/2) \approx 4.35\) dB above the Shannon floor. Equivalently, TurboQuant requires roughly \(0.72\) additional bits per dimension relative to the rate-distortion optimum. This gap arises from two compounding sources:

1. The scalar-vs-vector gap (the Panter-Dite penalty inherent to processing one coordinate at a time).
2. The approximation error in matching the post-rotation marginals to i.i.d. Gaussians (exact only as \(d \to \infty\)).

At \(b = 1\) the empirical gap is smaller (\(\approx 1.45\times\)) because the lower bound argument is not tight at very low rates.

7.3 Two-Stage Pipeline for Inner Product Preservation

MSE-optimal quantizers are biased: \(\mathbb{E}[Q(\mathbf{x})] \neq \mathbf{x}\) in general, so \(\mathbb{E}[\langle \mathbf{q}, Q(\mathbf{x})\rangle] \neq \langle \mathbf{q}, \mathbf{x}\rangle\). For KV-cache compression, inner product estimates feed directly into softmax attention and biased estimates degrade generation quality.

TurboQuant resolves this via a two-stage pipeline:

flowchart LR
    x["x ∈ S^(d-1)
input vector"]
    U["Random rotation
x_rot = Ux"]
    Q1["Stage 1: MSE quantizer
(b-1) bits per coord
x_hat = Q_mse(x_rot)"]
    R["Residual
r = x_rot - x_hat"]
    Q2["Stage 2: 1-bit QJL
s = sign(Ar)"]
    OUT["Decoder
combines x_hat + A^T s"]

    x --> U --> Q1 --> R --> Q2 --> OUT
    Q1 --> OUT

Stage 1. Apply the \((b-1)\)-bit MSE-optimal scalar quantizer to each coordinate of \(\tilde{\mathbf{x}} = \mathbf{U}\mathbf{x}\), giving \(\hat{\mathbf{x}} = Q_{\text{MSE}}^{(b-1)}(\tilde{\mathbf{x}})\).

Stage 2. Compute the residual \(\mathbf{r} = \tilde{\mathbf{x}} - \hat{\mathbf{x}}\) and apply a 1-bit JL sketch: \(\mathbf{s} = \mathrm{sign}(\mathbf{A}\mathbf{r})\) for \(\mathbf{A} \in \mathbb{R}^{d \times d}\) random Gaussian. The 1-bit sketch uses exactly \(d\) additional bits (1 bit per coordinate), so the total budget is \((b-1) + 1 = b\) bits per dimension.

Theorem 7.3 (Unbiasedness). The combined estimator \(\tilde{Q}(\mathbf{x}) = \hat{\mathbf{x}} + c\mathbf{A}^\top \mathbf{s}\) (for appropriate constant \(c\)) satisfies

\[\mathbb{E}[\langle \mathbf{q}, \tilde{Q}(\mathbf{x})\rangle] = \langle \mathbf{q}, \mathbf{x}\rangle\]

for all query vectors \(\mathbf{q}\).

This follows because \(\mathbb{E}[c\mathbf{A}^\top \mathrm{sign}(\mathbf{A}\mathbf{r})] = \mathbf{r}\) (up to a scaling constant) by the 1-bit JL unbiasedness property (Proposition 6.3), so the second stage corrects the bias of the MSE quantizer in expectation.

import numpy as np

def turboquant_ip(x: np.ndarray, q: np.ndarray, b: int,
                  U: np.ndarray, A: np.ndarray) -> float:
    d = x.shape[0]
    x_rot = U @ x

    # Stage 1: (b-1)-bit uniform scalar quantizer on each coordinate
    sigma = 1.0 / np.sqrt(d)
    clip = 3 * sigma
    levels = 2 ** (b - 1)
    delta = 2 * clip / levels
    x_hat = np.clip(np.round((x_rot + clip) / delta) * delta - clip, -clip, clip)

    # Residual
    r = x_rot - x_hat

    # Stage 2: 1-bit JL on residual
    s = np.sign(A @ r)  # (d,) binary

    # Scale factor for unbiasedness: E[A^T sign(Ar)] = c * r
    # c = sqrt(2/pi) * ||r|| / d approximately; use empirical scaling
    c = np.sqrt(2 / np.pi) * np.linalg.norm(r) / d
    x_reconstructed = U.T @ (x_hat + c * A.T @ s)

    return np.dot(q, x_reconstructed)

# Verify unbiasedness
d, n_trials = 128, 10_000
U, _ = np.linalg.qr(np.random.randn(d, d))
A = np.random.randn(d, d) / np.sqrt(d)
errors = []
for _ in range(n_trials):
    x = np.random.randn(d); x /= np.linalg.norm(x)
    q = np.random.randn(d)
    true_ip = np.dot(q, x)
    est_ip = turboquant_ip(x, q, b=4, U=U, A=A)
    errors.append(est_ip - true_ip)
print(f"Mean error: {np.mean(errors):.4f}, Std: {np.std(errors):.4f}")
# Expected: mean near 0, std small

References