Avraham.Adler at guycarp.com
2009-Jun-17 15:37 UTC
[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"
Ravi Varadhan
2009-Jun-17 16:34 UTC
[R] Matrix inversion-different answers from LAPACK and LINPACK
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.
Ravi Varadhan
2009-Jun-17 16:56 UTC
[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: 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.
Albyn Jones
2009-Jun-17 19:02 UTC
[R] Matrix inversion-different answers from LAPACK and LINPACK
As you seem to be aware, the matrix is poorly conditioned:> kappa(PLLH,exact=TRUE)[1] 115868900869 It might be worth your while to think about reparametrizing. albyn On Wed, Jun 17, 2009 at 11:37:48AM -0400, Avraham.Adler at guycarp.com wrote:> > 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. >