On Tue, Apr 26, 2011 at 6:36 PM, Jason Grout <jason-sage@creativetrax.com> wrote:
On 4/26/11 11:22 AM, Jason Grout wrote:
On 4/26/11 11:12 AM, Jason Grout wrote:
On 4/26/11 11:07 AM, Jason Grout wrote:
And indeed, I get a 0 row as the last row of the V**H matrix
I just double-checked things one last time and saw that I actually hadn't changed the first argument of zgesdd to "A" in the program that I actually ran. So with this change, I get a nonzero last row of the V**H matrix from the C call to zgesdd. So everything is consistent between the C call to zgesdd and the numpy svd call.
So now my remaining question is: if the Lapack docs only say that V**H is the full n-by-n matrix if M>=N, why is numpy returning it even if M<N?
One more post talking to myself...
I notice that the zgesvd routine docs guarantee that the V returned is unitary, regardless of the size of A. So this might be another argument for calling zgesvd instead of zgesdd.
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
If it is an Accelerate bug it got fixed on OS X 10.6. Below is the output with both current master and 1.5.1. It may be worth checking this on 10.5 with the binary Numpy installer from Sourceforge to make sure it's not the Sage build process causing this somehow. Ralf
import numpy import scipy.linalg A = numpy.array( [[1 - 1j, -3j, -2 + 1j, 1, -2 + 3j], ... [ 1 - 1j, -2 + 1j, 1 + 4j, 0, 2 + 1j], ... [ -1, -5 + 1j, -2 + 1j, 1 + 1j, -5 - 4j], ... [-2 + 4j, 2 - 1j, 8 - 4j, 1 - 8j, 3 - 2j]]) U, S, VH = scipy.linalg.svd(A) VH array([[-0.10378056+0.25534799j, 0.22371070-0.07765616j, 0.55736639-0.40276461j, -0.04874824-0.50177439j, 0.32025785-0.19584279j], [-0.28067375-0.0116824j , -0.49586800+0.1305711j , 0.14432786+0.08803575j, 0.21076793+0.04575091j, -0.19220732-0.73899332j], [ 0.10136055-0.1681669j , 0.13465050+0.69687595j, 0.19709231-0.03497225j, -0.42056091+0.37188099j, 0.28760069-0.14046187j], [-0.27289262+0.66077235j, 0.06177902+0.25464336j, -0.52050549+0.06307174j, 0.11847581-0.00883692j, 0.35771384-0.05719981j], [ 0.49800436+0.21785287j, -0.32014745+0.07789755j, 0.15924098-0.39775924j, 0.50925102+0.33273359j, 0.14247809+0.14849585j]]) numpy.dot(VH, VH.conj().T) array([[ 1.00000000e+00 +1.38777878e-17j, -8.32667268e-17 +2.84494650e-16j, 1.38777878e-16 -1.45716772e-16j, -1.94289029e-16 +1.38777878e-17j, -3.46944695e-17 +2.08166817e-17j], [ -8.32667268e-17 -2.77555756e-16j, 1.00000000e+00 +0.00000000e+00j, -7.56339436e-16 +1.38777878e-16j, 5.55111512e-17 -1.66533454e-16j, 2.11636264e-16 +1.11022302e-16j], [ 1.38777878e-16 +1.17961196e-16j, -7.56339436e-16 -1.28369537e-16j, 1.00000000e+00 +6.93889390e-18j, -3.88578059e-16 +4.16333634e-17j, -1.04083409e-16 +3.12250226e-17j], [ -1.94289029e-16 -1.04083409e-17j, 5.55111512e-17 +1.99493200e-16j, -3.88578059e-16 -2.08166817e-17j, 1.00000000e+00 -6.93889390e-18j, -9.71445147e-17 -8.15320034e-17j], [ -3.46944695e-17 +2.08166817e-17j, 2.11636264e-16 -1.00613962e-16j, -1.04083409e-16 -4.16333634e-17j, -9.71445147e-17 +6.93889390e-17j, 1.00000000e+00 +0.00000000e+00j]])