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."
1 Answers
Answers 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.
0 comments:
Post a Comment