矩阵乘法的特例
Posted
技术标签:
【中文标题】矩阵乘法的特例【英文标题】：special case of matrix multiplication 【发布时间】：2016-03-03 19:41:00 【问题描述】：我正在尝试将 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 中的东西往往运行得更快。这种“快速而肮脏”的方式赢得了数量级的胜利，并且只使用了基础。感谢您的建议（所有三个）
【问题讨论】：
您所做的基准测试可能并不是一个公平的比较。对于 meth4，该函数可能只需要一行“m1 %*% m2”。如果您在第一个实例中使用 Matrix() 而不是 matrix() 创建矩阵对象，则不需要将矩阵对象转换为 Matrix。最后类似地转换回矩阵几乎肯定是不必要的。如果您如您所说要使用非常大的矩阵运算，那么稀疏矩阵可以节省大量内存。 好的，很公平。不熟悉该包，因此希望对其余代码的干扰最小，例如考虑 data.table 如何对 data.frame 处理做有趣的事情的经验。但会更仔细地查看包裹 【参考方案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
【讨论】：
哇，谢谢。那里很深 - 我需要进一步研究您编写的这个 Rcpp，但即使第一个似乎也应该是一个更快的解决方案，在我进行的每次迭代中避免多次应用和中间的 ifelse 检查我的方法【参考方案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
【讨论】：
感谢您的提示。不知道这个包或稀疏矩阵。将检查时间。以上是关于矩阵乘法的特例的主要内容，如果未能解决你的问题，请参考以下文章