睡眠不足?! RSSフィード

2269098

2013-07-20

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

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

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

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

はじめに

std::vector<Object>のように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<int *>を構築する部分も時間に含める

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

#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に対する全ての結果を示す.表が見づらくなるので平均値のみを記載.(分散は結論に影響を与えない程度の大きさだった)

Nk (ratio)setunordered_setrandom_shuffle
10000.10.0000260.0000230.000042
10000.20.0000550.0000450.000042
10000.30.0000860.0000620.000043
10000.40.0001230.0000810.000043
10000.50.0001430.0000540.000021
10000.60.0001040.0000670.000021
10000.70.0001340.0000780.000020
10000.80.0001750.0000920.000021
10000.90.0002450.0001220.000021
100000.10.0001420.0000970.000215
100000.20.0003100.0002060.000205
100000.30.0005060.0002900.000217
100000.40.0007440.0004610.000205
100000.50.0010250.0005430.000207
100000.60.0013800.0006650.000212
100000.70.0018230.0008090.000205
100000.80.0024510.0011540.000211
100000.90.0035400.0013880.000212
1000000.10.00200.00120.0022
1000000.20.00460.00250.0022
1000000.30.00830.00390.0022
1000000.40.01230.00580.0022
1000000.50.01710.00730.0022
1000000.60.02320.00900.0022
1000000.70.03110.01210.0022
1000000.80.04210.01430.0022
1000000.90.06060.01770.0022
10000000.10.0324 0.0153 0.0265
10000000.20.0819 0.0364 0.0262
10000000.30.1547 0.0774 0.0264
10000000.40.2447 0.1005 0.0264
10000000.50.3583 0.1280 0.0263
10000000.60.5024 0.1990 0.0265
10000000.70.6924 0.2338 0.0267
10000000.80.9639 0.2818 0.0268
10000000.91.4319 0.3625 0.0264
100000000.10.7067 0.3181 0.6920
100000000.21.7642 0.7155 0.6909
100000000.33.0718 1.2310 0.6921
100000000.44.6572 1.5806 0.6909
100000000.56.5762 2.3766 0.6904
100000000.68.9766 2.7640 0.6929
100000000.712.1109 3.2506 0.6924
100000000.816.5619 3.9203 0.6911
100000000.924.2751 5.9600 0.6909
1000000000.111.9186 5.0229 9.4275
1000000000.228.4891 10.7473 9.4202
1000000000.348.9235 14.4396 9.4179
1000000000.473.7326 23.2700 9.4068
1000000000.5103.9370 27.3172 9.4198
1000000000.6141.7940 32.3755 9.4193
1000000000.7191.4050 38.9861 9.4176
1000000000.8262.2190 57.5724 9.4168
1000000000.9384.7380 71.8925 9.4186

手法毎にNとkの組み合わせについての3次元棒グラフは以下のとおり.上段が10^3~6まで.下段が10^3~5までの結果を全て同じ尺度で示している.

f:id:sleepy_yoshi:20130721001120p:image

実験より,以下の結果を確認した.

  • 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先生すごい!

Nk (ratio)setunordered_setrandom_shuffleknuth
10000.10.0000260.0000230.0000420.000025
10000.20.0000550.0000450.0000420.000030
10000.30.0000860.0000620.0000430.000034
10000.40.0001230.0000810.0000430.000038
10000.50.0001430.0000540.0000210.000040
10000.60.0001040.0000670.0000210.000118
10000.70.0001340.0000780.0000200.000036
10000.80.0001750.0000920.0000210.000032
10000.90.0002450.0001220.0000210.000028
100000.10.0001420.0000970.0002150.000253
100000.20.0003100.0002060.0002050.000300
100000.30.0005060.0002900.0002170.000184
100000.40.0007440.0004610.0002050.000186
100000.50.0010250.0005430.0002070.000195
100000.60.0013800.0006650.0002120.000189
100000.70.0018230.0008090.0002050.000170
100000.80.0024510.0011540.0002110.000150
100000.90.0035400.0013880.0002120.000134
1000000.10.00200.00120.00220.0012
1000000.20.00460.00250.00220.0014
1000000.30.00830.00390.00220.0017
1000000.40.01230.00580.00220.0018
1000000.50.01710.00730.00220.0019
1000000.60.02320.00900.00220.0019
1000000.70.03110.01210.00220.0017
1000000.80.04210.01430.00220.0015
1000000.90.06060.01770.00220.0013
10000000.10.03240.01530.02650.0121
10000000.20.08190.03640.02620.0144
10000000.30.15470.07740.02640.0166
10000000.40.24470.10050.02640.0184
10000000.50.35830.12800.02630.0192
10000000.60.50240.19900.02650.0187
10000000.70.69240.23380.02670.0170
10000000.80.96390.28180.02680.0150
10000000.91.43190.36250.02640.0129
100000000.10.70670.31810.69200.1222
100000000.21.76420.71550.69090.1451
100000000.33.07181.23100.69210.1675
100000000.44.65721.58060.69090.1860
100000000.56.57622.37660.69040.1943
100000000.68.97662.76400.69290.1887
100000000.712.11093.25060.69240.1729
100000000.816.56193.92030.69110.1529
100000000.924.27515.96000.69090.1368
1000000000.111.91865.02299.42751.2284
1000000000.228.489110.74739.42021.4638
1000000000.348.923514.43969.41791.6919
1000000000.473.732623.27009.40681.8807
1000000000.5103.937027.31729.41981.9683
1000000000.6141.794032.37559.41931.9175
1000000000.7191.405038.98619.41761.7639
1000000000.8262.219057.57249.41681.5684
1000000000.9384.738071.89259.41861.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<int> 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以外は前回の数字.

Nk (ratio)setunordered_setrandom_shuffleknuthrandom_shuffle_partial
10000.10.000026 0.000023 0.000042 0.000025 0.000005
10000.20.000055 0.000045 0.000042 0.000030 0.000007
10000.30.000086 0.000062 0.000043 0.000034 0.000010
10000.40.000123 0.000081 0.000043 0.000038 0.000013
10000.50.000143 0.000054 0.000021 0.000040 0.000015
10000.60.000104 0.000067 0.000021 0.000118 0.000019
10000.70.000134 0.000078 0.000020 0.000036 0.000021
10000.80.000175 0.000092 0.000021 0.000032 0.000024
10000.90.000245 0.000122 0.000021 0.000028 0.000026
100000.10.000142 0.000097 0.000215 0.000253 0.000113
100000.20.000310 0.000206 0.000205 0.000300 0.000059
100000.30.000506 0.000290 0.000217 0.000184 0.000087
100000.40.000744 0.000461 0.000205 0.000186 0.000112
100000.50.001025 0.000543 0.000207 0.000195 0.000091
100000.60.001380 0.000665 0.000212 0.000189 0.000081
100000.70.001823 0.000809 0.000205 0.000170 0.000096
100000.80.002451 0.001154 0.000211 0.000150 0.000113
100000.90.003540 0.001388 0.000212 0.000134 0.000128
1000000.10.0020 0.0012 0.0022 0.0012 0.0002
1000000.20.0046 0.0025 0.0022 0.0014 0.0003
1000000.30.0083 0.0039 0.0022 0.0017 0.0005
1000000.40.0123 0.0058 0.0022 0.0018 0.0008
1000000.50.0171 0.0073 0.0022 0.0019 0.0009
1000000.60.0232 0.0090 0.0022 0.0019 0.0011
1000000.70.0311 0.0121 0.0022 0.0017 0.0013
1000000.80.0421 0.0143 0.0022 0.0015 0.0014
1000000.90.0606 0.0177 0.0022 0.0013 0.0015
10000000.10.0324 0.0153 0.0265 0.0121 0.0020
10000000.20.0819 0.0364 0.0262 0.0144 0.0040
10000000.30.1547 0.0774 0.0264 0.0166 0.0068
10000000.40.2447 0.1005 0.0264 0.0184 0.0089
10000000.50.3583 0.1280 0.0263 0.0192 0.0109
10000000.60.5024 0.1990 0.0265 0.0187 0.0136
10000000.70.6924 0.2338 0.0267 0.0170 0.0155
10000000.80.9639 0.2818 0.0268 0.0150 0.0174
10000000.91.4319 0.3625 0.0264 0.0129 0.0191
100000000.10.7067 0.3181 0.6920 0.1222 0.0514
100000000.21.7642 0.7155 0.6909 0.1451 0.1022
100000000.33.0718 1.2310 0.6921 0.1675 0.1537
100000000.44.6572 1.5806 0.6909 0.1860 0.2010
100000000.56.5762 2.3766 0.6904 0.1943 0.2535
100000000.68.9766 2.7640 0.6929 0.1887 0.2976
100000000.712.1109 3.2506 0.6924 0.1729 0.3389
100000000.816.5619 3.9203 0.6911 0.1529 0.3756
100000000.924.2751 5.9600 0.6909 0.1368 0.4156
1000000000.111.9186 5.0229 9.4275 1.2284 0.6073
1000000000.228.4891 10.7473 9.4202 1.4638 1.2054
1000000000.348.9235 14.4396 9.4179 1.6919 1.7780
1000000000.473.7326 23.2700 9.4068 1.8807 2.3940
1000000000.5103.9370 27.3172 9.4198 1.9683 2.9615
1000000000.6141.7940 32.3755 9.4193 1.9175 3.5221
1000000000.7191.4050 38.9861 9.4176 1.7639 4.1683
1000000000.8262.2190 57.5724 9.4168 1.5684 4.7086
1000000000.9384.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倍の速さですね!

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証