Observations on SVD linpack errors, and a workaround

Up to R version 2.3.0, La.svd had an argument called "method" which allowed
one to choose between "dgesvd" and "dgesdd", but the later versions only use
"dgesdd", which is supposed to be faster for larger matrices.

I wonder if the convergence problem goes away when "dgesvd", which uses the
Lapack routine xBDSQR (uses QR decomposition of a bidiagonal matrix).  The
algorithm "dgesdd" uses the Lapack routine xBDSDC, which uses a divide and
conquer algorithm.

Also, how can I get the "badx" object from your website into my R session?

Ravi.
  

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

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

 

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


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Art Owen
Sent: Wednesday, October 17, 2007 2:07 AM
To: r-help at r-project.org
Subject: [R] Observations on SVD linpack errors, and a workaround


Lately I'm getting this error quite a bit:
Error in La.svd(x, nu, nv) : error code 1 from Lapack routine 'dgesdd'

I'm running R 2.5.0 on a 64 bit Intel machine running Fedora (8 I think).
Maybe the 64 bit platform is more fragile about declaring convergence. 
I'm seeing way more of these errors than I ever have before.

 From R-Help I see that this issue comes up from time to time.

I'm posting an observation that might help diagnose
the problem, and a workaround that improves the odds of success.

I have found that sometimes  svd(t(x)) will work when
svd(x) fails.  For example:

 > load("badx")
 > svd(badx)$d

Error in La.svd(x, nu, nv) : error code 1 from Lapack routine 'dgesdd'

 > svd(t(badx))$d

  [1] 1.572739e+02 9.614579e+01 7.719867e+01 7.127926e+01 6.490623e+01
 .... stuff deleted ....
[126] 8.889272e+00 8.738343e+00 8.447202e+00 8.290393e+00 1.338621e-11
[131] 1.590829e-12 6.154970e-13

badx was a residual matrix, hence the 3 small singular values.
I put the output of save(badx,file="badx") on the web if anybody wants to
play with it.  That matrix is 132 x 10270 entries and the file is
over 10Mb.  As I write this, it seems to be giving firefox a very bad
time loading it.  So proceed with caution (if at all) to the file badx
in the web page stat.stanford.edu/~owen/
There is also a smaller one, called badx2 which illustrates the much
rarer case where a skinny matrix makes svd choke, while its wide
transpose causes no trouble.  Also badx2 did not make firefox hang
so it might be a safer one to look at.

For now my workaround is to write a wrapper that first tries svd(x).
If that fails it then tries svd(t(x)).  In about 800 svds the first case
failed about 100 times.  But the combination never failed.

A simplistic wrapper is listed below.    If SVD failures get very
common for lots of people then a better solution would be to have
the svd function itself try both ways.  Another option is to have the
svd code try the Golub and Reinsch algorithm (or some other SVD)
on those cases where the Lapack one fails.


-Art Owen,  Dept Statistics, Stanford University


##
##   Wrapper function for SVD.  If svd(x) fails, try svd( t(x) ).
##   If both fail you're out of luck in the SVD department.
##   You might succeed by writing a third option based on
##   eigen().  That is numerically inferior to svd when the latter
##   works.
##   -Art Owen, October 2007
##
svdwrapper = function( x, nu, nv, verbose=T ){
#   Caution: I have not tested this much.
#   It's here as an example for an R-Help discussion.
  gotit = F
  try( {svdx = svd(x,nu,nv); gotit=T}, silent = !verbose )
  if( gotit )return(svdx)
  try( {svdtx = svd(t(x),nv,nu); gotit=T}, silent = !verbose )
  if( !gotit )stop("svd(x) and svd(t(x)) both failed.")
  if( verbose )print("svd(x) failed but svd(t(x)) worked.")
  temp    = svdtx$u
  svdtx$u = svdtx$v
  svdtx$v = temp
  svdtx
}