【问题标题】：Mystified by qr.Q(): what is an orthonormal matrix in "compact" form?被 qr.Q() 迷惑：什么是“紧凑”形式的正交矩阵？
【发布时间】：2010-06-13 05:37:53
【问题描述】：

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

【问题讨论】：

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

标签： r linear-algebra lapack blas

【解决方案1】：

简介

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/pi，其中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 中的 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 方式

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 是 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 中的顺序排列。

【讨论】：

【解决方案2】：

我已经研究了与 OP 要求相同的问题，但我认为这是不可能的。基本上，OP 问题是是否具有明确计算的 Q，是否可以恢复 H1 H2 ... Ht。如果不从头计算 QR，我认为这是不可能的，但我也很想知道是否有这样的解决方案。

我有一个与 OP 类似的问题，但在不同的上下文中，我的迭代算法需要通过添加列和/或行来改变矩阵 A。第一次，QR 是使用 DGEQRF 计算的，因此是紧凑的 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 格式。

【讨论】：