[Numpy-discussion] SVD does not converge on "clean" matrix

Lou Pecora lou_boog2000@yahoo....
Sun Aug 14 09:27:06 CDT 2011


Chuck wrote:


________________________________

Fails here also, fedora 15 64 bits AMD 940. There should be a maximum iterations argument somewhere...

Chuck 

---------------------------------------------------

   ***  Here's the "FIX":

Chuck is right.  There is a max iterations.  Here is a reply from a thread of mine in this group several years ago about this problem and some comments that might help you.

---- From Mr. Damian Menscher who was kind enough to find the iteration location and provide some insight:

Ok, so after several hours of trying to read that code, I found
the parameter that needs to be tuned.  In case anyone has this
problem and finds this thread a year from now, here's your hint:

File: Src/dlapack_lite.c
Subroutine: dlasd4_
Line: 22562

There's a for loop there that limits the number of iterations to
20.  Increasing this value to 50 allows my matrix to converge.
I have not bothered to test what the "best" value for this number
is, though.  In any case, it appears the number just exists to
prevent infinite loops, and 50 isn't really that much closer to
infinity than 20....  (Actually, I'm just going to set it to 100
so I don't have to think about it ever again.)

Damian Menscher
-- 
-=#| Physics Grad Student & SysAdmin @ U Illinois Urbana-Champaign |#=-
-=#| 488 LLP, 1110 W. Green St, Urbana, IL 61801 Ofc:(217)333-0038 |#=-
-=#| 1412 DCL, Workstation Services Group, CITES Ofc:(217)244-3862 |#=-
-=#| <menscher at uiuc.edu> www.uiuc.edu/~menscher/ Fax:(217)333-9819 |#=-


---- My reply and a "fix" of sorts without changing the hard coded iteration max:

I have looked in Src/dlapack_lite.c and line 22562 is no longer a line that sets a max. iterations parameter.  There are several set in the file, but that code is hard to figure (sort of a Fortran-in-C hybrid).  

Here's one, for example:

    maxit = *n * 6 * *n;   // Line 887

I have no idea which parameter to tweak.  Apparently this error is still in numpy (at least to my version). Does anyone have a fix?  Should I start a ticket (I think this is what people do)?  Any help appreciated.

I'm using a Mac Book Pro (Intel chip), system 10.4.11, Python 2.4.4.

============ Possible try/except =========================== 

#  A is the original matrix
try:
    U,W,VT=linalg.svd(A)
except linalg.linalg.LinAlgError:  # "Square" the matrix and do SVD

    print "Got svd except, trying square of A."
    A2=dot(conj(A.T),A)
    U,W2,VT=linalg.svd(A2)


This works so far.

---------------------------------------------------------------------------------------

I've been using that simple "fix" of "squaring" the original matrix for several years and it's worked every time.  I'm not sure why.  It was just a test and it worked.  

You could also change the underlying C or Fortran code, but you then have to recompile everything in numpy.  I wasn't that brave.


-- Lou Pecora, my views are my own.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/numpy-discussion/attachments/20110814/b8b649ec/attachment.html 


More information about the NumPy-Discussion mailing list