10

在 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 矩阵的秩？

4

3 回答 3

11

当然，只需去掉更多你不需要的东西（因为你知道值是什么），不要做任何检查，设置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可能很危险，因为您可以更改输入对象。

于 2012-08-30T17:13:40.313 回答
3

我们可以使用 RcppEigen 做得更好。

// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
using namespace Rcpp;
using   Eigen::Map;
using   Eigen::MatrixXd;
using   Eigen::FullPivHouseholderQR;
typedef  Map<MatrixXd>  MapMatd;

//calculate rank of a matrix using QR decomposition with pivoting 

// [[Rcpp::export]]
int rankEigen(NumericMatrix  m) {
   const MapMatd  X(as<MapMatd>(m));
   FullPivHouseholderQR<MatrixXd> qr(X);
   qr.setThreshold(1e-7);
   return qr.rank();
}

基准：

microbenchmark(rankEigen(X), qr3(X),times=1000)
Unit: microseconds
         expr   min    lq median    uq    max neval
 rankEigen(X) 1.849 2.465  2.773 3.081 18.171  1000
       qr3(X) 5.852 6.469  7.084 7.392 48.352  1000

但是，公差与 LINPACK 中的公差并不完全相同，因为公差定义不同：

test <- sapply(1:200, function(i) {
  Y <- matrix(c(10^(-i), 10^(-i), 10^(-i), 10^(-i)), ncol = 2, byrow=T)
  qr3(Y) ==  rankEigen(Y)
})

which.min(test)
#[1] 159

FullPivHouseholderQR 中的阈值定义为：

如果一个主元的绝对值严格大于 |pivot|≤ 阈值 * |maxpivot|，则枢轴将被视为非零 其中 maxpivot 是最大的枢轴。

于 2014-02-21T11:02:07.243 回答
2

如果这个功能在这种情况下缺乏一些预防措施，现在让我说，但它似乎相当快

myrank <- function(x)
  if(sum(x^2) < 1e-7) 0 else if(abs(x[1,1]*x[2,2]-x[1,2]*x[2,1]) < 1e-7) 1 else 2

microbenchmark(qr(X)$rank, qr2(X)$rank, qr3(X), myrank(X), times = 1000)
Unit: microseconds
         expr    min     lq median      uq      max
1   myrank(X)  7.466  9.333 10.732 11.1990   97.521
2  qr(X)$rank 52.727 55.993 57.860 62.5260 1237.446
3 qr2(X)$rank 30.329 32.196 33.130 35.4625  178.245
4      qr3(X) 11.199 12.599 13.999 14.9310  116.185

system.time(for(i in 1:10e5) myrank(X))
   user  system elapsed 
   7.46    0.02    7.85 
system.time(for(i in 1:10e5) qr3(X))
   user  system elapsed 
  10.97    0.00   11.85 
system.time(for(i in 1:10e5) qr2(X)$rank)
   user  system elapsed 
  31.71    0.00   33.99 
system.time(for(i in 1:10e5) qr(X)$rank)
   user  system elapsed 
  55.01    0.03   59.73
于 2012-08-30T22:39:32.487 回答