qiskitで愚直にQAOAをシミュレートしてみました．

QAOAとは

QAOAでは，次のようなユニタリ回路を作用させることでコストハミルトニアン$H_C$の基底状態を求めることを考えます．
$$
\left|\psi_p(\vec{\gamma}, \vec{\beta})\right\rangle=\left(\prod_{k=1}^p\left[e^{-i H_X \beta_k} e^{-i H_C \gamma_k}\right]\right)\left|\psi_{\mathrm{ref}}\right\rangle
$$
ここで，$H_X$は通常$\sum_i X_i$と取ります．また，初期状態も通常
$$
\left|\psi_{\mathrm{ref}}\right\rangle=|+\rangle^{\otimes n}
$$
と取ります．生成した状態$\left|\psi_p(\vec{\gamma}, \vec{\beta})\right\rangle$を用いて期待値$\langle H_C\rangle$を測定し，この値が小さくなるように古典コンピュータで$(\vec{\gamma}, \vec{\beta})$をアップデートしていきます．

設定

今回は論文¹に記載されている次のコストハミルトニアンを考えます²．
$$
H_C= \sum_{(i, j)\in E} w_{i, j}Z_i Z_j
$$
ここで，論文¹同様，グラフ$G=(V,E)$の頂点を$n=6$とし，次数が$D=3,5$の２パターンを考えます．また，$[0,1]$の中から一様ランダムにグラフのウェイト$w_{i,j}$を決定します．20種類のグラフをランダムに生成し，最適化の結果については，その平均と95%信頼区間を見ます．最適化手法はNelder-Meadを用います．

実装

まず，グラフを生成します．

import numpy as np

num_spin = 6

# 完全グラフと３正則グラフ
regular_5 = "regular_5"
regular_3 = "regular_3"

# 利用可能なグラフのリスト
kinds_graph = [regular_5, regular_3]

# 選択されたグラフの種類
chosen_graph = kinds_graph[0]

def make_graph(graph_name):
    # 重みを生成する関数
    def weight_func():
        return np.random.uniform(0, 1)

    # 考えるグラフに応じた結合とその重みの辞書
    weights = {}

    # グラフの種類ごとに結合を生成
    if graph_name == regular_5:
        # 完全グラフの場合
        for i in range(num_spin):
            for j in range(i + 1, num_spin):
                weights[(i, j)] = weight_func()
    elif graph_name == regular_3:
        # 3正則グラフの場合
        for i in range(num_spin):
            # 辺の生成条件に基づいて重みを設定
            if i < (i + 1) % num_spin:
                weights[(i, (i + 1) % num_spin)] = weight_func()
            if i < (i - 1) % num_spin:
                weights[(i, (i - 1) % num_spin)] = weight_func()
            if i < (i + num_spin // 2) % num_spin:
                weights[(i, (i + num_spin // 2) % num_spin)] = weight_func()

    print("edge : ", weights)

    return weights  # 生成された結合と重みの辞書を返す

次に，コストハミルトニアン$H_C$を実装します．ここで，以前にも述べた通り，$|00110\rangle_{43210}$というような順番でハミルトニアンを表現していることに注意してください．

from qiskit.quantum_info import SparsePauliOp

class QAOA:
    def __init__(self, num_spin, weights):
        """
        QAOAの初期化メソッド

        Args:
            num_spin (int): スピンの数
            weights (dict): グラフの重み付きエッジ
        """
        self.num_spin = num_spin
        self.weights = weights

    def cost_hamiltonian(self):
        """
        コストハミルトニアンを生成するメソッド

        Returns:
            hamiltonian (SparsePauliOp): コストハミルトニアン
        """
        # 期待値測定のための問題ハミルトニアンの生成
        operator_in_cost_hamiltonian = []

        for i, j in self.weights.keys():
            temp = ""
            for k in range(self.num_spin):
                if k == self.num_spin - i - 1 or k == self.num_spin - j - 1:
                    temp += "Z"
                else:
                    temp += "I"
            operator_in_cost_hamiltonian.append(temp)

        # SparsePauliOpを使ってハミルトニアンを構築
        hamiltonian = SparsePauliOp(operator_in_cost_hamiltonian, list(self.weights.values()))

        return hamiltonian

次に期待値を測定する関数も，このクラスの中に入れます．


class QAOA:
    # (前略)

    def expectation_hamiltonian(self, params):
        """
        ハミルトニアンの期待値を計算する関数

        Args:
            params (list): QAOAのパラメータ
            
        Returns:
            energy (float): コストハミルトニアンの期待値
        """
        qc = QuantumCircuit(self.num_spin)  # N量子ビットの量子回路を作成
        qc.h(range(self.num_spin))  # すべての量子ビットにアダマールゲートを適用 

        gamma = params[::2]
        beta = params[1::2]
        num_round = len(gamma)

        # 回路の構築
        for i in range(num_round):
            for qubits, coeff in self.weights.items():
                qc.rzz(2 * gamma[i] * coeff, qubits[0], qubits[1])  # γの計算: RZZゲートを適用

            qc.rx(2 * beta[i], range(self.num_spin))  # βの計算: Xの和を適用

        # ハミルトニアンの期待値を計算
        energy = Statevector(qc).expectation_value(self.cost_hamiltonian())

        return np.real(energy)

最後に，このクラスに古典の最適化を行う関数も入れておきます．

from scipy.optimize import minimize

class QAOA:
    # (前略)

    def classical_optimization(self, params):
        """
        古典最適化を実行する関数

        Args:
            params (list): QAOAの初期パラメータ

        Returns:
            res (scipy.optimize.OptimizeResult): 最適化の結果
        """
        
        # gammaとbetaを最適化
        res = minimize(self.expectation_hamiltonian, params, method="Nelder-Mead")

        return res

これで設定は終わりです．実際に最適化してみましょう．

from tqdm import tqdm 
from scipy.linalg import eigh

# サンプリング回数とスピンの数を定義
num_sampling = 20

# サンプリングごとに保存するエネルギーのリスト
sampled_normalized_optimized_energy = []

# サンプリングのループ
for _ in tqdm(range(num_sampling)):

    # 最適化するパラメータの初期値
    params = [0.01, 0.0]

    # グラフの生成
    weights = make_graph(chosen_graph)

    # QAOAインスタンスの作成
    QAOA_instance = QAOA(num_spin, weights)

    # 最大の数の繰り返し回数
    max_round = 9
    normalized_optimized_energy = []

    # エネルギーの最小値と最大値の計算
    eigenvalues, eigenvectors = eigh(QAOA_instance.cost_hamiltonian().to_matrix())
    min_energy = min(eigenvalues)
    max_energy = max(eigenvalues)

    # ラウンドごとの最適化のループ
    for num_round in range(1, max_round):

        gamma = params[::2]
        beta = params[1::2]

        # 最適化を実行
        res = QAOA_instance.classical_optimization(params)

        # 最適化結果を更新
        params[:2 * num_round] = res.x

        # 新しいパラメータを追加
        params.append(0.01)
        params.append(0.0)

        # 正規化された最適化エネルギーを計算して保存
        normalized_optimized_energy.append(np.abs(res.fun - min_energy) / (max_energy - min_energy))

    # 1回のサンプリング結果を保存
    sampled_normalized_optimized_energy.append(normalized_optimized_energy)

ここで，追加する$\gamma$の初期値が$0.01$なのは，論文¹に書かれていたためです．また，通常のQAOAと異なり，ここではラウンド数$P$で最適化されたパラメータを用いて，ラウンド数$P+1$での初期パラメータを決定しています（論文¹がおそらくそうしているため）．

以上の結果をプロットすると，以下のようになります（次数が$3$の場合）

．

論文¹と似た結果になっていることが分かります．

最後に

次回はADAPT-QAOAとの比較も行ってみようと思います．

Zhu, Linghua, et al. "Adaptive quantum approximate optimization algorithm for solving combinatorial problems on a quantum computer." Physical Review Research 4.3 (2022): 033029. ↩ ↩² ↩³ ↩⁴ ↩⁵
この論文自体はADAPT-QAOAの提案論文であり，ADAPT-QAOAと通常のQAOAの性能を比較している論文です． ↩