非復元抽出の高速かつ実装が簡単な方法を考える
※ @tomerun さんに書いてもらったコードとその検証結果を記事の最後に追記しました.(2013-07-21 2:00)
ふとしたきっかけで非復元抽出 (random sampling without replacement) を実装するときに気になったのでどんな実装がよいのか考えてみた.なお非復元抽出はビンゴのように,N個の要素の中からk個の異なる要素をランダムに選択するという意味である.
復元抽出については @unnonouno さんのブログなどに書いてあり,非復元抽出についてもリンクが張ってあったのだけれど,リンク先のブログ記事が読めない状態になっていていたのが残念.
はじめに
std::vector
今回は以下の3つの方法を考えた.
- (1) std::setを用いる方法
- (2) std::tr1::unordered_setを用いる方法
- (3) std::vector + std::random_shuffleを用いる方法
- ※ 本ブログ記事をほぼ書き終わったところでKnuth先生の方法を知ったのであとで追記 orz
(1),(2)は元の配列からランダムに選択してset/unordered_setにinsertして,コンテナの大きさがkになったら停止するというもの.(3)は元の配列の各要素へのポインタを保持したvectorを作成し,それをstd::random_shuffleを用いてシャッフルして先頭からk個取得すればよい.
何も考えずに実装すると,ランダムに要素を選択してsetに挿入する方法を思いつく.しかしsetの内部はmapと同様二分木 (赤黒木) で実装されているので,要素が増えるとコストが大きくなると予想される.というわけでハッシュ実装であるunordered_setの利用を思いつく.そしてNの要素を入れた配列をシャッフルすることができればkの数に依存せずに一定時間で任意の数の非復元抽出ができるじゃないかと思って3番目の方法を思いつく.
きっとここらへんが誰でも思いつくレベルかつ実装が簡単な方法ではないだろうか.というわけでこの3つの方法の速度の差を検証してみる.
実験
実験条件
3つの手法による処理時間を以下のNとkの組み合わせについて10回ずつ計測.今回はkを固定の数字ではなく比率とした.すなわちN=1000, k=0.1の場合には1000個の要素の中から100個の要素を非復元抽出する.
- N \in {10^3, 10^4,..., 10^8}
- k \in {0.1, 0.2, ..., 0.9}
- 測定方法
- gettimefoday(2)で処理時間を計測
- ※ random_shuffleの場合はstd::vector
を構築する部分も時間に含める
実験に利用したコードの一部は以下のとおり.
#define RAND (((double)rand()-1)/((double)RAND_MAX)) double test_set (int N, int k) { std::vector<int> vec(N); for (int i = 0; i < N; i++) { vec[ i ] = i; } double t1 = gettimeofday_sec(); std::set<int> s; while ((int)s.size() < k) { int idx = (int)(RAND * N); s.insert( vec[ idx ] ); } double t2 = gettimeofday_sec(); return (t2 - t1); } double test_hash (int N, int k) { std::vector<int> vec(N); for (int i = 0; i < N; i++) { vec[ i ] = i; } double t1 = gettimeofday_sec(); std::tr1::unordered_set<int> s; while ((int)s.size() < k) { int idx = (int)(RAND * N); s.insert( vec[ idx ] ); } double t2 = gettimeofday_sec(); return (t2 - t1); } double test_shuffle (int N, int k) { std::vector<int> vec(N); for (int i = 0; i < N; i++) { vec[ i ] = i; } double t1 = gettimeofday_sec(); std::vector<int*> shuffle_vec(N); for (int i = 0; i < (int)vec.size(); i++) { shuffle_vec[ i ] = &(vec[ i ]); } std::random_shuffle( shuffle_vec.begin(), shuffle_vec.end() ); double t2 = gettimeofday_sec(); return (t2 - t1); }
結果
まずNとkに対する全ての結果を示す.表が見づらくなるので平均値のみを記載.(分散は結論に影響を与えない程度の大きさだった)
N | k (ratio) | set | unordered_set | random_shuffle |
---|---|---|---|---|
1000 | 0.1 | 0.000026 | 0.000023 | 0.000042 |
1000 | 0.2 | 0.000055 | 0.000045 | 0.000042 |
1000 | 0.3 | 0.000086 | 0.000062 | 0.000043 |
1000 | 0.4 | 0.000123 | 0.000081 | 0.000043 |
1000 | 0.5 | 0.000143 | 0.000054 | 0.000021 |
1000 | 0.6 | 0.000104 | 0.000067 | 0.000021 |
1000 | 0.7 | 0.000134 | 0.000078 | 0.000020 |
1000 | 0.8 | 0.000175 | 0.000092 | 0.000021 |
1000 | 0.9 | 0.000245 | 0.000122 | 0.000021 |
10000 | 0.1 | 0.000142 | 0.000097 | 0.000215 |
10000 | 0.2 | 0.000310 | 0.000206 | 0.000205 |
10000 | 0.3 | 0.000506 | 0.000290 | 0.000217 |
10000 | 0.4 | 0.000744 | 0.000461 | 0.000205 |
10000 | 0.5 | 0.001025 | 0.000543 | 0.000207 |
10000 | 0.6 | 0.001380 | 0.000665 | 0.000212 |
10000 | 0.7 | 0.001823 | 0.000809 | 0.000205 |
10000 | 0.8 | 0.002451 | 0.001154 | 0.000211 |
10000 | 0.9 | 0.003540 | 0.001388 | 0.000212 |
100000 | 0.1 | 0.0020 | 0.0012 | 0.0022 |
100000 | 0.2 | 0.0046 | 0.0025 | 0.0022 |
100000 | 0.3 | 0.0083 | 0.0039 | 0.0022 |
100000 | 0.4 | 0.0123 | 0.0058 | 0.0022 |
100000 | 0.5 | 0.0171 | 0.0073 | 0.0022 |
100000 | 0.6 | 0.0232 | 0.0090 | 0.0022 |
100000 | 0.7 | 0.0311 | 0.0121 | 0.0022 |
100000 | 0.8 | 0.0421 | 0.0143 | 0.0022 |
100000 | 0.9 | 0.0606 | 0.0177 | 0.0022 |
1000000 | 0.1 | 0.0324 | 0.0153 | 0.0265 |
1000000 | 0.2 | 0.0819 | 0.0364 | 0.0262 |
1000000 | 0.3 | 0.1547 | 0.0774 | 0.0264 |
1000000 | 0.4 | 0.2447 | 0.1005 | 0.0264 |
1000000 | 0.5 | 0.3583 | 0.1280 | 0.0263 |
1000000 | 0.6 | 0.5024 | 0.1990 | 0.0265 |
1000000 | 0.7 | 0.6924 | 0.2338 | 0.0267 |
1000000 | 0.8 | 0.9639 | 0.2818 | 0.0268 |
1000000 | 0.9 | 1.4319 | 0.3625 | 0.0264 |
10000000 | 0.1 | 0.7067 | 0.3181 | 0.6920 |
10000000 | 0.2 | 1.7642 | 0.7155 | 0.6909 |
10000000 | 0.3 | 3.0718 | 1.2310 | 0.6921 |
10000000 | 0.4 | 4.6572 | 1.5806 | 0.6909 |
10000000 | 0.5 | 6.5762 | 2.3766 | 0.6904 |
10000000 | 0.6 | 8.9766 | 2.7640 | 0.6929 |
10000000 | 0.7 | 12.1109 | 3.2506 | 0.6924 |
10000000 | 0.8 | 16.5619 | 3.9203 | 0.6911 |
10000000 | 0.9 | 24.2751 | 5.9600 | 0.6909 |
100000000 | 0.1 | 11.9186 | 5.0229 | 9.4275 |
100000000 | 0.2 | 28.4891 | 10.7473 | 9.4202 |
100000000 | 0.3 | 48.9235 | 14.4396 | 9.4179 |
100000000 | 0.4 | 73.7326 | 23.2700 | 9.4068 |
100000000 | 0.5 | 103.9370 | 27.3172 | 9.4198 |
100000000 | 0.6 | 141.7940 | 32.3755 | 9.4193 |
100000000 | 0.7 | 191.4050 | 38.9861 | 9.4176 |
100000000 | 0.8 | 262.2190 | 57.5724 | 9.4168 |
100000000 | 0.9 | 384.7380 | 71.8925 | 9.4186 |
手法毎にNとkの組み合わせについての3次元棒グラフは以下のとおり.上段が10^3~6まで.下段が10^3~5までの結果を全て同じ尺度で示している.
実験より,以下の結果を確認した.
- k=0.1の場合に実行速度は unordered_set > random_shuffle > set
- kが0.2以上の場合に実行速度は random_shuffle > unordered_set > set
- なぜかN=1000,k=0.1-0.4のrandom_shuffleが遅い.(何回かやり直しても再現)
なおrandom_shuffle実装においてポインタ要素の配列を作成する処理を時間計測から外した場合においても上記の結果は変わらなかった.
まとめ
上記の結果より,以下の知見を得た.
- あらゆる状況においてsetは遅い.
- 非復元抽出を行う要素数が少ない場合には(2)の方法を用いる
- 復元抽出を行う要素数が抽出元サイズの2割以上の場合には(3)の方法を用いる
- ポインタ要素の配列作成のコストが小さいことから配列操作は高速ということがわかる.
やはりtreeはポインタをたどるのでキャッシュミスが発生しやすい->遅い,という印象がさらに強まる結果であった.N=1000,k=0.1-0.4のrandom_shuffleが遅いのはrandom_shuffleの実装を調べればわかるのだろうか? 配列の大きさが小さい場合には,ランダム性を保証するためにたくさんシャッフルするとか?
さよならset,君のことは忘れないよ.(あれ,どっかで聞いたことあるセリフ?)
発展手法について
どうやらKnuth先生の本に非復元抽出の高速なアルゴリズムが載っているらしい.
(参考: http://stackoverflow.com/questions/311703/algorithm-for-sampling-without-replacement)
今回の実験に合わせた表記にするとこんな感じ.
double test_knuth (int N, int k) { std::vector<int> vec(N); for (int i = 0; i < N; i++) { vec[ i ] = i; } double t1 = gettimeofday_sec(); std::vector<int> sample_vec(k); int t = 0; int m = 0; while (m < k) { double u = RAND; if ((N - t) * u >= (k - m)) { t++; } else { sample_vec[ m ] = vec[ t ]; t++; m++; } } double t2 = gettimeofday_sec(); return (t2 - t1); }
ざっと試してみたところ,どうやら一番高速そう.しかも実装も簡単.今回の実験で表とグラフを作ったあとにこの方法を知ったので相当げんなりしている.グラフを作り直す元気はなかったので,更新した表を最後に追記.
また今回は簡単な実装ということで効率的な実装について深追いしなかったけれど,論文になっている手法では以下のようなものがあるらしい.ざっと読んだけれどまったく理解できなかったので見なかったことにした.
- J. S. Vitter, "An efficient algorithm for sequential random sampling", ACM TOMS vol.13(1), pp.58-67, 1987
Appendix
Knuthの方法の結果を以下に追記.
結果は以下の解釈をしている.
- Nとkの数が共に小さい場合にはunordered_set
- それ以外の場合にはKnuth
結論: Knuth先生すごい!
N | k (ratio) | set | unordered_set | random_shuffle | knuth |
---|---|---|---|---|---|
1000 | 0.1 | 0.000026 | 0.000023 | 0.000042 | 0.000025 |
1000 | 0.2 | 0.000055 | 0.000045 | 0.000042 | 0.000030 |
1000 | 0.3 | 0.000086 | 0.000062 | 0.000043 | 0.000034 |
1000 | 0.4 | 0.000123 | 0.000081 | 0.000043 | 0.000038 |
1000 | 0.5 | 0.000143 | 0.000054 | 0.000021 | 0.000040 |
1000 | 0.6 | 0.000104 | 0.000067 | 0.000021 | 0.000118 |
1000 | 0.7 | 0.000134 | 0.000078 | 0.000020 | 0.000036 |
1000 | 0.8 | 0.000175 | 0.000092 | 0.000021 | 0.000032 |
1000 | 0.9 | 0.000245 | 0.000122 | 0.000021 | 0.000028 |
10000 | 0.1 | 0.000142 | 0.000097 | 0.000215 | 0.000253 |
10000 | 0.2 | 0.000310 | 0.000206 | 0.000205 | 0.000300 |
10000 | 0.3 | 0.000506 | 0.000290 | 0.000217 | 0.000184 |
10000 | 0.4 | 0.000744 | 0.000461 | 0.000205 | 0.000186 |
10000 | 0.5 | 0.001025 | 0.000543 | 0.000207 | 0.000195 |
10000 | 0.6 | 0.001380 | 0.000665 | 0.000212 | 0.000189 |
10000 | 0.7 | 0.001823 | 0.000809 | 0.000205 | 0.000170 |
10000 | 0.8 | 0.002451 | 0.001154 | 0.000211 | 0.000150 |
10000 | 0.9 | 0.003540 | 0.001388 | 0.000212 | 0.000134 |
100000 | 0.1 | 0.0020 | 0.0012 | 0.0022 | 0.0012 |
100000 | 0.2 | 0.0046 | 0.0025 | 0.0022 | 0.0014 |
100000 | 0.3 | 0.0083 | 0.0039 | 0.0022 | 0.0017 |
100000 | 0.4 | 0.0123 | 0.0058 | 0.0022 | 0.0018 |
100000 | 0.5 | 0.0171 | 0.0073 | 0.0022 | 0.0019 |
100000 | 0.6 | 0.0232 | 0.0090 | 0.0022 | 0.0019 |
100000 | 0.7 | 0.0311 | 0.0121 | 0.0022 | 0.0017 |
100000 | 0.8 | 0.0421 | 0.0143 | 0.0022 | 0.0015 |
100000 | 0.9 | 0.0606 | 0.0177 | 0.0022 | 0.0013 |
1000000 | 0.1 | 0.0324 | 0.0153 | 0.0265 | 0.0121 |
1000000 | 0.2 | 0.0819 | 0.0364 | 0.0262 | 0.0144 |
1000000 | 0.3 | 0.1547 | 0.0774 | 0.0264 | 0.0166 |
1000000 | 0.4 | 0.2447 | 0.1005 | 0.0264 | 0.0184 |
1000000 | 0.5 | 0.3583 | 0.1280 | 0.0263 | 0.0192 |
1000000 | 0.6 | 0.5024 | 0.1990 | 0.0265 | 0.0187 |
1000000 | 0.7 | 0.6924 | 0.2338 | 0.0267 | 0.0170 |
1000000 | 0.8 | 0.9639 | 0.2818 | 0.0268 | 0.0150 |
1000000 | 0.9 | 1.4319 | 0.3625 | 0.0264 | 0.0129 |
10000000 | 0.1 | 0.7067 | 0.3181 | 0.6920 | 0.1222 |
10000000 | 0.2 | 1.7642 | 0.7155 | 0.6909 | 0.1451 |
10000000 | 0.3 | 3.0718 | 1.2310 | 0.6921 | 0.1675 |
10000000 | 0.4 | 4.6572 | 1.5806 | 0.6909 | 0.1860 |
10000000 | 0.5 | 6.5762 | 2.3766 | 0.6904 | 0.1943 |
10000000 | 0.6 | 8.9766 | 2.7640 | 0.6929 | 0.1887 |
10000000 | 0.7 | 12.1109 | 3.2506 | 0.6924 | 0.1729 |
10000000 | 0.8 | 16.5619 | 3.9203 | 0.6911 | 0.1529 |
10000000 | 0.9 | 24.2751 | 5.9600 | 0.6909 | 0.1368 |
100000000 | 0.1 | 11.9186 | 5.0229 | 9.4275 | 1.2284 |
100000000 | 0.2 | 28.4891 | 10.7473 | 9.4202 | 1.4638 |
100000000 | 0.3 | 48.9235 | 14.4396 | 9.4179 | 1.6919 |
100000000 | 0.4 | 73.7326 | 23.2700 | 9.4068 | 1.8807 |
100000000 | 0.5 | 103.9370 | 27.3172 | 9.4198 | 1.9683 |
100000000 | 0.6 | 141.7940 | 32.3755 | 9.4193 | 1.9175 |
100000000 | 0.7 | 191.4050 | 38.9861 | 9.4176 | 1.7639 |
100000000 | 0.8 | 262.2190 | 57.5724 | 9.4168 | 1.5684 |
100000000 | 0.9 | 384.7380 | 71.8925 | 9.4186 | 1.3664 |
2013-07-21 2:00 追記 (Thanks to @tomerun さん)
ブログ記事upのつぶやきをしたら @tomerun さんから以下の返信を頂く.
ありがとうございます! @tomerun さんに書いてもらったコードの該当部分だけを転載します.
// Code by @tomerun さん // http://ideone.com/raF3O4 から引用 double test_shuffle_partial (int N, int k) { std::vector<int> vec(N); for (int i = 0; i < N; i++) { vec[ i ] = i; } double t1 = gettimeofday_sec(); std::vector<int> shuffle_vec; for (int i = 0; i < k; ++i) { int pos = rand() % (N-i); shuffle_vec.push_back(vec[i + pos]); swap(vec[i], vec[i + pos]); } double t2 = gettimeofday_sec(); return (t2 - t1); }
これを見ると,選択された要素以外からランダムに選ぶ.選択されたものを配列の先頭に移して,それ以外からランダムに要素を選ぶ,ということを繰り返してk個のアイテムを取得する.現実世界におけるビンゴのような非復元抽出手法をコードに落とし込んだ感じ.なるほど,こう実装すればよいのか! 目からうろこです.
なお,超細かいところで恐縮なのだけれどpush_backするとkが大きい場合に何度もvectorのreserve (realloc) が走るので,std::vector
さぁさぁ同じ環境で実行して以下の結果を得た.random_shuffle_partial以外は前回の数字.
N | k (ratio) | set | unordered_set | random_shuffle | knuth | random_shuffle_partial | |
---|---|---|---|---|---|---|---|
1000 | 0.1 | 0.000026 | 0.000023 | 0.000042 | 0.000025 | 0.000005 | |
1000 | 0.2 | 0.000055 | 0.000045 | 0.000042 | 0.000030 | 0.000007 | |
1000 | 0.3 | 0.000086 | 0.000062 | 0.000043 | 0.000034 | 0.000010 | |
1000 | 0.4 | 0.000123 | 0.000081 | 0.000043 | 0.000038 | 0.000013 | |
1000 | 0.5 | 0.000143 | 0.000054 | 0.000021 | 0.000040 | 0.000015 | |
1000 | 0.6 | 0.000104 | 0.000067 | 0.000021 | 0.000118 | 0.000019 | |
1000 | 0.7 | 0.000134 | 0.000078 | 0.000020 | 0.000036 | 0.000021 | |
1000 | 0.8 | 0.000175 | 0.000092 | 0.000021 | 0.000032 | 0.000024 | |
1000 | 0.9 | 0.000245 | 0.000122 | 0.000021 | 0.000028 | 0.000026 | |
10000 | 0.1 | 0.000142 | 0.000097 | 0.000215 | 0.000253 | 0.000113 | |
10000 | 0.2 | 0.000310 | 0.000206 | 0.000205 | 0.000300 | 0.000059 | |
10000 | 0.3 | 0.000506 | 0.000290 | 0.000217 | 0.000184 | 0.000087 | |
10000 | 0.4 | 0.000744 | 0.000461 | 0.000205 | 0.000186 | 0.000112 | |
10000 | 0.5 | 0.001025 | 0.000543 | 0.000207 | 0.000195 | 0.000091 | |
10000 | 0.6 | 0.001380 | 0.000665 | 0.000212 | 0.000189 | 0.000081 | |
10000 | 0.7 | 0.001823 | 0.000809 | 0.000205 | 0.000170 | 0.000096 | |
10000 | 0.8 | 0.002451 | 0.001154 | 0.000211 | 0.000150 | 0.000113 | |
10000 | 0.9 | 0.003540 | 0.001388 | 0.000212 | 0.000134 | 0.000128 | |
100000 | 0.1 | 0.0020 | 0.0012 | 0.0022 | 0.0012 | 0.0002 | |
100000 | 0.2 | 0.0046 | 0.0025 | 0.0022 | 0.0014 | 0.0003 | |
100000 | 0.3 | 0.0083 | 0.0039 | 0.0022 | 0.0017 | 0.0005 | |
100000 | 0.4 | 0.0123 | 0.0058 | 0.0022 | 0.0018 | 0.0008 | |
100000 | 0.5 | 0.0171 | 0.0073 | 0.0022 | 0.0019 | 0.0009 | |
100000 | 0.6 | 0.0232 | 0.0090 | 0.0022 | 0.0019 | 0.0011 | |
100000 | 0.7 | 0.0311 | 0.0121 | 0.0022 | 0.0017 | 0.0013 | |
100000 | 0.8 | 0.0421 | 0.0143 | 0.0022 | 0.0015 | 0.0014 | |
100000 | 0.9 | 0.0606 | 0.0177 | 0.0022 | 0.0013 | 0.0015 | |
1000000 | 0.1 | 0.0324 | 0.0153 | 0.0265 | 0.0121 | 0.0020 | |
1000000 | 0.2 | 0.0819 | 0.0364 | 0.0262 | 0.0144 | 0.0040 | |
1000000 | 0.3 | 0.1547 | 0.0774 | 0.0264 | 0.0166 | 0.0068 | |
1000000 | 0.4 | 0.2447 | 0.1005 | 0.0264 | 0.0184 | 0.0089 | |
1000000 | 0.5 | 0.3583 | 0.1280 | 0.0263 | 0.0192 | 0.0109 | |
1000000 | 0.6 | 0.5024 | 0.1990 | 0.0265 | 0.0187 | 0.0136 | |
1000000 | 0.7 | 0.6924 | 0.2338 | 0.0267 | 0.0170 | 0.0155 | |
1000000 | 0.8 | 0.9639 | 0.2818 | 0.0268 | 0.0150 | 0.0174 | |
1000000 | 0.9 | 1.4319 | 0.3625 | 0.0264 | 0.0129 | 0.0191 | |
10000000 | 0.1 | 0.7067 | 0.3181 | 0.6920 | 0.1222 | 0.0514 | |
10000000 | 0.2 | 1.7642 | 0.7155 | 0.6909 | 0.1451 | 0.1022 | |
10000000 | 0.3 | 3.0718 | 1.2310 | 0.6921 | 0.1675 | 0.1537 | |
10000000 | 0.4 | 4.6572 | 1.5806 | 0.6909 | 0.1860 | 0.2010 | |
10000000 | 0.5 | 6.5762 | 2.3766 | 0.6904 | 0.1943 | 0.2535 | |
10000000 | 0.6 | 8.9766 | 2.7640 | 0.6929 | 0.1887 | 0.2976 | |
10000000 | 0.7 | 12.1109 | 3.2506 | 0.6924 | 0.1729 | 0.3389 | |
10000000 | 0.8 | 16.5619 | 3.9203 | 0.6911 | 0.1529 | 0.3756 | |
10000000 | 0.9 | 24.2751 | 5.9600 | 0.6909 | 0.1368 | 0.4156 | |
100000000 | 0.1 | 11.9186 | 5.0229 | 9.4275 | 1.2284 | 0.6073 | |
100000000 | 0.2 | 28.4891 | 10.7473 | 9.4202 | 1.4638 | 1.2054 | |
100000000 | 0.3 | 48.9235 | 14.4396 | 9.4179 | 1.6919 | 1.7780 | |
100000000 | 0.4 | 73.7326 | 23.2700 | 9.4068 | 1.8807 | 2.3940 | |
100000000 | 0.5 | 103.9370 | 27.3172 | 9.4198 | 1.9683 | 2.9615 | |
100000000 | 0.6 | 141.7940 | 32.3755 | 9.4193 | 1.9175 | 3.5221 | |
100000000 | 0.7 | 191.4050 | 38.9861 | 9.4176 | 1.7639 | 4.1683 | |
100000000 | 0.8 | 262.2190 | 57.5724 | 9.4168 | 1.5684 | 4.7086 | |
100000000 | 0.9 | 384.7380 | 71.8925 | 9.4186 | 1.3664 | 5.2277 |
さて実験から以下の結果を得た.
- N<=10^6までは全ての手法でrandom_shuffle_partialが最速
- 全てのN,kにおいてrandom_shuffle_partial>random_shuffle
- Knuth vs. random_shuffle_parial 最速対決
Nとkが大きくなるとknuth先生に負けてしまう.おそらくkの絶対数が増えるとswapコストが効いてくるのだろうか.しかし,実験のためk=0.1-0.9を検証しているけれど,実際にはkはそんなに大きくない場合の方が多いので,@tomerun さん手法の方が高速な場面の方が多いと思われる.
改めて@tomerun さんに感謝.さすができる人は対応もコードも3倍の速さですね!