`qr.qy()` function in R

Question

Я не могу найти документацию относительно точной работы этой функции. У меня есть QR-разложение матрицы 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 

Функция qr(X) дает список:

$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" 

я выбираю диагональные элементы qr(X)$qr, которые я называю 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 

Пока, так хорошо. Теперь следующий звонок я не понимаю:

(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 

---

некоторый прогресс:

Таким образом, благодаря ответу и некоторые чтение, я думаю, что объект qr(X)$qr содержит R полностью в верхнем треугольник:

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

нижний треугольник qr(X)$qr содержит информацию о Q:

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

Как-то вызов qr.Q(qr(X)) возвращает Q, используя внутреннюю функцию qr.qy() с qr() и диагональной матрицей 1 в качестве входных данных.

Но как осуществляется эта операция? Как заполняется остальная часть верхнего верхнего угла Q? Я думаю, что он использует $qraux, но как ему добраться до:

 [,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 

Короче, Как qr.qy() работает конкретно?

Я только что нашел это: «qy.qr(): возвращает результаты умножения матрицы: Q% *% y, где Q - ортогональное (или унитарное) преобразование порядка (x), представленное qr «.

Answer 1

Матрица Q из QR-разложения является только неявной в возвращаемом значении функции qr. Элемент списка qr является компактным представлением Q-матрицы; содержится в нижней треугольной части элемента списка qr и в векторе qraux. Верхняя треугольная матрица R QR-декомпозиции представляет собой верхнюю треугольную часть элемента списка qr в возвращаемом значении.

Функция R qr.qy после некоторых промежуточных шагов в конечном итоге вызывает подпрограмму Lapack dormqr, которая НЕ генерирует матрицу Q явно. Он использует информацию, содержащуюся в элементах списка qr и qraux. См. http://www.netlib.org/lapack/explore-html/da/d82/dormqr_8f.html.

Так qr.qy НЕ преобразует компактную форму Q в фактический Q. Он использует компактную форму для вычисления Q %*% z.

Функция R qr.Q (что вы сделали) использует qr.qy с диагональной матрицей с 1 на диагоналях для генерации Q.

Почему это делается так? По соображениям эффективности.

С помощью следующего кода вы можете проверить это:

library(rbenchmark) 

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

с выходом

         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 

Урок: только генерировать Q, если это действительно необходимо в явном виде, и если вам нужно создать его много с разными входными матрицами. Если Q исправлена ​​и не изменяется, то вы можете использовать Q %*% z.