機械学習でのQAOA変分パラメータ推定
量子コンピュータで離散最適化問題を扱うための理論である量子断熱計算をパラメータ化したQAOAについてみてみたいとおもいます。通常QAOAのパラメータはランダムからスタートして収束させますが、参考文献でLSTMを利用して初期パラメータ位置を推定するのがありましたので、近いことをしたいと思います。
参考
Learning to learn with quantum neural networks
QAOA回路を確認する
blueqatではQAOAを自動的に計算してくれるのですが、ほぼブラックボックス化しています。
from blueqat import vqe
from blueqat.pauli import *
#式をハミルトニアンの形に
hamiltonian = -1.0*Z[0] + 1.0*Z[1] + 0.5*Z[0]*Z[1]
#ステップ数を決める
step = 1
#中では変分化したアルゴリズムが動いている
result = vqe.Vqe(vqe.QaoaAnsatz(hamiltonian, step)).run()
print(result.most_common(4))
(((0, 1), 0.757124296858263), ((1, 0), 0.22248154912991072), ((1, 1), 0.010197077005912936), ((0, 0), 0.010197077005912932))
自動的にパラメータ化をして最適化までしてくれるのでQAOAをあまり知らなくても実行できますが、今回はパラメータを最適化したいのでブラックボックス化しないようにします。
QAOAは時間発展演算子を利用して時間発展させます。
from blueqat import Circuit
a = [term.get_time_evolution() for term in hamiltonian]
time_evolution = Circuit()
angle = np.random.rand()*2*np.pi
for evo in a:
evo(time_evolution, angle)
time_evolution.run_with_ibmq(returns='draw', output='mpl')
<Figure size 327.252x204.68 with 1 Axes>
これに横磁場ハミルトニアンも時間発展させてみます。
transverse = X[0] + X[1]
b = [term.get_time_evolution() for term in transverse]
time_evolution_x = Circuit()
angle2 = np.random.rand()*2*np.pi
for evo in b:
evo(time_evolution_x, angle2)
time_evolution_x.run_with_ibmq(returns='draw', output='mpl')
<Figure size 267.052x204.68 with 1 Axes>
この二つをくっつければp=1のQAOAに近づきます。
(time_evolution + time_evolution_x).run_with_ibmq(returns='draw', output='mpl')
<Figure size 507.852x204.68 with 1 Axes>
量子断熱計算の初期状態は横磁場ハミルトニアンに対応した固有状態を取りますので、|+>状態を作って完了です。
(Circuit(2).h[:] + time_evolution + time_evolution_x).run_with_ibmq(returns='draw', output='mpl')
<Figure size 568.052x204.68 with 1 Axes>
今回は角度は適当なランダムな値を入れましたが、このアングルを最適化することによって量子断熱計算を変分計算として最適解に近い計算を導き出せます。
さきほどQAOAで実行した回路を見てみると、
result.circuit
Circuit(2).h[:].cx[0, 1].rz(-16.66276103353925)[1].cx[0, 1].rz(33.3255220670785)[0].rz(-33.3255220670785)[1].rx(4.517772793364023)[:]
このように問題設定のハミルトニアンに対応するパラメータと横磁場ハミルトニアンに対応するパラメータが見つかっています。これらのパラメータは角度をパラメータ化し、scipy minimizeで最適化問題を解いてもとまっています。最小化するのはハミルトニアンと状態ベクトルから得られる期待値です。
こちらをランダムからスタートしてどのように最適化されていくか最適化の過程での値をプロットしていきたいと思います。
変分化するパラメータをbeta, gammaとしてランダムスタートさせます
time_evolution = Circuit(2).h[:]
beta = np.random.rand()*2*np.pi
for evo in a:
evo(time_evolution, beta)
gamma = np.random.rand()*2*np.pi
for evo in b:
evo(time_evolution, gamma)
time_evolution.run_with_ibmq(returns='draw', output='mpl')
<Figure size 568.052x204.68 with 1 Axes>
この回路を実行した後の状態ベクトルは、
sv = time_evolution.run()
sv
array([ 0.30809098-0.44281209j, 0.21083871-0.43441041j,
-0.4292161 -0.02455541j, 0.30809098-0.44281209j])
こちらと、問題のハミルトニアンを使って、期待値を求めます。
#状態ベクトルから期待値を返す関数
def sv_expt(sv_tmp, h_param):
z0 = np.array([[1,0],[0,-1]])
i0 = np.eye(2)
h0 = h_param[0] * np.kron(z0, i0) + h_param[1] * np.kron(i0, z0) + h_param[2] * np.kron(z0, z0)
expt = sv_tmp.conj()@h0@sv_tmp
return expt.real
#試しにパラメータを入れてみます
h_init = [-1, 1, 0.5]
expt = sv_expt(sv, h_init)
print(expt)
-0.014666664541780838
パラメータを最適化していきます。
import copy
arr_expt = [expt] #期待値の履歴
arr_beta = [beta] #角度1の履歴
arr_gamma = [gamma] #角度2の履歴
step = 1000
#回路から期待値を返す関数
def te_circuit(te_a, te_b, angle_beta, angle_gamma, h_param):
te_tmp = Circuit(2).h[:]
for evo in te_a:
evo(te_tmp, angle_beta)
for evo in te_b:
evo(te_tmp, angle_gamma)
sv_tmp = te_tmp.run()
return sv_expt(sv_tmp, h_param)
#最適化のループ
for _ in range(step):
#微分用
h = 0.01
#学習率
e = 0.1
#期待値を求める
expt = te_circuit(a, b, beta, gamma, h_init)
#beta差分での期待値
beta_tmp = copy.copy(beta) + h
expt_delta_beta = te_circuit(a, b, beta_tmp, gamma, h_init)
#gamma差分での期待値
gamma_tmp = copy.copy(gamma) + h
expt_delta_gamma = te_circuit(a, b, beta, gamma_tmp, h_init)
#微分を使って角度を更新
beta -= (expt_delta_beta - expt) * e
gamma -= (expt_delta_gamma - expt) * e
#更新された値をチェック
expt = te_circuit(a, b, beta, gamma, h_init)
arr_beta.append(beta)
arr_gamma.append(gamma)
arr_expt.append(expt)
まず期待値の変遷を見てみます。
import matplotlib.pyplot as plt
plt.plot(arr_expt)
plt.show()
<Figure size 432x288 with 1 Axes>
角度の最適化されるのも見てみたいと思います。
plt.plot(arr_beta, label="beta")
plt.plot(arr_gamma, label="gamma")
plt.legend()
plt.show()
<Figure size 432x288 with 1 Axes>
つぎに今回はp=1でしたので、角度によってどのように極小値に向かうのかを図示してみました。
参考：https://pennylane.ai/qml/demos/learning2learn.html
grid = 30
beta_grid = np.linspace(0, 2*np.pi, grid)
gamma_grid = np.linspace(0, 2*np.pi, grid)
Beta, Gamma = np.meshgrid(beta_grid, gamma_grid)
data = np.zeros((grid, grid))
for i in range(grid):
for j in range(grid):
data[i][j] = te_circuit(a, b, beta_grid[i], gamma_grid[j], h_init)
#③等高線を描画
plt.contourf(Beta, Gamma, data)
plt.colorbar()
#角度を描画
plt.plot(arr_beta, arr_gamma, linestyle="--", color="red")
plt.plot(arr_beta[-1], arr_gamma[-1], linestyle="--", color="red", marker="x")
plt.xlabel(r"$\beta$")
plt.ylabel(r"$\gamma$")
plt.show()
<Figure size 432x288 with 2 Axes>
初期状態から少しずつ色の濃いエネルギーの低い状態へ移動しているのが分かります。
グラフ自体は3次元で見るとわかりやすいですね。
from mpl_toolkits.mplot3d import axes3d
ax = plt.figure().add_subplot(projection='3d')
ax.contour3D(Beta, Gamma, data, 200, cmap='rainbow')
plt.show()
<Figure size 432x288 with 1 Axes>
角度パラメータを数理最適を用いて最適化することでQAOAの解を求める様子が分かりました。次にこのコスト関数の角度の初期値を推定してなるべく解に近いところからスタートするように古典NNを学習させてみたいと思います。まずはデータを作ります。今回は2量子ビットのハミルトニアンを任意で作り、その収束した期待値をラベルとして学習します。
def make_data():
rand1 = np.random.rand()*2-1
rand2 = np.random.rand()*2-1
rand3 = np.random.rand()*2-1
h_make = [rand1, rand2, rand3]
hami1 = rand1*Z[0] + rand2*Z[1] + rand3*Z[0]*Z[1]
transverse2 = X[0] + X[1]
aa = [term.get_time_evolution() for term in hami1]
bb = [term.get_time_evolution() for term in transverse2]
beta = np.random.rand()*2*np.pi
gamma = np.random.rand()*2*np.pi
min1 = 0
min2 = rand1
min3 = rand2
min4 = rand1 + rand2 + rand3
minimum = min(min1, min2, min3, min4)
step = 1000
h = 0.01
e = 0.1
for _ in range(step):
expt = te_circuit(aa, bb, beta, gamma, h_make)
#beta差分での期待値
beta_tmp = copy.copy(beta) + h
expt_delta_beta = te_circuit(aa, bb, beta_tmp, gamma, h_make)
#gamma差分での期待値
gamma_tmp = copy.copy(gamma) + h
expt_delta_gamma = te_circuit(aa, bb, beta, gamma_tmp, h_make)
#微分を使って角度を更新
beta -= (expt_delta_beta - expt) * e
gamma -= (expt_delta_gamma - expt) * e
#更新された値をチェック
expt = te_circuit(a, b, beta, gamma, h_init)
if expt < minimum+0.5:
return [rand1, rand2, rand3], [beta, gamma]
else:
return 0
学習データを作り、訓練用と評価用に分けます。
from sklearn.neural_network import MLPRegressor
from sklearn.model_selection import train_test_split
X_data = []
y_data = []
for _ in range(1000):
get_data = make_data()
if get_data != 0:
X_data.append(get_data[0])
y_data.append(get_data[1])
print("*", end="")
else:
print("x", end="")
print("")
print(len(X_data))
X_train, X_test, y_train, y_test = train_test_split(X_data, y_data, test_size=0.2, random_state=0)
xx******xxxxxx******x*x*x*xx***x**xx***x*****x**xx**x**x*xx**xxx***x****x*x****x*x*xxxx**x*xx*x****x*xxx***xxx**xx*x*xx*xxx*x****x**xx*x**xxx*xx*x*x*x*xxx*x*xx**xxxxxxxxx*xx*x*x*x**x**x*x*xx**x***x***xxxxxx**xxxxxx*xxxx**x*xx*xx*x*xx**xx*x*x*xxxx**xxxxx***x*x*x*****xx***x***xx*x*xxxx*x*xxx**xxx*****x**xxxx***x******xxxx****x*xxxxx*xxxx*x***xx*x*xxx**xx**xx***x*xxx*xx*xxxxxx*x****x*xx***x*x***xx**x*x***x*x***x***xx**xxxx*x*xx****xxx*x****xx*xx**xx*xx*x*x*x**x*xx*xxx**xxx***xx****xx***xxxx*xxx*x**x*xxx*xx*xxxxx**xxx*xx****x**x*x**x*xxx*****x*xxxxxxx*xx**x**xxx*xx*x**x*xx*x*xx*x**xx**x*x*x*x**x**xxx**xx*******x*x*xx*xxxxx*x*x*xxxxxxx**x*xx*xxx**xx*xx*xx**x**x*xx*x*x*xxx*x***x*xxx***x**x***x***xxx*****xx*x**x***xx*xxx*x*x**xxxxx**xx*xx*xxx*xx*****x**x**x*x******xxx*xxxx*x*x*x*xxxxxxxxx*xx********xxxx*x*x*x*x**x*xx**x***x*xx*xx****x**xx**xxx*x*xxxxxxxxxxxxx***xx***x**x*x*x**xx**x**x*x**x*xxxx***xxx***x*xxxx**xx****x***xx*x**x**x*xx*xx*x**x*x*xxxx**x***xxxxx*xxx**x**xx*xx**x***x*xxxx***xxxxx
483
学習させたモデルを評価します。
model_nn = MLPRegressor(random_state=42)
model_nn.fit(X_train, y_train)
score = model_nn.score(X_train, y_train)
print("Training score:", score)
score2 = model_nn.score(X_test, y_test)
print("Valuation score:", score2)
Training score: 0.01067089406467775
Valuation score: -0.020423018316892194
/opt/conda/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:614: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (200) reached and the optimization hasn't converged yet.
warnings.warn(
実際に新しいハミルトニアンから変分パラメータを推定してみます。
#ランダムな係数からハミルトニアンを作成
rand1_pred = np.random.rand()*2-1
rand2_pred = np.random.rand()*2-1
rand3_pred = np.random.rand()*2-1
h_pred = [rand1_pred, rand2_pred, rand3_pred]
#角度を推定する
x_pred = model_nn.predict([h_pred])
x_pred = x_pred[0]
print(x_pred)
beta = x_pred[0]
gamma = x_pred[1]
#ハミルトニアンから時間発展を作成
hami_pred = rand1_pred*Z[0] + rand2_pred*Z[1] + rand3_pred*Z[0]*Z[1]
transverse = X[0] + X[1]
a_pred = [term.get_time_evolution() for term in hami_pred]
b_pred = [term.get_time_evolution() for term in transverse]
expt_pred = te_circuit(a_pred, b_pred, beta, gamma, h_pred)
print(expt_pred)
arr_expt = [expt_pred] #期待値の履歴
arr_beta = [beta] #角度1の履歴
arr_gamma = [gamma] #角度2の履歴
step = 1000
h = 0.01
e = 0.1
for _ in range(step):
expt = te_circuit(a_pred, b_pred, beta, gamma, h_pred)
#beta差分での期待値
beta_tmp = copy.copy(beta) + h
expt_delta_beta = te_circuit(a_pred, b_pred, beta_tmp, gamma, h_pred)
#gamma差分での期待値
gamma_tmp = copy.copy(gamma) + h
expt_delta_gamma = te_circuit(a_pred, b_pred, beta, gamma_tmp, h_pred)
#微分を使って角度を更新
beta -= (expt_delta_beta - expt) * e
gamma -= (expt_delta_gamma - expt) * e
#更新された値をチェック
expt = te_circuit(a_pred, b_pred, beta, gamma, h_pred)
arr_beta.append(beta)
arr_gamma.append(gamma)
arr_expt.append(expt)
[3.19008567 3.48471883]
-0.1294131295898297
次にランダムでも生成してみます。
beta = np.random.rand()*2*np.pi
gamma = np.random.rand()*2*np.pi
expt_pred = te_circuit(a_pred, b_pred, beta, gamma, h_pred)
print(expt_pred)
arr_expt2 = [expt_pred] #期待値の履歴
arr_beta2 = [beta] #角度1の履歴
arr_gamma2 = [gamma] #角度2の履歴
for _ in range(step):
expt = te_circuit(a_pred, b_pred, beta, gamma, h_pred)
#beta差分での期待値
beta_tmp = copy.copy(beta) + h
expt_delta_beta = te_circuit(a_pred, b_pred, beta_tmp, gamma, h_pred)
#gamma差分での期待値
gamma_tmp = copy.copy(gamma) + h
expt_delta_gamma = te_circuit(a_pred, b_pred, beta, gamma_tmp, h_pred)
#微分を使って角度を更新
beta -= (expt_delta_beta - expt) * e
gamma -= (expt_delta_gamma - expt) * e
#更新された値をチェック
expt = te_circuit(a_pred, b_pred, beta, gamma, h_pred)
arr_beta2.append(beta)
arr_gamma2.append(gamma)
arr_expt2.append(expt)
0.013404047639859196
ランダムでの生成との収束を比較してみます
import matplotlib.pyplot as plt
plt.plot(arr_expt, label="from NN")
plt.plot(arr_expt2, label="random")
plt.legend()
plt.show()
<Figure size 432x288 with 1 Axes>
収束についてランドスケープ上で軌跡を比較してみます。
grid = 30
beta_grid = np.linspace(0, 2*np.pi, grid)
gamma_grid = np.linspace(0, 2*np.pi, grid)
Beta, Gamma = np.meshgrid(beta_grid, gamma_grid)
data = np.zeros((grid, grid))
for i in range(grid):
for j in range(grid):
data[i][j] = te_circuit(a_pred, b_pred, beta_grid[i], gamma_grid[j], h_pred)
#③等高線を描画
plt.contourf(Beta, Gamma, data)
plt.colorbar()
#角度を描画
plt.plot(arr_beta, arr_gamma, linestyle="--", color="red")
plt.plot(arr_beta[-1], arr_gamma[-1], linestyle="--", color="red", marker="x")
#角度を描画2
plt.plot(arr_beta2, arr_gamma2, linestyle="--", color="blue")
plt.plot(arr_beta2[-1], arr_gamma2[-1], linestyle="--", color="blue", marker="x")
plt.xlabel(r"$\beta$")
plt.ylabel(r"$\gamma$")
plt.show()
<Figure size 432x288 with 2 Axes>
from mpl_toolkits.mplot3d import axes3d
ax = plt.figure().add_subplot(projection='3d')
ax.contour3D(Beta, Gamma, data, 200, cmap='rainbow')
plt.show()
<Figure size 432x288 with 1 Axes>
何度か試してみましたが、なんとなく気休めかもしれませんが、ある程度推定された方のパラメータのほうが全体的には収束しやすかった気がします。たまたまランダムのほうがいい値を取ることもあります。以上です。