この記事について
AWS Braket が公開されてから、一般の人にとっても量子コンピュータの実機に触れる機会が増え、それに伴い、私も量子コンピュータのお勉強をしています。
Braketではゲート型の最適化を行うにあたってQAOAというアルゴリズムを使用します。
AWS Braketのサンプル上ではGraph Coloring Problem という問題を解いており、「各ノードに赤や青の色をつけ、色の異なるノード間の各枝に点数をつけていきより良い点をしよう」という問題の求解を行っています。
こちらの問題ですが、以下の式で表されるコスト関数を定義し、（$z_n= {-1,1}$, $J$は枝間の重み）それを最小化する問題に帰着できます。
H_{C}=\sum_{i>j} J_{i,j} z_{i} z_{j},
このような問題に帰着できるものであれば、AWS Braketのサンプルを流用することで動かすことができますが、最適化問題をより一般的に表すには、以下の通り線形項も加味する必要があります。
H_{C}=\sum_{i>j} J_{i,j} z_{i} z_{j} + \sum_{i=1}^{n} h_{i} z_{i}
AWSのサンプルをそのまま動かすだけでは上記の問題を正確に解くことができません。(近似アルゴリズムなので別のアプローチでももちろん正確ではないが、そのまま使うよりは正確です。)
なのでこういった最適化問題を解くことができるような実装をAWS Braket上で行いました。
以前の記事では、AWS BraketのQAOAのサンプル実装の説明を行っていったのでそれを活用して、今回は少し応用的に実装していこうと思います。
まずはAWS Braket の解説を見たいという方は上記の記事を見てみてください。
また、量子コンピュータ関係の他の記事は、下記で紹介しています。
問題定義
今回扱う問題は以下のような形式になっている必要があります。
H_{C}=\sum_{i>j} J_{i,j} z_{i} z_{j} + \sum_{i=1}^{n} h_{i} z_{i}
最低一つの$i$において$h_i > 0$となっていないと今回やりたいことができません。
なのでとりあえずそれを実現するために何かしらの問題を定義します。
こういった形の問題は探せば結構あります。例えば巡回セールスマン問題や重みつきMax Cut問題等があります。
あまり難しいことはやりたくないので今回は直感的に理解しやすい巡回セールスマン問題を扱っていきます。
また、巡回セールスマン問題は0, 1バイナリ変数で扱うことが多いので今回も0, 1バイナリにマッピングしながらやっていこうと思います。
あくまで今回の焦点はQAOAなので問題設定は既存で定義されているほかの方の問題を採用させていただきたいと思います。
こちらの
D-Waveとblueqat.optで巡回セールスマン問題とmaxcut問題、クリーク問題、自然数分割問題
という記事に記載の問題と全く同じ問題を使用させていただきます。
なので説明は割愛させていただきますが、結果として以下のような行列で示すことができるのでそれを使っていこうと思います。
詳しくは上記の記事を参考にしてみてください。
J1 = np.array([
[-2,2,2,2,2,0,0,0,2,0,0,0,2,0,0,0],
[0,-2,2,2,0,2,0,0,0,2,0,0,0,2,0,0],
[0,0,-2,2,0,0,2,0,0,0,2,0,0,0,2,0],
[0,0,0,-2,0,0,0,2,0,0,0,2,0,0,0,2],
[0,0,0,0,-2,2,2,2,2,0,0,0,2,0,0,0],
[0,0,0,0,0,-2,2,2,0,2,0,0,0,2,0,0],
[0,0,0,0,0,0,-2,2,0,0,2,0,0,0,2,0],
[0,0,0,0,0,0,0,-2,0,0,0,2,0,0,0,2],
[0,0,0,0,0,0,0,0,-2,2,2,2,2,0,0,0],
[0,0,0,0,0,0,0,0,0,-2,2,2,0,2,0,0],
[0,0,0,0,0,0,0,0,0,0,-2,2,0,0,2,0],
[0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,2],
[0,0,0,0,0,0,0,0,0,0,0,0,-2,2,2,2],
[0,0,0,0,0,0,0,0,0,0,0,0,0,-2,2,2],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,2],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2],
])
J2 = np.array([
[0,0,0,0,0,2,0,2,0,1,0,1,0,3,0,3],
[0,0,0,0,2,0,2,0,1,0,1,0,3,0,3,0],
[0,0,0,0,0,2,0,2,0,1,0,1,0,3,0,3],
[0,0,0,0,2,0,2,0,1,0,1,0,3,0,3,0],
[0,0,0,0,0,0,0,0,0,4,0,4,0,2,0,2],
[0,0,0,0,0,0,0,0,4,0,4,0,2,0,2,0],
[0,0,0,0,0,0,0,0,0,4,0,4,0,2,0,2],
[0,0,0,0,0,0,0,0,4,0,4,0,2,0,2,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2],
[0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2],
[0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
])
J = J1 + 0.2*J2
これをAWSのisingハミルトニアンと置き換えるようにソースコードを変更し、
QAOAを動かしていきたいと思います。
実装
失敗してみる
失敗してみることが大切なのでまずは失敗してみます。
一応、イジングモデルを0, 1で定義しました。
なのでAWS Braketのサンプルコードのイジングモデル部分の定義、そして、{-1, 1}で表されているマッピングの部分だけソースコードを編集してみます。
イジングモデルの定義
# set Ising matrix
Jfull = nx.adjacency_matrix(G).todense()
Jfull = np.array(Jfull)
# get off-diagonal upper triangular matrix
J = np.triu(Jfull, k=1).astype(np.float64)
こちらのソースコードを先ほど紹介したモデルに変更します。
具体的には以下のように変更します。
J1 = np.array([
[-2,2,2,2,2,0,0,0,2,0,0,0,2,0,0,0],
[0,-2,2,2,0,2,0,0,0,2,0,0,0,2,0,0],
[0,0,-2,2,0,0,2,0,0,0,2,0,0,0,2,0],
[0,0,0,-2,0,0,0,2,0,0,0,2,0,0,0,2],
[0,0,0,0,-2,2,2,2,2,0,0,0,2,0,0,0],
[0,0,0,0,0,-2,2,2,0,2,0,0,0,2,0,0],
[0,0,0,0,0,0,-2,2,0,0,2,0,0,0,2,0],
[0,0,0,0,0,0,0,-2,0,0,0,2,0,0,0,2],
[0,0,0,0,0,0,0,0,-2,2,2,2,2,0,0,0],
[0,0,0,0,0,0,0,0,0,-2,2,2,0,2,0,0],
[0,0,0,0,0,0,0,0,0,0,-2,2,0,0,2,0],
[0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,2],
[0,0,0,0,0,0,0,0,0,0,0,0,-2,2,2,2],
[0,0,0,0,0,0,0,0,0,0,0,0,0,-2,2,2],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,2],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2],
])
J2 = np.array([
[0,0,0,0,0,2,0,2,0,1,0,1,0,3,0,3],
[0,0,0,0,2,0,2,0,1,0,1,0,3,0,3,0],
[0,0,0,0,0,2,0,2,0,1,0,1,0,3,0,3],
[0,0,0,0,2,0,2,0,1,0,1,0,3,0,3,0],
[0,0,0,0,0,0,0,0,0,4,0,4,0,2,0,2],
[0,0,0,0,0,0,0,0,4,0,4,0,2,0,2,0],
[0,0,0,0,0,0,0,0,0,4,0,4,0,2,0,2],
[0,0,0,0,0,0,0,0,4,0,4,0,2,0,2,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2],
[0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2],
[0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
])
J = J1 + 0.2*J2
print(J)
そうするとIsing matrixが以下のような定義になるのがわかるかと思います。
{0,1}への変更
{-1, 1}で表されているマッピングの部分の変更を行います。
util_qaoa.pyの以下の部分を変更します。
変更前
# convert results (0 and 1) to ising (-1 and 1)
meas_ising = result.measurements
meas_ising[meas_ising == 0] = -1
変更後
# convert results (0 and 1) to ising (-1 and 1)
meas_ising = result.measurements
#meas_ising[meas_ising == 0] = -1
実行
では実行してみます。
ValueError: Operator qubit count 2 must be equal to size of target qubit set QubitSet([Qubit(0)])
直訳すると「演算子のビットカウント2はターゲットビットセットの数と同じでなければなりません。」とのことです。
エラーの意味はあまりわかりませんが、
こちらの記事にも書いた通り、AWS Braketのサンプルコードでは{-1, 1}のイジングモデルの線形項や{0, 1}モデルの対角成分を扱うことができません。
なのでそこの部分を修正していこうと思います。
修正①：量子複製不能原理に対しての修正
まず量子コンピュータの基本的な原理として「量子複製不能原理」というものがあります。
古典回路と比べて述べると簡単に説明できるのですが、簡単に言うと
「古典回路ではファンアウト操作により、1つの入力を複製することが可能であったが、量子回路ではそれができない。」といったものになります。
実際にAWS Braketのソースコードを見るとわかるかと思います。
util_qaoa.pyのソースコードの以下のコードを見てみてください。
# helper function for evolution with cost Hamiltonian
def cost_circuit(gamma, n_qubits, ising, device):
"""
returns circuit for evolution with cost Hamiltonian
"""
# instantiate circuit object
circ = Circuit()
# get all non-zero entries (edges) from Ising matrix
idx = ising.nonzero()
edges = list(zip(idx[0], idx[1])) # ここで行列の比ゼロ要素を取得している
# apply ZZ gate for every edge (with corresponding interation strength)
for qubit_pair in edges:
# get interaction strength from Ising matrix
# qubit_pair[0]==qubit_pair[1]があった場合でもZZGateを適用してしまっている。
int_strength = ising[qubit_pair[0], qubit_pair[1]]
# for Rigetti we decompose ZZ using CNOT gates
if device.name == 'Rigetti':
gate = ZZgate(qubit_pair[0], qubit_pair[1], gamma*int_strength)
circ.add(gate)
# classical simulators and IonQ support ZZ gate
else:
gate = Circuit().zz(qubit_pair[0], qubit_pair[1], angle=2*gamma*int_strength)
circ.add(gate)
return circ
上記ソースコードに日本語のコメントで書いている通り、ising変数の対角要素が比ゼロの場合であってもZZGateが適用されるようになっています。
なので上記のfor文の中身の部分を以下に変更すれば量子複製不能原理に触れるエラーは回避できます。
for qubit_pair in edges:
if (qubit_pair[0] != qubit_pair[1]): #この行を追加
# get interaction strength from Ising matrix
int_strength = ising[qubit_pair[0], qubit_pair[1]]
# for Rigetti we decompose ZZ using CNOT gates
if device.name == 'Rigetti':
gate = ZZgate(qubit_pair[0], qubit_pair[1], gamma*int_strength)
circ.add(gate)
# classical simulators and IonQ support ZZ gate
else:
gate = Circuit().zz(qubit_pair[0], qubit_pair[1], angle=2*gamma*int_strength)
circ.add(gate)
実行
変更を行ったのでプログラムを実行してみようと思います。
少し時間はかかりますが、実際に実行すると以下のような結果が返ってきます。
（ランダムなので多少結果は異なります。）
Circuit depth hyperparameter: 3
Problem size: 16
Starting the training.
====================================================================
OPTIMIZATION for circuit depth p=3
Param "verbose" set to False. Will not print intermediate steps.
====================================================================
Optimization terminated successfully.
Current function value: 9.080000
Iterations: 5
Function evaluations: 606
Final average energy (cost): 9.080000000000002
Final angles: [3.84038054 3.6674896 3.95335633 1.96342744 1.12175435 1.33671581]
Training complete.
Code execution time [sec]: 142.866375207901
Optimal energy: -6.6
Optimal classical bitstring: [0 0 0 1 0 0 1 0 1 0 0 0 0 1 0 0]
「Optimal classical bitstring」は$[q_0, q_1, q_2, ..., q_{15}]$の順で並んでいるのでそれをもとに実際の巡回路に落とし込むと「C → D → B → A」となります。（参考：D-Waveとblueqat.optで巡回セールスマン問題とmaxcut問題、クリーク問題、自然数分割問題）
実際の経路の長さを求めると、7となり、明らかに最適解であることがわかります。
(*)式で（添え字を0~15に変更した状態をもとに戻して）計算しなおすと上記の「Optimal energy」になることがわかり、これが最適解であることがわかるかと思います。
実は上記変更は足りない部分があります。
（というより正直最適解が出なくてほしかった。。。のですが、しょうがない。）
なぜかというと今回行った変更を思い出してください。
同じ量子ビット同士の演算を飛ばしています。
なのでそこの部分を表現しないといけません。
と言いつつ、今回最適解を獲得できた理由は厳密にはわかりませんが、以下の二つがあるかと思っています。
（わかる方いたら教えていただけると幸いです。）
- 結局基底状態を求めるためのパラメータ最適化は量子回路の外側で行われています。なのでそれなりの今回の最適解を求めるのに最適なパラメータ設定が行われた。
- 今回の最適化における対角成分は結局のところすべての値が同じ値でした。つまり、(対角成分を表現しない基底状態) $\supset$ (対角成分を表現する基底状態) なので対角成分を表現しない基底状態を求めるパラメータ調整においても一番最適な状態に持っていくことができれば最適化状態となりうる？
といったようなことがあるのかな。と思っています。
（本当に何となくの考察なのでわかる方いたら教えていただけると幸いです。）
ということで対角成分を表現していきましょう。
修正②：対角成分の表現
対角成分を表すのに数式的にどのように表現するかはこちらをぜひ参考にしてください。
結論だけ述べると以下の通り、$R_Z$回転で表すことができます。
（各ビットは0,1で表されるので$q_iq_i = q_i$で一次項として表すことができることに注意）
e^{-i\gamma_{i} Z} = R_{Z}(2\gamma_{i})
実際のソースコードに起こしてみましょう。
ソースコードに起こす際もそれほど難しくなく、ZZGateを適用した後にそれぞれの量子ビットに対して$R_Z$回転を加えます。
# apply ZZ gate for every edge (with corresponding interation strength)
for qubit_pair in edges:
if (qubit_pair[0] != qubit_pair[1]):
# get interaction strength from Ising matrix
int_strength = ising[qubit_pair[0], qubit_pair[1]]
# for Rigetti we decompose ZZ using CNOT gates
if device.name == 'Rigetti':
gate = ZZgate(qubit_pair[0], qubit_pair[1], gamma*int_strength)
circ.add(gate)
# classical simulators and IonQ support ZZ gate
else:
gate = Circuit().zz(qubit_pair[0], qubit_pair[1], angle=2*gamma*int_strength)
circ.add(gate)
# 線形項に対応する量子ゲートを追加
for qubit_pair in edges:
if (qubit_pair[0] == qubit_pair[1]):
gate = Circuit().rz(qubit_pair[0], 2 * gamma)
circ.add(gate)
準備はできたので実際に動かしてみましょう。
実行
変更を行ったのでプログラムを実行してみようと思います。
少し時間はかかりますが、実際に実行すると以下のような結果が返ってきます。
（ランダムなので多少結果は異なります。）
Circuit depth hyperparameter: 3
Problem size: 16
Starting the training.
====================================================================
OPTIMIZATION for circuit depth p=3
Param "verbose" set to False. Will not print intermediate steps.
====================================================================
Optimization terminated successfully.
Current function value: 5.500000
Iterations: 2
Function evaluations: 292
Final average energy (cost): 5.500000000000001
Final angles: [5.66357773 4.31982407 4.72149979 1.63593414 1.0267231 2.03051354]
Training complete.
Code execution time [sec]: 74.82601428031921
Optimal energy: -6.6
Optimal classical bitstring: [0 0 0 1 0 0 1 0 1 0 0 0 0 1 0 0]
「Optimal classical bitstring」は$[q_0, q_1, q_2, ..., q_{15}]$の順で並んでいるのでそれをもとに実際の巡回路に落とし込むと「C → D → B → A」となります。（参考：D-Waveとblueqat.optで巡回セールスマン問題とmaxcut問題、クリーク問題、自然数分割問題）
これは修正①での結果と同じ結果でありますが、巡回順が変わった結果が出力されることや近似アルゴリズムなので最適解出ない結果が出ることももちろんありますが、今回は最適解であることがわかります。
以上のようにして対角成分を表した計算をAWS Braket上で行うことができました。
最後に
今回はAWS Braket上でより一般的に近い最適化問題を解いてみました。
今後もより柔軟な回路設計をして最適化問題を解いてみようかと思っているのでまた記事にできそうならばしてみようかと思います。