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
  • crossprodも%*%も直接BLASを呼び出しますがtはごにょごにょしますんでその差かと. -- なかま 2004-10-07 (木) 11:06:32
  • 中間さんのコメント参考にもう一つ実験。確かに 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には -- なかま 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のチェック処理は行われません。 つまるところ安全に計算させるのと,チューニングされるであろう関数の存在と言うところでしょうか。
  • いつも的確なコメントありがとうございます -> 中間さん。ソースを読まにゃあかんな。 -- 2004-10-08 (金) 07:39:44


トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Google
WWW を検索 OKADAJP.ORG を検索
Last-modified: 2015-03-01 (日) 01:15:59 (1720d)