Najszybszy sposób na obliczenie rangi matrycy 2 * 2?

Question

Zalecany sposób obliczyć rangę matrycy w R wydaje się być 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 

udało mi się poprawić efektywność modyfikując tę ​​funkcję na moim konkretnym przypadku:

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 

Korzystanie R , czy można jeszcze szybciej obliczyć rangę matrycy 2 * 2?

Answer 1

Jasne, po prostu pozbyć się więcej rzeczy nie trzeba (bo wiesz, jakie wartości są), nie rób żadnych czeków, ustaw DUP=FALSE, a jedynie zwrócić to, co chcesz:

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 

Nie jestem zwolennikiem usuwania czeków, ale one spowalniają działanie. x*1.0 i tol*1.0 zapewniają podwójne, więc jest to rodzaj czeku i dodaje trochę narzut. Zauważ też, że DUP=FALSE może być potencjalnie niebezpieczny, ponieważ możesz zmienić obiekt wejściowy (s).

Answer 2

Niech mi teraz, czy to funkcja brakuje pewnych środków ostrożności, w tym przypadku, ale wydaje się dość szybko

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

Możemy zrobić jeszcze lepiej przy użyciu 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(); 
} 

Wzorce:

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 

Jednakże tolerancja nie są dokładnie takie same jak w Linpacka powodu różnych definicji tolerancji:

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 

Próg w FullPivHouseholderQR jest zdefiniowany jako:

A pivot zostanie uznany za niezerowy, jeśli jego bezwzględna wartość jest ściśle równa większa niż | pivot | ≤ wartość progowa * | maxpivot | gdzie maxpivot jest największym osią .