SIZE(20){COLOR(magenta){Rコードの最適化例:行列のクロス積}} ([[Rコード最適化のコツと実例集]] に戻る)~

行列の積は順序が大事という例 (r-help 記事より、2004.10.07)~

失礼:順序の問題ではなく、クロス積を計算する二つの同値な関数 crossprod と 
%*% の違いというべきでした。実際、?crossprod によれば crossprod 関数の方が
早いと書いてありました。どうも outer , %o% と違い単なる別名になっていないようです。

失礼:順序の問題ではなく、クロス積を計算する二つの同値な関数 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
-crossprodも%*%も直接BLASを呼び出しますがtはごにょごにょしますんでその差かと. -- [[なかま]] &new{2004-10-07 (木) 11:06:32};
-中間さんのコメント参考にもう一つ実験。確かに COLOR(red){t(a)} を計算するオーバーヘッドがあるようですが、それだけでは説明できないような。
 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
-よくよくみるとsrc/main/array.cには -- [[なかま]] &new{2004-10-07 (木) 22:08:41};
 /* 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のチェック処理は行われません.
つまるところ安全に計算させるのと,チューニングされるであろう関数の存在
と言うところでしょうか。
などのチェック処理が存在するため,%*%の場合はdo_matprod->matprodと渡り、このチェック処理の為のオーバヘッドが入ります。
方やcrossprodの方はdo_matprod->crossprodと渡り、NANのチェック処理は行われません。
つまるところ安全に計算させるのと,チューニングされるであろう関数の存在と言うところでしょうか。
-いつも的確なコメントありがとうございます -> 中間さん。ソースを読まにゃあかんな。 --  &new{2004-10-08 (金) 07:39:44};

#comment

トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS