Schnellster Weg zur Berechnung des Rangs der 2 * 2 Matrix?

Question

Der empfohlene Weg, um den Rang einer Matrix in R zu berechnen scheint qr zu sein:

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 

Ich konnte die Effizienz verbessern, indem diese Funktion für meinen speziellen Fall ändern:

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 Verwenden Kann der Rang einer 2 * 2-Matrix noch schneller berechnet werden?

Answer 1

Sicher, nur mehr Sachen loswerden Sie nicht brauchen (weil Sie wissen, was die Werte sind), haben alle Kontrollen nicht tun, setzen DUP=FALSE, und nur zurückgeben, was Sie wollen:

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 

Ich bin kein Verfechter des Entfernens von Schecks, aber sie verlangsamen die Dinge. x*1.0 und tol*1.0 gewährleisten doppelte, so dass eine Art von Check und fügt ein wenig Overhead. Beachten Sie auch, dass DUP=FALSE potenziell gefährlich sein kann, da Sie die Eingabeobjekte ändern können.

Answer 2

mich jetzt lassen, wenn diese Funktion von einigen Vorsichtsmaßnahmen in diesem Fall fehlt, aber es scheint recht schnell

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

Wir tun sogar RcppEigen besser verwenden, können zu sein.

// [[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(); 
} 

Benchmarks:

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 

Jedoch ist die Toleranz nicht genau das gleiche wie in LINPACK, wegen der unterschiedlichen Toleranz Definitionen:

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 

Die Schwelle in FullPivHouseholderQR ist definiert als:

A Pivot wird als nicht Null betrachtet, wenn sein absoluter Wert streng größer als | pivot | ≤ threshold * | maxpivot | ist wo maxpivot der größte Pivot ist.