Sparse matrix rank: would a sparse-LU-based matrix_rank fit in scipy.sparse.linalg?

Contributor & Development Discussion
,
1

Disclosure: The code in the linked repository was partially written with the assistance of Claude, an AI assistant. I wrote this issue myself.

Problem

I was solving a Rank and Nullity problem (https://trentorch.com/ide/math-rank-and-nullity), and that is when I came to know that we can’t get the numerical rank of a sparse matrix without going through a dense process. The usual way, numpy.linalg.matrix_rank(A.toarray()), requires an N x N dense matrix, and its time complexity is about O(N^3) because of the SVD. For banded or low-fill matrices (such as finite-difference and finite-element discretizations), this usual way becomes impractical. It struggles well before N = 1e5 (100,000): at that size the matrix has shape 100,000 x 100,000 in float64 and needs 80 GB of memory.

Solution

I came across this problem and designed a solution: a function such as scipy.sparse.linalg.matrix_rank (open to suggestions) that works on square sparse matrices. It factorizes the matrix with splu. If the U factor (upper triangular matrix) is well-scaled and contains no nearly zero pivots, it reports full rank. Otherwise, the function falls back to dense SVD, so the result remains reliable.

I have created a prototype with tests and a benchmark script:

Running time of my prototype compared with dense matrix_rank, on a full-rank tridiagonal matrix on my machine (AMD64 Windows, Python 3.12.10, NumPy 2.5.4, SciPy 1.18.1):

  • N = 4,000: 1.3 ms vs 15.5 s
  • N = 10,000: 3.4 ms
  • N = 100,000: 38 ms (dense would require 80 GB)

Alternatives

  • Dense SVD via numpy.linalg.matrix_rank: correct, but it takes about O(N^3) time and O(N^2) memory.
  • Sparse QR factorization (e.g., SuiteSparseQR): it is rank-revealing and more robust, but SciPy does not have this functionality yet.
  • It is also possible to keep this as a separate package if its scope does not align with SciPy.

Limitations:

  • The high speed only applies to full-rank square sparse matrices. If the matrix is rank-deficient or not square, the function falls back to dense SVD and gives no performance benefit (1.31 s vs 1.25 s at N = 2000 in my case).
  • LU decomposition with partial pivoting is not rank-revealing in general. Therefore, the prototype uses LU factorization only to verify the “full-rank” condition when the pivots are well scaled. Otherwise, it falls back to dense SVD.
  • O(N) complexity applies to banded and low-fill sparse matrices. LU factorization of a general sparse matrix can result in a substantial amount of fill-in.

Questions for the maintainers:

  1. Is the full-rank test for sparse matrices too specialized to be useful in scipy.sparse.linalg, or is it just right?
  2. Would a rank-revealing algorithm (sparse QR) be preferred, and is there interest in adding a SuiteSparse dependency?
  3. Are there any test matrices it should be tested against?