R devel - Jun 2009 - Inverting a square... (PR#13762)

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 erroneous
> 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 at 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.dalgaard at biostat.ku.dk)              FAX: (+45) 35327907

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 Will 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 gives
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 1st 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  =3E 	h5 =3C- hilbert(5)
=3E  =3E 	hinv1 =3C- solve(qr(h5))
=3E  =3E 	hinv2 =3C- solve(qr(h5=2C LAPACK=3DTRUE))	
=3E  =3E 	all=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  =3E 	solve=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 =22LINPACK=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  =3E 	apply(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  =3E 	all=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
=3E    c/ /=27=5F --- Dept=2E of Biostatistics     PO Box 2099=2C 1014 Cph=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