矩阵乘法的特例

Question

我正在尝试将 R 中的矩阵相乘，但使用的是应用函数。在这种特殊情况下，我正在寻找处理 NA，为此我在 crossprod 中没有看到任何要处理的内容，或者使用 %*%

set.seed(3141)
mat1 <- c(1:50)
pos <- sample(c(1:50),14)
mat1[pos] <- NA
mat1 <- matrix(mat1,10,5)
mat2 <- matrix(sample(c(0,1),20,replace=T),5,4)

mat1:

       [,1] [,2] [,3] [,4] [,5]
  [1,]    1   11   NA   31   41
  [2,]   NA   12   NA   32   NA
  [3,]   NA   13   NA   NA   NA
  [4,]    4   14   24   34   44
  [5,]    5   15   25   NA   45
  [6,]    6   16   26   36   46
  [7,]    7   17   27   37   47
  [8,]    8   18   28   NA   NA
  [9,]    9   19   29   NA   49
 [10,]   10   20   NA   40   NA

mat2:

      [,1] [,2] [,3] [,4] 
 [1,]    0    0    0    1 
 [2,]    1    0    1    1 
 [3,]    0    1    0    0 
 [4,]    0    1    1    0 
 [5,]    1    1    1    1

所以 mat1 有一些 NA 被扔进去，而 mat2 就像一张旧的打孔卡，跟踪 mat1 的哪些元素保留在结果中（所以它不是真正意义上的完整乘法 - 打孔卡是真的是我所追求的，乘法似乎是一种获得它的方法）。使用 %*%，

mat3 <- mat1 %*% mat2

      [,1] [,2] [,3] [,4]
 [1,]   NA   NA   NA   NA
 [2,]   NA   NA   NA   NA
 [3,]   NA   NA   NA   NA
 [4,]   58  102   92   62
 [5,]   NA   NA   NA   NA
 [6,]   62  108   98   68
 [7,]   64  111  101   71
 [8,]   NA   NA   NA   NA
 [9,]   NA   NA   NA   NA
[10,]   NA   NA   NA   NA

到处都有 NA。第一次尝试处理它们：

    mat4 <- t(apply(mat1,1,function(x){apply(mat2,2,function(y){sum(x*y,na.rm=T)})}))

      [,1] [,2] [,3] [,4]
 [1,]   52   72   83   53
 [2,]   12   32   44   12
 [3,]   13    0   13   13
 [4,]   58  102   92   62
 [5,]   60   70   60   65
 [6,]   62  108   98   68
 [7,]   64  111  101   71
 [8,]   18   28   18   26
 [9,]   68   78   68   77
[10,]   20   40   60   30

这更好，但挑剔的复杂情况是我想删除任何试图从 mat1 中包含 NA 的结果，因此它不会对最终结果产生影响。

mat5 <- t(apply(mat1,1,function(x){
  apply(mat2,2,function(y){
    ifelse(is.na(sum(x[as.logical(y)])),
           0,
           sum(x*y,na.rm=T))
  })}))

      [,1] [,2] [,3] [,4]
 [1,]   52    0   83   53
 [2,]    0    0    0    0
 [3,]    0    0    0    0
 [4,]   58  102   92   62
 [5,]   60    0    0   65
 [6,]   62  108   98   68
 [7,]   64  111  101   71
 [8,]    0    0    0    0
 [9,]   68    0    0   77
[10,]    0    0    0    0

这就是我要去的地方，因为我只在 mat1 有 NA 时才抛出结果（即 mat2 有相应的 1，但如果没有，那么 NA 很好）。

问题是，这是一个有效的解决方案吗？我是否错过了基地中的某些东西，这会使这一切变得更快？ （缺乏并行化，因为我很遗憾在 Windows 上这样的事情不适合胆小的人）。这看起来很笨重，并且必须在多个阵列上执行数百万次，因此任何加速都是有用的。谢谢。

更新： 感谢您到目前为止的两个回复。我想我会在我的机器上进行时间比较，看看这些方法可能有何不同。不幸的是，我无法让 C++ 工作。我收到一条错误消息，指出构建共享库时出错。它建议我从 CRAN 下载兼容版本的 Rtools（我正在使用 R3.2.3），但我也在考虑这必须在需要额外安装等的其他计算机（比如我老板的）上运行使这项工作可能并不理想。包，我可以写入代码，但是如果代码抛出错误来修复它，访问一个站点以下载一些不属于标准安装的附加内容，这有点复杂。无论如何，对于其他人：

meth1 <- function(m1,m2){
  t(apply(m1,1,function(x){
    apply(m2,2,function(y){
      ifelse(is.na(sum(x[as.logical(y)])),
             0,
             sum(x*y,na.rm=T))
    })}))
}
meth2 <- function(m1,m2){
  m1[is.na(m1)] <- 10^20
  res <- m1 %*% m2
  res[abs(res) > 10^10] <- 0
  res
}

library(Matrix)
meth4 <- function(m1,m2){
  M1 <- Matrix(m1,sparse=TRUE)
  M2 <- Matrix(m2,sparse=TRUE)
  res <- M1 %*% M2
  res[is.na(res)] <- 0
  Matrix(res,sparse = F)
}

library(microbenchmark)
microbenchmark({meth1(mat1,mat2)},{meth2(mat1,mat2)},{meth4(mat1,mat2)},times=100)

屈服：

Unit: microseconds
                      expr      min       lq       mean   median       uq
 {     meth1(mat1, mat2) }  475.957  516.155  563.41297  535.826  568.754
 {     meth2(mat1, mat2) }    8.126    9.836   14.78396   15.609   18.816
 {     meth4(mat1, mat2) } 4535.489 4764.701 5016.47097 4901.331 5008.025
      max neval
 1763.565   100
   30.791   100
 9722.265   100

对 Rcpp 感到羞耻 - 我很欣赏它看起来付出了不小的努力，而且 C 中的东西往往运行得更快。这种“快速而肮脏”的方式赢得了数量级的胜利，并且只使用了基础。感谢您的建议（所有三个）

Answer 1

一个快速但肮脏的解决方案是将NA 替换为足够高的值，然后使用阈值来挑选零：

mat1[is.na(mat1)] <- 10^200
A <- mat1 %*% mat2
A[abs(A) > 10^100] <- 0
A
      [,1] [,2] [,3] [,4]
 [1,]   52    0   83   53
 [2,]    0    0    0    0
 [3,]    0    0    0    0
 [4,]   58  102   92   62
 [5,]   60    0    0   65
 [6,]   62  108   98   68
 [7,]   64  111  101   71
 [8,]    0    0    0    0
 [9,]   68    0    0   77
[10,]    0    0    0    0

或者您可以简单地使用 Rcpp 编写自己的代码：

library(inline)
library(Rcpp)
cppFunction(
    'NumericMatrix f(NumericMatrix mat1, NumericMatrix mat2) {
        double val;
        NumericMatrix X(mat1.nrow(), mat2.ncol());
        for (int i = 0; i < mat1.nrow(); ++i) {
            for (int j = 0; j < mat1.ncol(); ++j) {
                val = 0;
                for(int k = 0; k < mat1.ncol(); k++){
                    if(NumericVector::is_na(mat1(i, k))){
                        if( mat2(k, j) != 0) {
                            val = 0;
                            break;
                        }
                    } else val += mat1(i, k)*mat2(k, j);
                }
                X(i, j) = val;
            }
        }
        return X;
    }'
)

> f(mat1, mat2)
      [,1] [,2] [,3] [,4]
 [1,]   52    0   83   53
 [2,]    0    0    0    0
 [3,]    0    0    0    0
 [4,]   58  102   92   62
 [5,]   60    0    0   65
 [6,]   62  108   98   68
 [7,]   64  111  101   71
 [8,]    0    0    0    0
 [9,]   68    0    0   77
[10,]    0    0    0    0

Answer 2

最简单的方法可能是使用稀疏矩阵。

library(Matrix)
M1 <- Matrix(mat1,sparse=TRUE)
M2 <- Matrix(mat2,sparse=TRUE)
ans <- M1 %*% M2
ans
10 x 4 sparse Matrix of class "dgCMatrix"

 [1,] 52  NA  83 53
 [2,] NA  NA  NA NA
 [3,] NA  NA  NA NA
 [4,] 58 102  92 62
 [5,] 60  NA  NA 65
 [6,] 62 108  98 68
 [7,] 64 111 101 71
 [8,] NA  NA  NA NA
 [9,] 68  NA  NA 77
[10,] NA  NA  NA NA

如果您愿意，可以将 NA 替换为 0：

ans[is.na(ans)] <- 0
Matrix(ans,sparse = F)

 10 x 4 Matrix of class "dgeMatrix"
      [,1] [,2] [,3] [,4]
 [1,]   52    0   83   53
 [2,]    0    0    0    0
 [3,]    0    0    0    0
 [4,]   58  102   92   62
 [5,]   60    0    0   65
 [6,]   62  108   98   68
 [7,]   64  111  101   71
 [8,]    0    0    0    0
 [9,]   68    0    0   77
[10,]    0    0    0    0