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

Ravi Varadhan rvaradhan at jhmi.edu
Thu Oct 18 16:24:48 CEST 2007


Hi,

This is in response to Simon's observation about QR decomp being way too
slow for the "badx" matrix posted by Art Owen.  This is due to the use of
LINPACK routine DQRDC.  QR decomp is much faster when LAPACK routine is
used.  

> system.time(qr(badx, LAPACK=T))
[1] 1.11 0.03 1.14   NA   NA

> system.time(qr(badx))  # Simon's timing
   user  system elapsed
845.896   0.164 846.182

Since LAPACK is more recent and has better routines, I think that it should
be the default for QR decomp (as it is already for SVD).

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 Simon Wood
Sent: Wednesday, October 17, 2007 4:53 AM
To: r-help at r-project.org
Subject: Re: [R] Observations on SVD linpack errors, and a workaround

I get the same problem with this matrix with R 2.6.0 on a 64 bit (Intel
Xeon) machine running Suse 9, but the svd is fine on a 32 bit windows
machine ( R 2.6.0 again). 

`badx' has another weird property (64 bit linux version)... 
> x <- matrix(runif(prod(dim(badx))),nrow(badx),ncol(badx))
> system.time(qr(x)) ## x same dimensions as badx
   user  system elapsed
  1.209   0.031   1.241
> system.time(qr(badx))
   user  system elapsed
845.896   0.164 846.182
Something similar happens on 32 bit windows...  odd for a non-iterative
algorithm (qr.R(qr(badx)) has no NA entries, in case you are wondering, and
has the same singular values as badx, as it should).

best,
Simon

On Wednesday 17 October 2007 07:06, 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.
>
>  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.

-- 
> Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY 
> UK
> +44 1225 386603  www.maths.bath.ac.uk/~sw283

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



More information about the R-help mailing list