Clojure の mean を高速化できなかった

まとめ

Clojureの速度が気になったのでJuliaを使ってみたが、Clojureが忘れられず、partitionやreductionsをJuliaで書いていた。そこに手間をかけるのであれば、Clojureの高速化にひと手間かけても良いのではないかと考えた。

しかし、シーケンスの平均値を取る操作を長いベクタに対して行った際、その実行時間を満足いくまで縮めることができなかった。更に縮めるためには恐らくデータ構造を変える必要があり、それはClojureの良さを殺してしまうように思われるので、やはりJuliaを使おう思う。

5000万個のランダム値に対する平均を取る操作の比較が下記。Clojureが4.7秒に対し、Juliaは0.37秒。素のPythonが26秒で、Numpyは0.39秒。Common Lispは大体1.75秒。

;; Clojure
user=> (time (mean (for [_ (range 50000000)] (rand))))
"Elapsed time: 4762.215565 msecs"

# Julia
julia> @time mean([rand() for _ in 1:50000000])
  0.377605 seconds (52.42 k allocations: 383.969 MiB, 22.62% gc time, 9.38% compilation time)

# Python
from time import time
from random import random
from statistics import mean
>>> t = time(); mean([random() for _ in range(0, 50000000)]); time() - t
25.926154136657715

# Numpy
import numpy as np
>>> t = time(); np.random.rand(50000000).mean(); time() - t
0.3986964225769043

; Common Lisp
; 5000万個だと落ちるので、500万個
* (time (alexandria:mean     
    (let ((l nil))
      (dotimes (i 5000000 l)
        (push (random 1.0) l)))))
Evaluation took:
  0.175 seconds of real time
  0.172530 seconds of total run time (0.153032 user, 0.019498 system)
  [ Run times consist of 0.037 seconds GC time, and 0.136 seconds non-GC time. ]
  98.86% CPU
  469,180,265 processor cycles
  79,995,616 bytes consed

Clojure高速化メモ

REPL

大きな出力をさせようとすると、REPLがネックになることが多々ある。心当たりがある場合は、take等で出力範囲を区切る。

プロファイル

まずは、コードのどこが遅いかを見付ける必要がある。選択肢は下記。

time

時間を図りたい部分をtimeで括るだけ。結果の表示がうるさい/遅い場合は、do/nilも併用する。遅延評価部分は実行されないので、必要に応じてdoallも併用する。

(time (for [_ (range 10000)] (range 10000)))
; (out) "Elapsed time: 0.185629 msecs"
; REPLに結果を表示させようとするため、実際はしばらく待たないと結果が表示されず、REPLも使えない。

(time (do (for [_ (range 10000)] (range 10000)) nil))
; (out) "Elapsed time: 0.180818 msecs"
; (do ... nil) で括ると出力が無くなるため、実行後すぐにREPLを使える。

(time (do (doall (for [_ (range 10000)] (range 10000))) nil))
; (out) "Elapsed time: 2.578813 msecs"
; doallを使うと、遅延評価部分が現実化されるため、実行時間が長くなる。
; REPLに表示させる場合は、常にこの現実化+REPL用の処理時間がかかる。

下記のような関数にまとめてしまっても良いかもしれない。

(defn dotime [& something]
  (time (apply doall something))
  nil)

(dotime (for [_ (range 10000)] (range 10000)))
; (out) "Elapsed time: 3.307578 msecs"

tufte

参考記事ではcriterium (☆1.1k) が利用されているが、2021年11月から1年半ほど動きがないので、tufteを試す。こちらは2022年12月更新の、☆486。READMEを見る限り、こちらの方がcriteriumよりも細かく情報を取れそう。
ターミナルからREPLでcljしている場合は良いが、Conjureを利用している場合、Conjure出力画面ではなく裏で動いているターミナルに結果が表示されるので注意。

;; deps.edn
  :deps {
    com.taoensso/tufte {:mvn/version "2.4.5"}
  }

計測した部分を(p ...)で括って、最初の引数にキーワードを渡してあげるだけ。最後に全体を(profile {} ...)で括る。

(ns user
  (:require [taoensso.tufte :as tufte :refer [defnp p profile]]))

(tufte/add-basic-println-handler!
  {:format-pstats-opts {:columns [:n-calls :min :max :mean :mad :clock :total]}})

(profile
  {}
  (let [a (p :a (range 10000))]
    (p :body (range 10000))))
;; (省略)
;; pId           nCalls        Min        Max       Mean   MAD      Clock  Total
;; 
;; :a                 1     1.60μs     1.60μs     1.60μs   ±0%     1.60μs     1%
;; :body              1   387.00ns   387.00ns   387.00ns   ±0%   387.00ns     0%
;; 
;; Accounted                                                       1.99μs     1%
;; Clock                                                         180.80μs   100%

なお、profileの後ろの{}には、プロファイルのオプションを入れることができる。そしてtufteについても、関数を作っておく。

(defn doprof [& something]
  (profile {} (apply doall something) nil))

高速化

mean

Clojureは標準でmeanが実装されていないので、自分で簡単なものを書いていたが、これがボトルネックだった。なので、逐次的に値を求めるような書き方に変更をした。
下記の通りいくつかのパターンを試した (実際はそれぞれ何回かサンプルを取った) ところ、ライブラリのものが少し他より速かった。

;; これが遅い
(defn mean [coll]
  (/ (apply + coll) (count coll)))
; (time (mean (for [_ (range 1000000)] (rand))))
; (out) "Elapsed time: 106.225616 msecs"

;; reduceで一回だけcollを眺めれば良いように変更
(defn mean [coll]
  (let [r (reduce (fn [mem val] {:S (+ (mem :S) val) :n (inc (mem :n))})
                  {:S 0 :n 0}
                  coll)]
    (/ (r :S) (r :n))))
; (out) "Elapsed time: 70.677517 msecs"

;; recurバージョン
(defn mean 
  ([coll] (mean coll 0 0))
  ([coll s n]
    (if (empty? coll)
      (/ s n)
      (recur (rest coll) (+ s (first coll)) (inc n)))))
; (out) "Elapsed time: 79.161107 msecs"

;; recurベクタバージョン
(defn mean 
  ([coll] (mean (vec coll) 0 0))
  ([coll s n]
    (if (empty? coll)
      (/ s n)
      (recur (pop coll) (+ s (peek coll)) (inc n)))))
; (out) "Elapsed time: 207.029417 msecs"

;; ライブラリ
;; net.mikera/core.matrix {:mvn/version "0.63.0"}
(ns user
  (:require [clojure.core.matrix.stats :refer [mean]]))
; (out) "Elapsed time: 70.736877 msecs"

わざわざ書くほどでもないが、メモ。$M$=平均、$S$=合計。
$$
\begin{align*}
M_n &= \frac{a_1 + a_2 + \cdots + a_n}{n} = \frac{S_n}{n} \\
S_n &= S_{n-1}+a_n \\
M_n &= \frac{S_{n-1}+a_n}{n}
\end{align*}
$$

上記の速度はJuliaに比べると、十分に遅いように思われる。

# Julia
using Statistics
@time mean([rand() for _ in 1:1000000])
# 0.033709 seconds (52.42 k allocations: 10.131 MiB, 83.37% compilation time)