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

[Jason Grout]

On 4/25/11 12:57 PM, Pauli Virtanen wrote:
> The Numpy routine is a very thin wrapper of LAPACK's ZGESDD, and probably
> cannot have any bugs of this kind,

As noted in my other message, I've been digging through the ZGESDD docs 
to understand it better.  Here is the doc string for what becomes the 
V^T matrix:

*  VT      (output) COMPLEX*16 array, dimension (LDVT,N)
*          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
*          N-by-N unitary matrix V**H;
*          if JOBZ = 'S', VT contains the first min(M,N) rows of
*          V**H (the right singular vectors, stored rowwise);
*          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.

Notice that VT is supposed to contain the N-by-N V**H if M>=N (i.e., 
more rows than columns).  In our problem cases, we have M<N.  I looked 
at the numpy linalg code and it doesn't seem to check to see if M<N 
before returning the V matrix.  Is this a problem?

When I run an example C program that links to the OSX Accelerate 
framework, I get a V matrix that has a zero row at the bottom.  I got 
these results by:

1. Download the example found at 
http://software.intel.com/sites/products/documentation/hpc/mkl/lapack/mkl_lapack_examples/zgesdd_ex.c.htm

2. Change the zgesdd call to have first argument "A" (so it mimics the 
default numpy call)

3. Change

print_matrix( "Right singular vectors (stored rowwise)", m, n, vt, ldvt );

to

print_matrix( "Right singular vectors (stored rowwise)", n, n, vt, ldvt );

(this is to print out the full n-by-n matrix, instead of just the first 
m rows)

4. Compile the program with "gcc zgesdd_ex.c -framework Accelerate"

5. Run the program: "./a.out"

And indeed, I get a 0 row as the last row of the V**H matrix

However, when I do the SVD of this same matrix using numpy (which I 
*think* uses the Accelerate framework the way we compile it in Sage), I 
get a full V matrix:

In [13]: I=1j
In [14]: a=numpy.asarray([( -5.40+ I*  7.40), (  6.00+ I*  6.38), ( 
9.91+ I*  0.16), ( -5.28+ I* -4.16),   (  1.09+ I*  1.55), (  2.60+ I* 
0.07), (  3.98+ I* -5.26), (  2.03+ I*  1.11),   (  9.88+ I*  1.91), ( 
4.92+ I*  6.31), ( -2.11+ I*  7.39), ( -9.81+ I* 
-8.98)],dtype=complex).reshape(3,4)

In [15]: numpy.linalg.svd(a)Out[15]:
(array([[ 0.54742764+0.j        ,  0.76302168+0.j        , 
-0.34368721+0.j        ],
        [-0.03507684-0.15148438j,  0.27097680-0.22637514j,
          0.54572628-0.74386208j],
        [ 0.81299016+0.12325614j, -0.52311095-0.13956616j,
          0.13357577-0.1135282j ]]),
  array([ 21.75519279,  16.59545017,   3.97327576]),
  array([[ 0.23160531+0.20669796j,  0.36590896+0.3864613j ,
          0.24259328+0.32833854j, -0.56133932-0.37233547j],
        [-0.57911906+0.40329699j,  0.10921398+0.17242422j,
          0.59673801-0.27492812j,  0.15998810+0.05510835j],
        [ 0.60420072+0.12337134j, -0.18988750+0.29722068j,
          0.39210635+0.19697635j,  0.45451433+0.31015037j],
        [-0.08335334+0.13557583j,  0.69622499-0.25685378j,
         -0.15350399+0.43075518j,  0.46295826+0.02290113j]]))

That seems odd.  I must be missing something.  It also seems odd that 
numpy returns V as a full 4x4 matrix even when M<N.

Thanks,

Jason