AGC018: C – Coins を離散凸解析(L 凸関数最小化)で解く

最近,一部の競技プログラマの間で離散凸解析が流行っている(?)ようです.

そのブームに乗り,プログラミングコンテスト(AtCoder Grand Contest)で出題されたある問題の別解として,離散凸解析の L 凸関数最小化を使った解法を紹介します.記事の後半では,解法を理解・実装するために必要な離散凸解析の知識を説明します(分量的にはそちらがメインです).

想定読者は AtCoder 青色以上程度の方です.

問題概要

問題ページ: AtCoder Grand Contest 018: C - Coins

本質を見やすくするために,問題を数学的に書き直しておきましょう.

頂点集合が $U = \set{1,\dots,N},$ $V = \set{1,2,3}$ である完全二部グラフ $G = (U, V; E)$ を考える.以下が与えられる:

  • 正整数 $b_1, b_2, b_3$(ただし $b_1 + b_2 + b_3 = N$),
  • 各辺 $(i, j) \in E$ の正整数重み $c_{ij}$.

以下の $2$ 条件を満たすように辺部分集合 $S \subseteq E$ を選ぶとき,重みの和 $\sum_{(i, j) \in S} c_{ij}$ の最大値を求めよ:

  • 二部グラフ $(U, V; S)$ において各頂点 $i \in U$ の次数は $1$,
  • 二部グラフ $(U, V; S)$ において各頂点 $j \in V$ の次数は $b_j$.

入力に対する制約

  • $3 \leq N \leq 10^5$
  • $1 \leq c_{ij} \leq 10^9$

この形の問題は最大重み b-マッチング問題(最小費用流問題の特殊ケース)として知られており,多項式時間で解けます.しかし $N$ の上限が $10^5$ と大きいため,一般的なアルゴリズムでは $2$ 秒の制限に間に合いません.片側の頂点数が $|V| = 3$ ととても小さいことを利用した高速なアルゴリズムを設計する必要がありそうです.

解法

以下の 3 ステップで解きます:

  1. 問題を整数計画で定式化し,双対をとる.
  2. 変数の少ない問題に変形する.
  3. L 凸関数最小化アルゴリズムを適用する.

想定解法を $|V| \geq 4$ の場合に拡張するのは難しそうですが,ここで紹介する解法だと $|V| = 5$ 程度でも $2$ 秒以内に答えが計算できます.嬉しいですね!

整数計画による定式化とその双対

最大重み b-マッチング問題は整数計画として定式化できます:

\begin{align*} \max_{x \in \Z_+^{N \times 3}} \ & \sum_{i \in U,\ j \in V} c_{ij} x_{ij} \\ \text{s.t.}\ \ \ &\sum_{j \in V} x_{ij} = 1 \quad (i \in U),\\ & \sum_{i \in U} x_{ij} = b_j \quad (j \in V). \end{align*}

この最適化問題の変数は,$N \times 3$ の非負整数行列 $x = (x_{ij})$ です.非負整数と言っても,$1$ つ目の制約から $x_{ij}$ は $0$ か $1$ しかとれません.$x_{ij} = 1$ は辺 $(i, j)$ が $S$ に含まれることを表します.

最適化問題を見たらとりあえず双対をとりたくなりますよね.ここは知識パートですが,「上の問題を連続化し,双対をとって,離散化した問題」,つまり

\begin{equation}\label{prob-pq}\begin{aligned} \min_{p \in \Z^N,\ q \in \Z^3} & \sum_{i \in U} p_i - \sum_{j \in V} b_j q_j \\ \text{s.t.}\ \ \ \ & \ p_i - q_j \geq c_{ij} \quad (i \in U,\ j \in V)\end{aligned}\end{equation}

の最適値はもとの問題の最適値と等しいことが知られています.よってこの整数計画が解ければよいです.(この辺りの話は,室田一雄先生の講義資料「整数計画と LP の整数性」などにあります.)

変数の少ない問題に変形

問題 \eqref{prob-pq} が $N$ 次元ベクトル $p$ と $3$ 次元ベクトル $q$ についての最適化問題であることに注目し,よりサイズの小さいベクトルである $q$ を固定してみます.

問題 \eqref{prob-pq} の目的関数から $p_i$ の値はなるべく小さくしたいので,$q$ が固定された下で最適な $p$ は $p_i = \max_j \! \set{q_j + c_{ij}}$ と明示的に書けます.よって,問題 \eqref{prob-pq} は

\begin{equation}\label{prob-q} \min_{q \in \Z^3} \ f(q) := \sum_{i \in U} \! \max_{j \in V} \set{q_j + c_{ij}} - \sum_{j \in V} b_j q_j \end{equation}

という,たった $3$ 変数の問題と等価です! かなり解きやすくなりました.

L 凸関数最小化

問題 \eqref{prob-q} の目的関数 $f$ は L 凸関数というものになっているので,L 凸関数の最小化アルゴリズムが適用できます.これについては後の節で説明します.

解法全体の計算量は $O(N \max_{i,j} \log c_{ij} )$ となります(これも後で議論します).

実装例

L 凸関数最小化のために,後述の「スケーリング法」を使っています.実装は軽めで,

  • L 凸関数最小化ライブラリが 20 〜 30 行程度
  • 問題固有の部分は入出力を除くと 5 〜 10 行程度

です.

#include<bits/stdc++.h>
using namespace std;
using ll = long long;

// L 凸関数最小化ライブラリ
// O(2^n n^2 log (K/n)) * EVAL
template<class TD, class TR, TR (*f)(vector<TD>)>
struct Lconvex {
  int n;
  Lconvex(int _n) : n(_n) {}
  pair<TR, vector<TD>> min_scaling(vector<TD> x, TD alpha) const {
    TR f0, f1 = f(x);
    while (alpha) {
      f0 = f1, f1 = min_cube(x, f0, alpha);
      if (f0 == f1) alpha >>= 1;
    }
    return {f1, move(x)};
  }
private:
  TR min_cube(vector<TD> &x, TR f_best, const TD alpha) const {
    vector<TD> x0{x};
    for (int S = 1; S < (1 << n) - 1; S++) {
      vector<TD> xS{x0};
      for (int i = 0; i < n; i++) if (S >> i & 1) xS[i] += alpha;
      TR fS = f(xS);
      if (fS < f_best) f_best = fS, x = move(xS);
    }
    return f_best;
  }
};


int n;
array<int, 3> x;
vector<array<int, 3>> c;

ll f(vector<int> q) {
  ll res = 0;
  for (int i = 0; i < n; i++) res += max({c[i][0] - q[0], c[i][1] - q[1], c[i][2] - q[2]});
  for (int j = 0; j < 3; j++) res += ll(x[j]) * q[j];
  return res;
}

int main() {
  cin.tie(nullptr); ios::sync_with_stdio(false);
  for (auto &e : x) cin >> e;
  n = x[0] + x[1] + x[2];
  c.assign(n, array<int, 3>());
  for (auto &v : c) for (auto &e : v) cin >> e;

  Lconvex<int, ll, f> lc(3);
  auto [ans, q] = lc.min_scaling({0, 0, 0}, 1<<29);
  cout << ans << '\n';
}

ACコード(C++, 1407 Byte, 74 ms):
https://atcoder.jp/contests/agc018/submissions/23527593

L 凸関数

この節では AGC018: C - Coins の解法で現れた「L 凸関数」について説明し,次の節で L 凸関数最小化アルゴリズムについて説明します.

L 凸関数とは

L 凸関数の定義として同値なものがいくつか知られています.そのうちの一つを紹介します.以下では $\bar \R := \R \cup \set{\infty}$ とします.

定義: L 凸関数

関数 $f: \Z^n \to \bar \R$ が条件

  1. 任意の $x, y \in \Z^n$ に対し $\displaystyle f \Big( \Big\lfloor \frac{x+y}{2} \Big\rfloor \Big) + f \Big( \Big\lceil \frac{x+y}{2} \Big\rceil \Big) \leq f(x) + f(y)$
  2. ある実数 $r$ が存在し,任意の $x \in \Z^n$ に対し $f(x + \1) = f(x) + r$

をともに満たすとき,$f$ を L 凸関数という.

ここで,ベクトル $x$ に対する $\lfloor x \rfloor,$ $\lceil x \rceil$ は,要素ごとの床関数,天井関数を表します.また,太字の $\1$ は全成分が $1$ のベクトルです.

$f(x) = \infty$ となることも許していることに注意してください.ここでは $\infty + 3 = \infty$ などと形式的に計算することとし,$\infty = \infty$ や $\infty \leq \infty$ は正しい式だと考えます.

L 凸関数の定義には条件が 2 つありますが,条件 1 がより本質的です.連続変数の連続関数 $f: \R^n \to \bar \R$ が凸であることの必要十分条件は,\[ 2 f \Big( \frac{x+y}{2} \Big) \leq f(x) + f(y) \]が任意の $x,y \in \R^n$ で成り立つことです.条件 1 は,これを $\Z^n$ の世界に自然に翻訳したものと考えられます.このため,L 凸関数は名前の通り「凸っぽい」離散の関数と言えます.

条件 1 を満たすが条件 2 を満たすとは限らない関数 $f: \Z^n \to \bar \R$ は $\mathrm{L}^\natural$ 凸関数と呼ばれます(エル・ナチュラルと読みます).

L 凸関数の典型例

何らかの関数が得られたとき,それが L 凸かどうかを定義通りチェックするのは少々手間です.以下に挙げるような L 凸関数の典型例を知っておくと,「これは L 凸っぽいな」と分かるようになり便利です.

一次関数

一次関数\[f(x) = \sum_{i=1}^n a_i x_i + b\]は L 凸です.

差の凸関数

連続変数の凸関数 $g_{ij}: \R \to \bar \R$ $(1 \leq i,j \leq n)$ に対し\[ f(x) = \sum_{i,j} g_{ij}(x_i - x_j) \] は L 凸です.

L 凸関数の和

$2$ つの関数 $f_1, f_2: \Z^n \to \bar \R$ が L 凸ならば,その和 \[f(x) = f_1(x) + f_2(x)\] も L 凸です.これを繰り返し適用すれば,$3$ つ以上の L 凸関数の和も L 凸関数だと分かります.

L 凸関数の射影

$g: \Z^n \to \bar \R$ を L 凸関数とします.変数を $x = (x_1, x_2) \in \Z^{n_1} \times \Z^{n_2}$ と $2$ つのブロックに分割し,$x_2$ について最小化することで得られる $n_1$ 変数関数\[ f(x_1) = \min_{x_2 \in \Z^{n_2}} g(x_1, x_2)\] も L 凸です.

Coins の関数の L 凸性

AGC018: C - Coins の解法で登場した,式 \eqref{prob-q} の関数 $f$ が L 凸であることを確認しましょう.上で紹介した典型例を使って確認します.

$i \in U = \set{1,\dots,N},$ $j \in V = \set{1,2,3}$ に対し連続変数の関数 $h_{ij}: \R \to \bar \R$ を \[ h_{ij}(x) := \begin{cases} 0 & \text{if $x \geq c_{ij}$},\\ \infty & \text{otherwise} \end{cases} \]と定めると,これは凸関数です.よって $(N+3)$ 変数関数 \[ h(p, q) := \sum_{i \in U,\ j \in V} h_{ij}(p_i - q_i) \] は L 凸です(差の凸関数).

一次関数が L 凸であることと,L 凸関数の和が L 凸であることから,\[g(p, q) := h(p, q) + \sum_{i \in U} p_i - \sum_{j \in V} b_j q_j\]も L 凸です.この関数 $g$ は問題 \eqref{prob-pq} の制約を目的関数に組み込んだものになっています.この $h_{ij}$ や $h$ のように,条件を満たすなら $0$,満たさないなら $\infty$ をとる関数を用いて制約式を目的関数に組み込む方法は,数理最適化でよく用いられます.

式 \eqref{prob-q} の関数 $f$ は \[ f(q) = \min_{p \in \Z^N} g(p, q) \] として得られたものでした.これは L 凸関数の射影なので,確かに $f$ が L 凸関数だと分かりました!

上と同様にして,一般の有向グラフ上の最短路問題や最小費用流問題の双対問題(を離散化したもの)が,L 凸関数最小化問題の例になっていることが分かります.

L 凸関数最小化アルゴリズム

L 凸関数 $f: \Z^n \to \bar \R$ の最小化問題を考えます.

L 凸関数の定義の条件 2 より,ある $r$ に対し $f(x + \1) = f(x) + r$ が成り立ちますが,仮に $r \neq 0$ とすると $f$ は下に非有界になってしまう(最小値が存在しない)ので,以下では $r = 0$ を仮定します.もちろん式 \eqref{prob-q} の関数 $f$ でも $f(q + \1) = f(q)$ が成り立っています.

スケーリング法

L 凸関数 $f$ の最小化アルゴリズムはいくつか知られていますが,その中でも競技プログラミングと相性が良さそうなのがスケーリング法です.

L 凸関数最小化のためのスケーリング法
  1. $f(x) < \infty$ なる初期点 $x \in \Z^n$ と,$2$ 冪の整数 $\alpha$ を定める.
  2. $\alpha \geq 1$ である限り以下を繰り返す:
    1. $x + \set{0, \alpha}^n$ において $f$ を最小化する点を $y$ とする1
    2. $x \neq y$ ならば $x \gets y$ と更新する.$x = y$ ならば $\alpha \gets \alpha/2$ と更新する.
  3. 解 $x$ を出力する.

$x + \set{0, \alpha}^n$ は,点 $y \in \Z^n$ のうち,各 $i$ に対し $y_i \in \set{x_i, x_i + \alpha}$ であるようなものからなる集合を表します.これは $n$ 次元超立方体の頂点集合で,$2^n$ 個の点からなります.

アルゴリズム中の $\alpha$ は「歩幅(ステップサイズ)」のようなものです.つまり,最初は「大股」で解を探索し,大股ではより良い解が見つけられなくなったら歩幅を半分にする,ということを繰り返しています.直観的にも分かりやすいアルゴリズムですね.

このアルゴリズムでは,反復が進むごとに点 $x$ の各成分の値が広義単調に増加します.これで最適解が見つかるのは不思議に思うかもしれませんが,ここに $f(x + \1) = f(x)$ の仮定が効いています.つまり,$x^* \in \Z^n$ が最適解ならば,任意の $k \in \Z$ に対して $x^* + k \1$ も最適解なので,各成分の値を減らさずとも最適解が見つかるというわけです.

実は $\mathrm{L}^\natural$ 凸関数最小化も似たアルゴリズムで解くことができます.$\mathrm{L}^\natural$ の場合は $f(x + \1) = f(x)$ が使えないため,$x$ の成分の値が減少する方向にも探索する必要があります.

初期ステップサイズの選び方と計算量

スケーリング法の計算量($f$ の関数値評価の回数)を抑えるには,最初の歩幅 $\alpha$ の選び方が重要です.それに関する定理が知られています.

$f$ の実質的な定義域(実効定義域と呼ばれます)$\dom f$ を $\dom f := \set{ x \in \Z^n \mid f(x) < \infty }$ と定義し,$\dom f$ の大きさの指標を

\[ K := \max \set{ \|x - y\|_\infty \mid x, y \in \dom f, \ x_i = y_i \ (\exists i) } \]

と定めます.$\|x\|_\infty := \max_{i} |x_i|$ です.

定理: スケーリング法の反復回数

$K < \infty$ とする.スケーリング法において,初期ステップサイズを $\alpha = \Theta(K/n)$ となるように定めると,ステップ 2-1 の実行回数は $O(n^2 \log (K/n))$ となる.

スケーリング法において最も計算が重いのがステップ 2-1 です.例えば,ステップ 2-1 で $2^n$ 個の点を総当りすると,アルゴリズム全体で $f$ の関数値評価を $O(2^n n^2 \log (K/n))$ 回呼び出すことになります.詳細は割愛しますが,ステップ 2-1 で劣モジュラ関数最小化アルゴリズムを使い,$f$ の評価回数を $O(\mathrm{poly}(n) \log (K/n))$ に抑えることもできます.

この記事で紹介した実装例では,スケーリング法のステップ 2-1 で $2^n$ 個の点を総当りしています.

Coins の解法の計算量

AGC018: C - Coins の解法で現れた式 \eqref{prob-q} の関数 $f$ に対しては,$K = \infty$ となってしまっています.このままでは上の定理が使えませんが,この場合でも $f$ の実効定義域 $\dom f$ を適切に制限することで,$K < \infty$ にできることを確認します.

式 \eqref{prob-q} からずいぶん離れてしまったので再掲します:

\[ \min_{q \in \Z^3} \ f(q) := \sum_{i \in U} \max_{j \in V} \! \set{q_j + c_{ij}} - \sum_{j \in V} b_j q_j. \]

$U = \set{1,\dots,N},$ $V = \set{1,2,3}$ であることも再度書いておきます.

少々唐突ですが,ある $j, k \in V$ が存在し,任意の $i \in U$ に対し $q_j + c_{ij} > q_k + c_{ik}$ となるような $q$ は,最適解にはなりえません.そのような $q$ よりも $q_j$ を $1$ だけ小さくした解の方が,目的関数値が小さくなるからです.

よって,$q_j - q_k \leq \max_i \! \set{c_{ik} - c_{ij}}$ $(j, k \in V)$ を満たさない $q$ に対しては $f(q) = \infty$ であるとしても最小値は変わりません.このように $f$ を修正しても $f$ の L 凸性は保存されます2

こうすることで,$K < \infty$ になります.具体的には $K \leq 2 \max_{i,j} c_{ij}$ となることが示せます.よって,この問題の制約の下では $K = 2 \times 10^9$ と思ってスケーリング法を適用すればよいです.

$f$ の関数値評価の計算量は $O(N)$ なので,解法全体の計算量は $O(N \max_{i,j} \log c_{ij})$ となります.より正確には,$M := |V|$ として $O(2^M M^2 N \max_{i,j} \log c_{ij})$ です.

参考図書

この記事を書くにあたって参照した離散凸解析の教科書 4 冊を紹介します.

離散凸解析の考えかた

  • 離散凸解析が概観できる.
  • 図や例(反例)が豊富で読みやすい.最初の一冊にオススメ!
  • 第 13, 14 章の例が特に面白い.
  • 証明はあまり載っていない.
離散凸解析の考えかた 最適化における離散と連続の数理 | 室田 一雄 |本 | 通販 | Amazon
Amazonで室田 一雄の離散凸解析の考えかた 最適化における離散と連続の数理。アマゾンならポイント還元本が多数。室田 一雄作品ほか、お急ぎ便対象商品は当日お届けも可能。また離散凸解析の考えかた 最適化における離散と連続の数理もアマゾン配送商品なら通常配送無料。

離散凸解析と最適化アルゴリズム

  • 組合せ最適化の有名な問題とそのアルゴリズムが,離散凸解析の視点から学べる.
  • (離散凸解析というより)組合せ最適化の教科書としても使える.もちろん離散凸解析も学べる.
  • アルゴリズムがたくさん載っている.競技プログラマにオススメ!
離散凸解析と最適化アルゴリズム (数理工学ライブラリー) | 一雄, 室田, 昭義, 塩浦 |本 | 通販 | Amazon
Amazonで一雄, 室田, 昭義, 塩浦の離散凸解析と最適化アルゴリズム (数理工学ライブラリー)。アマゾンならポイント還元本が多数。一雄, 室田, 昭義, 塩浦作品ほか、お急ぎ便対象商品は当日お届けも可能。また離散凸解析と最適化アルゴリズム (数理工学ライブラリー)もアマゾン配送商品なら通常配送無料。

離散凸解析

  • 最も詳しい和書.
  • ほぼ全ての命題,定理に証明がついている.
  • 離散凸解析を腰を据えてじっくり学びたい人向け.
離散凸解析 (共立叢書 現代数学の潮流) | 室田 一雄 |本 | 通販 | Amazon
Amazonで室田 一雄の離散凸解析 (共立叢書 現代数学の潮流)。アマゾンならポイント還元本が多数。室田 一雄作品ほか、お急ぎ便対象商品は当日お届けも可能。また離散凸解析 (共立叢書 現代数学の潮流)もアマゾン配送商品なら通常配送無料。

Discrete Convex Analysis

  • 400 ページ近くある洋書.最も詳しい本.
  • ほぼ全ての命題,定理に証明がついている.
  • 離散凸解析に関連する研究をする人向け.
Amazon | Discrete Convex Analysis (Monographs on Discrete Mathematics and Applications, Series Number 10) | Murota, Kazuo | Calculus
Amazon配送商品ならDiscrete Convex Analysis (Monographs on Discrete Mathematics and Applications, Series Number 10)が通常配送無料。更にAmazonならポイント還元本が多数。Murota, Kazuo作品ほか、お急ぎ便対...
  1. $f$ を最小化する点が複数ある場合,そのうち極小なものをとります.ここでは,集合 $S \subseteq \Z^n$ の元 $a \in S$ に対し,$b \leq a$(成分ごとの不等式)となる $b \in S$ が $a$ 以外に存在しないとき,$a$ を $S$ の極小な点と呼んでいます.$x + \set{0, \alpha}^n$ において $f$ を最小化する点の集合を $S$ とすると,$S$ の極小な点がただ一つ存在することが知られています.
  2. 「差の凸関数」が L 凸であることと,L 凸関数の和が L 凸であることから分かります.
タイトルとURLをコピーしました