単純ランダムウォークの到達時間と全訪問時間のまとめ

概要

連結無向グラフ上の単純ランダムウォークの到達時間(hitting time)と全訪問時間(cover time)に関する様々な結果を紹介します. 前半で諸定義を書き後半は軽いサーベイになっています.

1. 用語と定義

この章は基本的な用語の定義だけなのでよく分からなければ読み飛ばしても大丈夫です. 本記事では有限集合$V$上の離散時間マルコフ連鎖$(X_t)_{t\geq 0}$をランダムウォークと呼ぶことにします. 便宜上, 集合$V$を頂点集合と呼びその要素$v\in V$を頂点と呼ぶことにします. 遷移確率が時刻に依存せず一定であれば, $(X_t)$は斉時的である (time-homogeneous) といいます. 以下では断りなく$(X_t)$はtime-homogeneousなランダムウォークを表すこととし, その遷移確率行列を$P$と表記します. 即ち, $P\in[0,1]^{V\times V}$であって$(X_t)$は

$$\Pr[X_{t+1}=v\mid X_t=u] = P(u,v)$$

の遷移規則を満たす確率変数列です. 数直線上のランダムウォークなどを考えると$V=\mathbb{Z}$は無限集合になることもあり, この時は$P$は行列ではなく線形作用素と見なすことになりますが, ここではそのようなケースは考慮しません.

$P$と初期点$X_0$を与えると$(X_t)$の軌道の分布が定まります. $X_t$の分布を$x_t\in[0,1]^V$とすると, $x_{t}=x_{t-1} P = \dots = x_0 P^t$を満たします. 大袈裟に言ってしまえば, $(X_t)$の分布を調べることは$P^t$の性質を調べることと同義なので, 線形代数が重要になってくることは想像に難くないと思います. 特に, $P$がその定常分布から誘導される計量(内積空間)において自己随伴的であることを仮定することが多く (これを可逆性; reversibility), この場合は特に実固有値を持ち, 特にCourant–Fischerの定理を使うことによって混交時間 (mixing time)などを抑えることができます. これらの基礎をちゃんと学びたい人は Reversible Markov Chains and Random Walks on Graphs や Markov Chains and Mixing Times などの教科書がオススメです.

例(単純ランダムウォーク)

無向グラフ$G=(V,E)$の頂点$u\in V$の次数を$\deg(u)$とします. 

\begin{align*} P(u,v) = \begin{cases} \frac{1}{\deg(u)} & \text{if }uv\in E,\\ 0 & \text{otherwise}\end{cases}\end{align*}

によって定義される$P$によって定まる$(X_t)$を単純ランダムウォーク(simple random walk; SRW) といいます. 言い換えると, 隣接頂点から一様ランダムに頂点を選びそこに遷移する過程で, 最も有名なランダムウォークだと思います. $X_t$の分布のベクトル$x_t\in[0,1]^V$に対し, $\lim_{t\to\infty}x_t$は収束するでしょうか?

結論から言うとNOとなる例が存在します. 無向二部グラフ上の単純ランダムウォークを考えてみましょう. すると左右の頂点集合に交互に移動することになるので, $x_t$は収束しないことがすぐわかります.

このようなランダムウォークの周期性 (periodicity) を除去するために, 確率$1/2$で今の位置にとどまるような以下の$P$で定まるランダムウォークを考えてみましょう:

\begin{align*} P(u,v) = \begin{cases} \frac{1}{2} & \text{if }u=v,\\ \frac{1}{2\deg(u)} & \text{if }uv\in E,\\ 0 & \text{otherwise}.\end{cases}\end{align*}

自己ループをくっつけたので, 先ほどの二部グラフのような周期性は除外され, 実は$x_t$の集束することが証明できます ($G$の無向性は必要). このようなランダムウォークは遅延単純ランダムウォーク (lazy simple random walk; LSRW) と呼びます (lazyの和訳でよく使われるものを知らないのでとりあえずここでは遅延と呼ぶことにします).

1.1. 到達時間, 全訪問時間, 混交時間

$V$上のランダムウォーク$(X_t)_{t\geq 0}$を考えます. 二つの頂点$u,v$に対し確率変数$$\tau_{\mathrm{hit}}(u,v)=\min\{t\geq 0\colon X_t=v,X_0=u\}$$

を考えます. これは頂点$u$からスタートして頂点$v$に初めて到達するまでに要する時間を意味します. 定義より明らかに$\tau_{\mathrm{hit}}(u,u)=0$です.

$$t_{\mathrm{hit}}=\max_{u,v\in V}\mathrm{E}[\tau_{\mathrm{hit}}(u,v)]$$

を到達時間(hitting time)と呼びます. 以降, ``$\tau$"が付くものは確率変数で``$t$"がつくものはそれの期待値だと思ってください.

次に, 全訪問時間(cover time)を定義します. 頂点$u$に対し確率変数

$$\tau_{\mathrm{cov}}(u)=\min\{t\geq 0\colon \{X_0,\dots,X_t\}=V,X_0=u\}$$

を定義します. これは$u$を始点としたランダムウォークが全ての頂点に少なくとも1回ずつ到達するのに要した時間を表す確率変数です.

$$t_{\mathrm{cov}}=\max_{u\in V}\mathrm{E}[\tau_{\mathrm{cov}}(u)]$$

を$(X_t)$の 全訪問時間(cover time) と呼びます.

1.2. SRWに対する有名なバウンド

この節では連結無向グラフ$G$上のSRWの到達時間と全訪問時間のバウンドに関する既存結果だけを簡単に紹介します. LSRWを考えても期待値的に値は2倍しか変わりません. $G$の頂点数を$n$, 辺数を$m$とします.

自己ループ付きの完全グラフ上のSRWの全訪問時間はクーポンコレクター問題になり, 古くから研究されており, $t_{\mathrm{cov}}=(1+o(1))n\log n$です.

$n$頂点のパスグラフは最悪な始点として端の頂点を考えればよく, この場合はもう一方の端の頂点に到達するまでに要する時間を考えれば$t_{\mathrm{hit}}=t_{\mathrm{cov}}=\Theta(n^2)$が示せます. $\Theta(n^2)$の直感としては, 無限の数直線$\mathbb{Z}$上の$T$ステップの単純ランダムウォークは$T$個の独立なRademacher確率変数 (一様に$\pm 1$のいずれかを出力する確率変数)の和として表すことができ, 中心極限定理より$T$ステップ後の座標は絶対値が大体$\Theta(\sqrt{T})$くらいになります. $n$だけずれるのに要する時間は$\Theta(\sqrt{T})=n$を$T$について解いて$T=\Theta(n^2)$を得ます.

明らかに到達時間は全訪問時間の下界になるため$t_{\mathrm{hit}}\leq t_{\mathrm{cov}}$なのですが, 実は$t_{\mathrm{cov}}=O(\log n\cdot t_{\mathrm{hit}})$が任意のグラフで成り立ちます. これはランダムウォークの文脈ではMatthew's bound [Matthew, Ann. Probab. 88]などと呼ばれる非常に有名なテクニックです.

Aleliunasら [Aleliunas, Karp, Lipton, Lovász, and Rackoff, FOCS1979]は任意の連結グラフに対し$t_{\mathrm{cov}}\leq 2m(n-1)$を証明しました. この証明は全域木に沿った軌道を考えそれぞれの辺に対するedge commute timeの総和を考えることで得られており, 1ページほどにまとめられていて面白いです.

Kahnら [Kahn, Linial, Nisan, Saks, J. Theor. Probab., 1989]は$\frac{m}{2d_{\min}}\leq t_{\mathrm{hit}} \leq t_\mathrm{cov}\leq 16mn/d_{\min}$ ($d_{\min}$は最小次数) を示しました. ここから系として, 正則グラフに対する$t_{\mathrm{cov}}=O(n^2)$というバウンドが得られます. この定理の証明はAleliunasらの証明における全域木を巧妙に構成することによって得られています. この上界は定数倍を除いてタイトで, t_{\mathrm{hit}}=\Omega(n^2)$を満たす3正則グラフが構成されています (cf. Aldous--Fill本).

Aleliunasらの上界がタイトになる例も知られており[Brightwell, Winkler, RSA1990]はロリポップグラフと呼ばれるグラフが$t_{\mathrm{hit}}=(4/27+o(1))n^3$になることを証明し, Feige [Feige, RSA1995]は任意の連結グラフで$t_{\mathrm{cov}}=(4/27)n^3+O(n^{2.5})$となることを証明しました.

hardcore lemmaとその証明

1. 導入

Impagliazzoのhardcore lemmaと呼ばれる非常に面白いメッセージ性を持つ平均計算量理論において有名な結果[Impagliazzo, FOCS95]を紹介し, その証明をします. 直感的に言えばこの結果は「任意の"まぁまぁ困難な"計算問題に対して"めちゃくちゃ難しい"インスタンスの集合$H$であって密なものが存在する」ということを主張します.

平均計算量理論については以前の記事で触りだけ紹介しています. 本記事では問題として判定問題$f\colon \{0,1\}^n\to\{0,1\}$を考え, 入力分布としてある有限集合$H$ (例えば$H=\{0,1\}^n$) 上の一様分布を考えます. また, 計算モデルとしては回路を考え, その効率性を回路サイズ (=ゲートの個数) で測ります. これを回路計算量と呼びます. アルゴリズムの時間計算量と回路計算量の関係性は非自明で以前の記事で触れていますが, 任意のアルゴリズムはその動作を回路で（少しのオーバーヘッドがあるものの）シミュレーションできるため, 「時間計算量が小さいならば回路計算量が小さい」は成り立ちます (逆は必ずしも成り立たちません).

2. Hardcore lemmaの主張

定義1 (関数のhardness).

関数$f\colon\{0,1\}^n\to\{0,1\}$は, 任意のサイズ$s$以下の回路$C$に対し$\Pr_{x\sim H}[C(x)=f(x)]\leq 1-\delta$を満たすとき, サイズ$s$に対し$H$上で$\delta$-hardであるという. ここで, $\Pr_{x\sim H}[...]$は$x\in\{0,1\}^n$が$H$上一様ランダムにサンプリングされた時の確率をとることを意味する.

$\delta$の値が大きいほど$f$は難しいと解釈できますが, 常に0を出力する回路と常に1を出力する回路を考えればどちらか一方は半分以上の入力$x$に対し$f(x)$を正しく計算できるので, $\delta>1/2$にはなりえません. 

定理1 (Impagliazzo's hardcore lemma)

任意の$\delta,\epsilon>0$に対し以下が成り立つ: $f\colon\{0,1\}^n\to\{0,1\}$がサイズ$s$に対し$\{0,1\}^n$上で$\delta$-hardであるならば、サイズ$\delta 2^n$以上のある集合$H$が存在して$f$はサイズ$O(\delta^2\epsilon^2 s)$に対し$H$上で$(1/2-\epsilon)$-hardである.

冒頭でも述べた通り, $f$が$\{0,1\}^n$上でまぁまぁ難しい (i.e., $0.01$-hard) ならば, $f$がめちゃくちゃ難しい(i.e., $0.499$-hard)ようなサイズ$0.01\times 2^n$以上のインスタンスの集合$H$が存在することになります.

3. 証明

主に二通りの証明方法 (ブースティングアルゴリズムに基づくものとNewman's min-max theoremに基づくもの) が知られています. 後者の証明はArora-Barakの本などにも紹介されているので, 本記事では前者の証明を紹介します.

3.1. Hardcore measure

まず, setの連続緩和としてmeasureという概念を定義し, hardcore setの存在性の代わりにhardcore measureの存在性を証明します. 

定義2 (measure)

関数$M\colon\{0,1\}^n\to\{0,1\}$をmeasureと呼ぶ (測度論における測度とは区別するのを注意してください). $M$のサイズを$|M|=\sum_{x\in\{0,1\}^n}M(x)$とし, 相対サイズを$\mu(M)=|M|/2^n$とする. $\mu(M)\geq \delta$であるとき, $M$は$\delta$-denseであるという.

定義3 (measureにおけるhardness)

$M$をmeasureとする. 関数$f\colon \{0,1\}^n\to\{0,1\}$は以下を満たすとき, サイズ$s$に対し$M$上で$\delta$-hardであるという: 任意のサイズ$s$以下の回路$C$に対し$\Pr_{x\sim M}[C(x)=f(x)]\leq 1-\delta$.

ここで, $\Pr_{x\sim M}[...]$は$M(x)$に比例する確率で$x$をサンプリングした時の確率を考える.

注釈: measureは集合$H\subseteq \{0,1\}^n$を"連続化したもの"と解釈できます. 例えば$M$として$H$の指示関数 ($M(x)=1$ if $x\in H$, $M(x)=0$ if $x\not\in H$) をとれば, $M$が$\delta$-dense $\iff$ $|H|\geq \delta 2^n$ となります.

まず, 定理1の代わりにそれのmeasure版となる以下の結果を証明します:

定理2 (hardcore measure)

任意の$\delta,\epsilon>0$に対し以下が成り立つ: $f$が$\{0,1\}^n$上でサイズ$s$に対し$\delta$-hardならば, ある$\delta$-denseな measure $M$が存在し, $f$は$M$上でサイズ$O(\delta^2\epsilon^2s)$に対し$(1/2-\epsilon)$-hardである.

(証明)

待遇 (i.e., 任意の$M$上でhardでないならば$\{0,1\}^n$上でhardでない) を示します. 具体的には, 任意の$\delta$-denseなmeasure $M$に対し, $f$は$M$上でサイズ$s'$に対し$(1/2-\epsilon)$-hardでないと仮定します. すなわち, 任意の$\delta$-denseな$M$に対しあるサイズ$s'$の回路$C_M$が存在して$\Pr_{x\sim M}[C(x)=f(x)]>1/2+\epsilon$とします. これらの回路$C_M$を組み合わせることによって, $\Pr_{x\sim \{0,1\}^n}[C'(x)=f(x)]>1/2+\delta$となるようなサイズ$O(\epsilon^{-2}\delta^{-2}s')$の回路$C'$を構成します.

最終的に得られる回路$C'$は$C'=\mathsf{maj}(C_1,\dots,C_t)$の形になります ($\mathsf{maj}(x_1,\dots,x_t)$は$\sum x_i\geq t/2$以上ならば$1$, そうでなければ$0$を返す関数で, この関数はサイズ$O(t)$の回路で計算できることが知られています.) 以下では$C_1,\dots,C_t$を定義します.

$\gamma:=\delta\epsilon$とします. まず, 定数measure $M_1\equiv 1$ から始めます. measure $M_i$に対し$C_i=C_{M_i}$とします (待遇の仮定からこのような回路が存在). 事象$Z$の符号付き指示関数を$[Z]$と定義し (i.e., $Z$が真なら$[Z]=1$, そうでないなら$[Z]=-1$), $R_i(x)=\sum_{j=1}^i [C_i(x)=f(x)]$とします. 次のステップにおけるmeasure $M_{i+1}$を,

$$M_{i+1}(x) = \begin{cases} 1 & \text{if }R_i(x)\leq 0,\\1-\gamma R_i(x) & \text{if }0<R_i(x)<1/\gamma,\\0 & \text{if }R_i(x)>1/\gamma\end{cases}$$

で定義します. イメージとしては, これまで得られた回路の過半数が間違えるような入力に対しては最大の重み$M_{i+1}(x)=1$を与え, これまでの回路のうち多くが正解するような入力に対しては重みを与えない ($M_{i+1}(x)=0$) ようなものです. これを$\mu(M_{i+1})\leq \delta$になるまで繰り返し, 得られた回路を$C_1,\dots,C_t$とします.

最終的な回路$C'$の性能を評価します. 入力$x$に対して$C'(x)\neq f(x)$ならば, $C_1,\dots,C_t$のうち過半数が違う値を計算してるので, $C_1(x)+\dots+C_t(x)<t/2$となり, $R_t(x)<0$となります. このとき上記の構成より$M_{t+1}(x)=1$です. したがって

$$\Pr_{x\sim \{0,1\}^n}[C'(x)\neq f(x)] \leq \Pr_{x\sim \{0,1\}^n}[M_{t+1}(x)=1] \leq \frac{1}{2^n}\sum_{x\in\{0,1\}^n}M_{t+1}(x)\leq \delta$$

となるので, $C'$は$\{0,1\}^n$上で$f$を$1-\delta$以上の確率で計算していることになります. また, 各$C_i$のサイズ$\leq s'$であることから$C'$のサイズは$O(ts')$です.

最後に$t=O(\epsilon^{-2}\delta^{-2})$を示せば証明が完了します. ポテンシャルとして

$$A_t:=\sum_{1\leq i\leq t}\sum_{x\in\{0,1\}^n}M_i(x)[C_i(x)=f(x)]$$

を考え, その上下界を求めます:

(1): $A_t\geq t2^{n+1}\gamma$.

回路$C_i$の$M_i$上の正答率の仮定から$$\sum_{x\in\{0,1\}^n}\frac{M_i(x)}{|M_i|}[C_i(x)=f(X)]=2\Pr_{x\sim M_i}[C_i(x)=f(x)]-1\geq 2\epsilon$$なので,

$$A_t=\sum_{1\leq i\leq t}|M_i|(2\Pr_{x\sim M_i}[C_i(x)=f(x)]-1)\geq t\delta 2^{n+1}\epsilon=t2^{n+1}\gamma.$$

(2): $A_t\leq 2^n(\gamma t/2+1/\gamma)$.

$x\in\{0,1\}^n$を固定し, $S:=\sum_{i=1}^tM_i(x)[C_i(x)=f(x)]$を上から抑えます. $a\in\mathbb{Z}$に対して

$$m(a) = \begin{cases} 1 & \text{if }a\leq 0,\\1-\gamma a & \text{if }0<a<1/\gamma,\\0 & \text{if }a>1/\gamma\end{cases}$$

と定めると, $M_i(x)=m(R_{i-1}(x))$と書けます. ここで, 次で定義される辺重み付き有向グラフ$G=(V,A)$を考えます: $V=\{R_i(x)\colon i=0,\dots,t\}$, $A=\{(R_i(x),R_{i+1}(x))\colon i=0,\dots,t-1\}$, 辺$(R_i(x),R_{i+1}(x))$の重みは$M_i(x)(R_i(x)-R_{i-1}(x))=m(R_{i-1}(x))(R_i(x)-R_{i-1}(x))$.

有向グラフ$G$

ここで, 辺のペア$(a,a+1)$と$(a+1,a)$を考えます. これらの辺は逆向きなので重みの符号は異なり, その大きさは$|M(a)-M(a+1)|\leq \gamma$を満たします. このような辺のペアは高々$t/2$本あり, 総和に対して高々$\gamma t/2$だけ寄与します. 次に, この辺のペアを$G$から除去すると残るグラフは右に行き続けるパスまたは左に行き続けるパスになります.

ペアを全て除去して残ったグラフ

上図のように左に行き続けるパス (つまり, 頂点$a\to a-1\to\dots \to a-\ell$を辿るパス) の辺重みは全て負なので, パス重みの総和には寄与しません. また, 右に行き続けるパス (頂点$a\to a+1\to \dots \to a+\ell$を辿るパス) の寄与は$m(a)+m(a+1)+\dots+m(a+\ell-1)$ですが, $b>1/\gamma$なる任意の$b$に対して$m(b)=0$なのでこの寄与は高々$1/\gamma$となります. 従って, パス重みの総和は高々$\gamma t/2 + 1/\gamma$です. 考えるべきパスは$2^n$個あるから, $A_t\leq 2^n(\gamma t/2+1/\gamma)$を得ます.

(1),(2)のバウンドより, 

$$t2^{n+1}\gamma \leq A_t \leq 2^n(\gamma t/2+1/\gamma)$$

となるので, $t\leq 2/(3\gamma^2) \leq (\delta\epsilon)^{-2}$を得ます.

3.2. Hardcore measure $\Rightarrow$ Hardcore set

3.1の議論により, hardcore measure $M$の存在性が言えました. randomized roundingとprobabilistic methodを組み合わせてhardcore set $H$の存在性を証明します. 整理すると, $M$は$\delta$-denseであり, かつ任意の小さい回路$C$に対して

$$\Pr_{x\sim M}[C(x)=f(x)]\leq \frac{1}{2}+\epsilon$$

を満たします.　この節では$|H|\geq \delta/2\cdot 2^n$であり, かつ任意の小さい回路$C$に対して

$$\Pr_{x\sim H}[C(x)=f(x)]\leq \frac{1}{2}+\epsilon$$

を満たす集合$H\subseteq\{0,1\}^n$の存在性を証明します.

Hardcore measure $M$ を考えます. 各$x\in\{0,1\}^n$に対して, 独立に確率$M(x)/|M|$で$x$を含むようなランダム集合を$H$とします. まず, $|H|$の期待値は$|M|$なので, Hoeffding boundより確率$1-2^{-\Omega(2^{n})}$で$|H|\geq 0.9|M|$を満たします (つまり$H$は高確率で密な集合). 次に$\Pr_{x\sim H}[C(x)=f(x)]=\frac{|\{x\in H\colon C(x)=f(x)\}|}{|H|}$を考えます ($H$がランダム集合なので, この値は$H$のランダムネスに起因する確率変数となる). 先ほど分母の下界を求めたので, 次は分子の上界を求めます. まず, 回路$C$を固定して$A_C=\{x\colon C(x)=f(x)\}$を回路$C$が正解する入力の集合とします. すると分子は$|H\cap A_C|=\sum_{x\in A_C}\mathbf{1}_{x\in H}$となり, これは($H$の構成法により)独立な確率変数の和となります. さらにその期待値は$\mathrm{E}[|H\cap A_C|]=|M|\sum_{x\in A_C}\frac{M(x)}{|M|}=|M|\Pr_{x\sim M}[C(x)=f(x)]\leq (1/2+\epsilon)|M|$となるので, 非常に高い確率 (確率$1-2^{-\Omega(2^n)}$) で$|H\cap A_C|\leq 1.001(1/2+\epsilon)|M|\leq (1/2+\epsilon')|M|$となります. 従って, 回路$C$に関するunion bound ($C$が小さいサイズ, 例えば多項式サイズならば, $C$の総数は$2^{\mathrm{poly}(n)}$以下になる) により, 

$$\Pr_H[\forall C,\Pr_{x\sim H}[C(x)=f(x)]\leq \frac{1}{2}+\epsilon'] \geq 1-2^{-\Omega(2^n)}\cdot 2^{\mathrm{poly(n)}}>0.$$

よって正の確率で$H$はhardcore setの条件を満たすので, 所望の$H$の存在性が証明できました (証明終).

4. 教師有り学習のブースティングとの関係

教師有り学習の文脈では, $f\colon\{0,1\}^n\to\{0,1\}$を二値分類タスクとみなし, $C_M$を$M$上で動作する弱い分類器とみなすことにより, $C'$の構成をブースティングの一種と解釈することができます (ただし上記の証明ではあくまでも$C'$の存在性にすぎず, その構成アルゴリズムを与えたものではないため, ブースティングアルゴリズムそのものは与えていません). 上記の証明は$M$を上手く更新しながら対応する分類器をとってきてそれらのmajorityをとるものになっており, 今回の証明以外にも$M$の更新方法は様々知られており, 例えばAdaBoostのようなmultiplicative updateに基づくもの (+Bregman projection)も知られており, $t=O(\epsilon^{-2}\log\delta^{-1})$まで改善されています (Barak et al. SODA09).

全点対最短経路問題の(+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]などにも応用される汎用的なアイデアなので, 非常に有用だと思います.

確率変数の中央値と期待値の差

確率集中不等式はしばし実数値確率変数$Z$に対して$\Pr[|Z-\mathbf{E}Z|\geq t]$の上界を与えますが, Talagrandの不等式に代表されるようなisoperimetric inequalityを考えると期待値$\mathbf{E}Z$の代わりに中央値$\mathbf{M}Z$を考えることがよくあります. ここで確率変数$Z$の中央値(median)とは, $\Pr[Z\geq \mathbf{M}Z]\geq 1/2$ かつ $\Pr[Z\leq\mathbf{M}Z]\leq 1/2$を満たす実数$\mathbf{M}Z$のことです. Isoperimetric inequalityを使うと$\Pr[|Z-\mathbf{M}Z|\geq t]$の上界を与えることができます.

平均値と中央値は一般に異なりますが, 例えば$Z$の分散が小さいとバラつきが小さいので $\mathbf{E}Z$と$\mathbf{M}Z$の間はそれほどギャップがないように思われます. この直感は実際正しく, 以下の事実が成り立ちます.

定理

$\mathbf{E}Z^2<\infty$を満たす確率変数$Z$に対し,

$$|\mathbf{M}Z-\mathbf{E}Z|\leq \sqrt{\mathrm{Var}Z}$$

本記事ではこの事実を証明していきます.

補題 (Cantelli's inequality)

任意の$t>0$に対し,

$$\Pr\left[\frac{Z-\mathbf{E}Z}{\sqrt{\mathrm{Var}Z}} \geq t\right] \leq \frac{1}{1+t^2}.$$

(補題の証明)

平行移動とスケーリングによって一般性を失わず$\mathbf{E}Z=0$, $\mathrm{Var}Z=1$と仮定してよい. すなわち, 平均0で分散1の確率変数$Z$に対して$\Pr[Z\geq t]\leq 1/(1+t^2)$を示せば良い.

マルコフの不等式より, 任意の$\lambda>0$に対して

$$\Pr[Z\geq t] \leq \Pr[(Z+\lambda)^2\geq (t+\lambda)^2] \leq \frac{\mathbf{E}Z^2+\lambda^2}{(t+\lambda)^2}=\frac{1+\lambda^2}{(t+\lambda)^2}.$$

関数$\lambda\mapsto \frac{1+\lambda^2}{(t+\lambda)^2}$は$\lambda=1/t$の時に最小化され, このとき$1/(1+t^2)$をとる (証明終).

(定理の証明)

補題において$t=1$を代入すると, $\Pr[Z\geq \mathbf{E}Z+\sqrt{\mathrm{Var}Z}] \leq 1/2$を得る. 中央値の定義より

$$\Pr[Z\geq \mathbf{E}Z+\sqrt{\mathrm{Var}Z}] \leq 1/2 \leq \Pr[Z\geq \mathrm{M}Z].$$

よって$\mathrm{M}Z\geq\mathbf{E}Z+\sqrt{\mathrm{Var}Z}$を得る.

一方, 新たな確率変数$-Z$に対して$t=1$として改めて補題を適用すると, $-1$倍しても分散は変わらないため, $\Pr[Z\leq \mathbf{E}Z-\sqrt{\mathrm{Var}Z}]\leq 1/2$となり, 同様にして$\mathrm{M}Z\leq\mathbf{E}Z-\sqrt{\mathrm{Var}Z}$を得る (証明終).

モーメント母関数のlogについて

実数値をとる確率変数 $X$ と適当な値$t\in\mathbb{R}$に対して, 確率$\Pr[X\geq t]$を抑えるときにモーメント母関数 (MGF: Moment Generating Function) を考えたくなります. 確率変数$X$のモーメント母関数は$M_X(c)=\mathrm{E}[e^{cX}]$で定義されます. これを使うと, 任意の$\lambda\geq 0$に対して

$\Pr[X\geq t] = \Pr[\lambda X \geq \lambda t] = \Pr[e^{\lambda X} \geq e^{\lambda t}] \leq e^{-\lambda t} M_X(\lambda)$

が成り立つ (最後の不等式ではMarkovの不等式を適用した) わけです. 一般に浸透された言い回しではありませんが, モーメント母関数のlogをとったものをここではlogMGFと呼ぶこととします. すなわち, 確率変数$X$のlogMGFは

$$\Psi_X(\lambda) = \log M_X(\lambda) = \log \mathrm{E}[e^{\lambda X}]$$

です. これを使うと任意の$\lambda\geq 0$に対して

$$\Pr[X\geq t] \leq \exp(-\lambda t + \Psi_X(\lambda))$$

となるので, 最適な$\lambda$をチューニングしたくなります. ここで$\Psi_X^*(t) = \sup_{\lambda\geq 0}\{\lambda t - \Psi_X(\lambda)\}$ を考えてやると,

$$\Pr[X\geq t] \leq \exp(-\Psi^*_X(t))$$

を得ますが, $\Psi_X^*(t)$をよくみると関数$\Psi_X(\lambda)$のルジャンドル変換っぽいことが分かります ($\sup$で動く$\lambda$の範囲が$\lambda\geq 0$でなく$\lambda\in\mathbb{R}$ならばルジャンドル変換と一致します. なお, $t\geq \mathrm{E}[X]$ならば, $\Psi_X^*(t)=\sum_{\lambda\in\mathbb{R}}\{\lambda t - \Psi_X(\lambda)\}$ となり, ルジャンドル変換と一致します).

$\Psi_X(\lambda)$の有名な性質として凸性が挙げられます. つまり, $S=\{\lambda>0\colon \Psi_X(\lambda)<\infty\}$とすると$\Psi_X(\lambda)$は$S$上で凸関数になっています. これは, $x,y\in S$を$0<x<y$を満たす実数として$p\in(0,1)$を任意に選んで$\Psi_X(px+(1-p)y)$をHölderの不等式を使って上から抑えることによって証明することができます.

直径の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のLinear Decision Tree Complexity

Linear Decision Treeと呼ばれる計算モデルについて解説します.

1. 決定木 (Decision Tree)

「ソートアルゴリズムは$\Omega(n\log n)$時間かかる」という話をよく聞くと思いますが, これは「ソートアルゴリズムは2要素の比較を$\Omega(n\log n)$回しなければならない」ことを意味します. 具体的には各ノードに「$a_i\leq a_j$?」と書かれた二分木があり, YESならば右の子ノード, NOならば左の子ノードに辿っていき葉に至ったらそこで終了しそこに書かれているラベル (例えば「$a_3\leq a_1 \leq a_2 \leq a_4$」のような要素の順位が書かれている)を出力するという計算モデルを考えています. このように二分木にYES/NOの質問が書かれており, その答えに対応する子ノードを辿っていき葉に至ったら終了するという計算モデルを決定木 (decision tree) といいます.

例えば上記の決定木では, どんな入力が与えられても質問に3回答えれば葉に至ります. このように. 決定木の効率性は主にその深さによって測ります (もちろん全体のノード数が小さい方が良いですが). ここで, 深さ$D$の決定木の葉の個数は$2^D$くらいです. つまり深さ$D$の決定木は高々$2^D$通りの出力しか区別することができないので, 出力しうるラベルの種類が$N$であるような問題を考えるとその問題を解く決定木の深さは$\lceil\log_2 N\rceil$必要であるという情報論的な下界が得られます. 例えば$n$要素のソートでは$n!= 2^{\Theta(n\log_2 n)}$通りのラベルを区別しなければならないため深さは$D=\Omega(n\log n)$となります. ある問題に対してその問題を解く決定木の深さの最小値を決定木複雑性 (decision tree complexity) といいます. ソートの決定木複雑性は$\Omega(n\log n)$ですが, 実際にマージソートなどは高々$O(n\log n)$回の数値比較でソートをするので, この下界は定数倍を除いてタイトとなります.

決定木の深さはちょうどクエリに回答する回数に対応するため, 計算時間に対応していると解釈できるわけです.

2. 線形決定木

一般に決定木のノードに書かれるラベルはどんな質問でも構いません. 極端な例で言うとノードのラベルにSATのインスタンスが書いてあってそれがYESインスタンスならば左の子に, そうでなければ右の子に辿るようなものを考えてもよいです. しかしながらそれでもなおソートの決定木複雑性は$\Omega(n\log n)$になります.

決定木のノードラベルの質問を線形不等式「$c_1a_1+c_2a_2+\dots+c_na_n\geq 0?$」の形に制限したものを線形決定木 (linear decision tree) といいます. ここで$a_1,\dots,a_n$は入力, $c_1,\dots,c_n$は適当な定数とします. また, 全てのノードのラベルにおける係数に含まれる非ゼロ要素が高々$k$個であるとき, 特に$k$-線形決定木 ($k$-LDT) といいます. 例えばソートの二要素の比較は$a_i-a_j\geq 0$と書き換えることができるので2-LDTとなります.

幾何的な解釈としては, 入力全体の空間を超平面で分割 (線形不等式クエリに対するYES/NOの応答に対応)していき「ラベル1を出力すべき入力の全体」「ラベル2を出力すべき入力の全体」...を区別していくことになります.

3. 決定木とアルゴリズムの違い: uniform vs nonuniform

一つの決定木は固定サイズの問題しか解くことができません. 従って決定木を使って問題を解く場合は決定木の族$(\mathcal{T}_n)_{n\in\mathbb{N}}$を考えます. 問題Pのサイズ$n$の任意の入力$x$に対し, $x$の答えを$\mathcal{T}_n$を使って計算できる時に$(\mathcal{T}_n)$は問題Pを解くといいます. つまり, 与えられる問題のサイズに応じて異なる決定木を使って良いということになります. 

一方でアルゴリズムを使って問題を解く時は入力サイズがどのような場合であっても必ず同一のアルゴリズムを使わなければなりません. このように, 決定木とアルゴリズム(正確にはチューリングマシン)の間には問題サイズによって異なるものを使って良いかどうかという重要な違いがあります. 決定木のように問題サイズごとに異なるものを使う計算モデルをnonuniform, チューリングマシンのようにどのような問題サイズであっても同じものを使う計算モデルをuniformといいます.

UniformとNonuniformな計算モデルのギャップは未解決な部分が多いですが, 少なくともギャップが存在することが分かっています. 具体的にはnonuniformな計算モデルを使うと停止性判定問題を解くことができます (以前の記事参照). ちなみにnonuniformな計算モデルの計算量として最もよく研究されているのは(私の主観ですが)回路計算量だと思います.

今回の記事では線形計算機とアルゴリズムのギャップを, 3SUMという問題に焦点を当てて見ていきます.

4. 3SUMのLDT複雑性

3SUMという問題は

$n$個の数字$A=\{a_1,\dots,a_n\}\in\mathbb{R}$が与えられた時に$a_i+a_j+a_k=0$となる$i,j,k$が存在するか?

を判定する問題です. $a_i+a_j$を格納したサイズ$n^2$の配列Xを用意してこれをソートした後, 各$a_k$に対して$X$が$-a_k$を含むかどうかを二分探索でチェックすれば$O(n^2\log n)$時間で解くことができます. また, 前回の記事で紹介したように$O(n^2)$時間アルゴリズムも存在します. しかしながら任意の定数$\epsilon>0$に対し$O(n^{2-\epsilon})$時間アルゴリズムが存在するかどうかは未解決で, 精緻な計算量 (fine-grained complexity)の分野では$n^{2-o(1)}$時間かかると予想されています. とりあえず詳しくない人は「3SUMを$n^{1.99999}$時間で解くアルゴリズムは知られていない」と思っておけば大丈夫です.

本記事では3SUMのLDT複雑性に関する以下の驚くべき結果の証明(概略)を説明します.

定理1 ([1])

3SUMを解く深さ$O(n^{1.5}\log n)$の4-LDTが存在する.

実は3SUMを解く任意の3-LDTは深さが$\Omega(n^2)$あることが知られていたため[2], 4-LDTにすると一気に深さが減らせるというのはLDT界隈の人にとっては驚くべき事実だったことは想像に難くありません. 何よりも「3SUMは$O(n^{1.5}\log n)$回の要素の比較で解ける」ということになるわけです.

しかし定理1の証明はかなり単純かつ非常に賢いアイデアに基づいています. 定理1の4-LDTは以下のステップからなります.

ステップ1. 入力$A=\{a_1,\dots,a_n\}$を昇順にソートする.

ステップ2. $A$の各要素を連続した$\sqrt{n}$個ずつに分割してそれらを$A_1,\dots,A_{\sqrt{n}}$とする.

ステップ3. $L:=\cup_{i\in[\sqrt{n}]}\{x-y\colon x,y\in A_i\}$をソートする.

ステップ4. 以下のしゃくとり法っぽい計算を行う. ただし文中に出てくる$A_{st}$は$A_{st}:= \{x+y\colon x\in A_s,y\in A_t\}$とする.

for a in A:

    $s \leftarrow 0$

    $t \leftarrow \sqrt{n}$

    while $s\leq t$:

        if $-a\in A_{st}$ then:                            (1)

            return True

        elif $-a > \max(A_s)+\min(A_t)$:            (2)

            $s\leftarrow s+1$

        else:                                                     (3)

            $t\leftarrow t-1$

return False

正当性 (イメージ).

(1)において$-a\in A_{st}$となるならば, ある$b\in A_s,c\in A_t$に対して$-a=b+c$, すなわち$a+b+c=0$となるので, 確かにYESインスタンスとなります. 逆に, 入力がYESインスタンスであるとします. このときwhile文以下のループ処理は以下のような幾何的な解釈ができます:

入力$A=\{a_1,\dots,a_n\}$に対し, 点集合$P=\{(a_i,a_j):a_i,a_j\in A\}$を$xy$平面にプロットする. 各$a\in A$に対して直線$x+y=-a$上に乗っているような$P$の点の有無を判定していきます. ここで, $\mathbb{R}^2$を図のよう長方形で区切っていき, 各グリッドが$P$の点を$n$個ずつ含むようにします (実際には図のようなグリッドにはならないとは思いますが, 説明の簡単のためグリッドで区切っています). すると, while文の最初の時点では$s=0,t=\sqrt{n}$で点$(b,c)\in A_s\times A_t$に着目するわけですが, これは最も左上にある長方形に着目することに対応します. (1)では着目している短形内に直線上の点があるかどうかを判定しています. (2)と(3)では短形内の右下にある点に着目し, この点と直線の位置関係によって右隣か下の短形どちらに進むべきかを判断しています.

3SUMのLDTのしゃくとり法のイメージ (右の吹き出しの$\sqrt{n}$はタイポで正しくは$n$です)

深さの解析.

- ステップ1はマージソートをそのまま実装すれば良く, 深さ$O(n\log n)$の2-LDTで実装できます.

- ステップ2は分割しただけです.

- ステップ3はステップ1と同様にソートするだけです. $L$の要素数は$n^{1.5}$なので, 深さは$O(n^{1.5}\log n)$となります. また, $L$の各要素は$a_i-a_j$の形で書けるため, 4-LDTとなります.

- ステップ4では, while文が$s\leq t$である限り続けますが1回の反復で$t-s$は1ずつ減少するため, 高々$\sqrt{n}$回の反復で終了します. 各反復では(1)の処理がボトルネックになります. もし各$A_{st}$があらかじめソートしてあるならば二分探索ができるので(1)は$O(\log n)$時間で行うことができます. しかしながら全ての$A_{st}$をソートしようとすると, $0\leq s\leq t\leq \sqrt{n}$に対し各$A_{st}$は要素数$n$を持つため, $O(n^2\log n)$時間または同程度の深さのLDTが必要になってしまいます.

ここで, ステップ3で$L$をソートしているという事実に着目しましょう. $L$の全ての要素の大小を比較しているため, 各$i,j$と$a_i,a'_i\in A_i$, $a_j,a'_j\in A_j$に対し

$$a_i-a'_i \leq a_j-a'_j$$

またはその逆方向の不等式どちらが成り立つかが分かっています. これを移項すると

$$a_i+a'_j \leq a'_i+a_j$$

となりますが, この両辺はともに$A_{ij}$の要素です. これら4つの要素は任意にとって良かったので, 

$L$の大小関係の情報が確定すれば各$A_{ij}$の全ての2要素の大小関係の情報は必ず確定します.

決定木でステップ1,2,3を経たノードにいるならば, $L$の要素の大小の順位の情報が確定しているということになるわけなので, 必然的に全ての$A_{ij}$の大小関係の情報が確定しており, 実質的に全ての$A_{ij}$がソートされているとみなすことができます. そこでこの情報を使ってステップ4で二分探索を行うことができるので, 全体の深さは$O(n^{1.5}\log n)$となります.

5. 何が肝なのか？

4章で紹介したLDTでは, ステップ3で計算した$L$の大小関係を使ってステップ4(1)で行う二分探索のために必要な$A_{st}$の要素の大小関係の情報を得られる, という部分が本質的に最も効いています. あくまでも「得られる」というわけであって「簡単に計算できる」というわけではないのがミソです。例えばソートアルゴリズムをそのままLDTに翻訳した場合は「$a_i-a_j\geq 0?$」というノードに対してインデックス$i,j$を簡単に計算して求められていたのですが, 今回考えたLDTではこのようなインデックスは存在性が保証されているだけであって簡単に計算できるとは限らないという点がポイントです。このように存在性と計算可能性のギャップがチューリングマシンとnonuniformモデルのギャップに繋がってくるのは面白い観点だと思いました。

参考文献

[1] Allan Grønlund and Seth Pettie, Threesomes, Degenerates, and Love Triangles, JACM, 65(4), pp.1--25, 2018 (preliminary version is in FOCS14).

[2] Jeff Erickson, Bounds for linear satisfiability problems, Chicago J. Theor. Comput. Sci.,1999(8), 1999 (preliminary version is in SODA95).