矩陣乘法的特例

我想在R中乘以矩陣，但使用了應用函數。在這種特殊情況下，我正在尋找處理來港定居，對此我什麼也沒看見在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中有一些NAs拋出，mat2就像老式的打卡一樣，跟蹤mat1中哪些元素保留在結果中（所以它不是完成 multiplic真正意義上的 - 事實上，打卡是我所追求的，乘法似乎是一種獲得它的方式）。使用％*％，

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

與NAs在各地。首次嘗試對付它們：

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不小的往往還是運行得更快。這種「快而骯髒」的種類贏得了數量級的一天，只使用基地。感謝您的建議（三個）

來源

2016-03-03 JasonD

你做基準未必真的是一個公平的比較。對於meth4，函數可能只需要一行「m1％*％m2」。如果使用Matrix（）在第一個實例中創建它們，而不是matrix（），則不需要將矩陣對象轉換爲Matrix。類似地在最後轉換回矩陣幾乎肯定是不必要的。如果你所說的是要處理非常大的矩陣運算，那麼稀疏矩陣可以節省大量的內存。 – dww

好，夠公平的。對包並不熟悉，因此希望對代碼的其他代碼造成最小的破壞，並考慮data.table如何處理data.frame處理的有趣事情。但將更密切地看待包裹 – JasonD

一個快速而骯髒的解決方案是一個suffciently高值來代替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

來源

2016-03-03 21:48:41

哇，謝謝。這非常深入 - 我需要進一步研究你編寫的這個Rcpp，但即使是第一個似乎它應該是一個更快的解決方案，避免了多次應用和中間的ifelse檢查，每次迭代都有我的方法 – JasonD

最簡單的方法可能是使用稀疏矩陣。

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

如果你願意，你可以用0代替NA：

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

來源

2016-03-04 01:06:36 dww

感謝您的提示。不知道這個包或稀疏矩陣。將檢查時間。 – JasonD

矩陣乘法的特例

回答

相關問題