[R] Problems understanding qr. family of functions

Douglas Bates bates at cs.wisc.edu
Wed Oct 13 19:40:35 CEST 2004


Kjetil Brinchmann Halvorsen wrote:
> Hola!
> 
> By my understanding of reading   ?qr
> the following shold result a 3 x 3 matrix:  (rw2000)
> 
> x <- matrix( rnorm(12), 4,3)
> y <- matrix( rnorm(12), 4,3)
> xqr <- qr(x)
> yqr <- qr(y)
> 
> qr.qty( xqr, qr.Q(yqr) ) # dim (3,3)
>            [,1]       [,2]      [,3]
> [1,]  0.16815466 -0.1970936 0.2351670
> [2,]  0.15667444  0.4939800 0.8443340
> [3,] -0.97096971  0.1029652 0.1456355
> [4,] -0.06629443 -0.8405570 0.4588976
> 
> # The following should be equal:
> t( qr.Q(xqr) ) %*% qr.Q(yqr)
>           [,1]       [,2]      [,3]
> [1,]  0.1681547 -0.1970936 0.2351670
> [2,]  0.1566744  0.4939800 0.8443340
> [3,] -0.9709697  0.1029652 0.1456355
> 
> but evidently is not.
> 
> There is at least one bug:
> 1) in R
> 2) in the help page
> 3) in my reading of the help page
> 
> Kjetil

You are assuming that qr.Q produces the complete Q matrix.  By default 
it produces only the initial columns of the Q matrix sufficient to 
generate the same size matrix as was decomposed.  Use qr.Q(qrx, complete 
= TRUE) to generate the complete Q matrix.


 > qrx <- qr(matrix(rnorm(12), 4, 3))
 > qr.Q(qrx)
            [,1]       [,2]        [,3]
[1,] -0.6477697 -0.2925012 -0.68404548
[2,]  0.6645340 -0.3817050 -0.33452846
[3,] -0.1567635  0.7073773  0.01126533
[4,] -0.3379560 -0.5180364  0.64810924
 > qr.qy(qrx, diag(4))  ## multiply Q %*% I
            [,1]       [,2]        [,3]      [,4]
[1,] -0.6477697 -0.2925012 -0.68404548 0.1640710
[2,]  0.6645340 -0.3817050 -0.33452846 0.5484401
[3,] -0.1567635  0.7073773  0.01126533 0.6891413
[4,] -0.3379560 -0.5180364  0.64810924 0.4442730
 > qr.Q(qrx, complete = TRUE)
            [,1]       [,2]        [,3]      [,4]
[1,] -0.6477697 -0.2925012 -0.68404548 0.1640710
[2,]  0.6645340 -0.3817050 -0.33452846 0.5484401
[3,] -0.1567635  0.7073773  0.01126533 0.6891413
[4,] -0.3379560 -0.5180364  0.64810924 0.4442730




More information about the R-help mailing list