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: rvarad...@jhmi.edu
Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml
----------------------------------------------------------------------------
--------
-----Original Message-----
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On
Behalf Of avraham.ad...@guycarp.com
Sent: Wednesday, June 17, 2009 11:38 AM
To: r-help@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@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.
