Merge pull request #2965 from charris/fix-qr-mode

Fix qr mode

3 files changed, 213 insertions, 34 deletions

diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py
index c66ab8c3a..ba242e7c6 100644
--- a/numpy/linalg/linalg.py
+++ b/numpy/linalg/linalg.py
@@ -16,6 +16,8 @@ __all__ = ['matrix_power', 'solve', 'tensorsolve', 'tensorinv', 'inv',
            'svd', 'eig', 'eigh','lstsq', 'norm', 'qr', 'cond', 'matrix_rank',
            'LinAlgError']
 
+import warnings
+
 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, \
@@ -539,7 +541,7 @@ def cholesky(a):
 
 # QR decompostion
 
-def qr(a, mode='full'):
+def qr(a, mode='reduced'):
     """
     Compute the qr factorization of a matrix.
 
@@ -548,24 +550,42 @@ def qr(a, mode='full'):
 
     Parameters
     ----------
-    a : array_like
-        Matrix to be factored, of shape (M, N).
-    mode : {'full', 'r', 'economic'}, optional
-        Specifies the values to be returned. 'full' is the default.
-        Economic mode is slightly faster then 'r' mode if only `r` is needed.
+    a : array_like, shape (M, N)
+        Matrix to be factored.
+    mode : {'reduced', 'complete', 'r', 'raw', 'full', 'economic'}, optional
+        If K = min(M, N), then
+
+        'reduced'  : returns q, r with dimensions (M, K), (K, N) (default)
+        'complete' : returns q, r with dimensions (M, M), (M, N)
+        'r'        : returns r only with dimensions (K, N)
+        'raw'      : returns h, tau with dimensions (N, M), (K,)
+        'full'     : alias of 'reduced', deprecated
+        'economic' : returns h from 'raw', deprecated.
+
+        The options 'reduced', 'complete, and 'raw' are new in numpy 1.8,
+        see the notes for more information. The default is 'reduced' and to
+        maintain backward compatibility with earlier versions of numpy both
+        it and the old default 'full' can be omitted. Note that array h
+        returned in 'raw' mode is transposed for calling Fortran. The
+        'economic' mode is deprecated.  The modes 'full' and 'economic' may
+        be passed using only the first letter for backwards compatibility,
+        but all others must be spelled out. See the Notes for more
+        explanation.
+
 
     Returns
     -------
     q : ndarray of float or complex, optional
-        The orthonormal matrix, of shape (M, K). Only returned if
-        ``mode='full'``.
+        A matrix with orthonormal columns. When mode = 'complete' the
+        result is an orthogonal/unitary matrix depending on whether or not
+        a is real/complex. The determinant may be either +/- 1 in that
+        case.
     r : ndarray of float or complex, optional
-        The upper-triangular matrix, of shape (K, N) with K = min(M, N).
-        Only returned when ``mode='full'`` or ``mode='r'``.
-    a2 : ndarray of float or complex, optional
-        Array of shape (M, N), only returned when ``mode='economic``'.
-        The  diagonal and the upper triangle of `a2` contains `r`, while
-        the rest of the matrix is undefined.
+        The upper-triangular matrix.
+    (h, tau) : ndarrays of np.double or np.cdouble, optional
+        The array h contains the Householder reflectors that generate q
+        along with r. The tau array contains scaling factors for the
+        reflectors. In the deprecated  'economic' mode only h is returned.
 
     Raises
     ------
@@ -580,8 +600,20 @@ def qr(a, mode='full'):
     For more information on the qr factorization, see for example:
     http://en.wikipedia.org/wiki/QR_factorization
 
-    Subclasses of `ndarray` are preserved, so if `a` is of type `matrix`,
-    all the return values will be matrices too.
+    Subclasses of `ndarray` are preserved except for the 'raw' mode. So if
+    `a` is of type `matrix`, all the return values will be matrices too.
+
+    New 'reduced', 'complete', and 'raw' options for mode were added in
+    Numpy 1.8 and the old option 'full' was made an alias of 'reduced'.  In
+    addition the options 'full' and 'economic' were deprecated.  Because
+    'full' was the previous default and 'reduced' is the new default,
+    backward compatibility can be maintained by letting `mode` default.
+    The 'raw' option was added so that LAPACK routines that can multiply
+    arrays by q using the Householder reflectors can be used. Note that in
+    this case the returned arrays are of type np.double or np.cdouble and
+    the h array is transposed to be FORTRAN compatible.  No routines using
+    the 'raw' return are currently exposed by numpy, but some are available
+    in lapack_lite and just await the necessary work.
 
     Examples
     --------
@@ -626,6 +658,20 @@ def qr(a, mode='full'):
     array([  1.1e-16,   1.0e+00])
 
     """
+    if mode not in ('reduced', 'complete', 'r', 'raw'):
+        if mode in ('f', 'full'):
+            msg = "".join((
+                    "The 'full' option is deprecated in favor of 'reduced'.\n",
+                    "For backward compatibility let mode default."))
+            warnings.warn(msg, DeprecationWarning)
+            mode = 'reduced'
+        elif mode in ('e', 'economic'):
+            msg = "The 'economic' option is deprecated.",
+            warnings.warn(msg, DeprecationWarning)
+            mode = 'economic'
+        else:
+            raise ValueError("Unrecognized mode '%s'" % mode)
+
     a, wrap = _makearray(a)
     _assertRank2(a)
     _assertNonEmpty(a)
@@ -653,26 +699,30 @@ def qr(a, mode='full'):
     lwork = int(abs(work[0]))
     work = zeros((lwork,), t)
     results = lapack_routine(m, n, a, m, tau, work, lwork, 0)
-
     if results['info'] != 0:
         raise LinAlgError('%s returns %d' % (routine_name, results['info']))
 
-    #  economic mode. Isn't actually economic.
-    if mode[0] == 'e':
-        if t != result_t :
-            a = a.astype(result_t)
-        return a.T
+    # handle modes that don't return q
+    if mode == 'r':
+        r = _fastCopyAndTranspose(result_t, a[:,:mn])
+        return wrap(triu(r))
 
-    #  generate r
-    r = _fastCopyAndTranspose(result_t, a[:,:mn])
-    for i in range(mn):
-        r[i,:i].fill(0.0)
+    if mode == 'raw':
+        return a, tau
 
-    #  'r'-mode, that is, calculate only r
-    if mode[0] == 'r':
-        return r
+    if mode == 'economic':
+        if t != result_t :
+            a = a.astype(result_t)
+        return wrap(a.T)
 
-    #  from here on: build orthonormal matrix q from a
+    #  generate q from a
+    if mode == 'complete' and m > n:
+        mc = m
+        q = empty((m, m), t)
+    else:
+        mc = mn
+        q = empty((n, m), t)
+    q[:n] = a
 
     if isComplexType(t):
         lapack_routine = lapack_lite.zungqr
@@ -684,20 +734,21 @@ def qr(a, mode='full'):
     # determine optimal lwork
     lwork = 1
     work = zeros((lwork,), t)
-    results = lapack_routine(m, mn, mn, a, m, tau, work, -1, 0)
+    results = lapack_routine(m, mc, mn, q, m, tau, work, -1, 0)
     if results['info'] != 0:
         raise LinAlgError('%s returns %d' % (routine_name, results['info']))
 
     # compute q
     lwork = int(abs(work[0]))
     work = zeros((lwork,), t)
-    results = lapack_routine(m, mn, mn, a, m, tau, work, lwork, 0)
+    results = lapack_routine(m, mc, mn, q, m, tau, work, lwork, 0)
     if results['info'] != 0:
         raise LinAlgError('%s returns %d' % (routine_name, results['info']))
 
-    q = _fastCopyAndTranspose(result_t, a[:mn,:])
+    q = _fastCopyAndTranspose(result_t, q[:mc])
+    r = _fastCopyAndTranspose(result_t, a[:,:mc])
 
-    return wrap(q), wrap(r)
+    return wrap(q), wrap(triu(r))
 
 
 # Eigenvalues
diff --git a/numpy/linalg/tests/test_deprecations.py b/numpy/linalg/tests/test_deprecations.py
new file mode 100644
index 000000000..13d244199
--- /dev/null
+++ b/numpy/linalg/tests/test_deprecations.py
@@ -0,0 +1,24 @@
+"""Test deprecation and future warnings.
+
+"""
+import numpy as np
+from numpy.testing import assert_warns, run_module_suite
+
+
+def test_qr_mode_full_future_warning():
+    """Check mode='full' FutureWarning.
+
+    In numpy 1.8 the mode options 'full' and 'economic' in linalg.qr were
+    deprecated. The release date will probably be sometime in the summer
+    of 2013.
+
+    """
+    a = np.eye(2)
+    assert_warns(DeprecationWarning, np.linalg.qr, a, mode='full')
+    assert_warns(DeprecationWarning, np.linalg.qr, a, mode='f')
+    assert_warns(DeprecationWarning, np.linalg.qr, a, mode='economic')
+    assert_warns(DeprecationWarning, np.linalg.qr, a, mode='e')
+
+
+if __name__ == "__main__":
+    run_module_suite()
diff --git a/numpy/linalg/tests/test_linalg.py b/numpy/linalg/tests/test_linalg.py
index 6750f059d..3f2d438f7 100644
--- a/numpy/linalg/tests/test_linalg.py
+++ b/numpy/linalg/tests/test_linalg.py
@@ -469,11 +469,115 @@ def test_reduced_rank():
 
 
 class TestQR(TestCase):
+
+
+    def check_qr(self, a):
+        # This test expects the argument `a` to be an ndarray or
+        # a subclass of an ndarray of inexact type.
+        a_type = type(a)
+        a_dtype = a.dtype
+        m, n = a.shape
+        k = min(m, n)
+
+        # mode == 'complete'
+        q, r = linalg.qr(a, mode='complete')
+        assert_(q.dtype == a_dtype)
+        assert_(r.dtype == a_dtype)
+        assert_(isinstance(q, a_type))
+        assert_(isinstance(r, a_type))
+        assert_(q.shape == (m, m))
+        assert_(r.shape == (m, n))
+        assert_almost_equal(dot(q, r), a)
+        assert_almost_equal(dot(q.T.conj(), q), np.eye(m))
+        assert_almost_equal(np.triu(r), r)
+
+
+        # mode == 'reduced'
+        q1, r1 = linalg.qr(a, mode='reduced')
+        assert_(q1.dtype == a_dtype)
+        assert_(r1.dtype == a_dtype)
+        assert_(isinstance(q1, a_type))
+        assert_(isinstance(r1, a_type))
+        assert_(q1.shape == (m, k))
+        assert_(r1.shape == (k, n))
+        assert_almost_equal(dot(q1, r1), a)
+        assert_almost_equal(dot(q1.T.conj(), q1), np.eye(k))
+        assert_almost_equal(np.triu(r1), r1)
+
+        # mode == 'r'
+        r2 = linalg.qr(a, mode='r')
+        assert_(r2.dtype == a_dtype)
+        assert_(isinstance(r2, a_type))
+        assert_almost_equal(r2, r1)
+
+
+
     def test_qr_empty(self):
         a = np.zeros((0,2))
         self.assertRaises(linalg.LinAlgError, linalg.qr, a)
 
 
+    def test_mode_raw(self):
+        a = array([[1, 2], [3, 4], [5, 6]], dtype=np.double)
+        b = a.astype(np.single)
+
+        # m > n
+        h1, tau1 = (
+                array([[-5.91607978,  0.43377175,  0.72295291],
+                      [-7.43735744,  0.82807867,  0.89262383]]),
+                array([ 1.16903085,  1.113104  ])
+                )
+        # m > n
+        h2, tau2 = (
+                array([[-2.23606798,  0.61803399],
+                       [-4.91934955, -0.89442719],
+                       [-7.60263112, -1.78885438]]),
+                array([ 1.4472136,  0.       ])
+                )
+
+        # Test double
+        h, tau = linalg.qr(a, mode='raw')
+        assert_(h.dtype == np.double)
+        assert_(tau.dtype == np.double)
+        old_assert_almost_equal(h, h1, decimal=8)
+        old_assert_almost_equal(tau, tau1, decimal=8)
+
+        h, tau = linalg.qr(a.T, mode='raw')
+        assert_(h.dtype == np.double)
+        assert_(tau.dtype == np.double)
+        old_assert_almost_equal(h, h2, decimal=8)
+        old_assert_almost_equal(tau, tau2, decimal=8)
+
+        # Test single
+        h, tau = linalg.qr(b, mode='raw')
+        assert_(h.dtype == np.double)
+        assert_(tau.dtype == np.double)
+        old_assert_almost_equal(h, h1, decimal=8)
+        old_assert_almost_equal(tau, tau1, decimal=8)
+
+
+    def test_mode_all_but_economic(self):
+        a = array([[1, 2], [3, 4]])
+        b = array([[1, 2], [3, 4], [5, 6]])
+        for dt in "fd":
+            m1 = a.astype(dt)
+            m2 = b.astype(dt)
+            self.check_qr(m1)
+            self.check_qr(m2)
+            self.check_qr(m2.T)
+            self.check_qr(matrix(m1))
+        for dt in "fd":
+            m1 = 1 + 1j * a.astype(dt)
+            m2 = 1 + 1j * b.astype(dt)
+            self.check_qr(m1)
+            self.check_qr(m2)
+            self.check_qr(m2.T)
+            self.check_qr(matrix(m1))
+
+
+
+
+
 def test_byteorder_check():
     # Byte order check should pass for native order
     if sys.byteorder == 'little':