私たちは誰でも、一日に何度も「 」問題を解きます。もちろん、意図せずに。仕事の途中で、インターネットを閲覧しているとき、または机の上で物を整理しているときに、この問題を解きます。 最短経路 ちょっと…やりすぎのように聞こえますか？調べてみましょう。 想像してください。あなたは友達と、まあ、カフェで会うことに決めたとしましょう。まず、カフェまでのルート（または経路）を見つける必要があるので、利用可能な公共交通機関（徒歩の場合）またはルートと道路（車の場合）を探し始めます。現在地からカフェまでのルートを探していますが、「任意の」ルート（最短または最速のルート）を探しているわけではありません。 もう一つの例を挙げます。これは前の例ほど明白ではありません。経路探索中に、脇道を通る「近道」を取ることにしました。なぜなら、それは「近道」であり、この道を通る方が「速い」からです。しかし、この「近道」の方が速いことをどうやって知ったのでしょうか。個人的な経験に基づいて「最短経路」問題を解き、脇道を通るルートを選択しました。 どちらの例でも、「最短」の経路は、ある場所から別の場所へ移動するために必要な距離または時間で決定されます。移動の例は、 、特に「最短経路」問題の非常に自然な応用です。しかし、机の上に物を並べる例はどうでしょうか。この場合、「最短」は、たとえば 1 枚の紙を手に入れるために実行する必要があるアクションの数または複雑さを表します。机を開け、フォルダーを開き、1 枚の紙を取り出すか、机から直接 1 枚の紙を取り出すかです。 グラフ理論 上記の例はすべて、グラフ内の 2 つの頂点間の最短経路 (2 つの場所間のルートまたはパス、ある場所から別の場所へ 1 枚の紙を取得するアクションの数または複雑さ) を見つける問題を表しています。このクラスの最短経路問題は と呼ばれ、これらの問題を解決するための基本アルゴリズムは、計算複雑度が です。 SSSP (Single Source Shortest Path) O(n^2) の Dijkstra アルゴリズム しかし、すべての頂点間の最短経路をすべて見つける必要がある場合もあります。次の例を考えてみましょう。 、 、 間の定期的な移動のマップを作成しているとします。このシナリオでは、 、 、 、 、 、 の 6 つのルートが作成されます (逆のルートは、一方通行などの理由で異なる場合があります)。 自宅 職場 劇場 work ⭢ home home ⭢ work work ⭢ theatre theatre ⭢ work home ⭢ theatre theatre ⭢ home マップに場所を追加すると、組み合わせが大幅に増加します。 論の によれば、次のように計算できます。 組み合わせ n の順列の公式 A(k, n) = n! / (n - m)! // where n - is a number of elements, // k - is a number of elements in permutation (in our case k = 2) これにより、4 つの場所に 12 のルート、10 の場所には 90 のルートが得られます (これは印象的です)。注: これは、場所間の中間点を考慮していない場合です。つまり、自宅から職場に行くには、4 つの道路を渡り、川に沿って歩き、橋を渡る必要があります。いくつかのルートに共通の中間点があると想像してください...そうですね...結果として、多数の頂点を持つ非常に大きなグラフが得られ、各頂点は場所または重要な中間点のいずれかを表します。グラフ内のすべての頂点ペア間の最短経路をすべて見つける必要がある問題のクラスは、 と呼ばれ、これらの問題を解決するための基本アルゴリズムは、 の計算複雑度を持つ です。 APSP (全ペア最短経路) O(n^3) Floyd-Warshall アルゴリズム これが今日実装するアルゴリズムです 🙂 フロイド・ワーシャルアルゴリズム フロイド・ワーシャルアルゴリズムは、グラフ内の頂点のペア間の最短経路をすべて見つけます。このアルゴリズムは、ロバート・フロイドによって [1] で公開されました (詳細については、「参考文献」セクションを参照してください)。同じ年に、ピーター・インガーマンは [2] で、3 つのネストされた ループの形式でアルゴリズムの最新の実装を説明しました。 for algorithm FloydWarshall(W) do for k = 0 to N - 1 do for i = 0 to N - 1 do for j = 0 to N - 1 do W[i, j] = min(W[i, j], W[i, k] + W[k, j]) end for end for end for end algorithm // where W - is a weight matrix of N x N size, // min() - is a function which returns lesser of it's arguments 行列形式で表現されたグラフを扱ったことがない場合は、上記のアルゴリズムが何を行うのか理解するのが難しいかもしれません。そこで、理解を深めるために、グラフを行列形式で表現する方法と、そのような表現が最短経路問題を解決するのになぜ役立つのかを見てみましょう。 下の図は、5 つの頂点を持つ グラフを示しています。左側のグラフは、円と矢印で構成された視覚的な形式で表示されています。各円は頂点を表し、矢印は方向のあるエッジを表します。円内の数字は頂点番号に対応し、エッジの上の数字はエッジの重みに対応します。右側では、同じグラフが重みマトリックスの形式で表示されています。重みマトリックスは の一種で、各マトリックス セルには「重み」、つまり頂点 (行) と頂点 (列) の間の距離が含まれます。重みマトリックスには、頂点間の「パス」に関する情報 ( から に移動するときに経由する頂点のリスト) は含まれず、これらの頂点間の重み (距離) のみが含まれます。 有向の重み付き 隣接マトリックス i j i j 重みマトリックスでは、セルの値が 頂点間の重みに等しいことがわかります。そのため、頂点 (行 ) からのパスを調べると、 から へのパスが 1 つしかないことがわかります。ただし、視覚的に表現すると、頂点 から頂点 と (頂点 経由) へのパスが明確にわかります。この理由は単純です。初期状態では、重みマトリックスには隣接する頂点間の距離のみが含まれます。ただし、この情報だけで残りを のに十分です。 隣接する 0 0 0 1 0 2 3 1 見つける どのように動作するか見てみましょう。セル に注目してください。その値は、頂点 から頂点 への重みが のパスがあることを示しています (簡単に言うと、 と記述できます)。この知識があれば、行 のすべてのセル (頂点 からのすべてのパスのすべての重みを含む) をスキャンし、この情報を行 にバックポートして、重みを ( の値) 増やすことができます。 W[0, 1] 0 1 1 0 ⭢ 1: 1 1 1 0 1 W[0, 1] 同じ手順を使用して、頂点 から他の頂点までのパスを見つけることができます。検索中に、同じ頂点につながるパスが複数あることがあり、さらに重要なことに、これらのパスの重みが異なる可能性があります。このような状況の例を下の図に示します。頂点 から頂点 までの検索で、重みの小さい頂点 へのパスがもう 1 つ見つかりました。 0 0 2 3 パスは 2 つあります。元のパス と、先ほど発見した新しいパス (重み行列にはパスが含まれていないため、元のパスにどの頂点が含まれているかはわかりません)。ここで、最短のパスを保持し、セル に 書き込むことを決定しました。 0 ⭢ 3: 4 0 ⭢ 2 ⭢ 3: 3 W[0, 3] 3 最初の最短経路を見つけたようです! 先ほど見た手順は、フロイド・ワーシャル アルゴリズムの本質です。もう一度、アルゴリズムの疑似コードを見てみましょう。 algorithm FloydWarshall(W) do for k = 0 to N - 1 do for i = 0 to N - 1 do for j = 0 to N - 1 do W[i, j] = min(W[i, j], W[i, k] + W[k, j]) end for end for end for end algorithm 上の 最も外側のサイクルでは、グラフ内のすべての頂点が反復処理され、各反復で、変数 頂点を表します。 上の 最も内側のサイクルでも、グラフ内のすべての頂点が反復処理され、各反復で、 頂点を表します。最後に、 上の 最も内側のサイクルでは、グラフ内のすべての頂点が反復処理され、各反復で、 これらを組み合わせると、次のようになります。各反復 で、すべての頂点 から頂点 を経由してすべての頂点 に至るパスを検索します。サイクル内では、パス ( で表される) とパス ( と の合計で表される) を比較し、最短のパスを に書き戻します。 k for k パスを検索している i for i パスを検索している j for j パスを検索している頂点を表します。 k i k j i ⭢ j W[i, j] i ⭢ k ⭢ j W[I, k] W[k, j] W[i, j] さて、仕組みを理解したら、アルゴリズムを実装する時です。 基本的な実装 ソース コードと実験グラフは、GitHub の で入手できます。実験グラフは、 ディレクトリにあります。この記事のすべてのベンチマークは、 ファイルに実装されています。 リポジトリ Data/Sparse-Graphs.zip APSP01.cs 実装に入る前に、いくつかの技術的な点を明確にする必要があります。 すべての実装は、線形配列の形式で表される重み行列で動作します。 すべての実装では整数演算を使用します。頂点間のパスが存在しないことは、特別な定数重み で表されます。 NO_EDGE = (int.MaxValue / 2) - 1 さて、その理由を考えてみましょう。 マトリックスについて話すとき、セルを「行」と「列」で定義します。このため、マトリックスを「正方形」または「長方形」の形で想像するのが自然なようです (図 4a)。 1 に関して。 しかし、これは必ずしも、想像力を働かせるために、行列を配列の配列の形式で表現しなければならないという意味ではありません (図 4b)。その代わりに、各セルのインデックスが次の式を使用して計算される線形配列の形式で行列を表現することができます (図 4c)。 i = row * row_size + col; // where row - cell row index, // col - cell column index, // row_size - number of cells in a row. 重み行列の線形配列は、フロイド・ワーシャル アルゴリズムを効果的に実行するための です。その理由は単純ではなく、詳細な説明には別の投稿、またはいくつかの投稿が必要です。ただし、現時点では、このような表現によって が大幅に向上し、実際にアルゴリズムのパフォーマンスに大きな影響を与えることに注意してください。 前提条件 データの局所性 ここで、前提条件として、この情報だけを念頭に置いて私を信じてほしいとお願いしていますが、同時に、時間をかけてこの問題について研究することをお勧めします。ちなみに、インターネット上の人々を信じないでください。 - 著者注 アルゴリズムの擬似コードをよく見ると、2 つの頂点間のパスの存在に関連するチェックは見つかりません。代わりに、擬似コードでは 関数が使用されています。理由は簡単です。元々、2 つの頂点間にパスがない場合、セルの値は に設定され、JavaScript を除くすべての言語では、すべての値は 未満です。整数の場合、「パスなし」の値として 使用したくなるかもしれません。ただし、 パスと パスの両方の値が に等しい場合、整数オーバーフローが発生します。そのため、 の半分より 1 小さい値を使用します。 2 について。 min() infinity infinity int.MaxValue i ⭢ k k ⭢ j int.MaxValue int.MaxValue おい！でも、計算を行う前にパスが存在するかどうかをチェックできないのはなぜでしょうか。たとえば、両方のパスを 0 と比較する（ゼロを「パスなし」の値と見なす場合）。 - 思慮深い読者 確かに可能ですが、残念ながらパフォーマンスが大幅に低下します。簡単に言うと、CPU は分岐評価結果の統計を保持します。たとえば、 文の一部が または に評価された場合などです。この統計を使用して、実際の 文が評価される前に「統計的に予測された分岐」のコードを事前に実行します (これを 呼びます)。これにより、コードの実行がより効率的になります。ただし、CPU の予測が不正確な場合、CPU は停止して正しい分岐を計算する必要があるため、正しい予測と無条件実行に比べてパフォーマンスが大幅に低下します。 if true false if 投機的実行と 各反復 で重みマトリックスの重要な部分を更新するため、コード パターンがなく、すべての分岐がデータのみに基づいているため、CPU 分岐統計は役に立たなくなります。そのため、このようなチェックでは、 が大量に発生します。 k 分岐予測の誤り ここで私は、（今のところは）私を信じていただき、その後、時間をかけてこのテーマを研究していただきたいと考えています。 うーん、理論的な部分は終わったようです。アルゴリズムを実装しましょう (この実装を と呼びます)。 Baseline public void Baseline(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { for (var i = 0; i < sz; ++i) { for (var j = 0; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } } } } 上記のコードは、1 つの例外を除いて、前述の疑似コードとほぼ同じコピーです。Math.Min の代わりに、必要な場合にのみセルを更新するために を使用します。 Math.Min() if おい！ちょっと待って、ここで if が良くない理由について長々と書いたのに、数行後には私たち自身が if を導入したのはあなたじゃなかったのか？ - 思慮深い読者 理由は簡単です。執筆時点では、JIT は と 両方の実装に対してほぼ同等のコードを出力します。詳細は で確認できますが、メイン ループ本体のスニペットを以下に示します。 if Math.Min sharplab.io // // == Assembly code of implementation of innermost loop for of j using if. // 53: movsxd r14, r14d // // compare matrix[i * sz + j] and distance (Condition) // 56: cmp [rdx+r14*4+0x10], ebx 5b: jle short 64 // // if matrix[i * sz + j] greater than distance write distance to matrix // 5d: movsxd rbp, ebp 60: mov [rdx+rbp*4+0x10], ebx 64: // end of loop // // == Assembly code of implementation of innermost loop for of j using Math.Min. // 4f: movsxd rbp, ebp 52: mov r14d, [rdx+rbp*4+0x10] // // compare matrix[i * sz + j] and distance (Condition) // 57: cmp r14d, ebx. // // if matrix[i * sz + j] less than distance write value to matrix // 5a: jle short 5e // // otherwise write distance to matrix // 5c: jmp short 61 5e: mov ebx, r14d 61: mov [rdx+rbp*4+0x10], ebx 65: // end of loop と どちらを使用しても、条件チェックが行われていることがわかります。ただし、 の場合は不要な書き込みはありません。 if Math.Min if 実装が完了したら、コードを実行して速度を確認してみましょう。 でテストを実行することで、コードの正確性を自分で確認できます。 リポジトリ 私は (バージョン 0.12.1) を使用してコードをベンチマークします。ベンチマークで使用されるすべてのグラフは グラフです。グラフ内のエッジの数は、有向非巡回グラフの場合、次のように計算できる最大値の約 80% です。 Benchmark.NET 、300、600、1200、2400、および 4800 頂点の有向非巡回 var max = v * (v - 1)) / 2; // where v - is a number of vertexes in a graph. ロックしましょう！ 以下は私の PC (Windows 10.0.19042、Intel Core i7-7700 CPU 3.60Ghz (Kaby Lake) / 8 論理プロセッサ / 4 コア) で実行したベンチマークの結果です。 方法 サイズ 平均 エラー 標準偏差 ベースライン 300 27.525 ミリ秒 0.1937 ミリ秒 0.1617ミリ秒 ベースライン 600 217.897 ミリ秒 1.6415 ミリ秒 1.5355 ミリ秒 ベースライン 1200 1,763.335 ミリ秒 7.4561 ミリ秒 6.2262 ミリ秒 ベースライン 2400 14,533.335 ミリ秒 63.3518 ミリ秒 52.9016 ミリ秒 ベースライン 4800 119,768.219 ミリ秒 181.5514 ミリ秒 160.9406 ミリ秒 結果からわかるように、計算時間はグラフのサイズに比べて大幅に増加します。300 頂点のグラフでは 27 ミリ秒、2400 頂点のグラフでは 14.5 秒、4800 頂点のグラフでは 119 秒、つまり かかります。 ほぼ 2 分 アルゴリズムのコードを見ると、計算を高速化するために何かできることがあるとは想像しにくいかもしれません... なぜなら、ループが 3 つあるからです。ループは 3 つだけです。 しかし、よくあることですが、可能性は細部に隠されています 🙂 データを知る – 疎グラフ フロイド・ワーシャル アルゴリズムは、特に グラフや グラフの場合に、すべてのペアの最短経路問題を解決するための基本アルゴリズムです (アルゴリズムはすべての頂点ペア間の経路を検索するため)。 密な 完全な ただし、私たちの実験では グラフを使用します。このグラフには、頂点 から頂点 へのパスがある場合、頂点 から頂点 へのパスは存在しないという素晴らしい特性があります。私たちにとって、これは、 から へのパスがない場合にスキップできる隣接していない頂点が多数あることを意味します (この実装を と呼びます)。 、有向非巡回 1 2 2 1 i k SpartialOptimisation public void SpartialOptimisation(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { for (var i = 0; i < sz; ++i) { if (matrix[i * sz + k] == NO_EDGE) { continue; } for (var j = 0; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } } } } 以下は、同じグラフ セットに対する以前の実装 ( ) と現在の実装 ( ) の結果です (最速の結果は で強調表示されています)。 Baseline SpartialOptimisation 太字 方法 サイズ 平均 エラー 標準偏差 比率 ベースライン 300 27.525 ミリ秒 0.1937ミリ秒 0.1617ミリ秒 1.00 部分最適化 300 12.399ミリ秒 0.0943 ミリ秒 0.0882 ミリ秒 0.45 ベースライン 600 217.897 ミリ秒 1.6415ミリ秒 1.5355 ミリ秒 1.00 部分最適化 600 99.122 ミリ秒 0.8230 ミリ秒 0.7698 ミリ秒 0.45 ベースライン 1200 1,763.335 ミリ秒 7.4561 ミリ秒 6.2262 ミリ秒 1.00 部分最適化 1200 766.675 ミリ秒 6.1147 ミリ秒 5.7197 ミリ秒 0.43 ベースライン 2400 14,533.335 ミリ秒 63.3518 ミリ秒 52.9016 ミリ秒 1.00 部分最適化 2400 6,507.878 ミリ秒 28.2317 ミリ秒 26.4079 ミリ秒 0.45 ベースライン 4800 119,768.219 ミリ秒 181.5514 ミリ秒 160.9406 ミリ秒 1.00 部分最適化 4800 55,590.374 ミリ秒 414.6051 ミリ秒 387.8218 ミリ秒 0.46 かなり印象的ですね！ 計算時間を半分に短縮しました。もちろん、グラフの密度が高くなるほど、速度の向上は少なくなります。ただし、これは、処理する予定のデータのクラスが事前にわかっている場合に便利な最適化の 1 つです。 もっとできるかな？ 🙂 ハードウェアを理解する – 並列処理 私の PC には、8 つの論理プロセッサ ( ) または テクノロジを備えた 4 つのコアを備えた プロセッサが搭載されています。コアが複数あるということは、作業に使える「手」が増えるようなものです。では、作業のどの部分を複数のワーカー間で効率的に分割できるかを確認し、それを並列化してみましょう。 HW ハイパースレッディング Inter Core i7-7700 CPU 3.60GHz (Kaby Lake) ループは、ほとんどの場合、反復が独立しており、同時に実行できるため、常に並列化の最も明白な候補です。アルゴリズムには、分析して並列化を妨げる依存関係があるかどうかを確認する必要がある 3 つのループがあります。 ループから始めましょう。各反復で、ループは頂点 を通るすべての頂点からすべての頂点へのパスを計算します。新しいパスと更新されたパスは重み行列に格納されます。反復を見ると、結果に影響を与えずに 0、1、2、3 または 3、1、2、0 の任意の順序で実行できることに気付くかもしれません。ただし、それらは順番に実行する必要があります。そうしないと、一部の反復では、新しいパスまたは更新されたパスが重み行列にまだ書き込まれていないため、それらを使用できなくなります。このような不一致は、間違いなく結果を台無しにします。そのため、引き続き調べる必要があります。 for of k k 次の候補は ループです。各反復で、ループは頂点 から頂点 へのパス (セル: )、頂点 から頂点 へのパス (セル: ]) を読み取り、次に から への既知のパス (セル: ) が パス (合計: + ) より短いかどうかを確認します。アクセス パターンを詳しく見ると、各反復で ループが 行からデータを読み取り、 行を更新していることに気付くかもしれません。これは基本的に、反復が独立しており、任意の順序で実行できるだけでなく、並列に実行できることを意味しています。 for of i i k W[i, k] k j W[k, j i j W[i, j] i ⭢ k ⭢ j W[i, k] W[k, j] i k i これは有望そうなので、実装してみましょう (この実装を と呼びます)。 SpartialParallelOptimisations も並列化できます。ただし、この場合、最も内側のサイクルの並列化は非常に非効率的です。 にいくつかの簡単な変更を加えることで、自分で検証できます。 for of j ソース コード - 著者注 public void SpartialParallelOptimisations(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { Parallel.For(0, sz, i => { if (matrix[i * sz + k] == NO_EDGE) { return; } for (var j = 0; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } }); } } 以下は、同じグラフ セットに対する 、 、および 実装の結果です (並列化は クラスを使用して行われます)。 Baseline SpartialOptimisation SpartialParallelOptimisations Parallel 方法 サイズ 平均 エラー 標準偏差 比率 ベースライン 300 27.525 ミリ秒 0.1937ミリ秒 0.1617ミリ秒 1.00 部分最適化 300 12.399ミリ秒 0.0943 ミリ秒 0.0882 ミリ秒 0.45 部分並列最適化 300 6.252 ミリ秒 0.0211ミリ秒 0.0187ミリ秒 0.23 ベースライン 600 217.897 ミリ秒 1.6415 ミリ秒 1.5355 ミリ秒 1.00 部分最適化 600 99.122 ミリ秒 0.8230 ミリ秒 0.7698 ミリ秒 0.45 部分並列最適化 600 35.825 ミリ秒 0.1003 ミリ秒 0.0837 ミリ秒 0.16 ベースライン 1200 1,763.335 ミリ秒 7.4561 ミリ秒 6.2262 ミリ秒 1.00 部分最適化 1200 766.675 ミリ秒 6.1147 ミリ秒 5.7197 ミリ秒 0.43 部分並列最適化 1200 248.801 ミリ秒 0.6040 ミリ秒 0.5043 ミリ秒 0.14 ベースライン 2400 14,533.335 ミリ秒 63.3518 ミリ秒 52.9016 ミリ秒 1.00 部分最適化 2400 6,507.878 ミリ秒 28.2317 ミリ秒 26.4079 ミリ秒 0.45 部分並列最適化 2400 2,076.403 ミリ秒 20.8320 ミリ秒 19.4863 ミリ秒 0.14 ベースライン 4800 119,768.219 ミリ秒 181.5514 ミリ秒 160.9406 ミリ秒 1.00 部分最適化 4800 55,590.374 ミリ秒 414.6051 ミリ秒 387.8218 ミリ秒 0.46 部分並列最適化 4800 15,614.506 ミリ秒 115.6996 ミリ秒 102.5647 ミリ秒 0.13 結果から、 ループの並列化により、以前の実装 ( ) と比較して計算時間が 2 ～ 4 倍短縮されたことがわかります。これは非常に印象的ですが、純粋な計算の並列化は利用可能なすべての計算リソースを消費し、システム内の他のアプリケーションのリソース不足につながる可能性があることを覚えておくことが重要です。 for of i SpartialOptimisation 並列化は注意して使用する必要があります。 ご想像のとおり、これで終わりではありません 🙂 プラットフォームを知る – ベクトル化 とは、単一の要素で動作するコードを、複数の要素を同時に動作するコードに変換することです。 ベクトル化 これは複雑な問題のように聞こえるかもしれませんが、簡単な例でどのように機能するかを見てみましょう。 var a = new [] { 5, 7, 8, 1 }; var b = new [] { 4, 2, 2, 6 }; var c = new [] { 0, 0, 0, 0 }; for (var i = 0; i < 4; ++i) c[i] = a[i] + b[i]; 方法では、CPU レベルでの上記の ループの反復 の実行は次のように表すことができます。 非常に単純化した for 0 ここで何が起こっているか見てみましょう。CPU はメモリから配列 と の値を CPU レジスタにロードします (1 つのレジスタに 保持できません)。次に、CPU はこれらのレジストリに対してスカラー加算演算を実行し、その結果をメイン メモリ (配列 に直接) に書き込みます。このプロセスはループが終了する前に 4 回繰り返されます。 a b は 1 つの値しか c ベクトル化とは、特別な CPU レジスタ（ベクトル レジスタまたは SIMD (単一命令複数データ) レジスタ）と、対応する CPU 命令を使用して、複数の配列値に対して一度に操作を実行することを意味します。 ここで何が起こっているか見てみましょう。CPU はメモリから配列 と の値を CPU レジスタにロードします (ただし、今回は 1 つのベクトル レジスタが 配列値を保持できます)。次に、CPU はこれらのレジストリに対してベクトル加算演算を実行し、結果をメイン メモリ (配列 に直接) に書き込みます。一度に 2 つの値を操作するため、このプロセスは 4 回ではなく 2 回繰り返されます。 a b 2 つの c .NET でベクトル化を実装するには、 と 型を使用できます ( 名前空間の型も使用できますが、現在の実装では少し高度なため、ここでは使用しませんが、自由に試してみてください)。 Vector 型 Vector<T> System.Runtime.Intrinsics public void SpartialVectorOptimisations(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { for (var i = 0; i < sz; ++i) { if (matrix[i * sz + k] == NO_EDGE) { continue; } var ik_vec = new Vector<int>(matrix[i * sz + k]); var j = 0; for (; j < sz - Vector<int>.Count; j += Vector<int>.Count) { var ij_vec = new Vector<int>(matrix, i * sz + j); var ikj_vec = new Vector<int>(matrix, k * sz + j) + ik_vec; var lt_vec = Vector.LessThan(ij_vec, ikj_vec); if (lt_vec == new Vector<int>(-1)) { continue; } var r_vec = Vector.ConditionalSelect(lt_vec, ij_vec, ikj_vec); r_vec.CopyTo(matrix, i * sz + j); } for (; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } } } } ベクトル化されたコードは少し奇妙に見えるかもしれないので、ステップごとに説明していきましょう。 パスのベクトルを作成します: 。結果として、ベクトルが 型の値を 4 つ保持でき、 パスの重みが 5 に等しい場合、4 つの 5 のベクトル ([5, 5, 5, 5]) が作成されます。次に、各反復で、頂点 から頂点 へのパスを同時に計算します。 i ⭢ k var ik_vec = new Vector<int>(matrix[i * sz + k]) int i ⭢ k i j, j + 1, ..., j + Vector<int>.Count プロパティ ベクター レジスタに収まる 型の要素の数を返します。ベクター レジスタのサイズは CPU モデルによって異なるため、このプロパティは CPU によって異なる値を返すことがあります。 Vector<int>.Count int - 著者注 計算プロセス全体は、次の 3 つのステップに分けられます。 重み行列からの情報をベクトル と に読み込みます。 ij_vec ikj_vec ベクトルと ベクトルを比較し、最小値を新しいベクトル に選択します。 ij_vec ikj_vec r_vec 重み行列に書き戻します。 r_vec と 非常に簡単ですが、 説明が必要です。前述のように、ベクトルでは複数の値を同時に操作します。したがって、一部の操作の結果は 1 つにはなりません。つまり、比較操作 複数の値を比較するため、 または 返すことができません。そのため、代わりに、ベクトル と の対応する値の比較結果を含む新しいベクトルを返します ( の値が の値より小さい場合は 、それ以外の場合は )。返されたベクトル (それ自体) はあまり役に立ちませんが、 ベクトル操作を使用して、 ベクトルと ベクトルから必要な値を抽出し、新しいベクトル にすることができます。この操作は、入力ベクトルの対応する 2 つの値のうち小さい方に等しい値を持つ新しいベクトルを返します。つまり、ベクトル の値が に等しい場合、操作は から値を選択し、それ以外の場合は から値を選択します。 1 3 は 2 は Vector.LessThan(ij_vec, ikj_vec) true false ij_vec ikj_vec ij_vec ikj_vec -1 0 Vector.ConditionalSelect(lt_vec, ij_vec, ikj_vec) ij_vec ikj_vec r_vec lt_vec -1 ij_vec ikj_vec の他に、説明が必要な部分がもう 1 つあります。2 番目の非ベクトル化ループです。これは、重み行列のサイズが 値の積ではない場合に必要です。その場合、メイン ループはすべての値を処理できず (ベクトルの一部をロードできないため)、2 番目の非ベクトル化ループが残りの反復を計算します。 2 番 Vector<int>.Count この実装とこれまでのすべての実装の結果は次のとおりです。 方法 サイズ 平均 エラー 標準偏差 比率 ベースライン 300 27.525 ミリ秒 0.1937 ミリ秒 0.1617ミリ秒 1.00 部分最適化 300 12.399ミリ秒 0.0943 ミリ秒 0.0882 ミリ秒 0.45 部分並列最適化 300 6.252 ミリ秒 0.0211ミリ秒 0.0187ミリ秒 0.23 部分ベクトル最適化 300 3.056ミリ秒 0.0301ミリ秒 0.0281 ミリ秒 0.11 ベースライン 600 217.897 ミリ秒 1.6415ミリ秒 1.5355 ミリ秒 1.00 部分最適化 600 99.122 ミリ秒 0.8230 ミリ秒 0.7698 ミリ秒 0.45 部分並列最適化 600 35.825 ミリ秒 0.1003 ミリ秒 0.0837 ミリ秒 0.16 部分ベクトル最適化 600 24.378 ミリ秒 0.4320 ミリ秒 0.4041 ミリ秒 0.11 ベースライン 1200 1,763.335 ミリ秒 7.4561 ミリ秒 6.2262 ミリ秒 1.00 部分最適化 1200 766.675 ミリ秒 6.1147 ミリ秒 5.7197 ミリ秒 0.43 部分並列最適化 1200 248.801 ミリ秒 0.6040 ミリ秒 0.5043 ミリ秒 0.14 部分ベクトル最適化 1200 185.628 ミリ秒 2.1240 ミリ秒 1.9868 ミリ秒 0.11 ベースライン 2400 14,533.335 ミリ秒 63.3518 ミリ秒 52.9016 ミリ秒 1.00 部分最適化 2400 6,507.878 ミリ秒 28.2317 ミリ秒 26.4079 ミリ秒 0.45 部分並列最適化 2400 2,076.403 ミリ秒 20.8320 ミリ秒 19.4863 ミリ秒 0.14 部分ベクトル最適化 2400 2,568.676 ミリ秒 31.7359 ミリ秒 29.6858 ミリ秒 0.18 ベースライン 4800 119,768.219 ミリ秒 181.5514 ミリ秒 160.9406 ミリ秒 1.00 部分最適化 4800 55,590.374 ミリ秒 414.6051 ミリ秒 387.8218 ミリ秒 0.46 部分並列最適化 4800 15,614.506 ミリ秒 115.6996 ミリ秒 102.5647 ミリ秒 0.13 部分ベクトル最適化 4800 18,257.991 ミリ秒 84.5978 ミリ秒 79.1329 ミリ秒 0.15 結果から明らかなように、ベクトル化により計算時間が大幅に短縮されました。非並列化バージョン ( ) と比較して 3 倍から 4 倍短縮されました。ここで興味深いのは、ベクトル化バージョンは、より小さなグラフ (2400 頂点未満) でも並列バージョン ( ) よりも優れていることです。 SpartialOptimisation SpartialParallelOptimisations そして最後に、ベクトル化と並列処理を組み合わせてみましょう。 ベクトル化の実際の応用に興味がある場合は、 ベクトル化した の一連の投稿を読むことをお勧めします (これらの結果は、後に のガベージ コレクターの更新に反映されました)。 Array.Sort Dan Shechter .NET 5 プラットフォームとハードウェアを理解しましょう – ベクトル化と並列処理! 最後の実装は並列化とベクトル化の取り組みを組み合わせたもので、正直に言うと個性に欠けています 🙂 なぜなら、 の 本体を のベクトル化されたループに置き換えただけだからです。 SpartialParallelOptimisations Parallel.For SpartialVectorOptimisations public void SpartialParallelVectorOptimisations(int[] matrix, int sz) { for (var k = 0; k < sz; ++k) { Parallel.For(0, sz, i => { if (matrix[i * sz + k] == NO_EDGE) { return; } var ik_vec = new Vector<int>(matrix[i * sz + k]); var j = 0; for (; j < sz - Vector<int>.Count; j += Vector<int>.Count) { var ij_vec = new Vector<int>(matrix, i * sz + j); var ikj_vec = new Vector<int>(matrix, k * sz + j) + ik_vec; var lt_vec = Vector.LessThan(ij_vec, ikj_vec); if (lt_vec == new Vector<int>(-1)) { continue; } var r_vec = Vector.ConditionalSelect(lt_vec, ij_vec, ikj_vec); r_vec.CopyTo(matrix, i * sz + j); } for (; j < sz; ++j) { var distance = matrix[i * sz + k] + matrix[k * sz + j]; if (matrix[i * sz + j] > distance) { matrix[i * sz + j] = distance; } } }); } } この投稿のすべての結果は以下に記載されています。予想どおり、並列処理とベクトル化を同時に使用すると最良の結果が得られ、最初の実装と比較して計算時間が最大 25 倍 (1200 頂点のグラフの場合) 短縮されました。 方法 サイズ 平均 エラー 標準偏差 比率 ベースライン 300 27.525 ミリ秒 0.1937ミリ秒 0.1617ミリ秒 1.00 部分最適化 300 12.399ミリ秒 0.0943 ミリ秒 0.0882 ミリ秒 0.45 部分並列最適化 300 6.252 ミリ秒 0.0211ミリ秒 0.0187ミリ秒 0.23 部分ベクトル最適化 300 3.056ミリ秒 0.0301ミリ秒 0.0281 ミリ秒 0.11 部分並列ベクトル最適化 300 3.008 ミリ秒 0.0075ミリ秒 0.0066ミリ秒 0.11 ベースライン 600 217.897 ミリ秒 1.6415 ミリ秒 1.5355 ミリ秒 1.00 部分最適化 600 99.122 ミリ秒 0.8230 ミリ秒 0.7698 ミリ秒 0.45 部分並列最適化 600 35.825 ミリ秒 0.1003 ミリ秒 0.0837 ミリ秒 0.16 部分ベクトル最適化 600 24.378 ミリ秒 0.4320 ミリ秒 0.4041 ミリ秒 0.11 部分並列ベクトル最適化 600 13.425ミリ秒 0.0319ミリ秒 0.0283 ミリ秒 0.06 ベースライン 1200 1,763.335 ミリ秒 7.4561 ミリ秒 6.2262 ミリ秒 1.00 部分最適化 1200 766.675 ミリ秒 6.1147 ミリ秒 5.7197 ミリ秒 0.43 部分並列最適化 1200 248.801 ミリ秒 0.6040 ミリ秒 0.5043 ミリ秒 0.14 部分ベクトル最適化 1200 185.628 ミリ秒 2.1240 ミリ秒 1.9868 ミリ秒 0.11 部分並列ベクトル最適化 1200 76.770 ミリ秒 0.3021 ミリ秒 0.2522 ミリ秒 0.04 ベースライン 2400 14,533.335 ミリ秒 63.3518 ミリ秒 52.9016 ミリ秒 1.00 部分最適化 2400 6,507.878 ミリ秒 28.2317 ミリ秒 26.4079 ミリ秒 0.45 部分並列最適化 2400 2,076.403 ミリ秒 20.8320 ミリ秒 19.4863 ミリ秒 0.14 部分ベクトル最適化 2400 2,568.676 ミリ秒 31.7359 ミリ秒 29.6858 ミリ秒 0.18 部分並列ベクトル最適化 2400 1,281.877 ミリ秒 25.1127 ミリ秒 64.8239 ミリ秒 0.09 ベースライン 4800 119,768.219 ミリ秒 181.5514 ミリ秒 160.9406 ミリ秒 1.00 部分最適化 4800 55,590.374 ミリ秒 414.6051 ミリ秒 387.8218 ミリ秒 0.46 部分並列最適化 4800 15,614.506 ミリ秒 115.6996 ミリ秒 102.5647 ミリ秒 0.13 部分ベクトル最適化 4800 18,257.991 ミリ秒 84.5978 ミリ秒 79.1329 ミリ秒 0.15 部分並列ベクトル最適化 4800 12,785.824 ミリ秒 98.6947 ミリ秒 87.4903 ミリ秒 0.11 結論 この投稿では、全ペア最短経路問題に少し踏み込んで、それを解決するために C# で Floyd-Warshall アルゴリズムを実装しました。また、データを尊重し、.NET とハードウェアの低レベル機能を利用するように実装を更新しました。 この投稿では、私たちは「オールイン」をしました。しかし、実際のアプリケーションでは、私たちだけではないということを覚えておくことが重要です。ハードコアな並列化は、システム パフォーマンスを大幅に低下させ、メリットよりもデメリットをもたらす可能性があります。一方、ベクトル化は、プラットフォームに依存しない方法で実行すれば、少し安全です。一部のベクトル命令は、特定の CPU でのみ使用できる場合があることに注意してください。 この記事を楽しんでいただき、たった 3 サイクルでアルゴリズムからパフォーマンスを少しだけ「絞り出す」楽しさを感じていただければ幸いです 🙂 参考文献 Floyd, RW アルゴリズム 97: 最短経路 / RW Floyd // Communications of the ACM. – 1962. – Vol. 5、№. 6. – P. 345-. Ingerman、PZ アルゴリズム 141: パス マトリックス / PZ Ingerman // Communications of the ACM. – 1962. – Vol. 5、№. 11. – P. 556。
