docs/tutorial/ja/machine_learning/qboost.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# アニーリングを用いたアンサンブル学習（QBoost)\n",
"\n",
"[QBoost](https://arxiv.org/abs/0811.0416)は量子アニーリングを用いたアンサンブル学習の一つです。\n",
"アンサンブル学習は弱い予測器を多数用意し、その予測器の各予測結果を組み合わせることで、最終的に精度の良い予測結果を得る手法です。\n",
"\n",
"## 概要と原理\n",
"\n",
"このQBoostは、入力信号$\\bm{x}$がどのような性質を持っているのかを精度良く識別することを目標としたアルゴリズムです。ここでは2つの値$\\pm 1$のどちらを入力信号に割り当てるべきかという問題を考えましょう。例として、$\\bm{x}$が画像データを表しており、その画像に写っているものが犬か猫かを識別するといったタスクを想像すると良いでしょう。アンサンブル学習では、複数の予測器を利用することで、より良い予測精度を達成すること(ブースティング)を目指します。ここでは、あまり性能の良くない予測器(弱い予測器)をたくさん用意します。性能が良くないという意味は、入力に対して正しい出力をしないことが多いことを意味します。これらの予測器の出力を$c_i (\\bm{x}) \\in \\{ -1, 1\\} \\ (i=0, 1, \\dots, N-1)$とします。いくつかの弱い予測器の出力の和を取ることで、より良い予測ができるというのが基本的な考え方です。これを数式で表すと\n",
"\n",
"$$\n",
"C(\\bm{x}) \n",
"= \\mathrm{sgn} \\left( \\sum_{i=0}^{N-1} w_i c_i (\\bm{x}) \\right) \\tag{1}\n",
"$$\n",
"\n",
"となります。ここで$w_i \\in \\{0, 1\\}$で、$i$番目の予測器を使うか使わないかを表します。どの予測器を用いると、できるだけ少ない数の弱い予測器でより良い性能が得られるかを明らかにしましょう。 \n",
"このために、教師あり学習を用いて最適な$\\{w_i\\}$の組を求めることにします。教師データを$(\\bm{x}^{(d)}, y^{(d)}) \\ (d= 0, 1, \\dots, D-1)$を多数用意します($D \\gg 1$)。それらをできるだけ忠実に再現するように$\\{w_i\\}$を調整します。 \n",
"この方針をより具体的に表すと、次のハミルトニアンを$\\{w_i\\}$について最小化することを目指せば良いとわかります。\n",
"\n",
"$$\n",
"H(\\bm{w}) \n",
"= \\sum_{d=0}^{D-1} \\left( \\frac{1}{N} \\sum_{i=0}^{N-1} w_i c_i (\\bm{x}^{(d)}) - y^{(d)}\\right)^2 + \\lambda \\sum_{i=0}^{N-1} w_i \\tag{2}\n",
"$$ \n",
"\n",
"このハミルトニアンの最小化を通して、教師データ$y^{(d)}$との差ができるだけ小さくなるようにします。式(1)の右辺をそのまま使うと、符号関数があるために$w_i$の2次形式にならず、イジング模型に帰着することができません。そのため、符号関数の引数$\\sum_i w_i c_i$の$1/N$倍と教師データ$y^{(d)}$との差の2乗を最小化する問題にしています。$1/N$の係数は、$\\sum_i w_i c_i(\\bm{x})$の最大値が$N$であるために$y^{(d)}= \\pm 1$との差が大きくなりすぎないのように調整するためのものです。$\\lambda (>0)$がかかった項は、あまり多くの$w_i$を1にせず、比較的少数の弱い予測器で効率良く構成するための項(正則化項)を表します。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## JijModelingによるモデル構築\n",
"\n",
"### 変数の定義\n",
"\n",
"式(2)で用いられている変数を、以下のようにして定義しましょう。"
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {},
"outputs": [],
"source": [
"import jijmodeling as jm\n",
"\n",
"# define variables\n",
"c = jm.Placeholder(\"c\", ndim=2)\n",
"y = jm.Placeholder(\"y\", ndim=1)\n",
"N = c.len_at(0, latex=\"N\")\n",
"D = c.len_at(1, latex=\"D\")\n",
"w = jm.BinaryVar(\"w\", shape=(N, ))\n",
"lamb = jm.Placeholder(\"lamb\", latex=\"\\lambda\")\n",
"i = jm.Element(\"i\", belong_to=(0, N))\n",
"d = jm.Element(\"d\", belong_to=(0, D))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`c = jm.Placeholder('c', ndim=2)`で式(2)の$c$を定義しています。\n",
"そのリストの要素数から、弱い予測器の数$N$と教師データ数$D$をそれぞれ`N, D`として定義しています。\n",
"それらを用いて、最適化に用いるバイナリ変数`w`と教師データのバイナリ値`y`を定義しています。\n",
"式(2)の$\\lambda$を`lamb`として定義し、最後に式(2)で用いられている添字を`i, d`のように表しています。\n",
"\n",
"### 目的関数の実装\n",
"\n",
"式(2)を実装しましょう。"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {},
"outputs": [],
"source": [
"# set problem\n",
"problem = jm.Problem(\"QBoost\")\n",
"# set objective function 1: minimize the sum of differences\n",
"obj1 = jm.sum(d, (jm.sum(i, w[i]*c[i, d])/N-y[d])**2)\n",
"# set objective function 2: minimize regularization term\n",
"obj2 = lamb * w[:].sum()\n",
"problem += obj1 + obj2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`w[:].sum()`とすることで、$\\sum_i w_i$を簡潔に実装することができます。 \n",
"ここまでの実装が正しくされているかを確認しましょう。"
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$$\\begin{array}{cccc}\\text{Problem:} & \\text{QBoost} & & \\\\& & \\min \\quad \\displaystyle \\sum_{d = 0}^{D - 1} \\left(\\left(\\sum_{i = 0}^{N - 1} w_{i} \\cdot c_{i, d} \\cdot N^{(-1)} - y_{d}\\right)^{2}\\right) + \\lambda \\cdot \\sum_{\\ast_{0} = 0}^{N - 1} w_{\\ast_{0}} & \\\\\\text{{where}} & & & \\\\& w & 1\\text{-dim binary variable}\\\\\\end{array}$$"
],
"text/plain": [
"<jijmodeling.Problem at 0x36a7680>"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### インスタンスの作成\n",
"\n",
"実際に実行するタスクなどを設定しましょう。今回は弱い予測器を[scikit-learn](https://scikit-learn.org/stable/)のdecision stump(決定株: 一層の決定木)を用います。\n",
"また用いるデータはscikit-learnの癌識別データセットを使用します。"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from sklearn import datasets\n",
"from sklearn.tree import DecisionTreeClassifier as DTC\n",
"\n",
"def prediction_from_train(N, X_train, y_train, X_test):\n",
" # set the number of ensembles to be taken out for one sample in bootstrap sampling\n",
" sample_train = 40\n",
" # set model\n",
" models = [DTC(splitter=\"random\", max_depth=1) for i in range(N)]\n",
" for model in models:\n",
" # extract randomly\n",
" train_idx = np.random.choice(np.arange(X_train.shape[0]), sample_train)\n",
" # make decision tree with variables\n",
" model.fit(X=X_train[train_idx], y=y_train[train_idx])\n",
" y_pred_list_train = []\n",
" for model in models:\n",
" # execute prediction with model\n",
" y_pred_list_train.append(model.predict(X_train))\n",
" y_pred_list_train = np.asanyarray(y_pred_list_train)\n",
" y_pred_list_test = []\n",
" for model in models:\n",
" # execute with test data\n",
" y_pred_list_test.append(model.predict(X_test))\n",
" y_pred_list_test = np.array(y_pred_list_test)\n",
" return y_pred_list_train, y_pred_list_test\n",
"\n",
"# load data\n",
"cancer_data = datasets.load_breast_cancer()\n",
"# set the number of train data\n",
"num_train = 200\n",
"# add noise to feature\n",
"noisy_data = np.concatenate((cancer_data.data, np.random.rand(cancer_data.data.shape[0], 30)), axis=1)\n",
"# convert (0, 1) label to (-1, 1)\n",
"labels = (cancer_data.target-0.5) * 2\n",
"# divide dataset to train and test\n",
"X_train = noisy_data[:num_train, :]\n",
"X_test = noisy_data[num_train:, :]\n",
"y_train = labels[:num_train]\n",
"y_test = labels[num_train:]\n",
"# set the number of classifer\n",
"inst_N = 20\n",
"# predict from train data using dicision tree classifier\n",
"y_pred_list_tarin, y_pred_list_test = prediction_from_train(inst_N, X_train, y_train, X_test)\n",
"# set lambda (coefficient of 2nd term)\n",
"inst_lamb = 10.0\n",
"instance_data = {'y': y_train, 'c': y_pred_list_tarin, 'lamb': inst_lamb, 'y_train': y_train, 'y_test': y_test, 'y_pred_list_test': y_pred_list_test}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"デモンストレーションのために、ノイズとなる特徴量を加えたものを実際のデータとして用います。\n",
"`prediction_from_train`関数を用いて弱い予測器の作成と、それらの予測器からの出力$c_i (\\bm{x}^{(d)})$を作成しています。\n",
"ここでは弱い予測器の数を20、教師データ数は200です。最後に式(2)の$\\lambda$の値を3.0としています。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### JijModeling transpilerによるPyQUBOへの変換\n",
"\n",
"ここまで行われてきた実装は、全てJijModelingによるものでした。\n",
"これを[PyQUBO](https://pyqubo.readthedocs.io/en/latest/)に変換することで、OpenJijはもちろん、他のソルバーを用いた組合せ最適化計算を行うことが可能になります。"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {},
"outputs": [],
"source": [
"import jijmodeling_transpiler as jmt\n",
"\n",
"# compile\n",
"compiled_model = jmt.core.compile_model(problem, instance_data, {})\n",
"# get qubo model\n",
"pubo_builder = jmt.core.pubo.transpile_to_pubo(compiled_model=compiled_model, relax_method=jmt.core.pubo.RelaxationMethod.AugmentedLagrangian)\n",
"qubo, const = pubo_builder.get_qubo_dict(multipliers={})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": []
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {},
"outputs": [],
"source": [
"import openjij as oj\n",
"\n",
"# set sampler\n",
"sampler = oj.SASampler()\n",
"# solve problem\n",
"result = sampler.sample_qubo(qubo, num_reads=100)"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [],
"source": [
"# decode a result to JijModeling sampleset\n",
"sampleset = jmt.core.pubo.decode_from_openjij(result, pubo_builder, compiled_model)\n",
"objectives = np.array(sampleset.evaluation.objective)\n",
"lowest_index = np.argmin(objectives)\n",
"w_indices = sampleset.record.solution[\"w\"][lowest_index][0][0]"
]
},
{
"cell_type": "code",
"execution_count": 74,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Accuracy of QBoost is 0.8888888888888888\n"
]
}
],
"source": [
"from sklearn.metrics import accuracy_score\n",
"\n",
"y_pred_list_train = instance_data['c']\n",
"y_train = instance_data['y_train']\n",
"y_test = instance_data['y_test']\n",
"y_pred_list_test = instance_data['y_pred_list_test']\n",
"accs_train_oj = []\n",
"accs_test_oj = []\n",
"for solution in sampleset.record.solution[\"w\"]:\n",
" idx_clf_oj = solution[0][0]\n",
" y_pred_train_oj = np.sign(np.sum(y_pred_list_train[idx_clf_oj, :], axis=0))\n",
" y_pred_test_oj = np.sign(np.sum(y_pred_list_test[idx_clf_oj, :], axis=0))\n",
" acc_train_oj = accuracy_score(y_true=y_train, y_pred=y_pred_train_oj)\n",
" acc_test_oj = accuracy_score(y_true=y_test, y_pred=y_pred_test_oj)\n",
" accs_train_oj.append(acc_train_oj)\n",
" accs_test_oj.append(acc_test_oj) \n",
"print('Accuracy of QBoost is {}'.format(max(accs_test_oj)))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.9.5 ('.venv': venv)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.5"
},
"vscode": {
"interpreter": {
"hash": "8e4ef42922154d6c53cb4e5074b09fe7c3d2526e8e95a3a7b0551dc780a36e7a"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}