From: <P.Dalgaard_at_biostat.ku.dk>

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

> Inverting a matrix with solve(), but using LAPACK=3DTRUE, gives erroneo=

us

> results:

> Here is an example:

> Here is a function that I wrote to correct this problem:

*> https://stat.ethz.ch/mailman/listinfo/r-devel
*

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Thu 18 Jun 2009 - 13:23:24 GMT

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

Refiling this. The actual fix was slightly more complicated. Will soon
be committed to R-Patched (aka 2.9.1 beta).

-p

rvaradhan_at_jhmi.edu wrote:

> Full_Name: Ravi Varadhan

*> Version: 2.8.1
**> OS: Windows
**> Submission from: (NULL) (162.129.251.19)
*

>=20 >=20

> Inverting a matrix with solve(), but using LAPACK=3DTRUE, gives erroneo=

us

> results:

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

Inside coef.qr, we have

coef[qr$pivot, ] <-

=2ECall("qr_coef_real", qr, y, PACKAGE =3D "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).

>=20

> Here is an example:

>=20

> hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }

> h5 <- hilbert(5)

> hinv1 <- solve(qr(h5))> hinv2 <- solve(qr(h5, LAPACK=3DTRUE))=09> all.equal(hinv1, hinv2) # They are not equal

>=20

> Here is a function that I wrote to correct this problem:

>=20> R-devel@r-project.org mailing list

> solve.lapack <- function(A, LAPACK=3DTRUE, tol=3D1.e-07) {

> # A function to invert a matrix using "LAPACK" or "LINPACK"

> if (nrow(A) !=3D ncol(A)) stop("Matrix muxt be square")> qrA <- qr(A, LAPACK=3DLAPACK, tol=3Dtol)> if (LAPACK) {> apply(diag(1, ncol(A)), 2, function(x) solve(qrA, x))=20> } else solve(qrA)> }

>=20

> hinv3 <- solve.lapack(h5)

> all.equal(hinv1, hinv3) # Now, they are equal

>=20

> ______________________________________________

--=20

O__ ---- Peter Dalgaard =C3=98ster 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 - 13:23:24 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:26 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.
*