qr.Q()神秘化:什么是"紧凑"形式的正交矩阵？

Question

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

Answer 1

介绍

R 默认使用LINPACK dqrdc例程,或者在指定时使用LAPACK DGEQP3例程来计算QR分解.两个例程都使用Householder反射来计算分解.mxn矩阵A被分解为mxn经济大小的正交矩阵(Q)和nxn上三角矩阵(R)作为A = QR,其中Q可以通过t Householder反射矩阵的乘积来计算,其中t是较小的m-1和n:Q = H ₁ H ₂ ... H _t.

每个反射矩阵H _i可以由长度 - (m-i + 1)向量表示.例如,H ₁需要长度为m的矢量用于紧凑存储.该向量的除了一个条目之外的所有条目都放置在输入矩阵的下三角形的第一列中(对角线由R因子使用).因此,每次反射都需要一个标量的存储空间,这是由一个辅助向量($qraux在R的结果中调用qr)提供的.

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

LINPACK方式

Householder反射计算为H _i = I - v _i v _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](即调用后A的下三角的一列qr)

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 _i v _i v _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](即调用后A的下三角的一列qr)

在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字段; 我们会回过头来看看.现在,我们生成的信息将q 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

再一次,上面计算的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上面的顺序进行置换.