![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
Rコードの最適化例：行列のクロス積 (Rコード最適化のコツと実例集 に戻る)
行列の積は順序が大事という例 (r-help 記事より、2004.10.07)
失礼：順序の問題ではなく、クロス積を計算する二つの同値な関数 crossprod と %*% の違いというべきでした。実際、?crossprod によれば crossprod 関数の方が早いと書いてありました。どうも outer , %o% と違い単なる別名になっていないようです。
二次形式の計算を様々に行ない比較する。結構計算速度が違うことがわかる。
> f1 <- function(a,X){ ignore <- t(a) %*% X %*% a }
> f2 <- function(a,X){ ignore <- crossprod(t(crossprod(a,X)),a) }
> f3 <- function(a,X){ ignore <- crossprod(a,X) %*% a }
> f4 <- function(a,X){ ignore <- (t(a) %*% X) %*% a }
> f5 <- function(a,X){ ignore <- t(a) %*% (X %*% a) }
> f6 <- function(a,X){ ignore <- crossprod(a,crossprod(X,a)) }
> a <- rnorm(100); X <- matrix(rnorm(10000),100,100)
> print(system.time( for(i in 1:10000){ a1<-f1(a,X)}))
[1] 3.98 0.02 4.09 0.00 0.00
> print(system.time( for(i in 1:10000){ a2<-f2(a,X)}))
[1] 2.40 0.00 2.45 0.00 0.00
> print(system.time( for(i in 1:10000){ a3<-f3(a,X)}))
[1] 1.31 0.00 1.34 0.00 0.00
> print(system.time( for(i in 1:10000){ a4<-f4(a,X)}))
[1] 4.01 0.00 4.01 0.00 0.00
> print(system.time( for(i in 1:10000){ a5<-f5(a,X)}))
[1] 4.06 0.00 4.09 0.00 0.00
> print(system.time( for(i in 1:10000){ a6<-f6(a,X)}))
[1] 1.19 0.00 1.19 0.00 0.00
> c(a1,a2,a3,a4,a5,a6) # 計算結果のチェック
[1] -56.73955 -56.73955 -56.73955 -56.73955 -56.73955 -56.73955
b <- t(a)
f0 <- function(a,b,X){ ignore <- b %*% X %*% a }
> print(system.time( for(i in 1:10000){ a1<-f0(a,b,X)}))
[1] 3.01 0.05 3.05 0.00 0.00/* Don't trust the BLAS to handle NA/NaNs correctly: PR#4582
* The test is only O(n) here
*/
for (i = 0; i < nrx*ncx; i++)
if (ISNAN(x[i])) {have_na = TRUE; break;}
if (!have_na)
for (i = 0; i < nry*ncy; i++)
if (ISNAN(y[i])) {have_na = TRUE; break;}
などのチェック処理が存在するため,%*%の場合はdo_matprod->matprodと渡り、このチェック処理の為のオーバヘッドが入ります。
方やcrossprodの方はdo_matprod->crossprodと渡り、NANのチェック処理は行われません。
つまるところ安全に計算させるのと,チューニングされるであろう関数の存在と言うところでしょうか。
