r - R 中的快速 QR 分解

Question

我有大量矩阵，我需要对其执行 QR 因式分解并存储生成的 Q 矩阵(归一化以使 R 矩阵在其对角线上具有正数)。除了使用 qr() 函数还有其他方法吗？

这是工作示例:

system.time({
  # Parameters for the matrix to be generated
  row_number <- 1e6/4
  col_number <- 4 
  
  # Generate large matrix of random numbers normally distributed. 
  # Basically it's a matrix containing 4x4 matrices for which I will perform the QR factorization:
  RND    <- matrix(data = rnorm(n=1e6 , mean  = 0,sd = 1),nrow = row_number, ncol = col_number)
  
  # Allocate a 0 matrix where the result will be entered 
  QSTACK <- matrix(0, nrow =  row_number, ncol = col_number) 
  
  number_of_blocks <- row_number/4 # The number of 4x4 matrices in RND => 62,500
  
  for (k in c(1:number_of_blocks)) {
    l1             <- 1 + col_number * (k-1)
    l2             <- col_number * k
    QR             <- qr(RND[l1:l2,]) # Perform QR factorization
    R              <- qr.R(QR) 
    QSTACK[l1:l2,] <- qr.Q(QR) %*% diag(sign(diag(R))) # Normalize such that R diagonal elements are positive
  }
})

# user  system elapsed 
# 3.04    0.03    3.07 

所以计算 62,500 个 QR 分解用了 3.07 秒。我想知道是否有更快的东西？

Answer 1

如果你想:

• 具有正对角元素的R因子

• 显式形成 Q 因子(而不是其连续的 Householder 向量格式)

您可以使用 Cholesky 分解作弊:

cheatQR <- function (X) {
  XtX <- crossprod(X)
  R <- chol(XtX)
  Q <- t(forwardsolve(R, t(X), upper.tri = TRUE, transpose = TRUE))
  list(Q = Q, R = R)
}

原始二维码:

rawQR <- function (X) {
  QR <- qr(X)
  Q <- qr.Q(QR)
  R <- qr.R(QR)
  sgn <- sign(diag(R))
  R <- sgn * R
  Q <- Q * rep(sgn, each = nrow(Q))
  list(Q = Q, R = R)
}

基准:

X <- matrix(rnorm(10000 * 4), nrow = 10000, ncol = 4)

ans1 <- rawQR(X)
ans2 <- cheatQR(X)
all.equal(ans1, ans2)
#[1] TRUE

library(microbenchmark)
microbenchmark(rawQR(X), cheatQR(X))
#Unit: microseconds
#       expr      min       lq     mean    median       uq      max neval
#   rawQR(X) 3083.537 3109.222 3796.191 3123.2230 4782.583 13895.81   100
# cheatQR(X)  828.025  837.491 1421.211  865.9085 1434.657 32577.01   100

为了进一步加速，通常建议我们将我们的 R 软件链接到优化的 BLAS 库，例如 OpenBLAS。但与您的上下文更相关，您正在为小矩阵计算大量 QR 因式分解，并行化您的 for 循环更值得。