Groverアルゴリズムでグラフ彩色問題を解く by Qmod

こんにちは。元々Groverを勉強しており、IBM2019のチャレンジ問題でグラフ彩色問題のコンストがあったこともあり気になっていたものの勉強したことなかったのでトイプロブレムでやってみます。数式は抜きでざっくりやります。

グラフ彩色問題

無向グラフに対して制約条件を満たしながら隣接する頂点は同じ色にならないように塗り分ける問題です。今回は以下の図の４色塗り分けのトイプロブレムを考えます。

図のように黒点４色問題を隣接ノード同士が同じ色にならない解をGroverで探します。

Groverアルゴリズム

Groverアルゴリズムは、量子計算における代表的なアルゴリズムの一つで、未整列データベースの探索問題を効率的に解くために設計されています。このアルゴリズムでは、特定の「ターゲット状態」を見つけるために、アダマールゲート (H)、オラクル (R)、およびグローバー拡散演算子 (D) を使います。

ブロック図（量子ゲート）を使って説明します。n個の量子ビットに対してアダマールゲートを作用させ全ての入力状態を最初に用意します。次に入力状態の中から制約を満たすものだけをオラクルRを使ってマーキングします。したがって、全ての入力に対して並列にマーキング処理をします。オラクルは問題ごとによって異なっており、今回の場合、グラフ彩色問題の制約条件を満たすようにユーザが設計する必要があります。

ただし、量子ではマーキングだけしても意味がなく（詳細な数式は省きます。）、そのためマーキングした状態に対して増幅をさせて他の入力状態と区別します。これはグローバー拡散演算子 (D)を使用することで達成できます。Dは初期状態で表記が可能で、設計方法は固定されています。このRとDを交互に作用させると増幅度合いが強まっていきますが、作用させすぎると精度が悪化します。必要な回数は$\sim \mathcal{O}(\sqrt{2^n/M})$で達成できます($M$は制約を満たす充足解の個数)。詳しい回路を知りたい方は例えば、Qiskitに詳細があります。

Qmodをつかって実装

今回はQiskitではなくQmodを使って実装してみます。以下のドキュメントつかってやってみます。

この文章を読むとQmodはgrover_searchというパッケージが用意されており、

grover_search(
        reps,
        oracle=lambda vars: phase_oracle(oracle, vars),
        packed_vars = vars,
    )

でGroverアルゴリズムが実装可能です。repsはGroverオペレータの繰り返し回数, oracleはRで, varsはオラクルを作用させる量子ビットになります。したがって、この入力状態を準備できればGroverアルゴリズムが実装するみたいです。やってみます。

グラフ彩色問題のオラクルを作成する

上記の割り当てされていない黒いノードに対して赤、青、オレンジ、緑を制約のもと割り当てます。ここで、色をプログラミングするために整数に置き換えます。

赤=0,
青=1,
オレンジ=2,
緑=3

とします。したがって色を表現するために最大で2量子ビットが必要です($3=11_2$)。以下のように各ノードに２量子ビットを割り当てます。

class PredicateVars(QStruct):
    #QNum[量子ビット数, 符号ありなし, 小数点]
    a: QNum[2, False, 0]
    b: QNum[2, False, 0]
    c: QNum[2, False, 0]
    d: QNum[2, False, 0]
    e: QNum[2, False, 0]

ここで

    QNum[num_qubits, False, digits]

はQmod独自の量子ビット変数で、直接整数エンコードで扱いが可能です。

次に、オラクルを作成します。オラクルの条件は各隣接ノードが同じ色を持ってはいけないという条件です。つまり、命題論理式のように記述が可能です。ここで色変数$C_i$とします。$i$はノード番号です。例えば、$C_A$の意味はノードAの色を表す変数です。

具体的にノードAの場合の制約条件を考えます。ノードAは3つの隣接ノードを持っているので以下の色の制約があります。
$ (C_A \neq 赤) \land (C_A \neq C_B) \land (C_A \neq 青)$
ノードBは
$ (C_B \neq 赤) \land (C_B \neq C_C) \land (C_B \neq C_E) \land (C_B \neq オレンジ)$
のように考えられます。$\land$はandです。

この考えで、単純に各ノードごとの条件をまとめて書くと、オラクルはQmodで以下のように書けます。

def classical_predicate(a, b, c, d, e):
    return logical_and(logical_and(logical_and(logical_and(logical_and(logical_and(logical_and(logical_and(logical_and(logical_and(logical_and(logical_and((a != 0), (a != 1)),(a != b)), (b != 0)),(b != c)),(b != e)), (b != 2)),(c != 0)),(c != d)),(c != e)),(d != 1)),(d != 2)),(e != 3))

@qfunc
def quantum_predicate(vars: PredicateVars, res: QBit):
    res ^= classical_predicate(vars.a, vars.b, vars.c, vars.d, vars.e)

ここで、classical_predicate箇所はqmodによるpython表現で, 直接オラクルを表現できます。quantum_predicateの箇所では, 古典表現したオラクル関数を満たすときに1,0（True, False）をres量子ビットに格納させます。ここで、

logical_and

は$\land$と同じ機能です。

オラクルがかけたので、最後に、

@qfunc
def main(vars: Output[PredicateVars]):
    allocate(vars.size, vars)
    grover_search(
        reps=3,
        oracle=lambda vars: phase_oracle(quantum_predicate, vars),
        packed_vars = vars,
    )

を書いて、グラフ彩色問題を解くGroverアルゴリズムがかけました。reps=3（繰り返し回数）とします。ここで、

allocate(vars.size, vars)

は量子ビット変数varsに対してvarsのサイズの量子ビットを割り当てるという意味みたいです。varsはここで、class PredicateVars(QStruct)の中身の合計の10量子ビットを割り当てています。
作成が完了したので、次に回路を自動生成します。コマンドを以下のようです。

qmod = create_model(main)
qprog = synthesize(qmod)
show(qprog)

得られた回路です。Hゲート, Grover opratorが３つ表示されています。

オラクルの箱の中身を見てみると、先ほど作成したclassiq_predicateのオラクルが自動で作成されています。

実行して答えを確認します。
QiskitのようにExecutionをします。

result = execute(qprog).result_value()

制約を満たす解が増幅されていそうです。

上記の量子回路のシミュレーション結果を古典計算で可視化します。以下のプログラミングを作成しました。

#Post processing
import numpy as np

def convert_color(int_val):
    if int_val == 0:
        return "red"
    elif int_val == 1:
        return "blue"
    elif int_val == 2:
        return "orange"
    elif int_val == 3:
        return "green"
    

NUM_SOLUTIONS = 16
NUM_VARIABLES = 5
solution_list = np.zeros((NUM_SOLUTIONS, NUM_VARIABLES),dtype=int)
for k in range(NUM_SOLUTIONS):
    parsed_result = result.parsed_counts[k].state["vars"]
    a, b, c, d, e = int(parsed_result["a"]), int(parsed_result["b"]), int(parsed_result["c"]), int(parsed_result["d"]), int(parsed_result["e"])
    solution_list[k,0] = a
    solution_list[k,1] = b
    solution_list[k,2] = c
    solution_list[k,3] = d
    solution_list[k,4] = e
    print("a =", a, ", b =", b, ", c =", c, "d =", d, "e =", e, ":", classical_predicate(a, b, c, d, e))

NUM_SOLUTIONSはGroverの最適な計算回数パラメータ$M$ですが、今回はわからないので実験的に探しました。
最後に可視化します。

import networkx as nx
import matplotlib.pyplot as plt

# グラフの作成
G = nx.Graph()

# ノードの追加
nodes = ["A", "B", "C", "D", "E", 0, 1, 2, 3]
G.add_nodes_from(nodes)

# エッジの追加
edges = [
    (0, "A"), (0, "B"), (0, "C"),
    ("A", "B"), ("A", 0), ("A", 1),
    ("B", "C"), ("B", "E"), ("B", 2),
    ("C", "D"), ("C", "E"),
    ("D", 1), ("D", 2),
    ("E", 3)
]
G.add_edges_from(edges)

# 元の座標
pos = {
    "A": (0, 2), "B": (-1, 1), "C": (0, 1), "D": (1, 1), "E": (0, 0),
    0: (-1, 2), 1: (1, 2), 2: (-1, 0), 3: (1, 0)
}


# 16個の異なるグラフを作成する
graphs = []
for i in range(NUM_SOLUTIONS):
    graphs.append(G)
    

fig, axes = plt.subplots(4, 4, figsize=(12, 12))  
# 各サブプロットに描画
for i, ax in enumerate(axes.flat): 
    nx.draw(
    G, pos, ax=ax, with_labels=True, node_size=400, nodelist= ["A", "B", "C", "D", "E"]
        ,node_color=[convert_color(solution_list[i,0]),convert_color(solution_list[i,1]),convert_color(solution_list[i,2]),convert_color(solution_list[i,3]),convert_color(solution_list[i,4])]
        , font_color="white"
)
    nx.draw_networkx_nodes(G, pos, ax=ax, nodelist=[0], node_color="red")
    nx.draw_networkx_nodes(G, pos, ax=ax, nodelist=[1], node_color="blue")
    nx.draw_networkx_nodes(G, pos, ax=ax, nodelist=[2], node_color="orange")
    nx.draw_networkx_nodes(G, pos, ax=ax, nodelist=[3], node_color="green")
    ax.set_title(f"solution {i+1}")

plt.tight_layout()
plt.show()

目視で確認してもあっていそうです。

以上です。