ENH: add hermitian=False kwarg to matrix_power

With a symmetric matrix, the more efficient `eigvalsh` method
can be used to find singular values.

2 files changed, 29 insertions, 7 deletions

diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py
index 4005fc93d..d2ae7befc 100644
--- a/numpy/linalg/linalg.py
+++ b/numpy/linalg/linalg.py
@@ -24,8 +24,8 @@ from numpy.core import (
     add, multiply, sqrt, maximum, fastCopyAndTranspose, sum, isfinite, size,
     finfo, errstate, geterrobj, longdouble, moveaxis, amin, amax, product, abs,
     broadcast, atleast_2d, intp, asanyarray, isscalar, object_, ones, matmul,
-    swapaxes, divide)
-
+    swapaxes, divide, count_nonzero
+)
 from numpy.core.multiarray import normalize_axis_index
 from numpy.lib import triu, asfarray
 from numpy.linalg import lapack_lite, _umath_linalg
@@ -1504,11 +1504,11 @@ def cond(x, p=None):
         return norm(x, p, axis=(-2, -1)) * norm(inv(x), p, axis=(-2, -1))
 
 
-def matrix_rank(M, tol=None):
+def matrix_rank(M, tol=None, hermitian=False):
     """
     Return matrix rank of array using SVD method
 
-    Rank of the array is the number of SVD singular values of the array that are
+    Rank of the array is the number of singular values of the array that are
     greater than `tol`.
 
     .. versionchanged:: 1.14
@@ -1526,6 +1526,12 @@ def matrix_rank(M, tol=None):
 
         .. versionchanged:: 1.14
            Broadcasted against the stack of matrices
+    hermitian : bool, optional
+        If True, `M` is assumed to be Hermitian (symmetric if real-valued),
+        enabling a more efficient method for finding singular values.
+        Defaults to False.
+
+        .. versionadded:: 1.14
 
     Notes
     -----
@@ -1589,12 +1595,15 @@ def matrix_rank(M, tol=None):
     M = asarray(M)
     if M.ndim < 2:
         return int(not all(M==0))
-    S = svd(M, compute_uv=False)
+    if hermitian:
+        S = abs(eigvalsh(M))
+    else:
+        S = svd(M, compute_uv=False)
     if tol is None:
         tol = S.max(axis=-1, keepdims=True) * max(M.shape[-2:]) * finfo(S.dtype).eps
     else:
-        tol = asarray(tol)[...,newaxis]
-    return (S > tol).sum(axis=-1)
+        tol = asarray(tol)[..., newaxis]
+    return count_nonzero(S > tol, axis=-1)
 
 
 # Generalized inverse
diff --git a/numpy/linalg/tests/test_linalg.py b/numpy/linalg/tests/test_linalg.py
index fa20cc5ea..8b3984883 100644
--- a/numpy/linalg/tests/test_linalg.py
+++ b/numpy/linalg/tests/test_linalg.py
@@ -1383,6 +1383,19 @@ class TestMatrixRank(object):
         # works on scalar
         yield assert_equal, matrix_rank(1), 1
 
+    def test_symmetric_rank(self):
+        yield assert_equal, 4, matrix_rank(np.eye(4), hermitian=True)
+        yield assert_equal, 1, matrix_rank(np.ones((4, 4)), hermitian=True)
+        yield assert_equal, 0, matrix_rank(np.zeros((4, 4)), hermitian=True)
+        # rank deficient matrix
+        I = np.eye(4)
+        I[-1, -1] = 0.
+        yield assert_equal, 3, matrix_rank(I, hermitian=True)
+        # manually supplied tolerance
+        I[-1, -1] = 1e-8
+        yield assert_equal, 4, matrix_rank(I, hermitian=True, tol=0.99e-8)
+        yield assert_equal, 3, matrix_rank(I, hermitian=True, tol=1.01e-8)
+
 
 def test_reduced_rank():
     # Test matrices with reduced rank