Quantum Kitchen Sinks

[1806.08321] Quantum Kitchen Sinks: An algorithm for machine learning on near-term quantum computers
が気になったので、そのメモ。

どんなもの?
Quantum Kitchen Sinks(QKS) は、一連の量子回路からランダムサンプリングし、各々の回路を用いて入力データを測定されたビット列へ非線形変換を実現するものである。このアルゴリズムによって量子回路のパラメータの最適化にかかるコストを排除することができる!!そのあと、連結された結果は、(量子ではない通常の)機械学習アルゴリズムを用いて処理する(量子回路の効果が分かるように、単純なものを用いる)。
f:id:puyokw:20180626230407p:plain:w400

技術や手法のキモはどこ?
Quantum Kitchen Sinks(QKS) は、Random Kitchen Sinks(RKS) にinspire されたもので、RKS では例えばcos 関数が使用されているが、これは単一の量子ビットにRabi 回転(下図では, RX) を適用することでそのような変換ができる。そのときのパラメータはランダムに決める。

f:id:puyokw:20180626233156p:plain:w400
上図の(a) とか(c) はうまくいくが、(b) はダメらしい。

入力データは実数だけど、どうやって量子回路に入れるの??
前処理:
 \Theta_{i,e} = \Omega_e u_i + \beta_e
ただし,  \Omega_e \mathcal N(0, \sigma^2) に,  \beta_e \mathcal U(0, 2\pi) に従う。
この \theta RX(\theta) に入れることになります。

メインの量子回路の処理を行う。

後処理:観測を行って, 出力される0, 1 を全エピソード数E 個分つなげていく(ここの部分は, 複数回観測した合計や平均などでも可) 。パラメータ付き量子回路のパラメータ数が4 個なら, 各トレーニングデータは長さ4E のベクトルになる。
E は大きい方が良い。

どうやって有効だと検証した?
(1) "picture frames" の二値分類
|- (a) logistic regression (scikit-learn) (~50%)
|- (b) QKS + logistic regression (> 99.9%)

(2) MNIST の数字の3 と5 の二値分類
|- (a) logistic regression (95.9%)
|- (b) QKS + logistic regression (98.6%)

ロジスティク回帰のみと量子回路+ロジスティク回帰を比較していて、QKS の段階で量子回路をたくさん作ることで精度を向上させることができる。
f:id:puyokw:20180626231513p:plain

もう一つの戦略として, qubit 数を増やすのもあるが, こちらでは途中から過学習している。
f:id:puyokw:20180626231538p:plain:w400
これは4 qubit の場合の例で, qubit 数を増やしたら精度が上がりそうに思うのですが…
f:id:puyokw:20180627012940p:plain:w200

議論はある?
・CNOT を含むある特定のネットワークが準最適である可能性がある。
・MNIST のデータセットでは、大きな回路のパラメータをチューニングするには不十分な大きさである。

今回試されていた回路は例えば以下のようなRX とCNOT を組み合わせたものでした。
例えば、4 qubit の場合だと以下のような構成。
RX(%x0) 0
RX(%x1) 1
RX(%x2) 2
RX(%x3) 3
CNOT 0 2
CNOT 1 3
CNOT 0 1
CNOT 2 3

感想:
量子回路の設計方針が分かりやすい!


TO DO
・再現実装をしてみる(まだ雰囲気程度)

from pyquil.quil import Program
from pyquil.api import QVMConnection
from pyquil.gates import *
from sklearn import datasets
from sklearn.linear_model import LogisticRegression
import numpy as np
import math
from tqdm import trange
from sklearn.metrics import log_loss, accuracy_score
from sklearn.cross_validation import StratifiedKFold
from sklearn.preprocessing import StandardScaler, MinMaxScaler
import pandas as pd

def rand():
	return np.random.rand()

def rand_u():
	return np.random.ranf()*2-1

def rand_n():
	return np.random.normal(0,1)

def make_pictureframe(num_data):
	data=np.zeros((num_data,2))
	teacher=np.zeros(num_data).astype(np.int32)
	for i in range(num_data):
		if rand()<=0.5:
			teacher[i]=0
		else:
			teacher[i]=1
	for i in range(num_data):
		tmp = rand()
		if tmp <= 0.25:
			data[i][0] = 1.0*rand_u() + rand_u()*0.1
			if rand() <= 0.5:
				data[i][1] = -1.0 + rand_u()*0.1
			else:
				data[i][1] = 1.0 + rand_u()*0.1
		elif tmp <= 0.5:
			data[i][0] = -1.0*rand_u() + rand_u()*0.1
			if rand() <= 0.5:
				data[i][1] = -1.0 + rand_u()*0.1
			else:
				data[i][1] = 1.0 + rand_u()*0.1
		elif tmp <= 0.75:
			data[i][1] = -1.0*rand_u() + rand_u()*0.1
			if rand() <= 0.5:
				data[i][0] = -1.0 + rand_u()*0.1
			else:
				data[i][0] = 1.0 + rand_u()*0.1
		else:
			data[i][1] = 1.0*rand_u() + rand_u()*0.1
			if rand() <= 0.5:
				data[i][0] = -1.0 + rand_u()*0.1
			else:
				data[i][0] = 1.0 + rand_u()*0.1
	for i in range(num_data):
		if teacher[i] == 1:
			data[i][0] *= 2
			data[i][1] *= 2
	return [data, teacher]


#np.random.seed(1)
data, target = make_pictureframe(100)

num_param=2
num_data=data.shape[0]
p=data.shape[1]
q=num_param # number of qubit
n_qubits=num_param # number of qubit
r=int(p/q)

num_episode=200
sigma=1

#np.random.normal(0, sigma)

Omega=[]
beta=[]
for e in range(num_episode):
	Omega.append( (np.concatenate( (np.random.randn(r, q), np.zeros((p-r,q))), axis=0)).T )
	beta.append( np.random.rand(q)*2*math.pi)

Omega=np.array(Omega)
beta=np.array(beta)

data_transform=[]
for i in trange(num_data):
	tmp=[]
	for e in trange(num_episode):
		u=data[i]
		theta = np.dot(Omega[e], u) + beta[e]
		qvm = QVMConnection()
		p1=Program()
		for j in range(q):
			p1.inst(RX(theta[j],j))
		
		p1.inst(CNOT(0, 1))
		p1.inst(MEASURE(0,0))
		p1.inst(MEASURE(1,1))
		tmp.append(np.array(qvm.run(p1, [j for j in range(q)], 10000)).mean(axis=0))
		
	tmp=[item for sublist in tmp for item in sublist]
	data_transform.append(tmp)

data_transform = pd.DataFrame(data_transform)
data_transform.columns=["col_"+str(i) for i in range(data_transform.shape[1])]
#data_transform.to_csv("trained.csv", index=False)

#train = pd.read_csv("trained.csv")
train=data_transform
train = np.array(train)

skf=StratifiedKFold(target, n_folds=5, shuffle=True)
for tr, te in skf:
	X_train, X_valid, y_train, y_valid = train[tr], train[te], target[tr], target[te]

lr = LogisticRegression()
lr.fit(X_train, y_train)

pred=lr.predict_proba(X_valid)
test_logloss=log_loss(y_valid, pred)
test_acc=lr.score(X_valid, y_valid)
pred=lr.predict_proba(X_train)
train_logloss=log_loss(y_train, pred)
train_acc=lr.score(X_train, y_train)
print("train acc: {}, test acc: {}".format(train_acc, test_acc))
print("train logloss: {}, test logloss {}".format(train_logloss, test_logloss))

pyquil は遅いので, 他のシミュレータで試した感じだと85% 程度しか精度でない....
ノイズがあった方が良いのかな??

Rigetti のpyquil のお試し

IBMQ 以外にも量子プログラミングができるプラットフォームとして Rigetti Forest があります。
Rigetti
Rigetti Forest のシミュレータは, メールアドレスを登録することで使用できます。実機での実行するためには, アカウントのアップグレードが必要なため, ここではシミュレーションの方を試します。python のpyquil を使用します。
GitHub - rigetticomputing/pyquil: A Python library for quantum programming using Quil.

pip install pyquil

でインストールできます。

メールアドレスを登録すると, key とuser_id が貰えますので, それを以下のような内容の .pyquil_config ファイルを作成します。

[Rigetti Forest]
key: <Rigetti Forest API key>
user_id: <Rigetti User ID>

Installation and Getting Started — pyQuil 2.0.0.dev0 documentation
の注意書きにあるように, Linux では /users/username に, Mac では /Users/Username に, Windows では C:\Users\Username にこのファイルを配置します。

以下では, 回路の作成, 回路の図示を行います。

from pyquil.quil import Program
from pyquil.api import QVMConnection
from pyquil.gates import *
qvm = QVMConnection()
p1=Program()
p1.inst(H(0),
        CNOT(0,1),
        H(1),
        MEASURE(0,0),
        MEASURE(1,1)) # 回路を作成
print(p1)
# 出力
#H 0
#CNOT 0 1
#H 1
#MEASURE 0 [0]
#MEASURE 1 [1]

pyquil では柔軟に回路を作成することができます。

p1.inst(H(0),
        CNOT(0,1),
        H(1),
        MEASURE(0,0),
        MEASURE(1,1)) 
の部分は
p1.inst(H(0))
p1.inst(CNOT(0,1))
p1.inst(H(1))
p1.inst(MEASURE(0,0))
p1.inst(MEASURE(1,1)) 
と書いても同じです。

また複数の回路を結合することもできます。

q1=Program()
q1.inst(Z(0),I(0))
q1.inst(X(1))
print(q1)
# 出力
# Z 0
# I 0
# X 1
# q1 の後にp1 を結合する
q1.inst(p1)
print(q1)
# 出力結果
# Z 0
# I 0
# X 1
# H 0
# CNOT 0 1
# H 1
# MEASURE 0 [0]
# MEASURE 1 [1]

シミュレーションの実行は次のようにして行います。

result=qvm.run(p1, [0,1], 10)
print(result) 
# [[1, 0], [0, 1], [1, 0], [0, 0], [0, 0], [0, 1], [1, 0], [0, 0], [1, 0], [0, 1]]

試行回数を増やした版で再度実験して集計します。

result=qvm.run(p1, [0,1], 8000)
# 各サブリストの中をまとめる
L=[]
for sublist in result:
    L.append(''.join([str(l) for l in sublist]))

# pandas を用いて集計する
import pandas as pd 
import matplotlib.pyplot as plt
df=pd.DataFrame({"state": L})
df['state'].value_counts().plot.bar()
plt.show()

f:id:puyokw:20180408040019p:plain

from pyquil.latex.latex_generation import to_latex 
tmpdir='tmp/'
filename='circuit1'
dc=to_latex(q1.inst()) # q1 の回路をtex 形式に変換
with open(tmpdir+filename+'.tex','w') as f:
    f.write(dc) 

回路をtex 形式に変換するまではpyquil でできますが, そこからpdflatex を用いてpdf に変換します。それをさらにjpeg などに変換して読み込みます。

import os
os.system("pdflatex -output-directory {} {}".format(tmpdir, filename+".tex"))
os.system("pdftocairo -jpeg {}.pdf {}".format(tmpdir+filename, tmpdir+filename))

from PIL import Image 
img = Image.open(os.path.join('./'+tmpdir+filename+'-1.jpg'))
img.show()

f:id:puyokw:20180408035143j:plain
他にもsympy のqasm を利用することでも回路図を表示することもできます。

from sympy.physics.quantum.qasm import Qasm
from sympy.physics.quantum.gate import IdentityGate
import matplotlib.pyplot as plt

# Qasm では, identity gate が含まれていないため, 継承して追加します
class QasmEx(Qasm):
    def identity(self, arg):
        self.circuit.append(IdentityGate(self.index(arg)))

# 回路を図示
def VisualizeCircuit(circuit):
    STR=str(circuit).split('\n')[:-1]
    SSTR=[]
    for s in STR:
        s_tmp=s.split(' ')
        if s_tmp[0]=='MEASURE':
            s_tmp=s_tmp[:-1]
        s_tmp2=s_tmp[0]
        if len(s_tmp)>1:
            s_tmp2+=' '
            s_tmp2+=','.join(list(map(lambda x: 'q'+x,s_tmp[1::])))
        SSTR.append(s_tmp2)
    SSTR=list(map(lambda x: x.replace("MEASURE","measure").replace("I","identity"),SSTR))
    SSTR=["qubit "+"q"+str(i) for i in range(2)]+SSTR
    qprog1=QasmEx()
    for i in range(len(SSTR)):
        qprog1.add(SSTR[i])
    qprog1.plot()

VisualizeCircuit(q1.inst())
plt.show()

f:id:puyokw:20180409023857p:plain

pyquil の印象としては, 回路を書くのは便利そうな感じがしました。

IBMQ のお試し

量子コンピュータの基本的な使い方のメモです。まず, IBMQ は以下のサイトから利用できます。
IBM Q Experience
簡単な回路であれば, GUI で試してみるのが良いと思います。
ドラッグ&ドロップで簡単に使うことができます。ただし, 使用するバックエンド(ibmqx2, ibmqx4 など)によって回路が異なるため, 制御NOT (CNOT)を使用するときにどのノードにおけるかは制限があります。
f:id:puyokw:20180407170318j:plain
simulate を押せばシミュレーションが実行でき, run を押せば実機で実行できます。

本格的に使いたい場合には, qiskit を使用するのが良さそうです。
GitHub - QISKit/qiskit-sdk-py: Quantum Software Development Kit for writing quantum computing experiments, programs, and applications.
python からibmqを使用することができます。インストールは

pip install qiskit

で, できます。
実機で実行するには, api token が必要となります。
github のinstallation の項に従って, IBM Q Experience のアカウントを作成して, api token を発行します。
(ユーザー名の部分をクリックして, My Account を選択。Advanced の項目の中にAPI Token が存在するので, そこのRegenerate をクリックして作成。)

Qconfig.py として以下のようなコードを準備します。(API token は先程取得したものを記載します。)

APItoken = 'PUT_YOUR_API_TOKEN_HERE'

config = {
    'url': 'https://quantumexperience.ng.bluemix.net/api',

    # If you have access to IBM Q features, you also need to fill the "hub",
    # "group", and "project" details. Replace "None" on the lines below
    # with your details from Quantum Experience, quoting the strings, for
    # example: 'hub': 'my_hub'
    # You will also need to update the 'url' above, pointing it to your custom
    # URL for IBM Q.
    'hub': None,
    'group': None,
    'project': None
}

if 'APItoken' not in locals():
    raise Exception('Please set up your access token. See Qconfig.py.')

これで準備が完了です。
以下では, 具体的に回路を作成して, 回路図の表示, シミュレーションの実行を行います。
f:id:puyokw:20180409035411j:plain
このようなtoy example を考えてみます。
回路自体は次のように書けます。

from qiskit import QuantumProgram
qp = QuantumProgram()

nq = 3  # number of qubits
q = qp.create_quantum_register("q", nq) # nq 量子ビット準備する
c = qp.create_classical_register("c", nq) # 測定用の領域を確保

circuits = ['testQ'] # 回路の名前を入力
testQ = qp.create_circuit(circuits[0], [q], [c]) # 回路の作成準備
# 1 列目
testQ.h(q[0]) # q0 にH を作用させる
testQ.h(q[1]) # q1 にH を作用させる
testQ.h(q[2]) # q2 にH を作用させる

# 2 列目
testQ.z(q[1]) # q1 にz を作用させる
testQ.cx(q[1],q[2]) # q1, q2 にcnotを作用させる

# 3 列目
testQ.h(q[2]) # q2 にH を作用させる

# 測定
testQ.measure(q[0], c[0]) # q0 を測定し, c0 に結果を格納
testQ.measure(q[1], c[1]) # q1 を測定し, c1 に結果を格納
testQ.measure(q[2], c[2]) # q2 を測定し, c2 に結果を格納

これをシミュレーションして, 結果をヒストグラムで表示します。

# 実行する
results = qp.execute(circuits, backend='local_qasm_simulator', shots=8192, seed=1) 
# 結果を図示する
from qiskit.tools.visualization import plot_histogram
plot_histogram(results.get_counts(circuits[0]))

f:id:puyokw:20180409035645p:plain
結果の状態(000や001など)では, 右側からq0, q1, q2 を示しています。この場合q2 は常に0 の状態です。

回路図を表示するのがやや面倒で, まずlatex_drawer を用いることで, tex ファイルに変換します。ここからは各自の環境に依存しますが, pdflatex を用いてpdf に変換して, さらにpdf2image もしくはpoppler のpdftocairo を用いて画像に変換して読み込みます。

import os 
from qiskit.tools.visualization import latex_drawer
from PIL import Image 
# ipython を使用する場合
# from IPython.display import Image 

filename='circuit'
tmpdir='tmp/'
if not os.path.exists(tmpdir):
    os.makedirs(tmpdir)

latex_drawer(testQ, tmpdir+filename+".tex", basis="h,z") # tex ファイルの作成
os.system("pdflatex -output-directory {} {}".format(tmpdir, filename+".tex")) # tex -> pdf 
os.system("pdftocairo -jpeg {}.pdf {}".format(tmpdir+filename, tmpdir+filename)) # pdf -> jpeg
img = Image.open(os.path.join('./'+tmpdir+filename+'-1.jpg'))
img.show()
#Image(filename=os.path.join(tmpdir+filename+'-1.jpg')) # ipython を使用する場合
#pdf2image.convert_from_path(tmpdir+filename+".pdf") # pdf2image が使える場合は, こちらでもOK 

上記のようなコードにより, 回路図が表示されます。
sympy のqasm を用いても回路図を表示できます。まず, 回路のコードを表示します。

QASM_source = qp.get_qasm('testQ')
print(QASM_source)
#OPENQASM 2.0;
#include "qelib1.inc";
#qreg q[3];
#creg c[3];
#h q[0];
#h q[1];
#h q[2];
#z q[1];
#cx q[1],q[2];
#h q[2];
#measure q[0] -> c[0];
#measure q[1] -> c[1];
#measure q[2] -> c[2];

次に, plot するために文字列を変換していきます。

from sympy.physics.quantum.qasm import Qasm

STR=str(QASM_source).replace(';','').replace('[','').replace(']','').replace('cx','cnot').split('\n')[4:-1]
SSTR=[]
for s in STR:
    s_tmp=s.split(' ')
    if s_tmp[0]=='measure':
        s_tmp=s_tmp[:-2]
    s_tmp2=s_tmp[0]
    if len(s_tmp)>1:
        s_tmp2+=' '
        s_tmp2+=','.join(s_tmp[1::])
    SSTR.append(s_tmp2)

SSTR=["qubit "+"q"+str(i) for i in range(nq)]+SSTR
print(SSTR)
# ['qubit q0', 'qubit q1', 'qubit q2', 'h q0', 'h q1', 'h q2', 'z q1', 'cnot q1,q2', 'h q2', 'measure q0', 'measure q1', 'measure q2'] 

これを用いて図示します。

import matplotlib.pyplot as plt
qprog1=Qasm()
for i in range(len(SSTR)):
    qprog1.add(SSTR[i])

qprog1.plot()
plt.show()

f:id:puyokw:20180409040419p:plain

他にも, 実機で実行する場合などのチュートリアルが以下に書かれています。
GitHub - QISKit/qiskit-tutorial: A collection of Jupyter notebooks using QISKit

python から実機での実行がまだ上手く行っていないので, そのあたりをもう少し調べていこうと思います。