当前位置： 文江博客 话题详情

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

发布于 2024-09-05 15:24:38 字数 319 浏览 11 评论 0 原文

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

原文

R has a qr() function, which performs QR decomposition using either LINPACK or LAPACK (in my experience, the latter is 5% faster). The main object returned is a matrix "qr" that contains in the upper triangular matrix R (i.e. R=qr[upper.tri(qr)]). So far so good. The lower triangular part of qr contains Q "in compact form". One can extract Q from the qr decomposition by using qr.Q(). I would like to find the inverse of qr.Q(). In other word, I do have Q and R, and would like to put them in a "qr" object. R is trivial but Q is not. The goal is to apply to it qr.solve(), which is much faster than solve() on large systems.

如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论， 你可以免费 注册 一个本站的账号。

评论（2）

你げ笑在眉眼 2024-09-12 15:24:38

简介

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

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

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

LINPACK 方式

Householder 反射计算为 Hi = I - vi viT/pi，其中I是单位矩阵，pi是<中的对应条目code>$qraux，vi如下：

  • vi[1..i-1] = 0,
  • vi sub>[i] = pi
  • vi[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。LAPACK

> 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 反射计算为 Hi = I - piviv iT，其中 I 是单位矩阵，pi 是 $qraux 中的相应条目，v i 如下：

  • vi[1..i-1] = 0,
  • vi[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 中的顺序排列。

Introduction

R uses the LINPACK dqrdc routine, by default, or the LAPACK DGEQP3 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 = H1H2...Ht.

Each reflection matrix Hi can be represented by a length-(m-i+1) vector. For example, H1 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).

The compact representation used is different between the LINPACK and LAPACK routines.

The LINPACK Way

A Householder reflection is computed as Hi = I - viviT/pi, where I is the identity matrix, pi is the corresponding entry in $qraux, and vi is as follows:

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

LINPACK Example

Let's work through the example from the QR decomposition article at Wikipedia in R.

The matrix being decomposed is

> 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

We do the decomposition, and the most relevant portions of the result is shown below:

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

This decomposition was done (under the covers) by computing two Householder reflections and multiplying them by A to get R. We will now recreate the reflections from the information in $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

Now let's verify the Q computed above is correct:

> 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

Looks good! We can also verify QR is equal to 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

The LAPACK Way

A Householder reflection is computed as Hi = I - piviviT, where I is the identity matrix, pi is the corresponding entry in $qraux, and vi is as follows:

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

There is another twist when using the LAPACK routine in R: column pivoting is used, so the decomposition is solving a different, related problem: AP = QR, where P is a permutation matrix.

LAPACK Example

This section does the same example as before.

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

Notice the $pivot field; we will come back to that. Now we generate Q from the information the Aqr.

> 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

Once again, the Q computed above agrees with the R-provided 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

Finally, let's compute QR.

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

Notice the difference? QR is A with its columns permuted given the order in Bqr$pivot above.

回复 原文
带刺的爱情 2024-09-12 15:24:38

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

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

  1. 无论是第一次计算时还是每次更新后都始终保持Q和R显式分离，以额外的触发器为代价，但保持所需的内存有很好的界限。
  2. 坚持紧凑的 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 格式。

I have researched for this same problem as the OP asks and I don't think it is possible. Basically the OP question is whether having the explicitly computed Q, one can recover the H1 H2 ... Ht. I do not think this is possible without computing the QR from scratch but I would also be very interested to know whether there is such solution.

I have a similar issue as the OP but in a different context, my iterative algorithm needs to mutate the matrix A by adding columns and/or rows. The first time, the QR is computed using DGEQRF and thus, the compact LAPACK format. After the matrix A is mutated e.g. with new rows I can quickly build a new set of reflectors or rotators that will annihilate the non-zero elements of the lowest diagonal of my existing R and build a new R but now I have a set of H1_old H2_old ... Hn_old and H1_new H2_new ... Hn_new (and similarly tau's) which can't be mixed up into a single QR compact representation. The two possibilities I have are, and maybe the OP has the same two possibilities:

  1. Always maintain Q and R explicitly separated whether when computed the first time or after every update at the cost of extra flops but keeping the required memory well bounded.
  2. Stick to the compact LAPACK format but then every time a new update comes in, keep a list of all these mini sets of update reflectors. At the point of solving the system one would do a big Q'*c i.e. H1_u3*H2_u3*...*Hn_u3*H1_u2*H2_u2*...*Hn_u2*H1_u1*H2_u1...*Hn_u1*H1*H2*...*Hn*c where ui is the QR update number and this is potentially a lot of multiplications to do and memory to keep track of but definitely the fastest way.

The long answer from David basically explains what the compact QR format is but not how to get to this compact QR format having the explicit computed Q and R as input.

回复 原文
~没有更多了~

关于作者

你在看孤独的风景

暂无简介

0 文章
0 评论
23 人气
发私信

相关话题

更多

热门标签

操作系统 程序设计 IT运维 Linux系统管理 JavaScript 服务器应用 solaris C/C++ PHP Shell BSD Vue.js aix Oracle Python HTML 系统管理 HTML5 CSS 前端
更多

推荐作者

胡图图

文章 0 评论 0

zt006

文章 0 评论 0

z祗昰~

文章 0 评论 0

冰葑

文章 0 评论 0

野の

文章 0 评论 0

天空

文章 0 评论 0

更多

友情链接

文江博客
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站，即表示您同意使用 Cookies 和您的相关数据。
    原文