[R] Observations on SVD linpack errors, and a workaround

[Art Owen]

Ravi Varadhan wrote:
> 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.
>   
Save it to a file called badx and then in R from the same
directory say load("badx").

-Art



> ----------------------------------------------------------------------------
> -------
>
> 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
> }
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>