ENH - added matrix_rank function to linalg

2 files changed, 84 insertions, 3 deletions

diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py
index ce319c17a..e0dbb9e90 100644
--- a/numpy/linalg/linalg.py
+++ b/numpy/linalg/linalg.py
@@ -11,13 +11,14 @@ zgetrf, dpotrf, zpotrf, dgeqrf, zgeqrf, zungqr, dorgqr.
 
 __all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv', 'inv',
            'cholesky', 'eigvals', 'eigvalsh', 'pinv', 'det', 'svd',
-           'eig', 'eigh','lstsq', 'norm', 'qr', 'cond', 'LinAlgError']
+           'eig', 'eigh','lstsq', 'norm', 'qr', 'cond', 'matrix_rank',
+           'LinAlgError']
 
 from numpy.core import array, asarray, zeros, empty, transpose, \
         intc, single, double, csingle, cdouble, inexact, complexfloating, \
         newaxis, ravel, all, Inf, dot, add, multiply, identity, sqrt, \
         maximum, flatnonzero, diagonal, arange, fastCopyAndTranspose, sum, \
-        isfinite, size
+        isfinite, size, finfo
 from numpy.lib import triu
 from numpy.linalg import lapack_lite
 from numpy.matrixlib.defmatrix import matrix_power
@@ -1376,6 +1377,66 @@ def cond(x, p=None):
     else:
         return norm(x,p)*norm(inv(x),p)
 
+
+def matrix_rank(M, tol=None):
+    ''' Return matrix rank of array using SVD method
+
+    Rank of the array is the number of SVD singular values of the
+    array that are greater than `tol`.
+    
+    Parameters
+    ----------
+    M : array-like
+        array of <=2 dimensions
+    tol : {None, float}
+       threshold below which SVD values are considered zero. If `tol` is
+       None, and ``S`` is an array with singular values for `M`, and
+       ``eps`` is the epsilon value for datatype of ``S``, then `tol`
+       set to ``S.max() * eps``.
+
+    Examples
+    --------
+    >>> matrix_rank(np.eye(4)) # Full rank matrix
+    4
+    >>> I=np.eye(4); I[-1,-1] = 0. # rank deficient matrix
+    >>> matrix_rank(I)
+    3
+    >>> matrix_rank(np.ones((4,))) # 1 dimension - rank 1 unless all 0
+    1
+    >>> matrix_rank(np.zeros((4,)))
+    0
+
+    Notes
+    -----
+    Golub and van Loan [1]_ define "numerical rank deficiency" as using
+    tol=eps*S[0] (where S[0] is the maximum singular value and thus the
+    2-norm of the matrix). This is one definition of rank deficiency,
+    and the one we use here.  When floating point roundoff is the main
+    concern, then "numerical rank deficiency" is a reasonable choice. In
+    some cases you may prefer other definitions. The most useful measure
+    of the tolerance depends on the operations you intend to use on your
+    matrix. For example, if your data come from uncertain measurements
+    with uncertainties greater than floating point epsilon, choosing a
+    tolerance near that uncertainty may be preferable.  The tolerance
+    may be absolute if the uncertainties are absolute rather than
+    relative.
+
+    References
+    ----------
+    .. [1] G. H. Golub and C. F. Van Loan, _Matrix Computations_.
+    Baltimore: Johns Hopkins University Press, 1996.
+    '''
+    M = asarray(M)
+    if M.ndim > 2:
+        raise TypeError('array should have 2 or fewer dimensions')
+    if M.ndim < 2:
+        return int(not all(M==0))
+    S = svd(M, compute_uv=False)
+    if tol is None:
+        tol = S.max() * finfo(S.dtype).eps
+    return sum(S > tol)
+
+
 # Generalized inverse
 
 def pinv(a, rcond=1e-15 ):
diff --git a/numpy/linalg/tests/test_linalg.py b/numpy/linalg/tests/test_linalg.py
index 01e96d6a2..b81561254 100644
--- a/numpy/linalg/tests/test_linalg.py
+++ b/numpy/linalg/tests/test_linalg.py
@@ -6,7 +6,7 @@ from numpy.testing import *
 from numpy import array, single, double, csingle, cdouble, dot, identity
 from numpy import multiply, atleast_2d, inf, asarray, matrix
 from numpy import linalg
-from numpy.linalg import matrix_power, norm
+from numpy.linalg import matrix_power, norm, matrix_rank
 
 def ifthen(a, b):
     return not a or b
@@ -315,5 +315,25 @@ class TestNormSingle(_TestNorm):
     dt = np.float32
     dec = 6
 
+
+def test_matrix_rank():
+    # Full rank matrix
+    yield assert_equal, 4, matrix_rank(np.eye(4))
+    # rank deficient matrix
+    I=np.eye(4); I[-1,-1] = 0. 
+    yield assert_equal, matrix_rank(I), 3
+    # All zeros - zero rank
+    yield assert_equal, matrix_rank(np.zeros((4,4))), 0
+    # 1 dimension - rank 1 unless all 0
+    yield assert_equal, matrix_rank([1, 0, 0, 0]), 1
+    yield assert_equal, matrix_rank(np.zeros((4,))), 0
+    # accepts array-like
+    yield assert_equal, matrix_rank([1]), 1
+    # greater than 2 dimensions raises error
+    yield assert_raises, TypeError, matrix_rank, np.zeros((2,2,2))
+    # works on scalar
+    yield assert_equal, matrix_rank(1), 1
+    
+    
 if __name__ == "__main__":
     run_module_suite()