こんにちは。ガチャって面白いですよね。量子コンピュータと量子アニーラを使って、ガチャを作りました。それぞれ作り方が違いますので、解説します。
イオントラップ型量子コンピュータ
原子核をイオン化し、レーザーを打ち込んで計算をするタイプです。今回はamazonのサービス経由でblueqatを組み合わせながらIonQというマシンで計算をしました。
量子ゲートボルンガチャ
ボルンの規則と呼ばれる測定をしたときの確率分布のルールがあります。この規則を利用して機械学習をするとき「ボルンマシン」と呼ばれることがあります。今回は状態ベクトルを使った「ボルンガチャ」を作りました。
利用機材
awsのbraketSDK、sagemaker、blueqatSDKを利用しました。余裕があればここにpytorchが入りますが、今回は小さいサイズでしたので、これで。sagemaker上でbraketSDKとblueqatSDKを動かしました。sagemakerを使わず、ローカルに入れても良いと思います。
インストール
sagemaker上ではbraketSDKはインストールしてるので、blueqatSDKを入れました。
# Check SDK version
!pip show amazon-braket-sdk | grep Version
Version: 1.0.0.post1
!pip install blueqat
無事これで入りました。
ツール読み込みblueqat
今回の学習はイジング定式化みたいに適当にやるにはパラメータが多いので、量子回路を作りますが、NISQ変分量子古典ハイブリッド回路を利用します。
from blueqat import Circuit
import matplotlib.pyplot as plt
import numpy as np
import time
%matplotlib inline
再現性を確保するために、乱数のシード固定し、初期の角度パラメータを5つ用意しました。
np.random.seed(39)
ainit = [np.random.rand()*np.pi*2 for i in range(5)]
print(ainit)
[3.436205910548335, 5.0133474071358375, 5.154737037760021, 0.7668619300107449, 3.7824901639645887]
In [68]:
量子回路を準備
今回の量子回路を準備します。全部で角度パラメータは5つ。任意回転4つとXXゲート1つです。初期の角度は適当に入れてみます。XXゲートがないマシンでは、CXで良いと思います。
def circ(a):
return Circuit().rx(a[0])[0].rz(a[1])[0].rx(a[2])[1].rz(a[3])[1].rxx(a[4])[0,1].run()
ガチャの出現確率を準備
今回のガチャのターゲットの確率分布は、
N 60%
R 30%
HR 8%
SSR 2%
としました。これをターゲットにしてみたいと思います。4要素の確率分布を作るには、2量子ビットあれば十分です。Nを00、Rを01、HRを10、SSRを11の量子状態に対応させ、出現確率の平方根を求めて、確率振幅に直します。
#target state vector
tgt = [np.sqrt(0.6),np.sqrt(0.3),np.sqrt(0.08),np.sqrt(0.02)]
print(tgt)
これによって、
[0.7745966692414834, 0.5477225575051661, 0.282842712474619, 0.1414213562373095]
In [72]:
こちらの状態ベクトルを作ることができれば、ガチャを回すことができます。今回は量子アニーラのようにイジングエネルギーでボルツマン分布を作るのに対して、直接状態ベクトルを触るので、「ボルンガチャ」と命名します。
早速学習
解説は省略しますが、初期パラメータからターゲットとの損失関数を求め、角度パラメータを微小変化させて微分をとり、勾配法で最適化をかけます。今回はちょっとズルをして、状態ベクトルを直接学習させてます。本来は、サンプリング手法を使って期待値計算をして学習をします。
#initial parameters
a = ainit.copy()
#result list
ar = []
h = 0.01
e = 0.01
#iterations
nsteps = 1000
start = time.time()
for i in range(nsteps):
r = np.random.randint(0,4)
c = circ(a)
loss = (np.abs(c[0]) - tgt[0])**2 + (np.abs(c[1]) - tgt[1])**2 + (np.abs(c[2]) - tgt[2])**2 + (np.abs(c[3]) - tgt[3])**2
ar.append(loss)
at = [0 for i in range(len(a))]
for j in range(len(a)):
aa = a.copy()
aa[j] += h
cc = circ(aa)
loss2 = (np.abs(cc[0]) - tgt[0])**2 + (np.abs(cc[1]) - tgt[1])**2 + (np.abs(cc[2]) - tgt[2])**2 + (np.abs(cc[3]) - tgt[3])**2
at[j] = a[j] - e*(loss2 - loss)/h
a = at
plt.plot(ar)
plt.show()
print(time.time() - start)
print(loss)
2.693882942199707
0.0008066113653701334
もっと頑張れますが、とりあえず良さそうです。上手い感じに収束し、学習してくれました。最適化された角度は、
print(a)
こうなりました。
[2.2016688877066573, 4.827924322644809, 4.311839683930082, 0.4078438118319338, 4.022940030803429]
In [105]:
確認してみる。
では、確認してみます。今学習した角度パラメータを導入した量子回路で、今度は状態ベクトルではなく、サンプリングをしてみます。
Circuit().rx(a[0])[0].rz(a[1])[0].rx(a[2])[1].rz(a[3])[1].rxx(a[4])[0,1].m[:].run(shots=1000)
これによって、
Counter({'00': 596, '10': 296, '11': 29, '01': 79})
In [85]:
どうでしょう。良い感じですね。実際にサンプリングをして、ちょうど良い感じのサンプルが取れました。これがガチャに対応します。
いよいよ量子コンピュータ
ここまではシミュレータで学習をしました。学習した内容を使って実機計算します。今回利用するのは、アメリカのメリーランド州にあるIonQというマシンです。
# Import Braket libraries
from braket.circuits import circuit, Circuit, Gate, Moments
from braket.circuits.instruction import Instruction
from braket.aws import AwsQuantumTask, AwsDevice
from braket.devices import LocalSimulator
import matplotlib.pyplot as plt
# magic word for producing visualizations in notebook
%matplotlib inline
import numpy as np
たくさんツールを読み込みます。そして、結果の格納のフォルダ指定します。
# Please enter the S3 bucket you created during onboarding in the code below
my_bucket = f"自分のバケット" # the name of the bucket
my_prefix = "simulation-output" # the name of the folder in the bucket
s3_folder = (my_bucket, my_prefix)
良い感じです。
量子回路の実装
blueqatで作った回路をbraketに移植します。
circ = Circuit();
circ.rx(0,a[0]).rz(0,a[1]).rx(1,a[2]).rz(1,a[3]).xx(0,1,a[4])
print(circ)
任意回転ゲートの定義がblueqatと同じなのかちょっと気になりましたが、時間がないのでそのままやりました。記述方法はちょっと違いますね。
T : | 0 | 1 | 2 |
q0 : -Rx(2.2)--Rz(4.83)--XX(4.02)-
|
q1 : -Rx(4.31)-Rz(0.408)-XX(4.02)-
T : | 0 | 1 | 2 |
できました。
実機準備
まずは、デバイスをionqに設定、shotsを今回は1000円ガチャですので、1shot=1円くらいなので、1000shots設定しました。回路、結果フォルダ、ショット、ポーリングのタイムアウトを設定して、タスクを投げます。
ionq = AwsDevice("arn:aws:braket:::device/qpu/ionq/ionQdevice")
ionq_task = ionq.run(circ, s3_folder, shots=1000, poll_timeout_seconds=5*24*60*60)
ionq_task_id = ionq_task.id
ionq_status = ionq_task.state()
print('Status of task:', ionq_status)
これで無事登録されました。うまくいくと、
Status of task: CREATED
がでます。
待つ
イオントラップはキューに格納されてから処理されるのが遅いです。待ちます。。。
# print status
status = ionq_task.state()
print('Status of (reconstructed) task:', status)
Status of (reconstructed) task: QUEUED
In [107]:
寝ます。
終了
寝て起きたら終わってました。
Created
Nov 11, 2020 15:02 (UTC)
Ended
Nov 11, 2020 15:43 (UTC)
今回はかなり早かったです。ジョブの投入から41分で処理されました。
結果の取得と表示
今回は結果の取得と表示は手元のpythonで確認します。
# recover task
task_load = AwsQuantumTask(arn="arn:aws:braket:us-east-1:hogehoge")
# print status
status = task_load.state()
print('Status of (reconstructed) task:', status)
if status == 'COMPLETED':
# get results
results = task_load.result()
# get all metadata of submitted task
metadata = task_load.metadata()
# example for metadata
shots = metadata['shots']
machine = metadata['deviceArn']
# print example metadata
print("{} shots taken on machine {}.".format(shots, machine))
# get measurement counts
counts = results.measurement_counts
print('Measurement counts:', counts)
# plot results: see effects of noise
plt.bar(counts.keys(), counts.values());
plt.xlabel('bitstrings');
plt.ylabel('counts');
plt.tight_layout();
plt.savefig('ionq.png', dpi=700);
elif status in ['FAILED', 'CANCELLED']:
# print terminal message
print('Your task is in terminal status, but has not completed.')
else:
# print current status
print('Sorry, your task is still being processed and has not been finalized yet.')
そうすると、
Status of (reconstructed) task: COMPLETED
1000 shots taken on machine arn:aws:braket:::device/qpu/ionq/ionQdevice.
Measurement counts: Counter({'00': 508, '10': 284, '01': 134, '11': 74})
良い感じでカウントが取れました。
Nは50.8%、Rは28.4%、HRは13.4%、SSRは7.4%となりました。結果としてはボチボチではないでしょうか。ボルンガチャが完成しました！是非やってみてください！1000連で1000円かかります！
次に量子アニーラを利用したボルツマンガチャ
上記はボルンの規則を使いました。量子アニーリングでは、ボルツマン分布を利用することで同じようなことができます。上記は量子回路を作りましたが、量子アニーリングでは回路を作らず式を作ります。
ボルツマン分布とQUBO
早速ガチャやります。今回はボルツマン分布を使ってみます。イジングモデルのエネルギー値に対して、確率分布である、
p(x) = 1/Z * exp(beta*E)
が成立します。これを利用し、QUBOを作り、エネルギー分布とガチャのレア度を対応させます。
ツール読み込み
今回は純粋にD-WaveのAdvantageを呼び出します。
import numpy as np
from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite
早速サンプリング
早速QUBOを作って、実機に投げてみます。出る確率分布は上記のエネルギーの値により、QUBOの形によります。今回はIonQと似た分布の式をあらかじめ作っておきました。tokenは自分のをあらかじめ設定します。
qubo = np.array([[0,-2,10],[0,0,-1],[0,0,-0.5]])
bqm = dimod.BQM.from_numpy_matrix(qubo)
response = EmbeddingComposite(DWaveSampler(endpoint='https://cloud.dwavesys.com/sapi', token=token, solver='Advantage_system1.1')).sample(bqm, num_reads=1000)
print(response)
作るQUBOによって分布が変わります。今回はAdvantage system 1.1に1000回投げました。
確率分布
結果はこのような感じででます。
0 1 2 energy num_oc. chain_.
0 1 1 0 -2.0 551 0.0
1 0 1 1 -1.5 375 0.0
2 0 0 1 -0.5 38 0.0
3 1 0 0 0.0 11 0.0
4 0 0 0 0.0 15 0.0
5 0 1 0 0.0 10 0.0
['BINARY', 6 rows, 1000 samples, 3 variables]
110が551回。011が375回、その他続くという感じです。背景にある確率分布は厳密に求めることができますが、今回はサンプリングですので、それに準じて出現回数が決まっています。
量子アニーラの特徴は応答が速いことです。量子ゲートマシンで結果を取得するのは時間がかかるので、今回はアナログのD-Waveの量子アニーラのサンプリングを使います。
理論的な話
実際に、背後の確率分布はそれぞれ、出ているエネルギー値-2.0や-1.5などによって変わります。サイズが小さい場合には、厳密な確率を求めることができますが、実際にはアナログマシンですので、マシンの性能や状態によっても分布が変わりますので、事前に実機で分布を取るのが望ましいです。
上記は、実機では、SSRが出る確率は、およそ1%となっています。この辺り調整するには、量子アニーラではイジング式をいじれば良いことになります。
ゲームは面白い
量子コンピュータや量子アニーラの不確実性をゲームに導入するのはとても面白いと思いますので、ぜひ新しい分野に挑戦してみてください。