外積は 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() を使う必要はなかった。