计算 2*2 矩阵秩的最快方法？

Posted 2023-02-16

【中文标题】计算 2*2 矩阵秩的最快方法？

【英文标题】：Fastest way for calculating rank of 2*2 matrix?

【发布时间】：2012-08-25 11:09:46

【问题描述】：

在 R 中计算矩阵秩的推荐方法似乎是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

使用 R，是否可以更快地计算 2*2 矩阵的秩？

【参考方案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 可能具有潜在危险，因为您可以更改输入对象。

【讨论】：

fortune(98) - 好吧，我想是 4 倍。

@BenBarnes：我用节省下来的时间看互联网上的笑笑猫。

我正在优化一个函数的性能，我需要在模拟中运行几百万次。 qr 在此函数的 while 循环内使用。因此，最终这些微秒数可能长达数小时。

【参考方案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 

【讨论】：

谢谢。您的函数比 Joshua 的要快，但在我的实际测试用例（我没有在这里给出）中使用时，似乎给出的结果与 qr(X)$rank 不完全相同。很难找出它何时以及为什么会给出不同的结果。由于你和Joshua的函数速度差别不大，我就拿他的函数。但我赞成你的回答。

@Roland，你是对的，我刚刚比较了我的函数和qr。 1e-7 是这里的问题：对于等级 0，我会说它应该是 == 0，然后等级 1 存在更多问题，因为即使所有条目都大约是 1e-300，qr 也会输出 2，这是正确的。但是这些条目的乘积在 R 中是 0，myrank 返回 1，所以这不再是一个有效的解决方案。划分行可能会起作用，但随后功能会变慢。

【参考方案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 是 最大的支点。