Dijkstra法による最小費用流を求める実装（Primal-Dual）とポテンシャルに関する個人的なメモ

前提

基本的な考えは、Bellman-Ford法による最小費用流の求め方と同じである。

以下の実装では高速化のためBellman-Ford法の実装個所をDijkstra法に置き換えており、同時に負の辺が存在するグラフに対してDijkstra法を適用するため、ポテンシャルを導入している。

計算量: O(flow|E|log|V|)（最大流flow、頂点数V、辺数Eとする）

ポテンシャルについて

※ 以下の説明は、niuezさんのブログ（参考1）を参考に自分の理解を整理するために書いたものです。

最短距離の数式化

単一始点最短経路を求め済みのグラフにおいて、頂点iの最短距離をd[i]、頂点(i,j)の辺のコストをdist(i,j)とすると、辺が存在する頂点(i, j)間には以下のことが成り立つ。（これが最短距離を数式として落とし込んだものになる）

• d[s] = 0 とする。（sは始点)
• d[i] + dist(i, j) >= d[j] （iが始点に近く、jが終点に近い頂点とする）

ポテンシャルを導入し、グラフ（辺のコスト）を再構築する

各頂点にポテンシャルhを決めて、辺の距離を dist'(i, j) = dist(i, j) + h[i] - h[j] に置き換えたグラフG'を考える。 G'における各頂点の最短距離をd'とする。G'にも、最短距離の数式が成り立ち、更に2行目の様に式変換ができる。

d'[i] + dist'(i, j) >= d'[j]
⇔ d'[i] + h[i] + dist(i, j) >= d'[j] + h[j]

2行目から、d[i] = d'[i] + h[i] と見れば、d[i]は最短距離の式の形をしており、d[s] = d'[s] + h[s] より各頂点iについて d'[i] + h[i] - h[s] を計算すれば元の最短距離が求められるとわかる。

各辺を非負とするポテンシャルの求め方

ここまでで、ポテンシャルを導入して元のグラフの最短距離を求められることが分かった。次に、Dijkstra法を適用できるようにするためdist'を非負にしたいので、ポテンシャルに適切な値を設定することを考える。

dist(i, j) + h[i] - h[j] >=0
⇔ h[i] + dist(i, j) >= h[j]

2行目から最短距離の式の形をしているため、始点からの最短距離がdist'を非負にするためのポテンシャルであるとわかる。つまり1度Bellman-Ford法等により、各頂点の最短距離を求めておけば、以降その値をポテンシャルとすることでDijkstra法によっても最短距離を求めることができるようになる。

まとめ

1. 一度Bellman-Ford法等によって最短距離を求められればポテンシャルも求まる
2. ポテンシャルから非負の辺のみからなるグラフG'を構築でき、これに対してDijkstra法を適用できる
3. 2で求めたG'における各頂点の最短距離d'から始点sのポテンシャルを減算すれば元のグラフにおける最短距離を求められる

ちなみに以下の最小費用流を求める実装においては、Dijkstra法によって1度目の最短距離を求めている。これは、初期状態のグラフにおいて移動可能な経路上に負の辺が存在しないことを前提としているため。

参考

参考1）ポテンシャルについて

niuez.hatenablog.com

参考2）負閉路が無い事とポテンシャルが存在することは同値であることの説明

検証用問題

onlinejudge.u-aizu.ac.jp

実装

template <typename capacity_t, typename cost_t>
class PrimalDual {
private:
  const cost_t INF = 1e9;
  struct Edge {
    int to, rev;
    capacity_t cap;
    cost_t cost;
    Edge(int to, int rev, capacity_t cap, cost_t cost) : to(to), rev(rev), cap(cap), cost(cost) {}
  };
  vector<vector<Edge>> G;
  vector<int> prev_v, prev_e;
  Edge &get_rev(Edge &edge) { return G[edge.to][edge.rev]; }

public:
  PrimalDual(int n) {
    G.resize(n);
    prev_v.resize(n);
    prev_e.resize(n);
  }
  void add_edge(int from, int to, capacity_t cap, cost_t cost) {
    G[from].push_back(Edge(to, (int)G[to].size(), cap, cost));
    G[to].push_back(Edge(from, (int)G[from].size() - 1, 0, -cost));
  }

  cost_t get_min_cost_flow(int s, int t, capacity_t flow) {
    int n = G.size();
    vector<cost_t> h(n); // ポテンシャル
    cost_t res = 0;
    while (flow > 0) {
      vector<cost_t> dist(n, INF);
      dist[s] = 0;
      priority_queue<pair<cost_t, int>, vector<pair<cost_t, int>>, greater<pair<cost_t, int>>> que;
      que.emplace(0, s);
      while (!que.empty()) {
        cost_t c;
        int v;
        tie(c, v) = que.top();
        que.pop();
        if (c > dist[v]) continue;
        for (int i = 0; i < G[v].size(); i++) {
          auto &edge = G[v][i];
          if (edge.cap > 0 and dist[edge.to] > dist[v] + edge.cost + h[v] - h[edge.to]) {
            dist[edge.to] = dist[v] + edge.cost + h[v] - h[edge.to];
            prev_v[edge.to] = v;
            prev_e[edge.to] = i;
            que.emplace(dist[edge.to], edge.to);
          }
        }
      }
      if (dist[t] == INF) return -1;
      for (int v = 0; v < n; v++) h[v] += dist[v];
      capacity_t d = flow;
      for (int v = t; v != s; v = prev_v[v]) {
        d = min(d, G[prev_v[v]][prev_e[v]].cap);
      }
      flow -= d;
      res += d * h[t];
      for (int v = t; v != s; v = prev_v[v]) {
        G[prev_v[v]][prev_e[v]].cap -= d;
        get_rev(G[prev_v[v]][prev_e[v]]).cap += d;
      }
    }
    return res;
  }
};