¿La forma más rápida para calcular el rango de la matriz 2 * 2?

Question

La forma recomendada para calcular el rango de una matriz en I parece ser 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 

pude mejorar la eficiencia mediante la modificación de esta función para mi caso concreto:

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 

con R , ¿es posible calcular el rango de una matriz 2 * 2 aún más rápido?

Answer 1

Claro, sólo deshacerse de más cosas que no es necesario (porque usted sabe cuáles son los valores), no hacen ninguna comprobación, ajuste DUP=FALSE, y sólo regresan lo que quiere:

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 

No soy partidario de eliminar cheques, pero sí reducen la velocidad. x*1.0 y tol*1.0 aseguran el doble, por lo que es una especie de control y agrega un poco de sobrecarga. También tenga en cuenta que DUP=FALSE puede ser potencialmente peligroso, ya que puede modificar el (los) objeto (s) de entrada.

Answer 2

mí Vamos ahora si esta función carece de algunas precauciones en este caso, pero parece ser bastante rápido

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

Podemos hacer aún mejor usando 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(); 
} 

Puntos de referencia:

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 

Sin embargo, la tolerancia no es exactamente el mismo que en LINPACK, debido a las diferentes definiciones de tolerancia:

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 

El umbral en FullPivHouseholderQR se define como:

A el pivote se considerará distinto de cero si su valor absoluto es estrictamente mayor que | pivote | ≤ umbral * | maxpivot | donde maxpivot es el pivote más grande .