计算2*2矩阵等级的最快方法？

Question

计算R中矩阵等级的推荐方法似乎是qr:

X <- matrix(c(1, 2, 3, 4), ncol = 2, byrow=T)
Y <- matrix(c(1.0, 1, 1, 1), ncol = 2, byrow=T)
qr(X)$rank
[1] 2
qr(Y)$rank
[1] 1

通过针对我的特定情况修改此函数,我能够提高效率:

qr2 <- function (x, tol = 1e-07) { 
  if (!is.double(x)) 
  storage.mode(x) <- "double"
  p <- as.integer(2)
  n <- as.integer(2)
  res <- .Fortran("dqrdc2", qr = x, n, n, p, as.double(tol),
                  rank = integer(1L), qraux = double(p), pivot = as.integer(1L:p), 
                  double(2 * p), PACKAGE = "base")[c(1, 6, 7, 8)]
  class(res) <- "qr"
  res}

qr2(X)$rank
[1] 2
qr2(Y)$rank
[1] 1

library(microbenchmark)
microbenchmark(qr(X)$rank,qr2(X)$rank,times=1000)
Unit: microseconds
         expr    min     lq median     uq      max
1  qr(X)$rank 41.577 44.041 45.580 46.812 1302.091
2 qr2(X)$rank 19.403 21.251 23.099 24.331   80.997

使用R,是否可以更快地计算2*2矩阵的秩？

Answer 1

当然,只需要摆脱更多你不需要的东西(因为你知道它们是什么),不要做任何检查,设置DUP=FALSE,只返回你想要的东西:

qr3 <- function (x, tol = 1e-07) {
  .Fortran("dqrdc2", qr=x*1.0, 2L, 2L, 2L, tol*1.0,
           rank = 0L, qraux = double(2L), pivot = c(1L,2L), 
           double(4L), DUP = FALSE, PACKAGE = "base")[[6L]]
}
microbenchmark(qr(X)$rank,qr2(X)$rank,qr3(X),times=1000)
# Unit: microseconds
#          expr    min      lq  median      uq     max
# 1  qr(X)$rank 33.303 34.2725 34.9720 35.5180 737.599
# 2 qr2(X)$rank 18.334 18.9780 19.4935 19.9240  38.063
# 3      qr3(X)  6.536  7.2100  8.3550  8.5995 657.099

我不是支持取消支票的倡导者,但它们确实减慢了速度. x*1.0并tol*1.0确保双打,所以这是一种检查,增加了一点开销.另请注意DUP=FALSE,由于您可以更改输入对象,因此可能存在危险.