容量スケーリング法のすゝめ

要約

最小費用流問題に対する最短路反復法(Successive Shortest Path, SSP)や Primal-Dual 法¹と呼ばれるアルゴリズムは, 少し弄ると弱多項式時間アルゴリズムにできる.
具体的には, $O(\U \cdot \mathrm{SP_+}(n, m, nC))$ が $O(m \log \U \cdot \mathrm{SP_+}(n, m, nC))$ になる.
$\mathrm{SP_+}(n, m, nC)$ は$n$ 頂点 $m$ 辺, 費用高々 $nC$ のグラフの一点から全点への最短路問題を解くのにかかる時間. 以下では実装をサボったダイクストラ法なので $O(m \log m)$ だが, $\exists k: m = O(n^k)$ を仮定すると合計 $O(m^2 \log \U \log n)$.

まえがき

ぼくの考えたさいきょうのフローライブラリでは, その速さは特に追求せず, SSP を用いたライブラリの設計例を紹介した. この記事では, 前回完成したライブラリを少し変更することで, 計算量を改善する.

よく教科書に乗っている容量スケーリング法とは, ちょっとだけ形を変えた実装になっている. 書いた後に参考文献を漁っていたら, この形の実装だと多項式時間な計算量が証明されていない とか見つけたので, もし僕の証明が間違っていたらごめんなさい. この形で多項式時間を示している文献を知っている方, もしくは間違えているところを見つけ方は, そっと報告してくれると喜びます. 正しいことが (僕以外の人の手によって) 知られている前処理方法についても, 最後に言及します.

実装のみを知りたい方は, 前半の章は読み飛ばして下さい.

背景知識

計算量解析で必要なものの準備を兼ねて, LP と容量スケーリング法のかかわりと, 容量スケーリング法の invariant を見る.

記号の準備

問題の定義すらまだだけれど, 先においておいた方が後で見やすいと思うので.

$\exE$ は $E$ に逆辺を追加した集合
$u_\vecf(e)$ は残余容量 (順辺に対しては $u_\vecf(e) := u_e - f_e$, 逆辺については $u_\vecf(\reve) := f_e - l_e$)
$c^\vecp(e)$ は簡約コスト(順辺 $e = (u, v)$ に対しては $c^\vecp(e) := c_e + p_u - p_v$, 逆辺については $c^\vecp(\reve) := -c^\vecp(e)$)
$G_\vecf$ は $\vecf$ に関する残余グラフ ($\exE$ のうち $u_\vecf(e) > 0$ なものを集めたもの. 容量は $u_\vecf(e)$.)
$G^\Delta_\vecf$ は $G_\vecf$ のうち容量が $\Delta$ 以上のもののみからなるグラフ
$\U := \max_e \{ u_e - l_e \}$
$\delta^+(v)$ は $v$ から出る $G$ の辺全体の集合, $\delta^-(v)$ は $v$ へ入る $G$ の辺全体の集合.
$e_\vecf(v) := b_v - (\sum_{e \in \delta^+(v)} f_e - \sum_{e \in \delta^-(v)} f_e)$ (余剰量)
$S_\vecf := \{ v | e_\vecf(v) > 0 \}$, $T_\vecf := \{ v | e_\vecf(v) < 0 \}$
$S^\Delta_\vecf := \{ v | e_\vecf(v) \ge +\Delta \}$, $T^\Delta_\vecf := \{ v | e_\vecf(v) \le -\Delta \}$

LP と Primal-Dual 法

最小費用流問題²は, その線型計画問題としての表現と強双対性を使うと, 次のような問題と捉えられる;

与えられた有向グラフ $G = (V, E)$ と, $\sum_v b_v = 0$ であるような $\vec{b} = \{b_v\}_{v \in V}$, $l_e \le u_e$ であるような $\vec{l} = \{l_e\}_{e \in E}$ と $\vec{u} = \{u_e\}_{e \in E}$, $\vec{c} = \{c_e\}_{e \in E}$ の組 $(G, \vec{b}, \vec{l}, \vec{u}, \vec{c})$ に対し, $\vec{f} = \{f_e\}_{e \in E}$, $\vec{p} = \{p_v\}_{v \in V}$ で, 次を満たすものを求めよ.

$l_e \leq f_e \leq u_e, \quad \forall e \in E$ (容量制約)
$\sum_{e \in \delta^+(v)} f_e - \sum_{e \in \delta^-(v)} f_e = b_v, \quad \forall v \in V$ (フロー整合性条件)
$f_e > l_e \Rightarrow c_e + p_u - p_v \leq 0, \quad \forall e = (u, v) \in E$ (相補性条件)
$f_e < u_e \Rightarrow c_e + p_u - p_v \geq 0, \quad \forall e = (u, v) \in E$ (相補性条件)

このフロー整合性条件は, 主問題が infeasible であった時のことを考慮し, 残余グラフ $G_\vecf$ と余剰量 $e_\vecf(v)$ を用いて次のように書き直しておく.

$G_\vecf$ 上に $S_\vecf$ の頂点から $T_\vecf$ の頂点へのパス(増大路)が存在しない. (カット条件³)

フロー整合性条件を満たすフローに対しては, 全頂点で $e_\vecf(v) = 0$ だからカット条件は成立する一方, フロー整合性条件を満たさないがカット条件を満たすフローが存在するならば, 主問題が infeasible であることが確かめられる.

また, 残余容量 $u_\vecf(e)$ と簡約コスト $c^\vecp(e)$, 逆辺を追加した辺集合 $\exE$ を用いると, 相補性条件は次の一つの式で表せる.

$u_\vecf(e) > 0 \Rightarrow c^\vecp(e) \geq 0, \quad \forall e = (u, v) \in \exE$ (相補性条件)

Primal-Dual 法は, 容量制約と相補性条件を満たす主双対ペア $(\vec{f}, \vec{p})$ を持ち, カット条件が満たされるまで解を改善するアルゴリズムである. 具体的には, 初期解を得た後, 次の2つのステップを交互に繰り返すアルゴリズムである.

$\vec{f}$ は固定しておき, 相補性条件を保ちつつ $\vec{p}$ を変更することで, 双対目的関数値を改善する.
$\vec{p}$ は固定しておき, 相補性条件を保ちつつ $\vec{f}$ を変更することで, フローをより整合させる. すなわち $G_\vecf$ 上での $S_\vecf$ 全体から $T_\vecf$ 全体へのカット容量を小さくする.

相補性条件の緩和

Primal-Dual 法は相補性条件を常に保っていた. 容量スケーリング法では, 一部を $\Delta$-緩和した, 次の3つの条件を考える.

$l_e \leq f_e \leq u_e, \quad \forall e \in E$ (容量制約)
$G^\Delta_\vecf$ 上に $S^\Delta_\vecf$ の頂点から $T^\Delta_\vecf$ の頂点へのパス($\Delta$-増大路)が存在しない ($\Delta$-緩和されたカット条件)
$u_\vecf(e) \geq \Delta \Rightarrow c^\vecp(e) \geq 0, \quad \forall e = (u, v) \in \exE$ ($\Delta$-緩和された相補性条件)

ここで, $G^\Delta_\vecf$ は $G_\vecf$ から容量が $\Delta$ 未満の辺を削除したものである.

容量が整数という条件の下, $\Delta = 1$ の時に元の問題と一致する. また, $\Delta$ が充分に大きいとき, 具体的には $\U := \max_e \{ u_e - l_e \}$ より大きいとき, 例えば $f_e = l_e, p_v = 0$ が上の条件を満たす.

容量スケーリング法は, 充分大きな $\Delta$ を取り, この $\Delta$-緩和された問題を解き, 得られた解を $\Delta/2$-緩和された問題の "充分に良い" 初期解に変形することを, $\Delta = 1$ になるまで繰り返すアルゴリズムである. $\Delta$-緩和された問題の初期解が与えられた時, これを解くアルゴリズムは, Primal-Dual 法と同様, 次の2つのステップを交互に繰り返す. これを $\Delta$-scaling phase と呼ぶ.

$\vec{f}$ は固定しておき, 緩和された相補性条件を保ちつつ $\vec{p}$ を変更することで, 双対目的関数値を改善する.
$\vec{p}$ は固定しておき, 緩和された相補性条件を保ちつつ $\vec{f}$ を変更することで, フローをより整合させる. すなわち $G^\Delta_\vecf$ 上での $S^\Delta_\vecf$ 全体から $T^\Delta_\vecf$ 全体へのカット容量を小さくする.

補足

上では主問題側の制約を $\Delta$-緩和した. 一方で双対問題側の制約を $\varepsilon$-緩和すると, 例えば

$u_\vecf(e) > 0 \Rightarrow c^\vecp(e) \geq -\varepsilon, \quad \forall e = (u, v) \in \exE$ ($\varepsilon$-緩和された相補性条件)

となる. これを用いると, コストスケーリング法 (Cost Scaling Algorithm) や平均長最小負閉路解消法 (Minimum Mean Cycle Cancel Algorithm) などが得られる.

具体的なアルゴリズム

引き続き, 以下の記号を使う.

$\exE$ は $E$ に逆辺を追加した集合
$u_\vecf(e)$ は残余容量 (順辺に対しては $u_\vecf(e) := u_e - f_e$, 逆辺については $u_\vecf(\reve) := f_e - l_e$)
$c^\vecp(e)$ は簡約コスト(順辺 $e = (u, v)$ に対しては $c^\vecp(e) := c_e + p_u - p_v$, 逆辺については $c^\vecp(\reve) := -c^\vecp(e)$)
$G_\vecf$ は $\vecf$ に関する残余グラフ ($\exE$ のうち $u_\vecf(e) > 0$ なものを集めたもの. 容量は $u_\vecf(e)$.)
$G^\Delta_\vecf$ は $G_\vecf$ のうち容量が $\Delta$ 以上のもののみからなるグラフ
$\U := \max_e \{ u_e - l_e \}$
$e_\vecf(v) := b_v - (\sum_{e \in \delta^+(v)} f_e - \sum_{e \in \delta^-(v)} f_e)$ (余剰量)
$S_\vecf := \{ v | e_\vecf(v) > 0 \}$, $T_\vecf := \{ v | e_\vecf(v) < 0 \}$
$S^\Delta_\vecf := \{ v | e_\vecf(v) \ge +\Delta \}$, $T^\Delta_\vecf := \{ v | e_\vecf(v) \le -\Delta \}$

今回紹介する容量スケーリング法は, 次のようなアルゴリズムである.

容量スケーリング法:

$\Delta = 2^{\floor{\log_2 \U} + 1}$ とする.
$\vec{p} = \veczero$, $\vec{f} = \vec{l}$ とする.
$\Delta = 1$ なら終了し, そうでないならば $\Delta := \Delta / 2$ と更新する.
$c^\vecp(e) < 0$ な $G^\Delta_\vecf$ の辺に対し, $u_\vecf(e) - \Delta$ より多く $u_\vecf(e)$ 以下だけ, つまり残余容量が $\Delta$ 未満になるように push する.
"$\Delta$-scaling phase" を行う.
3.に戻る.

ここで, "push する" は, $f_e$ の値を増やすの意味. これに伴い, 順辺と逆辺残余容量 $u_\vecf$, この辺の端点の余剰量 $e_\vecf$ が変更されることに注意. このアルゴリズムは, 次の不変条件を満たす.

2.の初期化が終わった時点以降, $\vecf$ は常に容量制約を満たす.
3.に(2. もしくは 6. から)突入した時点で, $(\vecf, \vecp)$ は $\Delta$-緩和されたカット条件を満たす
5.の(突入時, 脱出時を含めた)実行中, $(\vecf, \vecp)$ は $\Delta$-緩和された相補性条件を満たす.
5.の実行中, $S^\Delta_\vecf$, $T^\Delta_\vecf$ に頂点が新たに加わることはない.

例えば, 1. と 2. の初期化が終わった時点で, 任意の辺の残余容量は $\Delta$ 未満だから, この $(\vecf, \vecp)$ は $\Delta$-緩和されたカット条件と $\Delta$-緩和された相補性条件を満たす. また, 4. の後 $c^\vecp(e) < 0 \Rightarrow u_\vecf(e) < \Delta$ が満たされるが, これは $\Delta$-緩和された相補性条件の対偶である.

$\Delta$-scaling phase は, 次の $\Delta$-scaling dual step と$\Delta$-scaling primal step を, $G^\Delta_\vecf$ に $S^\Delta_\vecf$ から $T^\Delta_\vecf$ へのパスが無くなるまで, すなわち $\Delta$-緩和されたカット条件を満たすまで, 交互に行う.

$\Delta$-scaling dual step:

$G^\Delta_\vecf$ 上で $S^\Delta_\vecf$ から全点へのコスト $c^\vecp(e)$ による (多始点)最短路長 $\vec{d}$ を求める.
$\vecp := \vecp + \vec{d}$ と更新する.

これにより, $G^\Delta_\vecf$ 上での $S^\Delta_\vecf$ からの最短経路に入りうる辺全体で $c^\vecp(e) = 0$ となる. $S^\Delta_\vecf$ から $T^\Delta_\vecf$ へのパスが存在したから, その中で最短なものの辺が $c^\vecp(e) = 0$ を満たすようになった. 一方, この更新後も任意の $G^\Delta_\vecf$ の辺で $c^\vecp(e) \ge 0$ であるから, $(\vecf, \vecp)$ は $\Delta$-緩和された相補性条件を満たす.

$\Delta$-scaling primal step:

$S^\Delta_\vecf$ の頂点から $T^\Delta_\vecf$ の頂点へ $c^\vecp(e) = 0$ な $G^\Delta_\vecf$ の辺からなるパスを一つ取る.
そのパスに沿ってフローを流す. 見つかったのが $s$-$t$ path $P$ であったとすると, $\Delta$ 以上 $\min\{ e_\vecf(s), -e_\vecf(t), \min_{e \in P} u_\vecf(e) \}$ 以下だけ流す.

この更新は $c^\vecp(e) = 0$ な辺上で行われるから, $\Delta$-緩和された相補性条件を保存する. また, この更新により $S^\Delta_\vecf$ と $T^\Delta_\vecf$ に頂点が新たに加わることはない.

実装

ぼくの考えたさいきょうのフローライブラリから, 前処理, メインループ, primal step, dual step を変更する.

新しい前処理

前処理では, 容量制約を守る $\vecf$ を一つとり, $\Delta$ の初期値を決める.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
  std::pair<Status, Cost> solve() {
    potential.resize(n);
    for (auto &es : g) for (auto &e : es) {
      const Flow rcap = e.residual_cap();
      if (rcap < 0) {
        push(e, rcap);
        b[e.src] -= rcap;
        b[e.dst] += rcap;
      }
    }

    Flow inf_flow = 1;
    for (const auto &es : g) for (const auto &e : es) inf_flow = std::max(inf_flow, e.residual_cap());
    Flow delta = 1;
    while (delta <= inf_flow) delta *= SCALING_FACTOR;
    // 続く
  }

メインループと各フェーズの前処理

メインループは次のようになる.

1
2
3
4
    for (delta /= SCALING_FACTOR; delta; delta /= SCALING_FACTOR) {
      saturate_negative(delta);
      while (dual(delta)) primal(delta);
    }

ここで, saturate_negative は次のような関数である.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
  void saturate_negative(const Flow delta) {
    excess_vs.clear();
    deficit_vs.clear();
    for (auto &es : g) for (auto &e : es) {
      const Flow rcap = e.residual_cap();
      const Cost rcost = residual_cost(e.src, e.dst, e);
      if (rcost < 0 && rcap >= delta) {
        push(e, rcap);
        b[e.src] -= rcap;
        b[e.dst] += rcap;
      }
    }
    for (V_id v = 0; v < n; ++v) if (b[v] != 0) {
      (b[v] > 0 ? excess_vs : deficit_vs).emplace_back(v);
    }
  }

dual step

ほぼ, Primal-Dual におけるフロー型の 0 との比較を $\Delta$ との比較に変えただけ.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
  bool dual(const Flow delta) {
    dist.assign(n, unreachable);
    parent.assign(n, nullptr);
    excess_vs.erase(std::remove_if(std::begin(excess_vs), std::end(excess_vs),
                                   [&](const V_id v) { return b[v] < delta; }),
                    std::end(excess_vs));
    deficit_vs.erase(std::remove_if(std::begin(deficit_vs),
                                    std::end(deficit_vs),
                                    [&](const V_id v) { return b[v] > -delta; }),
                     std::end(deficit_vs));
    for (const auto v : excess_vs) pq.emplace(dist[v] = 0, v);
    farthest = 0;
    std::size_t deficit_count = 0;
    while (!pq.empty()) {
      const auto [d, u] = pq.top();
      pq.pop();
      if (dist[u] < d) continue;
      farthest = d;
      if (b[u] <= -delta) ++deficit_count;
      if (deficit_count >= deficit_vs.size()) break;
      for (auto &e : g[u]) {
        if (e.residual_cap() < delta) continue;
        const auto v = e.dst;
        const auto new_dist = d + residual_cost(u, v, e);
        if (new_dist >= dist[v]) continue;
        pq.emplace(dist[v] = new_dist, v);
        parent[v] = &e;
      }
    }
    pq = decltype(pq)(); // pq.clear() doesn't exist.
    for (V_id v = 0; v < n; ++v) {
      potential[v] += std::min(dist[v], farthest);
    }
    return deficit_count > 0;
  }

primal step

こちらも同様. $\Delta$ 以上しか流さないようにしていることに注意.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
  void primal(const Flow delta) {
    for (const auto t : deficit_vs) {
      if (dist[t] > farthest) continue;
      Flow f = -b[t];
      V_id v;
      for (v = t; parent[v] != nullptr && f >= delta; v = parent[v]->src) {
        f = std::min(f, parent[v]->residual_cap());
      }
      f = std::min(f, b[v]);
      if (f < delta) continue;
      for (v = t; parent[v] != nullptr;) {
        auto &e = *parent[v];
        push(e, f);
        const size_t u = parent[v]->src;
        parent[v] = nullptr;
        v = u;
      }
      b[t] += f;
      b[v] -= f;
    }
  }

計算量

primal step は $O(m)$, dual step は $O(m \log m)$ である.

$\U$ を $\max_e \{ u_e - l_e \}$ とすると, $\Delta$ の初期値は $O(\U)$ であり, 合計 $O(\log \U)$ 回の scaling phase がある.

$\Delta$-scaling phase で流すフロー量を $F_\Delta$ とする. フェーズ末尾を除く各 primal/dual step につき, 少なくとも $\Delta$ だけフローを流すから, 上の実装は $O ( \sum_\Delta (F_\Delta / \Delta) m \log m )$ になる. この $F_\Delta$ について考える.

大雑把に言うと, 次の2つを示す:

フェーズの前処理で追加されるフローの不整合度は $O(m \Delta)$ である.
前のフェーズから残っているフローの不整合度のうち, 解決出来るのは $O(m \Delta)$ である.

さて, $\Delta$ を一つ固定し, $\Delta$-scaling phase の前処理を行う直前 ($2\Delta$-scaling phase の直後又は初回) のフローを $\vecf'$, 前処理の直後のフローを $\vecf$ とし, (前処理でポテンシャルは不変なので両方に)対応するポテンシャルを $\vecp$ とすると,

$G^{2\Delta}_{\vecf'}$ 上に $S^{2\Delta}_{\vecf'}$ の頂点から $T^{2\Delta}_{\vecf'}$ の頂点へのパス($2\Delta$-増大路)が存在しない ($2\Delta$-緩和されたカット条件)
$c^\vecp(e) < 0 \Rightarrow u_{\vecf'}(e) < 2 \Delta, \quad \forall e = (u, v) \in \exE$ ($2 \Delta$-緩和された相補性条件の対偶)

が成り立つ. 各 primal step では $S^\Delta_\vecf$ の頂点から $T^\Delta_\vecf$ の頂点にフローを流す. この始点/終点になりうる頂点を4つに分割し,

$A = S^{2\Delta}_{\vecf'} \cap S^\Delta_\vecf$,
$B = T^{2\Delta}_{\vecf'} \cap T^\Delta_\vecf$,
$C = S^\Delta_\vecf \setminus (A \cup B)$,
$D = T^\Delta_\vecf \setminus (A \cup B)$

とすると, $F_\Delta$ は $A \cup C$ から $B \cup D$ に流れることになる.

$\Delta$-scaling phase の前処理は $c^\vecp(e) < 0$ な辺の残余容量が無くなるように各辺 $e$ に push するから, $(\vecf', \vecp)$ に対する $2\Delta$-緩和された相補性条件から, 各辺 $e$ で $|f_e - f'_e| < 2 \Delta$ であり, $\sum_v |e_\vecf(v) - e_{\vecf'}(v)| < 4 m \Delta$ である.

また, $s \in C$ に対し, $e_\vecf(s) \ge \Delta > 0$ であるが, $s \not \in S^{2\Delta}_{\vecf'}$, すなわち $e_{\vecf'}(s) < 2 \Delta$ であったから, $$|e_\vecf(s)| = e_\vecf(s) = e_\vecf(s) - e_{\vecf'}(s) + e_\vecf'(s) < |e_\vecf(s) - e_{\vecf'}(s)| + 2 \Delta$$ である. 同様に $t \in D$ に対しても $|e_\vecf(t)| < |e_\vecf(t) - e_{\vecf'}(t)| + 2 \Delta$ であるから, $C$ と $D$ に対する $|e_\vecf(v)|$ の合計は高々 $\sum_v |e_\vecf(v) - e_{\vecf'}(v)| + 2 \Delta |C \cup D| < 4m\Delta + 2n\Delta = O(m \Delta)$ である. 従って, $A$ から $D$, $C$ から $B$, $C$ から $D$ に流すフロー量の合計は $O(m \Delta)$.

一方, $(\vecf', \vecp)$ に対する $2\Delta$-緩和されたカット条件から, $A$ から $B$ は $G^{2 \Delta}_{\vecf'}$ 上にパスが存在しない. 従って, $G_{\vecf'}$ 上での最小 $A$-$B$ カットの容量は高々 $2 m \Delta$ である. $|f_e - f'_e| < 2 \Delta$ であったから, $G_{\vecf}$ での最小 $A$-$B$ カットの容量も $O(m \Delta)$ であり, $A$ から $B$ に流せるフロー量の合計は $O(m \Delta)$.

以上を纏めて, $F_\Delta = O(m \Delta)$, 従って全フェーズ合わせて $O ( \sum_\Delta (F_\Delta / \Delta) m \log m ) = O(\sum_\Delta m^2 \log m) = O(m^2 \log m \log \U)$ である.

scaling factor を $2$ ではなく $\alpha$ とすると, $O(m^2 \alpha \log m \log_\alpha \U)$ になる.

別の方針

実質的な内容は同じであるが, 人によっては次の方針が $A, B, C, D$ への分割を考えるより明瞭に思えるかもしれない.

$\vecf$ に対し, $V \cup \{s, t\}$ を頂点とするグラフ $H_\vecf$ で, $\overleftrightarrow{E}$ に残余容量で辺容量を入れた辺と, $s$ から $v$ に容量 $\max\{0, e_\vecf(v)\}$ の辺, $v$ から $t$ に容量 $\max\{0, -e_\vecf(v)\}$ の辺を持つものを考える.
$G_\vecf$ 上の $S^\Delta_\vecf$ から $T^\Delta_\vecf$ への $\Delta$-残余パスは $H_\vecf$ 上の $s$ から $t$ への $\Delta$-残余パスと一対一対応し, $F_\Delta$ は $\Delta$-scaling phase 開始時のフロー $\vecf$ に対応する $H_\vecf$ の最大 $s$-$t$ フロー量で上から抑えられる.
$\Delta$-scaling phase の前処理の直前のフロー $\vecf'$ に対応する $H_{\vecf'}$ の最小カットの容量は ($2\Delta$-カット条件から) $O(m \Delta)$.
$H_{\vecf'}$ と $H_\vecf$ で, 各辺の容量の差の絶対値の和は $O(m \Delta)$. これは元のグラフの辺容量だけでなく $e_\vecf$ の変化も含むことに注意.
従って, $H_\vecf$ の最小 $s$-$t$ カット容量も $O(m \Delta)$ で, これを言えばよかった.

注意

Z. Király and P. Kovács (2012) の Capacity Scaling (CAS) は恐らく上記のアルゴリズムとほぼ同じであるが, これはよく知られた Capacity Scaling Algorithm から前処理を少し削っている. そして "In this case, the polynomial running time bound is not proved" とあるので, 上の僕の証明は間違っているかもしれない... 不備を見つけたら, 連絡をいただけるとありがたいです.

さて, 教科書的な容量スケーリング法では, 任意の頂点対 $(u, v)$ に対し, $u$ から $v$ への充分大きな容量のパスがあることを仮定する (これは, 適当な頂点を一つ選び, 任意の頂点との間にコスト $\infty$ の辺を双方向に追加すれば容易に満たせる). また, $\Delta$ の倍数しか流さないようにすることが多いようだ. 充分容量の大きいパスがある仮定の下, 各 $\Delta$-scaling phase は $S^\Delta_\vecf$ か $T^\Delta_\vecf$ が空になることにより終了する事が示せる. これを用いると, 上でいう $A$ から $B$ に流れるフローを考える必要がなくなり, 証明がより簡単になる. 上の証明が正しく無さそうだと思う人は, この前処理を入れるとよいだろう.

参考文献

1.はちゃんと読んでいないが, 容量 $\infty$ なものについて扱い, その後輸送問題への変形を使って有限容量に対応しているようだ.
2.では 10.2 節で扱われている. ちゃんと書いていないが, $\Delta$-scaling phase の終了条件の正当性に Assumption 9.4 (上の注意で言った仮定)が効いている.
3.で上記の注意が出てくる. "証明されていないが, この実装の方が有意に速かったのでこっちを使った" という旨が書かれている.

謝辞

とこはるさんには事前に読んでいただき, 有用なコメントと共に安心感をいただきました. ありがとうございました!

ソースコード

長いので折りたたみ

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <numeric>
#include <queue>
#include <vector>

enum Objective {
    MINIMIZE = 1,
    MAXIMIZE = -1,
};
enum class Status {
    OPTIMAL,
    INFEASIBLE,
};

template <class Flow, class Cost, Objective objective = Objective::MINIMIZE, Flow SCALING_FACTOR = 2>
class MinCostFlow {
  using V_id = uint32_t;
  using E_id = uint32_t;

  class Edge {
    friend class MinCostFlow;

    V_id src, dst;
    Flow flow, cap;
    Cost cost;
    E_id rev;

  public:
    Edge() = default;

    Edge(const V_id src, const V_id dst, const Flow cap, const Cost cost,
         const E_id rev)
        : src(src), dst(dst), flow(0), cap(cap), cost(cost), rev(rev) {}

    [[nodiscard]] Flow residual_cap() const { return cap - flow; }
  };

public:

  class EdgePtr {
    friend class MinCostFlow;

    const MinCostFlow *instance;
    const V_id v;
    const E_id e;

    EdgePtr(const MinCostFlow *instance, const V_id v, const E_id e)
        : instance(instance), v(v), e(e) {}

    [[nodiscard]] const Edge &edge() const { return instance->g[v][e]; }

    [[nodiscard]] const Edge &rev() const {
      const Edge &e = edge();
      return instance->g[e.dst][e.rev];
    }

  public:
    [[nodiscard]] V_id src() const { return rev().dst; }

    [[nodiscard]] V_id dst() const { return edge().dst; }

    [[nodiscard]] Flow flow() const { return edge().flow; }

    [[nodiscard]] Flow lower() const { return -rev().cap; }

    [[nodiscard]] Flow upper() const { return edge().cap; }

    [[nodiscard]] Cost cost() const { return edge().cost; }

    [[nodiscard]] Cost gain() const { return -edge().cost; }
  };

private:
  V_id n;
  std::vector<std::vector<Edge>> g;
  std::vector<Flow> b;

public:
  MinCostFlow() : n(0) {}

  V_id add_vertex() {
    ++n;
    g.resize(n);
    b.resize(n);
    return n-1;
  }

  std::vector<V_id> add_vertices(const size_t size) {
    std::vector<V_id> ret;
    for (V_id i = 0; i < size; ++i) ret.emplace_back(n + i);
    n += size;
    g.resize(n);
    b.resize(n);
    return ret;
  }

  EdgePtr add_edge(const V_id src, const V_id dst, const Flow lower,
                   const Flow upper, const Cost cost) {
    const E_id e = g[src].size(), re = src == dst ? e + 1 : g[dst].size();
    assert(lower <= upper);
    g[src].emplace_back(Edge{src, dst, upper, cost * objective, re});
    g[dst].emplace_back(Edge{dst, src, -lower, -cost * objective, e});
    return EdgePtr{this, src, e};
  }

  void add_supply(const V_id v, const Flow amount) { b[v] += amount; }

  void add_demand(const V_id v, const Flow amount) { b[v] -= amount; }

private:
  // Variables used in calculation
  static Cost constexpr unreachable = std::numeric_limits<Cost>::max();
  Cost farthest;
  std::vector<Cost> potential;
  std::vector<Cost> dist;
  std::vector<Edge *> parent; // out-forrest.
  std::priority_queue<std::pair<Cost, int>, std::vector<std::pair<Cost, int>>,
                      std::greater<>>
      pq; // should be empty outside of dual()
  std::vector<V_id> excess_vs, deficit_vs;

  Edge &rev(const Edge &e) { return g[e.dst][e.rev]; }

  void push(Edge &e, const Flow amount) {
    e.flow += amount;
    g[e.dst][e.rev].flow -= amount;
  }

  Cost residual_cost(const V_id src, const V_id dst, const Edge &e) {
    return e.cost + potential[src] - potential[dst];
  }

  bool dual(const Flow delta) {
    dist.assign(n, unreachable);
    parent.assign(n, nullptr);
    excess_vs.erase(std::remove_if(std::begin(excess_vs), std::end(excess_vs),
                                   [&](const V_id v) { return b[v] < delta; }),
                    std::end(excess_vs));
    deficit_vs.erase(std::remove_if(std::begin(deficit_vs),
                                    std::end(deficit_vs),
                                    [&](const V_id v) { return b[v] > -delta; }),
                     std::end(deficit_vs));
    for (const auto v : excess_vs) pq.emplace(dist[v] = 0, v);
    farthest = 0;
    std::size_t deficit_count = 0;
    while (!pq.empty()) {
      const auto [d, u] = pq.top();
      pq.pop();
      if (dist[u] < d) continue;
      farthest = d;
      if (b[u] <= -delta) ++deficit_count;
      if (deficit_count >= deficit_vs.size()) break;
      for (auto &e : g[u]) {
        if (e.residual_cap() < delta) continue;
        const auto v = e.dst;
        const auto new_dist = d + residual_cost(u, v, e);
        if (new_dist >= dist[v]) continue;
        pq.emplace(dist[v] = new_dist, v);
        parent[v] = &e;
      }
    }
    pq = decltype(pq)(); // pq.clear() doesn't exist.
    for (V_id v = 0; v < n; ++v) {
      potential[v] += std::min(dist[v], farthest);
    }
    return deficit_count > 0;
  }

  void primal(const Flow delta) {
    for (const auto t : deficit_vs) {
      if (dist[t] > farthest) continue;
      Flow f = -b[t];
      V_id v;
      for (v = t; parent[v] != nullptr && f >= delta; v = parent[v]->src) {
        f = std::min(f, parent[v]->residual_cap());
      }
      f = std::min(f, b[v]);
      if (f < delta) continue;
      for (v = t; parent[v] != nullptr;) {
        auto &e = *parent[v];
        push(e, f);
        const size_t u = parent[v]->src;
        parent[v] = nullptr;
        v = u;
      }
      b[t] += f;
      b[v] -= f;
    }
  }

  void saturate_negative(const Flow delta) {
    excess_vs.clear();
    deficit_vs.clear();
    for (auto &es : g) for (auto &e : es) {
      const Flow rcap = e.residual_cap();
      const Cost rcost = residual_cost(e.src, e.dst, e);
      if (rcost < 0 && rcap >= delta) {
        push(e, rcap);
        b[e.src] -= rcap;
        b[e.dst] += rcap;
      }
    }
    for (V_id v = 0; v < n; ++v) if (b[v] != 0) {
      (b[v] > 0 ? excess_vs : deficit_vs).emplace_back(v);
    }
  }

public:
  std::pair<Status, Cost> solve() {
    potential.resize(n);
    for (auto &es : g) for (auto &e : es) {
      const Flow rcap = e.residual_cap();
      if (rcap < 0) {
        push(e, rcap);
        b[e.src] -= rcap;
        b[e.dst] += rcap;
      }
    }

    Flow inf_flow = 1;
    for (const auto &es : g) for (const auto &e : es) inf_flow = std::max(inf_flow, e.residual_cap());
    Flow delta = 1;
    while (delta <= inf_flow) delta *= SCALING_FACTOR;

    for (delta /= SCALING_FACTOR; delta; delta /= SCALING_FACTOR) {
      saturate_negative(delta);
      while (dual(delta)) primal(delta);
    }

    Cost value = 0;
    for (const auto &es : g) for (const auto &e : es) {
      value += e.flow * e.cost;
    }
    value /= 2;

    if (excess_vs.empty() && deficit_vs.empty()) {
      return { Status::OPTIMAL, value / objective };
    } else {
      return { Status::INFEASIBLE, value / objective };
    }
  }
};

template <class Flow, class Cost>
using MaxGainFlow = MinCostFlow<Flow, Cost, Objective::MAXIMIZE>;

使用例

Library-Checker: Minimum cost b-flow

問題はこちら

長いので折りたたみ

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
namespace methods_to_add {
  // 以下の3つの関数を MinCostFlow クラスに追加する
  template<class T>
  T get_result_value() {
    T value = 0;
    for (const auto &es : g) for (const auto &e : es) {
      value += (T)(e.flow) * (T)(e.cost);
    }
    value /= (T)2;
    return value / objective;
  }

  std::vector<Cost> get_potential() {
    // Not strictly necessary, but re-calculate potential to bound the potential values,
    // plus make them somewhat canonical so that it is robust for the algorithm chaneges.
    std::fill(std::begin(potential), std::end(potential), 0);
    for (size_t i = 0; i < n; ++i) for (const auto &es : g) for (const auto &e : es)
      if (e.residual_cap() > 0) potential[e.dst] = std::min(potential[e.dst], potential[e.src] + e.cost);
    return potential;
  }

  std::vector<size_t> get_cut() {
    std::vector<size_t> res;
    if (excess_vs.empty()) return res;
    for (size_t v = 0; v < n; ++v) {
      if (deficit_vs.empty() || (dist[v] < unreachable))
        res.emplace_back(v);
    }
    return res;
  }
}

#include <cstdint>
#include <cstdio>
#include <set>
#include <string>

#define REP(i, b, n) for (int i = (int)(b); i < (int)(n); ++i)
#define rep(i, n) REP(i, 0, n)
#define loop(n) rep(i##__COUNTER__, n)

int readI() {
  int n;
  scanf("%d", &n);
  return n;
}
long long readLL() {
  long long n;
  scanf("%lld", &n);
  return n;
}

template<class T>
std::string i2s(T value) {
  if (value < 0) return "-" + i2s(-value);
  if (value == 0) return "0";
  std::string s;
  while (value) {
    s += '0' + (value % 10);
    value /= 10;
  }
  std::reverse(s.begin(), s.end());
  return s;
}

int main(void) {
  using Flow = long long;
  using Cost = long long;
  using MCF = MinCostFlow<Flow, Cost>;
  const int n = readI();
  const int m = readI();
  MCF mcf;
  const auto vs = mcf.add_vertices(n);
  std::vector<Flow> original_bs(n);
  rep (v, n) {
    const Flow b = readLL();
    original_bs[v] = b;
    mcf.add_supply(vs[v], b);
  }
  std::vector<MCF::EdgePtr> edges;
  loop (m) {
    const int s = readI();
    const int t = readI();
    const Flow lower = readLL();
    const Flow upper = readLL();
    const Cost cost = readLL();
    edges.emplace_back(mcf.add_edge(s, t, lower, upper, cost));
  }
  const auto status = mcf.solve().first;
  if (status == Status::INFEASIBLE) {
    const auto cut_vec = mcf.get_cut();
    const std::set<size_t> cut_set(cut_vec.begin(), cut_vec.end());
    Flow left_sum = 0, right_sum = 0, cap_sum = 0;
    rep (v, n) (cut_set.count(v) ? left_sum : right_sum) += original_bs[v];
    for (auto &e : edges) {
      const auto sl = cut_set.count(e.src()) > 0;
      const auto tr = cut_set.count(e.dst()) == 0;
      if (sl != tr) continue;
      if (sl) {
        cap_sum += e.upper();
      } else {
        cap_sum -= e.lower();
      }
    }
    assert((left_sum > cap_sum) || (right_sum < -cap_sum));
    puts("infeasible");
  } else {
    const auto potential = mcf.get_potential();
    const auto result_value = mcf.get_result_value<__int128_t>();
    puts(i2s(result_value).c_str());
    for (const auto &v : vs) {
      puts(i2s(potential[v]).c_str());
    }
    for (const auto &e : edges) {
      puts(i2s(e.flow()).c_str());
    }
  }
}

厳密にはこの2つは異なる(例えば蜘蛛の巣本)が, この記事では深く言及しない. ↩
ここでは最小費用 $\vecb$-フロー問題を採用する. ↩
一般的にどう呼ばれるか知らないが, ここではカット条件と呼ぶことにする. ↩