# Re: [Rd] Inverting a square matrix using solve() with LAPACK=TRUE (PR#13762)

From: Peter Dalgaard <p.dalgaard_at_biostat.ku.dk>
Date: Thu, 18 Jun 2009 10:14:02 +0300

> 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-devel
```
Received on Thu 18 Jun 2009 - 07:17:44 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 - 13:36:25 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.