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