[Numpy-discussion] Zero row in SVD's unitary matrix on some Mac's

[Jason Grout]

On 4/26/11 11:49 AM, Pauli Virtanen wrote:
> Tue, 26 Apr 2011 11:36:19 -0500, Jason Grout wrote:
> [clip]
>> Okay, just one more data point.  Our people that are seeing the problem
>> with numpy returning a non-unitary V also see a non-unitary V being
>> returned by the test C call to zgesdd.  In other words, it really
>> appears that zgesdd follows the Lapack docs, and if rows<columns, the
>> returned V is not necessarily unitary, but may contain a zero row.  This
>> makes numpy's assumptions in using zgesdd false.
>>
>> You can see this report at
>> http://trac.sagemath.org/sage_trac/ticket/11248#comment:25
>
> What LAPACK promises is not very clear. Earlier on the the man page:
>
>         JOBZ    (input) CHARACTER*1
>                 Specifies options for computing all or part of the matrix U:
>                 = 'A':  all M columns of U and all N rows of V**H are  returned
>                 in the arrays U and VT; = 'S':  the first min(M,N) columns of U
>                 and the first min(M,N) rows of V**H are returned in the  arrays
>                 U and VT; = 'O':  If M>= N, the first N columns of U are over‐
>                 written in the array A and all rows of V**H are returned in the
>                 array VT; otherwise, all columns of U are returned in the array
>                 U and the first M rows of V**H are overwritten in the array  A;
>                 = 'N':  no columns of U or rows of V**H are computed.
>
> Looks to me the statement should be parsed as:
>
> 	return_n_rows = (jobz == 'A') or (jobz == 'O' and m>= n)
>
> So the current usage should be OK. (At least, as long as jobz == 'A',
> in the other cases, I don't think the wrapper does the right thing.)
>
> But apparently either there's a bug, or the LAPACK man page needs to
> be understood as you say.


Ah, you're right that it makes sense to parse their statement that way 
too, so I'm not so sure what Lapack really is saying anymore either.  If 
it's parsed the way you propose (which makes sense given the JOBZ 
description), I think it points to a bug in the Accelerate Lapack on the 
affected platforms, as we get the same zero row when we call the 
function directly from C, without numpy, python, or Sage in the middle.

The updated rwork calculation makes no difference with a 3x4 matrix 
(both the old calculation and the new calculation give 66 in the 3x4 
case), so I don't think that is affecting anything.

Jason