Thank you; that is perfect!
Should this bug be reported somewhere?
--Avraham
"Ravi Varadhan"
<rvarad...@jhmi.e
du> To
<avraham.ad...@guycarp.com>
06/17/2009 01:31 cc
PM <r-help@r-project.org>
Subject
RE: [R] Matrix inversion-different
answers from LAPACK and LINPACK
Avraham,
You can create a function to do the steps that I showed you, so that it is
easy to use:
solve.lapack <- function(A) {
# A function to invert a matrix using "LAPACK"
qrA <- qr(A, LAPACK=TRUE)
apply(diag(1, ncol(A)), 2, function(x) solve(qrA, x))
}
hess <- matrix(runif(100), 10, 10)
hess <- hess + t(hess)
hinv1 <- solve(qr(hess))
hinv2 <- solve.lapack(hess)
all.equal(hinv1, hinv2)
Hope this helps,
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: avraham.ad...@guycarp.com [mailto:avraham.ad...@guycarp.com]
Sent: Wednesday, June 17, 2009 1:10 PM
To: Ravi Varadhan
Cc: r-help@r-project.org
Subject: RE: [R] Matrix inversion-different answers from LAPACK and LINPACK
Thank you VERY much, that was fantastic. I wish I understood WHY your
suggestion works.
To extend that to an n-by-n square matrix, the proper procedure would be
(example, n=5 - most I would usually use (mixture of lognormals)):
Hinv <- matrix(NA, 5, 5)
Hinv[, 1] <- solve(qr(PLLH, LAPACK=TRUE), c(1,0,0,0,0)) Hinv[, 2] <-
solve(qr(PLLH, LAPACK=TRUE), c(0,1,0,0,0)) Hinv[, 3] <- solve(qr(PLLH,
LAPACK=TRUE), c(0,0,1,0,0)) Hinv[, 4] <- solve(qr(PLLH, LAPACK=TRUE),
c(0,0,0,1,0)) Hinv[, 5] <- solve(qr(PLLH, LAPACK=TRUE), c(0,0,0,0,1))
I'm glad that I know the size of the matrix a priori
Also, is this something of which the development list shou;d be made aware?
Once again, thank you very much!
--Avraham
"Ravi Varadhan"
<rvarad...@jhmi.e
du> To
<avraham.ad...@guycarp.com>,
06/17/2009 12:56 <r-help@r-project.org>
PM cc
Subject
RE: [R] Matrix inversion-different
answers from LAPACK and LINPACK
Avraham,
You can make LAPACK work by doing the following:
Hinv[, 1] <- solve(qr(PLLH, LAPACK=TRUE), c(1,0)) Hinv[, 2] <-
solve(qr(PLLH, LAPACK=TRUE), c(0,1))
Here is an example:
H <- matrix(runif(4), 2, 2)
H <- H + t(H)
Hinv <- solve(qr(H)) # this is the correct inverse from LINPACK
Hinv1 <- matrix(NA, 2, 2)
Hinv1[, 1] <- solve(qr(H, LAPACK=TRUE), c(1,0)) Hinv1[, 2] <- solve(qr(H,
LAPACK=TRUE), c(0,1))
Hinv2 <- solve(qr(H, LAPACK=TRUE)) # this won't work, as you found out!
all.equal(Hinv, Hinv1)
all.equal(Hinv, Hinv2)
Hope this helps,
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.
