Sparse Approximate Matrix Multiplication (SpAMM)

• Sparse Approximate Matrix Multiplication (SpAMM) is an algorithm that approximates the product of matrices by hierarchically pruning insignificant block contributions using Frobenius norm checks.
• It achieves near-linear computational scaling by focusing on significant sub-blocks, making it highly effective in electronic structure calculations and iterative solvers.
• Optimized data structures, parallel batching, and rigorous error control enable SpAMM to balance high accuracy with improved performance on modern hardware.

Sparse Approximate Matrix Multiplication (SpAMM) refers to a class of algorithms designed to accelerate matrix–matrix multiplication by exploiting the decay, sparsity, or structure in the input matrices. Rather than pursuing sparsity in the input or output spaces through element-dropping or traditional sparse routines, SpAMM directly prunes the product space by hierarchically omitting blocks whose contributions are provably negligible. This approach enables sub-cubic, often near-linear, scaling for large problem sizes, especially for matrices exhibiting exponential or algebraic decay, as commonly encountered in electronic structure theory and other areas of computational physics and data science (Challacombe et al., 2010, Bock et al., 2012, Artemov et al., 2020, Artemov, 2019, Liu et al., 2021).

1. Fundamental Principles and Definitions

The central idea of SpAMM is to replace exact evaluation of $C = AB$ by an approximation $\widetilde{C} \approx AB$, where only significant partial products are computed. In the hierarchical SpAMM variants, matrices $A$ and $B$ are recursively partitioned (typically in a quadtree or blockwise fashion). At every level, the Frobenius norms $\|A_{\text{block}}\|_F$ and $\|B_{\text{block}}\|_F$ are computed for pairs of sub-blocks. The product $A_{\text{block}}B_{\text{block}}$ is computed only if

$$\|A_{\text{block}}\|_F \cdot \|B_{\text{block}}\|_F \geq \tau$$

where $\tau$ is a user-specified or auto-tuned threshold. Otherwise, the product subtree is pruned and the result is set to zero. This selection criterion is a direct application of submultiplicativity of the norm; it ensures rigorous a priori error control for the omitted contributions (Challacombe et al., 2010, Bock et al., 2012, Artemov et al., 2020, Artemov, 2019).

The SpAMM approach does not require the matrices to be explicitly sparse; it is effective whenever large regions of the product tensor contribute negligibly due to decay or block-structured correlations.

2. Algorithmic Structure and Implementation Strategies

A standard SpAMM implementation proceeds as follows:

1. Matrix Representation: Matrices are stored as a hierarchy (e.g., quadtree), with each node annotated by the Frobenius norm of the encompassed block. Recursive partitioning is continued until a predetermined leaf size is reached (e.g., $4 \times 4$ or $16 \times 16$) (Bock et al., 2012, Artemov et al., 2020, Liu et al., 2021).
2. Recursive Pruning: At each interior node (block), test $$\|A_{\text{block}}\|_F \cdot \|B_{\text{block}}\|_F < \tau$$. If true, return zero. Otherwise, either recurse (if the block is not a leaf), or perform an explicit dense multiplication at the leaf level.
3. Hybrid Methods: Input matrices can be pre-truncated (element-dropping or block zeroing) to further reduce the problem size, with SpAMM applied subsequently in the product space (Artemov, 2019, Artemov et al., 2020).
4. Parallelism and Batching: Surviving block products are mapped to independent tasks. Effective parallel implementations batch these into arrays for BLAS/GPU kernels, exploiting both fine- and coarse-grained parallelism (Bock et al., 2012, Liu et al., 2021, Artemov, 2019, Artemov et al., 2020).
5. Optimized Data Structures: Highly optimized SpAMM codes use Morton-order (Z-order) or other space-filling curve-based schemes for block traversal and memory locality, efficient hash tables for storing block trees, and careful kernel unrolling for SIMD/SIMT execution (Bock et al., 2012, Liu et al., 2021).

The following table summarizes the recursive SpAMM kernel:

Step	Purpose	Typical Complexity
Frobenius norm checks	Prune insignificant block products	$O$(node visits)
Subdivision	Quadtree partitioning of blocks	$O$(log $n$) depth
Leaf multiplication	Batched small dense kernel (e.g. GEMM)	$O($num. active leaves$)$

3. Theoretical Analysis and Error Bounds

The approximation error due to SpAMM can be tightly bounded in both the Frobenius and max-norms:

• Each time a block product is omitted, the norm of the discarded term is less than $\tau$.
• The total number of omitted nodes $M$ scales as $O(n\log(1/\tau))$ for exponential (or fast algebraic) decay, so the total error

$$\|\widetilde{C} - C\| \leq O\left(n\log(1/\tau)\right)\tau$$

(Challacombe et al., 2010, Bock et al., 2012, Artemov et al., 2020, Artemov, 2019).

• For each matrix entry, $$|(AB)_{ij} - \widetilde{C}_{ij}| = O(n\tau\|A\|_F\|B\|_F)$$ for normalized matrices (Challacombe et al., 2015).
• In practical electronic structure calculations, the perturbation in the computed observable (e.g., total energy) can be rigorously bounded by the norm of the matrix error.
• For matrices with metric decay, the set of surviving products concentrates near low-dimensional manifolds (e.g., diagonals), further reducing effective complexity ("lensing" effect under iterative contractions) (Challacombe et al., 2015).
• For product sketches (low-dimensional approximations), the sparse co-occurring directions (SCOD) algorithm achieves high-accuracy operator-norm or Frobenius-norm error with optimal storage and near-linear time complexity for large, sparse matrices (Wan et al., 2020).

4. Computational Complexity and Performance Characteristics

• Dense Case: Classical matrix multiplication costs $O(n^3)$. For matrices with decay, SpAMM can achieve $O(n\log n)$ scaling in both arithmetic and data movement, dramatically reducing runtime for sufficiently large $n$ (Bock et al., 2012, Challacombe et al., 2010).
• Sparse and Blocked Cases: The actual speedup increases with sparsity/decay rate and is highest for systems with localized interactions (e.g., insulators in electronic structure) (Challacombe et al., 2010, Bock et al., 2012, Artemov, 2019, Artemov et al., 2020, Liu et al., 2021).
• Practical Benchmarks: On quantum chemistry matrices, SpAMM with single-precision achieves wall-clock speedups of $2$–$3\times$ over vendor-tuned DGEMM/SGEMM for $n\gtrsim 1000$ and can deliver better than native single-precision accuracy for moderate thresholds ($\tau\lesssim 2 \times 10^{-8}$) (Bock et al., 2012, Challacombe et al., 2010). On distributed-memory architectures, hybrid SpAMM reduces communication volume nearly $2\times$ relative to pure product-space or vector-space truncation (Artemov, 2019).
• GPU Acceleration: Implementations such as cuSpAMM achieve up to $10$–$50\times$ speedup over optimized cuBLAS/cuSPARSE for near-sparse, large matrices by offloading recursive pruning and block multiplies to multiple GPUs with fine-tuned data movement, blocking, and tensor core utilization (Liu et al., 2021). Scaling is linear to logarithmic with matrix size for decay-dominated problems.

5. Applications and Practical Use Cases

SpAMM is primarily deployed for large-scale scientific computing tasks where matrices exhibit decay:

• Electronic Structure Calculations: In Hartree–Fock, Kohn–Sham DFT, and density-matrix purification schemes, SpAMM accelerates the computation of Fock, overlap, and density matrices, reducing both arithmetic and memory requirements for large systems (e.g., proteins, water clusters, nanostructures) (Challacombe et al., 2010, Bock et al., 2012, Artemov et al., 2020, Artemov, 2019).
• Iterative Linear Algebra: SpAMM is effective in iterative square root, inverse square root, and $n$-body solvers, leveraging emergent “lensing” effects that concentrate computation along low-dimensional structures during contractive iterations (Challacombe et al., 2015).
• Task-based Parallel Frameworks: Integration into models such as Chunks & Tasks enables dynamic load balancing, fine-grain parallelism, and scalable weak scaling on distributed-memory clusters (Artemov, 2019, Artemov et al., 2020).
• Machine Learning: Efficient approximate product computation for kernel methods, neural network layers with weight decay, and streaming scenarios (via sparse SCOD) (Wan et al., 2020).

6. Variants, Extensions, and Comparisons

• Rank-One Sparsification: The Belabbas–Wolfe SpAMM variant constructs the $k$-sparse approximation to $AB$ by selecting and optimally reweighting rank-one terms, using a joint-Gram matrix and Schur complements for error analysis; greedy and randomized subset selection heuristics achieve near-optimal error decay (0707.4448).
• SCOD/Sketching Methods: Sparse co-occurring directions leverages randomized low-rank sketching to approximate $X Y^T$ within operator-norm error, exploiting input sparsity for fast, low-memory streaming and sketching applications, with strictly controlled errors (Wan et al., 2020).
• Hybrid Truncation: Combining vector-space and product-space truncations yields the best trade-off in terms of arithmetic cost, memory usage, and communication for large-scale parallel environments. The hybrid approach achieves error bounds and performance superior to either scheme individually (Artemov, 2019, Artemov et al., 2020).
• Comparison with Element Dropping: Pure element-dropping achieves $O(\varepsilon^2)$ error per pruned entry but retains more block products, leading to higher flop counts than SpAMM for matched global error tolerances (Challacombe et al., 2010, Bock et al., 2012).
• Parallelization and GPU Realization: Flat and hierarchical parallelization schemes, multi-level blocking, and fine-grained pruning ensure high utilization across diverse hardware, from CPUs to multi-GPU nodes; e.g., cuSpAMM achieves $5$–$13\times$ speedups over cuBLAS at various valid-ratios, with even greater benefits at large $n$ and tighter pruning (Liu et al., 2021).

7. Practical Considerations and Parameter Selection

• Threshold Selection: Choice of $\tau$ governs the trade-off between accuracy and computational savings. In practice, $\tau$ is chosen to achieve a target error, informed by the decay profile and system size. Schemes such as parameter sweeps (as in the CSE algorithm) rigorously bound all sources of omission error for full control (Artemov et al., 2020).
• Leaf Block Size: Performance is sensitive to the block size at which dense kernels are invoked; $16\times16$ or $4\times4$ blocks are common choices, balancing granularity and hardware efficiency (Bock et al., 2012, Challacombe et al., 2010).
• Scalability: Near-linear scaling with system size is observed for matrices with strong decay, both in arithmetic and memory requirements (Challacombe et al., 2010, Bock et al., 2012, Artemov, 2019), though communication can be a bottleneck in pure product-space truncation.
• Decaying and Ill-conditioned Systems: Regularization (e.g., Tikhonov shift) can be used to restore rapid decay in systems where it is lost due to poor conditioning, enabling SpAMM to remain effective via a nested “scoping product” representation (Challacombe et al., 2015).
• Empirical Choices: Tables of chosen $\tau$ versus target errors for various chemical basis sets, empirical timing, and error curves validate the algorithms' efficiency and reliability across real-world applications (Bock et al., 2012, Challacombe et al., 2010, Artemov et al., 2020).

• “Fast Multiplication of Matrices with Decay” (Challacombe et al., 2010)
• “An Optimized Sparse Approximate Matrix Multiply for Matrices with Decay” (Bock et al., 2012)
• “Approximate multiplication of nearly sparse matrices with decay in a fully recursive distributed task-based parallel framework” (Artemov, 2019)
• “Sparse approximate matrix-matrix multiplication for density matrix purification with error control” (Artemov et al., 2020)
• “Accelerating Sparse Approximate Matrix Multiplication on GPUs” (Liu et al., 2021)
• “A $N$-Body Solver for Square Root Iteration” (Challacombe et al., 2015)
• “On sparse representations of linear operators and the approximation of matrix products” (0707.4448)
• “Approximate Multiplication of Sparse Matrices with Limited Space” (Wan et al., 2020)