对 qr.Q() 感到困惑：什么是“紧凑”形式的正交矩阵？

2024-01-15

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

介绍

R uses the LINPACK http://www.netlib.org/linpack/ dqrdc http://www.netlib.org/linpack/dqrdc.f routine, by default, or the LAPACK http://www.netlib.org/lapack/ DGEQP3 http://www.netlib.org/lapack/double/dgeqp3.f routine, when specified, for computing the QR decomposition. Both routines compute the decomposition using Householder reflections. An m x n matrix A is decomposed into an m x n economy-size orthogonal matrix (Q) and an n x n upper triangular matrix (R) as A = QR, where Q can be computed by the product of t Householder reflection matrices, with t being the lesser of m-1 and n: Q = H₁H₂...H_t.

Each reflection matrix H_i can be represented by a length-(m-i+1) vector. For example, H₁ requires a length-m vector for compact storage. All but one entry of this vector is placed in the first column of the lower triangle of the input matrix (the diagonal is used by the R factor). Therefore, each reflection needs one more scalar of storage, and this is provided by an auxiliary vector (called $qraux in the result from R's qr).

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

LINPACK 方式

A Householder reflection is computed as H_i = I - v_iv_i^T/p_i, where I is the identity matrix, p_i is the corresponding entry in $qraux, and v_i is as follows:

v_i[1..i-1] = 0,
v_i[i] = p_i
v_i[i+1:m] = A[i+1..m, i] (i.e., a column of the lower triangle of A after calling qr)

LINPACK 示例

让我们通过QR 分解文章中的示例 http://en.wikipedia.org/wiki/QR_decomposition#Example_2在维基百科中的 R.

被分解的矩阵是

> 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

拉帕克之道

A Householder reflection is computed as H_i = I - p_iv_iv_i^T, where I is the identity matrix, p_i is the corresponding entry in $qraux, and v_i is as follows:

v_i[1..i-1] = 0,
v_i[i] = 1
v_i[i+1:m] = A[i+1..m, i] (i.e., a column of the lower triangle of A after calling qr)

在 R 中使用 LAPACK 例程时还有另一个转折：使用列旋转，因此分解正在解决一个不同的相关问题：AP = QR，其中 P 是置换矩阵 http://en.wikipedia.org/wiki/Permutation_matrix.

拉帕克示例

本节执行与前面相同的示例。

> 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场地;我们稍后再讨论这一点。现在我们根据信息生成 QAqr.

> 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 above.

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)