How to compute the rank of a large sparse matrix in MATLAB

Question

I am interested in computing the ranks of fairly large, the largest being of magnitude $10^6$ x $10^6$, sparse matrices whose entires are all 0, 1, or -1. I have been trying to use Matlab to accomplish this. In particular, my approach has been to use QR-factorization -- that is factoring my matrix $T$ as $TP=QR$ were $P$ is a permutation matrix, $Q$ is orthogonal, and $R$ is upper triangular -- and then counting the number of non-zero entries on the diagonal of $R$.

Using Matlab's built in QR function the following Matlab code computes the correct rank for some of my examples. (I have confirmed the rank for some of the smaller matrices against rank(full(T)).)

load(‘FILE_NAME.dat') 
T = spconvert(FILE_NAME);

tic

[Q,R,P]=qr(T);
[m,n] = size(T);
tol = max(size(T))*eps*abs(R(1,1));
r = 0;

while (r<min(m,n) & abs(R(r+1,r+1)) >= tol)
r = r+1;
end

tim=toc
r

That said even for some of my smaller matrices this code seems fairly slow and memory intensive. For example, on one of these matrices that is ~70,000 x ~50,000 the above code takes numerous hours and most of the memory on a 64 GB server. (One thing I noticed is that if I save the Q, R, and P computed above to a .mat file the file size seems to grow quickly. For example, the .mat file is ~1 GB when initial matrix is ~15,000 x ~15,000)

Is there a better (quicker or less memory intensive) way to approach computing the rank of large sparse matrices? For example, a more effective way to implement QR-factorization or even a method entirely different from QR?

Answer 1

There are two things you can do that may significantly reduce the computation time and memory used.

First, you aren't using the Q matrix, so don't ask MATLAB to compute it. Since it is dense, that will definitely save a lot of memory. You can do simply:

R=qr(T);

Second, in sparse factorizations, the ordering of the equations can have a dramatic effect on the amount of fill-in and computation time during factorization. This is true even for QR and LU factorizations where pivoting takes place to insure numerical stability.

An effective strategy in sparse QR factorizations is to first reorder the columns to preserve sparsity. Sometimes that will make a huge difference in both memory usage and factorization time-- depending on the matrix. That can be done, as follows:

[m,n]=size(T);
q=colamd(T, [n m]);
R=qr(T(:,q));

The details of the rank-revealing properties of this algorithm (and implementation) are discussed in this paper by Davis:

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.319.4864&rep=rep1&type=pdf

Answer 2

A rank revealing LU factorization such as LUSOL might be what you want:

http://web.stanford.edu/group/SOL/software/lusol/

Answer 3

If you know the rank will be small (say less than 100) and your matrix is square, use eigs(A,100) to get the 100 largest eigenvalues in magnitude. The eigenvalues are solved using sparse matrix techniques and would be much faster. In such a case, there is no need to create a full matrix.

On the other hand, if you know the rank will be close to the matrix size, then maybe try eigs(A,100,0) which should give you the 100 smallest eigenvalues.