[Rd] bug in qr.coef() and (therefore) in qr.solve (PR#8476)
sims@Princeton.EDU
sims at Princeton.EDU
Thu Jan 12 03:47:35 CET 2006
This is a multi-part message in MIME format.
--------------050206000203080003040803
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit
[I thought I'd submitted this bug report some time ago, but it's never showed up on the bug tracking system, so I'm submitting again.]
qr.solve() gives incorrect results when dealing with complex matrices or with qr objects that have been computed with LAPACK=TRUE, whenever the b argument has more than one column. This bug flows from qr.coef(), which has a similar problem. I believe the problem is the line
coef[qr$pivot, ] <- .Call("qr_coef_cmplx", qr, y, PACKAGE = "base")[1:p]
and the similar line in the LAPACK section of the qr.coef() code. As far as I can see, this line should read
coef <- .Call("qr_coef_cmplx", qr, y, PACKAGE = "base")[qr$pivot,]
With this change, qr.coef() gives correct results for my examples. In the examples, the qr.coeffCS() function is my version of qr.coef, with the two lines
in question changed as above (and no other modifications).
Examples:
> A <- matrix(rnorm(9),3,3)
> B <- matrix(rnorm(9),3,3)
> solve(A+1i*B,A+1i*B)
[,1] [,2] [,3]
[1,] 1.000000e+00+0.000000e+00i -1.853360e-17-1.199306e-17i 0+0i
[2,] 2.338819e-17-1.192988e-19i 1.000000e+00+1.155338e-20i 0+0i
[3,] -6.940188e-18+1.120842e-17i 5.188659e-17-3.226848e-17i 1+0i
> qr.solve(A+1i*B,A+1i*B)
[,1] [,2]
[1,] 1.000000e-00-2.583088e-16i 1.000000e-00-2.583088e-16i
[2,] -1.045057e-16+1.979352e-16i -1.045057e-16+1.979352e-16i
[3,] 3.966684e-16+7.360601e-16i 3.966684e-16+7.360601e-16i
[,3]
[1,] 1.000000e-00-2.583088e-16i
[2,] -1.045057e-16+1.979352e-16i
[3,] 3.966684e-16+7.360601e-16i
## Note: all columns the same, matrix should be the identity
> qr.coef(qr(A+1i*B),A+1i*B)
[,1] [,2]
[1,] 1.000000e-00-2.583088e-16i 1.000000e-00-2.583088e-16i
[2,] -1.045057e-16+1.979352e-16i -1.045057e-16+1.979352e-16i
[3,] 3.966684e-16+7.360601e-16i 3.966684e-16+7.360601e-16i
[,3]
[1,] 1.000000e-00-2.583088e-16i
[2,] -1.045057e-16+1.979352e-16i
[3,] 3.966684e-16+7.360601e-16i
> qr.coeffCS(qr(A+1i*B),A+1i*B)
[,1] [,2]
[1,] 1.000000e-00-2.583088e-16i -6.306738e-17-8.893088e-17i
[2,] -1.045057e-16+1.979352e-16i 1.000000e-00-1.921467e-17i
[3,] 3.966684e-16+7.360601e-16i 4.149193e-17-1.149122e-16i
[,3]
[1,] -7.350360e-17-6.340421e-17i
[2,] -2.685509e-17+3.225516e-17i
[3,] 1.000000e+00-8.587603e-18i
## Note: correct results from the function with two modified lines, whether
## LAPACK=TRUE or not and whether complex or not.
> qr.coeffCS(qr(A,LAPACK=TRUE),A)
[,1] [,2] [,3]
[1,] 1.000000e+00 0.000000e+00 0
[2,] 1.267908e-15 1.000000e+00 0
[3,] -1.889451e-15 4.074224e-17 1
> qr.coef(qr(A,LAPACK=TRUE),A)
[,1] [,2] [,3]
[1,] 1.000000e+00 1.000000e+00 1.000000e+00
[2,] 1.267908e-15 1.267908e-15 1.267908e-15
[3,] -1.889451e-15 -1.889451e-15 -1.889451e-15
> qr.solve(qr(A,LAPACK=TRUE),A)
[,1] [,2] [,3]
[1,] 1.000000e+00 1.000000e+00 1.000000e+00
[2,] 1.267908e-15 1.267908e-15 1.267908e-15
[3,] -1.889451e-15 -1.889451e-15 -1.889451e-15
> solve(A,A)
[,1] [,2] [,3]
[1,] 1.000000e+00 2.628787e-17 0
[2,] 1.899954e-17 1.000000e+00 0
[3,] 2.851722e-16 -6.860465e-17 1
> C <- solve(A,B)
> qr.solve(A,B)-C
[,1] [,2] [,3]
[1,] -8.881784e-16 1.332268e-15 4.440892e-16
[2,] -1.387779e-15 2.220446e-15 8.881784e-16
[3,] 2.664535e-15 -3.552714e-15 -1.110223e-15
## qr.solve() matches solve() with real, non-LAPACK argument
> qr.coef(qr(A),B)-C
[,1] [,2] [,3]
[1,] -8.881784e-16 1.332268e-15 4.440892e-16
[2,] -1.387779e-15 2.220446e-15 8.881784e-16
[3,] 2.664535e-15 -3.552714e-15 -1.110223e-15
> qr.coef(qr(A,LAPACK=TRUE),B)-C
[,1] [,2] [,3]
[1,] 1.110223e-15 3.380377 1.4697337
[2,] 2.164935e-15 1.812489 0.9821958
[3,] -3.552714e-15 -9.280706 -6.3293155
## qr.coef() gives different results with LAPACK=TRUE
> qr.coeffCS(qr(A,LAPACK=TRUE),B)-C
[,1] [,2] [,3]
[1,] 1.110223e-15 -1.776357e-15 -3.330669e-16
[2,] 2.164935e-15 -3.996803e-15 -6.661338e-16
[3,] -3.552714e-15 7.105427e-15 1.221245e-15
## the modified function gives the same results for LAPACK=TRUE as does
## qr.coef() in the LAPACK=FALSE case.
## lines below just show that qr.coeffCS() works in non-square cases and
## that the problem is there in ths same form in these cases for the original
## qr.coef()
> X <- matrix(rnorm(36),12,3)
> y <- matrix(rnorm(24),12,2)
> b <- qr.coef(qr(X),y)
> qr.coef(qr(X,LAPACK=TRUE),y)-b
[,1] [,2]
[1,] -5.551115e-17 0.2509164
[2,] 1.110223e-16 0.1846802
[3,] 0.000000e+00 1.0224349
> qr.coef(qr(X+0i),y)-b
[,1] [,2]
[1,] -5.551115e-17+0i 0.2509164+0i
[2,] 1.110223e-16+0i 0.1846802+0i
[3,] 0.000000e+00+0i 1.0224349+0i
> qr.coeffCS(qr(X+0i),y)-b
[,1] [,2]
[1,] -5.551115e-17+0i 2.775558e-17+0i
[2,] 1.110223e-16+0i 1.110223e-16+0i
[3,] 0.000000e+00+0i 5.551115e-17+0i
--------------050206000203080003040803
Content-Type: text/x-vcard; charset=utf-8;
name="sims.vcf"
Content-Transfer-Encoding: 7bit
Content-Disposition: attachment;
filename="sims.vcf"
begin:vcard
fn:Chris Sims
n:Sims;Chris
org:Princeton University;Department of Economics
adr:;;Fisher Hall;Princeton;NJ;08544-1021;USA
email;internet:sims at princeton.edu
tel;work:609 258 4033
tel;fax:609 258 6419
x-mozilla-html:FALSE
url:http://www.princeton.edu/~sims
version:2.1
end:vcard
--------------050206000203080003040803--
More information about the R-devel
mailing list