Refiling this. The actual fix was slightly more complicated. Will soon
be committed to R-Patched (aka 2.9.1 beta).
-p
rvarad...@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
> 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
> ______________________________________________
> R-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
--=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.dalga...@biostat.ku.dk) FAX: (+45) 35327907
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
