R devel - Jun 2009 - Inverting a square matrix using solve() with LAPACK=TRUE (PR#13762)

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

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