Warshall Floyd's algorithm ワーシャル フロイド

All Pairs Shortest Path (APSP) problem (全対最短経路問題) は、グラフ G(V, E) に対して、G に含まれるノードの全てのペアの最短経路（距離）を求める問題です。G に負の重みのあるエッジがなければ、この問題は各ノードを起点として Dijkstra's Algorithm を |V| 回行うことによって解くことができます。この方法のオーダーは O(|V|³) または O(|V|(|E| + |V|)log|V|) となります。

Floyd's algorithm は APSP 問題を O(|V|³) で解くことができます。しかも Floyd's algorithm は、G に負の閉路がない限り、G に負の重みを持つエッジが存在しても正しく動作します。負の閉路とは、その閉路を成すエッジの重みの合計が負となっている閉路です（負の閉路をもつとき G に最短経路は存在しません）。さらに、Floyd's algorithm は、G に負の閉路が存在するか否かを判定することができます。アルゴリズムが終了した時点で、G のあるノード v からノード v (それ自身) への最短コストが負になっていれば、G に負の閉路が存在します。

Floyd's algorithm は、ノードの全てのペア [i, j] に対する最短コストを与える配列 A[i, j] を以下の方法で求めます。まず、G(V, E) のノードを {1, 2, 3, ... |V|} とします。

A^k[i, j] をノード V^k = {1, 2, 3, ... k} のみを経由するノード i とノード j の最短コスト、その時の経路の１つを P^k[i, j] とします。Floyd's algorithm は A^k-1 を利用して A^k を求めるということを k = {1, 2, 3, ... |V|} に対して行い、A^|V| つまり A[i, j] (および P[i, j])を求めます（Dynamic Programming (動的計画法)）。Floyd's algorithm は以下のように数学的帰納法によって証明・説明することができます。

1.	A0 ( k = 0 ) の時の証明 A0[i, j] は、ノード i とノード j 以外に経由ノードが存在しないので、単にノード i とノード j を繋ぐエッジの重みになります。つまり、入力されたグラフ G に対応させます。G のノード i からノード j に向かって重み d のエッジがあれば A[i, j] = d とし、エッジがなければ A[i, j] = ∞ とます。これを全対の [i, j] に対して施します。ただし A[i, i] = 0 とします。
	この時点で、A0[i, j] がノード i からノード j への最短コストであることは明白なので、帰納法の基礎が確立できました。
2.	Ak ( k > 0 ) の時の証明 k = {1, 2, 3, ... |V|} について、Ak-1 に基づいて Ak が正しく求められるか？を示します。
	Pk[i, j] がノード k を経由しない場合とする場合を考えます。 (i) Pk[i, j] がノード k を経由しない場合、Pk[i, j] は端点であるノード i とノード j 以外では Vk-1 = {1, 2, 3, ... k-1} に属するノードしか通らないので、Pk-1[i, j] と等しいと考えることができます。よって、この場合は Ak[i, j] = Ak-1[i, j] となります。 (ii) Pk[i, j] がノード k を経由する場合、Pk[i, j] はノード k によって　i → k と k → j の２つの部分路に分けられます。この２つの部分路はどちらも Vk-1 = {1, 2, 3, ... k-1} に属するノードしか経由しません（下図参照）。従ってノード k を経由する最短経路の部分路は Pk-1[i, k] と Pk-1[k, j] となります。つまり、Ak[i, j] = Ak-1[i, k] + Ak-1[k, j] となります。
	(i)(ii) より、i, j の全ての対に対して、      Ak[i, j] = min( Ak-1[i, j], Ak-1[i, k] + Ak-1[k, j]) を行い Ak を求めます。

このアルゴリズムを実装する時に、A^k[i, j] を求めるために、３次元のメモリ空間を必要と思うかもしれません。つまり A^k[i, j] の計算中に A^k-1[i, j] の値が変化してしまうのでは？と考えてしまいます。しかし、これは以下の理由により問題はありません。

A^k[k, k] = 0 より、

     A^k[i, k] = min( A^k-1[i, k], A^k-1[i, k] + A^k-1[k, k] ) = A^k-1[i, k]
     A^k[k, j] = min( A^k-1[k, j], A^k-1[k, j] + A^k-1[k, k] ) = A^k-1[k, j]

よって、ノードの全てのペア [i, j] について A^k-1[i, j] を A^k[i, j] に上書きしても不都合はありません。

int N;
int A[MAX][MAX];

void floyd(){
    for ( int k = 0; k < N; k++ ){
        for ( int i = 0; i < N; i++ ){
            for ( int j = 0; j < N; j++ ){
                A[i][j] = min( A[i][j], A[i][k] + A[k][j] );
            }
        }
    }
}

Priority First Search 順位優先探索

Dijkstra's Algorithm のオーダーは、O(|V|^2) (ノードの数の２乗)です。これは、u を S に追加するためのループ（|V| 回）の中で、u を選択するためのループ（|V| 回）を行っていることから、簡単に導くことができます。このオーダーは、隣接行列で実装しても、隣接リストで実装しても同じですが、疎なグラフに対しては隣接リストで実装すべきです。（以下のテーブルを参照して下さい）。

タイプ	演算回数	オーダー
隣接行列	|V|^2 + |V|^2	|V|^2
隣接リスト	|V|^2 + |E|	|V|^2

Dijkstra's Algorithm は、隣接リストによる実装と、Priority Queue (Heap) を応用することによって、効率化することができます。このアルゴリズムは、Priority First Search (順位優先探索) とも呼ばれ、以下のように実装します。

	グラフ G(V, E) のノード全体の集合を V、 source をスタート地点、 S を「 source から S 内の各ノード v への最短経路が必ず S 内にあるような、最短経路の部分的な解の集合」とします（Dijkstra's Algorithm 参照）。
	各計算ステップにおいて、V - S の各ノード v について、S 内のノードのみを経由した source から v への最短経路のコスト(距離)を d[v] とします。また、各ノード間のコスト（距離） (x, y) は、行列 M[x][y] に記録されているものとします。
1.	初期状態で、S は空です。 source に対して    d[s] = 0 source 以外の V に属する全てのノード v に対して    d[v] = ∞ と初期化します。
	d[v] をキーとして、V - S に属するノードを minimum heap として構築します。これを heap とします。
2.	heap から d[u] が最小であるノード u を取り出します。   u = heap.extractMin()
	u を S に追加すると同時に、u に隣接しかつ V - S に属する全てのノード v に対する値を以下のように更新します。
	d[v] = min( d[v], d[u] + M[u][v] )     d[v] が更新された場合は、heap を更新します。         heap.update( v )　(v から upheap を施します)
	この処理を S = V となるまで繰り返します。

Priority First Search のオーダーは、O((|V| + |E|)log|V|) となります。u を heap から取り出すために |V|log|V|、d[v] を更新するために |E|log|V| の演算回数が必要です。

Priority First Search のプログラム例を以下に示します。

class Graph{
    public:
    vector<vector<int> > adjList;
    vector<vector<int>::iterator> pos;
    
    Graph(){}
    Graph( int size ){ resize( size );}
    
    void resize( int size ){
        adjList.resize( size ); pos.resize( size );
        for ( int i = 0; i < size; i++ ) adjList[i].clear();
    }
    
    int size(){ return adjList.size(); }
    
    void insert( int i, int j ){
        adjList[i].push_back( j );
    }
    
    int next( int i ){
        if ( pos[i] == adjList[i].end() ) return -1;
        return *pos[i]++;
    }
    
    void reset( int i ){ pos[i] = adjList[i].begin(); }
    void reset(){ for( int i = 0; i < size(); i++ ) reset(i); }
};

class Heap{
    public:
    int size;
    int H[MAX + 1]; // heap buffer
    int index[MAX + 1]; // index of v in the heap buffer
    int *key; // keys for heap condition
    
    Heap(){}
    Heap( int size, int d[] ): size(size){
        key = d;
    }
    
    void construct(){
        for ( int i = 0; i < size; i++ ){
            H[i + 1] = i; // heap is 1 base not 0
            index[i] = i + 1;
        }
        for ( int i = size / 2; i >= 1; i-- ) downheap( i );
    }
    
    int extractMin(){
        int v = H[1];
        H[1] = H[size];
        index[ H[size] ] = 1;
        size--;
        downheap(1);
        return v;
    }
    
    int update( int v ){
        upheap( index[v] );
    }
    
    private:
    void downheap( int k ){
        int j, v;
        v = H[k];
        while( k <= size/2 ){
            j = k + k;
            if ( j < size && key[ H[j] ] > key[ H[j+1] ] ) j++;
            if ( key[v] <= key[ H[j] ] ) break;
            H[k] = H[j];
            index[ H[j] ] = k;
            k = j;
        }
        H[k] = v;
        index[ v ] = k;
    }
    
    void upheap( int k){
        int v;
        v = H[k];
        while ( k > 1 && key[ H[k/2] ] >= key[v] ){
            H[k] = H[k/2];
            index[ H[k/2] ] = k;
            k = k/2;
        }
        H[k] = v;
        index[v] = k;
    }
};

int size;
int M[MAX][MAX];
Graph graph;

void dijkstra( int s, int d[] ){
    bool visited[MAX]; // S に属するノードは true
    for ( int i = 0; i < size; i++ ){
        d[i] = INT_MAX;
        visited[i] = false;
    }
    
    d[s] = 0; // 最初に s が u として選ばれる
    
    Heap heap = Heap( size, d );
    heap.construct();
    
    graph.reset();
    
    while( heap.size >= 1 ){
        int u = heap.extractMin();
        
        visited[u] = true; // u を S に追加
        
        int v;
        while ( ( v = graph.next( u ) ) != -1 ){
            if ( visited[v] ) continue;
            if ( d[u] + M[u][v] < d[v] ){
                d[v] =  d[u] + M[u][v];
                heap.update( v );
            }
        }
    }
}

このプログラム例ですと、自作の Heap を使っているので、プログラムが分かりづらくなってしまいました。私は、普段は自作の Heap を使うことはなく、以下のように既存の Priority queue に候補となるノードを突っ込んでいくアルゴリズムを実装しています。

class Graph{
    public:
    vector<vector<int> > adjList;
    vector<vector<int>::iterator> pos;
    
    Graph(){}
    Graph( int size ){ resize( size );}
    
    void resize( int size ){
        adjList.resize( size ); pos.resize( size );
        for ( int i = 0; i < size; i++ ) adjList[i].clear();
    }
    
    int size(){ return adjList.size(); }
    
    void insert( int i, int j ){
        adjList[i].push_back( j );
    }
    
    int next( int i ){
        if ( pos[i] == adjList[i].end() ) return -1;
        return *pos[i]++;
    }
    
    void reset( int i ){ pos[i] = adjList[i].begin(); }
    void reset(){ for( int i = 0; i < size(); i++ ) reset(i); }
};

class Node{
    public:
    int id, cost;
    Node(){}
    Node( int id, int cost ): id(id), cost(cost){}
    
    bool operator < ( const Node &n ) const{
        return cost > n.cost;
    }
};

int size;
int M[MAX][MAX];
Graph graph;

void dijkstra( int s, int d[] ){
    bool visited[MAX]; // S に属するノードは true
    for ( int i = 0; i < size; i++ ){
        d[i] = INT_MAX;
        visited[i] = false;
    }
    
    d[s] = 0;
    
    priority_queue<Node> PQ;
    PQ.push( Node( s, 0 ) );// 最初に s が u として選ばれる
    
    graph.reset();
    
    while( !PQ.empty() ){
        int u = PQ.top().id; PQ.pop();
        
        visited[u] = true; // u を S に追加
        
        int v;
        while ( ( v = graph.next( u ) ) != -1 ){
            if ( visited[v] ) continue;
            if ( d[u] + M[u][v] < d[v] ){
                d[v] =  d[u] + M[u][v];
                PQ.push( Node( v, d[v] ) );
            }
        }
    }
}

純粋に Heap を用いるのと比べて効率はやや劣りますが（なぜなら priority queue の中に余分なゴミが溜まっていくから）、オーダーがたかだか O((|V| + |E|)log|E|) になるだけなので、私の経験上、実用には問題ないようです。

Dijkstra's Algorithm ダイクストラのアルゴリズム

グラフ G(V, E) における Single Source Shortest Path を求めるためのアルゴリズムの１つが、 Dijkstra によって発見された Dijkstra's Algorithm (ダイクストラのアルゴリズム) です。

Dijkstra's Algorithm は Greedy Algorithm で、以下のように実装します。

	グラフ G(V, E) のノード全体の集合を V、 s をスタート地点 (source)、 S を「 s から S 内の各ノード v への最短経路が必ず S 内にあるような、最短経路の部分的な解の集合」とします（下図を参照して下さい）。
	各計算ステップにおいて、V - S の各ノード v について、S 内のノードのみを経由した s から v への最短経路のコスト(距離)を d[v] とします。また、各ノード間のコスト（距離） (x, y) は、行列 M[x][y] に記録されているものとします。
1.	初期状態で、S は空です。 s に対して     d[s] = 0 s 以外の V に属する全てのノード v に対して     d[v] = ∞ と初期化します。
2.	ノードの集合 V - S の中から、d[u] が最小であるノード u を選択します。
	u を S に追加すると同時に、u に隣接しかつ V - S に属する全てのノード v に対する値を以下のように更新します。
	d[v] = min( d[v], d[u] + M[u][v] )
	この処理を S = V となるまで繰り返します。

以下の図は、u が選ばれた直後、v を更新する様子を示しています。

2. の計算ステップ直後（つまり次に u を選ぶ直前）において、 d[v] には、s から S 内のノードのみを経由した v までの最短コストが記録されています。

すなわち、処理が終了した時点で、V に属する全てのノード d[v] には、s から v までの最短コスト（距離）が記録されています。

コスト（距離）だけでなく、経路も求めたい場合は、Shortest spanning tree における v の親である pi[v] を用意し、pi[s] = -1と初期化し、更新部分を以下のように書き換えます。

もし、d[u] + M[u][v] ＜ d[v] ならば、     d[v] = d[u] + M[u][v]     pi[v] = u;

ダイクストラのアルゴリズムのプログラム例です。

int N;
int M[MAX][MAX];

void dijkstra( int s, int d[] ){
    bool visited[MAX]; // S に属するノードは true
    for ( int i = 0; i < N; i++ ){
        d[i] = INT_MAX;
        visited[i] = false;
    }
    
    d[s] = 0; // 最初に s が u として選ばれる
    
    while( 1 ){
        int u; // 最適なノード u を選ぶ
        int mincost = INT_MAX;
        for ( int i = 0; i < N; i++ ){
            if ( !visited[i] && d[i] < mincost ){
                mincost = d[i]; u = i;
            }
        }
        
        // u が存在しなかったとき、つまり S がこれ以上増えないとき、終了
        if ( mincost == INT_MAX ) break;
        
        visited[u] = true; // u を S に追加
        
        for ( int v = 0; v < N; v++ ){
            // v が S に含まれる場合　または　エッジ(u, v)がない場合は　無視
            if ( visited[v] || M[u][v] == INT_MAX ) continue;
            
            d[v] = min( d[v], d[u] + M[u][v] );
        }
    }
}

経路を生成するために pi を加えたプログラムです。

001 int N;
002 int M[MAX][MAX];
003 
004 void dijkstra( int source, int d[], int pi[] ){
005     bool visited[MAX]; // S に属するノードは true
006     for ( int i = 0; i < N; i++ ){
007         d[i] = INT_MAX;
008         visited[i] = false;
009     }
010     
011     d[source] = 0; // 最初に s が u として選ばれる
012     pi[source] = -1;
013     
014     while( 1 ){
015         int u; // 最適なノード u を選ぶ
016         int mincost = INT_MAX;
017         for ( int i = 0; i < N; i++ ){
018             if ( !visited[i] && d[i] < mincost ){
019                 mincost = d[i]; u = i;
020             }
021         }
022         
023         // u が存在しなかったとき、つまり S がこれ以上増えないとき、終了
024         if ( mincost == INT_MAX ) break;
025         
026         visited[u] = true; // u を S に追加
027         
028         for ( int v = 0; v < N; v++ ){
029             // v が S に含まれる場合　または　エッジ(u, v)がない場合は　無視
030             if ( visited[v] || M[u][v] == INT_MAX ) continue;
031             if ( d[u] + M[u][v] < d[v] ){
032                 d[v] = d[u] + M[u][v];
033                 pi[v] = u;
034             }
035         }
036     }
037 }

最短経路

Shortest path problem (最短経路問題) とは、重みつきグラフ G(V, E) において、ある与えられたノード x, y を接続する path の中で、エッジの重みの総和が最小である path を求める問題です。この問題には、様々な応用がありますが、以下の２つが重要です。

	Single Source Shortest Path (SSSP) problem
	グラフ G において、与えられた１つのノード s から、他の全てのノードへの最短経路を求める問題。
	All Pairs Shortest Path (APSP) problem
	グラフ G において、全ての "ノードのペア" 間の最短経路を求める問題。

Single Source Shortest Path Problem (単一始点最短経路)

G(V, E) をエッジのコスト（長さ）が非負である重みつきグラフとし、s を、G のノードとします。ただし、s から G の全てのノードに対して経路があるものとします。このとき、s を根とし、s から G の全てのノードへの最短経路を包含する G の spanning tree T が存在します。このような tree を shortest path spanning tree といいます。以下に例を示します。