Qiskitで量子SVMを実装して性能評価してみた

はじめに

量子コンピュータ勉強中の機械学習エンジニアとして量子機械学習はぜひ学びたい。
おっ、量子Support Vector Machine(以下SVM)あるじゃん、ここから始めよう。
そう思ってひとまずQiitaを探したけど、まだこの辺は充実してない(2020/05時点)。
そもそもQiitaに限らず量子SVMの日本語の記事が少なく、あってもただ英語リソース翻訳したくらい。
実装して回してみた系の記事はほぼ皆無。

筆不精ながら一端のエンジニアとしてQiitaは書いてみようという思いはあったので、
書くなら少しでも喜んでくれる人がいてくれそうな、情報リソース少ない量子SVMでと本記事を書くに至りました。

読んだ論文

量子機械学習としてSupport Vector Machineの量子版をIBMの量子コンピュータ実機で実装したという論文です。2019年にNatureに掲載されました。
Supervised learning with quantum enhanced feature spaces
https://arxiv.org/pdf/1804.11326.pdf

紹介する論文では2量子ビットで量子特有の特徴量空間を作り二値分類を実現してます。

論文では量子教師あり学習として以下の２つ手法を紹介しています。
(1)Quantum Variation Classifier(変分量子回路)
　　　・特徴量空間マッピングを量子状態空間で実行
　　　・量子回路で予測まで行う(パラメータ最適化は古典)
(2)Quantum Kernel Estimation
　　　・カーネル関数を量子状態の内積として量子回路で再現
　　　・最終的な予測は古典SVMで実行

２つの手法のスキーム上の違いは以下のようになります。

前者の解説はQunaSysが運営するメディアの記事が分かりやすいです。

本記事では後者に数式の理解と実装の観点を加えて語ってみます。
ちなみにSVMを量子版にするというアイデア自体は2014年時点で提案はされていました(その論文のアルゴリズムはNISQ用ではなく、前述の論文とは若干とは若干異なります)。この実装のチュートリアルはIBM Qiskitのgithub上にもあります。

これらの実装は量子コンピュータ界隈でいうところの「NISQ(Noisy Intermediate Scale Quantum computer)デバイス」の応用になります。これは、各量子ビットに対する雑音を十分に取り除くことができない環境下での古典-量子ハイブリッドアルゴリズムということです。

古典(従来)のSVMおさらい

ラベルを分類する境界面を探すのがSVMです。

(注)
以下、簡単のためハードマージンを想定した説明です。ハードマージンとは境界面で完全にラベル毎に分離できる状況(線形分離可能)です。一方で、どんな境界面を作成しても片方のラベルのグループにもう一方のラベルが混じってしまう状況(線形分離不可能)のことをソフトマージンと呼びます。その場合はスナッグ変数という変数を導入して制約条件を緩めて分類をします。論文での量子SVMもソフトマージンを想定した問題ですが、本記事の説明ではハードマージンで量子カーネルの導入を説明するに留めます。

データセットのうち境界面と最も近い点$\boldsymbol{x}_i$と境界面の距離をマージンといい、このマージンを最大化することで境界面が決定します。

境界面はベクトル$\boldsymbol{w}$とスカラー$b$を用いて、
$$\boldsymbol{w}^T \cdot \boldsymbol{x}_i + b = 0 $$

と表現されます。
マージン$M$は、データ点を$\boldsymbol{x}_i$として次のように書けます(高校で習った点と直線の距離ですね)。
$$ M = \min_i \frac{|\boldsymbol{w}^T \cdot \boldsymbol{x}_i + b|}{|\boldsymbol{w}|}$$$\boldsymbol{w}$と$b$には定数倍の任意性があるので、$M=\frac{1}{|\boldsymbol{w}|}$となるように設定します。$$ \min_i |\boldsymbol{w}^T \cdot \boldsymbol{x}_i + b| = 1 \ \ \ (制約条件)$$
さらに、$t_i=\pm1$($\boldsymbol{x}_i$がラベルAに属するときに$t_i= 1$,ラベルBに属するときに$t_i=−1$)を導入します。このとき制約条件が次のように書き換えられます。
$$
\begin{aligned}
&\min_i |\boldsymbol{w}^T \cdot \boldsymbol{x}_i + b| = \min_i t_i (\boldsymbol{w}^T \cdot \boldsymbol{x}_i + b) = 1 \\
& \quad \iff t_i (\boldsymbol{w}^T \cdot \boldsymbol{x}_i + b) \ge 1 \quad \text{for} \ i = 1, 2,\dots,N
\end{aligned}
$$
この条件下で$M=\frac{1}{|\boldsymbol{w}|}$を最大化していきます。

「制約条件の下で目的関数を最大化/最小化する」というときはラグランジ未定係数法の出番です。

$$ L = \frac{1}{2}|\boldsymbol{w}|^2-\sum_i \alpha_i [t_i (\boldsymbol{w}^T\cdot \boldsymbol{x}_i+b)-1]
$$
と新たなパラメータ$\alpha_i$を用いてラグランジアンを定義します。
ラグランジアンによる変数変換で
$$\frac{\partial L}{\partial w}=\frac{\partial L}{\partial b}=0$$
という条件が加わり、$L$が次のように書き換えられます。
$$
L' = \sum_i \alpha_i -\frac{1}{2}\sum_{i,k} \alpha_i \alpha_k t_i t_k \ \boldsymbol{x}_i^T \cdot \boldsymbol{x}_k
$$
これで$\alpha_i \ge 0$と$\sum_i\alpha_i t_i = 0$の制約の下で$L'$を最適化する問題に帰結することができました。

ここまでが線形分離可能な問題での場合でした。線形分離できない問題は高次元空間に転写(特徴量マッピング)して考えます。そのためには上記の議論を$\boldsymbol{x} \rightarrow \phi(\boldsymbol{x})$として展開します。
すると$L'$は次のように書けます。
$$
L' = \sum_i \alpha_i -\frac{1}{2}\sum_{i,k} \alpha_i \alpha_k t_i t_k \ \phi(\boldsymbol{x}_i)^T \cdot \phi(\boldsymbol{x}_k)
$$
この$\phi(\boldsymbol{x}_i)^T \cdot \phi(\boldsymbol{x}_k)$が分かると最適化ができます。しかし$\phi(\boldsymbol{x})$の次元が大きいと計算するのが大変になります。(一般的にマッピングする前よりも次元が大きくなります。)
ここでデータと相性の良い特徴量マッピング$\phi(\boldsymbol{x}_i)$を探索してデータセット間の考えうる全組み合わせの内積を計算するのと、相性の良いカーネル$K(\boldsymbol{x}_i,\boldsymbol{x}_k)=\phi(\boldsymbol{x}_i)^T \cdot \phi(\boldsymbol{x}_k)$を考えることは等価です。
なので、カーネル関数を先に決めてしまえば巨大次元で特徴量の内積を直接計算しないで済みます。これをカーネルトリックと呼びます。

よく使われるカーネル関数として次のようなものがあります。

1. 多項式カーネル　$K(\boldsymbol{x}_i,\boldsymbol{x}_k)= (\boldsymbol{x}_i^T \cdot \boldsymbol{x}_k+c)^d$

　・d次元以下の全種類の単項式を各成分に持つ特徴量ベクトル$\phi$が対応します。
2. ガウスカーネル $K(\boldsymbol{x}_i,\boldsymbol{x}_k)= exp(-\frac{\left|\boldsymbol{x}_i - \boldsymbol{x}_k \right|^2}{2\sigma^2})$
　　・無限次元の特徴量ベクトル$\phi$を用いていることに等価
　　・その中に有用な特徴量も含まれていると期待できる
3. シグモイドカーネル $K(\boldsymbol{x}_i,\boldsymbol{x}_k)= tanh(b \boldsymbol {x}_i^T \cdot \boldsymbol{x}_k+c)$
・ニューラルネットワークと表層的に類似
　
次節で詳しく説明しますが、量子SVMでは量子カーネル関数 $K(\boldsymbol{x}_i,\boldsymbol{x}_k)=\left|\left.\langle \phi(\boldsymbol{x}_i)\right|\phi(\boldsymbol{x}_k)\rangle\right|^2$を使います。

量子SVMの概要

お待ちかねの量子SVMです。
古典では表現できないマッピングを考えれば、量子の優位性が期待されます。
そうして考えられたのがヒルベルト空間を利用したマッピングです。

(論文の図を引用)

*図は1量子ビットの場合で、その場合の対応するマッピングはBloch球になります。
量子ビットが増えるともっと複雑な特徴量空間になり、描画できません。

その特徴量マッピングのカーネルに対応するのが量子状態の内積$K(\boldsymbol{x}_i,\boldsymbol{x}_k)=\left|\left.\langle \phi(\boldsymbol{x}_i)\right|\phi(\boldsymbol{x}_k)\rangle\right|^2$です。

このアイデアを量子コンピュータとして実装するには、量子回路の初期状態$|0>$に何らかのゲート$U(\boldsymbol{x})$を作用させて$\phi(\boldsymbol{x})$を作る必要があります。例えば、特徴量が2つの場合、入力データ$\boldsymbol{x}=(x_1,x_2)$に対応する量子状態をゲート$U(\boldsymbol{x})$を用いて、
$$
U(\boldsymbol{x})|0\rangle = |\phi(\boldsymbol{x})\rangle
$$
と記述できます。
この$U(\boldsymbol{x})$は2次オーダーまでの展開で次のような形で書くことができます。

$$
U(\boldsymbol{x})=\left(\exp^{i\left(x_1Z_1+x_2Z_2+(\pi-x_1)(\pi-x_2)Z_1Z_2\right)}H^{\otimes2}\right)^l
$$
ここで$l$は回路の深さ(depth)と呼ばれる量で、増やすことで精度が向上するとされています。
回路で表現すると以下のような形になります。

ちなみに特徴量が3つの場合に対応する回路図は以下となります。

この回路で求まるカーネル関数$K(\boldsymbol{x}_i,\boldsymbol{x}_k)$の回路出力と理論値は論文によると概ね一致する結果となっています。

(論文の図を一部加工)

実装してみた

ここからはQiskitの量子SVMのモジュールを使って二値分類を実装してみます。
一連のまとまったjupyter notebookでのコードをみたい場合はこちらを参照ください。

(注意: バージョンによって以下コマンドを入れないとエラーが出ることがあります)

!pip install qiskit-aqua[cvx]==0.7.1

バージョン確認

まずqiskitは以下のバージョンで実行しました。

print(qiskit.__qiskit_version__)

{'qiskit-terra': '0.14.1', 'qiskit-aer': '0.5.1', 'qiskit-ignis': '0.3.0', 'qiskit-ibmq-provider': '0.7.1', 'qiskit-aqua': '0.7.1', 'qiskit': '0.19.2'}

同じバージョンを使いたい場合は

pip install qiskit==0.19.2

とするだけで十分です。

moduleのインポート

基本的なmoduleに加えてqiskitの以下moduleをインポートします。

import os
import time
import pickle
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import matplotlib as mpl
from sklearn import preprocessing
from tqdm.notebook import tqdm  #　jupyter用
# from tqdm import tqdm  # jupyter以外用

from qiskit import BasicAer
from qiskit.circuit.library import ZZFeatureMap
from qiskit.aqua import QuantumInstance
from qiskit.ml.datasets import ad_hoc_data, breast_cancer
from qiskit.aqua import aqua_globals
from qiskit.aqua.utils import split_dataset_to_data_and_labels, map_label_to_class_name
from qiskit.aqua.algorithms import SklearnSVM
from qiskit.aqua.algorithms import QSVM
from qiskit.circuit.library import ZZFeatureMap

パラメータ設定・関数定義

データセットの指定や、機械学習モデルのハイパーパラメータを設定します。

# 訓練・テストデータの設定
feature_dim = 2  # 特徴量の数
training_dataset_size = 20
testing_dataset_size = 10

# 予測データの設定
size = 50  # 50x50のメッシュを分類によって色分けし、擬似的に境界面を可視化する
mesh_list = [[2*i/size-1, 2*j/size-1] for i in range(size+1) for j in range(size+1)]

# 量子コンピュータのパラメータ設定
shots = 1024
seed = 10598

モデル評価を行うための関数を定義します。２値分類モデルの評価として精度(Accuracy)、適合率(Precision)、再現率(Recall)、特異度(Specificity)、F値(F1-Measure)を計算しました。それぞれの意味についてはこちらを参考にしてください。

def eval_model(pred_label, print_mode=True):
    test_label = ["A"]*testing_dataset_size + ["B"]*testing_dataset_size

    accuracy = sum([x == y for x, y in zip(pred_label, test_label)])/len(test_label)
    precision = sum([x == y for x, y in zip(pred_label, test_label) if x == "A"])/sum([x == "A" for x in pred_label])
    recall = sum([x == y for x, y in zip(pred_label, test_label) if y == "A"])/sum([y == "A" for y in test_label])
    specificity = sum([x == y for x, y in zip(pred_label, test_label) if y == "B"])/sum([y == "B" for y in test_label])
    f1 = 2*recall*precision/(recall + precision)

    eval_dict = {"accuracy": accuracy, "precision": precision, "recall": recall, "specificity": specificity, "F1-measure":f1}
    if print_mode:
        print("result: ", pred_label)
        print("truth : ", test_label)
        print(eval_dict)
    else:
        return eval_dict

結果を可視化する関数を定義します。
heatmap関数では50x50のメッシュをヒートマップとして色分けしたものをプロットします。
scatter_data関数では訓練データ・テストデータをヒートマップ上に色分けした上でプロットします。

def heatmap(pred_list, size=50):
    mat = np.flipud(pred_list.reshape(size+1, size+1, order='F'))
    centers = [-1, 1, -1, 1]
    dx, = np.diff(centers[:2])/(size)
    dy, = -np.diff(centers[2:])/(size)
    extent = [centers[0]-dx/2, centers[1]+dx/2, centers[2]+dy/2, centers[3]-dy/2]
    cmap = mpl.colors.ListedColormap(['orange', 'cyan'])
    # ヒートマップ表示
    plt.imshow(mat, interpolation='nearest', vmin=0, vmax=1, cmap=cmap, extent=extent)


def scatter_data(train_for_pred, test_for_pred, train_result, test_result,
                 yshift=-0.155, print_index=False):
    dataset_dict = {"train": train_for_pred, "test": test_for_pred}
    result_dict = {"train":train_result, "test":test_result}
    marker_dict = {"train": "o", "test": "s"}

    for data_type in ["train", "test"]:
        data_num_half = int(len(dataset_dict[data_type])/2) # ラベルA/Bのデータ数が1:1と仮定とする
        for label in ["A", "B"]:  
            if label == "A":
                (plot_data, color) = dataset_dict[data_type][:data_num_half], "red"
            elif label == "B":
                (plot_data, color) = dataset_dict[data_type][data_num_half:], "blue"
            plt.plot(plot_data[:,0], plot_data[:,1], marker_dict[data_type], color=color, markersize=10)

        # 誤分類を×マークでプロット
        for i, pred_label in enumerate(result_dict[data_type]):
            if (i < data_num_half and pred_label != 0)\
              or (i >= data_num_half and pred_label != 1): 
                if print_index:
                    plt.text(dataset_dict[data_type][i][0], dataset_dict[data_type][i][1], str(i), 
                            color="white", size=15, fontweight='bold')
                # ↓ x方向は自動、y方向は手動にてプロットの位置の微調整が現状ベスト
                plt.text(dataset_dict[data_type][i][0], dataset_dict[data_type][i][1] + yshift, "×", 
                         horizontalalignment='center', color="white", size=15, fontweight='bold')
        plt.axis('off')
        plt.title("Classification Boundary", size=15)

乳癌データセット(Breast Cancer Data)で機械学習

機械学習ではお馴染みの乳癌データセットを用いて分類器の構築を行います。
ここではQiskitのインターフェースを通じてscikit-learnの乳癌データセットを取得しています。
後述しますが、(ほぼ)同じ引数の与え方でQiskitから量子SVM用のデータセットを取得することができます。

sample_Total, training_input, test_input, class_labels = breast_cancer(
    training_size=training_dataset_size, test_size=testing_dataset_size,
    n=feature_dim, plot_data=True  # ,gap=0.3 # breast_cancer使用時はgapパラメータを削除
)

実行すると上記のような図がプロットされます。青色とオレンジ色でデータのラベルを表しています。

　古典SVM

以下を実行するとSVMのカーネルが表示され、テストデータの精度が表示されます。
ここで使用しているmoduleもQiskitを通してscikit-learnのSVMを呼び出しており、量子版と同じ引数の与え方で実行できるようになっています。

result = SklearnSVM(training_input, test_input, mesh_list).run()
print("kernel matrix during the training:")
kernel_matrix = result['kernel_matrix_training']
img = plt.imshow(np.asmatrix(kernel_matrix), interpolation='nearest', origin='upper', cmap='bone_r')
plt.show()

print("testing success ratio: ", result['testing_accuracy'])

ここで、訓練データとテストデータを古典SVMがどのように分類したか可視化するために、それぞれを予測データとしてSVMに分類させます(正解ラベルはマスキングしています)。
その上でSVMによるテストデータの分類能力をeval_model関数で評価します。

# モデル評価
# 誤分類チェックのためにtraining/testデータを予測データとして利用
(train_for_pred, _), _ = split_dataset_to_data_and_labels(training_input)
(test_for_pred, _), _ = split_dataset_to_data_and_labels(test_input)
train_result = SklearnSVM(training_input, test_input, train_for_pred).run()
test_result = SklearnSVM(training_input, test_input, test_for_pred).run()

eval_model(test_result["predicted_classes"])

result: ['B', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B']
truth : ['A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B']
{'accuracy': 0.85, 'precision': 1.0, 'recall': 0.7, 'specificity': 1.0, 'F1-measure': 0.82}

ヒートマップでSVMの境界面を可視化します。
赤/オレンジがA、青/水色がBのラベルとします。
●が訓練データで、■がテストデータです。
(Accuracyが0.85のわりに誤分類が目立つのは訓練データで誤分類が多いがテストデータをうまく分類できてしまったからです)

plt.figure(figsize=(7, 7))               
heatmap(result["predicted_labels"])
scatter_data(train_for_pred, test_for_pred, train_result["predicted_labels"], test_result["predicted_labels"],yshift=-0.014)
plt.show()
# 赤がA、青がBのラベル。●が訓練データで、■がテストデータ

量子SVM

同様の手順で量子SVMで分類器を作りデータセットを分類していきます。
回路の深さはデフォルト値の2で試しました。
量子コンピュータは実機でなくシミュレータ(qasm_simulator)を使用します(実機を使うと一般ユーザーは一撃で１日の使用上限に達します)。
ここの処理は私のjupyter環境で4分ほどかかりました。

%%time
backend = BasicAer.get_backend('qasm_simulator')
feature_map = ZZFeatureMap(feature_dim, reps=2)
svm = QSVM(feature_map, training_input, test_input, None)  # the data for prediction can be fed later.
svm.random_seed = seed
quantum_instance = QuantumInstance(backend, shots=shots, seed_simulator=seed, seed_transpiler=seed)
result = svm.run(quantum_instance)

CPU times: user 40.2 s, sys: 1.67 s, total: 41.9 s
Wall time: 3min 59s

SVMと同様のコードで、量子カーネル行列とモデル評価を表示できます。

result: ['A', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'A', 'B', 'B', 'B', 'B', 'A', 'A', 'B', 'B', 'B', 'B']
truth : ['A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B']
{'accuracy': 0.8, 'precision': 0.8, 'recall': 0.8, 'specificity': 0.8, 'F1-measure': 0.8}

量子SVMの境界面を可視化するには、数時間かかります。
jupyter上で実行するよりはターミナル上で実行する方が早いので必要な処理を.pyファイルに書き出します。
今回は一度に予測する点を17個、ループの数を51として、3つのpredict_xxx.pyファイルに分割します。
(51×51のメッシュを作ったので、51×51=3×51×17と分けました)

# 大量の予測データセットに対して実行する場合、terminal上で複数タブにて同時に実行する方が(若干)早い
# ただし、並列化するほど1つあたりの処理が大幅に遅くなるので、劇的な高速化はできない
(input_file_num, iter_num, input_size) = (3, 51, 17)
mesh_list_tmp = np.array(mesh_list).reshape(input_file_num, iter_num, input_size, feature_dim)


code_string = '''
import os
import pickle
from tqdm import tqdm

dir = os.path.dirname(os.path.abspath(__file__))
with open(dir+"/../model/model_qsvm_br_cn.pkl", "rb") as f:
    svm = pickle.load(f)

with open(dir+"/../input/input_INDEX.pkl", "rb") as f:
    input_array = pickle.load(f)

list_final = []
for epoch, epoch_array in enumerate(tqdm(input_array)):
    pred_tmp = svm.predict(epoch_array)
    print("epoch ", epoch, " has done")
    list_final.append(list(pred_tmp))
print(list_final)

with open(dir+"/../output/output_br_cn_INDEX.pkl", "wb") as f:
    pickle.dump(list_final, f)
'''

for path in ['./qsvm_terminal/' + folder for folder in ["input", "output", "script", "model"]]:
    os.makedirs(path, exist_ok=True)

with open("qsvm_terminal/model/model_qsvm_br_cn.pkl", "wb") as f:
    pickle.dump(svm, f)

for i in range(input_file_num):
    # inputはbreast_cancerとad_hoc_data共通でOK
    with open("qsvm_terminal/input/input_{}.pkl".format(i), "wb") as f:
        pickle.dump(mesh_list_tmp[i], f)

    with open("qsvm_terminal/script/predict_br_cn_{}.py".format(i), "w") as f:
        f.write(code_string.replace("INDEX", str(i)))

できたPythonファイルを、ターミナルで実行します。
ファイルはjupyterのディレクトリから見て、./qsvm_terminal/script/に書き込まれます。
ターミナル上で複数タブを開き各タブでpython predict_br_cn_1.pyのように最後の添字を変えて実行してください。

*ちなみに、どうしてもjupyter上で完結したい場合は、処理時間は増えますが以下のコードで実行可能です。
(Google ColabでGPUを使うと15分くらいで実行できます。)

# 時間かけてでもjupyter上で実行したい場合は以下のコードで実行可能
(iter_num, input_size) = (51*3, 17)
mesh_list_iter = np.array(mesh_list).reshape(iter_num, input_size, feature_dim)

mesh_predict_tmp = []
for i, iter_list in enumerate(tqdm(mesh_list_iter)):
    tmp_list = svm.predict(iter_list)
    mesh_predict_tmp.append(tmp_list)
    print("epoch ", i, " has done ",list(tmp_list))

print(mesh_result)
mesh_predict_result = np.array(mesh_predict_tmp).reshape(-1)

実行後はアウトプットデータを読み込み、jupyter上で可視化できるようにしていきます。

mesh_predict_tmp = []
for i in range(input_file_num):
    with open("qsvm_terminal/output/output_br_cn_"+str(i)+".pkl", "rb") as f:
        mesh_predict_tmp.append(pickle.load(f))

mesh_predict_result = np.array(mesh_predict_tmp).reshape(-1)

これで境界面を可視化することができます。

plt.figure(figsize=(7, 7))               
heatmap(mesh_predict_result)
scatter_data(train_for_pred, test_for_pred, train_result, test_result ,yshift=-0.014)
plt.show()

# 赤がA、青がBのラベル。●が訓練データで、■がテストデータ

量子SVM用に用意されたデータセット(Ad Hoc Data)で機械学習

基本的な処理に乳癌データの場合と同じですが、データセットの用意の仕方が若干違います。
ここでは、重複する処理は割愛し、データセット用意と結果のみ載せます。
コード全体を見たい方はこちらのリンクからjupyter notebookを確認ください。

ad_hoc_dataのロードはbreast_cancer使用時とoption引数の指定が違います。

aqua_globals.random_seed = seed
sample_Total, training_input_unnormalized, test_input_unnormalized, class_labels = ad_hoc_data(
    training_size=training_dataset_size, 
    test_size=testing_dataset_size, 
    n=feature_dim, gap=0.3, plot_data=True
)

このようなラベルの分布を持つデータセットを使用します。

breast_cancer使用時は始めからx軸・y軸が-1から+1の区間に規格化されていましたが、ad_hoc_dataは自身で規格化する必要があります。

# Breast_cancerと同じスケール(-1, 1)でプロットするためにdatasetを規格化
(train_for_pred, _), _ = split_dataset_to_data_and_labels(training_input_unnormalized)
(test_for_pred, _), _ = split_dataset_to_data_and_labels(test_input_unnormalized)
dataset_array = np.vstack([train_for_pred, test_for_pred])
min_array, max_array = dataset_array.min(), dataset_array.max()
training_input_normalized = {k:(v-min_array)/(max_array-min_array)*2-1 for k,v in training_input_unnormalized.items()}
test_input_normalized = {k:(v-min_array)/(max_array-min_array)*2-1 for k,v in test_input_unnormalized.items()} 

　古典SVM

分類するコードは全く同じなので、モデル性能と境界面だけ表示します。

result: ['A', 'B', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'A', 'B', 'B', 'A', 'B', 'A', 'A']
truth : ['A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B']
{'accuracy': 0.7, 'precision': 0.64, 'recall': 0.9, 'specificity': 0.5, 'F1-measure': 0.75}

　量子SVM

breast_cancerデータの時と同様の処理でモデル性能を得ることができます。

result: ['A', 'B', 'A', 'B', 'A', 'A', 'A', 'A', 'B', 'B', 'A', 'B', 'B', 'A', 'B', 'A', 'A', 'B', 'A', 'B']
truth : ['A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B']

境界面の作成においては、jupyter上から.pyファイルを作成する際のファイル名が微妙に違うのでその部分だけコード記載します。

# 大量の予測データセットに対して実行する場合、terminal上で複数タブにて同時に実行する方が(若干)早い
# ただし、並列化するほど1つあたりの処理が大幅に遅くなるので、劇的な高速化はできない
(input_file_num, iter_num, input_size) = (3, 51, 17) 
mesh_list_tmp = np.array(mesh_list).reshape(input_file_num, iter_num, input_size, feature_dim)


code_string = '''
import os
import pickle
from tqdm import tqdm

file_dir = os.path.dirname(os.path.abspath(__file__))
with open(file_dir+"/../model/model_qsvm_ad_hc.pkl", "rb") as f:
    svm = pickle.load(f)

with open(file_dir+"/../input/input_INDEX.pkl", "rb") as f:
    input_array = pickle.load(f)

list_final = []
for epoch, epoch_array in enumerate(tqdm(input_array)):
    pred_tmp  = svm.predict(epoch_array)
    print("epoch ", epoch, " has done")
    list_final.append(list(pred_tmp))
print(list_final)

with open(file_dir+"/../output/output_ad_hc_INDEX.pkl", "wb") as f:
    pickle.dump(list_final, f)
'''  

with open("qsvm_terminal/model/model_qsvm_ad_hc.pkl", "wb") as f:
    pickle.dump(svm, f)

for i in range(input_file_num):
    # inputはbreast_cancerとad_hoc_data共通でOK
    with open("qsvm_terminal/input/input_{}.pkl".format(i),"wb") as f:
        pickle.dump(mesh_list_tmp[i], f)

    with open("qsvm_terminal/script/predict_ad_hc_{}.py".format(i),"w") as f:
        f.write(code_string.replace("INDEX",str(i)))

ターミナル上で複数タブを開き各タブでpython predict_br_cn_1.pyのように最後の添字を変えて実行してください。
その後は同様の処理(読み出しファイルが違うことに注意)で境界面を表示できます。

# データ読み出し
mesh_predict_tmp = []
for i in range(input_file_num):
    with open("qsvm_terminal/output/output_ad_hc_"+str(i)+".pkl", "rb") as f:
        mesh_predict_tmp.append(pickle.load(f))

mesh_predict_result = np.array(mesh_predict_tmp).reshape(-1)

# 境界面可視化
plt.figure(figsize=(7, 7))               
heatmap(mesh_predict_result)
scatter_data(train_for_pred, test_for_pred, train_result, test_result ,yshift=-0.014)
plt.show()
# 赤がA、青がBのラベル。●が訓練データで、■がテストデータ

量子回路の深さ(量子SVMの概要の節を参照)をデフォルトのdepth=2ままでは、現状では両方のデータセットで古典SVMがモデル性能と計算時間において優っているという結果になりました。

次に、Ad Hoc Dataで回路の深さを変えた場合のモデル性能と計算時間の変化をみていきます。

回路の深さとモデル性能・計算時間の関係

今回は回路の深さを1から10まで変えて計算しました。
Jupyter上で計算するとメモリの問題かBrokenProcessPoolエラーが出てしまいました。
なのでdepth_variator.pyとして書き出し、Jupyter上でsubprocessコマンドで実行しています。

# jupyter上で回すとBrokenProcessPoolエラーが発生することがあるため、Terminalで実行。

# shotsやseedをjupyter上で変更する場合は注意
code_string_depth = '''
import os
import time
import pickle
import numpy as np
from qiskit import BasicAer
from qiskit.aqua.algorithms import QSVM
from qiskit.aqua import QuantumInstance
from qiskit.aqua import aqua_globals
from qiskit.circuit.library import ZZFeatureMap
from qiskit.ml.datasets import ad_hoc_data
from qiskit.aqua.utils import split_dataset_to_data_and_labels

# サンプリング設定
feature_dim = 2  # 特徴量の数
training_dataset_size = 20
testing_dataset_size = 10
shots = 1024
seed = 10598


# モデル評価
def eval_model(pred_label, print_mode=True):
    test_label = ["A"]*testing_dataset_size + ["B"]*testing_dataset_size

    accuracy = sum([x == y for x, y in zip(pred_label, test_label)])/len(test_label)
    precision = \
        sum([x == y for x, y in zip(pred_label, test_label) if x == "A"])/sum([x == "A" for x in pred_label])
    recall = \
        sum([x == y for x, y in zip(pred_label, test_label) if y == "A"])/sum([y == "A" for y in test_label])
    specificity = \
        sum([x == y for x, y in zip(pred_label, test_label) if y == "B"])/sum([y == "B" for y in test_label])
    f1 = 2*recall*precision/(recall + precision)
    eval_dict = {"accuracy": accuracy, "precision": precision, "recall": recall, "specificity": specificity, "F1-measure":f1}
    if print_mode:
        print("result: ", pred_label)
        print("truth : ", test_label)
        print(eval_dict)
    else:
        return eval_dict


# サンプル取得(ad_hoc_dataを使用)
aqua_globals. random_seed = seed
sample_Total, training_input_unnormalized, test_input_unnormalized, class_labels = ad_hoc_data(
    training_size=training_dataset_size, 
    test_size=testing_dataset_size, 
    n=feature_dim, gap=0.3, plot_data=False
)

# 規格化
(train_for_pred, _), _ = split_dataset_to_data_and_labels(training_input_unnormalized)
(test_for_pred, _), _ = split_dataset_to_data_and_labels(test_input_unnormalized)
dataset_array = np.vstack([train_for_pred, test_for_pred])
min_array, max_array = dataset_array.min(), dataset_array.max()
training_input_normalized = {k: (v-min_array)/(max_array-min_array)*2-1 for k, v in training_input_unnormalized.items()}
test_input_normalized = {k: (v-min_array)/(max_array-min_array)*2-1 for k, v in test_input_unnormalized.items()}

# 計算実行
dict_result = {}
for depth in range(1, 11):
    start_time = time.time()
    backend = BasicAer.get_backend('qasm_simulator')
    feature_map_depth = ZZFeatureMap(feature_dim, reps=depth)
    svm_depth = QSVM(feature_map_depth, training_input_normalized, test_input_normalized, test_for_pred)
    svm_depth.random_seed = seed
    quantum_instance = QuantumInstance(backend, shots=shots, seed_simulator=seed, seed_transpiler=seed)
    result_depth = svm_depth.run(quantum_instance)
    eval_metrics_dict = eval_model(result_depth["predicted_classes"], print_mode=False)
    print("depth: ", depth)
    eval_metrics_dict_print = {k: round(v, 2) for k, v in eval_metrics_dict.items()}
    print(eval_metrics_dict_print)
    print("--- %s seconds ---" % (round(time.time() - start_time)))
    dict_result[depth] = {"evaluation": eval_metrics_dict,
                          "time": time.time() - start_time,
                          "result": result_depth}

dir = os.path.dirname(os.path.abspath(__file__))
with open(dir+"/../output/depth_variator_output.pkl", "wb") as f:
    pickle.dump(dict_result, f)

'''

with open("qsvm_terminal/script/depth_variator.py", "w") as f:
        f.write(code_string_depth)
'''

with open("qsvm_terminal/script/depth_variator.py", "w") as f:
        f.write(code_string_depth)

実行はsubprocessコマンドでリアルタイムに取得しています。

cmd = ["python", "qsvm_terminal/script/depth_variator.py"]
proc = subprocess.Popen(cmd, stdout = subprocess.PIPE, stderr = subprocess.STDOUT)
for line in iter(proc.stdout.readline,b''):
    print(line.rstrip().decode("utf8"))

depth: 1
{'accuracy': 0.6, 'precision': 0.62, 'recall': 0.5, 'specificity': 0.7, 'F1-measure': 0.56}
--- 248 seconds ---
depth: 2
{'accuracy': 0.5, 'precision': 0.5, 'recall': 0.5, 'specificity': 0.5, 'F1-measure': 0.5}
--- 238 seconds ---
depth: 3
{'accuracy': 0.5, 'precision': 0.5, 'recall': 0.2, 'specificity': 0.8, 'F1-measure': 0.29}
--- 239 seconds ---
depth: 4
{'accuracy': 0.7, 'precision': 0.67, 'recall': 0.8, 'specificity': 0.6, 'F1-measure': 0.73}
--- 147 seconds ---
depth: 5
{'accuracy': 0.55, 'precision': 0.56, 'recall': 0.5, 'specificity': 0.6, 'F1-measure': 0.53}
--- 247 seconds ---
depth: 6
{'accuracy': 0.45, 'precision': 0.43, 'recall': 0.3, 'specificity': 0.6, 'F1-measure': 0.35}
--- 260 seconds ---
depth: 7
{'accuracy': 0.6, 'precision': 0.62, 'recall': 0.5, 'specificity': 0.7, 'F1-measure': 0.56}
--- 292 seconds ---
depth: 8
{'accuracy': 0.45, 'precision': 0.45, 'recall': 0.5, 'specificity': 0.4, 'F1-measure': 0.48}
--- 274 seconds ---
depth: 9
{'accuracy': 0.5, 'precision': 0.5, 'recall': 0.7, 'specificity': 0.3, 'F1-measure': 0.58}
--- 283 seconds ---
depth: 10
{'accuracy': 0.55, 'precision': 0.54, 'recall': 0.7, 'specificity': 0.4, 'F1-measure': 0.61}
--- 286 seconds ---

得られた結果をdataFrameにまとめておきます。

with open("qsvm_terminal/output/depth_variator_output.pkl", "rb") as f:
    result_dict = pickle.load(f)

records = [{**value["evaluation"], **{"time": value["time"]}} for value in result_dict.values()]
df = pd.DataFrame(records,index=range(1,11))
df

これで回路の深さごとにモデル性能を可視化することができます。

plt.figure(figsize=(10,6))
for col in df.columns:
    if col == "time": continue
    plt.plot(df[col],linestyle='-', marker='o',label=col)

plt.axis('auto')
plt.ylim(0,1)    
plt.xticks(range(1,11))
plt.rcParams["legend.edgecolor"] = 'black'
plt.ylabel("Model Performance Metrics",fontsize=18)
plt.xlabel("Circuit Depth",fontsize=18)
plt.legend(bbox_to_anchor=(1, 0), loc='lower right', borderaxespad=0, fontsize=11)
plt.show()

計算時間も同様に可視化しました。

plt.figure(figsize=(10,6))
plt.plot(df["time"],linestyle='-', marker='o')
plt.xticks(range(1,11))
plt.ylabel("Calclulation Time",fontsize=18)
plt.xlabel("Circuit Depth",fontsize=18)
plt.ylim(0,400)
plt.show()

モデル性能が比較的良くなっている回路の深さ(depth=4)では計算時間が抑えられるようです。
何度か実行して計算時間を計測しなおしてみても同様の傾向がみられました。

まとめ

この記事では量子SVMの解説と実行を行いました。
量子SVMはカーネル$K(\boldsymbol{x}_i,\boldsymbol{x}_k)=\left|\left.\langle \phi(\boldsymbol{x}_i)\right|\phi(\boldsymbol{x}_k)\rangle\right|^2$を利用するのがポイントです。
実装では量子シミュレータでモデル性能と計算時間を評価し、SVMの作る境界面を可視化しました。
量子SVMのハイパーパラメータで回路の深さをチューニングすることで、古典SVMと同等の精度が出ることを確認できました。

追記

この記事ではモデルの汎化性能は求めていません。
cross validationをするには自作で関数を定義する必要があり、単なるモデル作成にもかなり時間がかかることから省略しています。

修正にあたり、@med_quantum_menさんにアドバイスを頂きました。