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|≤threshold*|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 に答える