Julia でベクトル化プログラミングは有効か?

モンテカルロ法による円周率を求めるプログラムは，練習問題として言語を問わずよく取り上げられる。

言語によっては以下のようなプログラムを「まずいプログラム」とみなす場合があるが，Julia の場合は，実行効率が非常に高い。

function estimate_pi(n::Int)
    count = 0
    for i = 1:n
        x, y = rand(), rand()
        if x^2 + y^2 <= 1
            count += 1
        end
    end
    return 4 * count / n
end

estimate_pi (generic function with 1 method)

@time pi_est = estimate_pi(100000000)
println("Estimated pi: ", pi_est)

  0.254032 seconds (1 allocation: 16 bytes)
Estimated pi: 3.14138308

実行所要時間は平均的に 0.25 秒ほどである。

プログラマによっては以下のようにベクトル化するのを好む場合がある。しかし，Julia では上のプログラムに比べて，遅い!！

function estimate_pi2(n::Int)
    xs, ys = rand(n), rand(n)
    count0 = count(xs.^2 + ys.^2 .<= 1)
    return 4 * count0 / n
end

estimate_pi2 (generic function with 1 method)

@time pi_est = estimate_pi2(100000000)
println("Estimated pi: ", pi_est)

  0.666811 seconds (15 allocations: 3.737 GiB, 6.35% gc time)
Estimated pi: 3.14170676

実行所要時間はばらつきがあるがたいていは 1.6 秒ほどである。ときどき 0.6 秒ほどのときもあるが，そういうのは少ない。

ちなみに，ベクトル化といえば R の十八番であるが，R でも 1.4 秒台である。

system.time({n = 100000000;print(mean(runif(n)^2+runif(n)^2 < 1)*4)})

[1] 3.141476



   ユーザ   システム       経過  
     1.441      0.166      1.606 

考察するに，Julia の場合，ベクトル化を意識して書かれたプログラムは面倒くさい処理をするために複雑なコードを生成しているが，for ループで単純に書かれたプログラムは簡潔なコードを生成しているのであろうか。

一概にはいえないだろうが，まずは単純なプログラムを書いて，効果を見ながらブラッシュアップしていくというのがよい戦略ではないだろうか。