[Rd] bug in qr.coef() and (therefore) in qr.solve (PR#8476)

From: <sims_at_Princeton.EDU>
Date: Thu 12 Jan 2006 - 02:47:35 GMT


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@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--



R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Thu Jan 12 13:55:06 2006

This archive was generated by hypermail 2.1.8 : Mon 20 Feb 2006 - 03:21:39 GMT