Rigettiの量子ゲート量子コンピュータで組合せ最適化問題をやってみよう！

amazon braketを使って本日は簡単な組合せ最適化問題を解いてみましょう。

量子ゲート量子コンピュータで組合せ最適化問題を解くにはQAOAもしくはVQEを使うのが一般的となっています。QAOAは量子アニーリングのいとこみたいなものです。VQEは行列の固有値を探すための力業アルゴリズムです。

QAOAでいい答え出すの大変なので、シミュレータ後に実機に投げてみます。QAOAについて知りたい方はblueqat.comの中でいろいろ検索してもらうといろいろ出てきます。今回はある程度知ってる上でちゃっちゃか作ってしまいます。

from braket.circuits import Circuit
from braket.devices import LocalSimulator
import matplotlib.pyplot as plt
%matplotlib inline

まず組合せ最適化問題を解くためには、行列の固有値を求めればいいので、行列を設定します。

H = Z_0*Z_1 + Z_0

この式のZは値の範囲が-1から+1をとります。各項が-1-1=-2の時に最小値になると思われます。その時にZ0は-1、Z1は1になればOKです。
その際の量子ビットの値は、0番目の量子ビットが1、1番目の量子ビットが0です。

その際に式はほにゃほにゃして、パラメータを２つ用意します（QAOAにしたいばあいは）。
断熱時間発展化しますが、上記の係数は両方とも1なので、回路はRZZとRZになります。

|0> -H-RZZ-RZ--RX--M
|0> -H-RZZ-----RX--M

最初に固有状態を重ね合わせで準備し、RXはXの時間発展を用意し横磁場とします。
RZZとRZのパラメータは共通です。

今回はrigettiをつかってみます。ゲートセットを確認します。

from braket.aws import AwsDevice
device = AwsDevice("arn:aws:braket:us-west-1::device/qpu/rigetti/Aspen-M-2")
print(device.properties.action['braket.ir.jaqcd.program'].supportedOperations)

['cz', 'xy', 'ccnot', 'cnot', 'cphaseshift', 'cphaseshift00', 'cphaseshift01', 'cphaseshift10', 'cswap', 'h', 'i', 'iswap', 'phaseshift', 'pswap', 'rx', 'ry', 'rz', 's', 'si', 'swap', 't', 'ti', 'x', 'y', 'z', 'start_verbatim_box', 'end_verbatim_box']

さすがにRZZはないので、かわりにcnot-RZ-cnotで代用します。角度ですが、初期値はパラメータを二つ用意します。

import numpy as np

params = [np.random.rand(), np.random.rand()]

circuit = Circuit().h(0).h(1).cnot(0,1).rz(1,params[0]).cnot(0,1).rz(0,params[0]).rx(0,params[1]).rx(1,params[1])
print(circuit)

T  : |0|1|   2    |3|   4    |   5    |
                                       
q0 : -H-C----------C-Rz(0.81)-Rx(0.08)-
        |          |                   
q1 : -H-X-Rz(0.81)-X-Rx(0.08)----------

T  : |0|1|   2    |3|   4    |   5    |

これで出来上がりました。QAOAでは最小化したいのは期待値ですが、Z0Z1の期待値は同符号と異符号の回数を数えて簡単に求まります。またZ0の期待値は0と1の出る回数を数えればOKです。

device = LocalSimulator()
shots = 10000
result = device.run(circuit, shots=shots).result()

#計算結果の表示
counts = result.measurement_counts
print(counts)

#グラフを描画
# plot using Counter
plt.bar(counts.keys(), counts.values());
plt.xlabel('bitstrings');
plt.ylabel('counts')

Counter({'00': 2805, '11': 2566, '01': 2408, '10': 2221})

Text(0, 0.5, 'counts')

<Figure size 432x288 with 1 Axes>

expectation = ((counts['00']+counts['11']-counts['01']-counts['10']) + (counts['00']+counts['01']-counts['10']-counts['11']))/shots
print(expectation)

0.1168

これを使って最適化しましょう。最適化はScipyのCOBYLA法を使ってみます。

import scipy.optimize

shots = 100000

def f(para):
    circuit = Circuit().h(0).h(1).cnot(0,1).rz(1,-2*para[0]).cnot(0,1).rz(0,-2*para[0]).rx(0,-2*para[1]).rx(1,-2*para[1])
    result = device.run(circuit, shots=shots).result()
    counts = result.measurement_counts
    expectation = ((counts['00']+counts['11']-counts['01']-counts['10']) + (counts['00']+counts['01']-counts['10']-counts['11']))/shots
    return expectation

res = scipy.optimize.minimize(f, params, method='COBYLA')
print(res)

     fun: -0.55094
   maxcv: 0.0
 message: 'Optimization terminated successfully.'
    nfev: 31
  status: 1
 success: True
       x: array([ 2.19370847, -1.09945378])

circuit = Circuit().h(0).h(1).cnot(0,1).rz(1,-2*res.x[0]).cnot(0,1).rz(0,-2*res.x[0]).rx(0,-2*res.x[1]).rx(1,-2*res.x[1])
result = device.run(circuit, shots=shots).result()

#計算結果の表示
counts = result.measurement_counts
print(counts)

#グラフを描画
# plot using Counter
plt.bar(counts.keys(), counts.values());
plt.xlabel('bitstrings');
plt.ylabel('counts')

Counter({'10': 38645, '01': 26491, '11': 23699, '00': 11165})

Text(0, 0.5, 'counts')

<Figure size 432x288 with 1 Axes>

期待値は下がり切れていませんが、答えである10の出現確率が高くなっているのが確認できました。
再度Rigettiを呼び出してみます。

device = AwsDevice("arn:aws:braket:us-west-1::device/qpu/rigetti/Aspen-M-2")
print(device.properties.action['braket.ir.jaqcd.program'].supportedOperations)

['ccnot', 'cnot', 'cphaseshift', 'cswap', 'cy', 'cz', 'h', 'i', 'phaseshift', 'rx', 'ry', 'rz', 's', 'si', 'swap', 't', 'ti', 'v', 'vi', 'x', 'y', 'z', 'ecr', 'start_verbatim_box', 'end_verbatim_box']

circuit = Circuit().h(0).h(1).cnot(0,1).rz(1,-2*res.x[0]).cnot(0,1).rz(0,-2*res.x[0]).rx(0,-2*res.x[1]).rx(1,-2*res.x[1])
task = device.run(circuit) 
print(task.result().measurement_counts)

#そのまま描画
# plot using Counter
counts = task.result().measurement_counts
plt.bar(counts.keys(), counts.values());
plt.xlabel('bitstrings');
plt.ylabel('counts')

Counter({'10': 404, '11': 244, '01': 224, '00': 128})

Text(0, 0.5, 'counts')

<Figure size 432x288 with 1 Axes>

正解が出ました。実機でもQAOAができましたね！以上です。

Rigettiの量子ゲート量子コンピュータで組合せ最適化問題をやってみよう！

Yuichiro Minato