[Numpy-discussion] Matlab page on scipy wiki

[Fernando Perez]

Bill Baxter wrote:
> For what it's worth, matlab's rank function just calls svd, and returns the
> number singular values greater than a tolerance.  The implementation is a
> whopping 5 lines long.

Yup, and it would be pretty much the same 5 lines in numpy, with the same 
semantics.

Here's a quick and dirty implementation for old-scipy (I don't have new-scipy 
on this box):

def matrix_rank(arr,tol=1e-8):
     """Return the matrix rank of an input array."""

     arr = scipy.asarray(arr)

     if len(arr.shape) != 2:
         raise ValueError('Input must be a 2-d array or Matrix object')

     svdvals = scipy.linalg.svdvals(arr)
     return sum(scipy.where(svdvals>tol,1,0))


If you really hate readability and error-checking, it's a one-liner :)

matrix_rank = lambda arr,tol=1e-8: 
sum(scipy.where(scipy.linalg.svdvals(arr)>tol,1,0))


Looks OK (RA is RandomArray from Numeric):

In [21]: matrix_rank([[1,0],[0,0]])
Out[21]: 1

In [22]: matrix_rank(RA.random((3,3)))
Out[22]: 3

In [23]: matrix_rank([[1,0],[0,0]])
Out[23]: 1

In [24]: matrix_rank([[1,0],[1,0]])
Out[24]: 1

In [25]: matrix_rank([[1,0],[0,1]])
Out[25]: 2

In [26]: matrix_rank(RA.random((3,3)),1e-1)
Out[26]: 2

In [48]: matrix_rank([[1,0],[1,1e-8]])
Out[48]: 1

In [49]: matrix_rank([[1,0],[1,1e-4]])
Out[49]: 2

In [50]: matrix_rank([[1,0],[1,1e-8]],1e-9)
Out[50]: 2


Cheers,

f