如何对 Matrix 包中的“sparseMatrix”类对象进行 QR 分解？答案

【问题标题】：How do I do a QR decomposition on an object of class "sparseMatrix" in the Matrix package?如何对 Matrix 包中的“sparseMatrix”类对象进行 QR 分解？
【发布时间】：2012-04-06 01:50:36
【问题描述】：

我想在使用B<-as(A, "sparseMatrix") 创建的矩阵上使用Matrix:::qr() 函数进行QR 分解。我知道我可以用Matrix:::qr.R() 得到 R 矩阵。但是，我还需要 Q 矩阵。 Matrix 包中似乎没有 qr.Q() 函数。 Q矩阵如何获得？

【问题讨论】：

标签： r matrix

【解决方案1】：

Q 矩阵实际上存储在V 槽中。当前的 R Matrix 版本似乎包含一个错误 --- 它只是在进行 qr 分解之前神秘地将零行添加到矩阵 a 中。我希望开发人员能来解释一下。因此，以下代码可以帮助您恢复 R 和 Q：

gx.qr.Q <- function(a){
  if(is(a, "qr")){
    return(qr.Q(a, complete = T))
  }else if(is(a, "sparseQR")){
    Q <- diag(nrow = a@V@Dim[1], ncol = a@V@Dim[1])
    for(i in rev(1:a@V@Dim[2])){
      Q <- Q - (a@V[ ,i] * a@beta[i]) %*% (a@V[ ,i] %*% Q)
    }
    return(Q[order(a@p), ][1:a@Dim[1], 1:a@Dim[1]])
  }else{
    stop(gettextf("gx.qr.Q() fails on class '%s'", class(a)[1]))
  }
}

gx.qr.R <- function(a){
  if(is(a, "qr")){
    return(qr.R(a, complete = T)[ ,order(a$pivot)])
  }else if(is(a, "sparseQR")){
    if(length(a@q) == 0){
      return(a@R[1:a@Dim[1], ])
    }else{
      return(a@R[1:a@Dim[1] ,order(a@q)])
    }
  }else{
    stop(gettextf("gx.qr.R() fails on class '%s'", class(a)[1]))
  }
}

我已经通过随机设置矩阵大小和稀疏度进行了测试，它们工作顺利。然而，这是一种“不知道为什么就让它工作”的风格，在这里发布仅供讨论。因为我还没有破解“Matrix”包的实现细节。

【讨论】：

有人提醒我函数不处理“完成 = F”的情况。通过仔细删除一些行或列来解决问题。但是，我仍然对为什么 qr 分解与它的数学定义相比设计得如此复杂感到困惑。