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
<nwag...@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
>
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
