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:
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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
