🔖

Qiskitを使ってQAOAの経路最適化を試す

2025/08/20に公開

個人的な興味から「IBM Quantumで学ぶ量子コンピュータ」（秀和システム)の勉強を進めている。

アルゴリズムの理論面の理解はわかりやすく参考になる部分も大いにあるが、
実装をしているとライブラリのバージョンの関係でつまずいてしまう部分も多い。
これは著者の方もご自身のお仕事の関連サイトで書かれていて、ある程度のバージョンまではカバーしておられたようですが、永遠にはできないし無理もありません（なんなら技術書の宿命ですよね…）。

今回は書籍中の経路最適化をなんとか実装してみたのでそのメモとして書いておきます。

問題設定

最小化したい関数は次の通り（変数は q0〜q3 の4つ）：

f(q) = 4q_0 q_1 + 4q_0 q_2 + 8q_1 q_2 + 2q_1 q_3 + 2q_2 q_3 + 4q_0 + 4q_1 + 4q_2 + 4q_3

さらに制約として：

$q_0$ と $q_1$ は排他的に 1 つだけが 1
$q_2$ と $q_3$ も排他的に 1 つだけが 1

この制約を守るために、QAOA のミキサーとして 1/2(XX+YY) を使う。

実装

私は25/08中旬時点でのGoogle Colaboratoryで走らせています。
まずライブラリをインストールします。

pip install "qiskit>=1.0,<2" "qiskit-aer>=0.14,<0.15" "qiskit-algorithms>=0.3,<0.4" "qiskit-optimization>=0.6,<0.7"

このQiskitが結構クセモノで、ライブラリのバージョンごとに結構挙動が違って戸惑いました。

ライブラリのインポートをします。

from qiskit_optimization import QuadraticProgram
from qiskit.quantum_info import SparsePauliOp
from qiskit import QuantumCircuit
from qiskit.primitives import Sampler  # V1 Sampler（QAOAと互換性あり）
from qiskit_algorithms import QAOA
from qiskit_algorithms.optimizers import COBYLA
import matplotlib.pyplot as plt

1. QUBO定義とIsing変換

QUBOを定義した上でto_isingでイジング模型に変換します。
このQUBOの定義については書籍の内容に従うので省略。

qp = QuadraticProgram("pairwise-onehot")
for i in range(4):
    qp.binary_var(f"q{i}")

qp.minimize(
    linear=[4, 4, 4, 4],
    quadratic={(0,1):4, (0,2):4, (1,2):8, (1,3):2, (2,3):2}
)

cost_op, offset = qp.to_ising()

2. XYミキサー定義

今回はビット0,1と2,3の組み合わせごとにどちらかが0,どちらかが1の値を取るのでXYミキサーで得られる解の範囲を限定させることができます。

def xy_mixer_operator(num_qubits, pairs):
    labels, coeffs = [], []
    for (i, j) in pairs:
        lab = ['I'] * num_qubits
        lab[i] = lab[j] = 'X'
        labels.append(''.join(lab[::-1]))
        coeffs.append(0.5)
        lab = ['I'] * num_qubits
        lab[i] = lab[j] = 'Y'
        labels.append(''.join(lab[::-1]))
        coeffs.append(0.5)
    return SparsePauliOp.from_list(list(zip(labels, coeffs)))

n = 4
pairs = [(0,1), (2,3)]
mixer_op = xy_mixer_operator(n, pairs)

3. 可行初期状態

init = QuantumCircuit(n)
init.x(0); init.x(2)   # 例: |1010> （制約を満たす）

4. QAOA の実行

sampler = Sampler()
qaoa = QAOA(sampler=sampler, optimizer=COBYLA(maxiter=200), reps=3,
            mixer=mixer_op, initial_state=init)

result = qaoa.compute_minimum_eigenvalue(cost_op)
print("最良期待値:", result.eigenvalue.real + offset)

5. サンプリングと分布の可視化

実際は下のコードのqdで解を得ていますが、せっかくなのでこれをグラフィカルに可視化させてみます。

opt_circ = qaoa.ansatz.assign_parameters(result.optimal_point)
qd = sampler.run([opt_circ], shots=4096).result().quasi_dists[0]

def plot_qaoa_distribution(qd, n=4, order="q0..q3"):
    def label_of(k):
        s = format(k, f"0{n}b") if isinstance(k, int) else k
        return s[::-1] if order=="q0..q3" else s

    items = [(label_of(k), float(p)) for k,p in qd.items() if p>0]
    items.sort(key=lambda x: x[1], reverse=True)
    labels, probs = zip(*items)

    plt.figure(figsize=(max(6,len(items)*0.8),4))
    plt.bar(range(len(items)), probs)
    plt.xticks(range(len(items)), labels)
    plt.ylabel("Probability")
    plt.xlabel(f"bitstring ({order})")
    plt.show()
    return items[0]

best = plot_qaoa_distribution(qd, n=4)
print("最頻ビット列:", best)

実行結果の例

最良期待値は ~8 前後。
分布を見ると、q0=1, q1=0, q2=0, q3=1（つまり "1001"）が最頻ビット列として浮かび上がる。
制約どおり各ペアで 1 を１つにさせることができました。

ポイント

QuadraticProgram で QUBO を定義 → to_ising で Ising Hamiltonian に変換
QAOA 実行後は qd（quasi distribution）を可視化して最適解候補を確認させた

量子コンピューティング分野に関心があったので、勉強もかねてこの本を読んでみましたが、まだまだ勉強が必要そうです。
新しい書籍もいくつか出ているので、勉強しつつ最新の研究もフォローしていけたらと思います。