2 * 2 मैट्रिक्स के रैंक की गणना के लिए सबसे तेज़ तरीका?

Question

आर में एक मैट्रिक्स के पद गणना करने के लिए सुझाया गया तरीका 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 

आर का उपयोग करना , क्या 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 संभावित रूप से खतरनाक हो सकता है, क्योंकि आप इनपुट ऑब्जेक्ट को बदल सकते हैं। इस समारोह इस मामले में कुछ सावधानियों के अभाव है

Answer 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 

Answer 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 सबसे बड़ा पिवट है।