sims@Princeton.EDU
2006-Jan-12 02:47 UTC
[Rd] bug in qr.coef() and (therefore) in qr.solve (PR#8476)
[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
Reasonably Related Threads
- Inaccuracy in svd() with R ubuntu package
- solving cubic/quartic equations non-iteratively -- comparisons
- A limitation for polyroot ? (PR#880)
- Bug#403105: lomount not in path - lot of xen utils in /usr/lib/.../bin dir and therefore not accessible
- Error messages when start first time R: "You're using a non-UTF8 locale, therefore only ASCII characters will work."