`qr.qy()` function in R

Question

I can't find documentation as to the exact operation of this function. I have a QR factorization of a matrix X:

X = matrix(c(1,1,1,1,-1.5,-0.5,0.5,1.5,2.25,0.25,0.25, 2.25,-3.275,-0.125,0.125,3.375), 
nrow=4, byrow=F)

     [,1] [,2] [,3]   [,4]
[1,]    1 -1.5 2.25 -3.375
[2,]    1 -0.5 0.25 -0.125
[3,]    1  0.5 0.25  0.125
[4,]    1  1.5 2.25  3.375

The function qr(X) yields a list:

$qr (rounding output)
     [,1]       [,2]          [,3]          [,4]
[1,] -2.0          0          -2.5             0
[2,]  0.5     -2.236             0        -4.583
[3,]  0.5      0.447             2             0 
[4,]  0.5      0.894        -0.929        -1.341

$rank
[1] 4

$qraux
[1] 1.500000 1.000000 1.368524 1.341641

$pivot
[1] 1 2 3 4

attr(,"class")
[1] "qr"

I select the diagonal elements of qr(X)$qr, which I name z:

z = qr(X)$qr
z = Q * (row(Q) == col(Q))

     [,1]      [,2] [,3]      [,4]
[1,]   -2  0.000000    0  0.000000
[2,]    0 -2.236068    0  0.000000
[3,]    0  0.000000    2  0.000000
[4,]    0  0.000000    0 -1.341641

So far, so good. Now the next call I don't understand:

(raw = qr.qy(qr(X), z))

     [,1] [,2] [,3] [,4]
[1,]    1 -1.5    1 -0.3
[2,]    1 -0.5   -1  0.9
[3,]    1  0.5   -1 -0.9
[4,]    1  1.5    1  0.3

---

MAKING SOME PROGRESS:

So, thanks to the answer and some reading, I think that the object qr(X)$qr contains R completely in the upper triangle:

     [,1]       [,2]          [,3]          [,4]
[1,] -2.0          0          -2.5             0
[2,]          -2.236             0        -4.583
[3,]                             2             0 
[4,]                                      -1.341

The lower triangle of qr(X)$qr contains information about Q:

     [,1]       [,2]          [,3]          [,4]
[1,]                                            
[2,]  0.5                                       
[3,]  0.5      0.447                             
[4,]  0.5      0.894        -0.929      

Somehow calling qr.Q(qr(X)) returns Q using internally the function qr.qy() with qr() and a diagonal matrix of 1's as inputs.

But how is this operation carried out? How is the rest of the right upper corner of Q get filled? I think it makes use of $qraux, but how does it get to:

     [,1]       [,2] [,3]       [,4]
[1,] -0.5  0.6708204  0.5  0.2236068
[2,] -0.5  0.2236068 -0.5 -0.6708204
[3,] -0.5 -0.2236068 -0.5  0.6708204
[4,] -0.5 -0.6708204  0.5 -0.2236068

In short, How does qr.qy() work specifically?

I just found this: "qy.qr(): return the results of the matrix multiplications: Q %*% y, where Q is the order-nrow(x) orthogonal (or unitary) transformation represented by qr."

Answer 1

The matrix Q from a QR decomposition is only implicit in the return value of the function qr. The list element qr is a compact representation of the Q matrix; it is contained in the lower triangular part of list element qr and in the vector qraux. The upper triangular matrix R of a QR decomposition is the upper triangular part of the list element qr in the return value.

The R function qr.qy after some intermediate steps eventually calls the Lapack subroutine dormqr, which does NOT generate the Q matrix explicitly. It uses the information contained in list elements qr and qraux. See http://www.netlib.org/lapack/explore-html/da/d82/dormqr_8f.html .

So qr.qy does NOT transform the compact form of Q into the actual Q. It uses the compact form to compute Q %*% z.

The R function qr.Q (which you have done) uses qr.qy with a diagonal matrix with 1's on the diagonals to generate Q.

Why is it done like this? For efficiency reasons.

With the following code you can check this:

library(rbenchmark)

benchmark(qr.qy(XQR,z), {Q <- qr.Q(qr(X)); Q %*% z},  { Q %*% z},
          replications=10000, 
          columns=c("test","replications","elapsed","relative")  )

with output

                                     test replications elapsed relative
3                           {    Q %*% z}        10000   0.022    1.000
2       {    Q <- qr.Q(qr(X));   Q %*% z}        10000   0.486   22.091
1                           qr.qy(XQR, z)        10000   0.152    6.909

Lesson: only generate Q if you really need it in explicit form and if you need to generate it many times with different input matrices. If Q is fixed and doesn't change then you can use Q %*% z.