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

[Douglas Bates]

On 10/18/07, Ravi Varadhan <rvaradhan at jhmi.edu> wrote:
> Dear Prof. Ripley,

> You are right.  I failed to pay attention to the following remark in the
> help page:

> "Note that the storage used by DQRDC and DGEQP3 differs."

> Although this does not matter for solving a linear system (i.e. qr.solve()),
> the Q and R representations are indeed different.

I believe Brian was referring to the remarks about the rank.  The
LINPACK-based code implements a special pivoting scheme that is very
useful in decomposing model matrices.  It is actually a modified
version of the LINPACK code.  See the comments at the beginning of the
file src/appl/dqrdc2.f in the R source tree.  As stated further down
in the documentation

     Using LAPACK (including in the complex case) uses column pivoting
     and does not attempt to detect rank-deficient matrices.


> -----Original Message-----
> From: Prof Brian Ripley [mailto:ripley at stats.ox.ac.uk]
> Sent: Thursday, October 18, 2007 1:42 PM
> To: Ravi Varadhan
> Cc: 'Simon Wood'; r-help at r-project.org
> Subject: Re: [R] Observations on SVD linpack errors, and a workaround
>
> On Thu, 18 Oct 2007, Ravi Varadhan wrote:
>
> > 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).
>
> Please look more carefully at the help page: the default was not changed
> for a very good reason, that the LAPACK and LINPACK versions do not solve
> the same problem and the difference is important for typical R uses of
> qr().  (Hint: look at the value section.)
>
> --
> Brian D. Ripley,                  ripley at stats.ox.ac.uk
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford,             Tel:  +44 1865 272861 (self)
> 1 South Parks Road,                     +44 1865 272866 (PA)
> Oxford OX1 3TG, UK                Fax:  +44 1865 272595
>
> ______________________________________________
> 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.
>