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

Ralf Gommers ralf.gommers at googlemail.com
Tue Apr 26 13:19:08 EDT 2011


On Tue, Apr 26, 2011 at 6:36 PM, Jason Grout
<jason-sage at 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]])



More information about the NumPy-Discussion mailing list