量子ゲートブログ5:量子ゲートブログ1-4の技術を使い、量子アニーリング勢がQAOAに挑戦！Max-Cut問題を3つのアプローチで解き比べてみた

今回のブログは数式が多いので、基本lateXのコードで記号なども書いています。ご了承ください。

こんにちは！普段は量子アニーリングで組合せ最適化問題を解いていますが、今回は量子ゲート方式の代表的なアルゴリズムであるQAOA (Quantum Approximate Optimization Algorithm) に挑戦してみました。

題材は、量子計算の定番「最大カット問題（Max-Cut問題）」です。

今回は、実装の難易度や自由度を変えて、以下の3通りのパターンでコードを書いてみました。

ライブラリの自動作成機能を使ってハミルトニアンをQAOAに投げる
ハミルトニアンは自作して、QAOAの実行はライブラリに任せる
ハミルトニアンもQAOA回路もすべて自作する

単純にMax-Cut問題を解くと言っても、今回は量子アニーリング向けに開発された教科書の問題を量子ゲートで解きます(結局同じMax-Cut問題なので今までの変わりないのですが量子アニーリングとゲート両方で解けるという証明として)。

今回使用させていただいたのは日本量子コンピュティング協会(JQCA)が量子アニーリング検定で使用している「tytan」パッケージを開発された安田さん(QFのメンバーです）のgithubの教科書です。(湊さんと共同でtytanパッケージを作成されました)

https://github.com/tytansdk/tytan_tutorial

今回使用している教科書:

https://colab.research.google.com/drive/1Q39CJyYsUR3bbwWleLRgHdHUilgpO2KB?usp=sharing

1. Max-Cut問題とその定式化

まずは基本のおさらいです。Max-Cut問題とは、グラフの頂点集合を2つに分けたとき、その間をまたぐ「エッジの数」を最大にする問題です。

コスト関数（組合せ最適化）

頂点 $i$ が属するグループを $x_i \in \{0, 1\}$ とすると、最大化したいカット数 $C(x)$ は以下のように表せます。

$$C(x) = \sum_{(i,j) \in E} (x_i + x_j – 2x_ix_j)$$

イジング模型への変換

量子計算で扱いやすいように、変数を $x_i \to s_i \in \{-1, 1\}$ （スピン）に変換します。このとき、最大化問題から最小化問題（エネルギー最小化）に書き換えたハミルトニアン $H$ は以下のようになります。

$$H = \frac{1}{2} \sum_{(i,j) \in E} (Z_i Z_j – 1)$$

※ $Z_i$ はパウリZ演算子です。この $H$ の期待値を最小化することがQAOAの目的になります。

2. 実装パターン1：全自動アプローチ

最近のライブラリ（QiskitやcuTensorなど）は非常に優秀です。グラフを定義するだけで、内部で自動的にハミルトニアンを作り、QAOAのサーキットまで組み立ててくれます。

ポイント: アルゴリズムの構造を意識せず、まずは結果を見たい時に最適です。

#link = https://colab.research.google.com/drive/1Q39CJyYsUR3bbwWleLRgHdHUilgpO2KB?usp=sharing

qp = QuadraticProgram("tytan第3世代最大カット問題を量子ゲートで解く！")
edges = [(0, 1), (0, 2), (1, 3), (2, 3), (2, 4), (3, 4)]

for i in range(5):
    qp.binary_var(f"x{i}")

linear_terms = {}
quadratic_terms = {}

for i, j in edges:
    linear_terms[f"x{i}"] = linear_terms.get(f"x{i}", 0) + 1
    linear_terms[f"x{j}"] = linear_terms.get(f"x{j}", 0) + 1
    quadratic_terms[f"x{i}", f"x{j}"] = -2

qp.maximize(linear = linear_terms, quadratic = quadratic_terms)

print(qp.prettyprint())

sampler = Sampler()
optimizer = COBYLA(maxiter = 100)

qaoa = QAOA(sampler=sampler, optimizer=optimizer, reps = 2)

optimizer_wrapper = MinimumEigenOptimizer(qaoa)
result = optimizer_wrapper.solve(qp)

print(result.x)
print(result.fval)

3. 実装パターン2：ハミルトニアン自作＋ QAOAライブラリ

「問題の定義までは自分でやりたいけれど、QAOAの最適化ループ（古典計算とのやり取り）は既存のものを使いたい」というパターンです。

ポイント: 複雑な制約条件がある問題など、自分でイジングモデルを設計する必要がある実戦的な形式です。

Python

#tytan第3世代最大カット問題を量子ゲートで解く！
#link = https://colab.research.google.com/drive/1Q39CJyYsUR3bbwWleLRgHdHUilgpO2KB?usp=sharing

import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit import ParameterVector
from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector

from qiskit.primitives import StatevectorEstimator
from qiskit_algorithms import VQE
from qiskit_algorithms.optimizers import COBYLA


n = 5
edges = [(0, 1), (0, 2), (1, 3), (2, 3), (2, 4), (3, 4)]

pauli_list = []

for i, j in edges:
    pauli_str_list = ["I"]*n
    pauli_str_list[i] = "Z"
    pauli_str_list[j] = "Z"

    pauli_str_list.reverse()
    pauli_str = "".join(pauli_str_list)
    pauli_list.append((pauli_str, 1.0))

H = SparsePauliOp.from_list(pauli_list)

print(H)

hamiltonian = SparsePauliOp.from_list(pauli_list)
print("--- 構築されたハミルトニアン ---")
print(hamiltonian)

sampler = Sampler()
optimizer = COBYLA(maxiter=1000)
qaoa = QAOA(sampler=sampler, optimizer=optimizer, reps=3)
result = qaoa.compute_minimum_eigenvalue(hamiltonian)

print("\n最適解のビット列:", result.best_measurement['bitstring'])
print("最小エネルギー:", result.eigenvalue)

4. 実装パターン3：フルスクラッチ（全部自作）

最後は、ハミルトニアンの定義から、時間発展演算子、QAOAの量子サーキット、パラメータ最適化までをすべて自分で記述するパターンです。

ポイント: QAOAの内部で何が起きているか（コスト層とミキサー層の交互作用など）を完全に理解するために最も勉強になる方法です。

Python

#tytan第3世代最大カット問題を量子ゲートで解く！
#link = https://colab.research.google.com/drive/1Q39CJyYsUR3bbwWleLRgHdHUilgpO2KB?usp=sharing

import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit import ParameterVector
from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector

from qiskit.primitives import StatevectorEstimator
from qiskit_algorithms import VQE
from qiskit_algorithms.optimizers import COBYLA


n = 5
edges = [(0, 1), (0, 2), (1, 3), (2, 3), (2, 4), (3, 4)]

pauli_list = []
coeffs = []

for i, j in edges:
    pauli_str_list = ["I"]*n
    pauli_str_list[i] = "Z"
    pauli_str_list[j] = "Z"

    pauli_str_list.reverse()
    pauli_str = "".join(pauli_str_list)
    pauli_list.append((pauli_str, 1.0))

H = SparsePauliOp.from_list(pauli_list)

print(H)

p=2
gammas = ParameterVector("gamma", p)
betas = ParameterVector("beta", p)

qc = QuantumCircuit(n)

for i in range(n):
    qc.h(i)

for i in range(p):
    for k, j in edges:
        qc.cx(k, j)
        qc.rz(2*gammas[i], j)
        qc.cx(k, j)

    for t in range(n):
        qc.rx(2*betas[i], t)

print("\nQAOA ansatz circuit (p=2):")
print(qc.draw("text"))


##擬似量子コンピュータで計算

estimator = StatevectorEstimator()
optimizer = COBYLA(maxiter=100)

vqe = VQE(
    estimator=estimator,
    ansatz=qc,
    optimizer=optimizer,
)

result = vqe.compute_minimum_eigenvalue(H)

print("\n=== VQE (hand-made QAOA ansatz) result ===")
print("最小エネルギー:", result.eigenvalue)
print("最適パラメータ (gamma[0], gamma[1], beta[0], beta[1]):")
print(result.optimal_parameters)

# 1) 最適パラメータを回路に代入して状態ベクトルを作る
bound = qc.assign_parameters(result.optimal_parameters, inplace=False)
state = Statevector.from_instruction(bound)


# 2) 各計算基底状態の確率を計算
probs = state.probabilities()  # 長さ 2^4 = 16 の配列

# 最も確率の高い基底状態のインデックス
best_index = int(np.argmax(probs))

# 3) インデックス → 各量子ビットのビット値（q0〜q3）に変換
#    bits[0] が qubit0(=頂点1), bits[1] が qubit1(=頂点2) ... という対応
bits = [ (best_index >> i) & 1 for i in range(n) ]  # i=0がq0

print("\n=== Classical solution from QAOA state ===")
print("最も確率の高い基底状態 index:", best_index)
print("各qubitのビット列 [q0,q1,q2,q3] :", bits)

# 4) このビット列を「頂点のグループ分け」に解釈
#    ここでは 0 = グループA, 1 = グループB とする
group_A = [i+1 for i, b in enumerate(bits) if b == 0]  # 頂点番号は1始まりにして出力
group_B = [i+1 for i, b in enumerate(bits) if b == 1]

print("グループA (bit=0): 頂点", group_A)
print("グループB (bit=1): 頂点", group_B)

# 5) この切り方で何本の辺がカットされているかも計算
def cut_size(bits, edges):
    """bits: 長さ4の 0/1 リスト, edges: (i,j) のリスト"""
    cnt = 0
    for i, j in edges:
        if bits[i] != bits[j]:
            cnt += 1
    return cnt

cut = cut_size(bits, edges)
print("この切り方で切れている辺の本数 (cut size):", cut)

まとめ：やってみた感想

3通りの方法を試してみて、以下のことが分かりました。

手軽さ重視ならパターン1。ただ、局所解に陥りやすいので注意。
柔軟な問題設定をしたいならパターン2。比較的最適解を求められる。
アルゴリズムの研究やチューニングをしたいならパターン3。

アニーリングでは「いかにQUBOを作るか」に集中していましたが、ゲート方式では「いかに浅い回路で効率よく探索するか」というサーキット設計の視点が必要なのが新鮮でした。

皆さんも、自分の理解度に合わせて最適なアプローチを選んでみてください！

次(QAOAの勉強が終わったら)は、今回得られたパラメータ $\gamma, \beta$ の挙動を可視化して、アニーリングのスケジュールとの比較をしてみようと思います！