Calculating the rank of a huge sparse array

By virtue of the suggestion in my previous question, I constructed the sparse matrix whose size is $518400 \times 86400$, mostly filled with $0$ and $\pm 1$. Now I want to calculate its rank.

Since RowReduce requires a huge amount of memory to fill in $0$s, my first approach was to divide the matrix into many parts, e.g., each containing 144 rows ($144\times 86400$).

Module[{i, part, rank},
For [i = 0, 144 (i + 1) <= 518400, i++,
part = RowReduce[input[[144*i + 1 ;; 144 (i + 1)]]];
rank = MatrixRank[part];
If[i == 0, result = part[[1 ;; rank]],
result = Join[result, part[[1 ;; rank]]]]
];
]


Unfortunately, this does not work because all rows in each part are linearly independent so that the number of rows does not decrease.

Does anyone have such an experience? Should I export the matrix to the file and calculate its rank via other language / packages? If so, any suggestion?

I am not sure how robust this method is, but for a sparse, square matrix A, LinearSolve[A, Method-> "Multifrontal"] computes an LU-factorization (by using UMFPACK as the backend). As far as I can tell, the L-factor is always made invertible, so that the rank of A equals the rank of the U-factor. Since the U-factor is an upper-triangular matrix, its rank equals the number of nonzero entries on its diagonal.

The following function exploits this by filling up the input matrix with zeroes until it becomes a square matrix, calculates the LU-decomposition, and counts the nonzero diagonal elements of the U-factor.

SparseMatrixRank[A_SparseArray?MatrixQ] := Module[{a, d1, d2},
{d1, d2} = Dimensions[A];
a = If[d1 > d2,
Join[Transpose[A],
SparseArray[{}, {d1 - d2, d1}, A["Background"]]],
Join[A, SparseArray[{}, {d2 - d1, d2}, A["Background"]]]
];
Total[Unitize[Diagonal[Quiet@LinearSolve[a, Method -> "Multifrontal"]["getU"]]]]
];


Here a test suite for small-sized matrices

n = 1000;
d2 = 1200;
d1 = 500;
i = RandomInteger[{1, d1}, n];
j = RandomInteger[{1, d2}, n];
v = RandomReal[{-1, 1}, n];
A = SparseArray[Transpose[{i, j}] -> v, {d1, d2}];
SparseMatrixRank[A]
MatrixRank[A]


I haven't found any discrepancies yet. However, UMFPACK has its problems with random matrices with several million nonzero values. Actually, the topology of the graph behind a sparse matrix is much more important than its size or its number of nonzero values if it comes to factorizations. So it is hard to predict if it will work for your particular matrices.

• Fascinating. How'd you find "getU"? I hadn't realized a LinearSolveFunction had "Properties". – b3m2a1 Jun 15 '18 at 23:03
• @b3m2a1 Erm. If I am honest I cannot recall. Maybe from this post or this post by J.M.? – Henrik Schumacher Jun 15 '18 at 23:14
• Ah of course :) These interesting undocumented things always seem to come back to J.M. in the end – b3m2a1 Jun 15 '18 at 23:16