I ran into this problem as well a few months back. The reason for the empty residual array when M==N is that the LAPACK routine for Ax = b puts the solution for x in b. When M>N, the norm-squared is parceled out into the unused (M-N) points in the b array. When M==N, there's no room for the resids. Numpy could always return the optimal norm-squared, but at the expense of another matrix-matrix multiply. The other error (complex norms) can be fixed pretty easily this way. --- linalg.py (revision 6436) +++ linalg.py (working copy) @@ -1316,11 +1316,11 @@ 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) + resids = array([norm(ravel(bstar)[n:], 2)**2]) 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) + resids = array([norm(v[n:], 2)**2 for v in bstar]) st = s[:min(n, m)].copy().astype(_realType(result_t)) return wrap(x), wrap(resids), results['rank'], st Mike On Thu, Feb 19, 2009 at 12:51 PM, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
Hi all,
The summary of ticket 937 is incomplete. It should be "Complex matrices and lstsq".
http://projects.scipy.org/scipy/numpy/ticket/937
Nils _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion