[R] Matrix inversion-different answers from LAPACK and LINPACK

[Ravi Varadhan]

Hi Avraham,

I think this is a bug in solve() and qr.solve().

The structure of the QR object produced by LINPACK and LAPACK are different.


In fact, the help page for qr says:

"qr  a matrix with the same dimensions as x. The upper triangle contains the
R of the decomposition and the lower triangle contains information on the Q
of the decomposition (stored in compact form). Note that the storage used by
DQRDC and DGEQP3 differs."

The qr.coef() function, which is called by solve.qr(), seems to only handle
the QR object produced by LINPACK correctly.  It is incorrect when
LAPACK=TRUE is specified.  The bug might be in the following code segment in
qr.coef:

   z <- .Fortran("dqrcf", as.double(qr$qr), n, k, as.double(qr$qraux), 
        y, ny, coef = matrix(0, nrow = k, ncol = ny), info = integer(1), 
        NAOK = TRUE, PACKAGE = "base")[c("coef", "info")]
 
The Fortran routine "dqrcf" works correctly for QR object produced by
LINPACK, but not for that produced by LAPACK. 

Therefore, a different Fortran subroutine should be called when LAPACK =
TRUE. 


Best,
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_personal_pages/Varadhan.h
tml



----------------------------------------------------------------------------
--------


-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Avraham.Adler at guycarp.com
Sent: Wednesday, June 17, 2009 11:38 AM
To: r-help at r-project.org
Subject: [R] Matrix inversion-different answers from LAPACK and LINPACK


Hello.

I am trying to invert a matrix, and I am finding that I can get different
answers depending on whether I set LAPACK true or false using "qr". I had
understood that LAPACK is, in general more robust and faster than LINPACK,
so I am confused as to why I am getting what seems to be invalid answers.
The matrix is ostensibly the Hessian for a function I am optimizing. I want
to get the parameter correlations, so I need to invert the matrix. There are
times where the standard "solve(X)" returns an error, but "solve(qr(X,
LAPACK=TRUE))" returns values. However, there are times, where the latter
returns what seems to be bizarre results.

For example, an example matrix is PLLH (Pareto LogLikelihood Hessian)

                      alpha                    theta
alpha 1144.6262175141619082 -0.012907777205604828788
theta   -0.0129077772056048  0.000000155437688485563

Running plain "solve (PLLH)" or "solve (qr(PLLH))" returns:

                       [,1]                  [,2]
alpha    0.0137466171688024      1141.53956787721
theta 1141.5395678772131305 101228592.41439932585

However, running "solve(qr(PLLH, LAPACK=TRUE)) returns:

                      [,1]                  [,2]
[1,]    0.0137466171688024    0.0137466171688024
[2,] 1141.5395678772131305 1141.5395678772131305

which seems wrong as the original matrix had identical entries on the off
diagonals.

I am neither a programmer nor an expert in matrix calculus, so I do not
understand why I should be getting different answers using different
libraries to perform the ostensibly same function. I understand the
extremely small value of d²X/d(theta)² (PLLH[2,2]) may be contributing to
the error, but I would appreciate confirmation, or correction, from the
experts on this list.

Thank you very much,

--Avraham Adler



PS: For completeness, the QR decompositions per "R" under LINPACK and LAPACK
are shown below:

> qr(PLLH)
$qr
                          alpha                     theta
alpha -1144.6262175869414932095 0.01290777720653695122277
theta    -0.0000112768491646264 0.00000000987863187747112

$rank
[1] 2

$qraux
[1] 1.99999999993641619511209 0.00000000987863187747112

$pivot
[1] 1 2

attr(,"class")
[1] "qr"
>

> qr(PLLH, LAPACK=TRUE)
$qr
                           alpha                     theta
alpha -1144.62621758694149320945 0.01290777720653695122277
theta    -0.00000563842458249248 0.00000000987863187747112

$rank
[1] 2

$qraux
[1] 1.99999999993642 0.00000000000000

$pivot
[1] 1 2

attr(,"useLAPACK")
[1] TRUE
attr(,"class")
[1] "qr"
______________________________________________
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.