裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

外積の和

2020年09月02日 | ブログラミング

外積は R では outer(a, b, 演算子) で計算できる。

c(1, 2, 3, 4, 5) と c(2, 3, 4) で,演算子を "*" (積)とすると,以下のようになる。

> x = outer(1:5, 2:4, "*")
> x
     [,1] [,2] [,3]
[1,]    2    3    4
[2,]    4    6    8
[3,]    6    9   12
[4,]    8   12   16
[5,]   10   15   20

結果の行列の各要素の和は sum()  で求められる。

> print(sum(x))
[1] 135

さて,ここで,この結果を得るための最良の方法を探索しよう。

長さが共に n = 1e4 の正規乱数ベクトル a, b の sum(outer(a, b, "*")) を求めることとする。

outer() を使えば以下のようになる。答えは 2159.779 であり,0.309 秒ほどかかった。これでも十分速い。

> system.time({
+   x = outer(a, b, "*")
+   print(sum(x)) # 2159.779
+ })
[1] 2159.779
   ユーザ   システム       経過  
     0.309      0.150      0.459 

outer() がやっていることを R で素直(?)に書けば以下のようになる。

> system.time({
+   sum.of.outer.multiple = 0
+   for (i in seq_along(a)) {
+     for (j in seq_along(b)) {
+       sum.of.outer.multiple = sum.of.outer.multiple + a[i] * b[j]
+     }
+   }
+   print(sum.of.outer.multiple)
+ })
[1] 2159.779
   ユーザ   システム       経過  
     6.773      0.019      6.818

20  倍も遅い

内側の for ループは,a[i] * (b[1] + b[2] + ... + b[n]) を計算している訳だが,ループでは (b[1] + b[2] + ... + b[n]) はいつも同じなので,前もって計算しておいてそれを使えば,ループが省けるし計算量も削減できる。

> system.time({
+   sum.b = sum(b)
+   sum.of.outer.multiple = 0
+   for (a.i in a) {
+     sum.of.outer.multiple = sum.of.outer.multiple + a.i * sum.b
+   }
+   print(sum.of.outer.multiple)
+ })
[1] 2159.779
   ユーザ   システム       経過  
     0.003      0.000      0.003 

いやいや,待って待って。この for ループも (a[1] + a[2] + ... + a[n]) * sum.b を計算していますね。

な〜んだ。sum.a = a[1] + a[2] + ... + a[n] としておいて,sum.a * sum.b で答えが出ますね。

> system.time({
+   print(sum(a)*sum(b))
+ }) # 0.000
[1] 2159.779
   ユーザ   システム       経過  
         0          0          0 

瞬殺でした。outer() を使う必要はなかった。

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村