被qr.Q()困惑：什么是“紧致形式”的正交矩阵？

Question

被qr.Q()困惑：什么是“紧致形式”的正交矩阵？

18

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()要快得多。

- gappy

1

矩阵qr包含上三角矩阵中的因子R，包括对角线。R = qr[upper.tri(qr)]仅返回对角线上方的元素，并且它们不作为矩阵返回。要获取仅包含对角线的上三角形状的矩阵，一种选项是R = qr*upper.tri(qr, diag=TRUE)。 - David Alber

2个回答

0

我已经为与OP提出的同样问题进行了研究，我认为这是不可能的。基本上，OP的问题是：是否可以通过明确计算Q来恢复H1 H2 ... Ht。我认为这是不可能的，除非从头开始计算QR，但我也很想知道是否有这样的解决方案。

我在不同的情境下遇到了类似的问题，我的迭代算法需要通过添加列和/或行来改变矩阵A。第一次使用DGEQRF计算QR，因此使用紧凑的LAPACK格式。在矩阵A被改变后，例如添加新行，我可以快速构建一组新的反射器或旋转器，以消除现有R的最低对角线上的非零元素并构建一个新的R，但现在我有一组H1_old H2_old ... Hn_old和H1_new H2_new ... Hn_new（以及类似的tau），它们不能混合成单个QR紧凑表示。我有两种可能性，也许OP有相同的两种可能性：

始终保持Q和R明确分离，无论是在第一次计算还是在每次更新后，以额外的flops为代价，但保持所需内存的限制。
坚持使用紧凑的LAPACK格式，但每次有新的更新时，保留所有这些小型更新反射器的列表。在解决系统时，将执行一个大的Q'*c，即H1_u3*H2_u3*...*Hn_u3*H1_u2*H2_u2*...*Hn_u2*H1_u1*H2_u1...*Hn_u1*H1*H2*...*Hn*c，其中ui是QR更新编号，这可能需要进行大量的乘法运算和内存跟踪，但绝对是最快的方法。

David的长篇回答基本上解释了紧凑QR格式是什么，但没有说明如何将显式计算的Q和R作为输入转换为此紧凑QR格式。

- SkyWalker

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- David Alber · Accepted Answer

简介

R语言默认使用LINPACK dqrdc例程或指定时使用LAPACK DGEQP3例程来计算QR分解。这两个例程都使用Householder反射来计算分解。一个m x n矩阵A被分解为一个m x n经济尺寸的正交矩阵(Q)和一个n x n上三角矩阵(R)，满足A = QR，其中Q可以通过t个Householder反射矩阵的乘积来计算，其中t是m-1和n中较小的值：Q = H₁H₂...H_t。

每个反射矩阵H_i可以用一个长度为(m-i+1)的向量表示。例如，H₁需要一个长度为m的向量以实现紧凑存储。这个向量除了对角线（由R因子使用）外，其余所有条目都放在输入矩阵的下三角的第一列中。因此，每个反射需要多一个标量的存储空间，并且这由一个辅助向量提供（在R的qr结果中称为$qraux）。

LINPACK和LAPACK例程使用的紧凑表示方式是不同的。

LINPACK方法

计算Householder反射的公式为H_i = I - v_iv_i^T/p_i，其中I是单位矩阵，p_i是$qraux中相应的条目，v_i如下：

v_i[1..i-1] = 0，
v_i[i] = p_i
v_i[i+1:m] = A[i+1..m, i]（即在调用qr后，A的下三角的一列）

LINPACK示例

让我们用R语言来按照维基百科QR分解文章中的示例进行操作。

被分解的矩阵是

> 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的方法

一个Householder变换可以计算为H_i = I - p_iv_iv_i^T，其中I是单位矩阵，p_i 是 $qraux 中对应的条目，v_i 如下：

v_i[1..i-1] = 0,
v_i[i] = 1
v_i[i+1:m] = A[i+1..m, i] （即，在调用qr后A的下三角的一列）

当在R中使用LAPACK例程时，还有一个细节：使用了列置换，因此分解解决的是一个不同但相关的问题：AP = QR，其中P是置换矩阵。

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中列的排列顺序下的置换。