ANOVA君/高速化の試み - 井関龍太のページ

ホーム   編集 凍結 差分 バックアップ 添付 コピー 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS

ANOVA君/高速化の試み

Last-modified: 2013年09月28日 (土) 22:29:55 (1875d)
Top > ANOVA君 > 高速化の試み

平方和を計算する関数の高速化

ANOVA君の動作を速めるため,平方和の計算を担当する関数(ss.calc)の高速化を試みました。
anovakun 2.0.0までは,ss.calcでは,Rのlm関数とdrop1関数を使ってタイプⅡ平方和を計算していました。
しかし,この方法は,平方和の数だけ重回帰分析を行い,さらにその結果に修正を加える計算を行うことに等しく,かなりの計算時間を要します。
そこで,新バージョンでは,分析ごとに1つの計画行列を作成し,その部分行列に最小二乗法を繰り返すことによって計算を行うようにアルゴリズムを変更しました(具体的な計算手続きについては,Milliken & Johnson, 1992;高橋他, 1989が参考になりました)。
lm関数のように様々な指標を算出したり事後分析に対応したりする代わりに,平方和の計算のみに特化することで計算量を減らすことを試みました。

その結果を検討するために,前のバージョン(anovakun 2.0.0)と新バージョン(anovakun 2.1.0, anovakun 2.2.0)の計算速度を比べてみました。
さらに,anovakun 3.0.0では,デフォルトとしてタイプⅢ平方和を計算するようになりました。
以下のanovakun 3.0.0による結果は,すべてタイプⅢ平方和の計算時間です(anovakun 3.0.0でタイプⅡオプションを選択した場合の計算手順は,anovakun 2.2.0と同じになります)。
anovakun 3.1.0では,aov関数の手続きを利用して,誤差平方和の分解の部分を高速化しました。
この変更は,被験者内要因を含む分析にのみ影響します。
そのため,anovakun 3.1.0の結果は,被験者内要因を含むデータについてのみ示しています。

なお,以下の計算例の大部分は,ss.calcの計算速度を示したもので,ANOVA君の計算全体に要する時間ではありません。
また,実際の計算時間は,使用するマシンのスペック,実際のデータ等によって大きく変動するので,あくまで目安としてご覧ください。

被験者間計画の場合

Rのsystem.time関数を使って計算時間を調べます。

> system.time(計算時間を調べたい関数)

この関数を使うと,()内の関数の計算結果は表示されなくなる代わりに,計算に要した時間を返してくれます。
この関数の中でANOVA君のss.calc関数を使ってみます。

#anovakun 2.0.0
#3×2のABsデザイン(n = 10)
> system.time(ss.calc(full.elem, dat))
   ユーザ   システム       経過 
      0.09       0.00       0.09

#anovakun 2.1.0
#3×2のABsデザイン(n = 10)
> system.time(ss.calc(full.elem, dat))
   ユーザ   システム       経過 
      0.02       0.00       0.02 

#anovakun 2.2.0
#3×2のABsデザイン(n = 10)
> system.time(ss.calc(full.elem, dat))
   ユーザ   システム       経過 
      0.01       0.00       0.01 

#anovakun 3.0.0
#3×2のABsデザイン(n = 10)
> system.time(ss.calc(full.elem, dat))
   ユーザ  システム      経過 
      0.01       0.00       0.01 

単位は秒なので,それぞれ,0.09秒,0.02秒,0.01秒,0.01秒の計算時間ということになります。

ちなみに,以下のデータもそうですが,主に計算例用のデータなので,被験者数が少なかったり,水準数が多かったりします。
被験者数が少ないと差は現れにくいですが,それでも新バージョンの方が速くなっています。

混合要因計画の場合1

混合要因計画の場合も見てみます。

#anovakun 2.0.0
#2×3×2のABsCデザイン(n = 15)
> system.time(ss.calc(full.elem, dat))
   ユーザ   システム       経過 
      0.37       0.00       0.37 

#anovakun 2.1.0
#2×3×2のABsCデザイン(n = 15)
> system.time(ss.calc(full.elem, dat))
   ユーザ   システム       経過 
      0.07       0.00       0.07 

#anovakun 2.2.0
#2×3×2のABsCデザイン(n = 15)
> system.time(ss.calc(full.elem, dat))
   ユーザ   システム       経過 
      0.04       0.00       0.04 

#anovakun 3.0.0
#2×3×2のABsCデザイン(n = 15)
> system.time(ss.calc(full.elem, dat))
   ユーザ  システム      経過 
      0.03       0.00       0.03 

#anovakun 3.1.0
#2×3×2のABsCデザイン(n = 15)
> system.time(ss.calc(full.elem, dat))
   ユーザ  システム      経過 
      0.02       0.00       0.02 

この場合には,もう少し差が顕著になっていると思います。

混合要因計画の場合2

上の例よりもさらにデータ数の多いケースです。

#anovakun 2.0.0
#2×3×3のAsBCデザイン(n = 22)
> system.time(ss.calc(full.elem, dat))
   ユーザ   システム       経過  
      1.77       0.00       1.77 

#anovakun 2.1.0
#2×3×3のAsBCデザイン(n = 22)
> system.time(ss.calc(full.elem, dat))
   ユーザ   システム       経過 
      0.68       0.00       0.69 

#anovakun 2.2.0
#2×3×3のAsBCデザイン(n = 22)
> system.time(ss.calc(full.elem, dat))
   ユーザ   システム       経過 
      0.38       0.00       0.38 

#anovakun 3.0.0
#2×3×3のAsBCデザイン(n = 22)
> system.time(ss.calc(full.elem, dat))
   ユーザ  システム      経過 
      0.31       0.00       0.31 

#anovakun 3.1.0
#2×3×3のAsBCデザイン(n = 22)
> system.time(ss.calc(full.elem, dat))
   ユーザ  システム      経過 
      0.13       0.00       0.13 

違いがさらに顕著になっていると思います。

被験者内計画の場合

それでは,この新たなss.calc関数を実際にANOVA君に組み込んだ場合の違いを見てみます。
以下のデータは,一次の交互作用が多く,多重比較も必要になるために,平方和の計算を繰り返し必要とするケースです。

#anovakun 2.0.0
#下位検定の必要性が多いために時間のかかるケース(n = 12)
> system.time(anovakun(dat, "sABCD", 2, 2, 2, 3))
   ユーザ   システム       経過  
      35.78       0.53       36.60

#anovakun 2.1.0
#下位検定の必要性が多いために時間のかかるケース(n = 12)
> system.time(anovakun(dat, "sABCD", 2, 2, 2, 3))
   ユーザ   システム       経過 
      10.19       0.00       10.21 

#anovakun 2.2.0
#下位検定の必要性が多いために時間のかかるケース(n = 12)
> system.time(anovakun(dat, "sABCD", 2, 2, 2, 3))
   ユーザ   システム       経過 
       9.65       0.01       9.79 

#anovakun 3.0.0
#下位検定の必要性が多いために時間のかかるケース(n = 12)
> system.time(anovakun(dat, "sABCD", 2, 2, 2, 3))
   ユーザ  システム      経過 
       6.15       0.02       6.19 

#anovakun 3.1.0
#下位検定の必要性が多いために時間のかかるケース(n = 12)
> system.time(anovakun(dat, "sABCD", 2, 2, 2, 3))
   ユーザ  システム      経過 
       3.66       0.01       3.67 

このケースでは,もともと長い時間がかかっていたこともあり,かなり計算時間が短縮されました。
anovakun 3.0.0では,それまでのバージョンとは平方和の計算法が異なるため(タイプⅢ平方和),さらに時間が短くなっています(釣り合い型計画では,タイプⅡの場合と同じ結果になります)。
このスピードがanovakun 3.0.0でデフォルトの平方和を決定する際の決め手になりました。
また,誤差平方和の計算速度を速めたので,anovakun 3.1.0ではさらにスピードアップしています。

混合要因計画の場合3

しかし,新バージョンの計算法にも少々の工夫が必要なようです。

#anovakun 2.0.0
#3×4のAsBデザイン(n = 39)
> system.time(ss.calc(full.elem, dat))
   ユーザ   システム       経過 
      0.52       0.00       0.52 

#anovakun 2.1.0
#3×4のAsBデザイン(n = 39)
> system.time(ss.calc(full.elem, dat))
   ユーザ   システム       経過 
      1.26       0.00       1.26 

#anovakun 2.2.0
#3×4のAsBデザイン(n = 39)
> system.time(ss.calc(full.elem, dat))
   ユーザ   システム       経過 
      0.25       0.00       0.25 

#anovakun 3.0.0
#3×4のAsBデザイン(n = 39)
> system.time(ss.calc(full.elem, dat))
   ユーザ  システム      経過 
      0.22       0.01       0.23 

#anovakun 3.1.0
#3×4のAsBデザイン(n = 39)
> system.time(ss.calc(full.elem, dat))
   ユーザ  システム      経過 
      0.09       0.00       0.09 

これは,上の2つと同様に混合要因計画のデータですが,anovakun 2.1.0の方が計算に長くかかっています。
しかも,2要因の計画ですが,上の3要因の例よりも時間がかかっています。
どうやら要因の数が少なくても,水準数やデータ数が多いとかえって長くかかる場合があるようです。
いろいろ試してみたところ,混合要因計画で水準数が多いときに逆転しやすいようです。
おそらく,これは,このような場合に計画行列のサイズが大きくなりやすいことが原因のようです(行列演算は,サイズが大きいほど長くかかります)。

そこで,anovakun 2.2.0では,最小二乗解をより速く得られるオプションを使用し(LAPACK),さらに,最小二乗解の得られない部分を予め省略しておくことによって計算を速めました。
このおかげで,anovakun 2.2.0では,このケースでもanovakun 2.0.0よりも高速化することができました。
ただし,省略は混合要因計画でのみ有効であるため,被験者間計画や被験者内計画では,anovakun 2.1.0とanovakun 2.2.0の間にそれほどの差は出ないはずです。

anovakun 3.1.0では,これとは別な方法で誤差平方和を計算しているため,省略の必要はなくなりました。

おまけその1

上の被験者内計画のデータを別のPCで分析したものです。

#anovakun 2.0.0
#下位検定の必要性が多いために時間のかかるケース(n = 12)
> system.time(anovakun(dat, "sABCD", 2, 2, 2, 3))
   ユーザ   システム       経過  
     24.82       0.40      25.37 

#anovakun 2.1.0
#下位検定の必要性が多いために時間のかかるケース(n = 12)
> system.time(anovakun(dat, "sABCD", 2, 2, 2, 3))
   ユーザ   システム       経過  
      6.85       0.07       6.98 

#anovakun 2.2.0
#下位検定の必要性が多いために時間のかかるケース(n = 12)
> system.time(anovakun(dat, "sABCD", 2, 2, 2, 3))
   ユーザ   システム       経過  
      6.23       0.19       6.42  

#anovakun 3.0.0
#下位検定の必要性が多いために時間のかかるケース(n = 12)
> system.time(anovakun(dat, "sABCD", 2, 2, 2, 3))
   ユーザ   システム       経過  
      4.00       0.02       4.03 

#anovakun 3.1.0
#下位検定の必要性が多いために時間のかかるケース(n = 12)
> system.time(anovakun(dat, "sABCD", 2, 2, 2, 3))
   ユーザ   システム       経過  
      2.41       0.00       2.41 

マシンだけでなくRのバージョンも違うので(今回はR 2.6.0,他はすべてR 2.5.1)フェアな比較とはいえませんが,環境が違うとこのくらい違うことがあるという例としてご覧ください。

おまけその2

anovakun 3.1.0における誤差平方和の計算の高速化は,もちろん,タイプⅡ平方和の計算にも影響します。
以下の例は,上でも扱った被験者内計画のデータなので,計算結果はどちらも同じになります。

#anovakun 3.1.0;タイプⅢ平方和
#下位検定の必要性が多いために時間のかかるケース(n = 12)
> system.time(anovakun(dat, "sABCD", 2, 2, 2, 3))
   ユーザ  システム      経過 
       3.66      0.01      3.67 

#anovakun 3.1.0;タイプⅡ平方和
#下位検定の必要性が多いために時間のかかるケース(n = 12)
> system.time(anovakun(dat, "sABCD", 2, 2, 2, 3, type2 = T))
   ユーザ  システム      経過 
       3.87      0.00      3.91 

もともと,anovakun 3.1.0の更新は,以下のケースで非常に長い時間がかかってしまうことを解消するために行ったものでした。

#anovakun 3.0.0
#デザインが複雑でデータ数が多いため時間のかかるケース(n = 120)
> system.time(anovakun(dat, "ABsCD", 2, 2, 2, 2))
   ユーザ  システム      経過 
      33.41      0.60     34.34 

#anovakun 3.1.0
#デザインが複雑でデータ数が多いため時間のかかるケース(n = 120)
> system.time(anovakun(dat, "ABsCD", 2, 2, 2, 2))
   ユーザ  システム      経過 
       9.27      0.13      9.40 

もう少し速いとありがたいですが,とりあえずは,計算時間を短縮できたということで。

文献

Milliken, G. A., & Johnson, D. E. (1992). Analysis of messy data volume 1:Designed experiments. Chapman and Hall:New York.
高橋行雄・大橋靖雄・芳賀敏郎 (1989). SASによる実験データの解析. 東京大学出版会

TrackBack(0) | 外部リンク元 | このエントリーをはてなブックマークに追加