Rは本当に遅いのか？Juliaとの比較例

[このコード](https://gist.github.com/unaoya/fa7a9bfbf5f71981af8bd3fbc2b0f0f6)について、twitter上で「Juliaなら数十秒で終わるのにRだとめっちゃ時間かかったんだろうな…」的な発言が流れてきたのを見かけて、本当にそうなのか気になったので少し調べました。

## そもそもRのループは遅いのか？

たしかに、Rのforループが非常に遅いとされていた時代はありました。繰返し処理はなるべくベクトル化して書くべきものであって、forを書くというのは可能であれば避けるべき作法でした。

しかし、R 3.4.0からJITコンパイラが同梱されており、これがデフォルトでONになっている恩恵で、現代のRのforループの速度は相当に改善されています。場合によってはforを書いたほうが速度的に有利なことすらあります。

## やってみる

まあともかくやってみましょう。

### 実測とプロファイリング

まず「そもそも元のコードがどれだけ時間がかかるのか」を確認しておくと、私の手元のちょっと古くなってきたMac mini (Late 2012, 2.3 GHz Intel Core i7)だと4分くらいでした。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109348/9f3b5237-645f-8153-257d-7f5428c9d47d.png)


Juliaが数十秒で終わるのであれば（後で確認します）、確かに時間がかかっていると言えそうです。

`sample()`や`ifelse()`が実際には不要であるとか、多少のリファクタリングの余地はあるものの、おそらくこれは速度に対して本質的な部分ではありません。なので、まずはなるべく元のコードを尊重して話を進めることにします。

どこが遅いのかを確認するには、RStudioのProfile機能を使うと便利です。コードの調べたい部分を選択して「Profile > Profile Selected Line(s)」を実行するとプロファイリングが行われます。結果は次のようになります（調査のためループ回数は減らしてあります）。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109348/7e7c2744-191d-dc21-e169-3cca5250e235.png)


結果から明らかなように、`rnorm()`がかなり重たい処理のようです。逆に言えば、ここさえ高速化できれば速度改善が見込めそうです。

### dqrngパッケージを使う

[dqrngパッケージ](https://www.daqana.org/dqrng/)は、高速な疑似乱数の生成を可能にするパッケージです。このパッケージではpcg64、Xoroshiro、Threefryといった乱数生成器がサポートされており、デフォルトではこの中で最速のXoroshiro128+が採用されています。XoroshiroはXorshiftの変種で、統計的性質についてメルセンヌ・ツイスタよりも優れた部分もあるようです（cf. [xoshiro/xoroshiro generators and the PRNG shootout](http://xoshiro.di.unimi.it/)）。

各関数の利用方法はRの標準の乱数生成関数とほぼ同じです。「`rnorm()`に対応する関数は`dqrnorm()`である」といった具合に頭に`dq`が付いたものが`dqrng`版の関数になっています。ただ、現時点では対応している関数はまだ多くありません。

パッケージはCRANからもインストールできますが、今回は`sample()`の`dqrng`版である`dqsample()`も使いたいので、開発版をインストールします。開発版のインストールは`drat`パッケージを使う方法が推奨されているようです。

```r:dqrngの開発版をインストールする
if (!requireNamespace("drat", quietly = TRUE)) install.packages("drat")
drat::addRepo("daqana")
install.packages("dqrng")
```

で、元のコードを`dqrng`パッケージを使うように修正しました。具体的には、次の修正を行っただけです。軽微な変更なのでコードの提示は省略します。

- `library(dqrng)`を追加
- `rnorm`を`dqrnorm`に変更
- `sample`を`dqsample`に変更

結果、実行時間は次のようになりました。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109348/da65687f-c645-441d-3edc-cd3a197a8721.png)

単に`dqrng`パッケージを使うだけで、実行時間を約6割短くすることができました。しかしまだ100秒かかっているので「数十秒」と言うにはチョット辛いところです。

### 少しリファクタリングする

もう少し速度を追ってみたいと思うので、コードを整理してみます。具体的には、はじめに少し触れたように`sample()`と`ifelse()`が結果に対して不要なので、これを取り除きます。

変更箇所だけ示すとこうです。

```r
for (i in seq_len(n)) {
  M <- dqrnorm(N)
  F <- dqrnorm(N)
  f <- F[F > t]
  Mmean[i] <- mean(M[M > t])
  Fmean[i] <- mean(f[1:(length(f) * p)])
}
```

`result`という変数で計算していた部分は`sum(Mmean > Fmean, na.rm = TRUE)`で計算できます。

で、結果はこうです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109348/a2263601-7404-0936-25d5-dbce1d469544.png)

元の状態から8割程度の実行時間削減ができました。もう「数十秒」と言って良いような気もします。

## Juliaはどうなのか

問題はこちらです。

### Juliaで書いてみる

次のようなコードを書きました。

```julia
using Statistics
using DataFrames
using Gadfly

n = 1000000
N = 1000
p = 0.3
t = 2

function mfmean()
    M = randn(N)
    F = randn(N)
    f = F[F .> t]
    Mmean = mean(M[M .> t])
    Fmean = mean(f[1:Int(round(length(f)*p))])
    return DataFrame(
        category = ["Mmean", "Fmean"],
        value = [Mmean, Fmean]
    )
end

@time result = vcat([mfmean() for i in 1:n]...)
plot(result, x = "value", color = "category", Geom.histogram(bincount = 120, position = :dodge))
```

プロットはきれいに書けたので結果は間違ってないと思います。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109348/c7e27acb-b2dd-a3a2-99ef-166e577e3268.png)

あまり手を加えなくても格好良く出力されて良いですね。ちなみにAtomに[Juno](http://junolab.org/)を入れて書きました。使いやすかったです。

肝心の`@time`で計測した部分ですが、次のような結果になりました。

```
118.985376 seconds (91.30 M allocations: 28.143 GiB, 38.38% gc time)
```

遅いです。数十秒どころではないです。Rの完全勝利…？

ただ、この結果は結構怪しいと思っていて、なぜなら私はJuliaのコードを書き慣れていないからです。Juliaのコードは少し手を加えると劇的に高速化されることがよくあるので、上記のコードにもそういう部分がありそうな気がしています。詳しい方からのご指摘お待ちしております。

### 乱数生成とループ部分だけを比較する

私の適当なコードで速度比較するのも公平性に欠けるので、乱数生成とループの部分だけに注目して比較してみました。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109348/266e3c24-98c8-0c1a-837e-1fbf914a9c68.png)

Juliaの方が少し早いという結果になりました。ただ、Juliaもデフォルトでは乱数生成にメルセンヌ・ツイスタを使用しているようなので（[Random Numbers · The Julia Language](https://docs.julialang.org/en/v1/stdlib/Random/index.html)）、その点を加味すればやはりJuliaは相当早いということが言えそうです。Rで`rnorm()`を使った場合は次のような結果になります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109348/87b5a113-1116-66ca-8ed4-9aa59697fb4e.png)

やはりメルセンヌ・ツイスタ同士ではかなりの差が出ます。

きちんと調べていませんが、JuliaでXoroshiro128+を使うようにすれば相当な速度が出るのかもしれません。

### Julia再チャレンジ

（2019-04-18 21:00追記）

詳しい方からのご指摘が集まってきたので書き直してみました。

```julia
using Statistics
using DataFrames
using RandomNumbers.Xorshifts
const XOROSHIRO128PLUS = Xoroshiro128Plus()

function mfmeans(n)
    N = 1000
    p = 0.3
    t = 2
    Mmeans = Vector{Float64}(undef, n)
    Fmeans = Vector{Float64}(undef, n)
    for i = 1:n
        M = randn(XOROSHIRO128PLUS, N)
        F = randn(XOROSHIRO128PLUS, N)
        f = F[F .> t]
        Mmeans[i] = mean(M[M .> t])
        Fmeans[i] = mean(f[1:round(Int, length(f)*p)])
    end
    return DataFrame(
        category = [repeat(["Mmean"], n); repeat(["Fmean"], n)],
        value = [Mmeans; Fmeans]
    )
end

@time result = mfmeans(1000000)
```

コメント欄で頂いたコードと、[黒木さんのJupyter Notebook](https://nbviewer.jupyter.org/gist/genkuroki/35fad1963a5845a9fbee65bb1d818054)を参考にしました。

結果こうなりました。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109348/3a62db8f-28c5-7104-eea4-33f916d8000f.png)

多少ばらつきがありますが、概ね20秒前後です。やっぱりJuliaの方が数倍程度は早いようです。くやしみ。

## まとめ

- Rでは乱数生成に結構時間がかかっていました。
- 乱数生成部分を`dqrng`パッケージで高速化すると数十秒で計算が終わりました。
- ~~今回の処理はJuliaだとRの倍かかるという若干怪しい結果が出ました。~~
    - 詳しい方々の意見を参考に書き直したらJuliaのが3倍くらい早くなりました。
- 少なくとも「Rでも数十秒で終わる」ことは言えそうです。

## 参考

- [R 3.4.0 のJITバイトコンパイラってどんくらい早くなるの？ - でーたさいえんすって何それ食えるの？](http://masato-613.hatenablog.com/entry/2017/04/25/064632)
- [Performance improvements coming to R 3.4.0 (Revolutions)](https://blog.revolutionanalytics.com/2017/02/preview-r-340.html)
- [Rのforは遅いと誰が言った？ (data.frameの高速演算には列ごとならlapply、行ごとならReduceを使おうの補足) - Qiita](https://qiita.com/Atsushi776/items/c31f2345b9c698354c81)