Engineering Compressed Matrix Multiplication with the Fast Walsh-Hadamard Transform

TL;DR

This work adapts Pagh's compressed matrix multiplication to use the Fast Walsh-Hadamard Transform (FWHT) instead of the FFT, preserving unbiased estimation and variance guarantees for AB while enabling a simpler, memory-efficient, multithreaded CPU implementation. Through extensive synthetic experiments on dense, heavy-hitter–dominated products, the FWHT variant achieves up to 4x faster performance than the FFT version and up to 40x speedups over DGEMM under favorable sparsity and magnitude patterns. The authors provide a free, NumPy-compatible C++ implementation with Python bindings and APptainer containerization, and demonstrate robustness of the approach across multiple instance types and problem scales. The work suggests that randomized sketching for matrix multiplication can be a practical, competitive alternative in regimes where the product is sparse or has a few dominant entries, with future work extending to precision, sparsity formats, and GPU implementations.

Abstract

We present an implementation of Pagh's compressed matrix multiplication algorithm, a randomized algorithm that constructs sketches of matrices to compute an unbiased estimate of their product. By leveraging fast polynomial multiplication via the FFT, the algorithm achieves high performance when the product matrix is sparse or contains only a small number of entries with magnitudes significantly larger than the rest. We show empirically that the algorithm is practical and can outperform state-of-the-art DGEMM implementations when the product matrix has few nonzero entries or is otherwise dominated by a small subset of elements with large magnitude. As a minor theoretical contribution, we replace the FFT with the Fast Walsh-Hadamard Transform (FWHT) in sketched multiplication, preserving all correctness and variance guarantees of the original algorithm. Experiments with our carefully engineered multithreaded CPU implementation for dense double-precision matrices on 64-core CPU nodes across a range of synthetic benchmarks, exhibiting variable sparsity patterns, show that the FWHT variant is up to 4 times faster than the FFT-based version. Under favorable sparsity and magnitude patterns in the product matrix, our FWHT-based implementation achieves a speedup of up to 40 over DGEMM from Intel MKL, with low probability of error in the estimates. Our implementation is released as free software and comes with NumPy-compatible Python bindings.

Figures (5)

• Figure 1: Variance experiment results. Note the logarithmic axes.
• Figure 2: Scaling experiment results from Vera.
• Figure 3: Speedup results from Vera.
• Figure 4: Scaling experiment results from Minerva.
• Figure 5: Speedup results from Minerva.

Theorems & Definitions (2)

• Theorem 1
• proof