9
$\begingroup$

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?

Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Fascinating. How'd you find "getU"? I hadn't realized a LinearSolveFunction had "Properties". $\endgroup$
    – b3m2a1
    Jun 15, 2018 at 23:03
  • $\begingroup$ @b3m2a1 Erm. If I am honest I cannot recall. Maybe from this post or this post by J.M.? $\endgroup$
    – Henrik Schumacher
    Jun 15, 2018 at 23:14
  • 1
    $\begingroup$ Ah of course :) These interesting undocumented things always seem to come back to J.M. in the end $\endgroup$
    – b3m2a1
    Jun 15, 2018 at 23:16
  • 1
    $\begingroup$ @HenrikSchumacher The rank of an upper triangular matrix is not necessarily equal to the amount of non-zero diagonal elements. { {0,1}, {0,0} } has rank 1 yet it has only 0 elements on the diagonal. $\endgroup$
    – Gert
    Dec 9, 2019 at 14:41
  • $\begingroup$ @Gert Good point! I wrongly assumed that the pivoting would move ones onto the diagonal of the upper right block, but Mathematica actually complains in this case because it finds out that the matrix is singular =( Probably QRDecomposition (or SPQR for sparse matrices) are better suited for this task. $\endgroup$
    – Henrik Schumacher
    Dec 9, 2019 at 14:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.