[Numpy-discussion] Zero row in SVD's unitary matrix on some Mac's
Jason Grout
jason-sage at creativetrax.com
Tue Apr 26 12:07:56 EDT 2011
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
More information about the NumPy-Discussion
mailing list