diff options
author | Pauli Virtanen <pav@iki.fi> | 2010-10-10 21:29:27 +0200 |
---|---|---|
committer | Pauli Virtanen <pav@iki.fi> | 2010-10-10 21:29:27 +0200 |
commit | 99f79e3f90c09a6b3dc3d3b1a154a99dd1510e28 (patch) | |
tree | d2814f8a3d9d3dae61564531ff1b14430c2b22a8 | |
parent | 8630830bb05bd13aab0fcd869c8339c1515460de (diff) | |
download | numpy-99f79e3f90c09a6b3dc3d3b1a154a99dd1510e28.tar.gz |
BUG: linalg: lstsq should always return real residuals (#937)
-rw-r--r-- | numpy/linalg/linalg.py | 35 | ||||
-rw-r--r-- | numpy/linalg/tests/test_linalg.py | 83 |
2 files changed, 104 insertions, 14 deletions
diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py index eb8c8379a..c044176cf 100644 --- a/numpy/linalg/linalg.py +++ b/numpy/linalg/linalg.py @@ -1680,12 +1680,13 @@ def lstsq(a, b, rcond=-1): """ Return the least-squares solution to a linear matrix equation. - Solves the equation `a x = b` by computing a vector `x` that minimizes - the norm `|| b - a x ||`. The equation may be under-, well-, or over- - determined (i.e., the number of linearly independent rows of `a` can be - less than, equal to, or greater than its number of linearly independent - columns). If `a` is square and of full rank, then `x` (but for round-off - error) is the "exact" solution of the equation. + Solves the equation `a x = b` by computing a vector `x` that + minimizes the Euclidean 2-norm `|| b - a x ||^2`. The equation may + be under-, well-, or over- determined (i.e., the number of + linearly independent rows of `a` can be less than, equal to, or + greater than its number of linearly independent columns). If `a` + is square and of full rank, then `x` (but for round-off error) is + the "exact" solution of the equation. Parameters ---------- @@ -1706,7 +1707,7 @@ def lstsq(a, b, rcond=-1): Least-squares solution. The shape of `x` depends on the shape of `b`. residues : ndarray, shape (), (1,), or (K,) - Sums of residues; squared Euclidean norm for each column in + Sums of residues; squared Euclidean 2-norm for each column in ``b - a*x``. If the rank of `a` is < N or > M, this is an empty array. If `b` is 1-dimensional, this is a (1,) shape array. @@ -1772,6 +1773,7 @@ def lstsq(a, b, rcond=-1): if m != b.shape[0]: raise LinAlgError, 'Incompatible dimensions' t, result_t = _commonType(a, b) + result_real_t = _realType(result_t) real_t = _linalgRealType(t) bstar = zeros((ldb, n_rhs), t) bstar[:b.shape[0],:n_rhs] = b.copy() @@ -1811,16 +1813,27 @@ def lstsq(a, b, rcond=-1): 0, work, lwork, iwork, 0) if results['info'] > 0: raise LinAlgError, 'SVD did not converge in Linear Least Squares' - resids = array([], t) + resids = array([], result_real_t) if is_1d: x = array(ravel(bstar)[:n], dtype=result_t, copy=True) if results['rank'] == n and m > n: - resids = array([sum((ravel(bstar)[n:])**2)], dtype=result_t) + if isComplexType(t): + resids = array([sum(abs(ravel(bstar)[n:])**2)], + dtype=result_real_t) + else: + resids = array([sum((ravel(bstar)[n:])**2)], + dtype=result_real_t) else: x = array(transpose(bstar)[:n,:], dtype=result_t, copy=True) if results['rank'] == n and m > n: - resids = sum((transpose(bstar)[n:,:])**2, axis=0).astype(result_t) - st = s[:min(n, m)].copy().astype(_realType(result_t)) + if isComplexType(t): + resids = sum(abs(transpose(bstar)[n:,:])**2, axis=0).astype( + result_real_t) + else: + resids = sum((transpose(bstar)[n:,:])**2, axis=0).astype( + result_real_t) + + st = s[:min(n, m)].copy().astype(result_real_t) return wrap(x), wrap(resids), results['rank'], st def norm(x, ord=None): diff --git a/numpy/linalg/tests/test_linalg.py b/numpy/linalg/tests/test_linalg.py index 3a9584dd6..a672ed08a 100644 --- a/numpy/linalg/tests/test_linalg.py +++ b/numpy/linalg/tests/test_linalg.py @@ -33,6 +33,11 @@ class LinalgTestCase: b = array([2., 1.], dtype=double) self.do(a, b) + def test_double_2(self): + a = array([[1.,2.], [3.,4.]], dtype=double) + b = array([[2., 1., 4.], [3., 4., 6.]], dtype=double) + self.do(a, b) + def test_csingle(self): a = array([[1.+2j,2+3j], [3+4j,4+5j]], dtype=csingle) b = array([2.+1j, 1.+2j], dtype=csingle) @@ -43,6 +48,11 @@ class LinalgTestCase: b = array([2.+1j, 1.+2j], dtype=cdouble) self.do(a, b) + def test_cdouble_2(self): + a = array([[1.+2j,2+3j], [3+4j,4+5j]], dtype=cdouble) + b = array([[2.+1j, 1.+2j, 1+3j], [1-2j, 1-3j, 1-6j]], dtype=cdouble) + self.do(a, b) + def test_empty(self): a = atleast_2d(array([], dtype = double)) b = atleast_2d(array([], dtype = double)) @@ -70,6 +80,58 @@ class LinalgTestCase: self.do(a, b) +class LinalgNonsquareTestCase: + def test_single_nsq_1(self): + a = array([[1.,2.,3.], [3.,4.,6.]], dtype=single) + b = array([2., 1.], dtype=single) + self.do(a, b) + + def test_single_nsq_2(self): + a = array([[1.,2.], [3.,4.], [5.,6.]], dtype=single) + b = array([2., 1., 3.], dtype=single) + self.do(a, b) + + def test_double_nsq_1(self): + a = array([[1.,2.,3.], [3.,4.,6.]], dtype=double) + b = array([2., 1.], dtype=double) + self.do(a, b) + + def test_double_nsq_2(self): + a = array([[1.,2.], [3.,4.], [5.,6.]], dtype=double) + b = array([2., 1., 3.], dtype=double) + self.do(a, b) + + def test_csingle_nsq_1(self): + a = array([[1.+1j,2.+2j,3.-3j], [3.-5j,4.+9j,6.+2j]], dtype=csingle) + b = array([2.+1j, 1.+2j], dtype=csingle) + self.do(a, b) + + def test_csingle_nsq_2(self): + a = array([[1.+1j,2.+2j], [3.-3j,4.-9j], [5.-4j,6.+8j]], dtype=csingle) + b = array([2.+1j, 1.+2j, 3.-3j], dtype=csingle) + self.do(a, b) + + def test_cdouble_nsq_1(self): + a = array([[1.+1j,2.+2j,3.-3j], [3.-5j,4.+9j,6.+2j]], dtype=cdouble) + b = array([2.+1j, 1.+2j], dtype=cdouble) + self.do(a, b) + + def test_cdouble_nsq_2(self): + a = array([[1.+1j,2.+2j], [3.-3j,4.-9j], [5.-4j,6.+8j]], dtype=cdouble) + b = array([2.+1j, 1.+2j, 3.-3j], dtype=cdouble) + self.do(a, b) + + def test_cdouble_nsq_1_2(self): + a = array([[1.+1j,2.+2j,3.-3j], [3.-5j,4.+9j,6.+2j]], dtype=cdouble) + b = array([[2.+1j, 1.+2j], [1-1j, 2-2j]], dtype=cdouble) + self.do(a, b) + + def test_cdouble_nsq_2_2(self): + a = array([[1.+1j,2.+2j], [3.-3j,4.-9j], [5.-4j,6.+8j]], dtype=cdouble) + b = array([[2.+1j, 1.+2j], [1-1j, 2-2j], [1-1j, 2-2j]], dtype=cdouble) + self.do(a, b) + + class TestSolve(LinalgTestCase, TestCase): def do(self, a, b): x = linalg.solve(a, b) @@ -153,13 +215,28 @@ class TestDet(LinalgTestCase, TestCase): assert_equal(type(linalg.slogdet([[0.0j]])[0]), cdouble) assert_equal(type(linalg.slogdet([[0.0j]])[1]), double) -class TestLstsq(LinalgTestCase, TestCase): +class TestLstsq(LinalgTestCase, LinalgNonsquareTestCase, TestCase): def do(self, a, b): + arr = np.asarray(a) + m, n = arr.shape u, s, vt = linalg.svd(a, 0) x, residuals, rank, sv = linalg.lstsq(a, b) - assert_almost_equal(b, dot(a, x)) - assert_equal(rank, asarray(a).shape[0]) + if m <= n: + assert_almost_equal(b, dot(a, x)) + assert_equal(rank, m) + else: + assert_equal(rank, n) assert_almost_equal(sv, sv.__array_wrap__(s)) + if rank == n and m > n: + expect_resids = (np.asarray(abs(np.dot(a, x) - b))**2).sum(axis=0) + expect_resids = np.asarray(expect_resids) + if len(np.asarray(b).shape) == 1: + expect_resids.shape = (1,) + assert_equal(residuals.shape, expect_resids.shape) + else: + expect_resids = type(x)([]) + assert_almost_equal(residuals, expect_resids) + assert_(np.issubdtype(residuals.dtype, np.floating)) assert imply(isinstance(b, matrix), isinstance(x, matrix)) assert imply(isinstance(b, matrix), isinstance(residuals, matrix)) |