This paper explores the generalized problem of (approximately) computing the largest output entries, with an approximation error dependent solely on the smaller entries, from the viewpoint of sparse recovery, and shows that any sparse matrix multiplication algorithm with running time T(n, m_{in}, m_{out}) can be transformed into a robust algorithm running in time O(T(n, m_{in), k)$.
Abstract
In the seminal sparse matrix multiplication problem the goal is to compute the product of two $n \times n$ matrices when the matrices are sparse, i.e., when the number of nonzeros in the input matrices $m_{in}$ and/or the number of nonzeros in the output matrix $m_{out}$ are much smaller than $n^2$. In this paper, we explore the generalized problem of (approximately) computing the $k$ largest output entries, with an approximation error dependent solely on the smaller entries -- from the viewpoint of sparse recovery, this can be seen as a robust variant of sparse matrix multiplication. Despite the substantial research dedicated to sparse matrix multiplication, almost no existing algorithms are robust in this sense. The one exception is Pagh's algorithm in time $\widetilde O(m_{in} + nk)$ [ITCS'12], and it remained open whether other algorithms can be similarly made robust. Our principal contribution is a black-box reduction from robust sparse matrix multiplication to conventional sparse matrix multiplication with only polylogarithmic overhead. Specifically, we show that any sparse matrix multiplication algorithm with running time $T(n, m_{in}, m_{out})$ can be transformed into a robust algorithm running in time $\widetilde O(T(n, m_{in}, k))$. This reduction leverages an extensive toolkit from sparse recovery, and intriguingly, also involves solving a knapsack-type problem. By plugging in the state-of-the-art algorithm for sparse matrix multiplication by Abboud, Bringmann, Fischer, and K\"unnemann [SODA'24], we achieve significantly improved bounds such as $O((m_{in} + k)^{1.346})$. Notably, in the regime where $k \geq m_{in}^{1.762}$, our reduction culminates in an almost-optimal $k^{1+o(1)}$-time algorithm.
The symmetric rank-$k$ update (SYRK), $\C = \A\A^\top$, computes the dot product between each pair of rows of $\A$, producing the Gram matrix $\C$. Its sparse variant underpins similarity search in machine learning, graph analytics, and genomics, including Jaccard similarity on datasets too large for a single node. Yet, despite the symmetry in its inputs and outputs, existing distributed sparse matrix multiplication algorithms, such as Sparse SUMMA, treat sparse SYRK as generic multiplication, leaving performance untapped. In this paper, we present distributed sparse SYRK approaches that leverage symmetry. The approach partitions the off-diagonal blocks of the output between the upper and lower triangular portions of the process grid and computes only the lower-triangular part of each diagonal block, reducing per-process communication and computation compared with state-of-the-art distributed SpGEMM. A second variant reorders communication to avoid materializing $\A^\top$. On 32 nodes of the Perlmutter supercomputer, the algorithm achieves a $2\times$ speedup over an optimized Sparse SUMMA on matrices where local multiplication dominates the runtime; the advantage narrows on communication-bound inputs, a dependence the cost model predicts from the arithmetic intensity. Our variant rectifies this and consistently achieves superior scaling at high process counts. The approach is a drop-in replacement for any application computing $\C = \A \A^\top$ via a distributed SpGEMM routine, and its triangular output can be consumed directly by subsequent operations, reducing both compute and memory footprint.
Thomas McFarland, Julian Bellavita, G. Guidi· 0 citations
We describe a simple rejection-sampling-based algorithm to perform length-squared sampling on an $n \times n$ positive-semidefinite (psd) matrix: that is, to sample a column with probability proportional to its squared $\ell_2$-norm. The algorithm runs in just $O(n)$ expected time, which is significantly sublinear in the input matrix size. The runtime is optimal, even when the input is assumed to be diagonal. Our result has several applications. Length-squared sampling is used by a number of sublinear time algorithms for matrix problems, like low-rank approximation and eigenvalue approximation. Often, it is assumed that the algorithm is given access to the matrix column norms, and thus can perform length-squared sampling efficiently. Our result shows that, at least for psd matrices, we can remove this assumption. We also discuss an application to an asymptotically optimal algorithm for estimating the Frobenius norm of a psd matrix to relative error. Finally, we show that our sampling algorithm yields a very simple sublinear time algorithm for the robust psd low-rank approximation problem introduced by Bakshi et al. (FOCS, 2020), which nearly matches the more complex method developed there.
Rajarshi Bhattacharjee, Ethan N. Epperly, Cameron Musco et al.· 0 citations
We consider the probability that a discrete random matrix $M_n(\xi)$ is \emph{strongly non-singular}, meaning all its leading principal submatrices are non-singular. This property is equivalent to the existence of an LU factorization. We show that for any discrete random variable $\xi$ with finite support and $|\xi|_\infty<1$, there is a constant probability that $M_n(\xi)$ is strongly non-singular with a growth factor bounded by $n^{5/2+\delta}$. Furthermore, we provide a tight asymptotic lower bound for this probability as $|\xi|_\infty \to 0$. Finally, we provide exact counts for strongly non-singular binary matrices up to $n=9$ and use these to derive improved upper bounds for the Bernoulli case.
S. Mateo, John Urschel, Nicholas West· 0 citations
In this work, we consider the problem of sparse signal recovery known as compressed sensing using $\ell_1$-minimization.
We show how the $\ell_1$-minimization problem (also known as basis pursuit) can be transformed into an equivalent linear programming (LP) problem, and provide a proof of the equivalence of these two problems.
We conduct an experimental comparison of modern solvers (Gurobi, HiGHS, CPLEX, and Clarabel) for solving the LP problem on test data generated according to theoretical recovery guarantees for matrices with normally distributed elements. The results show that the open-source solver Clarabel is a competitive alternative to proprietary solvers in terms of speed.
We also propose a method for verifying the uniqueness of the obtained solution using an auxiliary quadratic programming problem with a strictly convex objective function.
A geometric interpretation of the uniqueness conditions is provided, and the application of the method is demonstrated on an example of a matrix with integer elements.
Anastasiia O. Storozhenko, P. Stetsyuk· KyivAcademUs2026· 0 citations
We study random sketching matrices with Khatri-Rao structure. In particular, we consider the Khatri-Rao product (i.e., column-wise tensor product) $A_1\odot\cdots\odot A_d \in \mathbb R^{(n_1 \cdots n_d) \times m}$ of random matrices $A_i \in \mathbb R^{n_i \times m}$ whose columns are isotropic, independent and sub-Gaussian (e.g., Gaussian matrices). Khatri-Rao sketching matrices are widely applied in randomized algorithms for linear algebraic computation and data analysis, when the input data has tensor structure that allows for fast multiplication with $A_1\odot\cdots\odot A_d$. However, existing theory is not able to fully explain their performance in practice. In particular, despite significant attention, our best bounds for the important \emph{oblivious subspace embedding} property with Khatri-Rao matrices lag behind what is achievable with standard unstructured matrices. For embedding a $k$-dimensional subspace to $(1\pm \epsilon)$ error, Bujanovi\'c et al. \cite{bujanovic2025subspace} prove that sketching dimension $m = O(k^{3/2}/\epsilon^2)$ suffices in the special case of $d = 2$. Their dependence on $k$ is weaker than the tight bound of $O(k/\epsilon^2)$ known for unstructured sub-Gaussian sketching matrices. In this work, we close this gap, showing that $m = \tilde O(k/\epsilon^2)$ suffices for subspace embedding with a Khatri-Rao sketching matrix with any fixed order $d$. Our proof is simple, leveraging just two basic properties of the Khatri-Rao sketching distribution: 1) the columns of $A_1\odot\cdots\odot A_d \in \mathbb R^{(n_1 \cdots n_d) \times m}$ are independent and isotropic, and 2) each column of $A_1\odot\cdots\odot A_d \in \mathbb R^{(n_1 \cdots n_d) \times m}$ satisfies a weak Johnson-Lindenstrauss type moment property.
Across test matrices with diverse singular-value decays, the proposed methods attain ranks close to the truncated-SVD optimum while meeting the prescribed tolerance with high probability.