如何获得 R 中稀疏矩阵的真 QR 分解？

Question

我有一个简单的问题：假设我在 R 中有一个稀疏矩阵 A。如果我编码

QR<-qr(A)
Q<-qr.Q(QR)
R<-qr.R(QR)

我没有获得 A 的实际 QR 分解，而是矩阵 PAP_=QR 其中 P 和 P_ 是一些置换矩阵（给定的）。 在大多数情况下，这可能不是问题，但出于我的目的，我相信我绝对需要“真实”矩阵 Q 和 R，使得 A=QR，其中 R 是上三角矩阵。据我所知，这不适用于稀疏矩阵？目前，我正在考虑将我的矩阵转换为稠密矩阵，然后应用 qr() 函数（因为对于稠密矩阵，似乎没有这样的排列）；然而，这似乎有点愚蠢。我想手动编码 QR 分解太复杂了，不是吗？ 任何人都可以帮助解决这个问题吗？

附言我是堆栈溢出的新手，我通常在数学论坛上；我不能在这里用 Latex 写吗？

编辑：我应该提到以上内容适用于 matrix 包中的 qr 分解。显然，我可以使用基本的 qr 函数，但在数字上这似乎不是很好。 （虽然这可能是我最好的选择？）

Answer 1

qr

类

dgCMatrix

的方法将

m

-by-

n

矩阵

A

（带有

m >= n

）分解为：

P1 * A * P2 = Q * R    <==>    A = P1' * Q * R * P2'

其中

P1

和

P2

是置换矩阵。在 C 级别，有一个选项可以禁用列旋转，因此

P2

是单位矩阵。然后你有

A = Q~ R

，其中

Q~ = P1' * Q

是正交的（因为

P1

和

Q

是正交的），你会得到

Q~

与

qr.Q

和

R

与

qr.R

。

不幸的是，在当前Matrix 1.5-4 版本中，该选项未在 R 级别公开。要禁用 R 中的列旋转，您必须直接调用 C 代码（我强烈反对在本示例之外这样做）：

library(Matrix)
set.seed(1)
x <- rsparsematrix(10L, 10L, 0.5)
stopifnot(is(x, "dgCMatrix"), validObject(x),
          packageVersion("Matrix") == "1.5-4")

qr.x <- qr(x)
all.equal(as(qr.Q(qr.x) %*% qr.R(qr.x), "matrix"), as(x, "matrix"))
## [1] "Mean relative difference: 1.059398"

qr.x <- .Call(Matrix:::dgCMatrix_QR, x, 0L, FALSE)
all.equal(as(qr.Q(qr.x) %*% qr.R(qr.x), "matrix"), as(x, "matrix"))
## [1] TRUE

碰巧我在上个月花了很多时间重构Matrix中的因式分解代码。在 Matrix 1.6-0（尚未发布）中，您 will 可以选择禁用列旋转，使用

qr(x, order = 0L)

.

副作用是所有矩阵分解类和方法的文档将得到很大改进。所以，人们可以期待...