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: [EMAIL PROTECTED]
>
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
> ----------------------------------------------------------------------------
> --------
>
>
> -----Original Message-----
> From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On
> Behalf Of Art Owen
> Sent: Wednesday, October 17, 2007 2:07 AM
> To: r-help@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@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.
>
______________________________________________
R-help@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.
