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:[EMAIL PROTECTED]
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--
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
