r - R 中的 `qr.qy()` 函数

Question

我找不到关于此功能的确切操作的文档。我有一个矩阵 X 的 QR 分解:

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)) 在内部使用函数 qr.qy() 返回 Q，其中 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 是由 qr 表示的 order-nrow(x) 正交(或酉)变换。”

Answer 1

来自 QR 分解的矩阵 Q 仅隐含在函数 qr 的返回值中。列表元素 qr 是 Q 矩阵的紧凑表示；它包含在列表元素 qr 的下三角部分和向量 qraux 中。QR分解的上三角矩阵R是返回值中列表元素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。