Processing math: 0%

2021年12月1日水曜日

全点対最短経路問題の(+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をそれぞれSTで表します. 即ち, 

SS\cap N(v)\neq \emptyset (\forall v\in H)

TT\cap N(w)\neq\emptyset (\forall w\in M)

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

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

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に対して以下の条件を満たすとき, 全域部分グラフSG(\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]などにも応用される汎用的なアイデアなので, 非常に有用だと思います.

2021年11月20日土曜日

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

確率集中不等式はしばし実数値確率変数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}を得る (証明終).



2021年11月19日金曜日

モーメント母関数の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 S0<x<yを満たす実数としてp\in(0,1)を任意に選んで\Psi_X(px+(1-p)y)をHölderの不等式を使って上から抑えることによって証明することができます.

2021年10月11日月曜日

直径の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}RS_vの頂点を一つ以上含みます. 従って, union boundにより確率1-1/nRは 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\ellH内のいずれかの頂点にたどり着くことができます. ここで\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_uuから近い順に頂点を取っていて, 距離\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.


2021年9月17日金曜日

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,ja_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).

2021年8月27日金曜日

小記事: 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)倍の改善はなされています.


2021年7月27日火曜日

講義を完走した感想

人生で初めて講義(教える側)をしたのですが、とても濃い経験になったので備忘録として残そうかと思います。大学で講義を「する」という経験は一般的に見ればなかなかない経験であまり想像がつかないと思いますので、どういう感じなのかを説明します。講義資料の準備で意識したことや悩んだことなども綴っていくので、これから講義をする機会のある方には参考になる部分もあるんじゃないかと思います。既に幾つも講義を受け持っているプロフェッショナルな方々には自身の経験を思い出して温かい目で見守っていただけたら幸いです。コロナ禍という特殊な環境ではあったものの初の対面講義で奮闘した一人の助教の記録です。


0. 背景

自分が担当した講義はプログラミングの発展的な内容を扱うもので、学部2年生向けの講義です。まだ学部2年生の前期なのでまだあまり専門分野を学んでいない状態です。とはいえプログラミングの基礎的な部分はある程度学んでいてPythonは少しは分かるという状態です。また、必修科目というわけでもないので、受講者はある程度プログラミングに興味を持っているという想定です。また、ORなどに興味を持ちそうな理系向けの講義になっています。


1. 講義内容と形式の決定

自分が最初に苦労したのはこの部分です。特に何をすべきかが厳密に定まっておらず、内容はほとんどを任されていましたが、昨年度までの資料を見る限りでは基礎的なアルゴリズムをやっていたためその方向性は守ろうと思いました。他の先生に話を聞くと「ダイクストラ法はやってほしい」みたいな声も聞いたので、この辺りはちゃんと取り込もうと思いました(もちろんダイクストラは元々やるつもりでしたが)。昨年度までは本当に幅広く(平面幾何、乱択アルゴリズム、組合せ最適化など)さまざまなアルゴリズムを軽く紹介するという感じだったので、組合せ最適化の内容を主に扱い、最後に自分の趣味(ランダムグラフ、ランダムウォークなど)を少し取り入れた感じにしようと考えていました。

私の所属する大学はクォーター制度を導入していたため、講義は全7回です (6月頭〜7月末)。週1回講義と演習の時間がそれぞれ100分ずつあります(なので計200分あります!)。また、2020年度はコロナ禍の影響でオンライン講義になっていたようですが、私が赴任した2021年度は演習ということもあり対面授業になるとのことでした。

私がD3だった2020年度はコロナ禍が始まった年ということもあり仕方のないことですが対面講義と比べるとオンライン講義はあまり楽しいものではありませんでした。なので対面で講義できるというのは嬉しく思いました。特に受講生の顔で反応を見ながらインタラクティブに教えたりするのは教育の楽しい部分だと個人的に思います。

とはいえ緊急事態宣言が発出されている状況なのでいつオンライン講義に変更になっても対応できるように、昨年度と同様に google colab (jupyter notebookをブラウザ上で動かすgoogleのサービス) のドキュメントを講義資料とすることにしました。特にgoogle colabだと学生のPCにpythonを動かすためのソフトを入れる必要もなく, googleアカウントさえあれば誰でも動かせるのもあって導入が非常に簡単です (今の時代だとほとんどの人がYouTubeを観たりするのでgoogleアカウントは持っています). ipynbファイルを講義資料として配布し、受講生は自分のPCを持参してもらえれば普通に講義を行うことができます。

2. 講義資料の準備 (第1回〜第3回)

4月に着任して早々に講義をどうするかが頭をよぎったので、かなり前もって準備をすることにしました。今思えばこの選択はかなり良かったです(何事も余裕を持って準備しておくことは非常に大切です)。

講義資料を作る前に自分の中ではグラフアルゴリズムをメインに据えるということを意識しました。グラフ理論やグラフアルゴリズムは自分が非常に得意とする分野であり、何度も実装したことがあるからです。そこで、当初は次の二つの事項を講義で扱おうと考えました:

1. Pythonを用いて基礎的なアルゴリズムを理解

  キーワード: 再帰, DFS, BFS, 貪欲法, 最短経路問題 (Dijkstra, Warshall—Floyd), DP.

2. 応用を知る: ネットワーク解析

  キーワード: ランダムウォーク, ランダムグラフ (BAモデル),  ネットワーク中心性 (PageRank, 媒介中心性) など.

ひとまず第1回の講義はPythonの復習を行い、その後は基礎的なアルゴリズムをひたすら紹介していき、最後の第7回講義で自分の趣味を交えた講義を行うという計画を立てました。

この方向で講義を組み立てて実際に講義資料を作り始めたのが4月の中旬 (講義開始2ヶ月前) でした。初回の講義資料が完成したので見返しました。初回はPythonの文法の復習だけやって「1から1000までの7の倍数の和を求めるコード」とかを1時間やるものだったので、純粋に「これはツマランな」と心配になりました。そこで少しでも退屈を紛らわそうと再帰関数の紹介を行いました。再帰を使ってフィボナッチ数列や階乗の計算のプログラムを書き、最後に非自明な例として分割数の漸化式を示しつつ分割数を計算するコードも紹介しました。

第1回の課題は「与えられた自然数が素数かどうか判定せよ」みたいな問題を出していましたが、なんとなく味気ないのでコラッツ予想を絡めた問題も出題してみました。こんな感じです。


この調子で第2回講義資料も取り掛かり、計算量とグラフの定義を行いました。グラフの隣接行列や隣接リストを紹介し、「グラフの次数列を計算せよ」「三角形(長さ3の閉路)の個数を求める関数を書け」といった課題を出していました。

課題を考える時に意識したのはデキる人でも退屈しない内容にしたいということです。講義そのものは流石に集中が散ったりするのでしょうがない部分もあるかもしれませんが、せめて競技プログラミングのように、問題を解く部分に楽しみを見出して楽しんでもらいたいと思っていました。第2回講義以降は毎週の課題の中に一つだけ学部2年生向けの講義にしてはありえないくらレベルの高い難問を設けることにしました (難問はチャレンジ課題枠としてあり、解けなくても成績に影響は出ないことを明言した上で出題しました)。

ちなみに第2回講義で出題した難問は以下の問題です:

要するに長さ4の閉路の判定をO(n^2)時間で行えという問題で、この記事にも紹介したように鳩の巣原理を使うと計算量が抑えられるという問題です。これを計算量の概念とグラフの定義を学んだばかりの学部2年生に出すのは正気の沙汰とは思えません。こんなん知らない状態から解けるやつおらんやん...

第1回と第2回の講義資料と課題の準備が一通り終わったので、他の先生や学生TAに難易度や分量を含めて確認してもらい、大丈夫そうとのことだったのでここでようやく一息つきました。ここまでの講義の内容は必須事項だったのであまり迷いはありませんでしたが、以降の構成に頭を悩ませることになりました。

ひとまず考えた計画としては残りの5回で

・可視化ライブラリ(matplotlibとnetworkx)の紹介、幅優先探索、深さ優先探索
・貪欲法、クラスカル法
・動的計画法
・Dijkstra、Warshall--Floyd、Bellman--Ford
・フロー (Ford--Fulkerson) 

をやるか、または最後のフローはやらずに代わりに自分の趣味(ランダムグラフなど)を紹介するか、どちらかにしようと決めました。

いずれにしてもグラフアルゴリズムをやることは決まっていたので、せっかくなのでアルゴリズムの可視化を行うことを思いつきました。ひとまずgoogle colab上でのアニメーションの方法を学び、試しにDFSやBFSのアニメーションを書いて遊びつつ講義資料の準備を進めました。

ちなみにこんな感じのアニメーションを作りました。受講生はgoogle colabで開いている講義資料にあるPythonコードを実行するとこのようなアニメーションが出てくるというようになっています。アルゴリズム自体を可視化しているのでグラフは自由に変えられます。また、ボタンを押した時に一コマだけ進ませるというようなことも可能です。


ちなみにこの動画のソースコードは本記事の5章で公開しています。

この調子でクラスカル法やダイクストラ法もアニメーションプログラムを作っていこうと決めていくうちに、とうとう初回講義の日がやってきました。

3. 対面講義初回


結論からいうと初回講義は機器の不調なども重なって非常に大変でした。そもそも自分の担当は講義100分+演習100分だったのですが、実は当時の自分は勘違いをしていて講義100分の中に演習が組み込まれているものと思っていました(なので実際の講義は想定の倍あったということになります)。

元々70分程度の想定の下で準備していた講義を緊張もあって少し早口で進めていき、めちゃ早く終わっちゃうと内心焦っていたら、突然PCの調子が悪くなり学生TAのPCを急遽借りて講義を再開して案の定めちゃくちゃ早く講義部分が終わって演習時間がものすごく長くなるなど、思い返しても大変な1日になってしまいました。結局初回は課題が簡単だったのでほとんどの学生がすぐに帰る感じになっていました。

今にして思えばもっとやりようはあっただろうとも思えますが、2ヶ月弱前から準備していたしやれることはやったのであまり気にせずとりあえず講義を進めていくことになります。

第2回講義は使用PCも変えたので第1回と比べるとかなりスムーズに終えられました。第2回課題には第1回とは違い難問も用意していましたが、これは予想以上に反応が良かったです。計算量を気にしてプログラムを書くということをあまり経験してこなかった学生がほとんどだと思うので四角形の高速な判定にどこまで取り組んでくれるか心配したのですが、10人弱の学生が18時くらいまで残って議論しながら取り組んでくれました。これには教員側からしてもとても嬉しかったです。


4. 講義の改善点の洗い出し


まだ2回ではあるもののここまで講義をやってきていくつか変更点が見えました。まず、そもそも講義時間を勘違いしたまま準備していたので単純問題として内容は倍にできる。どうするか?

また、学生の質問を受け付けていくうちに色々改善した方がいいかもなぁというものも多く出てきました。例えば

  • google colabのインデント幅の設定が人によって違う
  • 自身が書いたコードがちゃんと正しく動いてるかを確認するのが面倒(学生が自分でテストケースを用意して確認しなければならないが、グラフアルゴリズムだと自分で隣接リストなどを用意しなければならない)
  • printfデバッグみたいなことを教えた方がいいのか?
などです。printfデバッグは質問してきた学生に「こうしてprintをするとどこで詰まったか分かるよ」と直接教えることにしました。

インデントなど、全体に連絡すべき事項はちゃんと講義冒頭にスライドを使って説明するようにしました。

課題プログラムの自己チェックを可能にするために、Checkerというものを導入しました。これはつまるところ競技プログラミングのサンプルケースとその答えを事前にこちらで用意しておき、学生が書いたコードがこれらのサンプルケースで正しい答えを出すかどうかをチェックするプログラムを講義資料に入れておいたのです。こんな感じです。





*)完全に問題がAtCoder感があるかもしれませんが、こういう文章題は少なくて、普通の問題はグラフ理論の用語を使った問題文になっています。

また、課題の採点を楽にするためにテストケースを用意しました。基本的にはテストケースで正解していて大丈夫そうなコードを書いていたらその課題を満点にするという感じにしました。これのおかげで採点は劇的に楽になりました。

もちろん、サンプルケース, Checker, テストケースによる自動採点の準備はかなり大変でこれを講義資料の準備と並行して行うのは非常に時間がかかり大変でした。しかし、特にCheckerの導入は受講生にはかなり高評価だったようで、さらに「Checkerを走らせたらエラーが出た/Wrong answerだった」という類の質問がきてくれたことは良かったです(Checker導入以前はあまり質問をしてくれる学生がいなくて寂しかったです)。やはり、自分の書いたプログラムがちゃんと動いているかの確認が簡単にできるということは準備の大変さに目を瞑ればメリットしかありません。


5. 以降の講義資料の準備と講義


ここまで来ると講義自体も少しずつ慣れてきて、講義内容も固まってくるので特に悩むこともなくなってきます。以前立てた計画通り

・貪欲法、クラスカル法
・Dijkstra、Warshall--Floyd、Bellman--Ford
・動的計画法
・フロー (Ford--Fulkerson) 

の順番で講義資料を準備していきました。結局最終回はフローと最大二部マッチングをやることにしました。

ちなみに講義資料のために作ったアルゴリズムのアニメーションは

https://colab.research.google.com/drive/1ed0fAC6rUMcScJYxv0Kc0Ox-wDBf2FSR

で閲覧することができます(ダウンロードすれば誰でも編集とかができるので、講義などで使ってみたいという方ご自由にどうぞ。ソースコードの下から二行目のanim.save(...)の部分をコメントアウトすると動画をmp4ファイルで保存できます)。

ちなみに細かい講義内容は

第4回:

  • 貪欲法 (コイン支払い問題, 区間スケジューリング問題)
  • 最小全域木問題
  • クラスカル法


第5回

  • 最短経路問題
  • ダイクストラ
  • ベルマンフォード
  • ワーシャルフロイド


第6回

  • 動的計画法
  • 最長共通部分列問題
  • ナップザック問題
  • 重み付き区間スケジューリング問題
  • Held—Karp のTSP


第7回

  • フロー
    • 残余グラフ
    • Ford—Fulkerson
    • Edmonds—Karp
  • 二部マッチング
    • 割り当て問題

という感じになりました。難問枠としてmeet-in-the-middleなどを出題しました。

反省点は幾つもあります。
  • 第5回以外は意外と早く講義が終わってしまったというのがあります。アルゴリズムを広く浅く紹介して課題で問題を解きながら理解を深めてもらうというスタンスだったので、内容としてはちょうど良いのかもしれませんが講義自体は説明は割と早めに終わってしまいました。この辺りは説明をもっとゆっくりやればちょうど良くなりそうな気配があるので、経験を積むと感覚が掴めてくるのかもしれません。
  • 初回の講義についてはPythonの復習にたくさん時間を割くよりも、エラトステネスの篩、ユークリッドの互助法、二分探索などの軽いアルゴリズムの紹介などが出来た気がしています。復習に時間を割くよりもコードで出てくる度に「リストはこうやって要素を追加するよ」と言えば良さそうです。
  • 課題の提出についても注意喚起が足りなかった部分があります。課題を解いたのに保存で失敗して白紙の状態で提出してしまったという学生もいました。こういった不運な事故を防ぐ努力(注意喚起)はもうちょっとすべきだったと思います(流石に初めての講義でそこまで気が回らなかったのもしょうがない気がしますが。。。)

このように、色々と改善点はまだまだあるものの、初めてやったにしては上出来なんじゃないかと開き直ることにしています。

なお、講義資料の準備段階では他の先生の講義動画などが非常に参考になりました。説明の順序やその問題を考える動機などが自然に導入されており、流石だなぁと思いました(もはやこの動画が講義資料でいいじゃんと思っちゃうくらい)

6. 感想


元々自分は人に教えるのが好きだったので、とても楽しく完走することができました(もちろん準備はかなり大変でしたが。。。)来年度以降はゼロからアニメーションコードを開発したりまではしなくて良いのである程度苦労は減ることを期待していますが、今年度できなかった内容も取り入れたいし、実は少し楽しみにしています。例えば乱択アルゴリズムを取り入れて難問枠としてcolor codingを入れたりとかを妄想しています笑

一方でコロナ禍での対面授業ということもあり、感染拡大の防止にものすごく気をつけました。対面講義をするというのもあって、4月の緊急事態宣言以降は今もなお一度も人と一緒に外食をしておらず、そもそも食事はほとんど全てをテイクアウトか自炊で済ましています(講義が終わった今後も同様の生活を続けていくつもりです)

受講生自身からも対面授業が嬉しいというコメントがきている一方で彼らも感染防止にはかなり気を遣ってくれていて、もちろん完全に会話をゼロにするには流石に無理ですが、それでも常にマスクをつけて暑い中換気のために窓を開けるのに協力してくれたり、入室時には必ず手の消毒をしてくれていました。

私が担当した2年生は昨年度はコロナ禍でほとんど大学に来れなかったということもあり、対面授業で盛り上がって大声で話し出したらどうしようと危惧していて(実際そうなったとしたら流石に注意せざるを得ないが彼らの心情を考えると非常に心苦しい)悩んでいた時期もありましたが、基本的には近くの2~3人組で相談し、意欲のある学生は難問の相談を行う感じになってくれており、当初危惧していたようなこともなかったため非常に安堵しています。

私の周りや受講生で感染者が出たという話は今のところなく、おそらく2週間後くらいまでは完全に安心できないですが、ひとまず講義中にこのような問題が起きなかったことが一番ホッとしています。



2021年4月2日金曜日

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

単純無向グラフ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})時間で組合せ的 (高速行列積という飛び道具を使わず初等的に)解くことができます。HLに分解したあと、H内部の四角形を鳩の巣原理に基づくアルゴリズム (昔の記事参照) で判定し、Lの頂点を含む四角形は近傍探索を駆使して判定します。結構実装も簡単でシンプルかつ汎用性がありそうなアイデアなので、競技プログラミングで出ても良さそうな気が個人的にはします(というか既にある??)


参考文献


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

2021年3月31日水曜日

学生生活の振り返り

東京大学大学院 情報理工学系研究科 数理情報学専攻の博士課程を修了しました。自分は学部の4年間を東工大の理学部情報科学科で過ごし、その後に東大大学院に進み修士+博士の5年間を過ごしました。折角なので4+2+3年の学生生活を振り返りつつ、その過程で「もっとこれをすべきだった」や「これはしておいて良かった」と思う事などもまとめます。正直自分には特別なものはなく本当に多くのrejectをくらったりD3の1年間はコロナのせいで出張も行けず辛い思いを多くしてきたのですが、その一方で多くの論文を読んで勉強して研究し続けて理論計算機科学のトップ会議に何本か論文を通すことが出来たのは誇りに思っています。ただ、誰でもそうだとは思いますが、研究をしている最中は論文が通る保証はないわけで、その中で研究を続けていくことに不安を感じることも多々ありました。そんなごくごく普通の学生の吐露からなる記事ですが、興味があればぜひ読んでいってください。

1. D進の理由


自分はB4くらいからD進を考えていて、それ以外の道を考えたことがなかったので特に理由はありません。B4の最初に研究室所属をして初めて研究をした時にかなり楽しく研究できたことが大きかったと思います。学部の間に研究や勉強が楽しいと思えるのは非常に大事なポイントだと思います。今にして思えばD進に際してかなり視野が狭かった(就職に目を向けてインターンなどをしても良かった)気もしますが、逆に迷いなく研究に打ち込めたという面もあるのでどっちもどっちでしょうか。

2. 学生生活で良かったこと

  • とにかく勉強した。学生生活が終わると時間をとって勉強することが出来なくなるのではないかと漠然と思っていたので、将来後悔しないようにとにかくD1の1年間は研究はあまりせずとにかく勉強に打ち込みました。自分の場合は Arora & Barak の Computational Complexity や Frieze & Karonski の Introduction to Random Graphs を約1年くらいかけて読みました。特にランダムグラフの本は自分で手を動かしながら読んだのですが、そのおかげでランダムグラフに関する直感を自分で持つことができて今もたまにその直感が研究に役立つことがあるので、やってて良かったと思っています。
  • 色んな学会に参加した。ありきたりですが、全国の同年代の人たちと交流をもったり先生方と知り合いになれたのは本当に良かったと思っています。また色んな学会で発表をして自分を宣伝できれば、アカデミアの世界での就活で多少は有利に働くことになるんじゃないかと思います。
  • 共同研究をした。自分は博士の3年間はS先生とずっと研究をし、何本も論文を書きました。S先生は自分とは違うバックグラウンドを持つ方で、共同研究をしていく過程で自分一人では到達しえないアイデアを幾つも提案してくださり、結果としてより真理の深淵に近づけたと思います。確かに学会に参加すれば多くの方と交流する機会は得られますが、研究の技術的にディープな内容を討論することは出来ません。それを補いより研究を深めることが出来たというのは良かったと思います。

3. 辛かったこと

  • 「選ばれない」ことのショック。自分は幸運にも、受験などではこれまで「選ばれる」方だったのですが、そのために「選ばれない」経験をほとんどしてこなかったがためにDC1に不採択になった時は精神的ショックはかなり大きかったです。特に一時期は「周囲の知り合いが皆DC1を持っているのに自分だけ...」という部分に大きな負目を感じていました。自分の場合は、1ヶ月後くらいに「自分と他人を比較すべきではないから気にする必要はない。むしろ、これから凄まじい業績を出してDC1とった人たちを超えて当時の審査員を見返そう」と奮起してました。改めて見直すとちょっと引くくらいポジティブですが、なんにせよ失敗をひきずらない精神を持つのは重要だと思います。
  • 将来への不安。ありきたりで誰もが持つとは思います。自分の場合は「考えてもしょうがない。業績が出れば将来有利になるから今は研究に集中しよう」と考え、将来のことはほとんど考えないことにしました。

4. 後悔

  • 一つの究極的な目標をたててそれに邁進する研究をしたかった。自分は主に、個々の小問題を解決していくような研究をしてきました。それは学問における自分の研究の意義を薄くする要因になります。確かに問題を解くことも大事ですが、問題を解いたその後に何を見出すかもまた大事であり、この部分を蔑ろにしてしまったことへの後悔があります。自分が論文の中で新しく提案した部分はその分野においてどのように貢献するかを明白に意識して論文を書くべきでした。自分はこれが(ゼロではなかったにしても)まだまだ足りなかったためにD論でとても苦労しました。博士の3年間はこの力を養うのに最適な期間だと思うので、D進した方々はぜひ頑張ってください。また、将来のことを考えると業績をあげなければと焦るかもしれませんが、その時の流行のトピックにあやかって論文を書くよりも他の誰もやってないオリジナルの観点に基づいた研究を貫いた方が絶対強いと思います。
  • もっと海外に行きたかった。D1の1年間はとにかく勉強に邁進していたため、海外の学会などに参加しませんでした。D2になって国際学会に論文が通ったのですが、予期せぬアクシデント(千葉県を襲った巨大台風のせいで飛行機が飛ばず、代わりの飛行機が取れなかった)のためにこの時は海外に行けず、さらにD3では国際学会に3本論文が通ったのにCOVID19のために全てオンライン開催になってしまいました。そのため自分は修士の間に3回海外に行ったっきりで、Dの3年間は一度も海外に行ってません。今後もCOVID19の影響は続けば海外に行けないのも致し方ないですが、実際に海外のレベルの高い学会に参加してレベルの高い発表を聴いたりすると自分のモチベーションが上がったりするし良い機会になるのでやはり参加できるならばすべきだと思います。

5. 総括


ひたすら研究などに集中できたのは間違いなく周囲の環境のおかげであり、自分一人の力ではここまでやってこれないのは明白で、その環境を用意して下さった指導教員、研究室のメンバー、共同研究者、家族には感謝してもしきれません。

研究トピックを決めるのは本当に難しいと思いますが、これから進学していくみなさんには、単に問題を解くだけではなく、その分野における貢献を意識しながらトピックを策定していくことをぜひおすすめします(めちゃくちゃ難しいですが、、、)

2021年2月22日月曜日

The 123 Theorem

こんにちは. のぶしみです. 最近Mathlogという数学の記事共有サービスを使って記事を書くことにハマっておりまして, つい最近, エキスパンダーグラフの紹介記事を書きました. しかし今まで5年以上このブログをやってきているので, メジャーなトピックはmathlogで紹介し, マイナーなトピックは引き続きこのブログで紹介していきたいと思います.

そこで今回は123 Theoremと呼ばれる定理を紹介します. 変わった名前ですが, 次の定理です:

定理 (123 Theorem)
X,Yを同じ分布に従う実数上の二つの独立な確率変数とする. このとき
\Pr[|X-Y|\geq 2] \leq 3 \Pr[|X-Y|\geq 1].

123 Theoremという名前は不等式に出てくる三つの定数に由来します. 私がこの定理に初めて出会ったのは2021年1月に研究室の後輩とAlon & Spencer の The Probabilistic Method の輪読していた時に第1章の演習問題として出題されていた時でした. 私は基本的に本を読むときは演習問題は頑張って解く習慣があるので, 頑張って証明しようしたのですがかなり苦戦して, 結局丸一日を溶かして以下のようにして示すことが出来ました. 実は同じ証明が[1]に載っていてビックリしたので紹介します.


証明.
x_1,\ldots,x_m\in\mathbb{R}Xの分布に従って独立にサンプリングされたm個の点とし, T_r=\{(i,j):|x_i-x_j|\leq r\}とする ((i,i)\in T_rとする). まず

|T_2| \leq 3|T_1| ......(1)

を示す. この主張は定理の離散化みたいなもので, 実は(1)を示せば良いということが以下の議論から分かります:
\mathrm{E}[|T_2|]=m+m(m-1)\Pr[|X-Y|\leq 2]
\mathrm{E}[|T_1|]=m+m(m-1)\Pr[|X-Y|\leq 1]
および(1)より
m+m(m-1)\Pr[|X-Y|\leq 2] \leq 3m+3m(m-1)\Pr[|X-Y|\leq 1].
これを整理すると,
\Pr[|X-Y|\leq 2] \leq 3\Pr[|X-Y|\leq 1] + \frac{2}{m-1}.
これが任意のmに対して成り立つので, m\to\inftyとすれば,
\forall\epsilon>0,\,\Pr[|X-Y|\leq 2] \leq 3\Pr[|X-Y|\leq 1]+\epsilon.
即ち定理の主張が従う.

(1)の証明.
有向区間グラフの言葉で説明する. \{x_1,\ldots,x_m\}およびr>0に対し, 有向グラフG_rを, V(G_r)=[m]=\{1,\ldots,m\}, E(G_r)=T_rとする. 特に(i,j)\in E(G_r)ならば(j,i)\in E(G_r)でもあるので, 向きを付けることに本質的な意味はないが, 辺の本数=|T_r|として議論していきたいので有向グラフを考えることにする. つまり, |E(G_2)|\leq 3|E(G_1)|を示せば良い. 

二つのグラフG_1,G_2を考える. 頂点iのグラフG_rにおける出次数を\deg_r^+(i)とする. つまり, \deg_r^+(i)=|\{j\in[m]:(i,j)\in E(G_r)\}|であり, 自己ループも1としてカウントする.

主張(1)の証明はmに関する帰納法で行う. m=1ならT_r=\{(1,1)\}なので(1)は成り立つ. m\geq 2の場合を考える. i\in[m]G_1内で最大出次数を持つ頂点とする. このとき, \deg^+_2(i) \leq 3\deg^+_1(i)が成り立つことを言う. もし\deg^+_2(i)>3\deg^+_1(i)とすると, 数直線上の区間 [x_i+1,x_i+2][x_i-2,x_i-1] のどちらかに必ず>\deg^+_1(i)個の点がある (下図).

図. もし\deg_2^+(i)>3\deg_1^+(i)なら, より出次数の大きい点がとれる.


これは\deg^+_1(i)が最大であることに矛盾します. 従って \deg_2^+(i)\leq 3\deg_1^+(i)です. G_1G_2からそれぞれ頂点iを削除して得られるグラフをG'_1, G'_2とすれば, 帰納法の仮定より|E(G'_2)|\leq 3|E(G'_1)|であり, |E(G_r)|=|E(G'_r)|+2\deg^+_r(i)+1なので,

|E(G_2)|\leq 3|E(G'_1)|+6\deg^+_1(i)+1<|E(G_1)|

となり(1)が成り立つ.


補足


123 Theoremはタイトな例を作ることができ, 例えば\Pr[|X-Y|\leq 2] \leq 2.99\Pr[|X-Y|\leq 1]は一般には成り立ちません: \{2,4,6,\ldots,2n\}上の一様分布を考え, この分布に従う独立な二つの確率変数X,Yを考えます. このとき,
\Pr[|X-Y|\leq 2] = \frac{3}{n}-\frac{2}{n^2},
\Pr[|X-Y|\leq 1] = \Pr[X=Y]=\frac{1}{n}
より, \Pr[|X-Y|\leq 2] = 3\Pr[|X-Y|\leq 1] - \frac{2}{n^2}.


参考文献


[1] N. Alon and R. Yuster. The 123 Theorem and its extensions, Journal of Combinatorial Theory, Series A, 72(2)322--331, 1995.

2021年2月15日月曜日

ランダム正則グラフの sandwich conjecture

ランダム正則グラフ理論において有名な sandwich conjecture と呼ばれる予想があるのですが, 最近この予想に大きな進展[3,4]が見られたので紹介します. このトピックは私が修士の頃から論文を読んで追っていたトピックなので, 今回の進展をもたらした論文の登場は個人的にはかなり衝撃的でした.


1. ランダム正則グラフとErdős–Rényiグラフ


1.1. 定義

n 頂点 d-正則の頂点ラベル付きグラフ全体の集合から一様ランダムに取り出したグラフをランダムd-正則グラフと呼び, G_{n,d} で表します. また, n頂点の各頂点のペアを独立に確率pで辺で引いて得られるランダムなグラフを Erdős–Rényiグラフと呼び, G(n,p) で表します. G(n,p)の基本的な性質などは以前の記事を参照してください.


1.2. ランダムグラフの解析

ランダムグラフの解析をする上でそのランダムグラフの生成モデルに対する考察が不可欠です. たとえば, 

p=o(n^{-1}) に対して G(n,p) は確率 1-o(1) で三角形を含まない

という事実が知られていますが, これは, XG(n,p)内に含まれる三角形の個数と定めたときに, Markovの不等式より

\Pr[X\geq 1] \leq \mathbb{E}[X] = \binom{n}{3}p^3 \leq (np)^3 = o(1)

となることから従います. G(n,p)を考えると\mathbb{E}[X]の期待値が非常に計算しやすいためにこのような簡単に証明できるわけです.

それでは, G_{n,d} に含まれる三角形の個数 X に対して \mathbb{E}[X] を求める時はどのようにすれば良いでしょうか? G(n,p)の場合は各辺が独立に出現していたので簡単に \mathbb{E}[X] を求められましたが, G_{n,d} はそうとはいきません.

実は, G_{n,d} には configuration model と呼ばれる生成モデルが知られていて, これを使えば 2^{O(d^2)} nd 時間でランダム d-正則グラフを生成することができます. 技術的な詳細は省きますが, このconfiguration model に基づけば d=d(n)n に依存しない定数ならば G_{n,d} のさまざまな構造的性質を用意に解析することができます. 例えば, dが定数のときは n\to\infty の漸近で三角形の個数X はポアソン分布に従うことが証明できます. また, d=d(n)o(\sqrt{\log n}) くらいまでならなんとかできることもありますが, 例えば d=(1+o(1))n^{1/3} などの場合は難しくなります.

d=d(n) が大きいときのG_{n,d}の生成はそれ自体が一つの研究トピックになるほど難しい問題になっていて, 例えば d=o(n^{1/2}) に対して G_{n,d}O(nd^3) 時間でサンプリングする論文がFOCS15に採択されています[6]. このトピックは多くの論文がありますが, 大体は switching method と呼ばれる手法に基づいたものになっていて, この手法を考えることでG_{n,d} の解析を (時にはかなりアドホックな発想も必要になりつつも) 行うことができます ([2]の10章参照). しかしながら switching method を以てしても d=o(n^{1/2}) などの条件を課す必要があります. 

すなわち, d=d(n) が大きくなるにつれてG_{n,d}の生成と解析は困難なものになっていきます.


1.3. 構造的類似性

G(n,p) の平均次数 npnp=d=\omega(\log n) を満たすとき, G_{n,d}と近い構造を持つであろうことが以下の議論から推察できます: まず, G(n,p)の頂点 v を固定しその次数 \deg(v) を考えます. この値はG(n,p) がランダムに生成されるので確率変数となっていますが, \deg(v)=\sum_{u\in V\setminus \{v\}} \mathbb{1}_{uv\in E} と書けてこれは独立確率変数の和になっているため, 期待値 (n-1)p\approx np に集中します. 具体的には, Chernoff bound (補足の章の補題A.1参照) とUnion boundを組み合わせると, 任意のt>0に対し \Pr[\forall v\in V, |\deg(v)-(n-1)p|\geq t] \leq 2n\exp\left(-\frac{t^2}{2(n-1)p+2t/3}\right)が示せます. 特に, np=\omega(\log n) が成り立つときに t=10\sqrt{np\log n} を代入すると, 非常に高い確率で, 全ての頂点の次数が np\pm 10\sqrt{np\log n} の範囲におさまることが示せます. この議論から, G(n,p) は "ほぼ" (np)-正則グラフとみなすことができるので, np=dのときは G(n,p)G_{n,d}は似た構造的性質を持つことが予想されます. 実際, 彩色数や最大クリークのサイズなどは漸近的にほぼ同じ値を持つことが知られています.

(*) あくまでもこの類似性は d=np が大きい場合にのみ成り立つことに注意してください. 例えば d=np=3 のとき, G_{n,d} は確率 1-o(1) で連結である一方で G(n,p) は確率 1-o(1) で非連結です.


2. Sandwich conjecture


1章で述べた類似性をちゃんと議論した研究が幾つかあります [1,3,4,5]. 一番最初にこの類似性を研究したのは Kim & Vu [5] です. 彼らの定理を述べるために, \mathcal{G}(n,p)G(n,p) の分布, \mathcal{G}_{n,d}G_{n,d} の分布とし, X\sim \mathcal{G} と書いた時は確率変数 X は分布\mathcal{G}に従うことを意味します. また, 次の定理ではランダムグラフのカップリングについて考えていきます. カップリングの定義などについては補足A.2を参照ください.


定理1 (Informal; Theorem 2 of [5]).

d=d(n)\in\mathbb{N} を, d=\omega(\log n) かつ d=o(n^{1/3}/\log^2 n) を満たす関数とする.
p_1=p_1(n)=(1-O(\log n/d)^{1/3}) d/n,
p_2=p_2(n)=(1+o(1))d/n
とする. このとき, 次を満たす G_L\sim \mathcal{G}(n,p_1), G_d\sim \mathcal{G}_{n,d}, G_U\sim \mathcal{G}(n,p_2) のカップリング \pi が存在する:
(1) \Pr[E(G_L) \subseteq E(G_d)] = 1-o(1).
(2) 確率 1-o(1)E(G_U) \setminus E(G_d) のなすグラフの最大次数は o(\log n)




定理1の内容はinformalなもの (特に p_2 には実際にはもう少し条件が設定されている) であることに注意してください. 厳密なステートメントを知りたい方は元論文[5]を参照してください.

定理1から, G_{n,d}d が定理1の条件を満たす場合は p\approx d/n に対し G(n,p) を部分グラフとして含むということが従います. 例えばこの G(n,p) が直径\leq 10 であったならば, これを含む G_{n,d} も直径\leq 10 になることが直ちに従います. より一般的に, 辺の追加における単調性を持つような性質G_{n,d}が持つという証明をしたい場合は対応する G(n,p) で考えれば良いという議論になるので, ランダム正則グラフの構造的性質の証明が一気にしやすくなるわけです.

その一方で, (2)の結果では単に G(n,p_2)\setminus G_{n,d} の最大次数が抑えられるということしか主張しないので, G(n,p_2)直径>10であることが言えたとしてもそこからG_{n,d}が直径>10になることが従うわけではありません. つまり, 定理1には

d=d(n) に条件がついている (特に, d=o(n^{1/3})じゃないといけない).
・包含関係が片方しか示されていない.

という二つの問題点が残っていることになります. Kim と Vuは sandwich conjectureと呼ばれる次の予想を[5]で提示しています:

予想 (Sandwich conjecture, Conjecture 1 of [5]).

d=d(n)=\omega(\log n) とする. ある関数 p_1=(1-o(1))d/n, p_2=(1+o(1))d/n に対して, 以下を満たす G_L\sim \mathcal{G}(n,p_1), G_d\sim \mathcal{G}_{n,d}, G_U\sim \mathcal{G}(n,p_2) のカップリングが存在する:
\Pr[G_L \subseteq G_d \subseteq G_U] = 1-o(1).



3. Sandwich conjecture 解決への進展


Sandwich conjectureに対する進展を簡単に紹介していきます.

・Kim & Vu [5] では上にあげた二つの問題点を残していました.

・Dudek, Frieze, Ruciński, and Šileikis [1] は, d=\omega(\log n) かつ d=o(n) ならば, あるp_1=(1-o(1)) d/n に対して 定理1 の(1) の条件を満たすようなカップリングがあるということを示しています. 実際には[1]ではランダムグラフを一般化したランダムハイパーグラフを考えてカップリングの存在性を証明しています. ランダムグラフの本[2]の10章にこの結果の説明と証明が載っているので気になる方は参照ください. 実は, 私が2018年にSODAに通した論文も[1]の結果を用いています.

・Gao, Isaev, and McKay [4] は d=\omega(n/\sqrt{\log n}) に対して Conjecture 1 が真であることを主張しています. Kim and Vu の問題点だった逆方向の包含関係を (dが非常に大きいという仮定の下ではありますが) 解決したものになっています.

・Gao [3] は d=\Omega((\log n)^7) に対して Conjecture 1 は真であることを主張しています. [4]では d の仮定が非常に強いものでしたが, これをほぼ解決したということになります.


4. まとめ

ランダム正則グラフが研究に出てくるときは, p=d/n に対して G(n,p) を考えると良いです. もし G(n,p) が望ましい性質を持つならば, sandwich conjecture より, 対応する G_{n,d} も同様の性質を持つことが言えるかもしれません.



A. 補足


A.1. 集中不等式

補題A.1 (Chernoff bound; Theorem 21.6 of [2]).

X_1,\ldots,X_nを独立な確率変数で, 0\leq X_i\leq 1 とする (必ずしも同一である必要はない). S=\sum_{i=1}^n X_i とし, その期待値を \mu とする. このとき, 任意のt>0に対して以下が成り立つ:

\Pr[S\geq \mu+t] \leq \exp\left(-\frac{t^2}{2\mu + 2t/3}\right),

\Pr[S\leq \mu-t] \leq \exp\left(-\frac{t^2}{2\mu + 2t/3}\right).


A.2. カップリング


XYをそれぞれ分布 \mathcal{D}_X,\mathcal{D}_Y に従う確率変数とします. XYの同時分布 \pi であって, XYの周辺分布\pi_X, \pi_Y がそれぞれ\mathcal{D}_X, \mathcal{D}_Yと等しいものを, XYカップリングと呼びます. この記事ではG(n,p)G_{n,d}をそれぞれ確率変数とみなし, ランダムグラフ同士のカップリングを考えることがあります. また, 確率変数と分布を区別するために, G(n,p)G_{n,d}の分布をそれぞれ \mathcal{G}(n,p), \mathcal{G}_{n,d} と書くことにします.

例1: G(n,p)G(n,p) のカップリング.

\mathcal{G}(n,p)に従って得られた二つのグラフ G_1,G_2 が, G_1=G_2を満たすようなカップリングが存在します. 最初にXとしてG(n,p)を考え, 次にY=Xとするだけです. このカップリングにしたがって得られた二つのグラフ G_1,G_2は必ず G_1=G_2 となっており, G_1G_2の周辺分布はどちらも\mathcal{G}(n,p) になっています. 自明な例です.

例2: G(n,p_1)G(n,p_2) のカップリング.

p_1\leq p_2 とします. G_1\subseteq G_2 となるような G_1\sim \mathcal{G}(n,p_1)G_2\sim \mathcal{G}(n,p_2) のカップリングが存在します. まず, G_2\sim\mathcal{G}(n,p_2)を生成します. そして得られたグラフG_2の確辺を確率 1-p_1/p_2 で削除して得られたグラフをG_1とします. このように生成したグラフの組 (G_1,G_2)G_1\subseteq G_2の条件をみたします. G_2の周辺分布は明らかに\mathcal{G}(n,p_2)です. 最後にG_1の周辺分布を考えます. G_1の各辺は独立に現れ, その確率はp_2\cdot (1-(1-p_1/p_2))=p_1 となり, 確かにG_1の周辺分布は\mathcal{G}(n,p_1)となります.

この例で紹介したカップリングの存在生から, 例えば次の結果が用意に証明できます:

命題 A.2.

G(n,p_1)が確率1-o(1)で連結であるとする. このとき, 任意の p_2\geq p_1 に対して G(n,p_2)は確率1-o(1)で連結である.

証明:
例2のカップリング\piを考える. このとき
\Pr[G(n,p_2)\text{ is connected}] = \Pr_{(G_1,G_2)\sim \pi}[G_2\text{ is connected}]\geq \Pr_{(G_1,G_2)\sim \pi}[G_1\text{ is connected}]
となり, G_1\sim\mathcal{G}(n,p_1)は確率1-o(1)で連結であることから主張は直ちに従います. 上式の不等号はG_1\subseteq G_2から従います.

一般に, 連結性といった, 辺の追加でinvariantなグラフの性質ならどんなものでも命題A.2のような結果が例2のカップリングを使えばすぐに示せます.



参考文献

[1] Andrzej Dudek, Alan Frieze, Andrzej Ruciński, and Matas Šileikis. Embedding the Erdős–Rényi hypergraph into the random regular hypergraph and Hamiltonicity, JCTB, 2017.

[2] Alan Frieze, and Michał Karoński. Introduction to random graphs, Cambridge University Press, 2015

[3] Pu Gao, Kim-Vu's sandwich conjecture is true for all d=\Omega(\log^7n)arXiv, 2020. https://arxiv.org/abs/2011.09449

[4] Pu Gao, Mikhail Isaev, and Brendan D. McKay. Sandwiching random regular graphs between binomial random graphs. In Proceedings of SODA, 2020.

[5] Jeong Han Kim and Van H Vu. Sandwiching random graphs: Universality between random graph models. Advances in Mathematics, 2004.

[6] Pu Gao and Nicholas Wormald. Uniform Generation of Random Regular Graphs, In Proceedings of FOCS, 2015.

講義資料をNotionで書いてみた

 プログラミング応用という名前の講義を受け持っており, そこで組合せ最適化のベーシックな話題とそのPythonでの実装を教えているのですが, 資料をNotionで書いてみました. 講義資料をNotionで公開しているのでアルゴリズムの基礎とかNP困難性とかを勉強したい人はどう...