[R] Matrix inversion-different answers from LAPACK and LINPACK
Ravi Varadhan
RVaradhan at jhmi.edu
Thu Jun 18 00:50:38 CEST 2009
Avraham,
The hessian calculated by optim() can be inaccurate, since it uses an
inaccurate finite-difference approximation to the second partial
derivatives. A better approach is to use the hessian function, hessian(),
from "numDeriv" package, which uses an accurate, higher-order Richardson's
extrapolation scheme.
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 6:11 PM
To: Douglas Bates
Cc: dmbates at gmail.com; r-help at r-project.org
Subject: Re: [R] Matrix inversion-different answers from LAPACK and LINPACK
I will be the first one to admit I may be doing something stupid, which is
why I am asking here.
Yes, it is supposed to be a V-CoV matrix, but one found by numerical
iteration (a call to optim). I'm actually trying, in a sense, to reproduce
"hessian=TRUE" in cases where I know the analytic form of the first and
second partial derivatives of the distribution to which I am trying a
maximum likelihood fit. So I cannot guarantee that the result will be
positive semi-definite. I wanted to try the supposed increased speed and
precision obtained by passing these to optim using BFGS, for example, as
well as possibly being able to get parameter cross-correlations even in
cases where the simplex result of Optim returns a degenerate hessian. I'm
learning this as I go on using various texts (MASS, Crawley, etc.) and
internet webpages so it is more than likely that my ignorance is making
something go awry. If you can point me to any texts or sources which you
would consider helpful and educational, I would very much appreciate it.
Thank you,
--Avi
Douglas Bates
<bates at stat.wisc.
edu> To
Sent by: Albyn Jones <jones at reed.edu>
dmbates at gmail.com cc
Avraham.Adler at guycarp.com,
r-help at r-project.org
06/17/2009 05:55 Subject
PM Re: [R] Matrix inversion-different
answers from LAPACK and LINPACK
On Wed, Jun 17, 2009 at 2:02 PM, Albyn Jones<jones at reed.edu> wrote:
> 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.
Also, if it is to be a variance-covariance matrix then it must be positive
semidefinite so you should be considering a Cholesky decomposition, not a QR
decomposition.
I think we should insert code to print a warning
"Just because you found a formula involving the inverse of a matrix doesn't
mean that this is a good way to calculate the result - in fact it is usually
a very bad way."
whenever a user asks for
solve(x)
I was corresponding with Tim Davis, an renowned numerical analyst who wrote
the sparse matrix software used in the Matrix package, and he was horrified
that we even allow the one-argument form of the solve function. He said,
"But people will do very stupid things if you provide them with an easy way
of asking for a matrix inverse" and I said, "Yep".
I would amend
> fortune("rethink")
If the answer is parse() you should usually rethink the question.
-- Thomas Lumley
R-help (February 2005)
to say "parse() or solve(x)"
>
> 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.
>>
>
> ______________________________________________
> 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.
>
______________________________________________
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.
More information about the R-help
mailing list