From a1af647b13b33963eca3f9504f1d4b2603c3e4a3 Mon Sep 17 00:00:00 2001 From: Eric Wieser Date: Tue, 7 Nov 2017 20:05:50 -0800 Subject: [PATCH 1/5] MAINT: Remove similar branches from linalg.lstsq --- numpy/linalg/linalg.py | 47 ++++++++++++++++++++++++------------------ 1 file changed, 27 insertions(+), 20 deletions(-) diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py index 26254903d677..041b34b91559 100644 --- a/numpy/linalg/linalg.py +++ b/numpy/linalg/linalg.py @@ -2039,28 +2039,35 @@ def lstsq(a, b, rcond="warn"): 0, work, lwork, iwork, 0) if results['info'] > 0: raise LinAlgError('SVD did not converge in Linear Least Squares') - 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: - 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) + + # undo transpose imposed by fortran-order arrays + b_out = bstar.T + + # b_out contains both the solution and the components of the residuals + x = b_out[:n,:] + r_parts = b_out[n:,:] + if isComplexType(t): + resids = sum(abs(r_parts)**2, axis=-2) else: - x = array(bstar.T[:n,:], dtype=result_t, copy=True) - if results['rank'] == n and m > n: - if isComplexType(t): - resids = sum(abs(bstar.T[n:,:])**2, axis=0).astype( - result_real_t, copy=False) - else: - resids = sum((bstar.T[n:,:])**2, axis=0).astype( - result_real_t, copy=False) + resids = sum(r_parts**2, axis=-2) - st = s[:min(n, m)].astype(result_real_t, copy=True) - return wrap(x), wrap(resids), results['rank'], st + rank = results['rank'] + + # remove the axis we added + if is_1d: + x = x.squeeze(axis=-1) + # we probably should squeeze resids too, but we can't + # without breaking compatibility. + + # as documented + if rank != n or m <= n: + resids = array([], result_real_t) + + # coerce output arrays + s = s.astype(result_real_t, copy=True) + resids = resids.astype(result_real_t, copy=False) # array is temporary + x = x.astype(result_t, copy=True) + return wrap(x), wrap(resids), rank, s def _multi_svd_norm(x, row_axis, col_axis, op): From a311a8d45feaafa8a8ccc15d8e66fec94b7a002d Mon Sep 17 00:00:00 2001 From: Eric Wieser Date: Tue, 7 Nov 2017 22:49:08 -0800 Subject: [PATCH 2/5] MAINT: collect together type mangling --- numpy/linalg/linalg.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py index 041b34b91559..1f7f5f00af4f 100644 --- a/numpy/linalg/linalg.py +++ b/numpy/linalg/linalg.py @@ -1982,7 +1982,11 @@ def lstsq(a, b, rcond="warn"): ldb = max(n, m) if m != b.shape[0]: raise LinAlgError('Incompatible dimensions') + t, result_t = _commonType(a, b) + real_t = _linalgRealType(t) + result_real_t = _realType(result_t) + # Determine default rcond value if rcond == "warn": # 2017-08-19, 1.14.0 @@ -1997,8 +2001,6 @@ def lstsq(a, b, rcond="warn"): 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) bstar[:b.shape[0], :n_rhs] = b.copy() a, bstar = _fastCopyAndTranspose(t, a, bstar) From 6f83089e9a308bea4fd5346cbfe9c55c1590a23c Mon Sep 17 00:00:00 2001 From: Eric Wieser Date: Tue, 7 Nov 2017 23:25:28 -0800 Subject: [PATCH 3/5] DOC: Fix incorrect shape in documentation --- numpy/linalg/linalg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py index 1f7f5f00af4f..8357c92e6687 100644 --- a/numpy/linalg/linalg.py +++ b/numpy/linalg/linalg.py @@ -1915,7 +1915,7 @@ def lstsq(a, b, rcond="warn"): x : {(N,), (N, K)} ndarray Least-squares solution. If `b` is two-dimensional, the solutions are in the `K` columns of `x`. - residuals : {(), (1,), (K,)} ndarray + residuals : {(1,), (K,), (0,)} ndarray Sums of residuals; squared Euclidean 2-norm for each column in ``b - a*x``. If the rank of `a` is < N or M <= N, this is an empty array. From e3a50a964d64375f4058600006c59b9838c5284f Mon Sep 17 00:00:00 2001 From: Eric Wieser Date: Wed, 8 Nov 2017 23:50:29 -0800 Subject: [PATCH 4/5] MAINT: Avoid extra copies in linalg.lstsq This takes gh-5909 a little further. --- numpy/linalg/linalg.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py index 8357c92e6687..5cb381d3e5b7 100644 --- a/numpy/linalg/linalg.py +++ b/numpy/linalg/linalg.py @@ -2002,7 +2002,7 @@ def lstsq(a, b, rcond="warn"): rcond = finfo(t).eps * ldb bstar = zeros((ldb, n_rhs), t) - bstar[:b.shape[0], :n_rhs] = b.copy() + bstar[:m, :n_rhs] = b a, bstar = _fastCopyAndTranspose(t, a, bstar) a, bstar = _to_native_byte_order(a, bstar) s = zeros((min(m, n),), real_t) @@ -2066,9 +2066,9 @@ def lstsq(a, b, rcond="warn"): resids = array([], result_real_t) # coerce output arrays - s = s.astype(result_real_t, copy=True) - resids = resids.astype(result_real_t, copy=False) # array is temporary - x = x.astype(result_t, copy=True) + s = s.astype(result_real_t, copy=False) + resids = resids.astype(result_real_t, copy=False) + x = x.astype(result_t, copy=True) # Copying lets the memory in r_parts be freed return wrap(x), wrap(resids), rank, s From 3402dcf9b4acd32e66a9b8a02ec34a0e02052a5d Mon Sep 17 00:00:00 2001 From: Charles Harris Date: Thu, 9 Nov 2017 09:10:45 -0700 Subject: [PATCH 5/5] STY: Fix PEP8 vertical alignment violation. --- numpy/linalg/linalg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/numpy/linalg/linalg.py b/numpy/linalg/linalg.py index 5cb381d3e5b7..8bc1b14d3f60 100644 --- a/numpy/linalg/linalg.py +++ b/numpy/linalg/linalg.py @@ -2046,7 +2046,7 @@ def lstsq(a, b, rcond="warn"): b_out = bstar.T # b_out contains both the solution and the components of the residuals - x = b_out[:n,:] + x = b_out[:n,:] r_parts = b_out[n:,:] if isComplexType(t): resids = sum(abs(r_parts)**2, axis=-2)