From: <p.dalgaard_at_biostat.ku.dk>

Date: Thu, 18 Jun 2009 09:20:09 +0200 (CEST)

Date: Thu, 18 Jun 2009 09:20:09 +0200 (CEST)

rvaradhan_at_jhmi.edu wrote:

> Full_Name: Ravi Varadhan

*> Version: 2.8.1
**> OS: Windows
**> Submission from: (NULL) (162.129.251.19)
**>
**>
**> Inverting a matrix with solve(), but using LAPACK=TRUE, gives erroneous
**> results:
*

Thanks, but there seems to be a much easier fix.

Inside coef.qr, we have

coef[qr$pivot, ] <-

.Call("qr_coef_real", qr, y, PACKAGE = "base")[seq_len(p)]

which should be [seq_len(p),]

(otherwise, in the matrix case, the RHS will recycle only the 1st p elements, i.e., the 1st column).

*>
*

> Here is an example:

*>
**> hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }
**> h5 <- hilbert(5)
**> hinv1 <- solve(qr(h5))
**> hinv2 <- solve(qr(h5, LAPACK=TRUE))
**> all.equal(hinv1, hinv2) # They are not equal
**>
**> Here is a function that I wrote to correct this problem:
**>
**> solve.lapack <- function(A, LAPACK=TRUE, tol=1.e-07) {
**> # A function to invert a matrix using "LAPACK" or "LINPACK"
**> if (nrow(A) != ncol(A)) stop("Matrix muxt be square")
**> qrA <- qr(A, LAPACK=LAPACK, tol=tol)
**> if (LAPACK) {
**> apply(diag(1, ncol(A)), 2, function(x) solve(qrA, x))
**> } else solve(qrA)
**> }
**>
**> hinv3 <- solve.lapack(h5)
**> all.equal(hinv1, hinv3) # Now, they are equal
**>
**> ______________________________________________
**> R-devel_at_r-project.org mailing list
**> https://stat.ethz.ch/mailman/listinfo/r-devel
*

-- O__ ---- Peter Dalgaard ุster Farimagsgade 5, Entr.B c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard_at_biostat.ku.dk) FAX: (+45) 35327907 ______________________________________________ R-devel_at_r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-develReceived on Thu 18 Jun 2009 - 07:23:09 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Thu 18 Jun 2009 - 07:36:39 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel.
Please read the posting
guide before posting to the list.
*