計算2 * 2矩陣秩最快的方法是什麼？

推薦的方式來計算R A矩陣的秩似乎是qr：計算2 * 2矩陣秩最快的方法是什麼？

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

來源

2012-08-30 Roland

當然，剛剛擺脫更多的東西你不需要（因爲你知道的值），不做任何檢查，設置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-30 17:13:40

'幸運（98）' - 好吧，我想是4倍。 – BenBarnes

@BenBarnes：我用我保存的時間來看看互聯網上的lolcats。 –

我正在優化一個函數的性能，我需要在模擬中運行幾百萬次。在此函數中，'qr'在while循環中使用。所以，最後這些微秒有時甚至幾個小時。 – Roland

讓我現在如果此功能在這種情況下，一些預防措施的缺乏，但它似乎是相當快速

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-30 22:39:32 Julius

謝謝。你的函數比Joshua的函數更快，但在我的實際測試用例（我沒有在這裏給出）中使用時，似乎不會給出與'qr（X）$ rank'完全相同的結果。發現，何時以及爲什麼會給出不同的結果並不容易。由於你和Joshua功能之間的速度差別不大，我只是拿他的功能。但我贊成你的答案。 – Roland

@Roland，你說的對，我剛剛比較了我的函數和'qr'。 '1e-7'是這裏的問題：對於等級0，我認爲它應該是'== 0'，那麼等級1會出現更多問題，因爲即使所有條目都是關於1e-300 '，這是正確的。但是這樣的條目的產品在R中是0，並且'myrank'返回1，所以這不再是有效的解決方案。劃分行可能工作，但然後功能變得緩慢。 – Julius

我們甚至可以使用更好做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的閾值定義爲：

A如果其絕對值嚴格爲大於| pivot |≤閾值* | maxpivot |，則數據透視將被視爲非零值。其中maxpivot是最大的支點。

來源

2014-02-21 11:02:07 Roland

計算2 * 2矩陣秩最快的方法是什麼？

回答

相關問題