Goemans-Williamson Algorithmを理解したい

この記事はJij Advent Calendar 2024の6日目の記事です．

QAOAやQuantum Annealingでは，イジングモデルでかけるという理由で，例題でよくMaxcut問題を扱うことが多い．その時に，性能のベースラインとしてよく使われるのが，Goemans-Williamson Algorithmである．このアルゴリズムはそれ自体が非常に面白いので，この記事はGoemans-Williamson Algorithmについて自分が勉強したまとめである．

Maxcut問題

Maxcut問題は，無向グラフ $G(V,E)G(V,E)$とエッジの重み $w:E→R+w:E\to \R_{+}$が与えられたとき，ノードを $U⊂VU\subset V$と $W⊂VW \subset V$に分け，両端が $UU$に属するノードと $WW$に属するノードになっているエッジの合計を最大化する問題である．

例えば，以下の左の図のようなグラフが与えられたとする．
右の図のようにノードを塗り分ける(組を作る)ことで，ほとんどの隣り合うノードが別の組に入る分け方が実現できることがわかる．

Maxcut問題は，以下のように定式化することができる．
ノード $viv_i$が $UU$に属するとき $si=1s_i=1$， $WW$に属するとき， $si=−1s_i=-1$となる変数 $si∈{1,−1}s_i \in \{1,-1\}$を定義する．
このとき，エッジの両端 $(i,j)∈E(i,j)\in E$が同じ集合に含まれているとき， $sisj=1s_is_j = 1$であり，異なる集合に含まれているとき $sisj=−1s_is_j = -1$である．
この性質を用いると，次の数理モデルを作ることができる．

$$max⁡12∑(i,j)∈Ewij(1−sisj)(1) \tag{1} \begin{split} \max\quad& \frac{1}{2}\sum_{(i,j)\in E}w_{ij}(1-s_is_j)\\ \end{split}$$

ここで係数の$1/21/2$は，異なる組に入っているときに$1−sisj=21 - s_is_j = 2$となるので，それを補正するためにかけている．

Goemans-Williamson Algorithm

Goemans-Williamson Algorithm(GW)は，Maxcut問題に対する近似アルゴリズムの一つである[1]．このアルゴリズムについては，[2-4]に良い説明があるので，そちらを読む方が有意義かなとは思う．ここでは，[2-4]をベースにして,このアルゴリズムを説明しようと思う．

GWは，

• 問題(1)に対する適当な緩和問題を作る
• 緩和問題の解をうまく丸めて元の問題の解になるようにする

という二つのステップで構成されている．

緩和問題の作成

まずは，問題(1)をベクトル計画問題に緩和した問題を考える．
ベクトル計画問題は，変数がベクトル $vi∈Rn\bm{v}_i\in \R^n$で表現され，目的関数と制約条件がベクトルの内積 $vi⋅vj\bm{v}_i\cdot\bm{v}_j$に対して線形で表現される問題である．
問題(1)をベクトル計画問題に緩和した問題は，次のように記述することができる．

$$max⁡12∑(i,j)∈Ewij(1−vi⋅vj)s.t.vi⋅vi=1 ∀ivi∈Rn ∀i(2) \tag{2} \begin{split} \max\quad& \frac{1}{2}\sum_{(i,j)\in E} w_{ij}(1-\bm{v}_i\cdot\bm{v}_j)\\ \mathrm{s.t.}\quad &\bm{v}_i\cdot\bm{v}_i = 1\ \forall i\\ &\bm{v}_i \in \R^n\ \forall i\\ \end{split}$$

問題(1)のどんな実行可能解 $yiy_i$も$vi=(yi,0,…,0)v_i = (y_i,0,\dots,0)$とすれば， $vivi=1v_iv_i = 1$と $vivj=yiyjv_iv_j = y_iy_j$を満たし，ベクトル計画問題(2)の実行可能解になるので，ベクトル計画問題は元の問題の緩和問題になっていることがわかる．
また，この問題は $vi⋅vi=1\bm{v}_i\cdot\bm{v}_i = 1$という制約条件から， $vi\bm{v}_i$は $Rn\R^n$空間上の単位円 $SnS_n$上に解は分布することもわかる．

次に，上記のベクトル計画問題をさらに半正定値計画問題に変換する．
上記のベクトル計画問題では，変数は内積$vi⋅vj\bm{v}_i\cdot\bm{v}_j$の形でしか表れていないので，内積 $vi⋅vj\bm{v}_i\cdot\bm{v}_j$を半正定値計画の変数 $yijy_{ij}$に変換する．これによって，目的関数と制約条件は $yijy_{ij}$に対して線形となる．さらに， $yijy_{ij}$を要素にもつ行列 $Y∈Rn×nY\in \R^{n\times n}$は，簡単な計算から$YY$は半正定値行列であることがわかる．

以上から，ベクトル計画問題(2)を半正定値計画問題へと変換した結果得られる問題は，

$$max⁡12∑(i,j)∈Ewij(1−yij)s.t.yii=1 ∀iY⪰0Y∈Mn(3) \tag{3} \begin{split} \max\quad& \frac{1}{2}\sum_{(i,j)\in E}w_{ij}(1-y_{ij})\\ \mathrm{s.t.}\quad &y_{ii} = 1\ \forall i\\ &Y \succeq 0 \\ &Y \in M_n \end{split}$$

とかける．ここで， $MnM_n$は実 $n×nn \times n$行列の対称錐の集合であり，$vi⋅vi=1\bm{v}_i\cdot \bm{v}_i = 1$の制約条件から，$yii=1y_{ii} = 1$という制約条件をかけている．
この半正定値計画は内点法などを用いることで多項式時間で解くことができる．

実際に，上記の半正定値計画問題をソルバーで解きやすくするために標準形に持っていく．
これにはグラフラプラシアン$LL$を用いるのが便利である．
グラフラプラシアン$LL$は，隣接行列$AA$と次数行列$DD$を用いて，$L=D−AL = D - A$と書くことができ，このグラフラプラシアンを用いて，式(3)は

$$max⁡12L∙Ys.t.yii=1 ∀iY⪰0Y∈Mn(4) \tag{4} \begin{split} \max\quad& \frac{1}{2}L \bullet Y\\ \mathrm{s.t.}\quad &y_{ii} = 1\ \forall i\\ &Y \succeq 0 \\ &Y \in M_n \end{split}$$

と書き直すことができる．

ここでは，上で挙げた簡単な例を解いてみることにする．

import networkx as nx
G = nx.Graph()
G.add_nodes_from([0, 1, 2, 3, 4])
G.add_edges_from([(0, 1), (0, 4), (1, 2), (1, 3), (2, 3), (3, 4)])
pos = {0: (1, 1), 1: (0, 1), 2: (-1, 0.5), 3: (0, 0), 4: (1, 0)}
nx.draw(G, pos, with_labels=True)

グラフラプラシアンは

L = nx.laplacian_matrix(G, nodelist=list(range(5))).todense()

で得ることができる．

さて，これをソルバーで解いてみよう．
ここでは，cvxpyで実装して，それに付随しているソルバーで解くことにする．

import cvxpy as cp

Y = cp.Variable((5, 5), PSD=True)
constraints = [cp.diag(Y) == 1]
problem = cp.Problem(cp.Maximize(cp.trace(L @ Y) / 2.0), constraints)

problem.solve()

で解くことができる．

ラウンディングアルゴリズム

問題(3)の半正定値計画の解は得られたとして，次にこの解をどのように元の問題の解に戻すのかを考えよう．
得られた解 $YY$を特異値分解などを使えば，$Y=VTVY=V^TV$となるような $VV$を得ることができるので，ベクトル計画問題(2)の解を得ることができる．

U, singular_values, Vh = np.linalg.svd(Y.value)
V = (U @ np.diagflat(np.sqrt(singular_values))).T

さて，ベクトル計画問題(2)の最適解 $u1,…,un\bm{u}_1,\dots,\bm{u}_n$を得た．ここからどのように，(1)の解を得るかを考える．
二つのベクトル $ui,uj\bm{u}_i,\bm{u}_j$のなす角度を $θij\theta_{ij}$とすると， $∣u∣=1|{\bm{u}}| = 1$なので， $ui⋅uj=cos⁡θij\bm{u}_i\cdot\bm{u}_j = \cos\theta_{ij}$．
よって，ベクトル計画問題の目的関数は，

$$12∑(i,j)∈Ewij(1−cos⁡θij) \frac{1}{2}\sum_{(i,j)\in E}w_{ij}(1-\cos\theta_{ij})$$

となる．この目的関数を最大化するには， $θij\theta_{ij}$が $π\pi$に近い方がよい．つまり， $θij\theta_{ij}$が大きいときに $ui,uj\bm{u}_i,\bm{u}_j$を別の組になるようにすればいい．

そこで，以下の方法で，indexの集合を分類することにする．
初めに$nn$次元球面上から一様ランダムにベクトル $r\bm{r}$をサンプルする．このベクトルは各成分を正規分布 $N(0,1)\mathcal{N}(0,1)$からサンプルし，規格化することで得ることができる．
この$r\bm{r}$を用いて， $vi⋅r≥0\bm{v}_i\cdot \bm{r}\geq 0$となるベクトルのindexを $UU$に入れ，それ以外を $WW$に入れることで，全てのノードindexを $UU$と $WW$に分類することができる．

この方法についてもう少し詳しく述べよう．$ui\bm{u}_i$と $uj\bm{u}_j$が分離されるのは，ベクトル $r\bm{r}$に対して，

• $ui⋅r≥0\bm{u}_i\cdot \bm{r}\geq 0$かつ $uj⋅r≤0\bm{u}_j\cdot \bm{r}\leq 0$
• $ui⋅r≤0\bm{u}_i\cdot \bm{r}\leq 0$かつ $uj⋅r≥0\bm{u}_j\cdot \bm{r} \geq 0$

のどちらかの場合である．
$ui\bm{u}_i$と$uj\bm{u}_j$が張る2次元平面を考え， $r\bm{r}$をこの平面に射影することを考える．$r\bm{r}$は球面対称な分布にしたがっているので，この2次元射影した円に対しても一様に分布する．
ここで， $ui⋅r≥0\bm{u}_i\cdot \bm{r}\geq 0$かつ $uj⋅r≤0\bm{u}_j\cdot \bm{r}\leq 0$もしくは， $ui⋅r≤0\bm{u}_i\cdot \bm{r}\leq 0$かつ $uj⋅r≥0\bm{u}_j\cdot \bm{r} \geq 0$となるのは，$r\bm{r}$が$ui\bm{u}_i$と $uj\bm{u}_j$が作る弧の間つまり， $θij\theta_{ij}$の範囲にくる場合である．よって，その範囲に $r\bm{r}$がくる確率は $2θij/2π=θij/π2\theta_{ij}/2\pi = \theta_{ij}/\pi$である．
つまり，$θij\theta_{ij}$が$π\pi$に近いほど分離される確率が高くなることがわかる．

実際に，この方法で，(1)の解を得てみよう．

r = np.random.normal(size = (5,)) 
r /= np.linalg.norm(r)
cut = np.ones((5, 1))
for i in range(5):
    cut[i] = -1 if r @ V[:, i] <= 0 else 1

ちなみに，実際に得られた結果をプロットしたのが，冒頭に載せた右の図である．

近似率について

さて，最後に近似率について簡単に議論しておこう．

まず初めに，カットされた重みの期待値 $E[W]\mathbb{E}[W]$を計算してみる．$ui\bm{u}_i$と $uj\bm{u}_j$が分離される確率が$θij/π\theta_{ij}/\pi$であることを思い出すと，

$$E[W]=∑(i,j)∈Ewijθijπ \mathbb{E}[W] = \sum_{(i,j)\in E}w_{ij}\frac{\theta_{ij}}{\pi}$$

である．
ここで，天下り的だが，$α\alpha$が

$$α=2πmin⁡0≤θ≤πθ1−cos⁡θ≈0.87856… \alpha = \frac{2}{\pi}\min_{0\leq \theta\leq \pi}\frac{\theta}{1 - \cos \theta} \approx 0.87856\dots$$

で与えられるとしよう．
この定義から$α≤2πθ1−cos⁡θ\alpha \leq \frac{2}{\pi}\frac{\theta}{1 - \cos \theta}$なので，

$$θπ≥α(1−cos⁡θ)2 \frac{\theta}{\pi} \geq \frac{\alpha(1 - \cos \theta)}{2}$$

が成り立つ．これを代入すると，

$$E[W]≥α∑(i,j)∈Ewij(1−cos⁡θij)2 \mathbb{E}[W] \geq \alpha\sum_{(i,j)\in E}w_{ij}\frac{(1 - \cos \theta_{ij})}{2}$$

となる．右辺はベクトル計画問題(2)の最適値 $OPTv\mathrm{OPT}_v$であり，問題(2)は問題(1)の緩和問題であることを使うと，$OPTv≥OPT\mathrm{OPT}_v\geq \mathrm{OPT}$なので，

$$E[W]≥αOPT \mathbb{E}[W] \geq \alpha\mathrm{OPT}$$

が成り立つ．これはアルゴリズムの解の精度と最適値の比に対して，
よって，GWの近似率は$α≈0.878\alpha\approx 0.878$である．

最後に

Jijでは数理最適化・量子コンピュータに関連するソフトウェア開発を行っています。
OSSも複数開発、提供しています。Discordのコミュニティへぜひ参加してください！
https://discord.gg/Km5dKF9JjG

また，Jijでは各ポジションを絶賛採用中です！
ご興味ございましたらカジュアル面談からでもぜひ！ご連絡お待ちしております。
数理最適化エンジニア（採用情報）
量子コンピューティングソフトウェアエンジニア（採用情報）
Rustエンジニア(アルゴリズムエンジニア)（採用情報）
研究チームPMO（採用情報）
オープンポジション等その他の採用情報はこちら

参考文献

[1]Michel X Goemans and David P Williamson. “Improved approximation algorithms for maximum
cut and satisfiability problems using semidefinite programming”. In: Journal of the ACM (JACM)
42.6 (1995)
[2]Vijay V Vazirani. Approximation algorithms. Vol. 1. Springer, 2001.
[3]David P Williamson and David B Shmoys. The design of approximation algorithms. Cambridge
university press, 2011.
[4]Matsui Tomomi. “OR 研究の最前線 半正定値計画を用いた最大カット問題の.878 近似解法”. In:
オペレーションズ・リサーチ = Communications of the Operations Research Society of Japan: 経
営の科学 45.3 (Mar. 2000)