Groverのアルゴリズムで頂点彩色問題を解いてみる

# はじめに

Groverのアルゴリズムは量子コンピュータにおけるアルゴリズムの一種で、古典コンピュータにおける線形探索よりも高速に条件に合致するデータを抽出することができます。計算量で考えると、古典コンピュータの線形探索では$\mathrm{O}(n)$かかるところ、Groverのアルゴリズムでは$\mathrm{O}(\sqrt{n})$で探索ができます。

古典コンピュータの$\mathrm{O}(n)$の時点で十分に高速であり、量子回路にエンコーディングする手間を考えると、一般的なデータベースにおけるデータ検索で活用されるということはあまり想像できませんが、制約充足問題の解の探索のような、一般的に指数時間かかる問題の解の探索において活用できる可能性があります。当然ながら、二乗加速しても、指数時間かかることには変わりませんが、実用においては、最適化問題を解く上で、定数倍の加速にも価値があるケースは多く存在するため、二乗加速も十分に価値がある可能性があります。

QiskitのTextbook[[1]](#参考文献)には、Groverのアルゴリズムを用いて、単純化した数独を解くサンプルが用意されています。また、QiitaでもパズルをGroverで解く試みをされている記事がいくつかあります[[2, 3]](#参考文献)。

その他にも様々な試みがされており、制約充足問題を用途の一つとしては認識していましたが、自分で実際にオラクルの設計から取り組んだことはないため、この記事では、勉強も兼ねてグラフの頂点彩色問題をGroverのアルゴリズムで解くことに取り組んでみようと思います。頂点彩色問題はNP完全であることが知られています。

同様の取り組みは既に存在していますが、ここでは効率などは一旦おいておき、自分なりに実装を考えてみます。

# 頂点彩色問題

グラフの頂点彩色問題は、グラフ$G=(V,E)$と色の集合$C$が与えられた時、各頂点$v \in V$に対し、色$c \in C$を一つ決めます。このとき、隣接する頂点同士が異なる色になるような塗り分けが存在するかどうかを判定する問題です。

有名な例としては、地図の塗り分け問題があり、4色あれば、いかなる地図も塗り分けられることが知られています。これは、平面グラフにおける頂点彩色問題に相当します。また、数独も一種の頂点彩色問題と見なすことができます。

以下のグラフについて塗り分け方の一例を示します。

![元のグラフ](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2670216/3891e1b1-8265-9024-a59b-a39976043cf1.png)

![彩色後のグラフ](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2670216/0f86c6f5-7adc-9151-0399-684f29569709.png)

## 頂点彩色問題の定式化

頂点彩色問題を数式の形で表現してみます。とりあえずは一般的な数理最適化でよくありそうな形で定式化をしてみます。

変数
$x_{v, c}$：頂点$v \in V$を色$c \in C$で彩色する際に1、そうでない場合は0を取る変数

制約式
$ \sum_{c \in C} x_{v, c} = 1 \quad \forall v \in V $ : 各頂点は一色のみで塗られる
$ x_{u,c} + x_{v, c} \leq 1 \quad \forall (u, v)\in E, c \in C $ : 各枝の両端が同じ色で塗られていない（同時に1にならない）

変数に関しては、この形のまま、1変数1量子ビットで割り当てます。制約の実装方法については後述します。


# 実装

先ほどの例で示したグラフを元に、今回実装をします。

今回使った環境は以下の通りです。
```
python 3.11.3
networkx 2.8.4
qiskit 0.45.2
```

まずは必要なライブラリをインポートします。
```Python
import networkx as nx
import qiskit
from qiskit.circuit import QuantumCircuit
from qiskit import Aer
```

次にグラフの情報を定義します。
```Python
G = nx.Graph()

nodes = [0, 1, 2, 3]
edges = [(0, 1), (0, 2), (1, 2), (1, 3), (2, 3)]

G.add_nodes_from(nodes)
G.add_edges_from(edges)
```

可視化して確認してみます。

```Python
pos = nx.spring_layout(G, seed=0)
nx.draw_networkx(G, pos)
```

![元のグラフ](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2670216/3891e1b1-8265-9024-a59b-a39976043cf1.png)

Groverの実装を行うにあたり、必要なパラメータを設定します。

今回は、3色で塗り分けることを考えます。

```Python
n_nodes = len(nodes) # 頂点の数
n_edges = len(edges) # 枝の数
n_colors = 3 # 塗り分ける色の数
n_variables = n_nodes * n_colors # 変数の数
n_ancilla = n_nodes + n_edges # 実装に用いるアンシラの数
n_qr = n_variables + n_ancilla # 必要な量子ビット数

backend_sim = Aer.get_backend('qasm_simulator')
```

次に、制約を満たした状態をマークするためのオラクルを定義します。

まず、各ノードを一色のみで塗る制約についてです。三色の場合、(1, 0, 0), (0, 1, 0), (0, 0, 1)のいずれかの状態であればよいです。

今回は、マルチコントロールゲートを用いて、いずれかの状態であれば、アンシラが1になるように、以下のような実装をしました。

![頂点ごとの色の制約](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2670216/44fd145e-1f17-7971-39d5-1b923d166316.png)

そのほかには、任意の二色の組について、同時に1になっていないことを確認する、などの実装があり得ると思います。シミュレータ上では今回の実装が楽な気がしますが、実機で動かしやすいのはどちらかまでは考慮していません。


次に、エッジの両端の色が異なるようにする制約です。先ほどの定式化で示した形を直接表現することは難しそうです。同時に1になっている場合に、アンシラが0になるようにするため、以下のような実装をしました。

![枝ごとの色の制約](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2670216/eeb110a2-644e-a9b6-75a3-25b0cb4664ab.png)


これらを組み合わせて、以下にオラクルを定義します。

```Python
def oracle():
    circ = QuantumCircuit(n_qr)

    # 各ノードは一色で塗る
    for node in range(n_nodes):
        for i in range(n_colors):
            for j in range(n_colors):
                if i != j:
                    circ.x(node*n_colors + j)
            circ.mcx([node*n_colors + k for k in range(n_colors)], n_variables + node)
            for j in range(n_colors):
                if i != j:
                    circ.x(node*n_colors + j)
        circ.barrier()

    # 隣り合うノードは別の色にする
    ancilla_id = n_variables + n_nodes

    for u, v in edges:
        circ.x(ancilla_id)
        for i in range(n_colors):
            circ.ccx(u*n_colors + i, v*n_colors + i, ancilla_id)
        ancilla_id += 1
        circ.barrier()

    # 正解状態のときは位相反転
    circ.h(n_qr - 1)
    circ.mct(list(range(n_variables, n_qr - 1)), n_qr - 1)
    circ.h(n_qr - 1)
    circ.barrier()
    
    # アンシラの状態を元に戻す
    
    # 隣り合うノードは別の色にする
    ancilla_id = n_variables + n_nodes
    
    for u, v in edges:
        circ.x(ancilla_id)
        for i in range(n_colors):
            circ.ccx(u*n_colors + i, v*n_colors + i, ancilla_id)
        ancilla_id += 1
        circ.barrier()
            
    # 各ノードは一色で塗る
    for node in range(n_nodes):
        for i in range(n_colors):
            for j in range(n_colors):
                if i != j:
                    circ.x(node*n_colors + j)
            circ.mcx([node*n_colors + k for k in range(n_colors)], n_variables + node)
            for j in range(n_colors):
                if i != j:
                    circ.x(node*n_colors + j)
        circ.barrier()
        
    return circ
```

そして、振幅増幅をするための回路を作ります。この部分は問題に依存しない一般的な構造です。

```Python
def amplifier(n):
    circ = QuantumCircuit(n)
    circ.h(range(n))
    circ.x(range(n))
    circ.h(n - 1)
    circ.mct(list(range(n - 1)),n - 1)
    circ.h(n - 1)
    circ.x(range(n))
    circ.h(range(n))
    
    return circ
```


最後に、これらを結合して、一つの回路にまとめます。

操作の繰り返し回数ですが、ここでは10回で試してます。

```Python
n_loops = 10 # オラクルと振幅増幅の操作を何回繰り返すか

circ = QuantumCircuit(n_qr, n_variables)

for i in range(n_variables):
    circ.h(i)
circ.barrier()

for i in range(n_loops):
    circ = circ.compose(oracle())
    circ = circ.compose(amplifier(n_variables))

circ.measure(range(n_variables), range(n_variables))
```

```Python
n_shots = 10000
job_sim = backend_sim.run(qiskit.transpile(circ, backend_sim), shots=n_shots)
result_sim = job_sim.result()
counts = result_sim.get_counts(circ)
```

得られた結果から、出現回数が多いものを確認してみます。

```Python
sorted_counts = sorted(counts.items(), key=lambda x: x[1], reverse=True)[:10]
print(sorted_counts)
# [('010001100010', 908), ('010100001010', 867), ('100010001100', 865), ('001100010001', 858), ('100001010100', 857), ('001010100001', 819), ('011110100010', 7), ('100110011110', 7), ('000011100101', 7), ('110010010010', 6)]
```

上位6件の出現回数が突出しています。結論から言うと、塗り分けの仕方は$3!=6$通りなので、正しく動いていることが期待できそうです。

一方で、この出力を見ただけではうまくいってるのか少しわかりづらいので、グラフに可視化してみます。

```Python
color_dict = {0:"red", 1:"green", 2:"blue"}

def split_key(key):
    # 状態を各頂点に相当する部分ごとに切り分け、1になっている箇所を確認する
    color_node_dict = {i:[] for i in range(n_colors)} # 色をkeyに、頂点のリストを持つ辞書
    for i in range(0, n_variables, n_colors):
        color_string = key[i:i+n_colors] # 出力を切り分ける
        assert color_string.count("1") == 1 # 1色だけで塗られていることを確認する
        for j, c in enumerate(color_string):
            if c == "1": # 色を確認する
                color_node_dict[j].append(int(i / n_colors))
    return color_node_dict
                
def draw_color_graph(nodes, edges, node_color_dict):
    for key, value in color_dict.items():
        print(value, color_node_dict[key])
        nx.draw_networkx_nodes(G, pos, nodelist=color_node_dict[key], node_color=value)
        nx.draw_networkx_edges(G, pos)
```

出現回数が多かったものから順にいくつか確認してみます。
```Python
color_node_dict = split_key(sorted_counts[0][0]) # '010 001 100 010'
draw_color_graph(nodes, edges, color_node_dict)
```

![彩色後のグラフ 1](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2670216/0cb65f58-b5e0-319b-4b2d-426885a7f519.png)

```Python
color_node_dict = split_key(sorted_counts[1][0]) # '010 100 001 010'
draw_color_graph(nodes, edges, color_node_dict)
```

![彩色後のグラフ 2](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2670216/e2685f9e-9fbe-215f-1946-4cb0c155b9ed.png)

```Python
color_node_dict = split_key(sorted_counts[2][0]) # '100 010 001 100'
draw_color_graph(nodes, edges, color_node_dict)
```

![彩色後のグラフ 3](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2670216/2d1f9b78-d307-ab10-b62d-87c6eb675e13.png)

いずれも制約を守って正しく塗り分けできていそうです。

# まとめ
今回はグラフの頂点彩色問題を題材に、Groverのアルゴリズムの実装を試してみました。ひとまず動くものができたことには満足しています。

問題を数式の形で表現することは容易にできても、これをどのようにオラクルに落とし込むかは自明ではないように思え、このあたりを試行錯誤してみるのは面白かったです。今回のようなOne-Hot制約など、一般的によく現れる制約についてはある程度定番の形が定まっていそうな気はしますが、アンシラの数と回路の深さのトレードオフで色々と考えられそうな気はします。

また、今回の実装を試した後に、グラフの頂点彩色問題のGroverの実装を調べてみると、各頂点につき、色の塗り方をバイナリにエンコーディング(2ビットの場合、00〜11で4色)するやり方をしている方[[4]](#参考文献)もいたり、色々とやり方がありそうです。このやり方であれば、必要な量子ビット数が減りますし、頂点ごとのOne-Hot制約がいらなくなるのは強みです。通常の数理最適化のような定式化にとらわれず、Grover向きのやり方を考えるスキルも必要そうだと感じました。



# 参考文献

1. [Qiskit Textbook - Grover’s Algorithm](https://github.com/Qiskit/textbook/blob/main/notebooks/ch-algorithms/grover.ipynb)　（アクセス日：2024/09/18）

2. [@tax_nalgo, "量子コンピュータで論理パズルを解いてみた", Qiita, 2023.](https://qiita.com/tax_nalgo/items/48cf1ee47549fbf3b386)（アクセス日：2024/09/18）

3. [@mi_yuyu, "Groverのアルゴリズムで簡単なパズルを解く", Qiita, 2023.](https://qiita.com/mi_yuyu/items/7a5756c8462fad8eefc3)（アクセス日：2024/09/18）

4. [Mohammad Shoaib, "How to solve the Graph Coloring Problems using Qiskit’s Grover algorithm", Medium, 2023.](https://medium.com/@shoaib6174/how-to-solve-the-graph-coloring-problems-using-qiskits-grover-algorithm-e87eb52f203)（アクセス日：2024/09/18）