量子コンピュータを使った機械学習（量子機械学習）にはさまざまなアルゴリズムがあります。
代表的なものとして、量子サポートベクターマシン(QSVM)、量子ニューラルネットワーク(QNN)、量子GAN、量子ボルツマンマシンがありますが、ここでは最もアルゴリズムが理解しやすい QNN を使った検証実験を再現してみました。
再現する実験は論文
Quantum Machine Learning in Finance: Time Series Forecasting
に書かれているものです。
この文献では、BiLSTM という種類のニューラルネットワークと同等の性能を持つ QNN を構築しています。検証する問題としては、金融商品価格の予測（個別企業の株価やビットコインの価格変動予測）を行っています。
機械学習によって株価等を予測する問題は非常に難しく、評価も難しいので、本記事では論文の再現だけでなく、他の条件での実験も行っていきます。
実装方針
量子機械学習の実装には PyTorch と qiskit を採用しました。
qiskit には TorchConnector というクラスがあり、PyTorch の自動微分機能を使ってパラメータの学習ができるらしく、
それを試してみたいという理由です。
特に量子古典ハイブリッドな機械学習を行う場合は、PennyLane や TensorFlow Quantum の方が使いやすいと思われます。
ただし、TensorFlow Quantum は 公式に配布されているバイナリが、Python3.9 と TensorFlow 2.7 までしか対応しておらず、
開発が活発でなさそうに見えます。
また、株価やビットコイン価格のデータは yahoo-finance-api2 を使って取得します。
from yahoo_finance_api2 import share
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
データの取得と確認
論文にあわせて、米 Apple 社(コード:AAPL) の約 30 年分の日足データを取得します。
論文では 10000 データ点ほど取得できていますが、ここでは 7500 データ点ほど得られました。
恐らく、非営業日のデータが含まれていないためと思いますので、気にしないことにします。
stock = share.Share('AAPL')
symbol_data = stock.get_historical(share.PERIOD_TYPE_YEAR, 30, share.FREQUENCY_TYPE_DAY, 1)
df = pd.DataFrame(symbol_data)
df['timestamp'] = pd.to_datetime(df.timestamp, unit = 'ms')
print(len(df))
df.head(10)
7549
| timestamp | open | high | low | close | volume | |
|---|---|---|---|---|---|---|
| 0 | 1993-11-24 14:30:00 | 0.292411 | 0.299107 | 0.291295 | 0.294643 | 90440000 |
| 1 | 1993-11-26 14:30:00 | 0.292411 | 0.294643 | 0.287946 | 0.291295 | 43444800 |
| 2 | 1993-11-29 14:30:00 | 0.287946 | 0.290179 | 0.281250 | 0.283482 | 96712000 |
| 3 | 1993-11-30 14:30:00 | 0.283482 | 0.291295 | 0.281250 | 0.281250 | 112660800 |
| 4 | 1993-12-01 14:30:00 | 0.285714 | 0.287946 | 0.279018 | 0.281250 | 111216000 |
| 5 | 1993-12-02 14:30:00 | 0.283482 | 0.285714 | 0.276786 | 0.283482 | 100654400 |
| 6 | 1993-12-03 14:30:00 | 0.283482 | 0.285714 | 0.276786 | 0.281250 | 120467200 |
| 7 | 1993-12-06 14:30:00 | 0.281250 | 0.290179 | 0.279018 | 0.287946 | 156979200 |
| 8 | 1993-12-07 14:30:00 | 0.285714 | 0.287946 | 0.281250 | 0.287946 | 63851200 |
| 9 | 1993-12-08 14:30:00 | 0.285714 | 0.287946 | 0.281250 | 0.284598 | 39592000 |
raw_data = df['close'].to_numpy()
change_data = []
for i in range(len(raw_data) - 1):
change_data.append((raw_data[i+1] / raw_data[i] - 1))
fig = plt.figure(figsize = (10, 4))
ax1 = fig.add_subplot(1, 2, 1)
ax1.plot(df.timestamp, df.close, '-', markersize = 1)
ax2 = fig.add_subplot(1, 2, 2)
ax2.plot(df.timestamp[1:], change_data, 'o', markersize = 1)
データの前処理
グラフ右の成長率のデータについて、先頭から 4分の3 を学習データ、残りの 4分の1 を検証データとし、
学習データに MinMaxScaler を適用して、0.2〜0.8 の範囲の値になるようにします。
テストデータ、及び予測値は学習データの範囲を超えて、0〜1 の範囲を許容するものとします。（今回使用する QNN の出力値が 0〜1 の範囲になるように実装します）
train_len = int(len(change_data) * 0.75)
scaler = MinMaxScaler(feature_range=(0.2, 0.8))
train_data = scaler.fit_transform(np.reshape(change_data[:train_len], (-1, 1)))
test_data = scaler.fit_transform(np.reshape(change_data[train_len:], (-1, 1)))
モデルの入出力
学習データから連続する 17 日分のデータを抜き出し、1 日目から 16 日目までの値を入力として、
17 日目の値を出力するものとします。
テストも同様に、連続する 16 日分のデータをモデルに与えて、17 日目の値を予測させるようにします。
window_size = 16
# Training data
x_train, y_train = [], []
for i in range(window_size, len(train_data)):
x_train.append(train_data[i-window_size:i, 0])
y_train.append(train_data[i, 0])
np.shape(np.array(x_train))
# Training data
x_train, y_train = [], []
for i in range(window_size, len(train_data)):
x_train.append(train_data[i-window_size:i, 0])
y_train.append(train_data[i, 0])
# numpy arrayに変換
x_train, y_train = np.array(x_train), np.array(y_train)
# (batch, series, feature)
x_train = np.reshape(x_train, (x_train.shape[0], x_train.shape[1], 1))
y_train = np.reshape(y_train, (y_train.shape[0], 1))
# Test data
x_test, y_test = [], []
for i in range(window_size, len(test_data)):
x_test.append(test_data[i-window_size:i, 0])
y_test.append(test_data[i, 0])
# numpy arrayに変換
x_test, y_test = np.array(x_test), np.array(y_test)
# (batch, series, feature)
x_test = np.reshape(x_test, (x_test.shape[0], x_test.shape[1], 1))
y_test = np.reshape(y_test, (y_test.shape[0], 1))
BiLSTM
PyTorch の BiLSTM を使って、論文に記載の通りに予測モデルを構築します。
論文では 300 epochs も学習を行ったようですが、30 epochs 程度で十分収束していると判断しました。
パラメータ数は論文の値より少し多いですが、176,944 です。
import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F
num_epochs = 30
batch_size = 10
def make_dataset(x_train, y_train, batch_size):
x_batch = []
y_batch = []
for _ in range(batch_size):
idx = np.random.randint(0, len(x_train) - 1)
x_batch.append(x_train[idx])
y_batch.append(y_train[idx])
return torch.tensor(np.array(x_batch, dtype=np.float32)), torch.tensor(np.array(y_batch, dtype=np.float32))
class BiLSTM(nn.Module):
def __init__(self):
super(BiLSTM, self).__init__()
self.bd_seq = nn.LSTM(1, 128, bidirectional=True, batch_first=True)
self.bd_sin = nn.LSTM(128, 32, bidirectional=True, batch_first=True)
self.bd_1 = nn.LSTM(128, 1, bidirectional=True, batch_first=True)
self.bd_2 = nn.LSTM(32, 1, bidirectional=True, batch_first=True)
self.activation = nn.Tanh()
def forward(self, x):
bd_lstm, _ = self.bd_seq(x)
act_bd = self.activation(bd_lstm)
act_bd = act_bd[:, :, :128]
sin_lstm, _ = self.bd_sin(act_bd)
act_sin = self.activation(sin_lstm)
act_sin = act_sin[:, :, :32]
bd1_lstm, _ = self.bd_1(act_bd)
bd2_lstm, _ = self.bd_2(act_sin)
act_bd2 = self.activation(bd2_lstm)
out = torch.add(bd1_lstm, act_bd2)
return out[:, -1, 0].reshape([out.shape[0], 1])
lstm_model = BiLSTM()
lstm_optimizer = optim.Adam(lstm_model.parameters())
criterion = nn.MSELoss()
for epoch in range(num_epochs):
# training
running_loss = 0.0
training_accuracy = 0.0
for i in range((int)(x_train.shape[0] / batch_size)):
lstm_model.zero_grad()
x, y = make_dataset(x_train, y_train, batch_size)
output = lstm_model(x)
loss = criterion(output, y)
loss.backward()
lstm_optimizer.step()
running_loss += loss.item()
print('%d loss: %.6f' % (epoch + 1, running_loss))
予測
テストデータを使って、モデルに予測値を出力させます。評価指標として、予測と実際の値との間の平均二乗誤差を計算します。
with torch.no_grad():
y_pred = lstm_model(torch.tensor(x_test, dtype=torch.float32))
mse = F.mse_loss(y_pred, torch.tensor(y_test, dtype=torch.float32))
print(mse)
tensor(0.0053)
MinMaxscaler によって逆変換して、元の数値尺度で評価すると誤差の値が論文に記載の値に近づいて、0.0009 になりました。とりあえずこれを基準として QNN を使った場合との比較を行います。
（この予測精度の絶対評価は要求仕様に依存するので、ここでは深く言及しませんが、あまり良いものではないと思います。）
y_test_orig = np.ravel(scaler.inverse_transform(np.reshape(y_test, (-1, 1))))
y_pred_orig = np.ravel(scaler.inverse_transform(np.reshape(y_pred.numpy(), (-1, 1))))
mse_orig = F.mse_loss(torch.tensor(y_pred_orig), torch.tensor(y_test_orig))
print(mse_orig)
tensor(0.0009, dtype=torch.float64)
QNN
QNN は論文に合わせて、入力をエンコードする 16 量子ビットと、結果取得用の 1 量子ビットの合わせて 17 量子ビットを使用します。入力データ 1 個に対して 1 量子ビットを割り当て、Rx ゲートの回転角を使ってエンコードします。
ニューラルネットワーク部分は、以下の図（図は同論文から引用）のように XXゲート→ZZゲート→YYゲートからなる 3 層の量子回路を 2 個結合させた、6 層の量子回路となっています。
qiskit で EstimatorQNN を構築して、それを TorchConnector に渡すだけで PyTorch のモデルが出来上がります。パラメータ数は 96 です。
from qiskit.circuit import Parameter
from qiskit import QuantumCircuit
from qiskit.circuit import ParameterVector
from qiskit.quantum_info import SparsePauliOp
inputs = ParameterVector("input", window_size)
params_XX1 = ParameterVector("XX1", window_size)
params_ZZ1 = ParameterVector("ZZ1", window_size)
params_YY1 = ParameterVector("YY1", window_size)
params_XX2 = ParameterVector("XX2", window_size)
params_ZZ2 = ParameterVector("ZZ2", window_size)
params_YY2 = ParameterVector("YY2", window_size)
qc = QuantumCircuit(window_size + 1)
# エンコード
for i in range(window_size):
qc.rx(inputs[i], i+1)
# 推論
for i in range(window_size):
qc.rxx(params_XX1[i], 0, i+1)
for i in range(window_size):
qc.rzz(params_ZZ1[i], 0, i+1)
for i in range(window_size):
qc.ryy(params_YY1[i], 0, i+1)
for i in range(window_size):
qc.rxx(params_XX2[i], 0, i+1)
for i in range(window_size):
qc.rzz(params_ZZ2[i], 0, i+1)
for i in range(window_size):
qc.ryy(params_YY2[i], 0, i+1)
# 観測
observable = SparsePauliOp(["I"*window_size + "Z"])
from qiskit_machine_learning.neural_networks import EstimatorQNN
from qiskit_machine_learning.connectors import TorchConnector
estimator_qnn = EstimatorQNN(
circuit=qc, observables=observable, input_params=inputs,
weight_params=params_XX1.params + params_ZZ1.params + params_YY1.params + params_XX2.params + params_ZZ2.params + params_YY2.params
)
qnn_model = TorchConnector(estimator_qnn)
# Define optimizer and loss
qnn_optimizer = optim.Adam(qnn_model.parameters())
qnn_loss = nn.MSELoss()
# Start training
qnn_model.train()
for epoch in range(num_epochs):
# training
running_loss = 0.0
for i in range((int)(x_train.shape[0] / batch_size)):
qnn_model.zero_grad()
x, y = make_dataset(x_train, y_train, batch_size)
x *= np.pi
x = x.reshape([batch_size, window_size])
output = qnn_model(x)
# output は [-1, 1] なので、[0, 1] に変換してからlossを計算する
loss = qnn_loss((output + 1.0)/2.0, y)
loss.backward()
qnn_optimizer.step()
running_loss += loss.item()
print('%d loss: %.6f' % (epoch + 1, running_loss))
しかし、いつまで待っても学習が終わらない問題に直面しました。
デバッグをしてみたところ、backward() の計算に 1 回あたり数分もの時間がかかっていて、全く学習が進んでくれません。
PyTorch の機能を使いたかったのですが、残念な結果になりました。
TensorFlow Quantum で実装
今度は論文と同じく、TensorFlow Quantum を使った実装で実験を行ってみます。
こちらの実装では、学習データとテストデータをパラメータで与えるのではなく、固定の回転角を持つ量子回路を学習データの数だけ作成します。
Rx ゲートの回転角は入力データ$\times \pi$ とするのが望ましいと思われますが、ここでは $\pi$ を掛けていません。（その方がわずかに二乗誤差が小さかったです）
def convert_to_circuit(data):
"""データエンコード部"""
values = np.ndarray.flatten(data)
qubits = cirq.GridQubit.rect(4, 4)
circuit = cirq.Circuit()
for i, value in enumerate(values):
if value:
circuit.append(cirq.rx(value)(qubits[i]))
return circuit
x_train_circ = [convert_to_circuit(x) for x in x_train]
x_test_circ = [convert_to_circuit(x) for x in x_test]
x_train_tfcirc = tfq.convert_to_tensor(x_train_circ)
x_test_tfcirc = tfq.convert_to_tensor(x_test_circ)
続いて、QNN レイヤーを実装します。qiskit と比べて実装方法が分かりにくいです。
class CircuitLayerBuilder():
def __init__(self, data_qubits, readout):
self.data_qubits = data_qubits
self.readout = readout
def add_layer(self, circuit, gate, prefix):
for i, qubit in enumerate(self.data_qubits):
symbol = sympy.Symbol(prefix + '-' + str(i))
circuit.append(gate(qubit, self.readout)**symbol)
def create_quantum_model():
data_qubits = cirq.GridQubit.rect(4, 4)
readout = cirq.GridQubit(-1, -1)
builder = CircuitLayerBuilder(data_qubits = data_qubits,
readout=readout)
circuit = cirq.Circuit()
builder.add_layer(circuit, gate = cirq.XX, prefix='xx1')
builder.add_layer(circuit, gate = cirq.ZZ, prefix='zz1')
builder.add_layer(circuit, gate = cirq.YY, prefix='yy1')
builder.add_layer(circuit, gate = cirq.XX, prefix='xx2')
builder.add_layer(circuit, gate = cirq.ZZ, prefix='zz2')
builder.add_layer(circuit, gate = cirq.YY, prefix='yy2')
return circuit, cirq.Z(readout)
model_circuit, model_readout = create_quantum_model()
以降は、Keras の流儀に従ってモデルパラメータを学習していきます。
# Build the Keras model.
model = tf.keras.Sequential([
# The input is the data-circuit, encoded as a tf.string
tf.keras.layers.Input(shape=(), dtype=tf.string),
# The PQC layer returns the expected value of the readout gate, range [-1,1].
tfq.layers.PQC(model_circuit, model_readout),
])
model.compile(
loss='mse',
optimizer=tf.keras.optimizers.Adam(),
metrics='mse')
qnn_history = model.fit(
x_train_tfcirc, y_train,
batch_size=32,
epochs=10,
verbose=1,
validation_data=(x_test_tfcirc, y_test))
こちらの実装では、難なく 10 分程度で学習が進みました。
BiLSTM よりも少し誤差が大きく、0.0014 になっていますが、同程度の予測ができているとみなして良いのではないでしょうか。
y_pred = model.predict(x_test_tfcirc)
y_test_orig = np.ravel(scaler.inverse_transform(np.reshape(y_test, (-1, 1))))
y_pred_orig = np.ravel(scaler.inverse_transform(np.reshape(y_pred, (-1, 1))))
mse_orig = tf.keras.losses.MeanSquaredError()(y_test_orig, y_pred_orig).numpy()
print(mse_orig)
0.0013960395
他の条件で予測
これだけの実験では面白くないので、以下の条件・データでも予測を行ってみました。金額の推移を表すグラフと、予測誤差のみを示します。
-
量子回路の層を半分にして Apple 社の株価を予測
TorchConnector で学習を進めるための対処法が分からないので、回路の深さを半分（各量子ビットについて XX, ZZ, YY ゲート演算 1 回のみ）の小さい回路で学習を進めることにしました。それでも、30 epochs の学習に何日もの時間を要しました。
結果は、平均二乗誤差が 0.0243 と、あまり奮いませんでした。 -
米ドル / ビットコインのレートを予測
yahoo-finance-api2 での取得方法がわからなかったので、事前に 米 Yahoo Finance から csv を取得しました。
2014年9月以降のデータで、データ件数は 3352 件です。試した中では最も変化が激しいデータです。
-
円 / 米ドルのレートを予測
1996年11月以降のデータで、データ件数は 7017 件です。
-
日経平均株価を予測
1996年11月以降のデータで、データ件数は 7358 件です。
結果
使用したデータと予測結果をまとめると以下の表のようになります。表の平均二乗誤差は MinMaxScaler で変換する前の値です。
変化率の平均値が大きいアップルとビットコインの誤差が大きく、逆に変化率の平均が小さい米ドルと日経平均の誤差が小さくなっていることから、どのデータも同程度に予測できているのではないかと思われます。
| データ | アップル株価 | 米ドル / ビットコイン | 円 / 米ドル | 日経平均 |
|---|---|---|---|---|
| 件数 | 7549 | 3352 | 7017 | 7358 |
| 変化率の平均 | +0.12% | +0.20% | +0.0063% | +0.019% |
| 平均二乗誤差（BiLSTM） | 0.0009 | 0.0014 | 0.000085 | 0.00019 |
| 平均二乗誤差 (QNN) | 0.0014 | 0.0022 | 0.000172 | 0.00023 |
まとめ
本記事では QNN を使って、実データを使った時系列予測を行いました。
結果、96 パラメータの QNN で LSTM と同じような時系列予測モデルを、少ないパラメータで作ることができました。（LSTMのパラメータをもっと少なくしても同等の結果になるかもしれませんので、96 パラメータでも良かったという点だけ言及します。）
実際に試してみて、実装が容易で何らかの意味のある結果が出るという点で、量子コンピュータを使った実装事例として分かりやすく、将来の量子コンピュータの可能性を一般の方に認識してもらうのに適したアルゴリズムだと感じました。