demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Qboost: Binary Classification with Quantum Computer\n",
    "\n",
    "The D-Wave quantum computer has been widely studied as a discrete optimization engine that accepts any problem formulated as quadratic unconstrained  binary  optimization  (QUBO). In 2008, Google and D-Wave published a paper, [Training a Binary Classifier with the Quantum Adiabatic Algorithm](https://arxiv.org/pdf/0811.0416.pdf), which describes how the `Qboost` ensemble method makes binary classification amenable to quantum computing: the problem is formulated as a thresholded linear superposition of a set of weak classifiers and the D-Wave quantum computer is  used to optimize the weights in a learning process that strives to minimize the training error and number of weak classifiers.\n",
    "\n",
    "This notebook demonstrates and explains how the Qboost algorithm can be used to solve a binary classification problem. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Few Words on Ensemble Methods\n",
    "\n",
    "Ensemble methods build a strong classifier (an improved model) by combining weak classifiers with the goal of:\n",
    "\n",
    "* decreasing variance (bagging)\n",
    "* decreasing bias (boosting)\n",
    "* improving prediction (voting)\n",
    "\n",
    "![Boosting Algorithm](images/boosting.jpg)\n",
    "\n",
    "### Bagging, Boosting, and Voting\n",
    "\n",
    "The ensemble method produces new training data sets by random sampling with replacement from the original set. In _bagging_, any element has the same probability to appear in a new dataset; in _boosting_, data elements are weighted before they are collected in the new dataset. Another distinction is that bagging is parallelizable but boosting has to be executed sequentially. You can learn more about the differences between these methods here: https://quantdare.com/what-is-the-difference-between-bagging-and-boosting/.\n",
    "\n",
    "Voting operates on labels only. Unlike boosting, the aggeragated classification performance is not used to further polish each weak classifier. Voting has two typical requirements of its collection of  weak classifiers: that there be __many__ and that they be __diverse__.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Weak and Strong Classifiers\n",
    "For this reference example we chose the following four classifiers:\n",
    "    1. Adaboost\n",
    "    2. Decision Trees\n",
    "    3. Random Forest\n",
    "    4. Qboost\n",
    "Note that you can replace any of these with any commonly used classification model. Also, an ensemble method can use a strong classifier instead of a weak one, and in this example we embed the Qboost classifier itself, __QboostPlus__ in the following code, with the first three.  \n",
    "\n",
    "### Adaboost\n",
    "Adaboost combines a number of $N$ weak classifiers into a strong one as\n",
    "$$C(x) = sign\\left(\\sum_i^N w_i c_i(x)\\right),$$\n",
    "with $c_i(x) \\in [-1, +1]$ being the $i$-th weak classifier:\n",
    "\n",
    "$$c_i(x) = sign(w'*x + b)$$\n",
    "\n",
    "The loss function of Adaboost is defined as\n",
    "$$\n",
    "L = \\sum_{n=1}^N \\exp\\left\\{ - y_n \\sum_{s=1}^S w_sc_k(x_n)\\right\\}.\n",
    "$$\n",
    "\n",
    "The strong classifier $C(\\cdot)$ is constructed in an iterative fashion. In each iteration, one weak classifier\n",
    "is selected and re-learned to minimize the weighted error function. Its weight is adjusted and renormalized to make sure the sum of all weights equals 1. \n",
    "\n",
    "The final classification model will be decided by a weighted “vote” of all the weak classifiers. \n",
    "\n",
    "### Decision Trees\n",
    "A decision tree builds on a tree structure with non-leaf nodes encoding decision rules and leaf nodes encoding labels. You construct a decision tree by optimizing either entropic or information-theoretic metrics. Controlling the depth of a decision tree indirectly decides the sub-dimension of the dataset. \n",
    "\n",
    "Decision trees are often chosen as the weak classifiers in Adaboost because they are both simple to construct and fast to do inference. The `scikit-learn` package implements its Adaboost method with decision trees of depth 1, also known as _tree stumps_. This reference examples demonstrates an alternative implementation of boosting with a number of deeper decision trees.\n",
    "### Random Forest \n",
    "Random forest is an ensemble method that typically implements bagging on a set of decision trees. By introducing randomness in the selection of an optimized feature in the training of the underlying decision trees, the ensemble diversifies the weightings of its collection of weak classifiers, generally resulting in an improved model.    \n",
    "\n",
    "### Qboost\n",
    "To make use of the optimization power of D-Wave quantum annealer, we needs to formulate our objective function as a quadratic unconstrained binary optimization (QUBO) problem. Therefore, we replace the exponential loss as in Adaboost with the following quadratic loss\n",
    "$$\n",
    "w* = \\arg\\min_w\\left(\\sum_s \\left(\\frac{1}{N}\\sum_n^N w_nc_n(x_s) - y_s\\right)^2\\right) + \\lambda ||w||_0,\n",
    "$$\n",
    "where the regularization term is added to enable controlling of weight sparsity.\n",
    "\n",
    "Note in Qboost, the weight vector is binary.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training and Comparing Performance\n",
    "Now we define functions used in the following experiemnts to train our selected classifiers and provide metrics for comparing performance.\n",
    "First, let us import the required packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import necessary packages\n",
    "from sklearn import preprocessing, metrics\n",
    "from sklearn.ensemble import AdaBoostClassifier, RandomForestClassifier\n",
    "from sklearn.datasets.mldata import fetch_mldata\n",
    "from sklearn.datasets import load_breast_cancer\n",
    "from dwave.system.samplers import DWaveSampler\n",
    "from dwave.system.composites import EmbeddingComposite\n",
    "\n",
    "from qboost import WeakClassifiers, QBoostClassifier, QboostPlus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let us define the `metric` and `train_model` functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the functions required in this example\n",
    "def metric(y, y_pred):\n",
    "    \"\"\"\n",
    "    :param y: true label\n",
    "    :param y_pred: predicted label\n",
    "    :return: metric score\n",
    "    \"\"\"\n",
    "\n",
    "    return metrics.accuracy_score(y, y_pred)\n",
    "\n",
    "\n",
    "def train_model(X_train, y_train, X_test, y_test, lmd):\n",
    "    \"\"\"\n",
    "    :param X_train: training data\n",
    "    :param y_train: training label\n",
    "    :param X_test: testing data\n",
    "    :param y_test: testing label\n",
    "    :param lmd: lambda used in regularization\n",
    "    :return:\n",
    "    \"\"\"\n",
    "\n",
    "    # define parameters used in this function\n",
    "    NUM_READS = 1000\n",
    "    NUM_WEAK_CLASSIFIERS = 30\n",
    "    TREE_DEPTH = 2\n",
    "    DW_PARAMS = {'num_reads': NUM_READS,\n",
    "                 'auto_scale': True,\n",
    "                 'num_spin_reversal_transforms': 10,\n",
    "                 'postprocess': 'optimization',\n",
    "                 }\n",
    "\n",
    "    # define sampler\n",
    "    dwave_sampler = DWaveSampler()\n",
    "    emb_sampler = EmbeddingComposite(dwave_sampler)\n",
    "\n",
    "    N_train = len(X_train)\n",
    "    N_test = len(X_test)\n",
    "    print(\"\\n======================================\")\n",
    "    print(\"Train size: %d, Test size: %d\" %(N_train, N_test))\n",
    "    print('Num weak classifiers:', NUM_WEAK_CLASSIFIERS)\n",
    "\n",
    "    # Preprocessing data\n",
    "    imputer = preprocessing.Imputer()\n",
    "    scaler = preprocessing.StandardScaler()\n",
    "    normalizer = preprocessing.Normalizer()\n",
    "\n",
    "    X_train = scaler.fit_transform(X_train)\n",
    "    X_train = normalizer.fit_transform(X_train)\n",
    "\n",
    "    X_test = scaler.fit_transform(X_test)\n",
    "    X_test = normalizer.fit_transform(X_test)\n",
    "\n",
    "    ## Adaboost\n",
    "    print('\\nAdaboost')\n",
    "    clf1 = AdaBoostClassifier(n_estimators=NUM_WEAK_CLASSIFIERS)\n",
    "    clf1.fit(X_train, y_train)\n",
    "    y_train1 = clf1.predict(X_train)\n",
    "    y_test1 = clf1.predict(X_test)\n",
    "#     print(clf1.estimator_weights_)\n",
    "    print('accu (train): %5.2f'%(metric(y_train, y_train1)))\n",
    "    print('accu (test): %5.2f'%(metric(y_test, y_test1)))\n",
    "\n",
    "    # Ensembles of Decision Tree\n",
    "    print('\\nDecision tree')\n",
    "    clf2 = WeakClassifiers(n_estimators=NUM_WEAK_CLASSIFIERS, max_depth=TREE_DEPTH)\n",
    "    clf2.fit(X_train, y_train)\n",
    "    y_train2 = clf2.predict(X_train)\n",
    "    y_test2 = clf2.predict(X_test)\n",
    "#     print(clf2.estimator_weights)\n",
    "    print('accu (train): %5.2f' % (metric(y_train, y_train2)))\n",
    "    print('accu (test): %5.2f' % (metric(y_test, y_test2)))\n",
    "    \n",
    "    # Random forest\n",
    "    print('\\nRandom Forest')\n",
    "    clf3 = RandomForestClassifier(max_depth=TREE_DEPTH, n_estimators=NUM_WEAK_CLASSIFIERS)\n",
    "    clf3.fit(X_train, y_train)\n",
    "    y_train3 = clf3.predict(X_train)\n",
    "    y_test3 = clf3.predict(X_test)\n",
    "    print('accu (train): %5.2f' % (metric(y_train, y_train3)))\n",
    "    print('accu (test): %5.2f' % (metric(y_test, y_test3)))\n",
    "\n",
    "    # Qboost\n",
    "    print('\\nQBoost')\n",
    "    clf4 = QBoostClassifier(n_estimators=NUM_WEAK_CLASSIFIERS, max_depth=TREE_DEPTH)\n",
    "    clf4.fit(X_train, y_train, emb_sampler, lmd=lmd, **DW_PARAMS)\n",
    "    y_train4 = clf4.predict(X_train)\n",
    "    y_test4 = clf4.predict(X_test)\n",
    "    print(clf4.estimator_weights)\n",
    "    print('accu (train): %5.2f' % (metric(y_train, y_train4)))\n",
    "    print('accu (test): %5.2f' % (metric(y_test, y_test4)))\n",
    "\n",
    "    # QboostPlus\n",
    "    print('\\nQBoostPlus')\n",
    "    clf5 = QboostPlus([clf1, clf2, clf3, clf4])\n",
    "    clf5.fit(X_train, y_train, emb_sampler, lmd=lmd, **DW_PARAMS)\n",
    "    y_train5 = clf5.predict(X_train)\n",
    "    y_test5 = clf5.predict(X_test)\n",
    "    print(clf5.estimator_weights)\n",
    "    print('accu (train): %5.2f' % (metric(y_train, y_train5)))\n",
    "    print('accu (test): %5.2f' % (metric(y_test, y_test5)))\n",
    "\n",
    "    print(\"===========================================================================\")\n",
    "    print(\"Method \\t Adaboost \\t DecisionTree \\t RandomForest \\t Qboost \\t Qboost+\")\n",
    "    print(\"Train\\t %5.2f \\t\\t %5.2f \\t\\t %5.2f \\t\\t %5.2f \\t\\t %5.2f\"% (metric(y_train, y_train1),\n",
    "                                                                         metric(y_train, y_train2),\n",
    "                                                                         metric(y_train, y_train3),\n",
    "                                                                         metric(y_train, y_train4),\n",
    "                                                                         metric(y_train, y_train5),\n",
    "                                                                        ))\n",
    "    print(\"Test\\t %5.2f \\t\\t %5.2f \\t\\t %5.2f \\t\\t %5.2f \\t\\t %5.2f\"% (metric(y_test, y_test1),\n",
    "                                                                       metric(y_test, y_test2),\n",
    "                                                                       metric(y_test, y_test3),\n",
    "                                                                       metric(y_test, y_test4),\n",
    "                                                                       metric(y_test, y_test5)))\n",
    "    print(\"===========================================================================\")\n",
    "    \n",
    "    return [clf1, clf2, clf3, clf4, clf5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Experiments\n",
    "Now we're ready to run some experiments.\n",
    "First, import the required packages. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "#import os\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Experiment 1: Binary Classfication on the MNIST Dataset \n",
    "This example transforms the MNIST dataset (handwritten digits) into a binary classification problem. We assume all digits that are smaller than 5 are labelled as -1 and the rest digits are labelled as +1.\n",
    "\n",
    "First, let us load the MINIST dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mnist = fetch_mldata('MNIST original', data_home='data')\n",
    "\n",
    "idx_01 = np.where(mnist.target <= 9)[0]\n",
    "np.random.shuffle(idx_01)\n",
    "idx_01 = idx_01[:15000]\n",
    "idx_train = idx_01[:2*len(idx_01)//3]\n",
    "idx_test = idx_01[2*len(idx_01)//3:]\n",
    "\n",
    "X_train = mnist.data[idx_train]\n",
    "X_test = mnist.data[idx_test]\n",
    "\n",
    "y_train = 2*(mnist.target[idx_train] >4) - 1\n",
    "y_test = 2*(mnist.target[idx_test] >4) - 1\n",
    "\n",
    "print(\"Training data size: (%d, %d)\" %(X_train.shape))\n",
    "print(\"Testing data size: (%d, %d)\" %(X_test.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us visualize the digits: digits with class $+1$ are shown as images with a black background while digits with class $-1$  as images with a white background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(16):\n",
    "    if y_train[i] == 1:\n",
    "        COLORMAP = 'gray'\n",
    "    else:\n",
    "        COLORMAP = 'gray_r'\n",
    "    plt.subplot(4,4, i+1)\n",
    "    plt.imshow(X_train[i].reshape(28,28), cmap=COLORMAP)\n",
    "    plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now train the model and compare the results of the selected classifiers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# start training the model\n",
    "clfs = train_model(X_train, y_train, X_test, y_test, 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: for this cell, Graphviz executables must be on systems PATH \n",
    "# You can optionally visualize the decision trees by uncommenting the following code\n",
    "# import graphviz\n",
    "# from sklearn import tree\n",
    "# clf = clfs[0]\n",
    "# graph = graphviz.Source(tree.export_graphviz(clf.estimators_[0], out_file=None))\n",
    "# graph.render(None, view=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Experiment 2: Wisconsin Breast Cancer\n",
    "\n",
    "This example classifies tumors in scikit-learn's Wisconsis breast cancer dataset as either malignant or benign (binary classification).\n",
    "\n",
    "First, let us load the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wisc = load_breast_cancer()\n",
    "\n",
    "idx = np.arange(len(wisc.target))\n",
    "np.random.shuffle(idx)\n",
    "\n",
    "# train on a random 2/3 and test on the remaining 1/3\n",
    "idx_train = idx[:2*len(idx)//3]\n",
    "idx_test = idx[2*len(idx)//3:]\n",
    "\n",
    "X_train = wisc.data[idx_train]\n",
    "X_test = wisc.data[idx_test]\n",
    "\n",
    "y_train = 2 * wisc.target[idx_train] - 1  # binary -> spin\n",
    "y_test = 2 * wisc.target[idx_test] - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now train the model and compare the results of the selected classifiers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# train the model\n",
    "clfs = train_model(X_train, y_train, X_test, y_test, 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}