2010-09-04 185 views
56

我正在優化一個函數,我想擺脫慢循環。 我正在尋找一種更快的方式來乘以一個矢量矩陣的每一行。用矢量乘以矩陣的行嗎?

任何想法?

編輯:

我不是在尋找一個'古典'乘法。

例如,我有矩陣,有23列25行和長度爲23的矢量。 結果我想要矩陣25x23,每行乘以矢量。

回答

58

我認爲你正在尋找sweep()

> (mat <- matrix(rep(1:3,each=5),nrow=3,ncol=5,byrow=TRUE)) 
    [,1] [,2] [,3] [,4] [,5] 
[1,] 1 1 1 1 1 
[2,] 2 2 2 2 2 
[3,] 3 3 3 3 3 
> vec <- 1:5 
> sweep(mat,MARGIN=2,vec,`*`) 
    [,1] [,2] [,3] [,4] [,5] 
[1,] 1 2 3 4 5 
[2,] 2 4 6 8 10 
[3,] 3 6 9 12 15 

這是R的核心功能之一,雖然多年來已經有所改進。

31
> MyMatrix <- matrix(c(1,2,3, 11,12,13), nrow = 2, ncol=3, byrow=TRUE) 
> MyMatrix 
    [,1] [,2] [,3] 
[1,] 1 2 3 
[2,] 11 12 13 
> MyVector <- c(1:3) 
> MyVector 
[1] 1 2 3 

您可以使用兩種:

> t(t(MyMatrix) * MyVector) 
    [,1] [,2] [,3] 
[1,] 1 4 9 
[2,] 11 24 39 

或:

> MyMatrix %*% diag(MyVector) 
    [,1] [,2] [,3] 
[1,] 1 4 9 
[2,] 11 24 39 
-2

谷歌 「R矩陣multiplcation」 產量Matrix Multiplication,描述%*%操作,並表示「將兩個矩陣相乘,如果它們是一致的,如果一個參數是一個向量,它將被提升爲一個行或列矩陣來使這兩個參數一致,如果兩個都是向量,它將返回內部的p產品(作爲矩陣)「。

+2

問題不 – MHH 2013-12-23 03:34:36

21

其實,sweep是不是我的電腦上最快的選項:

MyMatrix <- matrix(c(1:1e6), ncol=1e4, byrow=TRUE) 
MyVector <- c(1:1e4) 

Rprof(tmp <- tempfile(),interval = 0.001) 
t(t(MyMatrix) * MyVector) # first option 
Rprof() 
MyTimerTranspose=summaryRprof(tmp)$sampling.time 
unlink(tmp) 

Rprof(tmp <- tempfile(),interval = 0.001) 
MyMatrix %*% diag(MyVector) # second option 
Rprof() 
MyTimerDiag=summaryRprof(tmp)$sampling.time 
unlink(tmp) 

Rprof(tmp <- tempfile(),interval = 0.001) 
sweep(MyMatrix ,MARGIN=2,MyVector,`*`) # third option 
Rprof() 
MyTimerSweep=summaryRprof(tmp)$sampling.time 
unlink(tmp) 

Rprof(tmp <- tempfile(),interval = 0.001) 
t(t(MyMatrix) * MyVector) # first option again, to check order 
Rprof() 
MyTimerTransposeAgain=summaryRprof(tmp)$sampling.time 
unlink(tmp) 

MyTimerTranspose 
MyTimerDiag 
MyTimerSweep 
MyTimerTransposeAgain 

這產生了:

> MyTimerTranspose 
[1] 0.04 
> MyTimerDiag 
[1] 40.722 
> MyTimerSweep 
[1] 33.774 
> MyTimerTransposeAgain 
[1] 0.043 

論是最慢的選擇上面,第二個選項達到內存限制(2046 MB)。然而,考慮到其餘的選擇,在我看來,雙轉換似乎比sweep好很多。


編輯

只是想較小的數據的次重複的次數:「你怎麼一個向量乘以矩陣」

MyMatrix <- matrix(c(1:1e3), ncol=1e1, byrow=TRUE) 
MyVector <- c(1:1e1) 
n=100000 

[...] 

for(i in 1:n){ 
# your option 
} 

[...] 

> MyTimerTranspose 
[1] 5.383 
> MyTimerDiag 
[1] 6.404 
> MyTimerSweep 
[1] 12.843 
> MyTimerTransposeAgain 
[1] 5.428 
+3

根據我的經驗,如果你向矩陣中投入一堆「NA」,「diag」所花費的時間似乎要經過屋頂。對於包含1E5'NA's的1E4x1E4墊,我得到:'MyTimerTranspose' = 0.014,'MyTimerSweep' = 0.042,'MyTimerDiag' = 67.738。我會複製,但我不耐煩...只是要記住。 – jbaums 2012-03-01 22:19:51

+0

我真的很喜歡雙轉換的答案,主要是因爲它顯示瞭如果我們用「列」替換「行」,答案是簡單的A * x,這是不明顯的,除非你真正理解R如何工作矩陣。 – MHH 2013-12-23 03:38:27