Pagh’s Compressed Matrix Multiplication

• The paper introduces a novel randomized algorithm that compresses the product of two matrices using hash-and-sign sketches combined with FFT-based polynomial multiplication.
• It achieves unbiased entry estimations with variance bounded by the tail energy, ensuring robust recovery of significant entries in spiky outputs.
• The approach supports exact sparse recovery in regimes with few nonzeros and has been extended via black-box reductions to modern robust sparse multiplication methods.

Searching arXiv for Pagh's compressed matrix multiplication and closely related follow-up work. Pagh’s compressed matrix multiplication is a randomized sketch-based method for representing the product $AB$ of two real $n\times n$ matrices in compressed form, so that individual entries can be estimated, large entries can be recovered, and sparse outputs can in some regimes be obtained exactly without materializing the full dense product (Pagh, 2011). In later terminology, especially in the sparse-multiplication literature, this problem is identified with “robust sparse matrix multiplication,” because the target guarantee is not classical Frobenius- or spectral-norm approximation but an entrywise additive guarantee governed by the tail of the output after removing its $k$ largest entries (Bringmann et al., 1 Jul 2026). The line of work originates with Rasmus Pagh’s 2011 paper “Compressed Matrix Multiplication” (Pagh, 2011), which compresses the matrix product via hash-and-sign sketches and FFT-based polynomial multiplication, and it has subsequently been reinterpreted as a sparse-recovery problem, engineered with alternative transforms such as the Fast Walsh–Hadamard Transform, and generalized by black-box reductions from exact sparse matrix multiplication (Andersson et al., 14 Jan 2026).

1. Historical position and terminology

Pagh introduced compressed matrix multiplication as a method for approximating and querying matrix products when the output is sparse or compressible in Frobenius mass, with motivating applications including sample covariance matrices and basis transforms that induce sparsity (Pagh, 2011). The central objective is not generic acceleration of dense exact matrix multiplication, but output-sensitive recovery of significant entries, together with efficient sparse exact multiplication when the true product has few nonzeros (Pagh, 2011).

Subsequent literature has clarified that “compressed matrix multiplication” in Pagh’s sense is distinct from both classical approximate matrix multiplication and exact algebraic compression schemes. Bringmann, Fischer, and Nakos state that the study of the robust problem “was initiated in 2013 by Pagh,” and note that Pagh “referred to the problem as compressed matrix multiplication,” whereas they prefer the term “robust sparse matrix multiplication” to align with sparse-recovery terminology (Bringmann et al., 1 Jul 2026). This terminological shift matters because the relevant guarantee is the matrix analogue of the $\ell_\infty/\ell_2$ sparse-recovery guarantee, not the usual projection- or sampling-based approximation guarantee depending on $\|AB\|_F$ or $\|A\|_F\|B\|_F$ (Bringmann et al., 1 Jul 2026).

A separate line, “Compressed Modular Matrix Multiplication,” uses deterministic $Q$-adic packing of residues modulo a small prime into machine words and is exact rather than randomized (0803.1975). That work is algebraic digit packing, not Pagh-style sketching. This distinction is essential because both use the phrase “compressed matrix multiplication,” but the underlying models, guarantees, and arithmetic assumptions are entirely different (0803.1975).

2. Core sketching construction in Pagh’s algorithm

Pagh’s method starts from the outer-product decomposition

$AB=\sum_{k=1}^n a_k b_k,$

where $a_k$ is column $k$ of $n\times n$0 and $n\times n$1 is row $n\times n$2 of $n\times n$3 (Pagh, 2011). The product is compressed into a polynomial

$n\times n$4

using random bucket hash functions $n\times n$5 and random sign functions $n\times n$6 (Pagh, 2011). Expanding the polynomial yields

$n\times n$7

so the coefficients aggregate signed contributions of product entries that collide under the pair hash

$n\times n$8

After reduction modulo $n\times n$9, the bucket coefficients are

$k$0

and the estimator for an entry is

$k$1

This is a CountSketch-style estimator over index pairs $k$2, with random signs canceling collision noise in expectation (Pagh, 2011).

The method is computationally viable because the sketch of each outer product is a convolution of two length-$k$3 vectors, so FFT computes it in near-linear time. Pagh’s dense-case running time is $k$4, and for sparse inputs with $k$5 total nonzeros the bound is $k$6 (Pagh, 2011). This suggests that the essential abstraction is not explicit formation of $k$7, but direct computation of a compressed representation from the factors.

3. Statistical guarantees and sparse-output exactness

The basic estimator is unbiased: $k$8 and has variance bounded by

$k$9

These guarantees arise because the target contribution is always present with squared sign $\ell_\infty/\ell_2$0, while off-target collisions have mean zero under the random signs and occur with probability about $\ell_\infty/\ell_2$1 under the hash family (Pagh, 2011).

Pagh also gives a stronger skew-aware entrywise guarantee in terms of the tail energy after removing large entries. Using $\ell_\infty/\ell_2$2 for the squared Frobenius norm after zeroing the $\ell_\infty/\ell_2$3 largest-magnitude entries, the median estimator over $\ell_\infty/\ell_2$4 repetitions satisfies

$\ell_\infty/\ell_2$5

with probability $\ell_\infty/\ell_2$6 (Pagh, 2011). This is the direct precursor of later robust sparse matrix multiplication guarantees phrased as $\ell_\infty/\ell_2$7-type tail bounds.

The algorithm also admits exact recovery in sufficiently sparse regimes. If $\ell_\infty/\ell_2$8 has at most $\ell_\infty/\ell_2$9 nonzero entries and $\|AB\|_F$0, median decoding returns every entry exactly with probability $\|AB\|_F$1 (Pagh, 2011). More strongly, if $\|AB\|_F$2 and $\|AB\|_F$3 have at most $\|AB\|_F$4 nonzero entries total and $\|AB\|_F$5 has at most $\|AB\|_F$6 nonzero entries, then a sparse representation of $\|AB\|_F$7 can be computed with high probability in time $\|AB\|_F$8 (Pagh, 2011). This exact sparse-output theorem is the point at which compressed matrix multiplication already overlaps with output-sensitive sparse matrix multiplication, rather than merely approximate estimation.

4. Recovery of heavy entries and the robust interpretation

A defining feature of Pagh’s framework is that it targets recovery of significant entries rather than low-rank or norm-relative approximation of the full dense product (Pagh, 2011). To identify large entries without scanning all $\|AB\|_F$9 positions, the algorithm augments sketching with error-correcting codes. Code bits are used to create masked sketches whose bucket behavior reveals the row and column indices of isolated heavy entries; decoding yields a candidate set of $\|A\|_F\|B\|_F$0 positions containing the large-magnitude entries above an appropriate threshold (Pagh, 2011).

Bringmann, Fischer, and Nakos place this guarantee in sparse-recovery language. For sparse inputs $\|A\|_F\|B\|_F$1 with $\|A\|_F\|B\|_F$2 and parameter $\|A\|_F\|B\|_F$3, they define robust sparse matrix multiplication by the requirement

$\|A\|_F\|B\|_F$4

where $\|A\|_F\|B\|_F$5 is obtained by zeroing out the $\|A\|_F\|B\|_F$6 largest entries of $\|A\|_F\|B\|_F$7 in absolute value (Bringmann et al., 1 Jul 2026). They state explicitly that Pagh’s “compressed matrix multiplication” is this robust problem under a later name, and that among known sparse matrix multiplication algorithms, Pagh’s algorithm with running time $\|A\|_F\|B\|_F$8 had been the only one known to be robust in this sense (Bringmann et al., 1 Jul 2026).

This interpretation resolves a common misconception. Pagh’s method is not merely an approximate matrix multiplication algorithm in the classical projection or sampling sense. Its central guarantee depends only on the tail of the output, not on the full Frobenius norm of the product, which is substantially stronger when the product is “spiky” and dominated by a small number of large entries (Bringmann et al., 1 Jul 2026).

5. Black-box robustification and the modern structural view

The principal modern development is the reduction of robust sparse matrix multiplication to exact sparse matrix multiplication. Bringmann, Fischer, and Nakos prove that if sparse matrix multiplication with $\|A\|_F\|B\|_F$9 nonzero inputs and $Q$0 nonzero outputs runs in time $Q$1, then robust sparse matrix multiplication runs in time

$Q$2

with high probability (Bringmann et al., 1 Jul 2026). Their central message is that Pagh’s guarantee is not a specialized consequence of his FFT-based polynomial encoding, but a black-box consequence of exact sparse multiplication up to polylogarithmic overhead.

The reduction is built around a query model for an unknown matrix $Q$3, instantiated later with $Q$4. One can query compressed products $Q$5, where $Q$6 has column sparsity $Q$7 and $Q$8 has row sparsity $Q$9; in the multiplication setting, such a query is answered by computing $AB=\sum_{k=1}^n a_k b_k,$0, $AB=\sum_{k=1}^n a_k b_k,$1, and then $AB=\sum_{k=1}^n a_k b_k,$2 (Bringmann et al., 1 Jul 2026). Sparse-recovery sketches are then applied to recover heavy entries from a small number of such compressed products. The technical toolkit includes an $AB=\sum_{k=1}^n a_k b_k,$3-tail estimator, a fast heavy-hitter recovery sketch, and a budget-allocation mechanism that assigns per-column recovery budgets $AB=\sum_{k=1}^n a_k b_k,$4 satisfying

$AB=\sum_{k=1}^n a_k b_k,$5

with high probability (Bringmann et al., 1 Jul 2026).

A notable element is a minimum-cost multiple-choice knapsack formulation used to estimate the global tail and derive feasible column budgets from sketch data (Bringmann et al., 1 Jul 2026). This suggests a reinterpretation of Pagh’s original problem: the difficult part is not polynomial compression per se, but allocation of sparse-recovery effort across a highly unbalanced matrix product. In this view, Pagh’s algorithm is the first concrete realization of a broader robustification principle.

6. Complexity improvements and implementation refinements

Instantiating the black-box reduction with the state-of-the-art sparse matrix multiplication algorithm of Abboud, Bringmann, Fischer, and Künnemann yields

$AB=\sum_{k=1}^n a_k b_k,$6

for robust sparse matrix multiplication (Bringmann et al., 1 Jul 2026). In the regime $AB=\sum_{k=1}^n a_k b_k,$7, the same reduction gives

$AB=\sum_{k=1}^n a_k b_k,$8

for any $AB=\sum_{k=1}^n a_k b_k,$9, with high probability (Bringmann et al., 1 Jul 2026). Since $a_k$0 time is necessary merely to write down $a_k$1 significant outputs, this is essentially optimal in that regime (Bringmann et al., 1 Jul 2026). Historically, this both subsumes and improves Pagh’s $a_k$2 bound.

On the implementation side, recent work has revisited Pagh’s original transform step. “Engineering Compressed Matrix Multiplication with the Fast Walsh-Hadamard Transform” replaces the FFT-based cyclic convolution by Walsh–Hadamard-based XOR convolution, using the pair hash

$a_k$3

instead of $a_k$4, while preserving unbiasedness and variance guarantees (Andersson et al., 14 Jan 2026). The paper states that all relevant correctness results from Pagh carry over to the FWHT variant, because the analysis depends on pairwise-independent bucket assignment and random signed collisions rather than specifically on cyclic convolution (Andersson et al., 14 Jan 2026).

Empirically, this engineering study reports that the FWHT variant is up to $a_k$5 faster than the FFT-based version, and under favorable sparsity and magnitude patterns in the product matrix can achieve a speedup of up to $a_k$6 over Intel MKL’s DGEMM, while retaining low failure probability in the estimates (Andersson et al., 14 Jan 2026). The implementation targets dense double-precision matrices on 64-core CPU nodes and provides NumPy-compatible Python bindings (Andersson et al., 14 Jan 2026). These results do not alter the algorithmic essence of Pagh’s estimator, but they demonstrate that the sketch-based approach can be practical when the output is sparse or dominated by a few heavy entries.

7. Relation to neighboring compressed-multiplication paradigms

Pagh’s compressed matrix multiplication belongs to a broad family of compressed-output methods, but it occupies a specific niche. It differs from classical approximate matrix multiplication in the Cohen–Lewis, Drineas–Kannan–Mahoney, and Sarlós traditions because those methods typically target Frobenius- or spectral-norm approximation of the whole product, whereas Pagh’s method targets entrywise additive recovery of heavy entries and exact sparse-output recovery when possible (Pagh, 2011, Bringmann et al., 1 Jul 2026). It also differs from exact deterministic packing methods such as $a_k$7-adic compressed modular multiplication, which encode several residues into one machine word and remain exact provided explicit arithmetic inequalities hold (0803.1975).

It is also distinct from recent compressed-sensing formulations of output-sparse multiplication over arbitrary rings. “Output-Sparse Matrix Multiplication Using Compressed Sensing” gives deterministic and randomized algorithms that compute $a_k$8 exactly under a sparsity promise $a_k$9, using explicit compressed sensing and matrix-multiplication verification rather than Pagh’s hash-and-sign polynomial sketch (Bennett et al., 14 Aug 2025). That work is conceptually allied to the “compress first, recover later” paradigm, but technically different: it uses deterministic sparse-recovery matrices over arbitrary rings and exact recovery, whereas Pagh’s original framework is randomized, real-valued, and organized around CountSketch-style bucket collisions (Bennett et al., 14 Aug 2025).

At a broader methodological distance are matrix-compression schemes that support direct exact multiplication in a structured compressed representation, such as grammar-compressed matrix–vector multiplication or adaptive multiplication of $k$0-matrices (Ferragina et al., 2022, Börm, 2023). Those lines of work treat “compressed matrix multiplication” as arithmetic directly on compressed operands rather than sketch-based recovery of an implicitly compressed output. The comparison helps isolate what is specific about Pagh’s contribution: not compression of stored matrices, but compressed sensing of the product itself without first forming the product.

In contemporary perspective, Pagh’s compressed matrix multiplication is best understood as the foundational sketch-based formulation of robust sparse matrix multiplication: a method that compresses the product $k$1 into a small randomized summary, supports unbiased entry queries with variance $k$2, recovers heavy entries through coding-based decoding, yields exact multiplication in sparse-output regimes, and—through later work—admits a black-box equivalence to exact sparse matrix multiplication up to polylogarithmic overhead (Pagh, 2011, Bringmann et al., 1 Jul 2026).