インテル® MKL クックブック

置換のない複数の単純なランダム・サンプリング

目的

サイズ N (1 ≤M≤N) の母集団から置換なしで K>>1 の単純なランダム長 M のサンプルを生成する。

ソリューション

問題の定義と詳細は、[SRSWOR] を参照してください。

部分 Fisher-Yates シャッフル・アルゴリズム [KnuthV2] の次の実装とインテル® MKL 乱数ジェネレーター (RNG) を使用して各サンプルを生成します。

部分 Fisher-Yates シャッフル・アルゴリズム

 A2.1: (初期化ステップ) PERMUT_BUF は自然数 1, 2, ..., N を含むものとする

 A2.2: for i from 1 to M do:

     A2.3: {i,...,N} で一様なランダムな整数を生成

     A2.4: PERMUT_BUF[i] と PERMUT_BUF[X] を交換

 A2.5: (コピーステップ) for i from 1 to M do: RESULTS_ARRAY[i]=PERMUT_BUF[i]

 End.

アルゴリズムを実装するプログラムは 11 969 664 の実験を行います。各実験 (1 から N までの M の一様なランダム自然数のシーケンスを生成) は、実際には N 要素の母集団から部分長 M をランダムにシャッフルします。 アルゴリズムのメインループは実際に抽選を行うため、プログラムで各実験は「N から M の抽選」と呼ばれています。

プログラムは M=6 および N=49 を使用して、結果のサンプル (長さ M のシーケンス) を単一配列 RESULTS_ARRAY に格納し、利用可能な並列スレッドをすべて使用します。

ソースコード: サンプル (http://software.intel.com/en-us/mkl_cookbook_samples) の lottery6of49 フォルダーを参照してください。

並列化

#pragma omp parallel
{
    thr = omp_get_thread_num(); /* スレッドのインデックス */
    VSLStreamStatePtr stream;
    /* このスレッドの RNG ストリームを初期化 */
    vslNewStream( &stream, VSL_BRNG_MT2203+thr, seed );
    ... /* (スレッド番号 thr で) 実験サンプルを生成 */
    vslDeleteStream( &stream );
}

コードは、OpenMP* #pragma parallel 宣言子を使用して、すべての CPU のすべての利用可能なプロセッサー・コアを利用します。 実験結果の配列 RESULTS_ARRAY は THREADS_NUM の部分に分割されます。ここで、THREADS_NUM は利用可能なスレッド数で、各スレッド (並列領域) は配列の個別の部分を処理します。

インテル® MKL 基本乱数ジェネレーターは、VSL_BRNG_MT2203 パラメーターを指定することにより、各スレッドで並列の独立したストリームをサポートします。

実験サンプルの生成

/* A2.1: 初期化ステップ */
/*  PERMUT_BUF は自然数 1, 2, ..., N を含むものとする */
 for( i=0; i<N; i++ ) PERMUT_BUF[i]=i+1; /* セット {1,...,N} を使用 */
 for( sample_num=0; sample_num<EXPERIM_NUM/THREADS_NUM; sample_num++ ){
     /* 次の抽選サンプルを生成 (ステップ A2.2、A2.3、A2.4): */
     Fisher_Yates_shuffle(...);
     /*  A2.5: コピーステップ */   
     for(i=0; i<M; i++)
         RESULTS_ARRAY[thr*ONE_THR_PORTION_SIZE + sample_num*M + i] = PERMUT_BUF[i];
 }

このコードは各スレッドで 部分 Fisher-Yates シャッフル・アルゴリズムを実装します。

多くの実験をシミュレーションする場合、初期化ステップは各実験の最初に 1 回のみ必要です。(実際の抽選のように) PERMUT_BUF 配列の自然数 1...N の順序は重要ではありません。

Fisher_Yates_shuffle 関数

/* A2.2: for i from 0 to M-1 do */
Fisher_Yates_shuffle (...)
{  
    for(i=0; i<M; i++) {
        /* A2.3: {i,...,N-1} からランダム自然数 X を生成 */
        j = Next_Uniform_Int(...);
        /* A2.4: PERMUT_BUF[i] と PERMUT_BUF[X] を交換 */
        tmp = PERMUT_BUF[i];
        PERMUT_BUF[i] = PERMUT_BUF[j];
        PERMUT_BUF[j] = tmp;
    }
}

ループ A2.2 の各反復は、実際に抽選を行います。残りのアイテム PERMUT_BUF[i], ..., PERMUT_BUF[N] を含む箱からランダムなアイテム X を取り出し、アイテム X を結果行 PERMUT_BUF[1],...,PERMUT_BUF[i] の最後に追加します。 長さ N の完全な置換ではなく、一部の長さ M のみを生成するため、アルゴリズムは部分的です。

注

アルゴリズムで説明した疑似コードと異なり、プログラムはゼロベースの配列を使用します。

説明

ステップ A2.3 で、プログラムは Next_Uniform_Int 関数を呼び出して次のランダム整数 X を生成し、{i, ..., N-1} で一様にします (詳細はソースコードを参照)。 インテル® MKL のベクトル化された RNG の能力を引き出しつつ、ベクトル化のオーバーヘッドを最小限に抑えるため、ジェネレーターは L1 キャッシュに収まるサイズ RNGBUFSIZE のベクトル D_UNIFORM01_BUF を生成する必要があります。 各スレッドは、独自のバッファー D_UNIFORM01_BUF およびそのバッファーで最後に使用された乱数の後を指すインデックス D_UNIFORM01_IDX を使用します。 Next_Uniform_Int 関数の最初の呼び出しでは (あるいはバッファーの乱数がすべて使用された場合は)、長さ RNGBUFSIZE とインデックス D_UNIFORM01_IDX をゼロにセットして vdRngUniform 関数を呼び出し、乱数バッファーを生成し直します。 

vdRngUniform( ... RNGBUFSIZE, D_UNIFORM01_BUF ... );

インテル® MKL は同じ分布の乱数値ジェネレーターのみ提供しますが、ステップ A2.3 では異なる範囲のランダム整数が必要なため、[0;1) で一様に分布された倍精度乱数をバッファーに格納した後、「整数のスケーリング」ステップで、これらの倍精度値を必要な整数範囲に収まるように変換します。 

number 0   distributed on {0,...,N-1}   = 0   + {0,...,N-1}
number 1   distributed on {1,...,N-1}   = 1   + {0,...,N-2}

...

number M-1 distributed on {M-1,...,N-1} = M-1 + {0,...,N-M}

(前の M ステップを繰り返す)

number M     distributed on: see (0)
number M+1   distributed on: see (1)

...

number 2*M-1 distributed on: see (M-1)
(前の M ステップを再度繰り返す)
...

続く

整数スケーリング

/* 整数スケーリング・ステップ */
for(i=0;i<RNGBUFSIZE/M;i++)
    for(k=0;k<M;k++)
        I_RNG_BUF[i*M+k] =
            k + (unsigned int)(D_UNIFORM01_BUF[i*M+k] * (double)(N-k));

ここで RNGBUFSIZE は M の倍数。

このコードに関するパフォーマンスの注意点は [SRSWOR] を参照してください。

使用するルーチン

タスク

ルーチン

RNG ストリームを作成して初期化する

vslNewStream

区間 [0;1) に一様に分散された倍精度数を生成する

vdRngUniform

RNG ストリームを削除する

vslDeleteStream

64 ビット境界でアライメントされたメモリーバッファーを割り当てる

mkl_malloc

mkl_malloc で割り当てられたメモリーを解放する

mkl_free

最適化に関する注意事項

インテル® コンパイラーは、互換マイクロプロセッサー向けには、インテル製マイクロプロセッサー向けと同等レベルの最適化が行われない可能性があります。これには、インテル® ストリーミング SIMD 拡張命令 2 (インテル® SSE2)、インテル® ストリーミング SIMD 拡張命令 3 (インテル® SSE3)、ストリーミング SIMD 拡張命令 3 補足命令 (SSSE3) 命令セットに関連する最適化およびその他の最適化が含まれます。インテルでは、インテル製ではないマイクロプロセッサーに対して、最適化の提供、機能、効果を保証していません。本製品のマイクロプロセッサー固有の最適化は、インテル製マイクロプロセッサーでの使用を目的としています。インテル® マイクロアーキテクチャーに非固有の特定の最適化は、インテル製マイクロプロセッサー向けに予約されています。この注意事項の適用対象である特定の命令セットの詳細は、該当する製品のユーザー・リファレンス・ガイドを参照してください。

改訂 #20110804

親トピック: インテル® MKL クックブック