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})$をアップデートしていきます．
設定
今回は論文1に記載されている次のコストハミルトニアンを考えます2．
$$
H_C= \sum_{(i, j)\in E} w_{i, j}Z_i Z_j
$$
ここで，論文1同様，グラフ$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$なのは，論文1に書かれていたためです．また，通常のQAOAと異なり，ここではラウンド数$P$で最適化されたパラメータを用いて，ラウンド数$P+1$での初期パラメータを決定しています（論文1がおそらくそうしているため）．
以上の結果をプロットすると，以下のようになります（次数が$3$の場合）
．
論文1と似た結果になっていることが分かります．
最後に
次回はADAPT-QAOAとの比較も行ってみようと思います．