被qr.q()所迷惑:“紧”形式的正交矩阵是什么?

内容来源于 Stack Overflow,并遵循CC BY-SA 3.0许可协议进行翻译与使用

  • 回答 (2)
  • 关注 (0)
  • 查看 (23)

r有aqr()函数,它使用Linpack或LAPACK执行QR分解。返回的主要对象是包含在上三角矩阵R中的矩阵“qr”。R=qr[upper.tri(qr)])。到目前为止还不错。QR的下三角部分包含Q“紧凑型”。可以从QR分解中提取q。qr.Q()...。我想找出qr.Q()...。换句话说,我有Q和R,想把它们放在“QR”对象中。r是平凡的,但q不是。目的是要应用于它。qr.solve(),这比solve()在大型系统上。

提问于
用户回答回答于

基本上,OP的问题是,是否有显式计算的Q,人们可以恢复H1 H2...。HT。

我有一个与OP类似的问题,但在不同的上下文中,我的迭代算法需要通过添加列和/或行来改变矩阵A。第一次,QR是用DGEQRF和紧凑的LAPACK格式计算的。在矩阵A发生突变之后,例如,对于新的行,我可以快速地构建一组新的反射器或旋转器,它将湮灭现有R的最低对角线中的非零元素,并生成一个新的R,但现在我有了一组H1。_老H2_老..。HN_老和H1_新H2_新的..。HN_新的(和类似的τ),不能混合成一个QR紧表示。我有两种可能性,也许OP有同样的两种可能性:

  1. 始终保持Q和R显式分开,无论是在第一次计算时还是在每次更新之后,代价是额外的失败,但要保持所需的内存保持良好的限制。
  2. 坚持紧凑的LAPACK格式,但是每次新的更新出现时,都要保留所有这些迷你更新反射器的列表。在解决这个系统的时候,一个人会做一个大的Q‘*c即H1_U3*氢_U3*...*HN_U3*H1_U2*氢_U2*...*HN_U2*H1_U1*氢_U1..。*HN_U1*H1*氢*...*HN*这里UI是QR更新号,这可能是很多乘法要做的事情,内存也要保持跟踪,但绝对是最快的方法。

David给出的长答案基本上解释了压缩QR格式是什么,但没有解释如何使用显式计算Q和R作为输入的压缩QR格式。

用户回答回答于

每个反射矩阵Hi可以用长度-(m-i+1)向量表示.。例如,h1需要长度-m矢量来压缩存储。除了一个项外,所有的向量都放在输入矩阵下三角形的第一列中(R因子使用对角线)。因此,每个反射需要多一个存储标量,这是由一个辅助向量(称为$qraux在R‘s的结果中qr)

使用的紧凑表示在Linpack和LAPACK例程之间是不同的。

Linpack方式

一个House Holder反射被计算为Hi=i-Vivit/pi,其中i是恒等矩阵,pi是$qraux,第六条如下:

  • 六1...I-1=0,
  • 六我=π
  • 六I+1:M=A我+1...我,我

Linpack示例

被分解的矩阵是

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

这个分解是通过计算两个House Holder反射并乘以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方式

一个House Holder反射被计算为hi=i-pivieT,其中i是恒等矩阵,pi是对应的$qraux,第六条如下:

  • 六1...I-1=0,
  • 六我=1
  • 六I+1:M=A我+1...我,我

在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菲尔德,我们会回到那个问题上来。现在我们从以下信息生成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上面。

扫码关注云+社区