[Rd] 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
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Inverting a square... (PR#13762)
Yes=2C Peter=2E I did look at it=2C but not carefully enought to catch = that=2E Thanks=2C Ravi=2E =5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F= =5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F= =5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F Ravi Varadhan=2C Ph=2ED=2E Assistant Professor=2C Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph=2E (410) 502-2619 email=3A rvaradhan=40jhmi=2Eedu - Original Message - From=3A Peter Dalgaard =3CP=2EDalgaard=40biostat=2Eku=2Edk=3E Date=3A Thursday=2C June 18=2C 2009 9=3A15 am Subject=3A Re=3A =5BRd=5D Inverting a square=2E=2E=2E (PR=2313762) To=3A rvaradhan=40jhmi=2Eedu Cc=3A R-bugs=40r-project=2Eorg =3E Refiling this=2E The actual fix was slightly more complicated=2E Wi= ll soon =3E be committed to R-Patched (aka 2=2E9=2E1 beta)=2E =3E = =3E -p =3E = =3E rvaradhan=40jhmi=2Eedu wrote=3A =3E =3E Full=5FName=3A Ravi Varadhan =3E =3E Version=3A 2=2E8=2E1 =3E =3E OS=3A Windows =3E =3E Submission from=3A (NULL) (162=2E129=2E251=2E19) =3E =3E = =3E =3E = =3E =3E Inverting a matrix with solve()=2C but using LAPACK=3DTRUE=2C g= ives erroneous =3E =3E results=3A =3E = =3E Thanks=2C but there seems to be a much easier fix=2E =3E = =3E Inside coef=2Eqr=2C we have =3E = =3E coef=5Bqr=24pivot=2C =5D =3C- =3E =2ECall(=22qr=5Fcoef=5Freal=22=2C qr=2C y=2C PACKAGE =3D =22base=22= )=5Bseq=5Flen(p)=5D =3E = =3E which should be =5Bseq=5Flen(p)=2C=5D =3E = =3E (otherwise=2C in the matrix case=2C the RHS will recycle only the 1= st p =3E elements=2C i=2Ee=2E=2C the 1st column)=2E =3E = =3E =3E = =3E =3E Here is an example=3A =3E =3E = =3E =3E hilbert =3C- function(n) =7B i =3C- 1=3An=3B 1 / outer(i -= 1=2C i=2C =22+=22) =7D =3E =3Eh5 =3C- hilbert(5) =3E =3Ehinv1 =3C- solve(qr(h5)) =3E =3Ehinv2 =3C- solve(qr(h5=2C LAPACK=3DTRUE)) = =3E =3Eall=2Eequal(hinv1=2C hinv2) =23 They are not equal =3E =3E = =3E =3E Here is a function that I wrote to correct this problem=3A =3E =3E = =3E =3Esolve=2Elapack =3C- function(A=2C LAPACK=3DTRUE=2C tol=3D1=2Ee= -07) =7B =3E =3E=23 A function to invert a matrix using =22LAPACK=22 or =22LIN= PACK=22 =3E =3E if (nrow(A) !=3D ncol(A)) stop(=22Matrix muxt be square= =22) =3E =3E qrA =3C- qr(A=2C LAPACK=3DLAPACK=2C tol=3Dtol) =3E =3E if (LAPACK) =7B =3E =3Eapply(diag(1=2C ncol(A))=2C 2=2C function(x) solve(qrA=2C x)) = =3E =3E =7D else solve(qrA) =3E =3E=7D =3E =3E = =3E =3E hinv3 =3C- solve=2Elapack(h5) =3E =3Eall=2Eequal(hinv1=2C hinv3) =23 Now=2C they are equal =3E =3E = =3E =3E =5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F= =5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F=5F= =5F =3E =3E R-devel=40r-project=2Eorg mailing list =3E =3E = =3E = =3E = =3E -- = =3E O=5F=5F Peter Dalgaard =D8ster Farimagsgade 5=2C= Entr=2EB =3Ec/ /=27=5F --- Dept=2E of Biostatistics PO Box 2099=2C 1014 C= ph=2E K =3E (*) =5C(*) -- University of Copenhagen Denmark Ph=3A (+45)= 35327918 =3E =7E=7E=7E=7E=7E=7E=7E=7E=7E=7E - (p=2Edalgaard=40biostat=2Eku=2Edk)= FAX=3A (+45) 35327907 =3E = =3E __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Error (or inaccuracy) in complex arithmetic (PR#13869)
Dear All,
I have been trying to compute "exact" derivatives in R using the idea of
complex-step derivatives. This is a really, really cool idea. It gives
"exact" derivatives by using a very small, complex step (e.g. 1.e-15i).
Unfortunately, I cannot implement this in R as the "complex arithmetic" in R
is inaccurate.
Here is an example:
#-- Classical Rosenbrock function in n variables
rosen <- function(x) {
n <- length(x)
x1 <- x[2:n]
x2 <- x[1:(n-1)]
sum(100*(x1-x2^2)^2 + (1-x2)^2)
}
x0 <- c(0.0094, 0.7146, 0.2179, 0.6883, 0.5757, 0.9549, 0.7136, 0.0849,
0.4147, 0.4540)
h <- c(1.e-15*1i, 0, 0, 0, 0, 0, 0, 0, 0, 0)
xh <- x0 + h
rx <- rosen(xh)
Re(rx)
Im (rx)
# rx = 190.3079796814885 - 12.13915588266717e-15 i # incorrect imaginary
part in R
However, the imaginary part of the above answer is inaccurate. The correct
imaginary part (from Matlab, S+, Scilab) is:
190.3079796814886 - 4.6677637664e-15 i # correct imaginary part
This inaccuracy is serious enough to affect the acuracy of the compex-step
gradient drastically.
I am wondering whether his problem might be related to the C++ complex.h
library.
I am using Windows XP operating system.
Thanks for taking a look at this.
Best regards,
Ravi.
---
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvarad...@jhmi.edu
Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml
[[alternative HTML version deleted]]
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
