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

[Albyn Jones]

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