arXiv CSopen access

Quantizing With Randomized Hadamard Transforms: Efficient Heuristic Now Proven

Quantizing With Randomized Hadamard Transforms: Efficient Heuristic Now Proven Ran Ben-Basat UCL and Broadcom

arXiv:2605.06014v1 [cs.LG] 7 May 2026

Amit Portnoy Microsoft

William Kuszmaul Carnegie Mellon University

Michael Mitzenmacher Harvard University

Shay Vargaftik VMware Research by Broadcom

Abstract Uniform random rotations (URRs) are a common preprocessing step in modern quantization approaches used for gradient compression, inference acceleration, KVcache compression, model weight quantization, and approximate nearest-neighbor search in vector databases. In practice, URRs are often replaced by randomized Hadamard transforms (RHTs), which preserve orthogonality while admitting fast implementations. The remaining issue is the performance for worst-case inputs. With a URR, each coordinate is individually distributed as a shifted beta distribution, which converges to a Gaussian distribution in high dimensions. Generally, one RHT is not suitable in the worst case, as individual coordinates can be far from these distributions. We show that after composing two RHTs on any d-sized input vector, the marginal distribution of every fixed coordinate of the normalized rotated vector is within O(d−1/2 ) of a standard Gaussian both in Kolmogorov distance and in 1-Wasserstein distance. We then plug these bounds into the analyses of modern compression schemes, namely DRIVE and QUIC-FL, and show that two RHTs achieve performance that asymptotically matches URRs. However, while two RHTs suffice for scalar quantization, they may be insufficient for Vector Quantization (VQ), which often requires weak correlation across fixedsize blocks of coordinates (as opposed to only marginal distribution convergence for single coordinates). We prove that a composition of three RHTs leads to decaying coordinate covariance. This ensures that any fixed, bounded, multi-dimensional VQ codebook optimized for URRs has the same expected error when using three RHTs, up to an additive term that asymptotically vanishes with the dimension. Finally, because practical inputs are rarely adversarial, we propose a linear-time O(d) check on the input’s moments to dynamically adapt the number of RHTs used at runtime to improve performance.

1

Introduction

A foundational technique in quantization is applying a uniform random rotation (URR) to the input vectors beforehand, which ensures the vector’s ℓ2 norm is evenly distributed across all coordinates. Moreover, because individual coordinates of a rotated vector rapidly converge to a Gaussian distribution as the vector’s dimension grows, rotations enable stable, nearly dimension-free scalar quantization designs. This robust statistical smoothing has recently found numerous applications, including Distributed Mean Estimation (DME) and gradient compression in federated learning (e.g., QUIC-FL [9]), distributed training (e.g., [22, 21]), privacy (e.g., ScionFL [10]), parameter quantization (e.g., Quartet-2 via MS-EDEN [29]) and activation quantization (e.g., HACK [41]), database vector search (e.g., RaBitQ [20]), KV cache quantization (e.g., TurboQuant [40] and EDEN [6]), and advanced post-training quantization frameworks (e.g., HIGGS [27]).

While URRs are conceptually and mathematically attractive, they take O(d2 ) operations to compute for d-dimensional vectors and thus are often too slow. This is why the literature repeatedly suggests replacing them with randomized Hadamard transforms (RHTs), which are structured random rotations that can be computed more efficiently, i.e., in O(d log d) operations. RHTs have been suggested for many quantization algorithms, e.g., [4, 5, 7, 10, 12, 18, 19, 24, 28, 32, 33, 35, 39, 40]. More broadly, RHTs are also standard structured surrogates for dense random matrices in fast Johnson-Lindenstrauss embeddings and in later kernel-approximation and distance-estimation results. The algorithmic advantage is immediate: the transform is orthogonal, fast, memory efficient, and GPU-friendly. At a theoretical level, however, the delicate issue is that RHTs are more difficult to work with, as they do not fully provide the same properties√as a URR. √ A single RHT may not be suitable for worst-case inputs; for example, the unit vector (1/ 2, 1/ 2, 0, 0, . . . , 0) after an RHT will have all coordinates p take on values in the set {0, ± 2/d}, rather than a (near-)Gaussian coordinate law. Intuitively, even with a smoother initial vector, there can be a gap in terms of readily obtaining results for applications that utilize the coordinate marginal distribution of a random rotation. In this paper, we show that, for scalar distributions, composing two successive RHTs eliminates such bad cases. Specifically, after one RHT, the resulting intermediate vector has a small expected third moment. Given this intermediate vector, each coordinate of the final vector becomes an ordinary Rademacher sum. We then apply the appropriate Berry-Esseen inequalities, yielding the necessary bounds for comparing these coordinate distributions to a Gaussian. We analyze this comparison using two metrics: the Kolmogorov distance, suited to thresholded and bounded-variation functionals, and the 1-Wasserstein distance, suited to Lipschitz functionals. With both estimates established, we then apply them to obtain new formal guarantees for recent scalar quantization (that quantize each entry separately) schemes, namely DRIVE [36] and QUIC-FL [9]. We expect similar guarantees to be derivable for other RHT-based scalar quantizers. Additionally, we turn to the problem of Vector Quantization (VQ), where blocks of entries are jointly quantized. In VQ, having a marginal distribution that is approximately Gaussian is insufficient, as we need weak correlation among different coordinates within each block. We show that this can be achieved by composing three successive RHTs. Finally, while composing RHTs restores strict worst-case guarantees, natural data distributions are rarely perfectly adversarial. A corollary of our general statements is that if an initial vector already exhibits suitably small values for specific norms, additional RHT steps are redundant. We establish that these requirements, specifically, the cubed ℓ3 norm dictating scalar marginals, and the squared ℓ∞ norm dictating multi-dimensional cross-correlation, can be explicitly computed in O(d) linear time. Thus, practical systems can dynamically check an input vector to determine whether 1, 2, or 3 RHTs are strictly required, preserving empirical speed without sacrificing theoretical bounds. Our presentation is intentionally modular. Specifically, we first state general results for composing two RHTs, without application-specific definitions. We then introduce the estimators and quantizers needed for DRIVE [36] and QUIC-FL [9]. We take a similar approach when analyzing compositions of three RHTs. This keeps the core probabilistic statements separate from the downstream subtleties of each method, which we hope provides an easier path for future applications. Our core contributions can be summarized as follows: • Universal Scalar Marginals (2-RHT): We prove that for every arbitrary input vector and every fixed coordinate, two RHTs provide a uniform O(d−1/2 ) normal additive approximation in both the Kolmogorov (dK ) and 1-Wasserstein (W1 ) metrics. • Application of 2-RHT Bounds to Scalar Quantizers: – Formal Guarantees for DRIVE [36] (Table 1a): We plug these scalar W1 bounds into the analysis of DRIVE, recovering the URR-based guarantees for both the biased and the unbiased estimator up to a decaying O(d−1/2 ) term. – Formal Guarantees for QUIC-FL [9] (Table 1b): We apply the dK bound to restore the exact URR-based transmission bandwidth threshold up to an additive O(d−1/2 ) decaying term, eliminating the 3.2× multiplicative penalty of a single RHT. • Decorrelation and VQ (3-RHT): We prove that two RHTs are insufficient for decorrelation of small coordinate blocks, whereas a composition of three RHTs provides the weak dependence required for VQ. This allows any fixed, bounded URR or Gaussian codebook to yield the same expected-error limit, up to a decaying additive term. 2

(a) DRIVE’s vNMSE (§ 4) URR 1-RHT 2-RHT (new)

(b) QUIC-FL’s guarantees (§5)

Biased

Unbiased

(E1 , E2 , E3 , E4 )

Sent frac.

1 − 2/π 0.5 1 − 2/π + εd

π/2 − 1 + εd none π/2 − 1 + εd

(1.520, 0.223, 0.044, 0.0098) (4.831, 0.692, 0.131, 0.0272) (1.520, 0.223, 0.044, 0.0098) + εd

p 3.2·p p + εd

Table 1: A summary of DRIVE’s [36] and QUIC-FL’s [9] guarantees with the known URR (O(d2 ) time) and 1-RHT bounds (O(d log d) time), and the new 2-RHT (O(d log d) time) bounds. Here, εd = O(d−1/2 ). In (b), Eb , for b ∈ {1, 2, 3, 4}, denotes QUIC-FL’s vNMSE at bit budget b. • Adaptive O(d) Linear-Time Check: We translate our theoretical constraints into a simple linear-time method that decides on the required number of RHTs to achieve theoretical guarantees based on the ℓ3 and ℓ∞ norms.

2

Notations and definitions

Randomized Hadamard Transform (RHT). Fix d = 2m . Let H be a Hadamard matrix [23]. Thus, e = √1 H. Thus, H e is orthogonal. H ∈ {±1}d×d and HH ⊤ = dId . Additionally, let H d e k HD e k−1 · · · HD e 1 , where D1 , D2 , . . . , Dk are diagonal matrices with i.i.d. For k ≥ 1, let Rk = HD ±1 Rademacher entries. We call Rk a composition of k RHTs. For convenience, we denote R0 = Id . Kolmogorov distance. Given two distributions whose CDFs are P , Q, their Kolmogorov distance is defined as dK (P, Q) = supx∈R |P (x) − Q(x)|. 1-Wasserstein distance (also known as the Earth Mover’s Distance). Utilizing the KantorovichRubinstein duality [3, 13, 37], we use an alternative definition of the 1-Wasserstein distance given by W1 (X, Y ) = supf ∈Lip(1) E[f (X)] − E[f (Y )] . Unit sphere. We denote by S d−1 the (d − 1)-dimensional unit sphere, i.e., the set of all vectors a ∈ Rd such that ∥a∥2 = 1. Coordinate. For a vector a ∈ Rd , we use ai and (a)i to denote its i’th Coordinate. Binary sign convention. We use sign(z) = 1 if z ≥ 0 and sign(z) = −1 otherwise, applied coordinatewise to vectors. x Normalized input For an input x ∈ Rd \ {0}, let x e = ∥x∥ ∈ S d−1 be its unit direction vector. 2

3

Composing Two RHTs

We start by considering the setting of composing two RHTs, under the Kolomogrov and 1-Wasserstein distance metrics. Specifically, we are interested in the properties of R2 x for any x ∈ Rd \ {0}. 3.1

Kolmogorov distance with two RHTs Pd−1 Let a ∈ S d−1 and define ρ3 (a) = i=0 |ai |3 . Intuitively, the ℓ3 norm cubed ρ3 (a) measures how “flat” the unit vector a is. If the coordinates of a all have size about d−1/2 , then ρ3 (a) is of order d−1/2 . If one coordinate is large, then ρ3 (a) is large. Thus, the analysis for two RHTs works by showing that the first RHT makes ρ3 sufficiently small on average, and then that the second RHT ensures that an individual coordinate is close to a Gaussian-distributed random variable when ρ3 (a) is small. We start with the following lemma. Lemma 1 (Average ℓ3 smoothing after one RHT). For every x e ∈ S d−1 and a random sign diagonal C3 3/4 e √ ≈ 2.2795. D, it holds that, Eρ3 (HDe x) ≤ d where C3 = 3 e x. Thus, a1 = √1 Pd−1 xi εi , where the εi are i.i.d. Rademacher ±1 signs (the Proof. Let a = HDe i=0 d signs of H are absorbed into the Rademacher signs). Since all coordinates have the same marginal 3

Pd−1 2 distribution, Eρ3 (a) = d E|a1 |3 . Let M = i=0 xi εi . Then EM 2 = ∥x∥2 = 1, and a direct fourthP 2 P P Pd−1 d−1 d−1 2 moment calculation gives EM 4 = i=0 x4i + 6 i<j x2i x2j = 3 − 2 i=0 x4i ≤ 3. i=0 xi By Hölder’s inequality, E|M |3 ≤ (EM 4 )3/4 ≤ 33/4 = C3 . Therefore, E|M |3 C3 1 Eρ3 (a) = d E|a1 |3 = d · E| √ M |3 = d · d−3/2 E|M |3 = √ ≤√ . d d d We also utilize a well-known Berry-Esseen bound for non-i.i.d variables derived in [34, Theorem 7] which we restate here for clarity. Theorem 7 of [34] (Berry–Esseen for non-i.i.d variables). Let X0 , . . . , Xd−1 be independent realPd−1 valued random variables with EXi = 0 and i=0 EXi2 = 1. If G ∼ N (0, 1), then ! d−1 d−1 X X sup P Xi ≤ t − P(G ≤ t) ≤ 0.5606 E|Xi |3 . t∈R

i=0

i=0

With these results at hand, we are now ready to prove the following theorem.

√ Theorem 1 (Kolmogorov distance with two RHTs). Let x ∈ Rd \ {0}. Define U (x) = d (R2 x e)1 .
√ C3 ≤ 1.28 √ , where G ∼ N (0, 1). Since the coordinates of R2 x are Then dK U (x), G ≤ 0.5606 d d √ identically distributed, the same estimate holds for every coordinate of d R2 x e. e 1x Proof. Let a = HD e = R1 x e ∈ S d−1 . Therefore, √ U (x) =

√ d (R2 x e)1 =

e 2 R1 x d (HD e)1 = (HD2 R1 x e)1 = (HD2 a)1 =

d−1 X

a i εi ,

i=0

where εi are D2 ’s i.i.d. Rademacher signs. By Theorem 7 of [34], dK U (x) | a, G ≤ 0.5606 ρ3 (a) . Taking expectation and applying Lemma 1 yields

0.5606 C3 √ dK U (x), G ≤ Ea [dK U (x) | a, G ] ≤ 0.5606 E[ρ3 (a)] ≤ . d 3.2

1-Wasserstein distance with two RHTs

We utilize and restate here for clarity the following lemma. Lemma 2.4 of [13](W1 Berry-Esseen Method1 ). Let X0 , . . . , Xd−1 be independent zero-mean Pd−1 random variables with i=0 E[Xi2 ] = 1. Let G ∼ N (0, 1) be a standard Gaussian. Then: ! d−1 d−1 X X   W1 Xi , G ≤ CW E |Xi |3 i=0

i=0

where CW ≤ 3 is an absolute constant. Now we move to bound the 1-Wasserstein distance. Theorem 2 (1-Wasserstein distance with two RHTs). Let x ∈ Rd \ {0}. Define again U (x) = √
CW d (R2 x e)1 . Then W1 U (x), G ≤ √dC3 , where G ∼ N (0, 1). Since the coordinates of R2 x are √ identically distributed, the same estimate holds for every coordinate of U (x) = d R2 x e. Pd−1 e 1x Proof. Again, let a = HD e = R1 x e ∈ S d−1 and U (x) = i=0 ai εi , where εi are D2 ’s i.i.d. Rademacher signs. Notice that conditioned on a, the random variables Xi = ai εi are independent, Pd−1 Pd−1 have zero mean, and satisfy i=0 E[(ai ϵi )2 ] = i=0 a2i = 1. 1 Also appears in [14, Section 3].

4

Applying Lemma 2.4 of [13] yields d−1 d−1 X X   3 W1 (U (x) | a, G) ≤ CW E |ai ϵi | | a = CW |ai |3 = CW ρ3 (a) . i=0

i=0

Taking expectation and applying Lemma 1 yields   W1 (U (x), G) = sup E E[f (U (x)) | a] − E[f (G)] f ∈Lip(1)

  CW C3 . ≤ Ea W1 (U (x) | a, G) ≤ CW E[ρ3 (a)] ≤ √ d

4

Applying the 2-RHT results to DRIVE

The DRIVE algorithm [36] provides a highly efficient method for biased and unbiased 1-bit (per entry) quantization. Then, in DRIVE, the input vector x ∈ Rd is rotated by Rx where R is a random rotation matrix. Then, the rotated vector is deterministically quantized to a binary vector by taking its sign. To reconstruct the vector, the estimate is x̂ = S · R−1 sign(Rx), for some scalar S ∈ R. Importantly, the representation of the quantized vector includes just the sign vector. R is not stored and is generated independently by the dequantizer using shared randomness [8]. This is a standard assumption [15, 26, 38] and can be implemented by using the same PRNG seed. The same applies when using RHT hereafter. In the biased configuration, the objective is to strictly minimize the vector’s Normalized Mean E[∥x−x̂∥22 ] Squared Error (vNMSE = ). DRIVE achieves this by calculating the optimal MSE∥x∥2 2

minimizing scale factor S = d1 ∥Rx∥1 . Because the randomized rotation forces the coordinates to behave similarly to Gaussian variables in high dimensions, this biased scaling yields a vNMSE that asymptotically converges to exactly 1 − π2 ≈ 0.363 for large d. ∥x∥2

∥x∥2

2 . Since the latter In the unbiased configuration, the scale factor is either S = ∥Rx∥21 or S = E∥Rx∥ 1 is more amenable to formal analysis and, as shown in [36, Appendix A.4], yields a tighter vNMSE bound, we use it. Both these specific scalar corrections explicitly guarantee that the reconstructed vector is unbiased, meaning E[x̂] = x. Imposing this constraint slightly increases the individual vNMSE, which asymptotically converges for both to π2 − 1 ≈ 0.571.

When quantizing low-dimensional vectors, applying URRs is computationally tractable and provides rigorously bounded theoretical guarantees. However, in practical implementations, high-dimensional data overwhelmingly favor RHT. Yet, as highlighted in the analysis of the DRIVE algorithm, this substitution is not mathematically flawless; while the RHT successfully mimics a uniform rotation for many natural inputs, it can structurally fail to uniformize certain worst-case vectors (such as sparse inputs). In particular, DRIVE [36] gives a weaker vNMSE ≤ 0.5 bound for the biased case when using RHT instead of URR, and no bound for the unbiased RHT-based variant. We now show that, after two RHTs, DRIVE recovers the URR vNMSE bounds for both the biased and unbiased variants, up to an additive term inversely polynomial in d. 4.1

The biased configuration

By Lemma 1 of the DRIVE, when using S = d1 ∥Rx∥1 , for any rotation matrix R, the vNMSE is: i 1 h x∥21 . (1) vNMSER = 1 − E ∥Re d For a URR, as mentioned, the error converges to exactly 1 − π2 ≈ 0.363 for large d. Our goal is to show that a matching asymptotic upper bound holds for R ≜ R2 . By Jensen’s inequality: i 1  2 1 h E ∥R2 x e∥21 ≥ E ∥R2 x e∥1 . d d   Pd−1 By linearity of expectation, E ∥R2 x e∥1 = i=0 E[|(R2 x e)i |] = d · E[|(R2 x e)1 |] . 5

Substituting this into our squared bound yields 2 i 1
2 1 h d E[|(R2 x e)1 |] = d · E[|(R2 x e)1 |] . E ∥R2 x e∥21 ≥ d d Substituting this into Equation (1) yields
2 vNMSER2 ≤ 1 − d · E[|(R2 x e)1 |] .

(2)

We have now expressed a bound on vNMSER2 in terms of a single coordinate. We now use Lemma 2 below, which uses the fact, similar to the URR, that any single coordinate is close to a Gaussian distribution. As our lower bound can be negative for small d, we define (z)+ = max(z, 0) to obtain: r  2  
2 2 2 CW C3 1 d · E[|(R2 x e)1 |] ≥ = −O √ − √ π π d + d Substituting this into Equation (2) yields the desired upper bound:  r  2 
2 2 CW C3 1 2 − √ vNMSER2 ≤ 1 − d · E[|(R2 x e)1 |] ≤ 1 − =1− +O √ . π π d + d √ x and define U (x) = d · (R2 x e)1 . Then, Lemma 2. Let x ∈ Rd \ {0} be the input vector, x e = ∥x∥ 2 r   2 CW C3 E |U (x)| ≥ − √ . π d Proof. For any u, v ∈ R, by the reverse triangle inequality, i.e., |u| − |v| ≤ |u − v| , the absolute value function f (z) = |z| is exactly 1-Lipschitz. Therefore, W1 bounds the L1 norm difference: E[|X|] − E[|Y |] ≤ W1 (X, Y )

(3)

Let G ∼ N (0, 1). Using Equation (3), rearranging, and applying Theorem 2 yields, r 2 CW C3 E[|U (x)|] ≥ E[|G|] − W1 (U (x), G) ≥ − √ . π d 4.2

The unbiased configuration

We now analyze the bias and variance of the unbiased 2-RHT-based DRIVE estimator for an arbitrary input vector x ∈ Rd \ {0}. As defined at the beginning of Section 4, to guarantee unbiasedness under ∥x∥2 a URR, the algorithm scales the quantized vector by S = E[∥RU 2x∥1 ] where RU is a random rotation √2 , where cd ≜ √1 E[∥RU x matrix. This simplifies to S = c∥x∥ e∥1 ] is the scaled L1 norm expectation d d d

for a uniformly distributed unit vector RU x e ∈ S d−1 . Recall the estimate: x̂ = SR2−1 sign(R2 x) =

∥x∥2 −1 √ R2 sign(R2 x e). cd d

Let us define the expected direction vector generated by the 2-RHT as µ = √1d E[R2−1 sign(R2 x e)]. 2 By factoring out the scalars, our expected estimator evaluates to E[x̂] = ∥x∥ cd µ. The estimation bias is the vector difference between this expectation and the input vector x:

B(x) ≜ E[x̂] − x =

∥x∥2 (µ − cd x e). cd

(4)

With this definition at hand, we continue with the following theorem. Theorem 3. Let x ∈ Rd \ {0} be an arbitrary vector, and define the 2-RHT DRIVE estimator with the √2 R−1 sign(R2 x). As d → ∞, the estimator’s normalized uniform-rotation unbiased scale by x̂ = c∥x∥ 2 d d squared bias decays to 0, and its normalized variance converges to the same value as for a URR: Var(x̂) π ∥B(x)∥22 = − 1 + O(d−1/2 ) (ii) ≤ O(d−1/2 ) (i) ∥x∥22 2 ∥x∥22 6

Proof intuition. To prove the bias decays to zero, we show that µ converges to the target vector cd x e. We achieve this by decomposing µ into a parallel component (its projection onto x e) and an orthogonal component (µ⊥ ). Analyzing these components directly against the algorithm’s constant cd involves cumbersome algebra. Instead,pour proof leverages a simpler proxy scale: the expected absolute value of a standard Gaussian, c ≜ 2/π. Because 2-RHT coordinates approximate standard Gaussians (Theorem 2), the scaled expected L1 norm of any unit vector under 2-RHT concentrates around c. We use this to bound both components: 1. Parallel component: The projection of µ onto x e equals the expected L1 norm of R2 x e, which is close to c. 2. Orthogonal component: We show that the perpendicular component µ⊥ is small, and we then bound its effect on the overall bias using the Pythagorean theorem. Finally, because the algorithm’s scale cd converges to c with an O(d−1 ) gap, we substitute c for cd at the end of the proof to establish the exact variance and bias bounds. The formal proof, along with the required asymptotic expansion of the scaling constant cd , is deferred to Appendix A. Implications for Distributed Mean Estimation (DME). When N arbitrary vectors x(c) are averaged PN −1 via their independent 2-RHT-based estimates x̂avg = N1 c=0 x̂(c) (each having a decaying bias), the Normalized Mean Squared Error (NMSE) decomposes into the aggregated variance and squared bias. Because vectors use independent 2-RHT rotations, the resulting overall variance scales down by a 1/N factor. Concurrently, by Jensen’s inequality, the squared magnitude of the aggregated bias vector is bounded by the average of the squared individual biases. Thus, the resulting error is close to π −1 E[∥x̂ −xavg ∥22 ] the one achieved via uniform rotations: NMSE ≜ 1 Pavg ≤ 2N + O(d−1/2 ) . N −1 ∥x ∥2 N

5

c=0

(c) 2

Applying the 2-RHT results to QUIC-FL

State-of-the-art Distributed Mean Estimation algorithms, such as QUIC-FL [9], rely heavily on Bounded Support Quantization (BSQ).2 In BSQ, the algorithm transmits coordinates that fall outside a predefined range at baseline precision, which often correspond to 32- or 16-bit values to avoid outlier errors, while stochastically quantizing the rest. Following the notation of [9], coordinates outside the range [−tp , tp ] are sent exactly (i.e., at baseline precision), where the subscript p represents that tp is chosen so that G ∼ N (0, 1) lies outside this range with probability p. This means that after a URR, each coordinate falls outside this range with probability approximately and asymptotically p. As discussed in [9], when QUIC-FL uses a (single) RHT instead of a URR, bounds can still be proven, but they are significantly weaker. In particular, the issue is that the expected fraction of items falling outside [−tp , tp ] can be substantially more than p. An inequality proven in [11] shows that for worst-case inputs, the tail probability of a coordinate obtained using a single RHT can be up to 3.1787p ≈ 3.2p; that is, around 3.2 times larger than that of a standard Gaussian. This multiplicative penalty in the bound on outliers cascades through the QUIC-FL analysis, leading to overestimates of both bandwidth provisioning and expected quantization error. Our results for composing two RHTs are directly applicable and remove this penalty. Specifically, using our results for the Kolmogorov metric (Theorem 1), we show that QUIC-FL can achieve essentially the same bounds when using two composed RHTs instead of a URR, up to terms asymptotically vanishing in the dimension. 5.1

Expected outlier fraction and bandwidth bounds

Recall that we let p = P(|G| > tp ) be the expected fraction of coordinates falling outside the range [−tp , tp ] for a standard Gaussian. With a single RHT, the expected √ fraction of exactly-sent coordinates (Theorem G.1 in [9]) can be bounded as follows: let Z = d(R1 x e)1 be the first coordinate after

one RHT; then P |Z| > tp ≤ 3.2 · P |G| > tp = 3.2·p. The bound on the expected number of coordinates requiring full precision is thus much higher than for URRs. 2 See also [16, 17, 25].

7

When composing two RHTs, our Kolmogorov bound (Theorem 1) guarantees that the expected fraction of outliers is asymptotically the same as for the Gaussian distribution. √ Theorem 4 (Expected Outlier Fraction with Two RHTs). As before, let U (x) = d(R2 x e)1 . Let tp > 0 be a threshold such that P(|G| > tp ) = p for G ∼ N (0, 1). Then, the probability that a coordinate falls outside [−tp , tp ] after two RHTs is bounded by:

2.56 P |U (x)| > tp ≤ p + √ = p + O d−1/2 . d √ Proof. Because the coordinates of d(R2 x e) are identically distributed, the expected fraction of outliers equals the probability that a single coordinate exceeds the threshold tp . We rewrite this tail probability in terms of the cumulative
distribution function and
apply the Kolmogorov bound from Theorem 1. Let FU (t) = P U (x) ≤ t and FG (t) = P G ≤ t denote the cumulative distribution functions of U (x) and G, respectively. Then:


P |U (x)| > tp = P U (x) > tp + P U (x) < −tp ≤ 1 − FU (tp ) + FU (−tp )
≤ 1 − FG (tp ) + FG (−tp ) + 2dK (U (x), G) = P |G| > tp + 2dK (U (x), G). √ by Theorem 1, substituting this bound yields the desired result. Since dK (U (x), G) ≤ 1.28 d

As d grows, the expected fraction of outliers converges to p, saving ≈ 70% of the bandwidth overhead if using the bound above for a single RHT. 5.2

Quantization error and bounded total variation

To bound the actual quantization error (vNMSE) of BSQ, we define the expected squared error function conditioned on the coordinate value: e(z) = E[(ẑ − z)2 | z]. For |z| ≤ tp , e(z) is defined by the stochastic quantizer. For |z| > tp , the values are sent exactly, meaning e(z) = 0. This creates a jump discontinuity at the thresholds ±tp . Because the error drops to zero, the function is non-Lipschitz. Consequently, the 1-Wasserstein metric used in Section 4 cannot be applied to bound this expected error. Instead, we use the Kolmogorov distance (dK ), which can bound the expectation of any function with bounded total variation (T V ). For a function Pn with jump discontinuities, the total variation is appropriately defined by T V (f ) = supz0 <···<zn j=1 |f (zj ) − f (zj−1 )| < ∞. Theorem 5 (Quantization Error with Two RHTs). Let e : R → R be right-continuous, of bounded variation and, with finite limits at ±∞ (representing the expected squared error function for BSQ). d Let T V (e) denote its total variation. For any √ input x ∈ R \ {0}, the expected quantization error of a scaled transformed coordinate U (x) = d(R2 x e)1 after two RHTs satisfies: E[e(U (x))] − E[e(G)] ≤ T V (e)dK (U (x), G) ≤

1.28 · T V (e) √ , d

where G ∼ N (0, 1). √ Proof intuition. According to Theorem 1, a scaled transformed coordinate U (x) = d(R2 x e)1 is close to Gaussian in Kolmogorov distance. Since Kolmogorov distance controls how much probability the two distributions assign below any threshold, and since T V (e) measures the total amount by which the quantization-error function can change, the expected value of e can change by at most their product. The complete proof is deferred to Appendix B. Refined worst-case error bounds. In [9, Theorem G.3], computing worst-case bounds for the quantization error Eb required partitioning the support into coarse intervals (e.g., [0, 1.5]) and penalizing the entire probability mass of that interval with the maximal error occurring at its boundary. This Riemann sum overestimation was necessary due to the 1-RHT multiplicative tail bound, leading to inflated theoretical bounds (e.g., producing 1-bit quantization error E1 ≤ 4.831 as d → ∞). By substituting 1-RHT with 2-RHT, the theoretical vNMSE bounds converge to the continuous values (identified numerically in [9, Section 3.5]). Transitioning from 1-RHT to 2-RHT thus recovers the desired theoretical error limits up to a decaying O(d−1/2 ) term. 8

Bit Budget 1-RHT Bound (Thm G.3) 2-RHT Bound Improvement b=1 E1 ≤ 4.831 E1 ≈ 1.520 3.17× tighter b=2 E2 ≤ 0.692 E2 ≈ 0.223 3.10× tighter b=3 E3 ≤ 0.131 E3 ≈ 0.044 2.97× tighter b=4 E4 ≤ 0.0272 E4 ≈ 0.0098 2.77× tighter Table 2: Comparing the b-bit vNMSE, i.e., Eb , as a function of the bit-budget b for 1-RHT and 2-RHT, adapting results from [9].

6

Beyond scalar quantization: Three RHTs for VQ

While two-RHTs are sufficient to ensure the one-dimensional marginal distributions required for scalar rotation-based quantization algorithms like DRIVE and QUIC-FL, some compression schemes require multi-dimensional VQ in which coordinates are partitioned into sets of size k and treated as k-dimensional vectors. This is particularly prominent in frameworks built for inference, KV cache quantization, and vector databases, e.g., [20, 27, 29]. Vector quantizers group multi-dimensional data into localized regions based on the nearest representative point. For standard codebooks to be efficient, the coordinates within a k-sized block must not only possess Gaussian marginals but must also be weakly correlated. If coordinates are highly correlated, the data clusters in a lower-dimensional subspace, leading to significant quantization error. In this section, we show that two RHTs are insufficient for VQ because they do not decorrelate sparse inputs. However, we establish that adding a third RHT yields a weak coordinate correlation sufficient for using Gaussian codebooks while maintaining their theoretical performance. Here, we consider Hadamard matrices constructed via the standard Sylvester’s construction. 6.1

The correlation bottleneck of Two RHTs

To understand why two RHTs may be insufficient for VQ, consider the conditional covariance between two coordinates of the transformed vector. Lemma 3 (2-RHT Correlation √ Bottleneck). There exist input vectors x e ∈ S d−1 and distinct coordie 1x nates i ̸= j such that, for U = d R2 x e and a = HD e, for every realization of the sign matrix D1 , CovD2 (Ui , Uj | a) = ±1 . Thus, the two coordinates are perfectly conditionally correlated or anti-correlated. Proof intuition. When evaluating the covariance between two output coordinates Ui and Uj (conditioned on the intermediate vector after the first RHT), the properties of the Hadamard matrix collapse the sum into pairwise products of the original input coordinates. If the input is highly sparse (e.g., only two non-zero coordinates), this sum reduces to a single surviving term. Because the marginal variance of each coordinate is 1, a conditional covariance of ±1 dictates that the coordinates are perfectly conditionally correlated or anti-correlated. The derivation is deferred to Appendix C.1. We note that the distinction between unconditional and conditional covariance is essential. Averaging over both RHTs may make two coordinates appear uncorrelated, but a vector quantizer operates on a block generated by one realized transform. As the lemma above shows, for sparse inputs under two RHTs, the conditional law of a coordinate pair can be supported on a one-dimensional subspace, even though the unconditional covariance vanishes after averaging over the first RHT. Thus, unconditional decorrelation is insufficient for using Gaussian VQ codebooks. Accordingly, we next show that a 3-RHT sequence bounds the conditional covariance itself in RMS, ruling out such bad cases for fixed coordinate pairs. 6.2

The 3-RHT conditional decorrelation guarantee

By using three RHTs, we can eliminate the conditional correlation patterns observed for two RHTs. √ Theorem 6 (3-RHT Decorrelation). Let x e ∈ S d−1 , y = R1 x e, U = dR3 x e, and for i ̸= j define the conditional covariance Ci,j (y, D2 ) = CovD3 (Ui , Uj | y, D2 ). Then, for every given y = R1 x e,: 9

  ED2 Ci,j (y, D2 ) = 0 , and the Root Mean Square (RMS) of this D3 -conditional over qcovariance   
 
1/2 1/2 log d D1 , D2 decays: ED1 ,D2 Ci,j (R1 x e, D2 )2 ≤ 2 ED1 ∥R1 x e∥2∞ ≤O . d Proof intuition. While the 2-RHT may fail, e.g., on sparse inputs, the first of a 3-RHT sequence makes the input to the final two RHTs flat on average, bounding its expected squared ℓ∞ norm by O( logd d ). For any fixed output y of the first RHT, we show that the final two RHTs decorrelate coordinates at a scale controlled by how flat y is. Since the first RHT makes y flat on average, the RMS conditional covariance decays as claimed. The complete proof is given in Appendix C.2. 6.3

Codebook universality: Matching the URR limit

In k-dimensional VQ, a codebook C = {c1 , . . . , cM } ⊂ Rk maps each k-dimensional block of the vector to its nearest centroid. The expected quantization error E(C, V ) is the expected squared distance between the block V and its selected centroid from the codebook C. For a URR, any k-dimensional block converges to a standard multivariate Gaussian Z ∼ N (0, Ik ) as d → ∞. Thus, the optimal codebook is derived via k-means clustering over this distribution. We show that a 3-RHT achieves the same quantization error up to a decaying term. √ Theorem 7 (3-RHT Codebook Universality). Let x e ∈ S d−1 and U = dR3 x e. Fix the block size k and a finite codebook C ⊂ Rk with radius B = maxc∈C ∥c∥2 . Then, the difference in expected quantization error satisfies: r  h i h i log d 2 2 E min ∥U0:k−1 − c∥2 − E min ∥Z − c∥2 ≤ Ok,B . c∈C c∈C d Proof intuition. We utilize the multivariate generalization of Stein’s Method to bound the 1Wasserstein distance between the 3-RHT block and a standard multivariate Gaussian. This distance is governed by two components: the deviation of the block’s covariance matrix from the identity Ik , and its expected third absolute moments. Our decorrelation guarantee (Theorem 6) bounds the covariance mismatch, while the 1-RHT smoothing property (Lemma 1) bounds the third moments. The only additional issue is that the VQ error function for vector v, minc∈C ∥v − c∥22 is quadratic and therefore not globally Lipschitz. Thus, we separate the quadratic part from the codebook-dependent part via
minc∈C ∥v − c∥22 = ∥v∥22 + minc∈C ∥c∥22 − 2⟨v, c⟩ . The second term is globally 2B-Lipschitz, and the first term cancels in expectation because both the 3-RHT block and Z ∼ N (0, Ik ) have expected squared norm k. The formal proof is deferred to Appendix C.3.

7

Adaptive linear-time verification

Because practical data distributions are often not adversarial, we can dynamically adapt the required number of RHTs. Our proofs show that the number of required RHT layers is governed by two x quantities of the normalized input x e = ∥x∥ . For the scalar marginal normal approximation, the 2 Pd−1 relevant quantity is the cubed ℓ3 mass ρ3 (e x) = r=0 |e xr |3 . For the RMS conditional correlation estimate used in the vector-quantization analysis, the relevant quantity is the squared norm ℓ∞ ∥e x∥2∞ . Thus, before applying the worst-case number of RHTs, a system can first check whether the input is already sufficiently flat. Scalar check. If ρ3 (e x) ≤ η3 , then the one-coordinate Berry–Esseen step used √in Theorems 1 and 2 e x)i , we have already applies to a single RHT. Namely, for every fixed coordinate Ui = d (HDe dK (Ui , G) ≤ 0.5606 η3 , W1 (Ui , G) ≤ CW η3 , where G ∼ N (0, 1). Therefore, if η3 = O(d−1/2 ), then 1-RHT already gives the same asymptotic guarantees as the worst-case 2-RHT theorem. Vector check. For VQ, scalar marginal normality is insufficient; the relevant quantity is the conditional covariance between coordinates in a fixed block. While Lemma 3 identifies a worst-case correlation bottleneck for 2-RHT, the analysis in Theorem 6 implies that if the input is sufficiently flat, ∥e x∥2∞ ≤ η∞ , then the decorrelation guarantees of the 3-RHT construction are already met by 10

√ only two RHT layers. More precisely, let U = d NHDb NHDa x e be the 2-stage RHT of x e. For distinct coordinates i ̸= j, define the conditional covariance as, Cij (e x, Da ) = CovDb (Ui , Uj | x e, Da ) . The covariance calculation then yields,  
1/2 √ EDa Cij (e x, Da )2 ≤ 2 η∞ .   Thus, if η∞ = O logd d , the preliminary smoothing RHT is unnecessary, and 2-RHT suffices to reach the same expected-error limit guaranteed by the worst-case 3-RHT construction in Theorem 6. Pd−1 2 Implementation. Both checks can be evaluated in one pass over x. Define S2 = r=0 xr , Pd−1 M2 S3 2 3 2 x∥∞ = S2 . Therefore, the S3 = r=0 |xr | , M2 = max0≤r<d xr . For x ̸= 0, ρ3 (e x) = 3/2 , ∥e S2

adaptive GPU-friendly verification step costs O(d) time and O(1) space. In summary, scalar quantization methods such as [9, 36] can use the ρ3 (e x) check to decide whether one RHT is already sufficient, while VQ methods can use the ∥e x∥2∞ check to decide whether the preliminary smoothing RHT in the three-RHT construction can be skipped.

8

Conclusion

The transition from dense URRs to fast, structured transformations like the Randomized Hadamard Transform (RHT) is a fundamental system-level optimization for modern quantization frameworks. Historically, however, this transition has carried a heavy theoretical penalty for worst-case inputs, yielding weak bounds and reducing theoretically optimal algorithms to practical heuristics. In this work, we demonstrated that this penalty is not an inherent limitation of structured matrices, but simply a consequence of using an insufficient number of transforms. By composing two RHTs, we established universal O(d−1/2 ) convergence bounds in both the Kolmogorov and 1-Wasserstein metrics for any arbitrary input vector. We demonstrated that these scalar guarantees seamlessly plug into the analyses of state-of-the-art scalar quantization schemes. By doing so, we successfully recovered the exact URR-based guarantees for DRIVE and the Bounded Support Quantization (BSQ) error integrals for QUIC-FL, without the prohibitive computational cost. Moving beyond scalar quantization, we identified a critical conditional correlation bottleneck on sparse inputs that structurally breaks 2-RHT configurations for VQ. We resolved this by proving that a 3-RHT sequence successfully decorrelates block coordinates, preserving the theoretical optimality of Gaussian VQ codebooks via Stein’s Method. Crucially, to ensure that these robust worst-case guarantees do not burden average-case execution, we established that simple linear-time O(d) checks on the input’s ℓ3 and ℓ∞ norms can dynamically dictate whether one, two, or three RHTs are strictly required at runtime. Our modular probabilistic framework opens several promising directions for future research. For example, while three RHTs are sufficient to guarantee weak dependence for standard fixed-size VQ blocks, exploring a four-RHT (4-RHT) composition presents a natural extension. Specifically, a 4-RHT sequence could provide enhanced independence structures through a fourth RHT step, which may be necessary to formally guarantee decorrelation across non-fixed-size blocks required by next-generation adaptive quantizers.

References [1] Milton Abramowitz and Irene A Stegun. Handbook of mathematical functions with formulas, graphs, and mathematical tables, volume 55. US Government printing office, 1948. [2] Nir Ailon and Bernard Chazelle. Approximate nearest neighbors and the fast johnsonlindenstrauss transform. In Proceedings of the thirty-eighth annual ACM symposium on Theory of computing, pages 557–563, 2006. [3] Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In International conference on machine learning, pages 214–223. Pmlr, 2017. 11

[4] Saleh Ashkboos, Amirkeivan Mohtashami, Maximilian L Croci, Bo Li, Pashmina Cameron, Martin Jaggi, Dan Alistarh, Torsten Hoefler, and James Hensman. Quarot: Outlier-free 4-bit inference in rotated llms. Advances in Neural Information Processing Systems, 37:100213– 100240, 2024. [5] Ran B Basat, Yaniv Ben-Itzhak, Michael Mitzenmacher, and Shay Vargaftik. Optimal and approximate adaptive stochastic quantization. Advances in Neural Information Processing Systems, 37:94265–94291, 2024. [6] Ran Ben-Basat, Yaniv Ben-Itzhak, Gal Mendelson, Michael Mitzenmacher, Amit Portnoy, and Shay Vargaftik. A note on turboquant and the earlier drive/eden line of work. arXiv preprint arXiv:2604.18555, 2026. [7] Ran Ben Basat, Yaniv Ben-Itzhak, Michael Mitzenmacher, and Shay Vargaftik. Better than optimal: Improving adaptive stochastic quantization using shared randomness. Proceedings of the ACM on Measurement and Analysis of Computing Systems, 9(3):1–44, 2025. [8] Ran Ben-Basat, Michael Mitzenmacher, and Shay Vargaftik. How to send a real number using a single bit (and some shared randomness). In ICALP 2021, pages 1–20. Schloss DagstuhlLeibniz-Zentrum für Informatik, 2021. [9] Ran Ben-Basat, Amit Portnoy, Gil Einziger, Yaniv Ben-Itzhak, Michael Mitzenmacher, et al. Accelerating federated learning with quick distributed mean estimation. In Forty-first International Conference on Machine Learning, 2024. [10] Yaniv Ben-Itzhak, Helen Möllering, Benny Pinkas, Thomas Schneider, Ajith Suresh, Oleksandr Tkachenko, Shay Vargaftik, Christian Weinert, Hossein Yalame, and Avishay Yanai. Scionfl: Efficient and robust secure quantized aggregation. In 2024 IEEE Conference on Secure and Trustworthy Machine Learning (SaTML), pages 490–511. IEEE, 2024. [11] Vidmantas Kastytis Bentkus and Dainius Dzindzalieta. A tight Gaussian bound for weighted sums of Rademacher random variables. Bernoulli, 21(2):1231 – 1237, 2015. [12] Sebastian Caldas, Jakub Konečny, H Brendan McMahan, and Ameet Talwalkar. Expanding the reach of federated learning by reducing client resource requirements. arXiv preprint arXiv:1812.07210, 2018. [13] Louis HY Chen. Stein’s method of normal approximation: Some recollections and reflections. The Annals of Statistics, 49(4):1850–1863, 2021. [14] Louis HY Chen, Larry Goldstein, and Qi-Man Shao. Normal approximation by Stein’s method. Springer Science & Business Media, 2010. [15] Xiaoqi Chen, Shay Vargaftik, and Ran Ben Basat. When ml training cuts through congestion: Just-in-time gradient compression via packet trimming. In Proceedings of the 23rd ACM Workshop on Hot Topics in Networks, pages 177–185, 2024. [16] Tim Dettmers, Mike Lewis, Younes Belkada, and Luke Zettlemoyer. Gpt3. int8 (): 8-bit matrix multiplication for transformers at scale. Advances in neural information processing systems, 35:30318–30332, 2022. [17] Tim Dettmers, Ruslan Svirschevski, Vage Egiazarian, Denis Kuznedelev, Elias Frantar, Saleh Ashkboos, Alexander Borzunov, Torsten Hoefler, and Dan Alistarh. Spqr: A sparse-quantized representation for near-lossless llm weight compression. arXiv preprint arXiv:2306.03078, 2023. [18] Ron Dorfman, Shay Vargaftik, Yaniv Ben-Itzhak, and Kfir Yehuda Levy. Docofl: Downlink compression for cross-device federated learning. In International Conference on Machine Learning, pages 8356–8388. PMLR, 2023. [19] Jianyang Gao, Yutong Gou, Yuexuan Xu, Yongyi Yang, Cheng Long, and Raymond ChiWing Wong. Practical and asymptotically optimal quantization of high-dimensional vectors in euclidean space for approximate nearest neighbor search. Proceedings of the ACM on Management of Data, 3(3):1–26, 2025. 12

[20] Jianyang Gao and Cheng Long. Rabitq: Quantizing high-dimensional vectors with a theoretical error bound for approximate nearest neighbor search. Proceedings of the ACM on Management of Data, 2(3):1–27, 2024. [21] Wenchen Han, Shay Vargaftik, Michael Mitzenmacher, and Ran Ben Basat. Dynamiq: Accelerating gradient synchronization using compressed multi-hop all-reduce. arXiv preprint arXiv:2602.08923, 2026. [22] Wenchen Han, Shay Vargaftik, Michael Mitzenmacher, Brad Karp, and Ran Ben Basat. Beyond throughput and compression ratios: Towards high end-to-end utility of gradient compression. In Proceedings of the 23rd ACM Workshop on Hot Topics in Networks, pages 186–194, 2024. [23] Kathy Horadam. Hadamard matrices and their applications. Princeton university press, 2012. [24] Jakub Konečnỳ, H Brendan McMahan, Felix X Yu, Peter Richtárik, Ananda Theertha Suresh, and Dave Bacon. Federated learning: Strategies for improving communication efficiency. arXiv preprint arXiv:1610.05492, 2016. [25] Changhun Lee, Jungyu Jin, Taesu Kim, Hyungjun Kim, and Eunhyeok Park. Owq: Outlieraware weight quantization for efficient fine-tuning and inference of large language models. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38, pages 13355–13364, 2024. [26] Minghao Li, Ran Ben Basat, Shay Vargaftik, ChonLam Lao, Kevin Xu, Xinran Tang, Michael Mitzenmacher, and Minlan Yu. THC: Accelerating Distributed Deep Learning Using Tensor Homomorphic Compression. In USENIX Symposium on Networked Systems Design and Implementation, 2024. [27] Vladimir Malinovskii, Andrei Panferov, Ivan Ilin, Han Guo, Peter Richtárik, and Dan Alistarh. Higgs: Pushing the limits of large language model quantization via the linearity theorem. In Proceedings of the 2025 Conference of the Nations of the Americas Chapter of the Association for Computational Linguistics: Human Language Technologies (Volume 1: Long Papers), pages 10857–10886, 2025. [28] Or Ordentlich and Yury Polyanskiy. Optimal quantization for matrix multiplication. IEEE Transactions on Information Theory, 2025. [29] Andrei Panferov, Erik Schultheis, Soroush Tabesh, and Dan Alistarh. Quartet ii: Accurate llm pre-training in nvfp4 by improved unbiased gradient estimation. arXiv preprint arXiv:2601.22813, 2026. [30] Martin Raič. A multivariate central limit theorem for lipschitz and smooth test functions. arXiv preprint arXiv:1812.08268, 2018. [31] Ananda Theertha, X Felix, H Brendan McMahan Yu, Sanjiv Kumar, et al. Distributed mean estimation with limited communication. In International Conference on Machine Learning, 2017. [32] Albert Tseng, Jerry Chee, Qingyao Sun, Volodymyr Kuleshov, and Christopher De Sa. Quip#: Even better llm quantization with hadamard incoherence and lattice codebooks. Proceedings of machine learning research, 235:48630, 2024. [33] Albert Tseng, Qingyao Sun, David Hou, and Christopher De. Qtip: Quantization with trellises and incoherence processing. Advances in Neural Information Processing Systems, 37:59597– 59620, 2024. [34] Ilya Tyurin. New estimates of the convergence rate in the lyapunov theorem. arXiv preprint arXiv:0912.0726, 2009. [35] Shay Vargaftik, Ran Ben Basat, Amit Portnoy, Gal Mendelson, Yaniv Ben Itzhak, and Michael Mitzenmacher. Eden: Communication-efficient and robust distributed mean estimation for federated learning. In International Conference on Machine Learning, pages 21984–22014. PMLR, 2022. 13

[36] Shay Vargaftik, Ran Ben-Basat, Amit Portnoy, Gal Mendelson, Yaniv Ben-Itzhak, and Michael Mitzenmacher. Drive: One-bit distributed mean estimation. Advances in Neural Information Processing Systems, 34:362–377, 2021. [37] Cédric Villani et al. Optimal transport: old and new, volume 338. Springer, 2009. [38] Ertza Warraich, Omer Shabtai, Khalid Manaa, Shay Vargaftik, Yonatan Piasetzky, Matty Kadosh, Lalith Suresh, and Muhammad Shahbaz. {OptiReduce}: Resilient and {TailOptimal}{AllReduce} for distributed deep learning in the cloud. In 22nd USENIX Symposium on Networked Systems Design and Implementation (NSDI 25), pages 685–703, 2025. [39] Yongyi Yang, Jianyang Gao, and Wei Hu. Raana: A fast, flexible, and data-efficient post-training quantization algorithm. arXiv preprint arXiv:2504.03717, 2025. [40] Amir Zandieh, Majid Daliri, Majid Hadian, and Vahab Mirrokni. Turboquant: Online vector quantization with near-optimal distortion rate. arXiv preprint arXiv:2504.19874, 2025. [41] Zeyu Zhang, Haiying Shen, Shay Vargaftik, Ran Ben Basat, Michael Mitzenmacher, and Minlan Yu. Hack: Homomorphic acceleration via compression of the key-value cache for disaggregated llm inference. In Proceedings of the ACM SIGCOMM 2025 Conference, pages 1245–1247, 2025.

A

Deferred analysis for Section 4.2

In this appendix, we first establish the asymptotic expansion of the scaling constant cd used by the unbiased DRIVE algorithm, and then provide the formal proof of Theorem 3. A.1

Asymptotic expansion of the scaling constant cd .

We now formalize the convergence gap between the scaling constant cd and the standard Gaussian q absolute expectation c =

2 π , where g ∼ N (0, 1).

For a URR matrix RU , define x eU ≜ RU x e. Then, x eU ∈ S d−1 is a point on the unit sphere, distributed uniformly at random. By [36, Lemma 9], its expected L1 norm evaluates via the Beta function to: E[∥e xU ∥1 ] =

2d
. (d − 1) · B 21 , d−1 2

xU ∥1 ], we rewrite the Beta function in terms of the Gamma Substituting this into cd = √1d E[∥e √ Γ(x)Γ(y) function, B(x, y) = Γ(x+y) . Evaluating Γ(1/2) = π yields: √ Γ(d/2) 2 d √ cd = . d − 1 πΓ((d − 1)/2) Applying the Gamma function recurrence relation Γ(z + 1) = zΓ(z) to the denominator provides (d − 1)Γ((d − 1)/2) = 2Γ((d + 1)/2). Substituting this simplifies the expression for cd to: r d Γ(d/2) cd = . (5) π Γ((d + 1)/2) To bound the resulting ratio of Gamma functions, we apply the Tricomi-Erdélyi expansion [1, Eq. 6.1.47, page 257]. For z → ∞ and bounded constants a and b, the expansion is stated strictly using Big-O notation as:    Γ(z + a) (a − b)(a + b − 1) 1 = z a−b 1 + +O . Γ(z + b) 2z z2 14

Letting z = d/2, a = 0, and b = 1/2, we obtain:  −1/2    d (−1/2)(−1/2) 1 Γ(d/2) = 1+ +O Γ((d + 1)/2) 2 2(d/2) d2 r    2 1 1 = . 1+ +O d 4d d2 Applying this exact expansion back into Equation (5) gives: r r    r    d 2 2 1 1 1 1 cd = = . 1+ +O 1+ +O 2 π d 4d d π 4d d2 Finally, using c =

p 2/π: r

|c − cd | =

2 − π

r

2 π



1 1+ +O 4d



1 d2

 =

1 √

2d 2π

 +O

1 d2

 .

Thus, the scaling discrepancy |c − cd | is bounded by O(d−1 ). A.2

Proof of Theorem 3

Proof. For clarity, we divide the proof into five steps. q Step 1: Proxy scale closeness. Let c = π2 . By Theorem 2 and the fact that the absolute value function is 1-Lipschitz, the expected absolute value of any individual 2-RHT coordinate differs from c by at most O(d−1/2 ). Recall that, for any arbitrary unit vector v ∈ S d−1 , all d coordinates of R2 v, share the same marginal distribution. Thus, its scaled expected L1 norm evaluates exactly to this single-coordinate absolute expectation. Thus, we establish a universal uniform bound δ defining the closeness of the true 2-RHT expectation to our Gaussian proxy value c: 1 √ E[∥R2 v∥1 ] − c = O(d−1/2 ). d v∈S d−1

δ ≜ sup

(6)

This guarantees that for any unit direction, the expected L1 norm after 2-RHT is in [c − δ, c + δ]. Step 2: The parallel component. We decompose the expected direction µ = √1d E[R2−1 sign(R2 x e)] into its parallel projection along x e and an orthogonal error vector µ⊥ . For the parallel component, using R2−1 = R2T , the identity ⟨R2⊤ a, b⟩ = ⟨a, R2 b⟩ and the property ⟨sign(y), y⟩ = ∥y∥1 , we have: 1 1 ⟨µ, x e⟩ = √ E[⟨sign(R2 x e), R2 x e⟩] = √ E[∥R2 x e∥1 ] . d d By applying Equation (6), the parallel projection is bounded, that is: ⟨µ, x e⟩ ∈ [c − δ, c + δ]. Step 3: The orthogonal component. For the orthogonal component µ⊥ , let u be the unit vector pointing in its direction, ensuring ⟨e x, u⟩ = 0 and ⟨µ, u⟩ = ∥µ⊥ ∥2 . Let t > 0 be a positive scalar and consider the perturbed vector R2 (e x + tu). The absolute value inequality dictates, |A + B| ≥ |A| + sign(A) · B . Applying this pointwise to every coordinate of the perturbed vector and summing them yields: ∥R2 (e x + tu)∥1 ≥ ∥R2 x e∥1 + ⟨sign(R2 x e), R2 tu⟩ = ∥R2 x e∥1 + t⟨sign(R2 x e), R2 u⟩. Taking expectations, multiplying by √1d , and applying the identity, √ E[⟨sign(R2 x e), R2 u⟩] = E[⟨R2⊤ sign(R2 x e), u⟩] = ⟨ dµ, u⟩ , bounds the orthogonal magnitude: 1 1 √ E[∥R2 (e x + tu)∥1 ] ≥ √ E[∥R2 x e∥1 ] + t⟨µ, u⟩ . d d 15

(7)

√ By the Pythagorean theorem, the length of the perturbed input vector is ∥e x + tu∥2 = 1 + t2 . Because the L1 norm √ expectation scales linearly with vector length, applying (6) upper-bounds the left side by (c + δ) 1 + t2 . We lower-bound the right side using c − δ + t∥µ⊥ ∥2 . Substituting these limits into (7) yields: p (c + δ) 1 + t2 ≥ c − δ + t∥µ⊥ ∥2 . √ 2 Using the standard scalar upper bound 1 + t2 ≤ 1 + t2 , expanding the left side, and subtracting c − δ results in: c+δ 2 2δ + t ≥ t∥µ⊥ ∥2 . 2 q δ If ∥µ⊥ ∥2 = 0, the bound holds trivially. Otherwise, dividing by t > 0 and setting t = 2 c+δ yields: p √ ∥µ⊥ ∥2 ≤ 2 δ(c + δ) = O( δ) = O(d−1/4 ). Step 4: Bridging to the algorithm’s scale. By the Pythagorean theorem, the squared Euclidean distance from µ to the proxy target ce x is bounded by: ∥µ − ce x∥22 = (⟨µ, x e⟩ − c)2 + ∥µ⊥ ∥22 ≤ δ 2 + O(δ) = O(δ) = O(d−1/2 ). The algorithm, however, uses cd . As derived in Appendix A.1, cd converges to c by |cd −c| = O(d−1 ). Thus, applying the triangle inequality yields: ∥µ − cd x e∥2 ≤ ∥µ − ce x∥2 + ∥(c − cd )e x∥2 ≤ O(d−1/4 ) + O(d−1 ) = O(d−1/4 ). Squaring this yields ∥µ − cd x e∥22 ≤ O(d−1/2 ). Recall that ∥B(x)∥22 =

∥x∥22 ∥µ − cd x e∥22 . Dividing by ∥x∥22 proves Equation (i): c2d

1 ∥B(x)∥22 = 2 ∥µ − cd x e∥22 ≤ O(d−1/2 ) . ∥x∥22 cd E[∥x̂∥2 ]−∥E[x̂]∥2

2 2 Step 5: Computing the variance. The relative variance we evaluate is Var(x̂) = . ∥x∥22 ∥x∥22 Because the 2-RHT matrix preserves lengths√due to its orthogonality, the uncentered ℓ2 norm of any d-dimensional ±1 sign vector is identically d. Thus, the uncentered squared norm is deterministic:  2 ∥x∥2 ∥x∥2 √ E[∥x̂∥22 ] = S 2 d = d = 22. cd cd d

Dividing by ∥x∥22 , the relative uncentered squared norm is exactly 1/c2d . For the squared mean, we ∥E[x̂]∥2

evaluate ∥x∥2 2 = c12 ∥µ∥22 . Since ∥µ∥22 = ⟨µ, x e⟩2 +∥µ⊥ ∥22 = (c±δ)2 +O(d−1/2 ) = c2 ±O(d−1/2 ), 2 d we subtract this from the relative uncentered norm: 1 Var(x̂) 1 1 − c2 = 2 − 2 (c2 ± O(d−1/2 )) = ± O(d−1/2 ). 2 ∥x∥2 cd cd c2d As established, cd = c + O(d−1 ), which directly implies c12 = c12 + O(d−1 ). We substitute this into d p the leading term, absorbing the O(d−1 ) difference into the O(d−1/2 ) error. Plugging in c = 2/π 2 1 π evaluates to 1−c c2 = c2 − 1 = 2 − 1, proving Equation (ii): Var(x̂) π = − 1 + O(d−1/2 ) . ∥x∥22 2

B

Deferred analysis for Section 5.2

Proof of Theorem 5 √ Proof. Let FU (z) = P(U (x) ≤ z) and FU (z − ) = P(U (x) < z) where U (x) = d · (R2 x e)1 be the CDF of the 2-RHT scaled coordinate, and let FG (z) = P(G ≤ z) = P(G < z) be the CDF 16

of the standard Gaussian G ∼ N (0, 1). Let de denote the signed Stieltjes measure induced by the right-continuous bounded-variation function e. We formulate the difference in the expected quantization error given by e(z) = E[(ẑ − z)2 | z] between U (x) and G using the Lebesgue-Stieltjes integral over the real line: Z ∞ Z ∞ E[e(U (x))] − E[e(G)] = e(z) dFU (z) − e(z) dFG (z) −∞ −∞ Z ∞
= e(z) d FU (z) − FG (z) . −∞

To isolate the CDF difference, we apply integration by parts for Lebesgue-Stieltjes integrals: Z ∞ Z ∞
h
i ∞
FU (z − ) − FG (z − ) de(z) . e(z) d FU (z) − FG (z) = e(z) FU (z) − FG (z) − −∞

−∞

−∞

We first evaluate the boundary terms. By the definition of Bounded Support Quantization, coordinates outside [−tp , tp ] are sent without quantization error, meaning e(z) = 0 for all |z| > tp . Consequently, the boundary evaluations vanish:

lim e(z) FU (z) − FG (z) = 0 and lim e(z) FU (z) − FG (z) = 0 . z→∞

z→−∞

This leaves only the integral term. Taking the absolute value yields: Z ∞
E[e(U (x))] − E[e(G)] = − FU (z − ) − FG (z − ) de(z) −∞ Z ∞ ≤ FU (z − ) − FG (z − ) |de|(z) . −∞

By definition, the Kolmogorov distance uniformly bounds the absolute difference between the CDFs at any real point: supz |FU (z)−FG (z)| = dK (U (x), G). The same bound holds for FU (z − )−FG (z − ). Therefore, we pull this supremum outside the integral: Z ∞ Z ∞ FU (z − ) − FG (z − ) |de|(z) ≤ dK (U (x), G) |de|(z) . −∞

−∞

R∞

The remaining integral term, −∞ |de|(z), corresponds to the Total Variation of the error function e(z), denoted T V (e). To properly account for the jump discontinuities at ±tp , the Total Variation is defined by the supremum over partitions: T V (e) =

sup z0 <···<zn

n X

e(zj ) − e(zj−1 ) .

j=1

If e(z) is piecewise absolutely continuous, this equals the sum of the integrals of |e′ (z)| on the smooth pieces plus the magnitudes of all jump discontinuities. Substituting T V (e) and applying the 2-RHT Kolmogorov bound derived in Theorem 1 gives the desired bound: 1.28 · T V (e) √ . E[e(U (x))] − E[e(G)] ≤ T V (e) dK (U (x), G) ≤ d

C

Deferred analysis for Section 6

C.1

Proof of Lemma 3

Proof. For clarity, we divide the proof into six steps. √ √ e 2 HD e 1x Step 1: Definitions. Let x e ∈ S d−1 e = dHD e. Define the intermediate √ and U = dR2 x e 1x e = H, we can write U = HD2 a. Conditioned on a, the vector a = HD e. Recalling dH Pd−1 randomness of D1 is fixed, and the i-th coordinate of U is given by Ui = r=0 Hir εr ar . 17

Step 2: Formulating the covariance. We seek the conditional covariance between two distinct coordinates Ui and Uj (i ̸= j) over the randomness of D2 . Since E[εr ] = 0, the expected value E[Ui | a] = 0, making the covariance equal to E[Ui Uj | a]: ! !# " X X Cov(Ui , Uj | a) = E[Ui Uj | a] = E Hir εr ar Hjs εs as . r

s

Because the signs ε are independent Rademacher variables, E[εr εs ] = 0 for r ̸= s, and E[ε2r ] = 1. The double sum thereby collapses to: d−1 X

Cov(Ui , Uj | a) =

Hir Hjr a2r .

r=0

Step 3: The XOR property. A fundamental property of the Sylvester-Hadamard matrix (using 0-based indexing) is that the product of two elements in the same column is given by the bitwise XOR of their row indices: Hir Hjr = Hi⊕j,r . Substituting this into our covariance yields: Cov(Ui , Uj | a) =

d−1 X

Hi⊕j,r a2r .

r=0

Step 4: Expanding the intermediate vector a. Recall ar = √1d

Pd−1

el . Squaring yields: l=0 Hrl sl x

d−1 d−1

a2r =

1XX Hrl Hrm sl sm x el xm . d m=0 l=0

Separating the diagonal (l = m) from the cross-terms (l ̸= m) and applying the XOR property Hrl Hrm = Hr,l⊕m , we get: d−1

a2r =

1X 2 1X x el + Hr,l⊕m sl sm x el xm . d d l=0

l̸=m

Step 5: Canceling the independent sum. We substitute a2r back into the covariance equation from Step 3. For the first term, non-zero row of a P because i ̸= j, we know i ⊕ j ̸= 0. The sum of anyP x e2l term. For the Hadamard matrix is 0 ( r Hi⊕j,r = 0), perfectly eliminating the independent d1 cross-terms, swapping the order of summation and using that H is symmetric yields: ! d−1 X 1X Cov(Ui , Uj | a) = sl sm x el xm Hi⊕j,r Hl⊕m,r . d r=0 l̸=m

The inner sum is the dot product of two Hadamard rows, which equals d if i ⊕ j = l ⊕ m, and 0 otherwise. This cancels the d1 and acts as an indicator function I(i ⊕ j = l ⊕ m), simplifying the covariance to: X Cov(Ui , Uj | a) = 2 sl sm x el xm . l<m : l⊕m=i⊕j

Step 6: The sparse input counterexample. Consider i = 0, j = 1, and e0 + e 1 x e= √ ∈ S d−1 . 2 Then i ⊕ j = 1. In the covariance formula from Step 5, the only nonzero pair satisfying l < m and l ⊕ m = i ⊕ j = 1 is (l, m) = (0, 1). Therefore, CovD2 (U0 , U1 | a) = 2s0 s1 x e0 x e1 = s0 s1 ∈ {±1}. Moreover, VarD2 (U0 | a) =

d−1 X

2 2 H0r ar = ∥a∥22 = 1

r=0

and similarly VarD2 (U1 | a) = 1. Hence, the conditional Pearson correlation is exactly ±1. 18

C.2

Proof of Theorem 6

Proof. For clarity, we divide the proof into five steps. √ √ e 3 HD e 2 (HD e 1x Step 1: Reduction to the 2-RHT case. Consider U = dR3 x e = dHD e). Let e y√ = HD1 x e = R1 x e be the output of the first RHT. If we condition on y, the final two RHTs e 3 HD e 2 y match the structure of the 2-RHT transformation analyzed in Lemma 3. Therefore, dHD we can apply the previous conditional covariance formula, replacing the original input x e with y, and using s to denote the random signs of the intermediate matrix D2 : X Ci,j (y, D2 ) = 2 sl sm yl ym . l<m : l⊕m=i⊕j

Step 2: The expected covariance is zero. We analyze the expected value of this conditional covariance over the randomness of D2 . The sum iterates over indices where l < m, where l ̸= m. Because the Rademacher signs of D2 are independent and zero-mean, ED2 [sl sm ] = E[sl ]E[sm ] = 0. Due to the linearity of expectation, the expression evaluates to zero for any given y: ED2 [Ci,j (y, D2 )] = 0. Step 3: Variance over disjoint pairs. To bound the magnitude of the conditional correlation, we compute its variance over D2 . Since the mean is zero, for any given y, the variance is VarD2 (Ci,j (y, D2 )) = ED2 [Ci,j (y, D2 )2 ] . When expanding the squared sum, we generate squared terms and cross-terms. The condition l ⊕ m = i ⊕ j partitions the d indices into exactly d/2 uniquely matched disjoint pairs. Because the pairs are entirely disjoint, the cross-terms consist of independent Rademacher products which vanish under expectation. We are left with the expected values of the squared terms. Since (sl sm )2 = 1: X 2 ED2 [Ci,j (y, D2 )2 ] = 4 yl2 ym . l<m : l⊕m=i⊕j

Step 4: Bounding with the maximum norm. Let α = i ⊕ j. Since i ̸= j, we have α ̸= 0. We can replace the strict l < m summation by a sum over all l ∈ {0, . . . , d − 1}, with m = l ⊕ α (which is equivalent to l ⊕ m = α). Since α ̸= 0, the map l 7→ l ⊕ α pairs each index with a distinct index. Therefore, the sum over all l counts each unordered pair twice: d−1 X

X

2 yl2 yl⊕α =2

l=0

2 yl2 ym .

l<m: l⊕m=α

Using the expression from Step 3, α = i ⊕ j and the above yields, ED2 [Ci,j (y, D2 )2 ] = 4

X

2 yl2 ym =2

l<m: l⊕m=α

d−1 X

2 yl2 yl⊕α .

l=0

We bound this sum using the maximum squared element of y, denoted by ∥y∥2∞ : ED2 [Ci,j (y, D2 )2 ] = 2

d−1 X

2 yl2 yl⊕α ≤ 2 max(yk2 ) k

l=0

d−1 X l=0

yl2 = 2∥y∥2∞

d−1 X

yl2 .

l=0

Because the first RHT preserves the ℓ2 norm of the input vector x e ∈ S d−1 , ∥y∥22 = 1. This yields ED2 [Ci,j (y, D2 )2 ] ≤ 2∥y∥2∞ . Step 5: Bounding via the First RHT. The vulnerability of the 2-RHT was the existence of highly sparse inputs. By passing x e through the first RHT R1 , standard  sub-Gaussian union bounds (e.g., see [31, Lemma 7] or [2]) dictate that ED1 ∥R1 x e∥2∞ ≤ O logd d . Consequently, taking the expectation over the randomness of the first RHT (i.e., over D1 ) and calculating the RMS, the conditional covariance between any pair of coordinates decays: ! r

log d 1/2 1/2 ED1 ,D2 [Ci,j (R1 x e, D2 )2 ] ≤ 2 ED1 [∥R1 x e∥2∞ ] ≤O . d 19

C.3 C.3.1

Proof of Theorem 7 Main proof.

Proof. For clarity, we divide the proof into five steps. Step 1: Deconstructing the matrix multiplication. Conditioned on the intermediate vector y = R2 x e, √ e 3 y = HD3 y. It is sufficient to analyze only the the last RHT applies the transformation U = dHD first k coordinates of U , which we denote as the block U0:k−1 . Let H0:k−1 be the sub-matrix consisting of only the first k rows of the Hadamard matrix. We can write our output block as: U0:k−1 = H0:k−1 (D3 y) . The vector D3 y is simply the column vector y with its coordinates multiplied by the random Rademacher signs εj from D3 . Therefore, the j-th element of D3 y is exactly εj yj . By definition, multiplying H0:k−1 by a column vector produces a linear combination of the columns of H0:k−1 . Let H (j) denote the j-th column of H0:k−1 (which is a deterministic vector of length k). The multiplication expands to: d−1 X U0:k−1 = (εj yj )H (j) . j=0

To separate the random signs ε from the deterministic values (i.e., conditioned on y), we group the deterministic parts together to define a set of d fixed “step vectors”, denoted as Vj : Vj = yj H (j) . Notice that the i-th coordinate of this step vector is exactly (Vj )i = Hi,j yj . Substituting this definition back into our sum, we have: U0:k−1 =

d−1 X

εj Vj .

j=0

Now, calculating the conditional covariance matrix Σy of this block becomes straightforward. By ⊤ definition, Σy = E[U0:k−1 U0:k−1 | y]. Because the signs ε are independent and zero-mean, all cross-terms E[εj εm ] vanish when j ̸= m. We are left only with the sum of the outer products of the individual step vectors: d−1 d−1 X X Vj Vj⊤ . Σy = E[ε2j ]Vj Vj⊤ = j=0

j=0

Looking at the specific entries of this k × k covariance matrix: P P • Diagonal entries (Variance): For coordinate i, this is j (Vj )2i = j (Hi,j yj )2 . Since the P 2 Hadamard matrix consists of ±1, Hi,j = 1, simplifying the sum to j yj2 = ∥y∥22 = 1. P • Off-diagonal entries (Covariance): For distinct coordinates i and l, this is j (Vj )i (Vj )l = P P 2 2 j Hi,j Hl,j yj . Using the Hadamard XOR property, this becomes exactly j Hi⊕l,j yj , which matches exactly the cross-correlation equation Ci,l (R1 x e) established in Theorem 6. Step 2: Bounding the covariance difference. Write the conditioned vector as e 2 y (1) , y = R2 x e = HD

y (1) = R1 x e.

The diagonal entries of Σy are equal to 1. For i ̸= ℓ, the off-diagonal entry satisfies (Σy )iℓ =

d−1 X

Hij Hℓj yj2 = Ci,ℓ (y (1) , D2 ),

j=0

20

where Ci,ℓ is the conditional covariance from Theorem 6. Therefore, by Jensen’s inequality and Theorem 6,  1/2 X    
  1/2 E ∥Σy − Ik ∥F ≤ E ∥Σy − Ik ∥2F = E Ci,ℓ (R1 x e, D2 )2  i̸=ℓ

p 
1/2  ≤O k ≤ 2k(k − 1) E ∥R1 x e∥2∞

r

log d d

! .

Step 3: Formalizing the Wasserstein bound via Stein’s method. To bound the distance between our block U0:k−1 and a standard Gaussian Z ∼ N (0, Ik ), we use a multivariate normal approximation lemma proved in C.3.2 that is based on [30, Theorem 2.15 and Eq. (3.5)]. For a sum of independent, Pd−1 Pd−1 zero-mean random vectors W = j=0 Xj with covariance Σ = j=0 E[Xj Xj⊤ ], its 1-Wasserstein distance to a standard Gaussian is bounded by:   d−1   X   W1 W, N (0, Ik ) ≤ Ak ∥Σ − Ik ∥F + E ∥Xj ∥32  j=0

where Ak is a constant that depend on k. Pd−1 In our setting, conditioned on the intermediate vector y, our block is the sum U0:k−1 = j=0 εj Vj . The summands Xj = εj Vj are independent because the Rademacher signs εj are independent. They are zero-mean (E[εj Vj | y] = 0), and their conditional covariance is Σy as derived in Step 1. Because εj ∈ {−1, 1}, the third absolute moment of each summand simplifies to:     E ∥εj Vj ∥32 | y = ∥Vj ∥32 E |εj |3 = ∥Vj ∥32 . qP √ k−1 2 We calculate the length of each step vector as ∥Vj ∥2 = k|yj |. Summing i=0 (Hi,j yj ) = these cubes over all d steps yields: d−1 X j=0

where ρ3 (y) =

Pd−1

j=0 |yj |

3

∥Vj ∥32 =

d−1 X

k 3/2 |yj |3 = k 3/2 ρ3 (y) ,

j=0

is the third absolute moment of the intermediate vector y.

Substituting the conditional covariance error and the third moments into the above bound gives the conditional Wasserstein distance:     W1 U0:k−1 | y, N (0, Ik ) ≤ Ak ∥Σy − Ik ∥F + k 3/2 ρ3 (y) . To obtain the unconditional bound, we take the expectation over the random vector y. By the joint convexity of the Wasserstein metric, the expected distance is bounded by the expectation of  the conditional Wasserstein distances, W1 (U0:k−1 , N (0, Ik )) ≤ Ey W1 (U0:k−1 | y, N (0, Ik )) . Applying this yields:        W1 U0:k−1 , N (0, Ik ) ≤ Ak E ∥Σy − Ik ∥F + k 3/2 E ρ3 (y) . q From Step 2, we have bounded the expected covariance difference as E[∥Σy − Ik ∥F ] ≤ O(k logd d ). Furthermore, by the 1-RHT smoothing property from Lemma 1, the expected third moment of y is bounded by E[ρ3 (y)] ≤ O( √1d ). Combining these terms gives our final multi-dimensional bound: " ! ! r r  3/2 #   log d k log d W1 U0:k−1 , N (0, Ik ) ≤ Ak O k +O √ = Ok . d d d Step 4: Separating the quadratic part of the VQ error. Define the VQ error function LC (v) = minc∈C ∥v − c∥22 . We decompose this error into a quadratic term and a codebook-dependent term: 21

LC (v) = ∥v∥22 + gC (v), where gC (v) = minc∈C ∥c∥22 − 2⟨v, c⟩ . Although LC is not globally Lipschitz because of the term ∥v∥22 , the function gC is globally Lipschitz since for any v, w ∈ Rk ,

|gC (v) − gC (w)| ≤ sup ∥c∥22 − 2⟨v, c⟩ − ∥c∥22 − 2⟨w, c⟩ c∈C

= 2 sup |⟨v − w, c⟩| ≤ 2B∥v − w∥2 . c∈C

Thus, gC is 2B-Lipschitz. Step 5: Final gap assembly via quadratic cancellation. Let Z ∼ N (0, Ik ). We first show that the quadratic terms agree in expectation. Conditioned on y = e, each coordinate of the 2x  RP d−1 2 2 2 block U0:k−1 = HD3 y has conditional second moment E Ui | y = j=0 Hij yj = ∥y∥22 = 1. P k−1 Therefore, E∥U0:k−1 ∥22 = i=0 E[Ui2 ] = k. On the other hand, since Z ∼ N (0, Ik ), E∥Z∥22 = k. Hence the quadratic parts cancel: E∥U0:k−1 ∥22 − E∥Z∥22 = 0. Using the decomposition from Step 4, we now obtain h i h i E min ∥U0:k−1 − c∥22 − E min ∥Z − c∥22 = |E[gC (U0:k−1 )] − E[gC (Z)]| . c∈C

c∈C

Since gC is 2B-Lipschitz, the definition of the multivariate W1 distance yields   |E[gC (U0:k−1 )] − E[gC (Z)]| ≤ 2B W1 U0:k−1 , N (0, Ik ) . Applying the Wasserstein estimate from Step 3 yields h i h i E min ∥U0:k−1 − c∥22 − E min ∥Z − c∥22 ≤ Ok,B c∈C

c∈C

r

log d d

! .

This concludes the proof. C.3.2

Auxiliary lemma for multivariate normal approximation via Stein’s method

We use the following estimate in the proof of Theorem 7. The point of the lemma is to allow the covariance of the summands to be close to Ik , rather than exactly equal to Ik . Lemma 4 (Multivariate W1 normal approximation with covariance mismatch). Fix k ≥ 1, and k let XP 1 , . . . , Xn be independent mean-zero random vectors in R with finite third moments. Define n S = j=1 Xj and Σ = Cov(S), and let Z ∼ N (0, Ik ). Then   n X W1 (S, Z) ≤ Ak ∥Σ − Ik ∥F + E∥Xj ∥32  , j=1

√ where Ak < ∞ depends only on k. Specifically, one may take Ak = 4Ck + 4 k + 1 where Ck = 11.1 + 0.83 log k. Proof. For clarity, we divide the proof into six steps. Step 1: The identity-covariance estimate. We first present the standard covariance case that follows from [30, Theorem 2.15 and Eq. (3.5)]. P  n Let Y1 , . . . , Yn be independent mean-zero random vectors in Rk , and suppose Cov Y j=1 j = Ik . Pn k Let T = j=1 Yj . For a test function f : R → R, write M1 (f ) for its Lipschitz constant with respect to the Euclidean norm. Then, the bound by [30, Theorem 2.15 and Eq. (3.5)] yields: n X   |Ef (T ) − Ef (Z)| ≤ M1 (f ) E ∥Yj ∥22 min {4.5, Ck ∥Yj ∥2 } , j=1

where Ck = 11.1 + 0.83 log k. Since min {4.5, Ck ∥Yj ∥2 } ≤ Ck ∥Yj ∥2 , we obtain |Ef (T ) − Ef (Z)| ≤ Ck M1 (f )

n X j=1

22

E∥Yj ∥32 .

Taking the supremum over all 1-Lipschitz functions f yields: W1 (T, Z) ≤ Ck

Pn

3 j=1 E∥Yj ∥2 .

Step 2: Notation for the general covarianceP case. We now return to the original summands n X1 , . . . , Xn . Define r = ∥Σ − Ik ∥F and B = j=1 E∥Xj ∥32 . We prove the desired estimate by splitting into two cases according to the size of r. Step 3: The small covariance-mismatch case. Assume that r ≤ 12 . Since Σ is a covariance matrix, it is symmetric positive semidefinite. Let λ1 , . . . , λk be the eigenvalues of Σ. Then, the eigenvalues Pk of Σ − Ik are λ1 − 1, . . . , λk − 1. Therefore, r2 = ∥Σ − Ik ∥2F = ℓ=1 (λℓ − 1)2 . Hence, for every ℓ = 1, . . . , k, !1/2 k X 1 |λℓ − 1| ≤ (λm − 1)2 =r≤ . 2 m=1 Consequently, 12 ≤ λℓ ≤ 32 well-defined.

for all ℓ = 1, . . . , k. In particular, Σ is positive definite, so Σ−1/2 is

−1/2 ej = Σ−1/2 Xj and Se = Pn X e Define the standardized summands X S. The random j=1 j = Σ e1 , . . . , X en are independent and mean zero. Moreover, vectors X

e = Σ−1/2 Cov(S)Σ−1/2 = Σ−1/2 ΣΣ−1/2 = Ik . Cov(S) Applying the identity-covariance estimate from Step 1 gives e Z) ≤ Ck W1 (S,

n X

E∥Σ−1/2 Xj ∥32 .

j=1

We now bound the third moments after this linear change of variables. For a matrix A, define its Euclidean operator norm by ∥A∥op = sup∥u∥2 =1 ∥Au∥2 . Equivalently, ∥A∥op is the smallest number L such that ∥Av∥2 ≤ L∥v∥2 for every v ∈ Rk . Indeed, if v ̸= 0, then u = v/∥v∥2 has unit norm, v and therefore ∥Av∥2 = ∥v∥2 A ∥v∥ 2

2

≤ ∥v∥2 ∥A∥op . The case v = 0 is trivial.

Applying this with A = Σ−1/2 and v = Xj , we obtain ∥Σ−1/2 Xj ∥2 ≤ ∥Σ−1/2 ∥op ∥Xj ∥2 . Cubing both sides and taking expectations gives E∥Σ−1/2 Xj ∥32 ≤ ∥Σ−1/2 ∥3op E∥Xj ∥32 . Summing over j Pn yields j=1 E∥Σ−1/2 Xj ∥32 ≤ ∥Σ−1/2 ∥3op B. Thus, e Z) ≤ Ck ∥Σ−1/2 ∥3op B. W1 (S, Since the eigenvalues of Σ lie in [1/2, 3/2], the eigenvalues of Σ−1/2 lie in

hq

2 3,

√ i 2 . Be-

cause Σ−1/2 is√symmetric positive definite, its operator norm is its largest eigenvalue. Therefore, ∥Σ−1/2 ∥op ≤ 2. Substituting this into the preceding bound gives √ e Z) ≤ 2 2 Ck B. W1 (S, Step 4: Returning from the standardized covariance to Ik . By the triangle inequality, W1 (S, Z) ≤ W1 (S, Σ1/2 Z) + W1 (Σ1/2 Z, Z). e We also use the following Lipschitz property of W1 : For the first term, use the identity S = Σ1/2 S. for any matrix A and any random vectors P, Q in Rk , W1 (AP, AQ) ≤ ∥A∥op W1 (P, Q). To see this directly from the dual definition of W1 , let f be any 1-Lipschitz function. Define g(v) := f (Av). Then g is ∥A∥op -Lipschitz, since |g(v) − g(w)| = |f (Av) − f (Aw)| ≤ ∥A(v − w)∥2 ≤ ∥A∥op ∥v − w∥2 . 23

g(v) If ∥A∥op = 0, then the desired inequality is trivial. Otherwise, set h(v) := ∥A∥ . Then h is op 1-Lipschitz. Hence, by the dual definition of W1 ,

|Ef (AP ) − Ef (AQ)| = |Eg(P ) − Eg(Q)| = ∥A∥op |Eh(P ) − Eh(Q)| ≤ ∥A∥op W1 (P, Q). Taking the supremum over all 1-Lipschitz functions f gives W1 (AP, AQ) ≤ ∥A∥op W1 (P, Q). e and Q = Z gives Applying this property with A = Σ1/2 , P = S, e Σ1/2 Z) ≤ ∥Σ1/2 ∥op W1 (S, e Z). W1 (S, Σ1/2 Z) = W1 (Σ1/2 S, Since the eigenvalues of Σ lie in [1/2, 3/2], the eigenvalues of Σ1/2 lie in q ∥Σ1/2 ∥op ≤ 32 . Combining this with the bound from Step 3 gives r √ 3 √ 1/2 W1 (S, Σ Z) ≤ · 2 2 Ck B = 2 3 Ck B ≤ 4Ck B. 2

h

√1 , 2

q i 3 2

. Thus,

For the second term, we compare Σ1/2 Z and Z using the same Gaussian   vector Z. For any 1Lipschitz function f , Ef (Σ1/2 Z) − Ef (Z) ≤ E f (Σ1/2 Z) − f (Z) ≤ E∥(Σ1/2 − Ik )Z∥2 . Taking the supremum over all 1-Lipschitz f yields, W1 (Σ1/2 Z, Z) ≤ E∥(Σ1/2 − Ik )Z∥2 . By Cauchy–Schwarz,  1/2 . E∥(Σ1/2 − Ik )Z∥2 ≤ E∥(Σ1/2 − Ik )Z∥22 Since Z ∼ N (0, Ik ), E∥(Σ1/2 − Ik )Z∥22 = ∥Σ1/2 − Ik ∥2F . Indeed, if M = Σ1/2 − Ik , then E∥M Z∥22 = E[Z ⊤ M ⊤ M Z] = tr(M ⊤ M E[ZZ ⊤ ]) = tr(M ⊤ M ) = ∥M ∥2F . Therefore,

W1 (Σ1/2 Z, Z) ≤ ∥Σ1/2 − Ik ∥F .

It remains to compare ∥Σ1/2 − Ik ∥F with ∥Σ − Ik ∥F . Since Σ is symmetric positive √ semidefinite, it has an orthonormal eigenbasis. In that basis, the eigenvalues of Σ1/2 − Ik are λℓ − 1, while the eigenvalues of Σ − Ik are λℓ − 1. For every λℓ ≥ 0, p |λℓ − 1| | λℓ − 1| = √ ≤ |λℓ − 1|. λℓ + 1 Hence, ∥Σ1/2 − Ik ∥F ≤ ∥Σ − Ik ∥F = r. Thus, W1 (Σ1/2 Z, Z) ≤ r. Combining the two terms, in the case r ≤ 1/2 we obtain W1 (S, Z) ≤ r + 4Ck B. Step 5: The large covariance-mismatch case. Assume that r > 21 . We use a crude bound that depends only on the second moments. Let f be any 1-Lipschitz function. Subtracting the constant f (0) does not change the difference of expectations, so |Ef (S) − Ef (Z)| = |E[f (S) − f (0)] − E[f (Z) − f (0)]| ≤ E|f (S) − f (0)| + E|f (Z) − f (0)| ≤ E∥S∥2 + E∥Z∥2 . Taking the supremum over all 1-Lipschitz f gives W1 (S, Z) ≤ E∥S∥2 + E∥Z∥2 . 24

By Jensen’s inequality, E∥S∥2 ≤ E∥S∥22

1/2

and E∥Z∥2 ≤ E∥Z∥22

1/2

.

Since ES = 0 and Cov(S) = Σ, E∥S∥22 = tr(Σ). Also, since Z ∼ N (0, Ik ), E∥Z∥22 = k. Therefore, p √ W1 (S, Z) ≤ tr(Σ) + k. We now bound tr(Σ). Since tr(Σ) = k + tr(Σ − Ik ), and since tr(Σ − Ik ) = ⟨Σ − Ik , Ik√ ⟩F , Cauchy’s inequality for the Frobenius inner product yields, tr(Σ − I ) ≤ ∥Σ − I ∥ ∥I ∥ = r k. k k F k F p p √ √ √ √ 2 Thus, tr(Σ) ≤ k + r k.√Consequently, tr(Σ) ≤√ k + r k ≤ k + r, because ( k + r) = √ √ √ k + 2r k + r2 ≥ k + r k. Hence, W1 (S, Z) ≤ 2 k + r. Since r > 1/2, 2 k ≤ 4 k r. Therefore, in the case r > 1/2,

√ W1 (S, Z) ≤ (4 k + 1)r.

Step 6: Combining the two cases. If r ≤ 1/2, Step 4 yields, √ W1 (S, Z) ≤ r + 4Ck B ≤ (4Ck + 4 k + 1)(r + B). If r > 1/2, Step 5 yields, √ √ W1 (S, Z) ≤ (4 k + 1)r ≤ (4Ck + 4 k + 1)(r + B). Thus, in all cases, √

W1 (S, Z) ≤ Ak (r + B),

where Ak = 4Ck + 4 k + 1. Substituting the definitions of r and B gives   n X W1 (S, Z) ≤ Ak ∥Σ − Ik ∥F + E∥Xj ∥32  . j=1

This proves the lemma.

25