量子コンピュータでグラフの最大クリークを発見する方法

Covid-19騒ぎに乗じて二週間ほど休みをもらいましたので、最初の3日ほどかけてUE4の再入門をやりました。
結果、最適化するゲームでも作りたいなと思い、Ocean-SDKとPython再入門しました。


# 目的
とにかくここ10年ほとんどまともにプログラミングしていないので、Python＋OceanSDKでなにも見ずに実装できるようになりたい。

# やったこと
Fixstarsさんの”組合せ最適化問題のイジング模型定式化”ページで問題をみてとりあえずおおよそ1日1個 計11個の定式化を実装していった。
どうしても解けない問題については、答えを見て実装した。

https://quantum.fixstars.com/techresouces/research/ising-model-formulation/clique-judgment/

# 結果
何も参考にせずに超基本は実装できるようになった。これだけでは腐れポエムで終わるので以下、実装例。

# 実装例

## 問題

無向グラフ$$G=(V,E)$$と自然数Kが与えられます。 Gには大きさKのクリークが存在するでしょうか？

なお、この問題の最適化問題版はNP困難な問題として知られています。 （無向グラフの大きさ最大のクリークを求める問題です。）

## 問題の定式化

### 解法検討

サイズKのクリークを探すイジングはよく見るので、一気に最大クリークを探し出す方法を検討する。その前に前提として、クリークのエッジ数とノード数の間には

$$クリークのエッジ数= 対象のノード数*(対象のノード数-1) / 2 $$

という関係ある。

i番目のノードが クリークを形成するノードとして選択された場合$x_i=1$ とした場合、目的関数と制約式は以下の通り。

#### 目的関数

 選択されたノード数の最大化
 $$ H_A = -\sum(x_i) $$

#### 制約式

 選択されたノードがクリークを形成していることは以下の2つの式が同値であることが条件となる。
 
 * 選択されたノード間に存在するエッジ数(高々以下のクリークである場合のエッジ数と同数の正の数)$\sum_{(u,v) \in V}( x_u * x_v)$ 
 * クリークである場合のエッジ数 $ \sum x_i( \sum x_i-1)/2$ 

よって、制約式は
$$ H_B =  \sum x_i( \sum x_i-1)/2 -\sum_{(u,v)\in V}( x_u * x_v) $$ 

## 実装

```python:問題作成
#まずはランダムな無効グラフを作成する
%matplotlib inline 
import networkx as nx 
import dwave_networkx as dnx 
from matplotlib import pyplot as plt 

NUM_NODES=11
G = nx.fast_gnp_random_graph(NUM_NODES,0.6,seed=1)

plt.figure(figsize=(8,8))
pos = nx.spring_layout(G,iterations=200,seed=1)
nx.draw_networkx(G,pos)
```
![download.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/121950/99f78c61-ae9c-3279-cf9b-08dca1ceab07.png)

先に結果の可視化関数を定義しておきます。
※もっとスマートに実装できる気もしますが…

```python:結果の可視化
#結果表示用ルーチン
import itertools
def show_result(sol):
    Target_node = []
    Target_edge = []

    #クリークの制約を満たしているか
    if len(sol[1]) != 0:
        print("Failed")
        print(sol[1])

    #可視化用のデータ作成
    solution = answer.first.sample
    for i in range(NUM_NODES):
        if sol[0]['x'][i] == 1:
            Target_node.append(i)

    for i in range(NUM_NODES):
        if sol[0]['x'][i] == 1:
            for edge in G.edges(i):
                if edge in itertools.permutations(Target_node,2): 
                    Target_edge.append(edge)

    plt.figure(figsize=(8,8))
    nx.draw_networkx(G,pos)
    nx.draw_networkx_nodes(G,pos,nodelist=Target_node,node_color='red')
    nx.draw_networkx_edges(G,pos,edgelist=Target_edge,edge_color='red')
```


”解法検討”で考えた内容をそのまま実装していきます。

```python:PyQUBOを用いた定式化
from pyqubo import Array,Sum,Constraint,Placeholder

x = Array.create('x',shape=NUM_NODES,vartype='BINARY')
H_A = -Sum(0,NUM_NODES,lambda i: x[i])

H_B = Constraint(Sum(0,NUM_NODES,lambda i: x[i])*(Sum(0,NUM_NODES,lambda i: x[i])-1)/2 - Sum(0,NUM_NODES,lambda i: Sum(0,NUM_NODES,lambda j: x[i]*x[j]*int(G.has_edge(i,j))))/2,label='Clique constraint was not satisfied')

H = H_A + Placeholder(label='balance') * H_B
model = H.compile()
```

上記で定式化しModelを作成しました。Placeholderを用いることで二つの式H_A、H_BのバランスをModel.compile()なしに調整できます。
何がうれしいかというと、複雑な式になるとcompile()で非常に時間がかかるのですが、それを回避できるため、パラメータチューニングで時間が無駄になりません。

Placeholderは以下のようにdict形式で定義できます。

```python:D-WaveのHybridSolverを使用して解く
from dwave.system import LeapHybridSampler,DWaveSampler
#Set Placeholder
feed_dict = {'balance': 1.2}
bqm = model.to_dimod_bqm(feed_dict)

sampler = LeapHybridSampler()
answer = sampler.sample(bqm,50) 
```

D-WaveのSamlerは処理結果をSampleSet型で返します。HybridSolverではあまり気にする必要がありませんが、せっかくなので便利なPyQUBOでデコードして、先に定義した可視化の関数をコールします。

```python:D-WaveのHybridSolverを使用して解く
solutions = model.decode_dimod_response(answer,topk=5,feed_dict=feed_dict)
show_result(solutions[0])
```

![download.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/121950/63cbfa9c-a221-d7ec-b576-da49123b5e54.png)

はい、とけました。
なんかいい感じになっているので、非常に気分がいいです！