Merge pull request #9582 from seberg/rcond-default

ENH: Warn to change lstsq default for rcond

2 files changed, 40 insertions, 2 deletions

diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py
index dd5ac6255..edd980dd5 100644
--- a/numpy/linalg/linalg.py
+++ b/numpy/linalg/linalg.py
@@ -1810,7 +1810,7 @@ def det(a):
 
 # Linear Least Squares
 
-def lstsq(a, b, rcond=-1):
+def lstsq(a, b, rcond="warn"):
     """
     Return the least-squares solution to a linear matrix equation.
 
@@ -1836,6 +1836,13 @@ def lstsq(a, b, rcond=-1):
         as zero if they are smaller than `rcond` times the largest singular
         value of `a`.
 
+        .. versionchanged:: 1.14.0
+           If not set, a FutureWarning is given. The previous default
+           of ``-1`` will use the machine precision as `rcond` parameter,
+           the new default will use the machine precision times `max(M, N)`.
+           To silence the warning and use the new default, use ``rcond=None``,
+           to keep using the old behavior, use ``rcond=-1``.
+
     Returns
     -------
     x : {(N,), (N, K)} ndarray
@@ -1909,6 +1916,20 @@ def lstsq(a, b, rcond=-1):
     if m != b.shape[0]:
         raise LinAlgError('Incompatible dimensions')
     t, result_t = _commonType(a, b)
+    # Determine default rcond value
+    if rcond == "warn":
+        # 2017-08-19, 1.14.0
+        warnings.warn("`rcond` parameter will change to the default of "
+                      "machine precision times ``max(M, N)`` where M and N "
+                      "are the input matrix dimensions.\n"
+                      "To use the future default and silence this warning "
+                      "we advise to pass `rcond=None`, to keep using the old, "
+                      "explicitly pass `rcond=-1`.",
+                      FutureWarning, stacklevel=2)
+        rcond = -1
+    if rcond is None:
+        rcond = finfo(t).eps * ldb
+
     result_real_t = _realType(result_t)
     real_t = _linalgRealType(t)
     bstar = zeros((ldb, n_rhs), t)
diff --git a/numpy/linalg/tests/test_linalg.py b/numpy/linalg/tests/test_linalg.py
index 97b2e328a..ab81fc485 100644
--- a/numpy/linalg/tests/test_linalg.py
+++ b/numpy/linalg/tests/test_linalg.py
@@ -793,7 +793,7 @@ class TestLstsq(LinalgSquareTestCase, LinalgNonsquareTestCase):
         arr = np.asarray(a)
         m, n = arr.shape
         u, s, vt = linalg.svd(a, 0)
-        x, residuals, rank, sv = linalg.lstsq(a, b)
+        x, residuals, rank, sv = linalg.lstsq(a, b, rcond=-1)
         if m <= n:
             assert_almost_equal(b, dot(a, x))
             assert_equal(rank, m)
@@ -814,6 +814,23 @@ class TestLstsq(LinalgSquareTestCase, LinalgNonsquareTestCase):
         assert_(imply(isinstance(b, matrix), isinstance(x, matrix)))
         assert_(imply(isinstance(b, matrix), isinstance(residuals, matrix)))
 
+    def test_future_rcond(self):
+        a = np.array([[0., 1.,  0.,  1.,  2.,  0.],
+                      [0., 2.,  0.,  0.,  1.,  0.],
+                      [1., 0.,  1.,  0.,  0.,  4.],
+                      [0., 0.,  0.,  2.,  3.,  0.]]).T
+
+        b = np.array([1, 0, 0, 0, 0, 0])
+        with suppress_warnings() as sup:
+            w = sup.record(FutureWarning, "`rcond` parameter will change")
+            x, residuals, rank, s = linalg.lstsq(a, b)
+            assert_(rank == 4)
+            x, residuals, rank, s = linalg.lstsq(a, b, rcond=-1)
+            assert_(rank == 4)
+            x, residuals, rank, s = linalg.lstsq(a, b, rcond=None)
+            assert_(rank == 3)
+            # Warning should be raised exactly once (first command)
+            assert_(len(w) == 1)
 
 class TestMatrixPower(object):
     R90 = array([[0, 1], [-1, 0]])