[Numpy-discussion] Summary of ticket 937

[M Trumpis]

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 at 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 at scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>