ざっくり解説: 焼きなまし法はなぜうまくいくのか？

概要

焼きなまし法は遺伝的アルゴリズムなどと並ぶ最適化問題の発見的解法の一つとして有名であり, 競技プログラミングの文脈ではマラソン(厳密解を求めるのが困難とされる一つの最適化問題に長時間取り組み最も最適値に近い解を得るという種目) においては常套手段の一つとして用いられています. 従ってこの文脈では様々な解説記事があるのですが, いわゆるマラソン競技特有の「コツ」や実装方法の説明がほとんどで, 焼きなまし法自体を数理的に深掘りした記事もなくはないのですが, やはり「なぜうまくいくのか」の観点で説明している記事はほとんど見たことがないので, なんとなくこの辺りが気になるという人向けにこの記事で解説していきます (公開時は把握してなかったのですが, 本記事の内容はこちらの記事とほぼ同じ内容になっていますので, 数理的な説明をする記事が全くないというわけではありません). 内容としてはメトロポリス・ヘイスティング法の説明になりますが, 数学的な厳密性は犠牲にして自分の理解で執筆していますので, 厳密な議論まで気になる方は論文を読んでください (もしかしたら他の説明の仕方があるかもしれません).

結論から述べると

• 焼きなまし法は, 遷移確率が時間とともに変動するランダムウォークを実用的な観点から修正したものとみなすことができる.
• このランダムウォークは, 遷移確率の変動が各時刻において十分に微小であれば, そのランダムウォークの収束先は最適解上の分布になることが保証される.
• しかし, 理論的に保証される収束は非常に遅いため, 実用の場面では, 適当なところでランダムウォークを打ち切ったり, より良い解を得るための(理論保証のない)修正を施す. これがマラソンなどの場面でよく見かける焼きなまし法である.

この記事にアクセスする人は焼きなまし法についてはなんとなく知ってるということを仮定します. ただ, 念の為, 初見の人にもアクセスできるよう, 焼きなまし法に関する説明を次の章で与えておきます. 焼きなまし法を知っている人はその次の2章までスキップして構いません.

1. 焼きなまし法とは?

焼きなまし法とは最適化問題に対し, 最適に近いであろう実行可能解を得る手法の一つで, 局所探索と呼ばれるカテゴリーに属します. 局所探索とは, 何かしらの実行可能解を一つ用意し, その解を逐次更新していくというアルゴリズムの総称であり, 焼きなまし法以外には山登り法などが知られています. ここでは以下の記法を用います.

• 実行可能解全体の空間は有限集合とし記号$S$で表します. 
• 実行可能解を$x,y,z\in S$などと表します
• 目的関数値を$f\colon S \to \mathbb{R}_{\ge 0}$で表し, 最小化を目指すこととします.
• 実行可能解$x\in S$の近傍集合を$N(x) \subseteq S$とし, その要素$y\in N(x)$を近傍と呼ぶ.

すなわち, $S$という有限集合と$S$上に非負値をとる関数$f$があり, $f(x)$が最小となる$x\in S$を見つけることが目的となります. 直感的には近傍$y \in N(x)$は$x$と「ちょっとだけ違う実行可能解」です. 実行可能解と目的関数の定義は考える最適化問題に依存しますが, 実行可能解はアルゴリズムの設計者が定義するため, 同じ最適化問題であっても様々な近傍を考えることができます.

山登り法では, 現在の実行可能解$x\in S$に対し, $f(y)$が最小となる近傍$y\in N(x)$を選び, それを$x$に置き換えて同じ処理を繰り返すというアルゴリズムです (または$N(x)$の元$y$を順番に一つずつ見ていき, $f(x)>f(y)$を満たす近傍$y$が見つかったらその瞬間に$x:=y$とするようなものも考えられる. これは, $|N(x)|$の要素数が非常に大きい時は$\min_{y\in N(x)} f(y)$を計算に時間がかかる際に有効な方法である).

このアルゴリズムだと谷底 (全ての近傍の目的関数値が今の解より真に良くならないような解; 局所最適解と呼ぶ) にいきつくとそれ以上動かなくなるのでアルゴリズムは終了します. 上図の例では偶然にも最適解に辿り着きますが, 一般に谷底は複数あり, どの谷底に至るかは一番最初に選んだ解に依存するため, 最適性の保証はありません (ただし空間が連続的で目的関数値が凸性を持つなど, 特定の状況下では最適性が保証できます).

一方, 焼きなまし法では, 現在の実行可能解$x\in S$に対し, 近傍$y\in N(x)$を順番に見ていき,

1. 目的関数値が改善, すなわち$f(x)>f(y)$を満たすとき, $x:=y$として次の反復へ移行
2. そうでない場合, $f(x),f(y)\in\mathbb{R}_{\ge 0}$, そしてこれまでの反復回数$t$のみから定まるある確率$p$で$x:=y$として次の反復へ移行

というアルゴリズムです. 例えば$p=0$ならば山登り法に一致します.

焼きなまし法では確率的に登ることがある.

ここで確率$p$は

$$ p = \exp\left(-\frac{f(y) - f(x)}{T}\right) \tag{1} $$

とします. ここで$T=T(t)>0$は温度と呼ばれるパラメータであり, これまでの反復回数$t$の関数であり, $t\to\infty$としたとき$T\to 0$となります. ここで, $f(y)\ge f(x)$を仮定しているため常に$p\in[0,1]$です.

確率$p$については以下の観察ができます:

• 解の改悪幅$f(y)-f(x)$が大きいほど$p$が小さい, すなわち$x\to y$の遷移が発生しにくい.
• 改悪幅を正の値で固定したとき, 温度$T$が$0$に近いほど$p$は小さくなります. つまり, $T$が大きいときは改悪($f(x)<f(y)$なる$y\in N(x)$への更新)が発生しやすくなるが, 反復を繰り返し$t\to\infty$としていくと$T\to 0$に近づくため, より山登り法に近い振る舞いをするようになる.

温度$T(t)$の定め方は「非常にゆっくり$T\to 0$に収束する関数」であることが理想ですが, 実際に焼きなまし法を走らせる際は手っ取り早く良い解を得たい場合がほとんどなので, あらかじめ反復回数を決定しておき, 最終的に$T\approx 0$になるように線形で下がっていくようなものを考えることもあります. ただし理論的な保証があるのは「非常にゆっくり」収束するという条件が必要なので, この時点で理論保証から外れてしまうことになります.

2. メトロポリス-ヘイスティング法

突然ですが, 有限集合$S$上の関数

$$ f\colon S \to \mathbb{R}_{\ge 0} $$

が定まっているとき, $f(x)$を最小にする$x\in S$を求めるにはどうすれば良いでしょうか? 唐突ですが, 温度と呼ばれるパラメータ$T>0$に対して$S$上の確率分布$\pi\in[0,1]^S$を以下で定めます:

$$ \pi (x) = \frac{\exp(- f(x)/T)}{Z}. \tag{2}$$

ここで$Z = \sum_{x\in S}\exp(-f(x)/T)$は$\pi$を分布に正規化する値 (様々な文脈で分配関数と呼ばれる値) です. また, この分布$\pi$を(エネルギー$f$に対する)ボルツマン分布と呼ぶらしいです.

分布$\pi$に従って生成される確率変数$X \sim \pi$を考えると, 以下の観察が芽生えます:

• 温度$T$が非常に大きい値, 極端に言えば$T=\infty$ のとき, $f(x)$の影響は無視されるため, $\pi$は$S$上で一様ランダムとなります.
• 温度$T$が非常に小さい値, すなわち$T\approx 0$ のとき, $\pi$は最適解$S^*=\{x\in S\colon f(x) = \min_{y\in S}f(y)\}$上の一様分布となります.

従って, 理想的には$T\approx 0$において$\pi$に従う確率変数を生成できれば最適化問題を解いたことになります!! ところが, 分配関数$Z$が簡単に計算できそうな見た目をしてないので, $\pi$に従う確率変数をサンプリングするということ自体が非自明なタスクになっています. メトロポリス-ヘイスティング法とは, 目標となる分布$\mu$があったとき, ランダムウォークを用いてこの分布に非常に近い分布に従う確率変数を生成する手法の一つです.

具体的には, $S$上の何らかのランダムウォークの遷移確率行列$P \in [0,1]^{S\times S}$が与えられているとします. ここで$P(x,y)$は$x$から$y$に遷移する確率を表しており, 例えば焼きなまし法の例でいうと, 一様ランダムな近傍に遷移する, というランダムウォークだと思えばOKです. 一般にランダムウォークはその遷移がとある条件 (例えば遷移の強連結性など) を満たせば, 定常分布と呼ばれるある分布に収束するという事実が知られています. 今考えているランダムウォーク$P$もその条件を満たしているとします. ただし$P$の定常分布は目標としている分布$\pi$とは異なっていることもあります.

メトロポリス-ヘイスティングス法のイメージ

メトロポリス-ヘイスティング法とは, 手元にあるランダムウォーク$P$に適当な重さの自己ループの追加などにより修正し, 目標分布$\pi$を定常分布とするような別のランダムウォーク$P'$をうまく定めるという手法です.

より詳細を述べると, 修正後のランダムウォーク$P'$は以下のような形をしています: ランダムウォーク$P$により, $x\to y$という遷移が発生したとき, 確率$a(x,y)$でその遷移を受理し実際に遷移し, 確率$1-a(x,y)$でその遷移を拒否し$x$にとどまる.

簡単のため$P(x,x)=0$とすると, $P'\in[0,1]^{S\times S}$の成分は

$$
P'(x,y) = \begin{cases}
P(x,y)\cdot a(x,y) & & \text{if }x\ne y,\\
1-\sum_{y} P(x,y)\cdot a(x,y) & & \text{otherwise}
\end{cases}
$$

となります. 肝心の受理確率$a(x,y)$の定め方ですが, $P'$が定常分布$\pi$を持つように, そして$a(x,y)$の値ができるだけ大きくなるように (自己ループの確率が少ない方が遷移する回数が多くなるため) 設計します. 詳細は省きますが釣り合い条件と呼ばれる条件から, $P'$は次の等式

$$ {}^{\forall} x,y\in S,\quad\pi(x) P'(x,y) = \pi(y) P'(y,x) $$

を満たしていなければなりません. これに$P'$の成分を代入すると

$$ \pi(x) a(x,y) P(x,y) = \pi(y) a(y,x) P(y,x) $$

となります. 

計算の簡単のため, $P$が対称行列だとしましょう. すると上の等式は$\pi(x)a(x,y) = \pi(y)a(y,x)$となります. ここで, できるだけ$a(\cdot,\cdot)$の値を大きくとると, 片方は$1$に設定できます. 従って

\begin{align*}
&a(x,y) = \min\left\{ \frac{\pi(y)}{\pi(x)},\,1\right\},\\
&a(y,x) = \min\left\{ \frac{\pi(x)}{\pi(y)},\,1\right\}
\end{align*}

を得ます. このようにして定まる受理確率$a(x,y)$は$\pi(x)$と$\pi(y)$の比のみに依存するという素晴らしい性質を持っています. 実際, $\pi$が(2)のような形になっているとき, その比は分配関数$Z$に依存せず,

$$ \frac{\pi(y)}{\pi(x)} = \exp\left( - \frac{f(y)-f(x)}{T} \right) $$

で簡単に計算できます. この値はまさに(1)で考えていた値と一致します. つまり, 焼きなまし方とは, メトロポリス-ヘイスティングス法の受理確率を模した確率的な遷移を許容する局所探索といえます.

なお, 実際の焼きなまし法では温度パラメータ$T$は遷移回数とともに変動します. 直感的にはこの変動は十分ゆっくりであれば局所的には$T$が固定された値とみなせるので解析できそうですが, この直感を厳密な形にするにはそれなりの苦労が必要のようです (例えばこの記事. この辺は私もまだわかっていない)

なお, メトロポリス-ヘイスティングス法で得られるランダムウォーク$P'$は, 収束先の定常分布が$\pi$であるということのみが保証されるため, その収束の速さについては何も保証がありません. しかし, 適当な仮定の下では何かしらの収束時間のバウンドを与えることは可能となります.

3. まとめ

冒頭にも述べましたが, 要するに

1. 最適化問題を解くというタスクをある目標分布に従うサンプリングに置き換える
2. このサンプリングタスクをメトロポリス-ヘイスティング法で解く
3. プロセス1,2の構造をできるだけ保ちつつ, 実用的に手っ取り早く良い解を得るための工夫を施して実装したものが焼きなまし法である

ということになります (もしもうちょっと深い理論の話があったら教えていただけると幸いです). 

なお, 最適化問題は$f$が凸性を持つときは「局所最適解は最適解の一つである」が成り立つため非常に解きやすくなるのですが, 同様に分布$\pi$は対数凹性を持つときはサンプリングしやすいという事実があるようです. 実際, $f$が凸のとき, 対応するボルツマン分布$\pi$は$\log \pi(x) = -\frac{f(x)}{T}$より, 対数凹性を持ちます.

緩和の順番を工夫してBellman-Ford法は改善できる

久々にブログのタイトル通りにアルゴリズム(Bellman-Ford法)に関する記事を書きます. それなりに丁寧に書いたつもりなので, Bellman-Ford法を知らない人でも楽しめると思います. プログラムの実行時間はアルゴリズムに加えて実装や実行環境に依存するため, アルゴリズムの計算量では原則として定数倍を無視して議論するわけですが, この記事で紹介する手法はそのような実行環境に依存せず, 本質的に演算回数を定数倍減らしているので競プロ的な(ナイーブな手法を改善していくという)面白さがあると思います.

概要. Bellman-Ford法は$O(mn)$時間で単一始点最短経路問題を解きます. このアルゴリズムは「各辺に対して緩和と呼ばれる操作を行う」ことを$n-1$回繰り返すというアルゴリズムです. 実は緩和を行う辺の順番を工夫することによって反復回数を$n-1$から$n/2$に改善することが可能です [Yen, 1970]. 本稿ではこの単純かつ素晴らしいアイデアを紹介します.

1. Bellman-Ford法とその正当性の証明

フォーマルに書くとBellman-Ford法は以下で定まるアルゴリズムです. まず距離を格納するリスト$\mathrm{dist}[u]$を,

\begin{align*}
\mathrm{dist}[u] = \begin{cases}
0 & \text{if $s=u$},\\
\infty & \text{otherwise}.
\end{cases}
\end{align*}
で初期化し, その後, 「各辺$e\in E$に対して $\mathrm{relax}(e)$ を行う」を$n-1$回繰り返します. ここで$\mathrm{relax}(e)$とは, 有向辺$e=(u,v)$に対し, その重みを$w_e$とすると,
\[
\mathrm{dist}[v] \leftarrow \min\{\mathrm{dist}[v], \mathrm{dist}[u]+w_e\}
\]
で更新する操作である.

Bellman-Ford法では見る辺の順番に関しては特に言及されておらず, どのような順番で辺を見たとしても高々$n-1$回の反復で(負閉路がない場合の)単一始点最短経路が正しく解けることが保証されます. その正当性として以下の「よくある」証明が知られています.

「よくある」Bellman-Ford法の証明.

各$i=0,\dots,n-1$に対して, $i$回目の反復が終了した時点での配列$\mathrm{dist}$を$\mathrm{dist}_i$と書くことにします. ただし$\mathrm{dist}_0$は初期化によって得られる$\mathrm{dist}$で, $\mathrm{dist}_{n-1}$はBellman-Ford法の出力です. 次の主張を$i$に関する帰納法で示します: 各$u\in V$に対して$\mathrm{dist}_i[u] = su$-経路のうち辺数が高々$i$であるもののうち最短の長さである. 実際, $i=0$の場合はこの主張は正しいです. また, $i\leq j$まで主張が正しいと仮定すると, $j+1$回目の反復について考えると, 漸化式
\begin{align*}
\mathrm{dist}_{j+1}[u] = \min_{v \in N^-(u) \cup \{u\}}\{ \mathrm{dist}_j[v] + w_{vu} \}
\end{align*}
を得ます (ただし$N^-(u)=\{v\in V\colon (v,u)\in E\}$は頂点$u$に入っていく辺のもう一方の端点からなる集合で, $w_{uu}=0$と定義している). 帰納法の仮定より各$\mathrm{dist}_j[v]$は辺数$j$本以下の中での最短$sv$経路の長さなので, $\mathrm{dist}_{j+1}[u]$は所望のものになっており, 主張が証明されました. さて, 元のグラフに負閉路がない場合は$s$から辿り着ける全ての頂点$u$に対して最短$su$経路は辺数が高々$n-1$なので, アルゴリズムの出力$\mathrm{dist}_{n-1}[u]$は確かに頂点$s$を始点とした最短経路長を格納した配列になっています. (証明終)

上記の帰納法に基づく証明により, Bellman-Ford法は始点$s$から近い順に最短経路長を決めていることがわかります. イメージとしては, $s$を根とした最短経路木(最短経路に属する辺からなる木. 最短経路が複数あるときは任意に一つだけ選んでその辺を木に追加していって構成できる) に対して, 根から順番に(幅優先探索の要領で)最短経路長が決定していくことがわかります. 特に, 各反復では最短経路木に含まれる辺のうち, 少なくとも1本が最短経路として同定されていきます. 木の辺数は$n-1$であることからもBellman-Ford法の反復回数が高々$n-1$であることが納得できると思います.

2. Yen(1970)の方法

Bellman-Ford法の解析では 「各反復では最短経路木に含まれる辺のうち, 少なくとも1本が最短経路として同定されていく」ことを示したわけですが, それでは「各反復では最短経路木に含まれる辺のうち, 少なくとも2本が最短経路として同定されていく」ならば, 反復回数は高々$\lceil n/2\rceil$となります. Yen(1970)による定数倍改善はこの観察に基づきます. Yen(1970)は, 以下のアルゴリズムによって定まる反復に基づくBellman-Ford法の反復回数は高々$\lceil n/2 \rceil$回であることを主張します:

1. 頂点に適当に番号づけて, $V=\{v_1,\dots,v_n\}$, 始点を$s=v_1$とする. 各辺$e=(v_i,v_j)$に対し, $i<j$ならば$e$を上昇辺, $i>j$ならば$e$を下降辺と呼ぶ.

2. 頂点を$v_1,\dots,v_n$の順番に見ていき, 各$v_i$に対して$v_i$を始点とする全ての上昇辺に対して緩和を実行する.

3. 頂点を$v_n,\dots,v_1$の順番に見ていき, 各$v_i$に対して$v_i$を始点とする全ての下降辺に対して緩和を実行する.

以下にこのアルゴリズムの正当性の証明を載せます.

Yen(1970)の手法の正当性の証明

頂点$u\in V\setminus \{s\}$を固定し, $su$-最短経路を任意に一つ選び, その経路を頂点列とみなし$P=(w_0,w_1,\dots,w_\ell)$とします (ただし$s=w_0,u=w_\ell$). このパス$P$は上昇辺からなる「増加部分」と下降辺からなる「減少部分」の交互列に分解できます. 例えば$P$が

$v_1 \to v_3 \to v_5 \to v_4 \to v_2 \to v_7 \to v_6$

という頂点列の場合, 

• $v_1 \to v_3 \to v_5$
• $v_5 \to v_4 \to v_2$
• $v_2 \to v_7$
• $v_7 \to v_6$
  

のように, 番号に関する極大な増加列と極大な減少列が交互に現れるように分解できます (始点を$s=v_1$にしているので最初は必ず増加列です). この分解を$P=(U_1,D_1,U_2,D_2)$と表すことにします. ここで$U_1=(v_1,v_3,v_5), U_2 = (v_2,v_7)$は増加列, $D_1=(v_5,v_4,v_2), D_2 = (v_7,v_6)$は減少列です. 一般に増加部分を$U_1,\dots,U_p$, 減少部分を$D_1,\dots,D_q$とし, 連続する上昇列と下降列を繋げて得られる$(U_i,D_i)$という$P$の部分列を山折りと呼ぶことにします. 各山降りは少なくとも2本の辺を含むため, 最短経路$P$が含む山折りの個数は高々$\lceil n/2\rceil$です. Yen(1970)のステップ2では小さい頂点から順番に上昇列を緩和していっているため, $P$の上昇列を順番に同定していきます. 同様にステップ3では$P$の下降列を順番に同定していきます. ステップ2-3を交互に繰り返すため, 一回の反復では一つの山折りに対応する部分列$(U_i,D_i)$の最短経路が順番で同定されていきます. 従って, 山折りの個数, すなわち高々$\lceil n/2\rceil$回の反復で最短経路が求まります. (証明終)

3. その後の進展

[Bannister and Eppstein, 2012]によって, 乱択を用いた反復回数のさらなる定数倍改善$n/3$が得られています (著者のブログ記事). このアルゴリズムも非常に単純なテクニックに基づいたものになっています. 具体的には, Yen(1970)では任意に選ばれた頂点順番に基づいて上昇辺と下降辺が定義されていましたが, これをランダムな順番にすると期待値の意味で反復回数が$n/3$になる, ということを示しています. ここではラフな証明を載せておきます. Yen(1970)の上記の解析により, 任意の最短経路$P$を同定するのに必要な反復回数はそれに含まれる山折りの個数で抑えることができます. 従って, 頂点番号がランダムに振られたときの$P$に含まれる山折りの個数の期待値について考察します. 最短経路の頂点列を$(v_0,\dots,v_\ell)$とし, 各$v_i$はランダムに振られた番号を表すとすると, 山折りの個数の期待値はおおよそ$\ell \cdot \sum_{i=1}^n (1/n) \cdot (i/n)^2 \approx \ell \cdot \int_0^1 x^2 \mathrm{d}x \leq n/3$となります (この式における$(1/n)\cdot (i/n)^2$というのは, パス中の頂点を一つ固定したときにその頂点の番号が$i$になり, かつ両隣の頂点の番号がどちらも$i$以下となる確率を意味します). ざっと論文を目を通したくらいなので詳しいことはわかりませんが, おそらく山折りの個数の集中性を示して(集中不等式の文脈ですでに知られている結果がほぼブラックボックスに使えるっぽい), 終着頂点に関するunion boundをとることによって証明がなされています.

なお, 執筆時(2023年12月時点)では負辺ありの単一始点最短経路問題に対するアルゴリズムのオーダーの改善も知られており, 特にFOCSのbest paperをとったブレイクスルーの結果[Bernstein, Nanongkai, Wulff-Nilsen, 2022]によって, 負辺を含むグラフの単一始点最短経路問題が$O(m\log^8(n)\log W)$時間で解けることが示されました. ここで$W>0$は, $w_e\geq -x$を満たす最小の$x$に対して$W=\max\{2,x\}$で定義される値です. ($2$とのmaxをとっているのはlogをとったときに$0$にならないようにするため.) さらに翌年, [Bringman, Cassis, Fischer, 2023]によって$O(m\log^2(n)\log(nW)\log\log n)$時間アルゴリズムが発表されました.

少し余談になりますが, 同じ会議(FOCS2023)では最短経路に関する別の論文[Duan, Mao, Shu, Yin, FOCS23]もあって, この論文ではcomparison-addition modelと呼ばれるモデルにおいて負辺なしの最短経路問題をダイクストラ法より高速に$O(m\cdot \sqrt{\log n\cdot \log\log n})$時間で解く乱択アルゴリズムを提案しています. comparison-addition modelとは, Bellman-Ford法の緩和操作のように, 数値の加算と比較のみが定数時間で行えるような環境を考えるモデルです. 彼らのアルゴリズムのアイデアは根幹の部分ではダイクストラ法と似たようなことをやっています. 通常のダイクストラ法では全ての頂点をヒープに入れて近い順に頂点を取り出してその周りの辺を緩和するというものでしたが, Duanらのアルゴリズムではパラメータ$k=\sqrt{\log n / \log\log n}$に対して, $n/k$個の頂点をランダムに選び, それらをまずFibnacci heapに入れ, 順番に頂点$u$を取り出して, $u$に"近い"頂点を緩和していくようなことをしています (arXiv版のSection 3参照).

4. 感想

講義のために最小費用流や最短経路問題の色々なアルゴリズムを勉強してみたのですが, マトロイドや離散凸といった一般的な枠組みを追求していく面白さとは違った面白さがあり, シンプルで賢いアルゴリズムは勉強していてとても楽しいですね.

5. 参考文献

• "An algorithm for Finding Shortest Routes from all Source Nodes to a Given Destination in General Network", Jin Y. Yen, Quarterly of Applied Mathematics, 27(4), pp.526--530, 1970.
• "Randomized Speedup of the Bellman–Ford Algorithm", Michael J. Bannister and David Eppstein, Proceedings of ANALCO, 2012.
• "Negative-Weight Single-Source Shortest Paths in Near-linear Time", Aaron Bernstein, Danupon Nanongkai, and Christian Wulff-Nilsen, Proceedings of FOCS, 2022.
• "Negative-Weight Single-Source Shortest Paths in Near-Linear Time: Now Faster!", Karl Bringmann, Alejandro Cassis, Nick Fischer, Proceedings of FOCS, 2023.

全点対最短経路問題の(+2)近似アルゴリズム

本記事はデータ構造とアルゴリズム Advent Calender 2021の1日目の記事です. Dor, Halperin, Zwick(2000)による重みなしグラフの全点対最短経路の近似を高速に計算するアルゴリズムを紹介します.

1. 全点対最短経路問題

与えられたグラフ$G=(V,E)$の各頂点対に対してその間の最短経路長を求める問題を全点対最短経路問題 (APSP; All-Pairs Shortest Paths) といいます. この問題の計算量は入力のグラフが辺重みを持つかどうかで変わってくると信じられています. 

この記事では$O(nm)$時間のアルゴリズムであっても計算量は$O(n^3)$と表記します (つまり$m$への依存は考慮せず, 最悪ケースの$m=n^2$だけを考える).

1.1. 辺重み有りの場合:

負閉路が存在しなければWarshall-Floyd法により$O(n^3)$時間で解くことができます. Warshall-Floyd法以外にも様々なアルゴリズムが提案されており, 以下の表のようななっています.

2021年12月現在ではWilliams(2014)による$n^3/\exp(\Omega(\sqrt{\log n}))$時間が(オーダーの意味で)最速ですが, $\exp(\sqrt{\log n})=n^{1/\sqrt{\log n}}=n^{o(1)}$なので, 任意の定数$\epsilon>0$に対して(重み付きグラフの)APSPを$O(n^{3-\epsilon})$時間で解くアルゴリズムは知られていません. つまり, およそ60年にわたってWarshall-Floyd法は$n^{o(1)}$倍の改善はあれど, $n^{\epsilon}$倍の改善はなされていません! このことから2010年代からは

任意の定数$\epsilon>0$に対して重み付きグラフのAPSPは$O(n^{3-\epsilon})$時間で解けない

という命題 (APSP予想と呼ばれる) が正しいという強い仮定の下で他の様々な問題の"精緻な"計算量下界が明らかにされてきました. これはちょうど, 「$\mathrm{P}\neq\mathrm{NP}$の下では最大カットは多項式時間で計算できない」などといったNP困難性の議論と同じ流れです (仮定する命題は大きく異なる). このアプローチは精緻な計算量 (fine-grained complexity) と呼ばれ, 以前の記事でも紹介しているので興味のある方は読んでみてください.

1.2. 辺重みなしの場合

辺重みなしグラフのAPSPは, 以前の記事で紹介したSeidelのアルゴリズムによって高速行列積に基づいた高速なアルゴリズムが知られており, 二つの$n\times n$行列の積の計算が$O(n^\omega)$時間でできるならば重みなしグラフのAPSPは$O(n^\omega\log n)$時間で計算できます. 2021年12月現在では$\omega < 2.37286$であることが知られています [Alman and Williams, 2021].

しかしながら, 高速行列積のアルゴリズムは複雑で, 実用上では素朴な$O(n^3)$時間アルゴリズムの方が高速であることが多いとされています. 従って高速行列積を使わない"組合せ的な"アルゴリズムの存在性が研究の議題になることもあります (e.g., [Bansal and Williams, 2012]).

2. 重みなしAPSPの近似アルゴリズム [DHZ00]

行列積が$O(n^2)$時間で計算できるならば重みなしグラフ上のAPSPも$O(n^2\log n)$時間で解けるのですが, 現状ではまだまだ程遠いです. 一方で行列積に基づかないAPSPアルゴリズムというのは自明な$O(n^3)$時間アルゴリズムとほぼ同程度の計算量しか達成されておらず, 例えば$O(n^{2.999})$時間のものは知られていません.

本記事では, 行列積を使わずにAPSPの(+2)近似を$O(n^{7/3}\log n)$時間で計算する乱択アルゴリズムを紹介します [Dor, Halperin, Zwick, 2000]. ここで, APSPの(+2)近似は全ての頂点対$u,v$に対して, $uv$間の最短経路長を$d(u,v)$と表すとして, $d(u,v)\leq D(u,v) \leq d(u,v)+2$を満たす$D(u,v)$を計算することを意味します. また, このアルゴリズムはランダムシードに依らず必ず$O(n^{7/3}\log n)$時間で終了しますが高確率(確率$1-O(1/n)$)で(+2)近似を計算します.

一般に近似精度が理論保証されている近似アルゴリズムというのはNP困難な問題に対して最適値の$c$倍以内の解を出力する (つまり, 近似比の保証がある) 多項式時間アルゴリズムを指すことが多いですが, この記事では$n^3$時間で解ける問題に対してadditive errorの保証がなされている近似アルゴリズムを考えているのに注意してください.

なお, $7/3\approx 2.333... < \omega$ なので, 現在知られている高速行列積アルゴリズムよりも高速な近似アルゴリズムになっています. ちなみに論文のタイトルは All Pairs Almost Shortest Paths というオシャレなタイトルで, 著者の1人であるUri Zwickはcolor codingの元論文の著者の1人です. なお, APSPの(+2)近似を$O(n^2)$時間で計算するアルゴリズムなどは分かってないです.

まず, 頂点集合を次数の大きさによって以下の三つの部分集合に分割します:

$L=\{v\in V\colon \deg(v)<n^{1/3}\}$

$M=\{v\in V\colon n^{1/3}\leq \deg(v)\leq n^{2/3}\}$

$H=\{v\in V\colon \deg(v)>n^{2/3}\}$

($L,M,H$はそれぞれlow, medium, high degreeだと思ってください). 各最短経路を

1. $L$内の頂点のみからなる最短経路
2. $H$の頂点を少なくとも一つ含む最短経路
3. $H$の頂点は一つも含まないが, $M$の頂点を少なくとも一つ含む最短経路

の三ケースに分けて考察します (下図).

2.1. $L$内の最短経路の計算

端点も含め全ての頂点が$L$に含まれるような最短経路を考えます. この最短経路は誘導部分グラフ$G[L]$に含まれるので, $G[L]$内で全頂点からBFSを行うことにより計算できます. このようにして得られた距離を$d_1(\cdot,\cdot)$とします (つまり, $u,v\in L$に対し$d_1(u,v)$は$G[L]$上の$uv$最短経路長).

$G[L]$に含まれる辺の本数は, $L$の次数制約より$O(n\cdot n^{1/3})$なので, このBFSは全体で$O(n^{7/3})$時間で完了します.

2.2. ランダムサンプリングによるHitting Setの計算

二つ目のケースを考える準備として, 別の記事で紹介したhitting setの概念を導入します. 頂点$v\in V$の隣接点の集合を$N(v)$と表すことにします. 頂点部分集合族$\{N(v)\colon v\in H\}$と$\{N(w)\colon w\in M\}$のhitting setをそれぞれ$S$と$T$で表します. 即ち, 

$S$は$S\cap N(v)\neq \emptyset$ ($\forall v\in H$)

$T$は$T\cap N(w)\neq\emptyset$ ($\forall w\in M$)

をそれぞれ満たします. $S$は以下のような簡単なランダムサンプリングで計算できます:

$S=\emptyset$からスタートし, 各頂点$u\in V$に対して独立に確率$p:=2n^{-2/3}\log n$で$u$を$S$に追加する.

$H$に属する各頂点$v$は$\deg(v)\geq n^{2/3}$を満たすので, $N(v)$が一つも$S$に追加されない確率を考えると

$$\Pr[N(v)\cap S=\emptyset] \leq (1-p)^{n^{2/3}}\leq \exp(-n^{2/3}p) \leq n^{-2}.$$

よってunion boundより, $\Pr[\exists v\colon N(v)\cap S=\emptyset]\leq n^{-1}$, つまり$S$は確率$1-n^{-1}$で hitting set になる. さらに, $|S|$は独立な確率確率変数(indicator)の和になるので, Chernoff boundにより, 確率$1-n^{-1}$で$|S|=\mathbf{E}[|S|]=\Theta(n^{1/3}\log n)$を満たします.

同様のランダムサンプリングにより, $T$はサイズ$\Theta(n^{2/3}\log n)$を満たすようにとることが高確率でできます.

2.3. $H$を通過する最短経路の近似

先ほど計算したhitting set $S$の各頂点$s\in S$に対し, $s$を始点としたBFSを計算することによって, $\{d(s,w)\colon s\in S,w\in V\}$を$O(|S|n^2)=O(n^{7/3}\log n)$時間で計算できます. ここで, 全ての頂点$u,v\in V$に対して $d_2(u,v)=\min_{s\in S}\{d(u,s)+d(s,v)\}$ とします. 

補題1.

$uv$間の最短経路であって$H$の頂点を含むものが存在するならば, $d(u,v)\leq d_2(u,v)\leq d(u,v)+2$.

(証明).

$h\in H$を通る$uv$最短経路を$P$とする. $S$はhitting setなので, $S\cap N(h)\in h'$とする(下図参照). このとき, 三角不等式より

$$d(u,v)\leq d_2(u,v)\leq d(u,h')+d(h',v)\leq d(u,h)+d(h,v)+2 =d(u,v)+2.$$

2.4. $H$を通らないが$M$を通る最短経路の近似 (方針)

このパートでは方向性だけ説明します.

$H$を通る最短経路は2.3で計算済みなので, 残りの最短経路は与えられたグラフから$H$を全て削除しても影響がありません. 残った辺は$L\cup M$に接続しており, $L\cup M$に属する頂点の次数は全て$\leq n^{2/3}$なので, 残った辺は$O(n\cdot n^{2/3})$本しかありません. つまりこのグラフはある程度は疎です! 

残ったグラフ上で2.2で計算したhitting set $T$に属する各頂点$t$に対し, $t$を始点としたBFSを行い, $\{d_{G'}(t,v)\colon t\in T,v\in V\}$ を$O(|T|n^{5/3})=O(n^{7/3}\log n)$時間で計算します. ここで$d_{G'}(\cdot,\cdot)$は$G'$上の最短経路を意味します. 

各頂点$u,v$に対し, $\min_{t\in T}\{d_{G'}(u,t)+d_{G'}(t,v)\}$を考えます. 補題1と同様の議論及び本パート冒頭の議論 ($H$間の辺の削除をしても, $H$を通過しない最短経路に影響はない) より, 以下が成り立ちます.

補題2.

元のグラフ上において, $uv$間の最短経路が$H$を通らないが$M$を通るとする. このとき, $d(u,v)\leq \min_{t\in T}\{d_{G'}(u,t)+d_{G'}(t,v)\}\leq d(u,v)+2$. 

(図による証明)

残念ながら, $\min_{t\in T}\{d_{G'}(u,t)+d_{G'}(t,v)\}$を愚直に計算すると$O(n^2|T|)=O(n^{8/3})$時間かかってしまいます. 以下では別のグラフを新たに構成することによってこの計算を工夫します.

2.5. $d_3(u,v)$の計算アルゴリズム

各頂点$x\in V$に対し, 辺重みつきグラフ$G_x=(V,E_x)$を以下のように構成します: $E_x=\emptyset$から始めます.

1. $L$に接続する辺を全て$E_x$に追加します. これらの重みは全て1とします. ここで追加される辺数は$L$の次数の条件より$O(n\cdot n^{1/3})$です.
2. 各$t\in  T$に対して$E_x$に辺$\{x,t\}$を追加します. この辺の重みを$d_{G'}(x,t)$とします. これらの重みは2.4で計算済みです. ここで追加される辺数は$O(|T|)=O(n^{2/3}\log n)$です.
3. $T$-$M$間で元のグラフに含まれている辺を全て重み1の辺として$E_x$に追加する. すなわち, 各$m\in M$と全ての$z\in T\cap N(m)$に対し ($T$はhitting setなのでこのような$z$は必ず存在します), 重み1の辺$\{m,z\}$を$E_x$に追加する. ここで追加される辺数を考えて見ましょう. 各$t\in T$の次数は高々$n^{2/3}$ ($H$の頂点は削除したから) なので, 辺数は高々$$|T|\cdot n^{2/3}=O(n^{4/3}\log n)$です.

上記の操作1~3で追加した辺を便宜上それぞれパターン1~3の辺と呼ぶことにします.

$G_x$上で$x$を始点としたDijkstra法を行うことにより, $\{d_{G_x}(x,v)\colon v\in V\}$を計算します. なお, 上記1~3の議論から$|E_x|=O(n^{4/3}\log n)$なので, このDijkstra法はFibonacci heapを用いた実装を使えば$O(|E_x|+n\log n)=O(n^{4/3}\log n)$時間で終わります (Fibonacci heapを使わなくても$O(n^{4/3}\log^2 n)$時間で終わります). 全ての$x\in V$に対して実行しても全体の計算時間は$O(n^{7/3}\log n)$です.

$d_3(u,v)=d_{G_u}(u,v)$とします. 以下が成り立ちます.

補題3.

元のグラフ上において, $uv$間の最短経路が$H$を通らないが$M$を通るとする. このとき, $d(u,v)\leq d_3(u,v)\leq d(u,v)+2$. 

(証明) 元のグラフよりショートカットするような辺を追加していないので, $d_{G_u}(u,v)\geq d(u,v)$は明らか. $d_{G_u}(u,v)\leq d(u,v)+2$を示す. 元のグラフ$G$上において, $L\cup M$内の$uv$最短経路であって$m\in M$を通るものを$P$とする. また, $P$が通過する$M$の頂点のうち一番最後に通過するものを$m$とする. すると, $m$以降の経路に含まれる辺は全て$L$に接続しているため, $G_u$においてパターン1で追加された辺からなる. すなわち$d(m,v)=d_{G_u}(m,v)$が成り立つ.

さらに, $G_u$上において, 

$$d_{G_u}(u,m)\leq d_{G_u}(u,t)+1\leq d_{G'}(u,t)+1\leq d_{G'}(u,m)+2=d(u,m)+2$$

が成り立つ. ここで, 最初の不等号における$t$はパターン3で追加された辺を利用した三角不等式で, 2つめの不等号はパターン2の辺の付け方の定義より成り立ち, 3番目の不等号は$G'$上の三角不等式から成り立ち, 最後の等号は$P$の性質 ($H$の頂点を通過しない) ことから得られます. 以上より

$$d_{G_u}(u,v)\leq d_{G_u}(u,m)+d_{G_u}(m,v)\leq d(u,m)+2+d(m,v)=d(u,v)+2$$

となり補題3を得ます.

2.6. アルゴリズムの概要

2.1〜2.5で計算した三つの距離$d_1,d_2,d_3$に対し, $\min\{d_1(u,v),d_2(u,v),d_3(u,v)\}$を出力します. $d_3$は必ずしも対称性 ($d_3(u,v)=d_3(v,u)$) が成り立つとは限らないので, 対称性を保証するために$d_3(u,v):=\min\{d_{G_u}(u,v),d_{G_v}(v,u)\}$を計算しても良いです. これまで議論してきたように, $d_1$から$d_3$は高々$O(n^{7/3}\log n)$時間で計算でき, しかも補題1~3により(+2)近似であることが保証されます.

3. 余談 (Spanner)

密な重みなしグラフ$G=(V,E)$が与えられた時, そのグラフの全域部分グラフ$S=(V,E_S)$のうち頂点間の距離をできるだけ保つようなものをspannerといいます. もう少しフォーマルな定義をいうと, 二つの正のパラメータ$\alpha,\beta$に対して以下の条件を満たすとき, 全域部分グラフ$S$は$G$の$(\alpha,\beta)$-spannerであるといいます.

任意の二頂点$u,v\in V$に対して$d_G(u,v)\leq d_S(u,v) \leq \alpha d_G(u,v)+\beta$.

なので, もし辺数$|E_S|=m_S$の$(\alpha,\beta)$-spannerが$T(n)$時間で計算することが出来るならば, APSPに対する$T(n)+O(nm_S)$時間近似アルゴリズムであって, multiplicative と additive の error がそれぞれ$\alpha,\beta$となるようなものを与えることができます. ちなみに$(1+1/k,0)$-spannerを$O(n^{1+1/k})$時間で計算するアルゴリズムが構成できます. 別の記事でも紹介したように, spannerは理論計算機科学でもかなり活発に研究されている分野で, 組合せ論におけるErdős Girth Conjectureと深いつながりがあります.

3. まとめ

重みなしグラフ上のAPSPに対して$O(n^{7/3}\log n)$時間の(+2)近似アルゴリズムを紹介しました. このアルゴリズムは現状の高速行列積よりもオーダー的に高速で, それほど複雑でもないため苦労せず実装できるものになっています.

アイデアとしては次数によって頂点集合を分割し, 高次数の頂点を経由する最短経路はhitting setを使って近似するというものになっています.

これら二つのアイデアはそれぞれ

1. 次数分離: スパースなグラフで役立つアルゴリズムのアイデア
   
2. 直径の1.5近似アルゴリズム
   

でも紹介しているので, 興味のある方は是非読んでみてください. 特にhitting setのアイデア (hitting set argument) は他にもgirthの近似[Ducoffe 2021]などにも応用される汎用的なアイデアなので, 非常に有用だと思います.

直径の1.5近似アルゴリズム

概要

$n$頂点$m$辺グラフの直径の1.5近似を$\tilde{O}(m\sqrt{n}+n^2)$時間で計算するAingworthら[1]のアルゴリズムを紹介します ($\tilde{O}(\cdot)$は$\mathrm{polylog}(n)$ factorを無視). 次にこのアルゴリズムの計算時間を改良したRoditty and Williams[2]による$\tilde{O}(m\sqrt{n})$時間乱択アルゴリズムを紹介します.

1. 導入

重みなし無向グラフ $G=(V,E)$ に対し二頂点$u,v$間の距離 $\mathrm{dist}(u,v)$ を最短経路長とします. このとき, $G$の直径 $\mathrm{diam}(G)$ は $\mathrm{diam}(G):=\max_{u,v\in V}\mathrm{dist}(u,v)$ で定義されます.

本記事のトピックは$n$頂点$m$辺グラフ$G$が与えられた時にその直径$\mathrm{diam}(G)$を計算する問題の計算量です. なお, 計算モデルは$O(\log n)$-word RAMを仮定します.

各頂点からBFSをすることによって$O(mn)$時間で計算することができます. 従ってこれよりも高速に計算できるかが焦点になります.

グラフの直径の計算は精緻な計算量(fine-grained complexity)においてかなり活発に議論されているトピックで, 直径3以下のグラフが与えられた時にその直径が$2$であるか$3$であるかを$O(m^{2-\epsilon})$時間で判別することはSETHと呼ばれる仮定の下では不可能であることが証明されています ($\epsilon>0$は任意の定数)[2]. すなわち, SETHを仮定すると直径の1.49近似を$O(m^{2-\epsilon})$時間で計算できないので, 例えば疎なグラフ ($m=O(n)$) では各点BFSによる$O(n^2)$時間アルゴリズムが最適であることを意味します. ここで, 直径の$c$近似 ($c>1$) というのは直径の$1/c$倍以上の値を出力することを意味します (例えば最大カットなどの文脈では$c^{-1}$近似の言い回しをよく使いますが, 直径近似研究の文脈では$c$近似と言うことが多く, 本記事もそれにならいます). ですので1.49近似アルゴリズムがあれば, 直径3のグラフに対してそのアルゴリズムを走らせると出力値は$\geq 3/1.49>2$となるため, 直径2と直径3のグラフを区別することができます.

それでは, SETHの下での近似困難性はタイトなのでしょうか? 具体的には, 疎なグラフに対し$O(n^{1.99})$時間で直径の1.5近似が計算できるかどうかが気になります. 

まず, 与えられるグラフ$G$が木である場合を考えてみます. このときはよく知られたdouble sweepと呼ばれる方法を使えば$O(m)$時間で直径の厳密値を計算できます. また, 一般グラフに対して適用すると, 直径の2近似が$O(m)$時間で計算できます.

本記事ではまず, Aingworthら[1]による$\tilde{O}(m\sqrt{n}+n^2)$時間で直径の"ほぼ"1.5近似を計算するアルゴリズムを紹介します. 具体的には彼らのアルゴリズムは $\lfloor 2/3 \mathrm{diam}(G) \rfloor \leq D \leq \mathrm{diam}(G)$を満たす値$D$を出力します. なお, 乱択を許しちょっと修正することによってAingworthらのアルゴリズムと同じ条件を満たす値$D$を$\tilde{O}(m\sqrt{n})$時間で出力する乱択アルゴリズムが知られています[2].

2. Aingworthらのアルゴリズム

Aingworthらによるアルゴリズム[1]を紹介します. double sweepを知っておくと感覚を若干ながら掴みやすいかもしれませんが, 知らなくても大丈夫です. また, 以下の説明ではBFSが出てきますが, これをDijkstra法に置き換えるとweighted graphに対して良い感じで動くと思うのですが, 私はあまり深く考えたことがないのでこの記事では割愛します. 以下, $O(m+n)$と書くのが面倒なので, $m\geq n$を仮定し, 単一始点BFSの時間計算量は$O(m)$と書いていきます.

アルゴリズム

1. 各頂点$v\in V$に対しBFSを行い$v$に近い順に頂点を$\sqrt{n}$個選び$S_v$とする (BFSを行い$\sqrt{n}$個の頂点に訪問したら途中で打ち切るようにする). なお, $v\in S_v$とする.

2. ステップ1で得られた集合族$(S_v)_{v\in V}$に対し, サイズ$O(\sqrt{n}\log n)$の hitting set $H\subseteq V$ を計算する (つまり, $H$は全ての$v\in V$に対し$H\cap S_v\neq\emptyset$を満たし, かつ$|H|=O(\sqrt{n}\log n)$を満たす集合. 存在性は後で示します).

3. $H$の各頂点$h\in H$に対し, $h$を始点としたBFSを行い, 全ての$h\in H,v\in V$に対し$\mathrm{dist}(h,v)$を計算する. その後, $D_1:= \max_{h\in H,v\in V}\mathrm{dist}(h,v)$とする.

4. $H$から最も遠い頂点を$u$とする. つまり$u$は$\min_{h\in H}\mathrm{dist}(h,v)$を最大にする$v$である.

5. ステップ4で得られた頂点$u$に対し, 各$x\in S_u$を始点としてそれぞれBFSを行い, 得られた距離の中での最大値を$D_2$とする. すなわち, $D_2:= \max_{x\in S_u,v\in V}\mathrm{dist}(x,v)$. 最後に$\max\{D_1,D_2\}$を出力.

計算時間:

- ステップ1では各頂点に対してBFSを行うので, 各$v\in V$に対し$O(|E(G[S_v])|)=O(n)$時間かかるので計$O(n^2)$時間です.

- ステップ2を考えます. まず存在性を議論します. 各頂点を独立に確率$q:=\frac{2\log n}{\sqrt{n}}$で選んで得られるランダムな頂点部分集合を$R$とします. このとき$\mathrm{E}[|R|]=2\sqrt{n}\log n$であって, Chernoff boundにより$|R|$は高確率 (具体的には確率$1-n^{-1}$で) $|R|=O(\sqrt{n}\log n)$となります. また, 固定した頂点$v\in V$に対して「$S_v$の全ての頂点が選ばれない」確率は

$(1-q)^{|S_v|}\leq \exp(-q|S_v|)=\exp(-2\log n)=n^{-2}$

なので, 確率$1-n^{-2}$で$R$は$S_v$の頂点を一つ以上含みます. 従って, union boundにより確率$1-1/n$で$R$は hitting set になっています. 以上より, $R$は高確率でサイズ$O(\sqrt{n}\log n)$の hitting set になってます. もし条件を満たすhitting setが存在しないならばこの確率は$0$にならなければいけないので, 存在性を示せたことになります. また, ランダムにとってくればほぼ確実に条件を満たすhitting setになっているので, 乱択を許せば$O(n)$時間でステップ2を計算できます. また, greedyにやれば ($H=\emptyset$からスタートして$S_{v_1}$の頂点を一つ選んで$H$に追加し, 暫定の$H$がhitしていないような$S_{v_2}$を選んで...を繰り返す) 条件を満たすようにできます.

-ステップ3は$O(|H|m)=O(m\sqrt{n}\log n)$時間で計算できます.

-ステップ4はステップ3の情報を使えば計算できます.

-ステップ5は単一始点BFSを$\sqrt{n}$回行うので$O(m\sqrt{n})$時間で計算できます.

以上より, 全体の時間計算量は$O(m\sqrt{n}\log n+n^2)=\tilde{O}(m\sqrt{n}+n^2)$です.

近似率の解析

アルゴリズムの出力を$D=\max\{D_1,D_2\}$とします. 単純のため, $\mathrm{diam}(G)$は3の倍数とします. 直径を達成する二頂点を$a,b\in V$とします (つまり$\mathrm{dist}(a,b)=\mathrm{diam}(G)$).

二つのケースを考えます:

(1). ある $h\in H$, $v\in V$ が存在して $\mathrm{dist}(h,v)\geq \frac{2}{3}\mathrm{diam}$ がなりたつならば, $D_1$が直径の1.5近似になっているので大丈夫です.

(2). 全ての$h\in H$, $v\in V$に対して $\mathrm{dist}(h,v)<\frac{2}{3}\mathrm{diam}$が成り立つとします. ステップ4で計算した頂点$u$に対し, $\ell:=\min_{h\in H}\mathrm{dist}(u,h)$とすると, $u$の定義より任意の$v\in V$は距離$\leq\ell$で$H$内のいずれかの頂点にたどり着くことができます. ここで$\mathrm{dist}(a,h)\leq \ell$を満たす頂点$h\in H$をとります. すると三角不等式より

$\mathrm{diam}\leq\mathrm{dist}(a,h)+\mathrm{dist}(h,b)<\ell+\frac{2}{3}\mathrm{diam}$,

すなわち $\ell>\frac{\mathrm{diam}}{3}$ を得ます (下図). また, $H$ は $(S_v)_{v\in V}$の hitting set であるため, 特に $H \cap S_u \neq \emptyset$ です. 従って (*)$u$ から距離$\leq\frac{\mathrm{diam}}{3}$ なる頂点は全て$S_u$に属しているはずです (行間補足: $S_u$ は $u$から近い順に頂点を取っていて, 距離$\ell>\frac{\mathrm{diam}}{3}$ だけ離れている頂点 ($H$との共通部分) も取っているのでそれより近い頂点は全て$S_u$に入っていなければならない).

ここで, $D_2 < \frac{2}{3}\mathrm{diam}$と仮定します. このとき, $D_2$の定義より, 全ての頂点は$u$から距離$<\frac{2}{3}\mathrm{diam}$となっています. 次に, $\mathrm{diam}=\mathrm{dist}(a,b)\leq \mathrm{dist}(a,u)+\mathrm{dist}(u,b) < \mathrm{dist}(a,u)+\frac{2}{3}\mathrm{diam}$ より, まとめると $\frac{1}{3}\mathrm{diam}<\mathrm{dist}(u,a) < \frac{2}{3}\mathrm{diam}$ を得ます.

最短$ua$経路上の頂点であって, $u$からの距離がちょうど$\mathrm{diam}/3$となる頂点を$p$とします. $\mathrm{dist}(a,u)=\mathrm{dist}(a,p)+\mathrm{dist}(p,u)<\frac{2}{3}\mathrm{diam}$となるため, $\mathrm{dist}(a,p)<\frac{1}{3}\mathrm{diam}$です. さらに, $p$の取り方から主張(*)より, $p\in S_u$です.

三角不等式より,

$$\mathrm{diam}=\mathrm{dist}(a,b)\leq \mathrm{dist}(a,p)+\mathrm{dist}(p,b)<\mathrm{dist}(p,b)+\frac{1}{3}\mathrm{diam} $$

よって, $D_2\geq \mathrm{dist}(p,b) > \frac{2}{3}\mathrm{diam}$ となってしまい仮定に矛盾します. よって$D_2\geq \frac{2}{3}\mathrm{diam}$ となります (証明終).

ちょっとした改善

最後に, Roditty と Williams [2] がによるちょっとした改善を紹介します. 基本は上記のAingworthらのアルゴリズムと同じですが計算時間を$\tilde{O}(m\sqrt{n}+n^2)$ から $\tilde{O}(m\sqrt{n})$にしています. つまりグラフが疎であれば高速になります.

アイデアは単純で, ステップ1を行わず, ステップ2で得られた hitting set $H$ を, 全頂点を独立に確率$2\log n/\sqrt{n}$で選んで得られるランダムな頂点部分集合とします. そしてステップ5では$u$に対して$S_u$を計算するだけです. Aingworthらのアルゴリズムではステップ1の計算で$O(n^2)$時間がかかっていましたが, そもそも$H$が $(S_v)_{v\in V}$ の hitting setであれば良いだけなので, 計算時間の解析でも紹介したようにランダムにとるだけで良いので省くことができます.

感想

賢すぎて草

参考文献

[1] D. Aingworth, C. Chekuri, P. Indyk, and R. Motwani. Fast estimation of diameter and shortest paths (without matrix multiplication). SIAM J. Comput., 28(4):1167–1181, 1999 (Preliminary version is in SODA96).

[2] L. Roditty and V. V. Williams, Fast Approximation Algorithms for the Diameter and Radius of Sparse Graphs, STOC13.

小記事: 3SUMのO(n^2)時間アルゴリズム

 3SUMと呼ばれる次の判定問題を考えます:

入力 : $n$個の整数値 $A=(a_1,\dots,a_n)$ ただし$a_i\in[-n^c,n^c]$.

出力 : YES iff 相異なる三つ組 $(i,j,k)$ が存在して $a_i+a_j+a_k=0$

この問題は

for i in range(n):

    for j in range(n):

        for k in range(n):

            if a[i]+a[j]+a[k]==0: return YES

return NO

というコードによって$O(n^3)$時間で解けますが, 少し工夫して次のようにすれば $O(n^2\log n)$ 時間で解けます.

L = []

for i in range(n):

    for j in range(n):

        L.append(a[i]+a[j])

L.sort()

for k in range(n):

    if -a[k] is in L:     // 二分探索をして-a[k]がLにあるかどうかを判定

        return True

return False

しかし尺取り法に基づいた次のアルゴリズムを用いると$O(n^2)$時間で解くことができます (簡単のため, 入力配列$A$の中身は全て相異なる要素であるとします).

A.sort()

for i in range(n):

    l = 0

    u = n-1

    while True:

        if u==l: return False

        if A[i]+A[u]+A[l]==0: return True

        elif A[i]+A[u]+A[l]>0: u-=1

        elif A[i]+A[u]+A[l]<0: l+=1

            

3SUMという問題は計算幾何学や精緻な計算量理論において非常に重要な問題の一つとされていて, 任意の定数$\epsilon>0$に対して$O(n^{2-\epsilon})$時間では解けないと予想されています.

一方で$\mathrm{polylog}(n)$倍の改善はなされています.

次数分離: スパースなグラフで役立つアルゴリズムのアイデア

単純無向グラフ$G=(V,E)$上で長さ$k$の閉路を探したり直径を求めたくなったりする時があります。そして時には$G$はスパースになっていることを利用したくなります。今回の記事では、スパース性を利用した単純ながら面白く賢いトリック (次数分割: high-degree low-degree technique)を紹介し、これに基づいてグラフに含まれる三角形 (長さ3の閉路) の個数の数え上げを行う非自明なアルゴリズムを紹介します。次数分離という和名は私が勝手に名付けたものであり世間で浸透しているものではないので注意してください (もし世間で浸透している名前をご存知の方がいらっしゃいましたら教えていただけると嬉しいです...!!）。このアイデアはAlon, Yuster, Zwick [1] によるものです。

$|E|=m$, $|V|=n$とします。一般には$m=O(n^2)$になりえますが、今回は$G$がスパースであることを仮定してるので$m=O(n)$を想定します。ただ次数分離は一般の$m$に対しても適用でき、例えば$m=O(n^{1.5})$とかでも似たような解析で計算時間を見積もることができます。

次数分割の基本アイデア

頂点$v$の次数を$\deg(v)$とします。$\alpha\in[0,1]$をパラメータとし、

$$H=\{v\in V: \deg(v)>n^\alpha\},$$

$$L=V\setminus H$$

とします ($H$は high degree、$L$は low degree のイメージです）。すると、$2m=\sum_v \deg(v)$ なので

$$2m \geq \sum_{v\in H} \deg(v) > |H|n^\alpha$$

より$|H|<2m/n^\alpha$となります。特に$m=O(n)$ならば$|H|=O(n^{1-\alpha})$となり、$n$と比べて小さくなります。これらの頂点部分集合は$O(m)$時間で得られます (計算モデルは$O(\log n)$bit-Word-RAM モデルを仮定)。

応用. 三角形(長さ3の閉路$C_3$)の数え上げ

入力としてグラフ$G=(V,E)$が与えられ、$G$に含まれる長さ$3$の閉路$C_3$と同型な部分グラフを数え上げる問題を解いてみましょう。以下の二つの素朴なアルゴリズムを考えます。

・アルゴリズム1:

$G$の隣接行列の三乗を計算し対角成分の総和を見ると、長さ3の閉じたウォークの個数となり、これは三角形の個数の$3!=6$倍になっています。従って$O(n^{\omega})$時間で数え上げられます。 ($\omega<2.373$は高速行列積の指数)。

・アルゴリズム2:

各頂点に対して隣接点を列挙して隣接点のペアで隣接しているものがあるかどうかを見れば$O(\sum_{v\in V}\deg(v)^2)$時間で数え上げが解けます。$(2m)^2=\left(\sum_{v\in V}\deg(v)\right)^2\geq \sum_{v\in V}\deg(v)^2$より、このアルゴリズム2は$O(m^2)$時間なので、例えば$m=O(n)$なら$O(n^2)$となりアルゴリズム1以上に高速になります ($\omega=2$なら同等のオーダーになる)。

ここで、冒頭に説明したアイデアを考えます。適当なパラメータ$\alpha\in[0,1]$に対して上で定義した頂点部分集合$H,L$を得ます。入力グラフ$G$に含まれる三角形$C_3$は以下の二つのケースに分類されます:

ケース1. 三頂点全てが$H$に含まれる場合

ケース2. 一つ以上の頂点が$L$に含まれる場合

ケース1は誘導部分グラフ$G[H]$内に含まれる三角形の個数を数え上げれば良く, $|H|=O(n^{1-\alpha})$より素朴なアルゴリズム1を用いれば$O(n^{\omega(1-\alpha)})$時間で数え上げることができます。

ケース2は$L$に含まれる全ての頂点に対し素朴なアルゴリズム2を用いることで列挙できるので、数え上げることができます。このとき、$k$頂点が$L$に含まれるような三角形 $(k=2,3)$ は$k$回列挙されるため適当な値で割ったりする必要があるので気をつけてください。列挙にかかる時間は、$L$に属する頂点の次数が高々$n^\alpha$であるので $\sum_{v\in L}\deg(v)^2 \leq n^{1+2\alpha}$ となります。

以上より、$m=O(n)$の場合は任意の$\alpha\in[0,1]$に対して三角形の数え上げを

$$O(n^{\omega(1-\alpha)}+n^{1+2\alpha})$$

時間で解くことができます。$\alpha=\frac{\omega-1}{\omega+2}$とおくとこの計算時間は$O(n^{3-\frac{6}{\omega+2}})$時間になります。現状では$\omega<2.373$なのでこの計算時間は$O(n^{1.628})$となり, $O(n^2)$よりも高速になっています。もし将来$\omega=2$になることがあれば$O(n^{1.5})$となります。

上の解析では$m=O(n)$を仮定していましたが一般の$m$に対しても同様の解析で$O(m^{3-\frac{6}{\omega+2}})$時間で計算できることを示せます ($|H|=O(m^{1-\alpha})$となるので)。

まとめ

スパースなグラフに対して、ある閾値を設定しそれより高い次数を持つ頂点集合とそれ以外の頂点集合に分割するテクニックを紹介し、その応用として三角形の数え上げ問題を取り上げました。本質的には、このような分離を行い、低次数の頂点の探索は近傍を見るコストが低いことを利用し、高次数の頂点の探索はまとめて処理できる賢いアルゴリズムを利用する、というだけのことです。

三角形と同様の解析を用いれば、例えば入力グラフが四角形(長さ4の閉路)を部分グラフとして持つかどうかの判定問題を$O(m^{4/3})$時間で組合せ的 (高速行列積という飛び道具を使わず初等的に）解くことができます。$H$と$L$に分解したあと、$H$内部の四角形を鳩の巣原理に基づくアルゴリズム (昔の記事参照) で判定し、$L$の頂点を含む四角形は近傍探索を駆使して判定します。結構実装も簡単でシンプルかつ汎用性がありそうなアイデアなので、競技プログラミングで出ても良さそうな気が個人的にはします（というか既にある？？）

参考文献

[1] N. Alon, R. Yuster, and U. Zwick, Finding and counting given length cycles, Algorithmica, 1997.

閉路検出問題の計算複雑性

グラフ $H$ を固定し, $H$-Detection という問題を考えます. この問題はグラフ $G$ が入力で与えられたときに $G$ が $H$ を部分グラフとして含むかどうかを判定する問題です.

例えば, $H$ が長さ$|V(G)|$ の閉路だとするとこれはハミルトン閉路問題になるのでNP-困難です. 一方で $H$ が辺数 $|V(G)|/2$ のマッチングならば完全マッチングを含むかどうかの判定問題になるので, 多項式時間で解けます. それでは, どのような$H$に対して $H$-Detection は多項式時間で解けるのでしょうか？

この話はグラフアルゴリズムの分野でよく研究されている分野であり, 結論から言うと未解決なのですが, $H$のサイズをパラメータとしたときに$H$-Detectionは「木幅がunbounded ならばW[1]-困難である」という予想が有ります[1] (逆に木幅がboundedならばcolor codingでFPT時間で解けます). 本記事では $H$を様々な長さの閉路としたときどのような結果が知られているかについてまとめていきます. 基本的に$n=|V(G)|$, $m=|E(G)|$ とします.

$k$-cycle

一般の$k$に対しては color coding というアルゴリズムによって, $2^{O(k)}\cdot n^{\omega}$時間で解けます ($n\times n$の行列積が$O(n^{\omega})$ 時間でできるとします. 2019年6月1日現在では, $\omega\approx 2.373$ です).

$k$が定数だとしましょう. 実は, 現状では$k$の偶奇によって計算量が変わります. $k$が奇数のときはこれまで知られている最速のアルゴリズムは $O(n^{\omega})$ 時間です. 例えば $k=3$ の場合は隣接行列の3乗を計算して体格成分だけを見れば良いですね. 一方で$k$が偶数のときは, $O(n^2)$で解くアルゴリズムが知られています[2]. 例えば$k=4$の時は以下の単純なアルゴリズムで四角形を見つけることができます. ただし $N(u)$は頂点$u$の隣接頂点の集合とします.

このアルゴリズムでは配列 $\mathrm{T}[i][j]$ は二頂点$ij$間の長さ2のパスをカウントしていて, もし$\mathrm{T}[i][j]\geq 2$ なる$ij$が見つかったらその瞬間に四角形が有ると判定してアルゴリズムを終了します. 配列$\mathrm{T}$の各要素は高々2回しかアクセスされないはず (2回アクセスしたらその瞬間にアルゴリズムはtrueを返して終了する) なので, 計算時間は$O(n^2)$になるわけです.

また, 辺数$m$に関しては近年にもう少し早いアルゴリズムが知られていたり[2]して, 自分が把握する限りをまとめるとこんな感じになります.

ハミルトン閉路

ハミルトン閉路はTSPで使うDPを考えると$O(2^n n^2)$時間で解くことができるのですが, それでは$O(1.99^n)$時間で解けるのでしょうか？ (乱択を許せば)できることが知られています[4]. アルゴリズムではまず, 入力のグラフ$G$を使ってある多変数多項式 $P$ を構成し, $P$が恒等的に$0$かどうかを乱択で判定するというものです. ここで構成する多項式$P$は, 「$P$が恒等的に0 iff グラフがハミルトン閉路を持つ」という性質を満たすようなものを上手く構成します. したがって多項式が恒等的に0かどうかの判定が出来る決定的アルゴリズムがあればハミルトン閉路の検出も決定的に出来るようになります. CyganらによるFPT algorithmの本[5]の10.4.2に分かりやすく解説されているので興味のある方はどうぞ.

[1]. M. Grohe, The complexity of homomorphism and constraint satisfaction problems seen from the other side, J. ACM, 2007.

[2]. R. Yuster and U. Zwick, Finding even cycles even faster, SIAM J. Discrete Math, (1997).

[3] S. Dahlgaard, M. B. T. Knudsen, M. Stöckel, Finding even cycles faster via capped k-walks, In Proc. of STOC2017.

[4]. A Björklund, Determinant sums for undirected Hamiltonicity, SIAM J. COMPUT., 2014.

[5]. M. Cygan, F. V. Fomin, Ł. Kowalik, D. Lokshtanov, D. Marx, M. Pilipczuk, M. Pilipczuk, S. Saurabh, Parameterized Algorithms, Springer Publishing Company, 2015.

マトロイドの基の数え上げの決定的多項式時間近似アルゴリズム

$\mathcal{M} = (E,\mathcal{B})$ をマトロイドとします. ここで $E$ は台集合で$|E|=n$を満たし, $\mathcal{B}$ は基の全体とします. マトロイドの独立性オラクルが与えられたとき, $|\mathcal{B}|$ を求めるという問題を考えます.

マトロイドの数え上げは例えばグラフの全域木の総数やグラフのreliability (分配関数みたいなもの)の計算を求める問題を含んでおり, かなり多くの研究結果があります. 今回の記事では

Log-Concave Polynomials I: Entropy and a Deterministic Approximation Algorithm for Counting Bases of Matroids

という論文の内容について, 簡略に紹介したいと思います. ここで載せたリンクはarXivのものですが, この論文はFOCS 2018に採択されています[1]. 論文そのものはlog-concave polynomialという凄い道具について解析をしており, その一つの応用としてマトロイドの数え上げ近似をするという内容なのですが, 凄い道具の解析は凄い難しく私も咀嚼しきれていないため, 本記事では凄い道具をどのように使って数え上げ近似をするのかについて焦点を絞ります.

論文の概要: log-concave polynomialに関する理論を築いた. その帰結として, ランク $r$ に対して, 独立性オラクルの下でマトロイドの基の総数の$2^{O(r)}$-近似を出力する多項式時間決定的アルゴリズムを与えた. 二つのマトロイドの共通基の個数も同様の近似比のアルゴリズムを与えた. 近似解はマトロイドの基多面体上のある最適化問題を解くことによって得られる.

0. マトロイドの基の数え上げについて

マトロイドの基の数え上げは割と長く研究されていて, これまでは Balanced matroid と呼ばれるクラスのマトロイドに対してはMCMCに基づく多項式時間乱択近似アルゴリズムが知られていました. この辺の話題はNIPS2018のチュートリアルでも関連した話題が取り扱われており, 注目の高さが伺えます.

グラフの全域木の総数は行列木定理により簡単に数えることが出来ます. ではこれを拡張してマトロイドの基は数えられるのでしょうか？

結論から言うとマトロイドの基の数え上げはとても難しいと信じられています. 具体的に言うとバイナリマトロイドの基の数え上げは#P-困難であることが知られています. 従って方向性としては近似アルゴリズムを考えるのが自然です. しかし近似率にも下界が知られており, 独立性オラクルモデルを多項式回呼び出すどんな決定的アルゴリズムも近似比は (P≠NPとは独立に) $2^{\Omega(n/(\log n)^2)}$になることが知られています [2]. この結果の系として, ランク $r\gg\log n$であれば近似比の下界

$2^{\Omega(r/(\log n)^2)}$          ...(1)

が導出されます.

冒頭で述べた通り, balanced matroidと呼ばれるクラスのマトロイドに対してはMCMCに基づく近似アルゴリズムが知られています[3]. 用語を知っている人向けに言うと, このアルゴリズムは
Balanced matroid ならば $\mathcal{M}$ とそのいかなるマイナーに対してもnegative correlationが成り立つので, base exchange graph上のランダムウォークにより基の(近似的な)一様サンプリングが可能になる. そこで Jerrum, Valiand and Vaziraniによる有名な結果 (サンプリングと数え上げ近似の関係) を使うと基の数え上げの近似が出来る.

という流れになっています. 噛み砕いて言うとbalanced matroidだと base exchange graph ($\mathcal{B}$が頂点集合で, 基交換公理によって行き来できる基同士を辺で結んだ巨大なグラフ) がエクスパンダーになっていて定常分布（ここでは一様分布）に早く収束し, 基の一様サンプリングが出来ます. そして実は一様サンプリングが出来れば数え上げで良い近似が出来るという結果が元からあるのでそれを使って数え上げの近似を行う, という流れです. 一般のマトロイドに対して base exchange graph がエクスパンダーかどうかはMihail-Vazirani 予想と呼ばれ, 長年未解決でしたが, 近年解決されたそうです (2019年1月時点ではarXiv上ですが).

重要な用語の定義を述べていきます. サイズ $n$ の台集合$E$上の集合族 $\mathcal{S}$ 及び $\mathcal{S}$上の確率測度 $\mu$ に対して

$g_\mu (x_1,\ldots,x_n) := \sum_{F\in \mathcal{S}} \mu(F)\prod_{e\in F} x_e$

という関数を定めます. この関数を $\mu$ の generating polynomial と呼びます. 以下では特に$\mathcal{S}$に属する任意の集合がサイズ$r$であるような場合のみを考えます (例えばマトロイドの基族). これまで知られている結果として, この関数が $\mathrm{Im}(x_i)>0$ for every $i=1,\ldots,n$ なる任意の点 $\mathbf{x}=(x_1,\ldots,x_n)$に対して $g(\mathbf{x})\neq 0$ を満たす (この性質を real stable と呼びます) ならば, $\mu$に対して自然に定義される$\mathcal{S}$上のランダムウォークは定常分布$\mu$に収束するというのがあります. これによって数え上げの近似っぽいことが出来ます.
マトロイドの基上の一様分布に対するgenerating polynomialがreal stableであればそのマトロイドがbalancedであるということが知られています.

ここまでで述べたように, 基の数え上げに対するアルゴリズムとしてはMCMCに基づくものが基本形でした (特定のマトロイドに関してはad-hocに数え上げるアルゴリズムは存在する). しかし今回紹介する論文では, 決定的なアルゴリズムで近似比が $2^{O(r)}$ となるものを提案します. これは近似比の下界(1)と比較すると右肩にちょっとだけギャップがありますが, とても面白いアルゴリズムです.

1. エントロピーを用いた$|\mathcal{B}|$の近似

情報理論や物理学ではエントロピーという概念がよく知られています. ここでは, 基族上の確率測度 $\mu:\mathcal{B}\to\mathbb{R}_{\geq 0}$ に対して

$H(\mu)=\sum_{F\in\mathcal{B}}\mu(B)\log\frac{1}{\mu(B)}$

をエントロピーと呼びます. もし$\mu$が$\mathcal{B}$上の一様分布ならば, $H(\mu)=\log|\mathcal{B}|$ が成り立ちます. もしもランク$r$のマトロイドの基族$\mathcal{B}$上の一様分布$\mu$に対し,
$H(\mu)-r \leq \beta \leq H(\mu)$
を満たす $\beta$ が計算できたとすると, $\exp(\beta)$ は $|\mathcal{B}|$ の$2^{O(r)}$-近似になっています. この$\beta$を計算することを考えます.

台集合$E$の要素 $e\in E$ に対し, $\mu_e:=\mu(\{F\in\mathcal{B}:F\ni e\})$ と定めます. つまり$\mu_e$はアイテム$e$を含む確率 (marginal probability) に対応しています. 論文中では次の結果を証明しています:

Theorem 1.
もしも$\mu$の generating polynomial $g_\mu$ が log-concaveならば, そのエントロピー $H(\mu)$ について, $H(\mu)\geq \sum_{e\in E} \mu_e\log\frac{1}{\mu_e}$ が成り立つ.

Theorem 2.
マトロイドの基族$\mathcal{B}$上の一様分布$\mu$に対応する generating polynomial は log-concave である.

また, エントロピーの劣加法性より, 以下の事実が簡単に確認できます:

Proposition 3.
$H(\mu) \leq \sum_{e\in E}\left(\mu_e\log\frac{1}{\mu_e}+(1-\mu_e)\log\frac{1}{1-\mu_e}\right)=\sum_{e\in E}H(\mu_e)$.

ここで, 関数 $f:x\mapsto (1-x)\log\frac{1}{1-x}$ は $0\leq x\leq 1$に対して $f(x)\leq x$が成り立ちます. 従って, ランク$r$のマトロイドの基族$\mathcal{B}$を考え, $x=\mu_e$を代入して$e\in E$について総和をとると,

$\sum_{e\in E}\mu_e\log\frac{1}{1-\mu_e}\leq \sum_{e\in E}\mu_e \leq \sum_{F\in\mathcal{B}}r\cdot \mu(F)=r$

が得られます. ここで Theorem 1 と Proposition 3 を組み合わせると,

$\sum_{e\in E}(1-\mu_e)\log\frac{1}{\mu_e}\geq H(\mu)-r$

となるので, 系として

Corollary 4.

ランク$r$のマトロイドの基族$\mathcal{B}$上の分布$\mu$に対し, その generating polynomial が log-concave ならば,

$H(\mu)\geq \sum_{e\in E}\mu_e\log\frac{1}{\mu_e}\geq H(\mu)-r$

が成り立つ.

2. 提案アルゴリズム

前節の結果より, marginal probability $\mu_e$ を計算出来れば良いのですが, それはそもそも「アイテム$e$を含む基の個数は？」という問題の答えになっているため, 基の数え上げよりも難しそうです. そこで$\mu_e$を計算せずに済む方法を考えます.

マトロイド$M$の基族 $\mathcal{B}$ 上の一様分布 $\mu$ に対して, $n$次元ベクトル $(\mu_e)_{e\in E}=(\mu_1,\ldots,\mu_n)$ を考えると, 基多面体 $\mathcal{P}_M$ に属する点の一つになっていることがわかります. ここで, 基多面体とは, それぞれの基の indicator vector の凸包として定義されます.

$(\mu_1,\ldots,\mu_n) = \frac{1}{|\mathcal{B}|}\sum_{F\in \mathcal{B}} 1_F$

より確かに marginal distribution のなすベクトルは基多面体$\mathcal{P}_M$内の点に含まれています. そこで, 次の最適化問題を考えます:

maximize: $\sum_{i=1}^n H(p_i)$
s.t.
     $(p_1,\ldots,p_n)\in \mathcal{P}_M$.

Proposition 3より, この問題の最適値 $\beta$ は求めたいエントロピー$H(\mu)=\log|\mathcal{B}|$の上界になっています. しかも, 目的関数は$p$について凹関数になっており, マトロイドの独立性オラクルがあれば$\mathcal{P}_M$上で分離オラクルが構成でき, 楕円体法を用いて多項式時間で解くことができます.

次に, 最適解 $p\in\mathcal{P}_M$ を考えます. 実は(論文中で示されているのですが) marginal distributionが$p$となるような$\mathcal{B}$上の分布$\hat{\mu}$で, generating polynomial $g_{\hat{\mu}}$が log-concave となるものが構成出来ます (つまり, $p_e=\hat{\mu}_e$が成り立つ). この$\hat{\mu}$に対して Corollary 4を適用すると,

$H(\hat{\mu})\geq \sum_{e\in E}\hat{\mu}_e\log\frac{1}{\hat{\mu}_e}=\sum_{e\in E}p_e\log\frac{1}{p_e}\geq \sum_{e\in E}H(p_e)-r\geq H(\mu)-r$.

更に, $H(\hat{\mu})\leq H(\mu)$ が成り立つ (一様分布$\mu$がエントロピーを最大化する) ので, 結果として, 最適化問題の最適値$\beta$に対して

$\beta-r \leq H(\mu)=\log|\mathcal{B}|\leq \beta$

が成り立つので, 近似比 $2^{O(r)}$ が成り立ちます.

3. まとめ

アルゴリズムは単純ですが, 近似比の保証において「generating polynomial $g_\mu$ が log-concaveならば〜」というステートメントが多く出ています. ここで紹介した論文の貢献としては log-concave function の解析にあって, 提案アルゴリズムはその応用先の一つでしかありません (とはいっても重要な結果です). 著者たちは更に別の続編とも言うべき論文 (こちらは2019年1月時点ではarXivに上がっているだけですが) で, Mihail-Vazirani 予想を解決しマトロイドのbase の数え上げに対してMCMCに基づく FPRAS を提案しており, 今後も離散の世界における log-concave 性が発展していくことになるかと思います (連続の世界ではlog-concave はすでに注目されていていろんなことが知られています).


[1]. Log-concave polynomials, entropy, and a deterministic approximation algorithm for counting bases of matroids, N. Anari, S. O. Gharan, and C. Vinzant, In Proceedings of the 59th Annual Symposium on Foundations of Computer Science (FOCS) , 2018.

[2]. On the problem of approximating the number of bases of a matroid, Y. Azar, A. Z. Broder, and A. M. Frieze, Information Processing Letters Volume 50, Issue 1, 1994.

[3]. Balanced matroids, T. Feder, and M. Mihail, In Proceedings of the 24th Annual ACM Symposium on Theory of Computing (STOC), 1992.

color-codingと脱乱択化

乱択アルゴリズムとは, ランダムネスの力を上手く活用したアルゴリズムのことで, ここでは「必ず多項式時間で終わるが間違った解を出力するかもしれない」というタイプのアルゴリズムを考えます（モンテカルロアルゴリズム）.

本記事では, パラメタ化計算量の分野で有名な乱択アルゴリズムとしてcolor coding[4]を紹介し, その「脱乱択化」の方法 (perfect hash family) を紹介します.

1. color coding

1.1. 背景 (k-PATH)

グラフGが与えられたとき, それが長さkのパスを持つかを判定する問題(k-PATH)を考えます. パスは同じ頂点を二度通らないことに注意してください. $k=n-1$ のとき, ハミルトンパスの判定問題を表現するのでこの問題は $k$ を入力としたときはNP完全となります. 一方で $k=O(1)$ としたときは $n^k$ 時間, すなわち多項式時間で解くことが出来ます. それでは, $k=O(\log n)$ の時は多項式時間で解けるでしょうか? $k=\sqrt{n}$ のときはどうでしょうか? 結論から言うと, $k=O(\log n)$ならばk-PATHは多項式時間で解けますが, $k=\omega(\log n)$の場合は(ETHの下で)多項式時間で解けないだろうと予想されています. color codingを使うと, k-PATHを $2^{O(k)}\cdot n^{O(1)}$ 時間で解くことができます. 以下のそのアルゴリズムを紹介していきますが, 同じ内容は吉田 悠一さんによるこちらの記事にも非常に分かりやすく載っているので紹介させていただきます.

1.2. アルゴリズム

1. 与えられたグラフ $G=(V,E)$ の各頂点をランダムに $k+1$ 色で塗ります. つまり, 各頂点に対して$\{1,\ldots,k\}$から一様ランダムに一つの値を選択しその値の色で塗ります.




2. 色付きグラフ $G=(V,E)$ がカラフルな k-path を含むかどうかを判定します. カラフルなk-path とは, 長さkのパスであって各頂点の色が全て異なる (i.e., $k+1$ 種類の色全てを含む) ようなものです. もしカラフルな k-path を含んでいたら, 元のグラフ $G$ は k-path を含むので終了し, そうでなければステップ1に戻ります.

1.3. アルゴリズムの解析

与えられたグラフ $G=(V,E)$ が k-path を含んでいたとして, その頂点を $v_1,\ldots,v_{k+1}$ としましょう. これらの $k+1$ 個の頂点がカラフルになる確率は

$\frac{k!}{k^k} \approx \sqrt{2\pi k} \exp(-k) \leq 3^{-k}$

となります (スターリングの近似を用いました). したがって上述のステップ1-2を $9^k$回繰り返したとき, パス$(v_1,\ldots,v_{k+1})$が毎回カラフルにならない確率は

$(1-3^{-k})^{9^k} \leq \exp(-3^{-k}\cdot 3^{2k}) = \exp(-3^k)$

となり, 非常に小さくなります (ここでは $1+x\leq \mathrm{e}^x$を用いました). つまり非常に高確率で, $9^k=2^{O(k)}$回の繰り返しの中で一回は k-path はカラフルになり, ステップ2で検出されるので高確率でk-PATHを解くことができることになります.

1.4. カラフルなk-パスの検出

動的計画法を使ってカラフルなk-パスの検出を $O(2^k\cdot n^3)$ 時間で解くことが出来ます. 二頂点$u,v\in V$および色の集合$S\subseteq \{1,\ldots,k+1\}$に対して配列として

$\mathrm{T}[u][v][S]:=$頂点$u$から$v$への, $S$-カラフルなパスが存在すればtrue, そうでなければfalse.

と定めます. ここで $S$-カラフルなパスとは, $|S|$個の頂点からなるカラフルなパスで, 各$c\in S$に対して色$c$で塗られた頂点を一個だけ含むようなものです (下図参照).

色つきグラフ$G$において$N(v)$で頂点$v$の隣接点の集合を表し, $\chi(v)\in \{1,\ldots,k+1\}$で$v$の塗られた色を表すこととします.

すると, $\mathrm{T}[u][v][S]$に対して以下の漸化式が成り立ちます.

$\mathrm{T}[u][v][S]=\bigvee_{w\in N(v)}\mathrm{T}[u][w][S-\{\chi(v)\}]$.

basic caseとして $\mathrm{T}[v][v][\{\chi(v)\}]=\mathrm{true}$ と定めます. また, 自明ですが$\{\chi(u),\chi(v)\}\not\subseteq S$ならばT[u][v][S]=falseです. これにより, カラフルk-パスの判定は$O(2^kn^3)$時間で解けるので, アルゴリズム全体では$9^k\cdot n^3=2^{O(n)}\cdot n^{O(1)}$時間で解けます ($n$の右肩の値はおそらく$2$になると思いますが, 大雑把に見積もっても高々$3$なのは明らかなので今回はそう記しました).

1.5. その他の問題への応用

k-PATH以外にもk-CYCLE (長さkの閉路の検出問題) といった問題も color coding を応用したアルゴリズムによって $2^{O(k)}\cdot n^{O(1)}$時間で解くことができます. 一般に, 入力で与えられたグラフ$G$が頂点数$k$の別のグラフ$H$を含むかどうかは 全探索により$O(n^{k+O(1)})$時間で解けますが, color codingを使うと $2^{O(k)}\cdot n^{O(\mathrm{tw})}$ 時間で解くことが出来ます (twは$H$の木幅です). したがってtwが定数ならば$2^{O(k)}\cdot n^{O(1)}$時間で解けます.

1.6. k-PATHに対するより早いアルゴリズム (advanced topic)

P≠NP予想より強いETHという予想を仮定すると, k-PATH問題は$2^{o(k)}\cdot n^{O(1)}$時間では解けないことが知られています[1]. 今回紹介したcolor coding では$9^k\cdot n^3$時間でしたが, 解析をちゃんとやると $O(\mathrm{e}^kn^3\log n)$時間でwith high probability (i.e., 適当な定数$\epsilon>0$に対し確率 $\geq 1-n^{-\epsilon}$) で正解するアルゴリズムになるはずです. では, 例えば$2.3^k\cdot n^{O(1)}$時間で解けるのでしょうか?

結論から言うと, $2^{k}\cdot n^{O(1)}$時間の乱択アルゴリズムが存在します[2]. このアルゴリズムは

1. 「グラフがk-pathを含まない iff Pは恒等的に0」となる多項式Pを構成する (但しガロア体$\mathrm{GF}(2^{\lceil 4k \rceil})$上の多項式).

2. Schwartz-Zippelの補題を用いてPが恒等的に0かどうかを乱択で判定する (その際に包除原理を用いる).

という流れです. Schwartz-Zippelの補題については相馬 輔 さんによる こちらのページで証明も含めて紹介されています. 更に, 現在では$2^{0.75k}\cdot n^{O(1)}$ 時間アルゴリズムも知られています[3].

この$2^{k}\cdot n^{O(1)}$時間アルゴリズムは様々な道具が登場し非常に面白いものなので, 機会があれば紹介したいと思います.

2. 脱乱択化 (derandomization)

脱乱択化とは文字通り「ランダムネスを排除する」ことです. これまで考えていたアルゴリズムは高確率で正解を計算するものでしたが, 脱乱択化をすると必ず正解を計算するアルゴリズムになります. 計算量理論では, 多項式時間の乱択モンテカルロアルゴリズムで(高確率で)解ける問題のクラス (BPP) が多項式時間で決定的に解ける問題のクラス(P) と等しいかどうかは未解決ですが, 多くの研究者はP=BPPだと信じているようです. P=BPPならば多項式時間乱択アルゴリズムは多項式時間のまま脱乱択化できる, ということになります. ここでは, それに関連したポジティブな結果の一つとして, 先ほど紹介したcolor codingの脱乱択化の方法について紹介します. キーワードはperfect hash familyです.

2.1. color coding is revisited

color codingに基づくアルゴリズムは

(i) グラフを$k$色で塗る.

(ii) カラフルな部分グラフを動的計画法で$2^{O(k)}\cdot n^{O(1)}$時間で見つける.

(iii) 見つけたい部分グラフが(i)でカラフルになる確率が$\geq 2^{-O(k)}$となるので, $2^{O(k)}$回(i)と(ii)を繰り返す.

というものでした.

即ち, 「見つけたい部分グラフが(i)でカラフルになるような塗り分け」をランダムネスを使わず決定的な方法で見つけられたら脱乱択化になります. 

2.2. Perfect hash family

入力で与えられるグラフを$G=(V,E)$とし, その頂点数を$n=|V|$とします. そして塗り分けの集合 $\mathcal{F}$ を考えます. $[m]:=\{1,\ldots,m\}$と表すことにすると, 各要素 $f\in \mathcal{F}$は写像$f:[n]\to [k]$となり, $f(v)$は頂点$v$の色を意味します. つまり$\mathcal{F}$は写像の集合になります. 

定義.

$\mathcal{F}$ が以下の条件(＊)を満たすとき, $\mathcal{F}$は $(n,k)$-perfect hash familyと呼びます.

(＊) サイズ$k$の任意の部分集合$S\subseteq [n]$ of $|S|=k$ に対して, ある$f\in\mathcal{F}$が存在して, $f(S)=[k]$となる (i.e., 塗り分け$f$の下では$S$はカラフルになる).

$\mathcal{F}$ を $(n,k)$-perfect hash family としましょう. グラフ$G=(V,E)$が入力で与えられたとき, $\mathcal{F}$を構成して全ての$f\in \mathcal{F}$のそれぞれに対してcolor codingの動的計画法のアルゴリズムを使ってカラフルな部分グラフを判定すれば, $T(n,k)+|\mathcal{F}|\cdot 2^{O(k)}\cdot n^{O(1)}$時間で決定的に解くことが出来ます. ここで$T(n,k)$は$\mathcal{F}$の構成にかかる時間です.

$(n,k)$-perfect hash familyについては以下の結果が知られています[4]:

定理.

$(n,k)$-perfect hash family $\mathcal{F}$で

$|\mathcal{F}|\leq \mathrm{e}^kk^{O(\log k)}\log n$

を満たすものを $\mathrm{e}^kk^{O(\log k)}n\log n$ 時間で構成することができる.

これにより, color codingの脱乱択化を$O(\log n)$時間程度のロスでやることが出来ます.

参考文献

[1] D. Lokshtanov, D. Marx, and S. Saurabh. Lower bounds based on the Exponential Time Hypothesis. Bulletin of the EATCS, 84, pp.41–71, 2011.

こちらのサーベイ論文は(題目通り)ETHやSETHの基づいたlower boundとその証明テクニックを幅広く紹介しており, 大変良い論文です. SETHについては以前の記事で紹介しています.

[2] R. Williams, Finding a path of length $k$ in $O^*(2^k)$ time, Journal Information Processing Letters, 109-6, pp. 315-318, 2009.

color codingではk-PATHを$\mathrm{e}^k n^{O(1)}$時間でしたが, こちらの論文はSchwartz-Zippelの補題を用いて$O(2^{k}n^{O(1)})$時間でk-PATHを解くアルゴリズムを提案しています.

[3] M. Cygan, F. V. Fomin, L. Kowalik, D. Lokshtanov, D. Marx, M. Pilipczuk, M. Pilipczuk S. Saurabh, Parameterized Algorithms, Springer International Publishing, 2015.

パラメタ化計算量に関するアルゴリズムを幅広く紹介しています. FPTといった計算量クラス基本的な定義やW-階層といった困難性のクラスやそれらに基づくlower boundについても載っています. 各chapterの後ろには豊富な演習問題とBibliographic notesが載っており, ここには計算量の改善の歴史などが載っています. 誤植もあまりなく, 英語も読みやすいため個人的には非常にオススメの本です.

[4] N. Alon, R. Yuster, U. Zwick, Color-coding, J. ACM, 42-4, pp. 844-856, 1995.

color codingを提案した非常に有名な論文です. 論文中では$k$-CYCLE (長さ$k$の閉路の検出) を用いてcolor codingのアルゴリズムを説明していますが, それのみならず, 木幅が定数のグラフの検出, perfect hash familyを用いた脱乱択化についても書かれています.

高速な全点対最短経路アルゴリズム

1. 全点対最短経路問題の計算量の外観

　グラフ $G$ が与えられたときに全点対最短経路(APSP : all pair shortest paths)を求めたいときがあります. そんなときはほとんどの人がワーシャル・フロイド法または各点ダイクストラなどを使うかと思います. しかしこれらのアルゴリズムは頂点数 $n$ に対して最悪 $O(n^3)$ 時間かかってしまうので, $n=1000$くらいでもちょっと厳しい感じがあります.

　実は辺重み付きグラフのAPSPの計算量には多くの歴史があり、ワーシャルフロイド法以降様々なアルゴリズムが提案されてきました.

APSPのアルゴリズムの歴史

　しかし、任意の $\epsilon>0$ に対して $O(n^{3-\epsilon})$ 時間でAPSPを解くアルゴリズムはこれまで知られていません. なので現在ではそのようなアルゴリズムは存在しないのではないかと予想されています. ちなみにSETHの仮定の下で計算量の下界が示されている問題が多くありますが, APSPのこの下界をSETHの仮定の下で示すのは難しいだろう, という根拠が2016年にCarmosino, Gao, Impagliazzoらによって示されました.

2. 重み無しグラフの全点対最短経路

　以上のことから, 辺重み付きグラフのAPSPを$O(n^{2.99})$で解くのは絶望的であると言えます. では, 辺重み無しグラフのAPSPについてはどうなのでしょうか？以降は「グラフ」は辺重み無しグラフを指すこととします.

　頂点数$n$, 辺数$m$のグラフのAPSPは各頂点から幅優先探索を行うことによって $O(nm)$ 時間で解くことができます. したがってグラフが密で, $m=\Omega(n^2)$ であるならば計算時間は $O(n^3)$ だけかかってしまいます. 

　結論から言うと辺重み無しグラフのAPSPは$m$に依存せず $O(n^{\omega}\log n)$ 時間で解くアルゴリズムが存在します. そのアルゴリズムは[Seidel, 1992]で初めて与えられました. 非常にシンプルなので, 本記事ではこのアルゴリズムについて解説します.

3. Seidel のアルゴリズム

　前提として入力として与えられるグラフ$G$は連結であるとします. したがってその直径は高々$n$です. 更に, 入力において$G$は隣接行列として与えられるとします. そしてアルゴリズムの出力は距離行列$D$とします. ただし$D_{ij}=ij$間の最短経路長と定義します.

　グラフ $G$ に対し, 任意の辺$uv,vw \in E(G)$ に対して 新たな辺$uw$を$G$に追加して得られるグラフを$\tilde{G}$ と定めます. 気持ちとしては$\tilde{G}$ は「全ての長さ2のパスに対してその端点を結ぶ」ことを出来るだけ行って得られるグラフです.

追加した辺を赤く塗っています. 三角形がめちゃくちゃ増えます.

$A$を$G$の隣接行列とします. すると$\tilde{G}$の隣接行列$\tilde{A}$は

になります ($i\neq j$に対し$A^2$の$ij$成分は$ij$間の長さ2のパスの本数を表すからです). したがって, $A^2$を$O(n^{\omega})$時間で計算することで$\tilde{A}$が計算できるので, $\tilde{G}$を得ることが出来ます. 

　この$\tilde{G}$に対して再帰によってAPSPの距離行列$\tilde{D}$を得ます. $\tilde{G}$の直径は$G$の直径の半分くらいになっているので, 再帰を繰り返した結果, 最終的には$G$の直径が$1$(=完全グラフ)となり, この時は(再帰のベーシックケースとして)出力$D$は非対角成分が全て1の行列となります. 更に, 再帰の深さは$\log n$くらいです. というのも$G$の直径が高々$n$で, 一回の再帰で直径が半分くらいになるからです.

　さて, 再帰の出力$\tilde{D}$を使って$G$の距離行列$D$を得ましょう. 下図のように, $\tilde{G}$の$ij$間の最短経路を考えます. 明らかに, 元のグラフ$G$における$ij$間の最短経路長$d_{ij}$が奇数であれば, $\tilde{G}$における$ij$間の最短経路長$\tilde{d}_{ij}$に対して $d_{ij} = 2\tilde{d}_{ij}-1$ が成り立ちます (下図(a)参照. $\tilde{G}$では長さ2のパスを1ホップで通過出来るので, 出来るだけ新しい辺を利用した方が経路長が短くなる). 一方で$\tilde{d}_{ij}$が偶数であれば, $d_{ij}=2\tilde{d}_{ij}$が成り立ちます (下図(b)). 以上より, $d_{ij}$のパリティが分かれば$\tilde{d}_{ij}$を用いて$d_{ij}$を復元できます.

元の最短経路長のパリティが分かれば復元できる.

　ここで, $\tilde{G}$上で ① $ij$間に(a)の図のような最短経路が存在する場合, $d_{ij}=2\tilde{d}_{ij}-1$となります. 一方で ② そのような最短経路が存在しない場合(つまり任意の$ij$間の最短経路が(b)のように赤い辺のみで構成されている場合), $d_{ij}=2\tilde{d}_{ij}$となることに注意します. つまり$d_{ij}$を復元するためには①と②の場合を判別する必要があります.

　そのために行列$X:=\tilde{D}A$を考えます. これは行列積なので$O(n^{\omega})$時間で計算できます. 行列積の定義より$X_{ij}$は

$X_{ij}=\sum_{k=1}^n \tilde{d}_{ik}\cdot A_{kj} = \sum_{k \in N(j)} \tilde{d}_{ik}$

となります. ここで$N(k)$は$G$における$k$の隣接点の集合, $\tilde{d}_{ik}$は$\tilde{G}$における$ik$間の最短経路長です. 

　②を仮定します. すると任意の$k\in N(j)$に対し, $\tilde{d}_{ik} \geq \tilde{d}_{ij}$ が成り立ちます (そうでないとすると, $\tilde{G}$上の$ik$間の最短経路に辺$kj$を足したものが$ij$間の最短経路になりますが, それは②に反するからです. 下図参照).

②で考えるケース. どの$k\in N(j)$も$\tilde{G}$上の$ij$間の最短経路上に含まれていない.

したがって

$\sum_{k \in N(j)} \tilde{d}_{ik} \geq \tilde{d}_{ij}\cdot \mathrm{deg}(j)$.

となります. ここで$\mathrm{deg}(j)=|N(j)|$は頂点$j$の$G$における次数です.

　次に①について考えます. この時はある$k\in N(j)$に対して$k$を含むような$\tilde{G}$上の$ij$間最短経路が存在します (下図).

　この$k$に対して $\tilde{d}_{ik} = \tilde{d}_{ij}-1$ が成り立ちます. 更に他の$k'\in N(j)$ に対しても $\tilde{d}_{ik'} \leq \tilde{d}_{ik}+1 = \tilde{d}_{ij}$が成り立ちます ($k$と$k'$間は$\tilde{G}$の辺で繋がっています. 下図参照).

したがって

$\sum_{k \in N(j)} \tilde{d}_{ik} < \tilde{d}_{ij}\cdot \mathrm{deg}(j)$.

が成り立ちます.

以上より, 

① $\iff$ $X_{ij}=\sum_{k \in N(j)} \tilde{d}_{ik} < \tilde{d}_{ij}\cdot \mathrm{deg}(j)$.

② $\iff$ $X_{ij}=\sum_{k \in N(j)} \tilde{d}_{ik} \geq \tilde{d}_{ij}\cdot \mathrm{deg}(j)$.

が成り立つので, 

として$d_{ij}$を求めることが出来ます. これでめでたく出力すべき$D=(d_{ij})$を構成することが出来ました.

　アルゴリズム全体の流れとしては

1. $G$の直径が$1$ならば自明な解を出力.
2. そうでないときは, $A^2$を計算することで$\tilde{A}$を計算し, 再帰に投げることで$\tilde{D}$を得る.
3. $X=\tilde{D}A$を計算し, $X$と$\tilde{D}$を用いて各$ij$に対して$d_{ij}$を計算する.
4. $D=(d_{ij})$を出力.

という感じになります. まず再帰の深さは$O(\log n)$で, ステップ2, 3においては行列積を計算しているので計算時間は$O(n^{\omega})$となり, ステップ4では$O(n^2)$かかるので, 全体の計算時間は$O(n^{\omega}\log n)$となります.

4. まとめ

　辺に重みがある場合と無い場合でAPSPの計算量には大きなギャップがあります. 辺重み付きグラフでは$O(n^{2.99})$では解けなさそうである一方, そうでないグラフに対しては$O(n^{\omega}\log n)=O(n^{2.38})$時間で解くアルゴリズムが知られています. もちろん辺重みが無いスパースなグラフに対しては依然として幅優先探索を用いた$O(nm)$の方が高速ですが, 密な場合はSeidel のアルゴリズムの方が高速であると言えます.
　では, 重み付きグラフに対して$O(n^{2.99})$時間でAPSPを計算するアルゴリズムは存在するのでしょうか？これは$P\neq NP$ほどではないですが, 計算量理論においてかなりホットな話題なので, みなさんもぜひ挑戦してみてください.

問題を解きやすいグラフクラスとその拡張

NP-困難な問題はP≠NPの仮定のもとでは多項式時間で解けません. しかしNP-困難な問題を解きたいという場面が少なからずあります. そのようなときには近似を考えたり入力インスタンスのクラスを制限する、ということを考えることが多々あります.

今回のトピックはグラフに対するNP-困難問題に対し, 「グラフクラスを制限することによって多項式時間で解く」アプローチについて広く浅く列挙します.

- 樹状幅を制限したグラフクラス (bounded tree-width)

　樹状幅とはもともと, Robertson・Seymour によるグラフマイナー定理の証明に用いられた重要な概念であり, グラフの「木っぽさ」を表す指標になっています. グラフクラス(有限グラフの無限個の集合)として樹状幅がある定数で抑えられるようなクラスのことをbounded tree-widthと呼びますが, このクラスに対してはクラスカルの定理を用いることでマイナー順序がwell-quasi-orderであることを示し, そうでないクラスについては「樹状幅がboundできない」という情報からグラフの構造上の特徴を示し, 曲面に「だいたい」埋め込んだ後に色々して証明している(らしい)です.
　一方でNP-困難な問題は木に対しては効率よく解けるということがよくあり, 例えば最大独立点集合問題は良い例です. したがって「木っぽい」グラフクラスであるbounded tree-widthに対しても最大独立点集合問題を効率よく解けるんじゃないか、という観点から研究が流行りました.
　樹状幅が高々kであるようなグラフに対しては $O(2^k n)$ 時間で幅が小さい樹状分解を求めることができるので、$k$を定数とみなすとこれは線形時間になっています. これを用いて樹状分解上での動的計画法や分割統治法を用いて様々な問題を解きます. このような(樹状幅などのパラメータを定数とみなしたときに多項式時間で解けるかどうかを議論する)枠組みのことを FPT (Fixed-Parameter Tractability) といいます.
　また, NP-困難な問題以外でも「グラフの直径の計算」「グラフの平均距離の計算」といったクラスPの問題(この２つの問題は $O(n^3)$ で解けるが $O(n^{2.999})$ 時間で解けるアルゴリズムが知られていない) に対しても同様のアプローチを考えることができて, 平均距離の計算を $O(2^{k \log (k)}\cdot n^{1+o(1)})$ 時間で行うアルゴリズムが2009年に提案されています. また, 2016年のSODAで FPT in P の概念が提唱され, (僕の中で)話題になりました.

- 平面的グラフ

　グラフクラスを平面的グラフに限ると高速になるという話がたくさんあります. 例えばグラフの最大カット問題は一般にはNP-困難ですが, 平面グラフに対しては多項式時間で解くことができます.
　一方でグラフクラスを平面的グラフに限ってもハミルトン閉路問題や最大独立点集合問題はなおNP-困難です. 平面的グラフ全体のクラスはbounded tree-width ではないので, 樹状幅によるFPTが通用しません. 実際, $k\times k$の格子グラフの樹状幅は$k$であることを示すことが出来ます.
　平面的グラフの辺数は最大で$3n - 6$なので, スパースです. 辺数$m$のグラフの直径を$o(m^2)$ 時間で計算するアルゴリズムがSETHの仮定のもとで存在しないことが知られています. つまり, スパース(i.e. $m=O(n)$)なグラフの直径には(SETHのもとで)$\Omega(n^2)$時間が必要であることになります. しかし, 平面的グラフの直径・平均距離を $O(n^{11/7} \cdot \mathrm{poly}\log n)$ 時間で行うアルゴリズムが2017年に提案されています(この論文はSODA2017でBest Paper Awardを受賞).

- bounded local tree-width

　bounded local tree-width というグラフクラスに触れます. このクラスに属すグラフ$G$ は次の条件を満たしています:

ある関数 $f(r)$ が存在して, 任意の頂点 $v$ に対して, $v$ から半径 $r$以内にある頂点から誘導される $G$の部分グラフの樹状幅 が $f(r)$ で上から抑えられている.

bounded local tree-width のイメージ.

　局所的な部分の樹状幅が functionally bounded であるということです. 実は平面的グラフは bounded local tree-width です. 実際, $f(r) = 3r-1$ が成り立ちます (この事実は全然自明じゃないです). また, bounded local tree-width なグラフ$G$に対しては $\mathrm{tw}(G) \leq f(\mathrm{diam})$ が成り立ちます. (ここで$\mathrm{diam})$は $G$ の直径です).

したがって,  bounded local tree-width というグラフクラスは bounded tree-width を拡張して平面的グラフを含むようにしたというようなグラフクラスになっています. このクラスに対するFPTのような枠組みは僕はあまり見かけたことはありませんが, 面白い研究題材になるんじゃないかと思っています.

余談 : bounded tree-width と bounded local tree-width のマイナーによる特徴づけ.

グラフマイナー理論において, bounded tree-width をマイナーの言葉で特徴づけしている次の定理が知られています:

グラフクラス $\mathcal{C}$ が bounded tree-width でないことの必要十分条件は, 任意の平面的グラフ $H$ に対して$H$をマイナーとして含むグラフ $G\in \mathcal{C}$ が存在することである.

　この定理は片方の向き (右から左) は自明です. 左の主張を仮定します. $H$ を $n\times n$ の格子グラフであるとすると, $H$をマイナーとして含む $G\in \mathcal{C}$ が存在しますが, $H$の樹状幅は$n$なので, $\mathrm{tw}(G)\leq n$ となります. $n$は任意にとれるので, $\mathcal{C}$は bounded tree-width ではありません.
　この定理の凄いところは左から右です. 実は, 「任意の$n$に対してある$k$が存在して, 樹状幅$k$以上の任意のグラフは$n\times n$の格子グラフをマイナーとして含む」という定理が知られていて, 更に「任意の平面的グラフ$H$に対して$H$をマイナーとして含むような格子グラフ$H'$が存在する」ということが簡単に言えます. この２つの事実を組み合わせることによって左から右の主張を証明することができます. 左の主張を仮定します. 任意の平面的グラフ$H$に対してある$n$が存在して, $n\times n$ 格子グラフ$H'$は$H$をマイナーとして含みます. この $n$に対してある$k$が存在して, 樹状幅$\geq k$の任意のグラフは$H'$, すなわち$H$をマイナーとして含みます. 左の主張の仮定より, $\mathcal{C}$は樹状幅$\geq k$のグラフ$G$を必ず持つので, 右の主張における所望のグラフをこの$G$が得られます.

　では, bounded local tree-width をマイナーの言葉で特徴づけすることはできるのでしょうか？ 実は, 次の定理が知られています[Eppstein, 2000]:

minor-closed なグラフクラス $\mathcal{C}$ が bounded local tree-width でないことの必要十分条件は, 任意のApex graph $H$ に対して$H$をマイナーとして含むグラフ $G \in \mathcal{C}$ が存在することである.

Apex graph とは, 「頂点を一つ取り除くことによって平面的になるようなグラフ」のことです. Apex graph の全体は平面的グラフを含みます. 例えば $n\times n$格子グラフに新たな頂点を一つ追加し, 新たな頂点を格子グラフの$n^2$個の頂点に接続して得られるグラフ$H$を考えましょう. このグラフは直径が2ですが, 樹状幅$=n$なので, bounded local tree-width ではありません. このことから定理の片方の向き (右から左) が言えます.

全点間最短経路の計算をちょっとの工夫でちょっと速くする話

概要

最近研究で重みなし連結正則グラフの全点間の最短経路を求めることが多く、ちょっと工夫したらちょっとだけ早くなったのでそれについて簡単に書きます。

素朴な方法

一番素朴な方法としては幅優先探索があります。
全ての頂点から幅優先探索をして全点間の最短経路を全て計算する方法です。

アルゴリズム

ところである点から幅優先探索を行うと図のように三角形みたいな図が出てくると思います。



　　　　　　　　　　　　　　

　幅優先探索をやっているとキューに頂点を取り出しては突っ込んでを繰り返して最短距離を計算していくわけですが、よく考えるとこの図の最下層にいる頂点たちはキューから取り出されてもその後何もせずに放り投げられてしまうことになります。したがって、キューに突っ込んだ頂点の個数がグラフの頂点数とちょうど等しくなったらキューを空にして幅優先を終了し、別の頂点からまた幅優先する」というアルゴリズムでちょっとだけ計算が速くなることが期待されます。
　考え方としては簡単でセコイように見えますが、グラフの形状によっては紫の部分の頂点がすごく多かったりすることがあるので、密なグラフに対してはそこそこ高速になるんじゃないかと期待できることと思われます。

数値実験

ということで幾つかの種類のグラフに対し「全点間最短経路長」を計算する以下に示す２つのプログラムを用意し、比較しました。
A. 素朴な幅優先探索
B. 上に述べた工夫をしたもの

用意したグラフは次のモデルによって生成されたグラフです。
1. Erdős–Rényi model (最も代表的なランダムグラフのモデル)
2. Configuration Model (ランダム正則グラフの代表的なモデル)

　　　　　　　　　　　

Erdős–Rényi modelによって生成されたグラフに対するプログラムの実行時間の比較です。
横軸はモデルのパラメータpの値で、右に行くほど辺の本数が多くなっています。頂点数は10000で固定しています。 縦軸は実行時間で、単位は[ms]です。紫の線が素朴な方法(A)で、緑の線が工夫したもの(B)です。


Configuration Modelによって生成されたグラフでの比較です。横軸は頂点の次数を表しているため、右に行くほど線形で辺の本数が増えていきます。縦軸は先ほどと同様に実行時間で単位は[ms]です。頂点数は10000です。なお、１番目の実験で生成したグラフよりも辺数は圧倒的に少ないためこちらの方が実行時間はかなり短くなっています。紫が素朴な方法で、緑が工夫した方法です。


二つの結果を見ると、意外とこの工夫によって早くなることが分かると思います。

考察・まとめ

・ 空間計算量を変えずにちょっとした工夫でBFSを早くするという話でした。

・ 恐らく緑の方では実行時間はグラフの直径によって大きく左右されるのでほぼ横ばいみたいな感じになっていると思います。
・ 直径をDだとすると、#{ 最短距離がDとなるような頂点ペア } が多いほど工夫が効果的になると思います。
・ したがって例えばグラフが木だったりするとぜんぜん効果がないと思います。

・ 実データとかでは実験してないのでその辺を含めてもうちょっと調べてみたいと思います。
・ さすがにこれくらいの工夫は誰か思いついてるでしょ...