[Numpy-discussion] SVD error in Numpy. NumPy Update reversed?

Lou Pecora lou_boog2000@yahoo....
Wed Mar 19 12:30:45 CDT 2008

I recently had a personal email reply from Damian
Menscher who originally found the error in 2002.  He


I explained the solution in a followup to my own post:
-- in short, find the dlasd4_ routine (for the current
1.0.4 version
it's at numpy/linalg/dlapack_lite.c:21902) and change
the max
iteration count from 20 to 100 or higher.

The basic problem was that they use an iterative
method to converge on
the solution, and they had a cutoff of the max number
of iterations
before giving up (to guard against an infinite loop or
cases where an
unlucky matrix would require an excessive number of
iterations and
therefore CPU).  The fix I used was simply to increase
the max
iteration count (from 20 to 100 -- 50 was enough to
solve my problem
but I went for overkill just to be sure I wouldn't see
it again).  It
*may* be reasonable to just leave this as an infinite
loop, or to
increase the count to 1000 or higher.  A lot depends
on your preferred
failure mode:
  - count too low -> low cpu usage, but "SVD did not
converge" errors
somewhat common
  - very high count -> some matrices will result in
high cpu usage,
non-convergence still possible
  - infinite loop -> it will always converge, but may
take forever

NumPy was supposedly updated also (from 20 to 100, but
you may want to
go higher) in bug 601052.  They said the fix made it
into CVS, but
apparently it got lost or reverted when they did a
release (the oldest
release I can find is v1.0 from 2006 and has it set to
20).  I just
filed another bug (copy/paste of the previous one) in
hopes they'll
fix it for real this time:



I looked at line 21902  of dlapack_lite.c, it is,

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

Indeed the upper limit for iterations in the
linalg.svd code is set for 20.  For now I will go with
my method (on earlier post) of squaring the matrix and
then doing svd when the original try on the original
matrix throws the linalg.linalg.LinAlgError.  I do not
claim that this is a cure-all.  But it seems to work
fast and avoids the original code from thrashing
around in a long iteration. 

I would suggest this be made explicit in the NumPy
documentation and then the user be given the option to
reset the limit on the number of iterations.  

-- Lou Pecora,   my views are my own.

Never miss a thing.  Make Yahoo your home page. 

More information about the Numpy-discussion mailing list