シミュレーション、統計推定、および再サンプリングのために乱数配列をすばやく生成します

現在、 モンテカルロ法を使用したシミュレーションは、外部から受信したデータの分析に基づいて複数の決定が必要とされる運用活動のほぼすべての分野で使用されています。 同時に、スペシャリストが解決する実際の問題の特徴を抽象メソッドに与えるために使用される乱数ジェネレーターの品質、パフォーマンス、および可用性が重要な役割を果たし始めています。 私が最近発見したように、この問題は並列プログラミングへの移行において決定的な役割を果たし始めます...あなたもこの問題に遭遇し、Windowsで必要な分布を持つ乱数の配列をすばやく取得する方法を知りたいですか？



/ *この記事ではかなりの数の英語の引用を使用しています。読者に十分なレベルの語学トレーニングを行っていると思います* /



私の勤務サービスアプリケーションの1つは、データベースからの要因の長いシーケンスでシミュレーションモデルを実行します。 詳細に入ることなく、将来の結果を予測する際の不確実性を減らすために、実際の実装の履歴からランダム性の強い要素を持つ繰り返しバイナリ（2ソース）プロセスの確率的特性を研究するというタスクが最近発生しました。 プロセスの実装は互いに独立しており、均等に（通常）分散していると見なされ、プロセスは固定されていると見なされます。 目標-可能な限りN個の連続したプロセス実装の結果の最も正確な予測-は、長さNの期間にわたる肯定的なプロセス実装の平均および最低予想頻度の評価を伴います。



ここで、確率変数の実現がその標準偏差のkを超える値によって期待値から逸脱しない信頼レベルの計算で標準統計的アプローチを使用できます。 このアプローチは、確率への周波数収束の概念に基づいています（ ベルヌーイの定理 ）。 もちろん、データベースに記録された数十万のプロセス実装に基づいて平均値を計算するのは非常に魅力的ですが、サブサンプルNの長さ（数十のオーダー）に制限があるため、制限されたサンプルでのプロセスのランダムな性質は、平均。 だから、与えられたNに対して、予想される平均周波数からの可能な最大の偏差が「定められる」ことを、一定の確実性で知る必要があるのです。



統計的なブートストラップ （ブートストラップ）の方法で必要な推定値を取得することにしました。これは、主な例のセットからモンテカルロ法を使用して、N個の実現からサンプルを複数生成します。 この方法の本質は、元の要素セットのランダムな組み合わせで構成される利用可能なサンプルから十分に多数の擬似サンプルを形成することです（その結果、一部の初期要素は1つの擬似サンプルで数回発生する場合がありますが、他の完全に存在しない場合もあります）ばらつきを調べるために、分析された統計特性の値を決定します。



予備的な推定によると、数千のパラメーターサブグループのそれぞれに対して、最大1,000万個の擬似サンプルを生成する必要があります。 データ量やその他の要因を考えると、数十時間のサーバーコンピューティング時間について話しています。 Visual Basic 6では、プロトタイプアプリケーションシミュレータが作成され、テストされました。 ただし、マルチスレッドの戦闘バージョンをコンパイルした後、複数のスレッドでシミュレーションを起動したときに、CPUコアあたりの負荷が予想される100％ではなく5％であることが突然明らかになりました。



Intel Vtune Amplifierプロファイラーの基本的なホットスポット分析では、異常なことは何も示されませんでした。 コードには算術演算とランダムインデックスによるメモリアクセス以外のものが含まれていないため、シミュレータはすべてのスレッドに作業を完全にロードする必要があると確信しました。

For k = 1 To NumIters Randomize For i = 1 To NumPockets NumFired = 0 For j = 1 To chainLength ind = d + Int(B * Rnd) NumFired = NumFired + FiredOrNot(ind) Next j Tmp = NumFired / chainLength If Tmp < MinProb Then MinProb = Tmp If Tmp > MaxProb Then MaxProb = Tmp Values(i) = Tmp TotalFired = TotalFired + NumFired Next i AvgProb = TotalFired / (k * NumPockets) / chainLength Next k
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

このコードはメモリ帯域幅に非常に結びついていますか？！ そして、なぜ数メガバイトのプロセッサキャッシュが必要なのか、なぜ対応できないのか！ 問題は数日間私を困惑させました。 遅れて開始された同時実行性分析は、VBシステムライブラリの腸内のどこかに大きなドローダウンを明らかにしましたが、これらの高価なシステム関数がどこから呼び出されるかを示すことができませんでした。







最後に、実験的に、問題はRnd演算子にあることがわかりました！ どうやら、Rndの並列使用は開発者によって提供されていません。 おそらく、Microsoftの実装では、何らかの種類の内部保留ロックを使用しています。 実際、アプリケーションを高速化する代わりにRndを並列化すると、何千もの速度が低下しました。 同様の動作は、VBAを使用するMS Office製品では一般的だと思います。 ランダムな整数インデックスを迅速に生成するという問題がこれまでにないほど重要になりました！ 乱数の配列を生成するための既製のライブラリを探すことにしました。



統計関数パッケージを備えたインテル®マスカーネルライブラリが選択されました。





MKLはテストモードで使用できます。 現在（バージョン11.02）、11の基本的なランダム変数ジェネレーター、12の連続、11の離散出力分布が利用可能です。 私にとって非常に重要であることが判明したのは、開発者がライブラリのスレッドセーフを保証し、さらにマルチスレッドモードで使用した場合の製品の利点を指摘していることです... MKLのドキュメントでは、次のセクションに注意が向けられました。

非決定的乱数ジェネレーター（RDRANDベースのジェネレーターのみ）[AVX]、[IntelSWMan]。 基になるハードウェアでサポートされている場合にのみ、非決定的乱数ジェネレーターを使用できます。 Intel CPUが非決定的乱数ジェネレーターをサポートしているかどうかを検出する方法については、たとえば、[AVX]の第8章：32nm後のプロセッサ命令または[BMT]の第4章：RdRand命令の使用を参照してください。

この「鉄」機能を試してみたかったのですが、残念ながら、私のプロセッサにはハードウェア乱数ジェネレータが装備されていないことがわかりました。 上記に加えて、MKLパッケージの重要な機能は、ほとんどのプログラミング言語で演算子呼び出しごとに1つの乱数値を取得する通常の慣行とは対照的に、乱数ストリームのブロック生成のパラダイムです。



Xeon Phiコプロセッサアーキテクチャを使用したシミュレーションにMKL VSLの乱数を使用した実用例、たとえばShuo Liのケーススタディ：Stepwise Optimization Frameworkを使用したMonte Carlo European Optionでの高性能の実現は、ネットワーク上ですでに利用可能です。 たぶん誰かがこの記事を翻訳したいですか？ ちなみに、乱数の連続生成から「ストリーミング」への移行中に、適切な速度の増加が記録されます。

R NGをC ++ TR1からインテルMKLに変更すると、元のコードが5.53倍改善されるだけでなく、コンパイラーがインライン関数呼び出しを実行し、内部ループでコードをベクトル化することも可能になりました...

MICアーキテクチャ用に再コンパイルすると、2つの8コアSandy Bridge @ 2.6 Hz Xeonsと比較してさらに大きく成長します。

ホストシステムで完了するのに44秒かかった同じプログラムが、わずか5秒未満で完了するようになりました。これは8.82倍の改善です。

しかし、記事から私はまだ著者がどのモデルPhiを使用したのか理解していませんでした、明らかに、4120ドルで7120X（16GB、1.238 GHz、61コア）。 ナイツランディングを待たずに、今までに時間がなかったファイを連れて行く時間でしょう。;-)戦闘で試して「ファイを言う」のは面白いでしょう。



いずれにせよ、乱数ブロックを生成するMKLのアプローチに興味があります。 まず、vslNewStreamを呼び出すと、指定された特性を持つ乱数ストリームジェネレーターが作成されます。たとえば、ベースジェネレーターの名前VSL_BRNG_SFMT19937（SIMD指向のFast Mersenne Twister擬似乱数ジェネレーター）は魅力的です。 次に、このジェネレーターを（および異なるストリームから）使用して、割り当てられたメモリブロックに、希望する分布とデータ型の乱数列を直接入力します。 たとえば、virnguniform（メソッド、ストリーム、n、r、a、b）などの呼び出しは、ストリームnジェネレーターに基づいて、均一に分散されたaからbまでのランダムな値で配列rを満たします。 作業後、vslDeleteStreamを呼び出してストリームを削除する必要があります。



擬似乱数ジェネレーターの有限期間により、乱数の並列使用/生成には重要な側面があります。これは、「 モンテカルロアプローチをさらに簡単かつ高速にする 」という記事で、同胞のセルゲイメイダノフとアンドレイナライキンが美しく描いたものです。 好奇心のために、抜粋を引用します（翻訳なし）：

MKLライブラリーでは、Intelの優秀なスタッフが、サブシーケンスに分割するための先読みアプローチと跳躍アプローチの両方を実装していることを言わなければなりません。 ただし、私のタスクの一部として、1つの基本的なジェネレーターからの複数のストリームが使用されますが、それらは独立してデータの異なるサブグループで使用されるため、互いからのストリームの絶対的な統計的独立性の要件はありません。 したがって、私にとっての主な技術的困難は、CおよびFortranによって強化されたMKLインターフェイスをVisual Basicプロジェクトに接続することです。 Microsoft * Office Excelのインテル®MKLを使用する記事では、MKL自体を使用してラッパーdllを生成する方法を説明し、 ExcelからMLK関数（もちろんstdcall規則）の呼び出しの例を示します。



概念は次のとおりです。MKLBuilderに、指定された呼び出し規約に従って、必要なMKL関数の任意のリストをエクスポートするdllを自動的に生成するように依頼できます。 もちろん、VC ++で独自のdllプロジェクトを作成できますが、自動生成はより普遍的な方法であり、必要なファイルをより高速に取得するのに役立つように思えました。 記事が古くなっていることがすぐに明らかになり、テストコマンドnmake libia32 threading =シーケンシャルインターフェイス= stdcall export = user_example_list name = mkl_customも多くの理由ですぐにクラッシュしました。 msvcrt.libファイルをmkl \ tools \ builder \フォルダーに追加して、問題の1つを解決しました。 複数のメッセージ、 未解決の外部シンボル___security_cookieを別々に楽しみました。



マイクロソフトの推奨事項と半日をかけての殺害を使用して、bufferoverflowU.libパラメーターを後にリンカーコマンドラインに割り当てることにより、ビルダーフォルダーのメイクファイルを変更する必要があるという結論に達しました：$（CB_NAME）.dll。 さて、スクラッチのあるメカニズムが機能し、 nmake libia32 threading = parallel interface = stdcall export = vml_vsl_stdcall_example_list name = mkl_customを発行することで、ベクター統計ライブラリのすべての機能を含む38メガバイトのダイナミックライブラリの幸運な所有者になりました。 私自身の健全な資本主義欲を克服し、ライブラリに必要な関数を3つだけ残しました。その後、サイズは数メガバイトに減少しました。



次に、呼び出しコードについて。 VB6およびVBAでは、必要な関数は次のように宣言する必要があります。

 Public Declare Function vslNewStream Lib "mkl_custom.dll" (ByRef StreamPtr As Long, ByVal brng As Long, ByVal seed As Long) As Long Public Declare Function vslDeleteStream Lib "mkl_custom.dll" (ByRef StreamPtr As Long) As Long Public Declare Function viRngUniform Lib "mkl_custom.dll" (ByVal method As Long, ByVal StreamPtr As Long, ByVal n As Long, ByRef Data As Long, ByVal A As Long, ByVal B As Long) As Long
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

テストとして、ベースジェネレーターVSL_BRNG_MCG31（31ビットの乗算合同ジェネレーター）を使用して、1〜10,000（10,000を含まない）の500万個のランダムな整数インデックスの配列を100回生成してみましょう。

 Dim i&, theStream&, v&(), res&, t! Dim theStream As Long ReDim v(1 To 5000000): res = vslNewStream(theStream, VSL_BRNG_MCG31, 777) t = Timer For i = 1 To 100 res = viRngUniform(VSL_RNG_METHOD_UNIFORM_STD, theStream, 5000000, v(1), 1, 10000) Next i Debug.Print Timer - t res = vslDeleteStream(theStream)
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

ランタイム-1.19秒。 ところで、SSEの使用方法を知っているVSL_BRNG_SFMT19937ジェネレーターは、0.99秒で同じタスクを実行できます。 さて、必要な操作の呼び出しをメインコードに追加し、プロセッサの負荷を監視します。







ビンゴ、フルダウンロード！ cuRANDパッケージを使用してGPGPUを実装するというアイデアがありましたが、その利点はNvidiaによって明確に示されています。



インテルMKLと比較した非常に高速なCURANDパフォーマンス





ただし、最初の目標はMKLで達成され、Cudaバージョンをテストする自由時間はまだありません。 :-)



要約すると、乱数の並列生成におけるいくつかの問題について学び、インテルMKLを使用してそれらを解決する方法に関する実用的な推奨事項を受け取りました。 この強力なライブラリは、最新のVisual Studioの任意の言語（C ++、C＃、VB、およびVisual Fortran）から使用できます。また、音声技術を使用して、外部DLL（VB6、VBA、LabView、Matlabなど）を呼び出すことができる環境からも使用できますその他。



ご清聴ありがとうございました！