超入門で、Hello Many Worldsがあったので、それに沿って補足しながら進めていきます。実際には量子回路とNISQの知識がないと厳しいですが、なるべくその辺りをフォローしながら進めます。

https://www.tensorflow.org/quantum/tutorials/hello_many_worlds?hl=ja

Tensorflow Quantumとは？

Tensorflowを使って量子機械学習をするためのライブラリで、主にGoogleの提供するTensorflowという深層学習むけとCirqという量子コンピュータ向けをつなげています。

突如として発表されました。

https://www.tensorflow.org/quantum?hl=ja

インストールです

Google colabの場合には、tensornetwork2を使うようにして、

try:
%tensorflow_version 2.x
except Exception:
pass

#=>TensorFlow 2.x selected.

インストールします。

pip install -q tensorflow-quantum

ツールを読み込みます。

ツールは一番上がtensorflow関連

その次が量子回路シミュレータのcirqとsympyやnumpyなどの数値ライブラリ、

最後にmatplotlibで表示関連をまとめます。

import tensorflow as tf
import tensorflow_quantum as tfq

import cirq
import sympy
import numpy as np

# 表示用
%matplotlib inline
import matplotlib.pyplot as plt
from cirq.contrib.svg import SVGCircuit

目的は「1量子ビットの任意の量子状態を操作する」というのがモチベーションになっています。

しかし、実機マシンではこの量子状態を取得することができませんので、間接的に期待値計算を行います。

本来は人間がこれまでは頭で行ったりしていましたが、機械学習ではかなり複雑な量子状態を使うことがあるので、任意の量子状態人間が考えて作り出すのは現実的でないとのことで、ニューラルネットワークを使って狙った量子状態に持っていけますよというのを示しているのが今回のチュートリアルになります。

量子回路をパラメータ化する

今回はtensorflow-quantumに収録されているGoogleの量子回路シミュレータのCirqを利用します。

https://github.com/quantumlib/Cirq

Cirqでは量子回路を角度を使ったパラメータで準備しますが、通常はパラメータを決めてから実行をすることもありますが、Cirqでは、Sympyのシンボルを利用して角度パラメータをあとで代入できるように準備しています。

今回は量子回路の中で利用される角度パラメータをaとbの二つ準備します。

a, b = sympy.symbols('a b')

量子回路を最初に準備します。今回は２量子ビットを使うので、それを明示します。Cirqでは量子回路の形状を決められますが、今回は2量子ビットなのであまり深く考えなくて平気そうです。

q0にrxゲートを角度aで、q1にryゲートを角度bで適用します。そして、そのあとにq0をコントロールビットに、q1をターゲットビットにしてCNOT回路を適用します。

これによって、２つの量子ビットをもつれさせながら、量子状態を決めていくことができます。

# ２量子ビットの量子回路を準備します。
q0, q1 = cirq.GridQubit.rect(1, 2)

# rxゲートにaの角度をryゲートにbの角度を入れて、それぞれq0とq1の量子ビットに適用します。そのあとに、cnotを適用
circuit = cirq.Circuit(
cirq.rx(a).on(q0),
cirq.ry(b).on(q1), cirq.CNOT(control=q0, target=q1))

SVGCircuit(circuit)

回路が見れました。上記の入力回路が実装されてます。

ここで基本的な確認

量子ビットは|0>か|1>の状態を取り、ベクトルで表されます。

また、量子状態は下記のように重ね合わせで表現されます。

２量子ビット以上はテンソル積を利用します。

後ほど出てくるパウリ行列は下記です。

回転パウリは、

あと、念のためにアダマール行列も、こいつは今回直接は出てきません。

上記の量子回路は念のためおさらいをしておくと、パラメータを使ってRxとRyとCXゲートで表現できます。状態ベクトルはq0とq1共に|0>からスタートしますので、

のテンソル積をとって、

となります。

RxゲートとRyゲートをそれぞれかける際には、これもゲートのテンソル積をとって、それをかければいいです。そのあとにCNOTゲートもきていますので、上の状態ベクトル|00>にそれぞれ行列をかけると、

で計算できます。ちなみに、ゲート同士のテンソルは、

極座標表記

量子ビットの極座標の表記も確認します。

1量子ビットの量子状態で、

で、極座標で表示をすると、

ブロッホ球は二つの直交する純粋状態の重ね合わせで表現できる量子状態を単位球面上に表す表記法。

https://ja.wikipedia.org/wiki/%E3%83%96%E3%83%AD%E3%83%83%E3%83%9B%E7%90%83

上記の角度によって、極座標表示で量子状態が表現できると書いてあります。

任意の状態を実現するには回転ゲートと呼ばれる回転操作に対応する量子ゲートを利用します。X,Y,Z,H,T,Sなどです。回転ゲートに相当する操作はユニタリ行列となっており、|0>から上記の状態を作るためのユニタリ操作として考えられるのが、2*2のユニタリ行列として、

があります。状態の∣0⟩

∣0⟩はベクトル表記をすると、

ですから、上記の行列をかけると、

無事なりました。

上記を直接実現するゲート操作はないので、分解をして実際のゲートから作る必要があります。任意回転のゲートを組み合わせてみましょう。

任意の回転を表現する

基本的にはRzゲートを利用するようです。

マイナスはどちらでも大丈夫だと思います。これを利用することで、任意の量子状態を実現するためのユニタリ行列は、

こちらは後ほど出てきます。

状態ベクトルの確認

上記はsympyのシンボルを使っているので、パラメータを実際に代入して状態ベクトルを確認したい時に見てみます。

cirq.ParamResolver オブジェクトを使ってaとbに角度をそれぞれ代入します。上記の量子回路ではaとbが決まれば量子状態は決まります。

cirq.Simulator を利用して状態ベクトルを取得して表示してみます。

# a=0.5とb=-0.5の時の状態ベクトル
resolver = cirq.ParamResolver({a: 0.5, b: -0.5})
output_state_vector = cirq.Simulator().simulate(circuit, resolver).final_state
output_state_vector
array([ 0.9387913 +0.j , -0.23971277+0.j ,
0. +0.06120872j, 0. -0.23971277j], dtype=complex64)

上記は状態ベクトルと呼ばれる量子状態を表すベクトルですが、一般的にこの量はシミュレータでしか取得できませんので、実際の量子コンピュータで利用したい場合には、状態ベクトルから得られる結果の形を使います。

具体的にはある行列に対する固有値を取得しますが、行列自体はパウリ行列で表現され、上記のパウリ行列を組み合わせたエルミート行列の形になります。

例えば、下の例ではpauliZを利用してZZ0の期待値を求めたり、1/2*Z0+X1の数値を求めたりします。パウリ行列の右下の数は対応する量子ビットの数です。

z0 = cirq.Z(q0)

qubit_map={q0: 0, q1: 1}

z0.expectation_from_wavefunction(output_state_vector, qubit_map).real
0.8775825500488281

このように得られましたが、実際にちょっとこれを計算してみましょう。得られる固有値は、

Ax=λx

の形でかけます。量子コンピュータの場合には、

H∣ψ⟩=λ∣ψ⟩

のようになります。左側から⟨ψ∣をかけると、

⟨ψ∣H∣ψ⟩=⟨ψ∣λ∣ψ⟩=λ⟨ψ∣∣ψ⟩=λ

となりますので、

⟨ψ∣H∣ψ⟩

を計算すれば良いことになります。古典表現をすると、

実際に計算してみます。先ほどの状態ベクトルを取り出し、

v_a = output_state_vector

複素共役を取ります。

v_b = np.conjugate(a)
v_b
array([ 0.9387913 -0.j , -0.23971277-0.j ,
0. -0.06120872j, 0. +0.23971277j], dtype=complex64)

求めたいのは、v_b@Z@v_aになりますので、素直に計算します。ZとIを準備して、

Z = np.array([[1,0],[0,-1]])
I = np.eye(2)

ZとIのテンソル積をとった行列をbとaで挟んで計算します。

v_b@np.kron(Z,I)@v_a
(0.8775825505048687+0j)

先ほどとほぼ同じ値が出ました。

続いて、1/2*Z0+X1を見ると、

z0x1 = 0.5 * z0 + cirq.X(q1)

z0x1.expectation_from_wavefunction(output_state_vector, qubit_map).real
-0.04063427448272705

こちらも素直に計算をすると、pauliXを準備して、

X = np.array([[0,1],[1,0]])

上記のハミルトニアンは線型結合されているので、分解でき、

の計算に落とし込めます。

(v_b@np.kron(Z,I)@v_a)/2 + v_b@np.kron(I,X)@v_a
(-0.04063427246870743+0j)

このようにうまくいきました。下記のようにまとめて計算もできますが、実際の量子回路で合わせるためバラバラにしてます。

v_b@(np.kron(Z,I)/2 + np.kron(I,X))@v_a
(-0.04063427246870744+0j)

量子回路のテンソル化

テンソルフローとオブジェクトを共有するためにTensorflow Quantum (TFQ)は、

tfq.convert_to_tensorの関数を提供していて、Cirqの量子回路をtensorオブジェクトに変換します

# 階数1のテンソルに（ベクトル）
circuit_tensor = tfq.convert_to_tensor([circuit])

print(circuit_tensor.shape)
print(circuit_tensor.dtype)
(1,)
<dtype: 'string'>
# こちらは２つのパウリ演算子を持つ階数1のテンソルに（ベクトル）
pauli_tensor = tfq.convert_to_tensor([z0, z0x1])
pauli_tensor.shape

できました。

TensorShape([2])

こちらは後ほど量子回路の入力値として利用します。

量子回路のバッチ処理

TFQでは量子回路のバッチ処理を提供しています。基本的には今後、量子回路から得られた結果の期待値を計算しますが、cirq.ParamResolversで指定されたパラメータを使ってシミュレートします。

aとbの角度パラメータをバッチ処理で準備します。ここでは、aとbにランダム角度を指定し、それを５セット準備します。

batch_vals = np.array(np.random.uniform(0, 2 * np.pi, (5, 2)), dtype=np.float32)
batch_vals
array([[5.8682933, 5.061834 ],
[3.7570376, 2.8654807],
[4.493578 , 2.910096 ],
[2.7371805, 1.9324272],
[1.313978 , 3.5935752]], dtype=float32)

そして、Cirqの回路に上記の角度パラメータをループで処理します。

cirq_results = []
cirq_simulator = cirq.Simulator()

for vals in batch_vals:
resolver = cirq.ParamResolver({a: vals[0], b: vals[1]})
final_state = cirq_simulator.simulate(circuit, resolver).final_state
cirq_results.append(
[z0.expectation_from_wavefunction(final_state, {
q0: 0,
q1: 1
}).real])

print('cirq batch results: \n {}'.format(np.array(cirq_results)))
cirq batch results:
[[ 0.91515988]
[-0.81651652]
[-0.21706918]
[-0.91933388]
[ 0.25400463]]

期待値の計算はたびたび出るので、TFQで簡略化されます。

tfq.layers.Expectation()(circuit,
symbol_names=[a, b],
symbol_values=batch_vals,
operators=z0)
<tf.Tensor: shape=(5, 1), dtype=float32, numpy=
array([[ 0.9151599 ],
[-0.81651634],
[-0.21706916],
[-0.91933393],
[ 0.2540047 ]], dtype=float32)>

量子古典ハイブリッド最適化

これで期待値の求め方がわかったので、ハイブリッド計算に入れます。量子ニューラルネットのチュートリアルをやってみます。

１量子ビットの量子回路を使って最適化計算をしてみます。

tfq.layers.ControlledPQCでパラメータ化された量子回路を使ってやってみます。

チュートリアルの実装で、全体のアーキテクチャは３つの部分に分かれます。

１、最初の３つのゲートでデータの入力

２、次の３つのゲートで量子状態の操作

３、古典のニューラルネットワークでパラメータの調整

をします。

まずは回路の実装

上記の回路の準備をします。

# 古典ニューラルネットで操作する角度パラメータを準備します。３つあります。
control_params = sympy.symbols('theta_1 theta_2 theta_3')

# そして、１量子ビットを準備し、上記のパラメータをrz,ry,rx回路に当てはめます。
qubit = cirq.GridQubit(0, 0)
model_circuit = cirq.Circuit(
cirq.rz(control_params[0])(qubit),
cirq.ry(control_params[1])(qubit),
cirq.rx(control_params[2])(qubit))

SVGCircuit(model_circuit)

古典コントローラの設定

次にコントローラネットワークを実装します。

# こちらは古典の通常のネットワークです
controller = tf.keras.Sequential([
tf.keras.layers.Dense(10, activation='elu'),
tf.keras.layers.Dense(3)
])

初期状態では値はランダムなので役には立ちません。

controller(tf.constant([[0.0],[1.0]])).numpy()
array([[ 0. , 0. , 0. ],
[-0.36744112, 0.13772306, -0.4180936 ]], dtype=float32)

古典コントローラを量子回路に接続

今回はkerasのモデルとして古典NNのコントローラを量子回路に接続します。

入力は２種類あり、古典コントローラに入れる入力値と量子回路に最初に入れるランダム回路です。 このランダム回路の入力を古典NNが学び、調整をしてZの期待値を学びます。

入力２つを準備します。

前者の入力は量子回路が入力されますが、テンソル化されて、tf.stringで入ってくることに注意します。

# 量子回路に準備する古典NNが修正する量子状態の入力
circuits_input = tf.keras.Input(shape=(),
# The circuit-tensor has dtype `tf.string`
dtype=tf.string,
name='circuits_input')

#古典NNに入力する0か1の数値。
commands_input = tf.keras.Input(shape=(1,),
dtype=tf.dtypes.float32,
name='commands_input')

ハミルトニアンZの期待値を求める処理を記述します。

dense_2 = controller(commands_input)

# Zの期待値を量子回路から求めます
expectation_layer = tfq.layers.ControlledPQC(model_circuit,
# Observe Z
operators = cirq.Z(qubit))
expectation = expectation_layer([circuits_input, dense_2])

あとは、これを全て tf.keras.Model としてまとめます

# The full Keras model is built from our layers.
model = tf.keras.Model(inputs=[circuits_input, commands_input],
outputs=expectation)

アーキテクチャを確認します。 graphviz package　が必要になります。

tf.keras.utils.plot_model(model, show_shapes=True, dpi=70)

入力は古典NNへのコマンド入力が0か1になり、入力１つに対して、出力が３つの角度パラメータに。一方量子回路も入力があります。

データセット

このモデルでは、古典への入力値0か1によって、pauliZの期待値+1か-1が出力されます。

# 古典のニューラルネットに対する入力は0か1
commands = np.array([[0], [1]], dtype=np.float32)

量子回路から出力されるハミルトニアンZの期待値は1か-1

expected_outputs = np.array([[1], [-1]], dtype=np.float32)

このほかに量子回路としての入力があります。

入力回路の準備

入力回路は今回ランダムで準備され、この入力回路を補正するようにpauliZの期待値を求めます。

random_rotations = np.random.uniform(0, 2 * np.pi, 3)
datapoint_circuits = tfq.convert_to_tensor([
cirq.Circuit(
cirq.rx(random_rotations[0])(qubit),
cirq.ry(random_rotations[1])(qubit),
cirq.rz(random_rotations[2])(qubit))
] * 2) # Make two copied of this circuit
datapoint_circuits.shape
TensorShape([2])

トレーニング

それぞれのコマンドに対してのデータポイントの値を確認できます。

model([datapoint_circuits, commands]).numpy()
array([[-0.3738031 ],
[ 0.32185417]], dtype=float32)

早速これをトレーニングしてみます。

optimizer = tf.keras.optimizers.Adam(learning_rate=0.05)
loss = tf.keras.losses.MeanSquaredError()
model.compile(optimizer=optimizer, loss=loss)
history = model.fit(x=[datapoint_circuits, commands],
y=expected_outputs,
epochs=30,
verbose=0)
plt.plot(history.history['loss'])
plt.title("Learning to Control a Qubit")
plt.xlabel("Iterations")
plt.ylabel("Error in Control")
plt.show()

これによって、NNが量子ビットの状態の制御をできていることがわかります。

異なる演算子の固有状態を求める

上記ではpauliZの固有状態は1もしくは0に対しては任意で設定できたので、入力1に対して、+Zの期待値、入力0に対して今度は-Xの固有状態を設定もできるはずです。

ここで、ZとXに対してはハミルトニアンの固有状態の測定回路が異なります。古典NNで測定を含むZとXの期待値を求めるようにトレーニングできます。

今回も期待値を求めますが、パラメータは１つ増えます。

新しいモデルの作成

モデルを新しく作り直します。pauliZとpauliXが同居しているので、測定演算子を規定します。

# Define inputs.
commands_input = tf.keras.layers.Input(shape=(1),
dtype=tf.dtypes.float32,
name='commands_input')
circuits_input = tf.keras.Input(shape=(),
# The circuit-tensor has dtype `tf.string`
dtype=tf.dtypes.string,
name='circuits_input')
operators_input = tf.keras.Input(shape=(1,),
dtype=tf.dtypes.string,
name='operators_input')

コントローラネットワークはほぼ同じです。

# Define classical NN.
controller = tf.keras.Sequential([
tf.keras.layers.Dense(10, activation='elu'),
tf.keras.layers.Dense(3)
])

接続します。

dense_2 = controller(commands_input)

# Since you aren't using a PQC or ControlledPQC you must append
# your model circuit onto the datapoint circuit tensor manually.
full_circuit = tfq.layers.AddCircuit()(circuits_input, append=model_circuit)
expectation_output = tfq.layers.Expectation()(full_circuit,
symbol_names=control_params,
symbol_values=dense_2,
operators=operators_input)

# Contruct your Keras model.
two_axis_control_model = tf.keras.Model(
inputs=[circuits_input, commands_input, operators_input],
outputs=[expectation_output])

データセット

今回はXとZで異なる測定演算子を準備します。

# The operators to measure, for each command.
operator_data = tfq.convert_to_tensor([[cirq.X(qubit)], [cirq.Z(qubit)]])

The command input values to the classical NN.

commands = np.array([[0], [1]], dtype=np.float32)

The desired expectation value at output of quantum circuit.

expected_outputs = np.array([[1], [-1]], dtype=np.float32)

トレーニング

トレーニングをします。

optimizer = tf.keras.optimizers.Adam(learning_rate=0.05)
loss = tf.keras.losses.MeanSquaredError()

two_axis_control_model.compile(optimizer=optimizer, loss=loss)

history = two_axis_control_model.fit(
x=[datapoint_circuits, commands, operator_data],
y=expected_outputs,
epochs=30,
verbose=1)
Train on 2 samples
Epoch 1/30
2/2 [==============================] - 1s 322ms/sample - loss: 0.4949
Epoch 2/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.2670
Epoch 3/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.1252
Epoch 4/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.0595
Epoch 5/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.0293
Epoch 6/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.0112
Epoch 7/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.0027
Epoch 8/30
2/2 [==============================] - 0s 3ms/sample - loss: 4.6870e-04
Epoch 9/30
2/2 [==============================] - 0s 4ms/sample - loss: 3.7340e-04
Epoch 10/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.0013
Epoch 11/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.0040
Epoch 12/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.0087
Epoch 13/30
2/2 [==============================] - 0s 4ms/sample - loss: 0.0133
Epoch 14/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.0156
Epoch 15/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.0151
Epoch 16/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.0126
Epoch 17/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.0092
Epoch 18/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.0059
Epoch 19/30
2/2 [==============================] - 0s 2ms/sample - loss: 0.0033
Epoch 20/30
2/2 [==============================] - 0s 3ms/sample - loss: 0.0016
Epoch 21/30
2/2 [==============================] - 0s 2ms/sample - loss: 6.9725e-04
Epoch 22/30
2/2 [==============================] - 0s 2ms/sample - loss: 2.7992e-04
Epoch 23/30
2/2 [==============================] - 0s 3ms/sample - loss: 1.1527e-04
Epoch 24/30
2/2 [==============================] - 0s 3ms/sample - loss: 6.2470e-05
Epoch 25/30
2/2 [==============================] - 0s 3ms/sample - loss: 5.8352e-05
Epoch 26/30
2/2 [==============================] - 0s 3ms/sample - loss: 8.7101e-05
Epoch 27/30
2/2 [==============================] - 0s 3ms/sample - loss: 1.5469e-04
Epoch 28/30
2/2 [==============================] - 0s 3ms/sample - loss: 2.7070e-04
Epoch 29/30
2/2 [==============================] - 0s 3ms/sample - loss: 4.3545e-04
Epoch 30/30
2/2 [==============================] - 0s 2ms/sample - loss: 6.3206e-04
plt.plot(history.history['loss'])
plt.title("Learning to Control a Qubit")
plt.xlabel("Iterations")
plt.ylabel("Error in Control")
plt.show()

しっかり学習できています。

スタンドアローンでコントローラを確認できます。入力値0と1に対応する出力パラメータは、

controller.predict(np.array([0,1]))
array([[-0.42406332, -1.2505993 , -0.08478868],
[ 0.10877109, -0.05627206, -0.11060521]], dtype=float32)

以上です

今回は簡単なモデルに対して学んでいきました。かなり基本的な内容なのでこれを大きくしてより深いモデルも学習していきましょう。