Fast QR Factorization in R-CodePudding

I have a large number of matrices for which I need to perform a QR factorization and store the resulting Q matrices (normalized such that the R matrix has positive number in its diagonal). Is there another method than using the qr() function?

Here is the working example:

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

So that took 3.07 seconds to compute 62,500 QR factorization. I'm wondering if there is something faster?

CodePudding user response：

If you want:

the R factor to have positive diagonal elements
to explicitly form the Q factor (rather than its sequential Householder vectors format)

you can cheat by using Cholesky factorization:

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)
}

The raw QR:

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)
}

Benchmark:

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

For further speedup, it is often advised that we link our R software to an optimized BLAS library, like OpenBLAS. But more relevant to your context, where you are computing large number of QR factorizations for small matrices, it is more worthwhile to parallelize your for loop.