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

[Avraham.Adler at guycarp.com]

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"