PaperlayerView exact revision on arXiv

Source-linked AI summary

Subzero matrix completion for sparse data analysis: large-scale learning of latent low-rank structure

Lawrence K. Saul, Ningyuan Huang, Dennis Bollweg, Jeff Soules, Diana C. Halikias

arXiv:2608.21607v1cs.LGstat.ML

TL;DR

The paper asks when sparse nonnegative matrices can be represented by much lower-rank real-valued matrices after negative entries are zeroed. It develops a stochastic block-based alternating least-squares method with sparse and CUDA optimizations, and applies it to large sparse datasets including a 139,255-by-139,255 fruit-fly connectome. The resulting approach reveals latent low-rank structure and supports cell-category prediction, while its stochastic updates can diverge for poorly chosen batch or step sizes.

  • Problem

    Sparse nonnegative matrices may contain latent low-rank structure, but prior algorithms require a dense, full-size auxiliary matrix and nonlinear decompositions are harder to interpret than linear ones.

  • Method

    The paper develops geometric analyses and a stochastic alternating least-squares algorithm that learns low-rank factors from smaller blocks, with sparse optimizations and customized CUDA kernels.

  • Results

    Customized CUDA kernels yield speedups of up to 4x in Python and over 10x in MATLAB, while the fruit-fly connectome exhibits latent low-rank structure predictive of cell categories.

  • Takeaways & Limitations

    Subzero matrix completion provides a scalable way to analyze latent low-rank structure in sparse matrices with hundreds of thousands of rows and columns.

  • Takeaways & Limitations

    Stochastic updates can diverge for poorly chosen batch sizes or step sizes, and the paper notes that a formal convergence proof remains desirable.

Abstract

from arXiv · show

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 $\textit{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.

1 Introduction

The paper studies a nonlinear connection between sparsity and low rank, then develops scalable algorithms for discovering latent low-rank structure in very large sparse matrices.

  • Motivation: Sparse nonnegative matrices can often be approximated as S ≈ max(0, L), where L has much lower rank and negative entries replace zeros.The shared positive entries make L a low-rank encoding of S; greater sparsity can permit more negative replacements and lower rank.
  • Motivation: Existing subzero matrix-completion algorithms use alternating minimization but generally require storing a dense, full-size auxiliary matrix.This creates a practical scaling challenge for large problems.
  • Approach: The paper combines geometric analysis of the decompositions with a stochastic alternating least-squares algorithm that operates on smaller blocks.It also provides GPU-compatible Python and MATLAB implementations with optimized CUDA kernels.
  • Applications: The method targets sparse matrices with hundreds of thousands of rows and columns, including a female fruit-fly connectome recording synaptic counts between neurons.The connectome has over 139,000 rows and columns.
  • Results: Customized CUDA kernels provide further speedups of up to 4x in Python and over 10x in MATLAB on large problems.The paper reports that stochastic optimizers are much faster than deterministic ones on problems of this size.

2 Motivation and examples

The paper studies when sparse nonnegative matrices admit much lower-rank real-valued completions after negative entries are thresholded to zero, and develops a geometric account of this relationship. Analytical constructions show when rank reduction is possible or limited, while empirical datasets exhibit slowly decaying spectra yet may retain latent low-rank structure.

  • Definitions and motivation: A subzero matrix completion represents S by thresholding a real-valued matrix L elementwise, with rank(L) prescribed to be low.The defining relation is S_ij = max(0, L_ij).
  • Geometric constructions: Thresholded similarity matrices generated by vectors in R^d have subzero completions of rank at most d+1.Nonzero entries correspond to row and column vectors whose cosine exceeds the threshold τ, with larger τ producing greater sparsity.
  • Geometric constructions: The geometric picture associates low-dimensional vectors with rows and columns, so nonzero entries encode pairs that are most nearly parallel.The converse proposition shows that a rank-r completion also induces a thresholded-similarity representation in R^(r+1).
  • Graphical and matrix examples: Permutation matrices, including the identity, have subzero matrix completions of rank 3 or less, despite potentially full rank themselves.For the identity completion, the rank gap comes with a nuclear norm that grows cubically with n.
  • Counterexample: The relative-prime matrix is a counterexample: its rank and every completion’s rank grow at nearly linear rate with matrix size.Its rank equals Π(n), while every completion has rank at least π(n), and both scale asymptotically as n/log n.
  • Sparse datasets: Three empirical matrices are large, sparse, and nonnegative, with fewer than 1% positive entries; all have very slowly decaying singular-value spectra.The datasets include MNIST similarities, a fruit-fly connectome with over 139K brain cells, and Wikipedia bigram statistics.

3 Algorithm

The paper formulates subzero matrix completion as approximating a sparse nonnegative matrix by the elementwise rectification of a low-rank matrix, then develops alternating optimization strategies for scalable computation.

  • Problem formulation: Subzero completion seeks a rank-r matrix L whose rectification max(0, L) approximates sparse nonnegative S.The prescribed rank r is typically much smaller than the matrix dimensions.
  • Objective function: RMSE and weighted Jaccard distance measure reconstruction error, but both ignore how strongly negative entries in L correspond to positive entries in S.Their derivatives also vanish for negative elements, limiting their suitability for hill-climbing methods.
  • Objective function: An auxiliary matrix Z yields an objective that is zero only when Z=L and S=max(0,L), while penalizing strongly negative predictions more severely.The objective also provides an upper bound on RMSE at feasible points.
  • Alternating minimization: Alternating minimization updates Z and the rank-r approximation L in turn, decreasing the objective monotonically.For fixed Z, L is obtained by a truncated SVD retaining the leading r components; for fixed L, Z is updated from the rectified values of L.
  • Factorized optimization: Factorizing L as AB⊤ replaces the full low-rank matrix with tall, skinny factors A and B and enables three-way updates over Z, A, and B.The updates can incorporate momentum, regularization, and extrapolation, and can be rewritten using the difference Δ=Z−AB⊤.

4 Scaling to large problems

The paper scales subzero matrix completion beyond dense-memory limits by combining blockwise and stochastic least-squares updates with GPU kernel optimizations. These strategies trade monotonic objective decrease against faster convergence and can be complemented by sparse computations, although their interaction with GPU acceleration is problem-dependent.

  • Memory management: Dense matrices Z or ∆ can exceed available memory, motivating updates that operate on smaller matrix blocks.The bottleneck arises when the dense intermediate cannot be stored all at once.
  • Parallel least-squares: Blockwise least-squares updates subsample rows of A and corresponding blocks of ∆, while analogous B updates use matching columns of S, Ŝ, and ∆.These updates can be parallelized across cores or use the largest blocks that fit in memory on a single core.
  • Stochastic least-squares: Stochastic least-squares conserves memory by introducing randomness into alternating minimization and using mini-batches constrained by the available storage for ∆.Momentum can further accelerate the stochastic updates.
  • Complementarity: Parallel updates monotonically decrease the objective without step sizes, whereas stochastic updates converge faster when different submatrices of S have similar statistics.For very large problems, the two strategies can be combined by subsampling both rows and columns of ∆.
  • Kernel fusion: Fused CUDA kernels avoid writing and rereading the dense intermediate AB^⊤, addressing the memory-bandwidth bottleneck of separate GPU launches.The fused operation combines multiplication, rectification, and the subsequent multiplication by the pseudoinverse.
  • Computational throughput: Custom fused TF32 kernels reach up to 75% of peak throughput on H100 and 82% on RTX Pro 6000 Blackwell GPUs.The comparison is against a two-kernel PyTorch baseline using epilogue fusion; Table 2 reports achieved throughputs in TFLOP/s.
  • Arithmetic intensity: Over 100x higher arithmetic intensity is achieved by the fused kernel for a representative rank-10 problem with m=n=10^5 and block size n_block=10^3.The paper states that fused kernels maintain much higher arithmetic intensity when m,n≫r, including on smaller blocks.
  • Exploiting sparsity: Sparse storage and range-search optimizations may contend with dense CUDA optimizations, so their combined benefit depends empirically on problem size, sparsity, hardware, and optimization stage.The paper does not claim that sparse and GPU-based speedups always reinforce one another.

5 Experimental results

The experiments show that SUMAC combines stochastic least-squares updates with custom GPU kernels to scale subzero matrix completion, while retaining competitive accuracy and revealing latent structure. On the fly connectome, SALSA is faster initially and in wall-clock time, whereas ADAM performs slightly better late in optimization.

  • Optimization design: Random-block stochastic least-squares updates combined with custom CUDA kernels formed the most effective SUMAC configuration for the tested problem sizes.The stochastic updates conserved memory, while GPU speedups were easier to realize than fast similarity search.
  • Optimization design: On a single H100 GPU, 1000 iterations took 9 minutes with PyTorch ADAM versus 2–3 minutes with SUMAC in Python and MATLAB.The comparison used the fly connectome experiments shown in Figure 5.
  • Optimization design: SALSA decreased the objective more quickly at the outset, while ADAM performed similarly or marginally better in the tail.Both optimizers used the same initialization, mini-batch size, and momentum; the objective was averaged over five initializations for rank 16.
  • Matrix approximation: SUMAC revealed latent low-rank structure that purely linear models could not, with crossover values below 10 for digits and connectome matrices and below 100 for bigram counts.These crossover values were also commensurate with the rank of the largest positive submatrix identified in Table 1.
  • Matrix approximation: In some severely underparameterized models, such as normalized bigram counts with rank r ≤32, SUMAC errors exceeded those from truncated SVD.The authors associate this behavior with shallow local minima and instability of leading subspaces across mini-batches.
  • Cell classification: Truncated SVD embeddings produced error rates far below random guessing, indicating that cell labels correlate with synaptic connectivity patterns.The linear low-rank structure also summarized these patterns more compactly than the full connectome matrix.
  • Cell classification: Connectome subzero completions achieved comparable or lower 1nn error rates with ranks 8–16 than truncated SVD embeddings evaluated up to rank 1024.The subzero models therefore used much lower-dimensional embeddings for cell-category classification.
  • Cell classification: The most accurate subzero classification results came from models optimized for only five epochs, despite crude RMSE, weighted Jaccard, and objective values.The authors describe this result as positive but not fully understood.

6 Related work

The paper builds on matrix factorization, sparse-completion, rectified low-rank, manifold-learning, and sign-rank perspectives. These connections situate subzero matrix completion within broader approaches to representing sparse or incomplete data with constrained low-dimensional structure.

  • Matrix factorization and completion: Low-rank factorization approximates a large real-valued matrix with two smaller factors, while truncated SVD applies when all elements are observed.Low-rank completion addresses settings with missing elements.
  • Matrix factorization and completion: Matrix factorization has been generalized with binary, discrete, bounded, nonnegative, sparse, stochastic, and low-norm constraints on factors or encoded matrices.These variants broaden the types of structure that factorization methods can represent.
  • Subzero matrix completion: The paper extends earlier rectified low-rank encodings of sparse nonnegative matrices and prior alternating-minimization methods, including expectation-maximization, coordinate descent, and ADMM.Subsequent work targeted the objective with increasingly varied optimization procedures.
  • Related geometric connections: Manifold-learning methods also connect sparse representations of nearest neighbors or local reconstructions with low-rank embeddings.This provides a related example of sparse and low-rank structure serving complementary roles.
  • Related geometric connections: For sparse Boolean matrices, sign rank is the minimum rank of a real matrix whose Heaviside threshold reproduces the Boolean pattern.The identity matrix of size n≥3 has sign rank 3.

7 Discussion

The paper characterizes when sparse nonnegative matrices admit lower-rank real-valued representations after negative entries are zeroed, and presents a scalable discovery algorithm for very large matrices.

  • The study analyzes when sparse nonnegative matrices can be recovered from much lower-rank real-valued matrices by zeroing negative entries.
  • The proposed stochastic updates operate on random blocks of factors A and B, enabling scaling to matrices with hundreds of thousands of rows and columns.
  • Further speedups may come from adaptive momentum and step-size selection, selective sparse optimizations, and randomized SVD methods.

A Details on sparse datasets

The appendix constructs sparse matrices from image-neighbor, connectome, and text-bigram data, illustrating diverse sparse nonnegative datasets used for empirical analysis.

  • Digits: The MNIST matrix contains 70K rows and columns, with nonzeros marking each image’s 16 nearest neighbors under cosine distance.
  • Digits: The MNIST construction is over 99.9% sparse, has exactly 16 nonzeros per row, and is asymmetric.
  • Connectome: The connectome matrix records synapse counts among 139255 cells, with labels available at superclass, class, subclass, and cell-type levels.
  • Bigrams: The bigram matrix is a stochastic matrix derived from counts for the 250K most frequent Wikipedia tokens, with nonzero rows normalized to sum to one.

B Experimental settings

The experiments generally use stochastic optimization with fixed default hyperparameters, while convergence requires adjustments for heavily underparameterized models and the largest bigram matrix.

  • Experimental settings: Most models were trained for 1000 epochs with step size η = 1, momentum γ = 0.9, and 100 mini-batches per epoch.
  • Experimental settings: Table 7 lists the top 10 labels for each connectome cell category, while Table 8 lists model hyperparameter settings.
  • Experimental settings: Smaller step sizes and larger mini-batches are required for the smallest, most heavily underparameterized models to converge.
  • Experimental settings: 4000 rather than 1000 iterations are generally required for convergence on the largest sparse bigram matrix, with nearly 250K rows and columns.

C GPU optimizations

The GPU implementation accelerates fused matrix computations through Tensor Core kernels, tiled intermediate processing, asynchronous copy/compute overlap, autotuning, and JIT compilation.

  • GPU optimizations: The fused kernel avoids the split-kernel computation’s 2wmn-byte memory round-trip for the matrix product AB⊤.
  • GPU optimizations: Two implementations use warp-synchronous mma instructions for A100 GPUs and warp-group asynchronous wgmma instructions for H100 GPUs.
  • GPU optimizations: Tiling computes local intermediates Lij = AiBj⊤ and accumulates Yi += max(0, Lij)Cj without writing the intermediate tile to global memory.
  • GPU optimizations: Two staging buffers support asynchronous copy/compute overlap, with producers reusing a completed buffer for tile j+2.
  • GPU optimizations: Autotuning selects tile sizes, warp or warpgroup row allocations, pipeline stages, operand modes, and instruction shapes for each input shape.
  • GPU optimizations: JIT compilation avoids precompiling many parameter combinations, and plain FP32 fused kernels are provided alongside TF32 kernels.
Loading 2608.21607v1…