[Numpy-tickets] [NumPy] #162: problems with numpy.linalg.svd

[NumPy]

#162: problems with numpy.linalg.svd
--------------------------+-------------------------------------------------
 Reporter:  tovrstra      |       Owner:  somebody
     Type:  defect        |      Status:  new     
 Priority:  normal        |   Milestone:          
Component:  numpy.linalg  |     Version:  0.9.8   
 Severity:  major         |    Keywords:          
--------------------------+-------------------------------------------------
 Hi,

 This simple example should always work, since svd is (supposed to be) a
 very stable algorithm.

 {{{
 import numpy
 a = numpy.random.uniform(-1.0, 1.0, (160,260))
 print numpy.linalg.svd(a)
 }}}

 - For numpy-0.9.8 with LAPACK (atlas-3.6.0), the output is
 {{{
  ** On entry to DGESDD parameter number 12 had an illegal value
 }}}
 and the script crashes. For the latest development version of numpy
 (revision 2707), also compiled with LAPACK (atlas-3.6.0), the problem is
 the same

 - When numpy is compiled without LAPACK
 {{{
 export ATLAS=None
 export LAPACK=None
 python setup.py install
 }}}
 version 0.9.8 hangs on the example program with 100% CPU usage. I could
 not interrupt the test program with ctrl-C. The latest development version
 (without LAPACK) works perfect.

 - The problem vanishes for smaller array sizes.

 - To verify the resulting matrices, one can use this code
 {{{
 import numpy

 a = numpy.random.uniform(-1.0, 1.0, (160,260))
 u,s,vh = numpy.linalg.svd(a, full_matrices=True)

 if a.shape[1] > a.shape[0]:
     error1 = sum((numpy.dot(u*s,vh[:a.shape[0]]) - a).ravel()**2)
 else:
     error1 = sum((numpy.dot(u[:,:a.shape[1]]*s,vh) - a).ravel()**2)
 print "error1", error1

 error2 = sum((numpy.dot(u.transpose(), u) - numpy.identity(a.shape[0],
 float)).ravel()**2)
 print "error2", error2

 error3 = sum((numpy.dot(vh.transpose(), vh) - numpy.identity(a.shape[1],
 float)).ravel()**2)
 print "error3", error3
 }}}
 where the three errors should be very small (<1e-20). This verification
 works for the latest development release compiled without lapack.

 - There seems to be a small typo in the svd documentation. The line
 {{{
 a == dot(u,dot(S,vh))
 }}}
 should be replace by
 {{{
 if a.shape[1] > a.shape[0]:
     a == numpy.dot(u*s,vh[:a.shape[0]])
 else:
     a == numpy.dot(u[:,:a.shape[1]]*s,vh)
 }}}

-- 
Ticket URL: <http://projects.scipy.org/scipy/numpy/ticket/162>
NumPy <http://projects.scipy.org/scipy/numpy>
The fundamental package needed for scientific computing with Python.