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

[Charanpal Dhanjal]

Thanks very much Lou for the information. I tried delving into the C 
code and found a line in the dlasd4_ routine which reads:

for (niter = iter; niter <= MAXITERLOOPS; ++niter) {

This is apparently the main loop for this subroutine and the value of 
MAXITERLOOPS = 100. All I did was increase the maximum number of 
iterations to 200, and this seemed to solve the problem for the matrix 
in question. Let this matrix be called A, then

>>> P0, o0, Q0 = numpy.linalg.svd(A, full_matrices=False)
>>> numpy.linalg.norm((P0*o0).dot(Q0)- A)
1.8558089412794851
>>> numpy.linalg.norm(A)
4.558649005154054
>>> A.shape
(1952, 895)

It seems A has quite a small norm given its dimension, and perhaps this 
explains the error in the SVD (the numpy.linalg.norm((P0*o0).dot(Q0)- A) 
bit).

To investigate a little further I tried finding the SVD of A*1000:

>>> P0, o0, Q0 = numpy.linalg.svd(A*1000, full_matrices=False)
>>> numpy.isfinite(Q0).all()
False
>>> numpy.isfinite(P0).all()
False
>>> numpy.isfinite(o0).all()
False

and hence increasing the number of iterations does not solve the 
problem in this case. That was about as far as I felt I could go with 
investigating the C code. In the meanwhile I will try the squaring the 
matrix solution.

Incidentally, I am confused as to why numpy calls the lapack lite 
routines - when I call numpy.show_config() it seems to have detected my 
ATLAS libraries and I would have expected it to use those.

Charanpal




On Sun, 14 Aug 2011 07:27:06 -0700 (PDT), Lou Pecora wrote:
> 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
> |#=-
> -=#|  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.