r - 由qr.Q()迷惑:什么是“紧凑”形式的正交矩阵？

Question

R具有qr()函数，该函数使用LINPACK或LAPACK执行QR分解（以我的经验，后者快5％）。返回的主要对象是一个矩阵“ qr”，该矩阵包含在上三角矩阵R中（即R=qr[upper.tri(qr)]）。到目前为止，一切都很好。 qr的下部三角形部分包含“紧凑形式”的Q。可以使用qr.Q()从qr分解中提取Q。我想找到qr.Q()的反函数。换句话说，我确实有Q和R，并希望将它们放在“ qr”对象中。 R是微不足道的，但Q不是。目标是将其应用于qr.solve()，它在大型系统上比solve()快得多。

Answer 1

介绍

R默认情况下使用LINPACK dqrdc例程，或在指定时使用LAPACK DGEQP3例程来计算QR分解。这两个例程都使用Householder反射来计算分解。将mxn矩阵A分解为mxn经济规模正交矩阵（Q）和nxn上三角矩阵（R），如A = QR，其中Q可以通过t家用反射矩阵的乘积来计算，t较小m-1和n的乘积：Q = H1H2 ... Ht。

每个反射矩阵Hi可以由长度-（m-i + 1）向量表示。例如，H1需要长度为m的向量以进行紧凑存储。此向量中除一个外的所有条目都放置在输入矩阵下三角的第一列中（R因子使用对角线）。因此，每个反射都需要一个以上的标量存储，这由一个辅助矢量（在R的$qraux结果中称为qr）提供。

LINPACK和LAPACK例程之间使用的紧凑表示形式有所不同。

LINPACK方式

根据Hi = I-viviT / pi来计算Householder反射，其中I是单位矩阵，pi是$qraux中的相应条目，vi如下所示：


vi [1..i-1] = 0，
vi [i] = pi
vi [i + 1：m] = A [i + 1..m，i]（即，调用qr后A的下三角的列）


LINPACK示例

让我们研究R中的Wikipedia中的example from the QR decomposition article。

分解的矩阵是

> A <- matrix(c(12, 6, -4, -51, 167, 24, 4, -68, -41), nrow=3)
> A
     [,1] [,2] [,3]
[1,]   12  -51    4
[2,]    6  167  -68
[3,]   -4   24  -41

我们进行分解，结果的最相关部分如下所示：

> Aqr = qr(A)
> Aqr
$qr
            [,1]         [,2] [,3]
[1,] -14.0000000  -21.0000000   14
[2,]   0.4285714 -175.0000000   70
[3,]  -0.2857143    0.1107692  -35

[snip...]

$qraux
[1]  1.857143  1.993846 35.000000

[snip...]

通过计算两个Householder反射并将其乘以A得到R（在幕后）进行分解。我们现在将根据 $qr中的信息重新创建反射。

> p = Aqr$qraux   # for convenience
> v1 <- matrix(c(p[1], Aqr$qr[2:3,1]))
> v1
           [,1]
[1,]  1.8571429
[2,]  0.4285714
[3,] -0.2857143

> v2 <- matrix(c(0, p[2], Aqr$qr[3,2]))
> v2
          [,1]
[1,] 0.0000000
[2,] 1.9938462
[3,] 0.1107692

> I = diag(3)   # identity matrix
> H1 = I - v1 %*% t(v1)/p[1]   # I - v1*v1^T/p[1]
> H2 = I - v2 %*% t(v2)/p[2]   # I - v2*v2^T/p[2]

> Q = H1 %*% H2
> Q
           [,1]       [,2]        [,3]
[1,] -0.8571429  0.3942857  0.33142857
[2,] -0.4285714 -0.9028571 -0.03428571
[3,]  0.2857143 -0.1714286  0.94285714

现在，让我们验证上面计算的Q是否正确：

> qr.Q(Aqr)
           [,1]       [,2]        [,3]
[1,] -0.8571429  0.3942857  0.33142857
[2,] -0.4285714 -0.9028571 -0.03428571
[3,]  0.2857143 -0.1714286  0.94285714

看起来挺好的！我们还可以验证QR等于A。

> R = qr.R(Aqr)   # extract R from Aqr$qr
> Q %*% R
     [,1] [,2] [,3]
[1,]   12  -51    4
[2,]    6  167  -68
[3,]   -4   24  -41

LAPACK方式

户主反射的计算公式为Hi = I-piviviT，其中I是单位矩阵，pi是 $qraux中的相应条目，vi如下所示：


vi [1..i-1] = 0，
vi [i] = 1
vi [i + 1：m] = A [i + 1..m，i]（即，调用 qr后A的下三角的列）


在R中使用LAPACK例程时，存在另一种扭曲：使用了列透视，因此分解解决了另一个相关问题：AP = QR，其中P是 permutation matrix。

LAPACK示例

本节与之前的示例相同。

> A <- matrix(c(12, 6, -4, -51, 167, 24, 4, -68, -41), nrow=3)
> Bqr = qr(A, LAPACK=TRUE)
> Bqr
$qr
            [,1]        [,2]       [,3]
[1,] 176.2554964 -71.1694118   1.668033
[2,]  -0.7348557  35.4388886  -2.180855
[3,]  -0.1056080   0.6859203 -13.728129

[snip...]

$qraux
[1] 1.289353 1.360094 0.000000

$pivot
[1] 2 3 1

attr(,"useLAPACK")
[1] TRUE

[snip...]

注意 $pivot字段；我们将回到这一点。现在，我们从 Aqr信息生成Q。

> p = Bqr$qraux   # for convenience
> v1 = matrix(c(1, Bqr$qr[2:3,1]))
> v1
           [,1]
[1,]  1.0000000
[2,] -0.7348557
[3,] -0.1056080


> v2 = matrix(c(0, 1, Bqr$qr[3,2]))
> v2
          [,1]
[1,] 0.0000000
[2,] 1.0000000
[3,] 0.6859203


> H1 = I - p[1]*v1 %*% t(v1)   # I - p[1]*v1*v1^T
> H2 = I - p[2]*v2 %*% t(v2)   # I - p[2]*v2*v2^T
> Q = H1 %*% H2
           [,1]        [,2]       [,3]
[1,] -0.2893527 -0.46821615 -0.8348944
[2,]  0.9474882 -0.01602261 -0.3193891
[3,]  0.1361660 -0.88346868  0.4482655

再次，上面计算的Q与R提供的Q一致。

> qr.Q(Bqr)
           [,1]        [,2]       [,3]
[1,] -0.2893527 -0.46821615 -0.8348944
[2,]  0.9474882 -0.01602261 -0.3193891
[3,]  0.1361660 -0.88346868  0.4482655

最后，让我们计算QR。

> R = qr.R(Bqr)
> Q %*% R
     [,1] [,2] [,3]
[1,]  -51    4   12
[2,]  167  -68    6
[3,]   24  -41   -4

注意区别吗？ QR是A，按上面 Bqr$pivot中的顺序排列其列。