非復元抽出の高速かつ実装が簡単な方法を考える

c++

※ @tomerun さんに書いてもらったコードとその検証結果を記事の最後に追記しました．(2013-07-21 2:00)

ふとしたきっかけで非復元抽出 (random sampling without replacement) を実装するときに気になったのでどんな実装がよいのか考えてみた．なお非復元抽出はビンゴのように，N個の要素の中からk個の異なる要素をランダムに選択するという意味である．

復元抽出については @unnonouno さんのブログなどに書いてあり，非復元抽出についてもリンクが張ってあったのだけれど，リンク先のブログ記事が読めない状態になっていていたのが残念．

はじめに

std::vectorのようにN個のオブジェクトを格納した配列からk個の異なるオブジェクトをランダムに選択したいというような状況を考える．このような処理を何回も行うため，可能な限り高速に処理したい．

今回は以下の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を構築する部分も時間に含める
  • 実験環境
    • gcc 4.4.5 (-O3オプション)
    • Linux kernel 2.6.32-5
    • Core i7 950 (3.07GHz)

実験に利用したコードの一部は以下のとおり．

#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);
}

ざっと試してみたところ，どうやら一番高速そう．しかも実装も簡単．今回の実験で表とグラフを作ったあとにこの方法を知ったので相当げんなりしている．グラフを作り直す元気はなかったので，更新した表を最後に追記．

また今回は簡単な実装ということで効率的な実装について深追いしなかったけれど，論文になっている手法では以下のようなものがあるらしい．ざっと読んだけれどまったく理解できなかったので見なかったことにした．

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 shuffle_vec(k)として，push_backの代わりにshuffle_vec[ i ] = vec[ i + pos ];とするとkが大きいときに少し速くなる．N=10^8,k=0.9で5.2277->4.972程度．まぁやればできる系の高速化なのでtrivialですが．

さぁさぁ同じ環境で実行して以下の結果を得た．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=10^7,k<0.4においてrandom_shuffle_partial>Knuth>その他
    • N=10^7,k>=0.4においてKnuth>random_shuffle_partial>その他
    • N=10^8,k<0.3においてrandom_shuffle_partial>Knuth>その他
    • N=10^8,k>=0.3においてKnuth>random_shuffle_partial>その他

Nとkが大きくなるとknuth先生に負けてしまう．おそらくkの絶対数が増えるとswapコストが効いてくるのだろうか．しかし，実験のためk=0.1-0.9を検証しているけれど，実際にはkはそんなに大きくない場合の方が多いので，@tomerun さん手法の方が高速な場面の方が多いと思われる．

改めて@tomerun さんに感謝．さすができる人は対応もコードも3倍の速さですね!