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

[Douglas Bates]

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.
>