ENH: A very fast algorithm for computing matrix rank

[touqir14]

🚀 Feature

I have devised an LUP based matrix rank computing algorithm that outperforms significantly the standard Numpy's matrix_rank function and PyTorch's GPU accelerated matrix_rank algorithm for most cases. For instance, for UAR generated 15000 x 15000 matrices, it offers a speedup of about 25 times over Numpy's matrix_rank. See the benchmarks. Unlike the LU factorization, LUP is stable in practice and I have thoroughly tested the correctness of my implementation.

Motivation

Computing the matrix rank using Numpy's matrix_rank function can be rather slow, particularly if one deals with large matrices and needs to run multiple times for at least tens of times. This motivated me to come up with a much faster approach than the SVD method used in Numpy's matrix_rank.

Pitch

I would like to consider the prospect of implementing this as a Scipy function. My implementation uses scipy's LUP function and a couple of Cython optimized functions along with the standard Numpy matrix operations.

Let me know if you need to know more about the algorithm's specifics and testing/benchmarks.

[ilayn]

Hi @touqir14 that's great that you are working on this. The immediate next step is actually to validate the results based on precision and not so much about the speed itself.

The main issue with matix-rank and for all other upper/lower-semicontinuous functions for that matter is that they have to be correct at exactly very precise points. Hence while reaping the speed benefits of LU as opposed to an SVD decomposition it should not loose precision which is I am a bit skeptical about since their comparison has been pretty well-established in terms of precision elsewhere.

[rgommers]

Also, is there any literature on this method?

[touqir14]

Hi @ilayn , I agree with your point. I am not trying to propose a method that matches the precision of SVD or even beats it. Having said that, I have carried various tests on this algorithm and the one case where this method has some likelihood of failing is when the spectrum of the singular values of the input matrix is quite large. For instance, I have been randomly generating rank-deficient matrices which have the ratio of the largest singular value to the smallest non-zero singular value to be at least as large as 10^11 and thats pretty much the breaking point. Below 10^11, at 10^10, there are cases of failure that are rather rare. For a matrix of shape 100 x 100, after carrying out 1000 runs I could notice a couple of failures (SVD rank and LUP rank do not match up). Also, the mismatches were usually by a gap of 1. For a larger matrix of size 1000 x 1000, there does not seem to be any failure from 1000 runs. Below 10^10, the results seemed to be fine. I believe this algorithm is appropriate for real world cases where it is unlikely to come across matrices with a very large singular value spectrum. The SVD matrix_rank method's breaking point is probably close to 10^13. Do let me know if you are looking for results on some particular tests that you think may add validity. I have also carried out tests on random matrices from various distributions and they all were fine.

Hi @rgommers , from my literature survey it seems like the literature on rank revealing LUP is quite small. Smaller than other rank revealing techniques like RRQR. Here is an algorithm https://core.ac.uk/download/pdf/82771231.pdf that uses full pivoting. My algorithm uses partial pivoting and from what I know, partial pivoting offers the same precision as full pivoting while being faster than the later. To the best of my knowledge, my algorithm is not present in the literature. The principle behind my algorithm, briefly speaking, is that the partial pivoted LU factorization breaks the input matrix A into PA = LU. The rank of PA is the same as rank of A. Since L is a unit diagonal matrix, it is a full rank matrix. This leaves us to U which reflects the rank deficiency of A. Thus the problem gets reduced to computing the rank of U. U is not necessarily in a Row Reduced Echelon (REF) form. The crux of the algorithm is to perform successive LUP factorizations on the sub-matrices of U in a particular way to convert it into its REF form and from there the rank can be easily computed.

[rgommers]

Ref the same proposal to PyTorch: pytorch/pytorch#41622

[touqir14]

The GPU accelerated algorithm uses PyTorch's LUP function and my benchmarks show that this algorithm significantly outperforms PyTorch's GPU accelerated matrix rank function.

[rgommers]

How about adding a graph on precision like percentage failures as a function of ratio of smallest to largest singular value, for np.linalg.matrix_rank and your LUP-based function?

[touqir14]

That would be good. I will post that.

[touqir14]

I have slightly modified the algorithm for better precision. I have written a post comparing the success rates of LUP-rank-computer and numpy.linalg.matrix_rank

[ilayn]

Can you please also add how you generate these matrices?

One thing I might warn about these benchmarks is that it is very hard to get rank-deficient matrices through random generators (other than the obvious A @ A.T for rectangular A). You are not necessarily generating problematic rank-deficient matrices but mostly ill-conditioned matrices.

The ill-conditioned matrices are kind of OK for SVD and not-so-much for LU (there is much literature about this) but matrices in the form of

A = np.array([[1., 1.],
              [1+eps, 1.]])

are the actual enemies as eps approaches 0. If you even use this matrix you would get an idea of precision differences.

Moreover the success/fail should be predictable and extendable. In other words, there must be a way to shrink/increase tolerance values as in np.linalg.matrix_rank. I don't think this possible for LU based approaches since the elimination is hidden from user.

[touqir14]

I have added a subsection in that post detailing how I created those random matrices.

I can see that the actual rank may not be the same as the intended rank but I think thats only the case for a very large target ratio of the largest singular value to the smallest non-zero singular value, like at least 1e13. For those cases I can definitely test by creating random rectangular matrices of condition number = t^0.5 whereas t is the target ratio of the final matrix since doing A @ A.T squares the condition number of A.

Your example highlights the weakness of LU factorization when used for solving systems of equations. The LUP matrix computer works fine upto eps=1e-12 while numpy.linalg.matrix_rank works fine upto eps=1e-14. This accuracy gap between SVD's rank algorithm and my LUP based rank computer has already been highlighted in the experiments since eps=1e-12 corresponds to a condition number almost 1e12 while eps=1e-14 corresponds to 1e14.

The tolerance value that I am using as of now is not chosen based on theoretical grounds as in numpy.linalg.matrix_rank and I am in search for one. Can the bound (page 22, theorem 9.3) in https://who.rocq.inria.fr/Laura.Grigori/TeachingDocs/UPMC_Master2/Slides_Spr2016/DenseLU.pdf be used to compute a tolerance?

[touqir14]

@ilayn , changing the random matrix generation process to A @ A.T does degrade the success rate of LUP method for those ranges of condition number where it previous had a 100% success rate such that for most condition numbers, at least a couple of cases out of a 1000 pop up where the LUP method fails. The success rate of SVD remains almost unchanged over all condition numbers. Slightly changing the tolerance level of the LUP method allows it to achieve a 100% success rate for those condition numbers where it achieved a 100% success rate using the previous random matrix generation method. I think the focus should go towards finding a principled way of tuning the tolerance level.

[ilayn]

I don't have too much hope to get equivalent performance from LUP to be honest. Even if there is a sensible tolerance machinery as long as it cannot match the SVD results it would still be inferior.

The performance should also worsen as the matrices gets taller instead of squares.

[touqir14]

Even if it is stays inferior to SVD, don't you think it still has utility as long as it works for most cases? I was not hoping to match the SVD results in the first place. The original motivation was to provide an algorithm that works quite fast compared to SVD based matrix_rank algorithm when the later (possibly over multiple runs) can take hours for large matrices. For most of the failed cases, the matrix rank was off by 1. There are algorithms that even require the use of approximate rank.

[ilayn]

Oh don't get me wrong I think this is still a valuable result but matrix rank is notoriously problematic and expensive hence often we try to avoid computing it and try to resort to other means to work around it. Many commercial software do the same.

To gain more speed and to approximate the rank, often people resort to iterative methods or randomized methods. For dense arrays and algebraic methods, a wrong result is not acceptable unfortunately. For the latter approaches that's perfectly fine. But then your opponents are https://arxiv.org/abs/1503.07157 and alike.

[touqir14]

Can you recommend me approximate rank computational algorithms that come with scipy? I know the interpolative decomposition comes with an approximate rank computational method in scipy and it might be a good candidate to compare with. Also, can the LUP rank algorithm be considered an addition to scipy as an approximate rank computational algorithm like the ID's approximation rank computational algorithm if it is a better competitor than the existing approximate rank computational algorithms in scipy?

[touqir14]

On second thoughts, ID's rank computational method will be a poor competitor. See #12593

[touqir14]

@ilayn , I have added experimental results comparing the LUP-rank-computer with the https://arxiv.org/abs/1503.07157 approximate rank computational algorithm here. The LUP-rank-computer seems to be superior. I don't think they have mentioned how the tolerance parameter can be set.

[ilayn]

I will check as soon as I can spare some time but another thing to note

One such case is when the input matrix has a very wide spectrum of singular values i.e when the ratio of the largest singular value to the smallest singular value is quite big.

That's true for every rank-deficient matrix. Because sigma_min is always 0. for a rank deficient matrix. That's actually what you are trying to figure out. Hence sigmamax/sigmamin comparison is not an indicator of anything. If sigmamax = 100 (which is pretty common) and sigmamin is 1e-15 (which is roughly the size for rank deficiency) you get an artificial 1e17 which is just noise.

[touqir14]

Actually, I meant the smallest non-zero singular value.

[ilayn]

That's a problem for every algorithm so LUP is no exception.