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

[Peter Dalgaard]

Art Owen wrote:
> 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.
>   
Nitpicking:

(a) This is precisely *not* a linpack issue (set LINPACK=TRUE and it
goes away).
(b) If you don't know which Fedora you are running, it is probably not
v.8, which is still in testing...

Apart from that, I also see the symptoms on SuSE 10.2 / 64bit. It seems,
however, that svd(badx,0,0)$d does work, i.e. the issue lies in the
calculation of left and right eigenvectors (however both svd(badx,1,0)
and svd(badx,0,1) fail).

>  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.
>   


-- 
   O__  ---- Peter Dalgaard             Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)                  FAX: (+45) 35327907