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_la... 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