Subzero matrix completion for sparse data analysis: large-scale learning of latent low-rank structure
Authors: Lawrence K. Saul, Ningyuan Huang, Dennis Bollweg, Jeff Soules, Diana C. Halikias
Organizations: Center for Computational Mathematics, Flatiron Institute, New York, NY 10010 · Scientific Computing Core, Flatiron Institute, New York, NY 10010 · Courant Institute, NYU, New York, NY 10012
We investigate when a sparse nonnegative matrix can be recovered from a real-valued matrix of much lower rank by zeroing out its negative elements. The potential for such decompositions suggests a mathematical connection between sparsity and rank; we analyze a number of sparse matrices with this latent low-rank structure and use them to illustrate the geometric origins of this connection. Previous algorithms have discovered these decompositions via an alternating minimization over the factors of a low-rank matrix, but to do so, they have also needed to compute and store another matrix, neither sparse nor low-rank, that is the size of their product. We develop a stochastic, alternating least-squares algorithm that operates on smaller blocks of this dense matrix and scales as a result to much larger problems. We also show how to further accelerate this algorithm with sparse optimizations and customized CUDA kernels. As one example, we use the algorithm to analyze the sparse matrix of synaptic weights for the recently published Drosphilia connectome. The nonzero elements of this matrix, with 139,255 rows and columns, record the number of synapses between cells in the nervous system of a female fruit fly. Despite a slowly decaying spectrum of singular values, this matrix exhibits a latent low-rank structure that is predictive of cell categories across multiple levels of specificity.
Figures & tables
Figure 1: A sparse nonnegative matrix S of full rank ( left ) can be transformed into a matrix L of lower rank ( right ) by replacing zeros in S with negative values in L (shown in red). Since S and L share the same positive entries (shown in blue), the latter provides a low-rank encoding of the former via S=max(0,L) . In this case we say that L is a subzero matrix completion of S .
Figure 2: Graphs whose signless Laplacians have subzero matrix completions of much lower rank than the number of nodes ( left , middle ) or the number of connected components ( right ).
S
m
n
ρ+
ε10
ε100
ε1000
m+
n+
rank( S+ )
digits
70,000
70,000
0.00023
0.998
0.987
0.923
11
9
9
connectome
139,003
138,955
0.00102
0.964
0.866
0.662
61
59
59
bigrams
249,992
249,987
0.00088
0.479
0.379
0.307
218
214
214
Table 1: Sparse nonegative matrices analyzed numerically in this paper. For each matrix S , the table shows the number of rows and columns (m,n) with at least one nonzero element, the fraction of nonzero elements ρ+ , and the relative errors εk in Equation 22 from truncated singular value decompositions of rank k∈{10,100,1000} . The table also shows the number of rows and columns (m+,n+) and rank of the largest positive submatrix S+ found by a randomized greedy procedure.
Figure 3: The sparse similarity matrix for MNIST digit images in Table 1 exhibits an approximately block diagonal structure if its rows and columns are appropriately permuted. Above, shown side-by-side, are the upper and lower halves of this reconstituted matrix.
Figure 4: Roofline plots showing the theoretical peak FP32 and TF32 throughput ceilings for the H100 and RTX Pro 6000 Blackwell GPUs as a function of arithmetic intensity. The dots show the achieved performance of our fused TF32 kernels (in red) and a two-kernel epilogue-fused PyTorch baseline (in blue); see Table 2 for details.
H100
PCIe
80GB
Table 2: Achieved throughputs in TFLOP/s on H100 ( left ) and RTX GPUs ( right ) for the fused TF32 kernels versus the two-kernel PyTorch implementation with an option for TF32 acceleration.
Figure 5: Objective in Equation 26 versus the number of epochs, averaged over five different random seeds, for models of subzero matrix completion of rank 16 by SALSA and ADAM on the fly connectome matrix. SALSA converges more rapidly at the outset (and in shorter wall-clock times), while ADAM performs slightly better in the tail; see text for details.
Figure 6: Comparison of matrix approximation errors from truncated SVD and SUMAC with stochastic alternating least-squares updates. Left: normalized root-mean-squared errors from Equation 24 . Right: weighted Jaccard distances from Equation 25 . The dashed lines show the rank of the largest positive submatrix found by a randomized greedy search.
category
cardinality
labeled
missing
superclass
9
139,255
0
class
27
106,969
32,286
subclass
96
94,946
44,309
cell type
8,547
137,677
1,578
Table 3: Numbers of labeled cells in the connectome for categories of differing specificity.
SVD
1nn error
rate (%)
Table 4: 1nn error rates (%) from row and column embeddings obtained from different truncated SVDs of the connectome matrix.
SUMAC
1nn error
rate (%)
Table 5: 1nn error rates (%) from row and column embeddings obtained from different subzero matrix completions of the connectome matrix. Each subzero matrix completion was optimized for only five iterations. The table shows the mean and standard deviation of the error rates over ten different random seeds.
SUMAC
superclass
1nn error
rate (%)
Table 6: 1nn superclass error rate (%) from row and column embeddings obtained from different subzero matrix completions of the connectome matrix.
Appendix figures & tables3 assets
Supplementary material from the paper’s appendix.
Appendix
superclass
class
subclass
cell type
optic
77849
optic-lobe-intrinsic
76928
unlabeled
44309
R1-6
8525
central
32388
unlabeled
32286
transmedullary
17947
KCg-m
2190
sensory
16981
visual
11469
photo-receptor
9777
T2a
1774
visual-projection
7665
Kenyon-Cell
5177
medulla-intrinsic
7859
Tm3
1756
ascending
2362
CX
2869
lamina-monopolar
7584
T4c
1710
descending
1303
mechanosensory
2648
distal-medulla
6538
T3
1676
Appendix
Table 7: Top 10 labels per cell category in the connectome data set.
dataset
rank
step-size
iterations
mini-batches
digits
22,23,…,28
1
1000
100
connectome
24,25,…,28
1
1000
100
bigrams
24,25,…,29
0.1
4000
100
connectome
22
1
1000
25
connectome
23
1
1000
50
bigrams
22
0.01
1000
50
Appendix
Table 8: Hyperparameter settings for the models of subzero matrix completion in Section 5.2 .
Figure 7: Schematic asynchronous copy/compute pipeline overlap for the fused TF32 kernels with two staging buffers. Box widths and offsets are schematic and do not represent actual timings. For each tile, the first MMA block forms Lij , the ReLU block rectifies Lij and converts it back to TF32, and the second MMA block accumulates Yi . In the WGMMA kernel, the dashed vertical line marks completion of the stage holding tile j , after which the producer reuses that buffer for tile j+2 .
We study nonconvex methods for matrix completion, the problem of recovering a low-rank matrix from a subset of its entries. Convex methods achieve sample complexity linear in the matrix dimension and the rank, up to logarithmic factors, whereas global guarantees for commonly used nonconvex methods require a higher polynomial dependence on the rank. We close this gap by analyzing Riemannian gradient descent (RGD) and Riemannian Gauss--Newton (RGN) methods. For an n×n matrix of rank r with incoherence parameter μ and condition number κ, the two methods achieve exact recovery with high probability from O(μnrlognlog(nκ)) and O(μnrlognlog(2μrκ)) observations, respectively. The methods use a multiscale residual initialization, while the analysis simultaneously controls the spectral error and incoherence. The resulting RGD iterates converge linearly, whereas RGN eventually converges Q-quadratically.
Jian-Feng Cai, Xiliang Lu, Juntao You
Department of Mathematics, Hong Kong University of Science and Technology. · School of Mathematics and Statistics, Hubei Center for Applied Mathematics, and Hubei Key Laboratory of Computational Science, Wuhan University, Wuhan 430072, China. · School of Artificial Intelligence, Hubei Center for Applied Mathematics, and Hubei Key Laboratory of Computational Science, Wuhan University, Wuhan 430072, China.
Compressing a matrix for reconstruction still defaults to the truncated SVD, approximating the data with a single low-rank structure. It is common to reduce the residual further by adding an overlapping row-sparse component, but methods that solve this joint problem often require iterative solvers and tuning of regularization parameters. We propose the Disjoint Row-Sparse plus Low-Rank (D-SLR) decomposition, a closed-form drop-in for the truncated SVD that improves or exactly matches it. D-SLR restricts rows to either being stored verbatim or approximated by the low-rank fit, never both. Under squared error this restriction costs nothing: the joint optimum is attainable disjointly with fewer parameters at every non-trivial rank and stored row count (shape). With zero stored rows D-SLR reduces to the truncated SVD, so it never does worse at equal cost. The algorithm scores the entire error-versus-parameters tradeoff, and the solution is chosen afterwards by a supplied error target or parameter count, or by a selection rule. The grid and solution together cost three SVDs, with no tuning or regularization. We derive an assumption-free, a-posteriori lower bound on the error at every shape, giving each solution a computable certificate on the potential gain of any other choice of rank and stored rows. Experiments on synthetic and real data (LLM embedding tables, network traffic, hyperspectral images) confirm the gains and quantify the certificate.
Vincent Szolnoky
Department of Mathematical Sciences Chalmers University of Technology
Rearranging the rows or columns of a sparse matrix using an appropriate ordering can significantly reduce fill-ins, i.e., new nonzeros introduced during matrix factorization, decreasing memory usage and runtime. However, finding an ordering that minimizes fill-ins is NP-complete. Existing approaches, including graph-theoretic and deep learning methods, rely on surrogate objectives without theoretical guarantees. The Fill-Path Theorem reveals a direct and intrinsic relationship between fill-in generation and the sparse structure of the matrix as path triplet inequalities. Here we first employ a multigrid graph network to capture structural information for each vertex. We then derive a triplet sampling strategy based on inequalities. Finally, we introduce an end-max chain loss function to reduce the number of triplets whose predicted scores satisfy these inequalities. Experimental evaluations on the publicly available SuiteSparse matrix collection demonstrate the superiority of the proposed method in terms of both fill-in reduction and speedup in LU factorization time.
Ziwei Li, Tao Yuan, Fangfang Liu +3
Institute of Software, Chinese Academy of Sciences, Beijing, China