最大流問題をプリフロー・プッシュ (Push-Relabel) 法 で解く

はじめに

画像処理でグラフカットを使う場面があったのですが，どうせなら最大流/最小カットのアルゴリズムについても軽くおさえておくかという事で，本記事を書いています．
バグ・間違い等ありましたらご指摘頂けるとうれしいです．

この記事で扱う事

• プリフロー・プッシュ法のざっくりとした説明
• Pythonによるプリフロー・プッシュ法の実装

この記事で扱わない事

• 難しい理論
• 性能（速度）の追求

最大流問題

有向グラフ $G=( \{ s \} \cup \{ t \} \cup V, E)$ を考える．各辺 $e \in E$ にはフローを流すことのできる容量 $c$ が決められており，$c$ 以上のフローを流すことはできない．ソース $s$ からシンク $t$ への流量を最大にするようなフローを求める．

今回の実装はフローの総流量のみを出力しています．

プリフロー・プッシュ法

Push-Relabel アルゴリズムとも呼ばれます．

用語定義

これらの定義は，「劣モジュラ最適化と機械学習」を参考にさせて頂きました．

• 残存量 (excess)
  全ての頂点は，入ってきたフローを貯めておくタンクのようなものを持っています．このタンクに入っているフローの量を 残存量と呼びます．
  つまり，(頂点 $e$ の残存量) = (頂点 $e$ に入ってくるフロー総量) - (頂点 $e$ から出ていくフロー総量)　が成立します．

• プリフロー
  残存量の存在を許容しているフローをプリフローと呼びます．プリフローは，最大流問題の制約を満たしていません．

• 活性頂点 (active vertex)
  残存量が正の頂点を活性頂点と呼びます．

• 距離ラベル（"高さ"とも呼ばれる）
  全ての頂点には"高さ"のような概念が存在します．これは，各頂点が「フローを流す(Push)」という操作をするための条件に関係し，反復的に更新されます．また，ソースの距離ラベルは $|V|+2$ に，シンクの距離ラベルは $0$ に設定されています．

• 可能辺 (admissible arc)
  次の条件を満たす辺を可能辺と呼びます．
  (辺 $e$ の始点の距離ラベル) = (辺 $e$ の終点の距離ラベル) + 1

アルゴリズム概要

とてもざっくりと，アルゴリズムの流れについて追っていきます．

1. 初期化
   ソースから出ている全ての辺に，容量いっぱいまでフローを流します．さらに，ソースとシンクを除く頂点の距離ラベルを $0$ に設定します．
2. 活性頂点 $i$ の選択
   活性頂点が無ければ，プリフローを最大流として出力してアルゴリズムを終了します．
3. 可能辺 $e$ の選択
   頂点 $i$ を始点とする可能辺を選択します．可能辺がなければ Relabel操作 を行った後，ステップ1に戻ります．
4. Push操作
   選択した可能辺 $e$ に対して Push操作 を行った後，ステップ1に戻ります．

グラフの実装

今回の実装では，以下のような形でグラフを実装しています．
辺は「始点」「終点」「容量」「逆辺」の4つを属性として保持しており，「逆辺」は残余ネットワークにおける辺を表現しています．また，今回は「流量」と「容量」を兼ねているため，破壊的な実装になります．（アルゴリズム前にコピーしておけば保持はできます．）
また，Network.addNode()のような一部のメソッドに関しては，今回の記事では使用しないため省略しています．

class Edge:
    def __init__(self, from_node, to_node, capacity):
        self.from_node = from_node
        self.to_node = to_node
        self.capacity = capacity

    def __repr__(self):
        return '{}->{} capacity:{}'.format(self.from_node, self.to_node, self.capacity)

class Network:
    def __init__(self):
        # 隣接リスト(key: Node Label, value: Edge List)
        self.adj_list = {}

    def addEdge(self, from_node, to_node, capacity=0):
        new_edge = Edge(from_node, to_node, capacity)
        new_edge_rev = Edge(to_node, from_node, 0)
        new_edge.rev = new_edge_rev
        new_edge_rev.rev = new_edge
        self.adj_list[from_node].append(new_edge)
        self.adj_list[to_node].append(new_edge_rev)

    def getEdge(self, from_node, to_node=None):
        if to_node is None:
            # to_nodeが指定されていない場合は，全てのedgeを返す
            return self.adj_list[from_node]
        else:
            # to_nodeが指定されている場合は，該当するedgeを返す
            for edge in self.adj_list[from_node]:
                if edge.to_node == to_node:
                    return edge

    def getNodeLabels(self):
        return [label for label in self.adj_list.keys()]

Push 操作

ある辺 $e$ が以下の条件を満たしている時，プリフローを更新する操作が可能になります．これを Push と呼びます．

• $e$ の始点が活性頂点である
• $e$ が可能辺である
• ($e$ の容量) $ > 0$

対象の辺 $e$ の容量がいっぱいになるまで，プリフローを流します．
(辺 $e$ のプリフロー) $ \leftarrow $ (辺 $e$ のプリフロー) $+ \min($ 辺 $e$ の始点の残存量, 辺 $e$ の残り容量 $)$

もちろん，残余ネットワークのフロー流量も，この操作で更新します．

def can_push(g, from_node, to_node, d_label):
    # Push操作が可能であるかの判定
    return d_label[from_node] == d_label[to_node] + 1 and g.getEdge(from_node, to_node).capacity > 0

def push(g, from_node, to_node, d_label, excess_flow):
    edge = g.getEdge(from_node, to_node)
    flow = min(excess_flow[from_node], edge.capacity)
    edge.capacity -= flow
    edge.rev.capacity += flow
    excess_flow[edge.from_node] -= flow
    excess_flow[edge.to_node] += flow

Relabel 操作

距離ラベルの更新操作を Relabel と呼びます．ある頂点 $u$ について，残存量がまだ存在するのにも関わらず，自身よりも低い頂点が近傍にない場合に限り実行できます．
更新は，以下のように行います．
頂点 $u$ の距離ラベル $\leftarrow \min($近傍にある頂点の距離ラベル$) + 1$

def relabel(g, node, d_label):
    d_label[node] = min([d_label[edge.to_node] for edge in g.getEdge(node) if edge.capacity > 0]) + 1

全体の処理

全体の処理になります．ここまでに説明した処理をwhileループ内で行っています．
また，活性頂点の選択・可能辺の選択方法については色々とあるのですが，今回は半ば思考停止的に，先頭から調べていき最初に当たった頂点・辺を選択するような実装にしています．

def preflow_push(g, source=None, target=None):
    # Set source node and target node
    if source is None:
        source = 's'
    if target is None:
        target = 't'

    # ステップ1
    # 距離ラベル(高さ)の初期化
    d_label = {node_label: 0 for node_label in g.getNodeLabels()}
    d_label[source] = len(g.getNodeLabels())
    # 過剰フロー(残存量)の初期化
    excess_flow = {node_label: 0 for node_label in g.getNodeLabels()}
    excess_flow[source] = np.inf

    initialize(g, source, d_label, excess_flow)

    while True:
        # 活性頂点の選択 (ステップ2)
        node_label = chooseActiveNode(g, source, target, excess_flow)
        if node_label is None:
            return sum([edge.rev.capacity for edge in g.getEdge(source)])

        # 可能辺の選択 (ステップ3)
        edge = chooseActiveEdge(g, node_label, d_label)

        if edge is None:
            # Relabel操作 (ステップ4に行けなかった場合)
            relabel(g, node_label, d_label)
        else:
            # Push操作 (ステップ4)
            push(g, edge.from_node, edge.to_node, d_label, excess_flow)


def chooseActiveNode(g, source, target, excess_flow):
    # 活性頂点の選択
    for node_label in g.getNodeLabels():
        if node_label == source or node_label == target:
            continue
        if excess_flow[node_label] > 0:
            return node_label
    return None


def chooseActiveEdge(g, node_label, d_label):
    # 可能辺の選択
    for edge in g.getEdge(node_label):
        if can_push(g, edge.from_node, edge.to_node, d_label):
            return edge
    return None


def initialize(g, source, d_label, excess_flow):
    # 初期化（ソースからフローを容量いっぱいまで流す）
    for edge in g.getEdge(source):
        push(g, source, edge.to_node, d_label, excess_flow)

おわりに

長くなりましたが，説明は以上です．
最近(?)，Qiitaで競プロの話題をよく見かけるような気がしますが，アルゴリズムの復習・実装のためにもまったり始めてみたいなーとか思うようになりました．参考文献にも1冊挙げていますが，グラフを実装するのに，非常に分かりやすく読みやすかったです．

参考文献

以下の3冊を参考にさせて頂きました．
- 「アルゴリズムイントロダクション 第2巻」
- 「プログラミングコンテストチャレンジブック」
- 「劣モジュラ最適化と機械学習」