フォード・ファルカーソンのアルゴリムで最大流-最小切断問題を解く
PRMLの8.3節では、マルコフ確率場を応用した事例として画像のノイズ除去がでてきますが、PRMLで紹介されている「グラフカットアルゴリズム」を試すには最大流問題を解く必要があるのでアルゴリズムの勉強もかねて実装してみました。
最大フロー問題または最大流問題とは、各枝に容量(capacity)と流量(flow)が設定された有効グラフにおいて、ある始点（ソース）から終点（シンク）へのフローで最大となるフローを求める問題です。
この問題を解くアルゴリズムとしてフォード・ファルカーソンのアルゴリズムが一般的に良く知られていて、僕の持っているアルゴリズムC・新版—基礎・データ構造・整列・探索にも載っています。
アルゴリズムCから引用すると、フォード・ファルカーソンのアルゴリズムは「辺の流量がすべてゼロである状態から開始し、ソースからシンクに向かう道で、飽和した前向きの辺や流量ゼロの後ろ向きの辺を含まないものに沿って流量を増やし、ネットワークにそのような道がなくなるまで繰り返す」というように要約できます。「後ろ向きの辺」が何かというと、元の向きと逆向きの辺を考えて、その辺を流れるときは逆に流量を減らすというものです。この逆向きの辺を考えないと、最初に変な道の選び方をしてしまうと、もう最大フローにはたどり着かなくなってしまいます。
また、フォード・ファルカーソンのアルゴリズム自体は「ソースからシンクに向かう道で、飽和した前向きの辺や流量ゼロの後ろ向きの辺を含まないもの」を探すアルゴリズム自体を規定しているわけではないので、このアルゴリズムは自分で探索アルゴリズムを組む必要があります。アルゴリズムCでは最良優先探索を使用しています。
以下では、アルゴリズムCにならって最良優先検索で最大流および最小切断を求めてみます。
対象となる最大流問題は、アルゴリズムCの図33.2のネットワークを対象としています。
（図がダウンロードできなかったので、貼付けられず。。ご容赦。。）
#include <iostream> #include <vector> #include <math.h> using namespace std; namespace GraphCut { //NodeとEdgeで相互参照するため、あらかじめ宣言しておく。 class Node; class Edge { public: Node* s_Node_; Node* e_Node_; Edge (Node* s_Node, Node* e_Node, double size); double AddFlow (double flow); void set_reverse_edge (Edge* re); double get_flow (); void set_flow (double flow); double get_size (); private: Edge* reverse_edge_; double flow_; double size_; }; class Node { public: vector<Edge*> vecFlowPath; Node* dad_; string name_; Node (string name); void AddFlowPath (Node* e_node, double size); double get_val (); void set_val (double val); void visit (); bool get_is_visited (); Edge* dad_path_; void reset (); private: double val_; bool is_visited_; }; Edge::Edge (Node* s_Node, Node* e_Node, double size) : s_Node_(s_Node), e_Node_(e_Node), size_(size), flow_(0) {} double Edge::AddFlow (double flow) { flow_ += flow; reverse_edge_->set_flow(-flow); return flow; } double Edge::get_flow () { return flow_; } void Edge::set_flow (double flow) { flow_ = flow; } void Edge::set_reverse_edge (Edge* re) { reverse_edge_ = re; } double Edge::get_size () { return size_; } Node::Node (string name) : name_(name), val_(0), is_visited_(false) {} void Node::AddFlowPath (Node* e_node, double size) { Edge* e = new Edge(this, e_node, size); vecFlowPath.push_back(e); Edge* re = new Edge (e_node, this, -size); re->set_reverse_edge(e); e->set_reverse_edge(re); e_node->vecFlowPath.push_back(re); } double Node::get_val () { return val_; } void Node::set_val (double val) { val_ = val; } void Node::visit () { is_visited_ = true; } bool Node::get_is_visited () { return is_visited_; } void Node::reset () { is_visited_ = false; val_ = 0; } } using namespace GraphCut; //ノードnを含む最小切断を再帰的に求めてmincutに追加する void SearchMinCut (Node* n, vector<Node*>& mincut) { mincut.push_back(n); n->visit(); vector<Edge*> vecAdjEdges = n->vecFlowPath; for (int j=0; j<vecAdjEdges.size(); j++) { Edge* adjEdge = vecAdjEdges[j]; Node* adjNode = adjEdge->e_Node_; if (adjNode->get_is_visited()) continue; //飽和した前向きの辺や、流量ゼロの後ろ向きの辺ではない場合 if ((adjEdge->get_size() > 0 && adjEdge->get_flow() != adjEdge->get_size()) || (adjEdge->get_size() < 0 && adjEdge->get_flow() != 0)) { SearchMinCut(adjNode, mincut); } } } //全ノードvecNodesを順位優先で検索して、流量の増加の最も多い経路を見つける。 bool PFS (vector<Node*> vecNodes, Node* source, Node* sink) { printf("start priority-first search\n"); //初期化 for (int j=0; j<vecNodes.size(); j++) { Node* n = vecNodes[j]; n->reset(); } source->set_val(INT_MAX); Node* v = source; bool findSource2Sink = false; bool findMax = true; while (findMax) { findMax = false; v->visit(); double val = v->get_val(); //新たに訪問済になったノードvの隣接ノードのpriority（流量）の更新 vector<Edge*> vecAdjEdges = v->vecFlowPath; for (int j=0; j<vecAdjEdges.size(); j++) { Edge* adjEdge = vecAdjEdges[j]; Node* adjNode = adjEdge->e_Node_; if (adjNode->get_is_visited()) continue; //priorityの更新 //vを経由した場合の、adjNodeへの流量を求める double priority = -(adjEdge->get_flow()); double size = adjEdge->get_size(); if (size > 0) priority += size; if (priority > val) priority = val; //新たに求めたadjNodeへの流量が、既に設定された流量を上回る場合、adjNodeへの経路と流量を更新する if (adjNode->get_val() < priority) { adjNode->set_val(priority); adjNode->dad_ = v; adjNode->dad_path_ = adjEdge; } } //流量最大のノードを求める double maxval = 0; int max = -1; int nodeNum = vecNodes.size(); for (int j=0; j<nodeNum; j++) { Node* n = vecNodes[j]; if (n->get_is_visited()) continue; if (n->get_val() > maxval) { maxval = n->get_val(); max = j; } } //流量最大のノードが見つかった場合 if (max >= 0) { findMax = true; v = vecNodes[max]; if (v == sink) { //ソースからシンクまでの経路が見つかった場合 findSource2Sink = true; break; } } } return findSource2Sink; } int main (int argc, char * const argv[]) { //ノードと辺の初期化 vector<Node*> vecNodes; Node* a = new Node("a"); Node* b = new Node("b"); Node* c = new Node("c"); Node* d = new Node("d"); Node* e = new Node("e"); Node* f = new Node("f"); vecNodes.push_back(a); vecNodes.push_back(b); vecNodes.push_back(c); vecNodes.push_back(d); vecNodes.push_back(e); vecNodes.push_back(f); a->AddFlowPath(b,6); a->AddFlowPath(c,8); b->AddFlowPath(d,6); b->AddFlowPath(e,3); c->AddFlowPath(d,3); c->AddFlowPath(e,3); d->AddFlowPath(f,8); e->AddFlowPath(f,6); Node* source = a; Node* sink = f; //ソースからシンクまでの経路が尽きるまで、流量最大経路の探索を繰り返す while (PFS(vecNodes, source, sink)) { //ソースからシンクまでの流量が増加する経路が発見された場合、シンクから逆に発見された経路をたどり、流量を増加させる。 Node* n = sink; while (n != source) { Edge* e = n->dad_path_; e->AddFlow(sink->get_val()); printf("add flow to %s-%s by %f\n", e->s_Node_->name_.c_str(), e->e_Node_->name_.c_str(), sink->get_val()); n = n->dad_; } } //最小切断を求めるために全ノードを初期化 for (int j=0; j<vecNodes.size(); j++) { Node* n = vecNodes[j]; n->reset(); } //再帰的に最小切断を求める vector<Node*> mincut; SearchMinCut(source, mincut); cout << "display mincut" << endl; for (int j=0; j<mincut.size(); j++) { Node* n = mincut[j]; printf("%s\n", n->name_.c_str()); } return 0; }
実行結果
start priority-first search add flow to d-f by 6.000000 add flow to b-d by 6.000000 add flow to a-b by 6.000000 start priority-first search add flow to e-f by 3.000000 add flow to c-e by 3.000000 add flow to a-c by 3.000000 start priority-first search add flow to e-f by 3.000000 add flow to b-e by 3.000000 add flow to d-b by 3.000000 add flow to c-d by 3.000000 add flow to a-c by 3.000000 start priority-first search display mincut a c
というわけで、最小切断がa,cだということが分かりました。
PFS内で、流量最大のノードを探索している箇所をヒープソートにしたりすれば多分もっと早くなります。
ちなみにですが、石川先生のグラフカットのサーベイ論文によれば、最大流問題を解くアルゴリズムとして、フォード・ファルカーソンのアルゴリズム以外にpush-relabelアルゴリズム(PDF)が知られていて、フォード・ファルカーソンが計算量O(n^3)に対し、push-relabelはO(nmlog(n^2/m))（nは頂点の数、mは辺の数）で、push-relabelの方が速いそうです。
ですが、次回はとりあえずこのままフォード・ファルカーソンのアルゴリズムで、グラフカットによるエネルギー最小化問題を解いてみることにします。
