ML 3.2 - Support Vector Machine - from scratch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Internet Resources:\n",
    "\n",
    "[Python Programming.net - machine learning episodes 20-33](https://pythonprogramming.net/support-vector-machine-intro-machine-learning-tutorial/)  \n",
    "[MIT OpenCourseWare 16 Machine Learning: Support Vector Machines](https://www.youtube.com/watch?v=_PwhiWxHK8o)  \n",
    "[Victor Lavrenko - IR20.6 Sequential minimal optimization (SMO)](https://www.youtube.com/watch?v=I73oALP7iWA)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAUfUlEQVR4nO3db2xb1eHG8cfYCWkgJLHdJktaB1LSTSDRKWpV1K0wVq8vGJMAoSqtACE0ITVjSGPDUG1jTF01z6xkKikqL6rAuheBSavQpv1BHhqV1kmUuN2yjvWfioMUksy5TQm0aWPn/l6U+hfLbpPWN74+yffzCh/f3jw9Lo9ujq99PLZt2wIAGOs6twMAAIpDkQOA4ShyADAcRQ4AhqPIAcBwFDkAGM7n1g8eHBzMGwsGg0qlUi6kuXZknnum5ZXIXCoLLXNTU1PBca7IAcBwFDkAGI4iBwDDUeQAYDiKHAAM58hdK3/4wx/0zjvvyOPxaNmyZers7FRlZaUTpwbmHe/AgGpiMXmHhpRpbNR4JKJMKOR2LBis6CK3LEt/+tOf1NXVpcrKSr300ks6cOCAvva1rzkQD5hfvAMD8nd0qCKZzI5VJBKyenspc1wzR5ZWpqamdOHCBWUyGV24cEH19fVOnBaYd2pisZwSl6SKZFI1sZhLiTAfFH1F7vf79a1vfUtbtmxRZWWlVq5cqZUrV+YdF4/HFY/HJUnRaFTBYDA/jM9XcLyckXnumZZXunxmn2UVPL7Kslz/O86neS5nc5G56CL/9NNPdfDgQe3atUvV1dV66aWXtH//ft111105x4XDYYXD4ezjQp9sWmif0nKLaZlNyytdPnOd36/qAsdP+P0ac/nvOJ/muZyV5Sc7+/v7tWTJEt10003y+Xxas2aNjh07VuxpgXlpPBLRZEtLzthkS4vGIxGXEmE+KPqKPBgM6vjx4zp//rwqKyvV39+v5cuXO5ENmHcyoZCs3t6Ld60MDyvT0MBdKyha0UXe1tamO++8U88++6y8Xq9uvvnmnCUUALkyoZDGurvdjoF5xJH7yDdu3KiNGzc6cSoAwFXik50AYDiKHAAMR5EDgOEocgAwHEUOAIajyAHAcBQ5ABiOIgcAw1HkAGA4ihwADEeRA4DhHPmuFaCcsUcm5juKHPMae2RiIWBpBfMae2RiIaDIMa95h4YKjw8PlzgJMHcocsxrmcbGwuMNDSVOAswdihzzGntkYiHgzU7Ma+yRiYWAIse8xx6ZmO9YWgEAw1HkAGA4ihwADEeRA4DhKHIAMBxFDgCGo8gBwHAUOQAYjiIHAMNR5ABgOIocAAxHkQOA4Rz50qzPPvtMu3fv1kcffSSPx6MtW7ZoxYoVTpwaADADR4q8p6dHX/7yl/X9739f6XRa58+fd+K0AIBZKHpp5ezZs/rggw/09a9/XZLk8/l0ww03FB0MADA7Htu27WJO8OGHH+rVV1/V0qVLlUwm1draqscee0xVVVU5x8XjccXjcUlSNBrVhQsX8s7l8/mUTqeLiVNyZJ57puWVyFwqCy1zZWVlwfGii/zkyZP64Q9/qG3btqmtrU09PT1atGiROjo6rvjnBgcH88aCwaBSqVQxcUqOzHPPtLwSmUtloWVuamoqOF700kogEFAgEFBbW5sk6c4779SpU6eKPS0AYJaKLvK6ujoFAoHsFXZ/f7+WLl1adDAAwOw4ctfK448/rp07dyqdTmvJkiXq7Ox04rQAgFlwpMhvvvlmRaNRJ04FALhKfLITAAxHkQOA4ShyADAcRQ4AhqPIAcBwFDkAGI4iBwDDUeQAYDiKHAAMR5EDgOEocgAwnCPftQIAThkY8CoWq9HQkFeNjRlFIuMKhTJuxyprFDmAsjEw4FVHh1/JZEV2LJGoUG+vRZlfAUsrAMpGLFaTU+KSlExWKBarcSmRGShyAGVjaMhbcHx4uPA4LqLIAZSNxsbCyycNDSyrXAlFDqBsRCLjammZzBlraZlUJDLuUiIz8GYngLIRCmXU22spFqvR8LBXDQ3ctTIbFDmAshIKZdTdPeZ2DKOwtAIAhqPIAcBwFDkAGI4iBwDDUeQAYDiKHAAMR5EDgOEocgAwHEUOAIajyAHAcBQ5ABjOsSKfmppSJBJRNBp16pRAWfEODKjuyScVeOgh1T35pLwDA25HAiQ5+KVZf/zjH9Xc3Kxz5845dUqgbHgHBuTv6FBFMpkdq0gkZPX2KhMKuZgMcOiKfHR0VIlEQuvXr3fidEDZqYnFckpckiqSSdXEYi4lAv6fI1fkr732mh5++OErXo3H43HF43FJUjQaVTAYzA/j8xUcL2dknnvlkNdnWQXHqyyLf8suIvPn5yz2BH19faqtrVVra6uOHDly2ePC4bDC4XD2cSqVyjsmGAwWHC9nZJ575ZC3zu9XdYHxCb9fY/xbds1Cy9zU1FRwvOgiP3r0qN5//30dOnRIFy5c0Llz57Rz50499dRTxZ4aKBvjkYgqEomc5ZXJlhaNRyIupgIuKrrIN2/erM2bN0uSjhw5ot///veUOOadTCgkq7dXNbGYvMPDyjQ0aDwS4Y1OlAW2egNmKRMKaay72+0YQB5Hi/z222/X7bff7uQpAQAz4JOdAGA4ihwADEeRA4DhKHIAMBxFDgCGo8gBwHAUOQAYjiIHAMNR5ABgOIocAAxHkQOA4ShyADAcRQ4AhqPIAcBwFDkAGI4iBwDDUeQAYDiKHAAMR5EDgOEocgAwHEUOAIajyAHAcBQ5ABiOIgcAw1HkAGA4ihwADEeRA4DhKHIAMBxFDgCGo8gBwHAUOQAYzlfsCVKplHbt2qWxsTF5PB6Fw2Hde++9TmQDAMxC0UXu9Xr1yCOPqLW1VefOndNzzz2nO+64Q0uXLnUiHwBgBkUvrdTX16u1tVWStGjRIjU3N8uyrKKDAQBmp+gr8ulGRkZ06tQp3XrrrXnPxeNxxeNxSVI0GlUwGMwP4/MVHC9nZJ57puWVyFwqZL7IY9u27cSJJiYm9JOf/EQPPvig1qxZM+Pxg4ODeWPBYFCpVMqJOCVD5rlnWl6JzKWy0DI3NTUVHHfkrpV0Oq0dO3Zo3bp1sypxAIBzii5y27a1e/duNTc367777nMiEwDgKhS9Rn706FHt379foVBIzzzzjCRp06ZNam9vLzocAGBmRRf5l770Jb355ptOZAEAXAM+2QkAhqPIAcBwFDkAGI4iBwDDUeQAYDiKHAAMR5EDgOEocgAwHEUOAIajyAHAcBQ5ABjO0Y0lAJhnYMCrWKxGluWT31+nSGRcoVDG7Vi4ChQ5sIANDHjV0eFXMlnx+Ui1EokK9fZalLlBWFoBFrBYrGZaiV+UTFYoFqtxKRGuBUUOLGBDQ96C48PDhcdRnihyYAFrbCy8fNLQwLKKSShyYAGLRMbV0jKZM9bSMqlIZNylRLgWvNkJLGChUEa9vdbnd61Uye+f4K4VA1HkwAIXCmXU3T2mYDCoVGrM7Ti4BiytAIDhKHIAMBxFDgCGo8gBwHAUOQAYjiIHAMNR5ABgOIocAAxHkQOA4ShyADAcRQ4AhnPku1YOHz6snp4eTU1Naf369br//vudOG0O78CAamIxeYeGlGls1Hgkokwo5PjPAQDTFF3kU1NT2rNnj370ox8pEAho69atWrVqlZYuXepEPkkXS9zf0aGKZDI7VpFIyOrtpcwBLHhFL62cOHFCjY2NamhokM/n09q1a3Xw4EEnsmXVxGI5JS5JFcmkamIxR38OAJio6Ctyy7IUCASyjwOBgI4fP553XDweVzwelyRFo1EFg8H8MD5f4XHLKvizqyyr4PGldLnM5cy0zKbllchcKmT+/JzFnsC27bwxj8eTNxYOhxUOh7OPU6lU3jEXvw85f7zO71d1gZ894fdrrMDxpXS5zOXMtMym5ZXIXCoLLXNTU1PB8aKXVgKBgEZHR7OPR0dHVV9fX+xpc4xHIppsackZm2xp0Xgk4ujPAQATFV3ky5cv18cff6yRkRGl02kdOHBAq1atciJbViYUktXbq7MPPKDza9fq7AMP8EYnAHyu6KUVr9erxx9/XNu3b9fU1JTuueceLVu2zIlsOTKhkMa6ux0/LwCYzpH7yNvb29Xe3u7EqQAAV4lPdgKA4ShyADAcRQ4AhqPIAcBwFDkAGI4iBwDDUeQAYDiKHAAMR5EDgOEocgAwHEUOAIZz5LtWAGBgwKtYrEZDQ141NmYUiYwrFMq4HWtBoMgBFG1gwKuODr+SyYrsWCJRod5eizIvAZZWABQtFqvJKXFJSiYrFIvVuJRoYaHIARRtaMhbcHx4uPA4nEWRAyhaY2Ph5ZOGBpZVSoEiB1C0SGRcLS2TOWMtLZOKRMZdSrSw8GYngKKFQhn19lqKxWo0POxVQwN3rZQSRQ7AEaFQRt3dY27HWJBYWgEAw1HkAGA4ihwADEeRA4DhKHIAMBxFDgCGo8gBwHAUOQAYjiIHAMNR5ABgOIocAAxHkQOA4Yr60qy9e/eqr69PPp9PDQ0N6uzs1A033OBUNgCYNy7taWpZPvn9dY5+O2RRRX7HHXdo8+bN8nq9+s1vfqN9+/bp4YcfdiQYAMwX+XuaVju6p2lRSysrV66U13txK6cVK1bIsqyiAwHAfDPXe5o69n3k77zzjtauXXvZ5+PxuOLxuCQpGo0qGAzmh/H5Co6XMzLPPdPySmQuFVMyW1bhqrWsKkfyz1jk27Zt09hY/pfFd3R0aPXq1ZKk3/3ud/J6vVq3bt1lzxMOhxUOh7OPU6lU3jHBYLDgeDkj89wzLa9E5lIxJbPfXyepusD4hFKp2W/G0dTUVHB8xiL/8Y9/fMXn//a3v6mvr0/PP/+8PB7PrAMBwEIRiYwrkajIWV5xck/TopZWDh8+rLfeeks//elPdf311zsSCADmm+l7mlpWlfz+ifK5a2XPnj1Kp9Patm2bJKmtrU1PPPGEI8EAYD65tKfpxeUgZ/c2LarIX375ZadyAACuEZ/sBADDUeQAYDiKHAAMR5EDgOE8tm3bbocAAFy7sroif+6559yOcNXIPPdMyyuRuVTIfFFZFTkA4OpR5ABgOO8LL7zwgtshpmttbXU7wlUj89wzLa9E5lIhM292AoDxWFoBAMNR5ABgOMd2CLoWXV1dGhwclCSdPXtW1dXVevHFF/OO+853vqOqqipdd9118nq9ikajpY6a9eabb+qvf/2rbrrpJknSpk2b1N7ennfc4cOH1dPTo6mpKa1fv173339/qaNmzXaTbLfneaY5s21bPT09OnTokK6//np1dna6uj6aSqW0a9cujY2NyePxKBwO695778055siRI4rFYlqyZIkkac2aNXrooYfciJs10+tcbvM8ODiorq6u7OORkRFt3LhR3/zmN7Njbs/zK6+8okQiodraWu3YsUOS9Omnn6qrq0v/+9//tHjxYn3ve9/TjTfemPdnHekKu0y8/vrr9m9/+9uCz3V2dtpnzpwpcaLC3njjDfutt9664jGZTMZ+8skn7aGhIXtyctL+wQ9+YH/00UclSpjv8OHDdjqdtm3btvfu3Wvv3bu34HFuzvNs5qyvr8/evn27PTU1ZR89etTeunWrK1kvsSzLPnnypG3btn327Fn7qaeeysv873//2/75z3/uRrzLmul1Lrd5ni6Tydjf/va37ZGRkZxxt+f5yJEj9smTJ+2nn346O7Z371573759tm3b9r59+wr+f+dUV5TF0opt2/rHP/6hr3zlK25HccSJEyfU2NiohoYG+Xw+rV27VgcPHnQtjwmbZM9mzt5//33ddddd8ng8WrFihT777DOdPn3apcRSfX199kp10aJFam5uLsu5vVrlNs/T9ff3q7GxUYsXL3Y7So7bbrst72r74MGDuvvuuyVJd999d8EOcKorXF1aueSDDz5QbW2tvvCFL1z2mO3bt0uSvvGNb+Ts/emGv/zlL9q/f79aW1v16KOP5r2AlmUpEAhkHwcCAR0/frzUMQuaaZNst+Z5NnNmWVbORrWBQECWZam+vr5kOS9nZGREp06d0q233pr33LFjx/TMM8+ovr5ejzzyiJYtW+ZCwlxXep3LeZ7//ve/X/aCr9zm+cyZM9k5q6+v1yeffJJ3jFNdMedFPpvNm6/04lw6h9/v15kzZ/Szn/1MTU1Nuu2221zJvGHDhuza2xtvvKFf//rX6uzszDnOLnBH51zvZ+rEJtmlnufpZjNnbszrbExMTGjHjh167LHHVF2du8HuLbfcoldeeUVVVVVKJBJ68cUXtXPnTpeSXjTT61yu85xOp9XX16fNmzfnPVeO8zwbTs31nBf5TJs3ZzIZvffee1d8Y83v90uSamtrtXr1ap04cWJOC2amzJesX79ev/jFL/LGA4GARkdHs49HR0fn/GrGiU2ySz3P081mzgKBQM6O6aWY15mk02nt2LFD69at05o1a/Ken17s7e3t2rNnjz755JPsm+VumOl1Lsd5lqRDhw7plltuUV1dXd5z5TjPtbW1On36tOrr63X69OmCWZzqCtfXyPv7+9XU1JTz68V0ExMTOnfuXPa///WvfykUCpUyYo7pa4XvvfdewV/fli9fro8//lgjIyNKp9M6cOCAVq1aVcqYOS5tkv3ss89edpNst+d5NnO2atUq7d+/X7Zt69ixY6qurna1YGzb1u7du9Xc3Kz77ruv4DFjY2PZq64TJ05oampKNTU1pYyZYzavc7nN8yVX+s293OZZujiP7777riTp3Xffzf5mPJ1TXeH6Jzt37dqltrY2bdiwITtmWZZeffVVbd26VcPDw/rlL38p6eLV+1e/+lU9+OCDbsXVyy+/rA8//FAej0eLFy/WE088ofr6+pzMkpRIJPT6669rampK99xzj6uZv/vd7yqdTmfX8i9tkl1u81xozt5++21J0oYNG2Tbtvbs2aN//vOfqqysVGdnp5YvX17SjNP997//1fPPP69QKJT9LWfTpk3Zq9kNGzboz3/+s95++215vV5VVlbq0Ucf1Re/+EXXMl/udS7neZak8+fPa8uWLeru7s5efU/P7PY8/+pXv9J//vMfjY+Pq7a2Vhs3btTq1avV1dWlVCqlYDCop59+WjfeeOOcdIXrRQ4AKI7rSysAgOJQ5ABgOIocAAxHkQOA4ShyADAcRQ4AhqPIAcBw/weq/rL58aqlVgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import style\n",
    "import math\n",
    "\n",
    "style.use(\"ggplot\")\n",
    "\n",
    "# we create our own data\n",
    "\n",
    "# features\n",
    "samples_x = [\n",
    "    [1,7],\n",
    "    [2,4],\n",
    "    [3,8],\n",
    "    [5,1],\n",
    "    [6,-1],\n",
    "    [7,5],\n",
    "    [10,-2],\n",
    "    [-7,0]\n",
    "]\n",
    "\n",
    "# classes/labels\n",
    "samples_y = [1,1,1,-1,-1,-1, -1, 1]\n",
    "\n",
    "[plt.scatter(xi[0], xi[1], color = \"r\" if yi == 1 else \"b\") for xi, yi in zip(samples_x, samples_y)]\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quick explaination of a SVM:\n",
    "======\n",
    "A support vector machine is a supervised machine learning algorithm that can be used for classification and regression. In this example we'll stick to classifaction.\n",
    "It works by finding the hyperplane (a plane with dimensions n-1 in a space of dimension n) which best divides the training dataset into its different classes while maximizing the distance of the dataset to the plane. The left graph shows a bunch of different hyperplanes, while the right picture shows the optimal hyperplane for this datset. The dotted lines are called support vectors.\n",
    "<table><tr><td><img src=https://miro.medium.com/max/300/0*9jEWNXTAao7phK-5.png width=\"350\"></td><td><img src=https://miro.medium.com/max/300/0*0o8xIA4k3gXUDCFU.png width=\"350\"></td></tr></table>\n",
    "\n",
    "So essentially this is a constrained optimization problem where we want to maximze the width between the support vectors.   \n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Math behind the SVM:\n",
    "===\n",
    "In this problem we will only consider linearly separable data with two classes: Positive(+) & Negative(-)\n",
    "<br>\n",
    "<br>\n",
    "### Descision Rule ###\n",
    "The normal equation of the hyperplane is:  \n",
    "$$\\vec{x} \\cdot \\vec{w} = \\vec{b} \\iff \\vec{x} \\cdot \\vec{w} - \\vec{b} = 0$$\n",
    "where $\\vec{w}$ is the normal vector of the plane, $\\vec{x}$ is some point in the plane and $\\vec{b}$ is some offset (otherwise the plane would always go through the origin).  \n",
    "Given some query $\\vec{x_*}$, we can easily calculate on which side of the hyperplane it is i.e. which class it belongs to:\n",
    "$$\n",
    "\\begin{align}\n",
    "\\vec{x_*} \\cdot \\vec{w} &- \\vec{b} > 0 \\implies \\vec{x_*} \\in + \\\\\n",
    "\\vec{x_*} \\cdot \\vec{w} &- \\vec{b} < 0 \\implies \\vec{x_*} \\in -\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Support Vectors ###\n",
    "For each sample $\\vec{x_i}$ in our training data we demand:  \n",
    "$$ \\mid \\vec{x_i} \\cdot \\vec{w} - \\vec{b} \\mid \\, \\geq 1 $$\n",
    "This ensures that ${\\left\\|w\\right\\|}$ is related to the margin, that is the size of the gap between support vectors. Furthermore for the the data points that lie on the support vectors i.e. the points nearest to the plane we want\n",
    "$$ \\mid \\vec{x_i} \\cdot \\vec{w} - \\vec{b} \\mid \\, = 1 $$\n",
    "which actually holds two equations\n",
    "$$ \n",
    "\\begin{align}\n",
    "\\vec{x_i} \\in + \\implies \\vec{x_i} \\cdot \\vec{w} & - \\vec{b} = 1 \\\\\n",
    "\\vec{x_i} \\in - \\implies \\vec{x_i} \\cdot \\vec{w} & - \\vec{b} = -1 \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "To make this more easier to deal with we'll introduce another vairable $y_i$ which will just hold the value of whatever class $\\vec{x_i}$ is:\n",
    "$$y_i = \n",
    "\\begin{cases}\n",
    "+1,& \\text{if } \\vec{x_i} \\in + \\\\\n",
    "-1,& \\text{else if } \\vec{x_i} \\in -\n",
    "\\end{cases}\n",
    "$$\n",
    "With this we can express the above in one equation\n",
    "$$\n",
    "\\begin{align}\n",
    "y_i(\\vec{x_i} \\cdot \\vec{w} &- \\vec{b}) -1 \\, \\geq 0 \\;\\; \\text{for all } \\vec{x_i} \\\\\n",
    "y_i(\\vec{x_i} \\cdot \\vec{w} &- \\vec{b}) -1 \\, = 0 \\;\\; \\text{for support vectors} \\tag{1}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Margin ###\n",
    "The Margin, can be calculated as:\n",
    "$$ \n",
    "\\text{Margin} = \\left(\\vec{x_+}-\\vec{x_-}\\right) \\cdot \\frac{\\vec{w}}{\\left\\|w\\right\\|} = \n",
    "\\frac{\\vec{x_+} \\cdot \\vec{w} - \\vec{x_-} \\cdot \\vec{w}}{\\left\\|w\\right\\|}\n",
    "$$\n",
    "where $\\vec{x_+} \\in +$ and  $\\vec{x_-} \\in -$ are samples in the training set that lie on the support vectors and $\\frac{\\vec{w}}{\\left\\|w\\right\\|}$ is the unit normal vector.  \n",
    "With (1) evaluated for $\\vec{x_+}$ and $\\vec{x_-}$\n",
    "$$\n",
    "\\begin{align}\n",
    "&y_+(\\vec{x_+} \\cdot \\vec{w} - \\vec{b}) -1 = \\vec{x_+} \\cdot \\vec{w} - \\vec{b} -1 = 0 \\implies \\vec{x_+} \\cdot \\vec{w} = 1+b \\\\\n",
    "&y_-(\\vec{x_-} \\cdot \\vec{w} - \\vec{b}) -1 = -\\vec{x_-} \\cdot \\vec{w} + \\vec{b} -1 = 0 \\implies \\vec{x_-} \\cdot \\vec{w} = -1+b\n",
    "\\end{align}\n",
    "$$\n",
    "we have \n",
    "$$ \\text{Margin} = \\frac{1+b-(-1+b)}{\\left\\|w\\right\\|}=\\frac{2}{\\left\\|w\\right\\|}$$\n",
    "So if we want to maximize the margin, we need to minimize $\\left\\|w\\right\\|$.  \n",
    "<img src=https://upload.wikimedia.org/wikipedia/commons/thumb/7/72/SVM_margin.png/1024px-SVM_margin.png width=\"350\">\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem formulation ###\n",
    "This is a constrained optimization problem:\n",
    "We want to minimze \n",
    "$$ \\frac{1}{2} \\left\\|\\vec{w}\\right\\|^2$$\n",
    "such that \n",
    "$$ y_i(\\vec{x_i} \\cdot \\vec{w} - \\vec{b}) -1 \\geq 1 \\iff  y_i(\\vec{x_i} \\cdot \\vec{w} - \\vec{b}) \\geq 0$$  \n",
    "The reason why we are saying we want to minimze $ \\frac{1}{2} \\left\\|\\vec{w}\\right\\|^2$ instead of $\\left\\|w\\right\\|$ is that this makes things nicer to compute later on.  \n",
    "<br>\n",
    "Using the method of [lagrange multipliers](https://en.wikipedia.org/wiki/Lagrange_multiplier) gives us the **primal form** of this problem\n",
    "(if you dont know what the lagrangian or lagrange multipliers are, [this](https://www.khanacademy.org/math/multivariable-calculus/applications-of-multivariable-derivatives/lagrange-multipliers-and-constrained-optimization/v/constrained-optimization-introduction) video series from will give you an intuitive introduction) \n",
    "$$ \\underset {w, b}{\\operatorname {arg\\,max} } \\ L_p(w,b)=\\frac{1}{2} \\left\\|\\vec{w}\\right\\|^2 - \\sum \\alpha_i \\left[ y_i \\left(\\vec{x_i} \\cdot \\vec{w} - \\vec{b} \\right) -1 \\right]$$\n",
    "<br>\n",
    "where the $\\alpha_i$ are lagrange multipliers. Setting $\\nabla L = 0$ gives us  \n",
    "<br>\n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{\\partial L}{\\partial w} &= \\vec{w}- \\sum \\alpha_i y_i \\vec{x_i} = 0  \\iff \\vec{w} = \\sum \\alpha_i y_i \\vec{x_i} \\\\\n",
    "\\frac{\\partial L}{\\partial b} &= \\sum \\ \\alpha_i y_i = 0 \n",
    "\\end{align}\n",
    "$$\n",
    "Plugging these back in to $L_p$  \n",
    "$$\n",
    "\\begin{align}\n",
    "L_p(w,b) &=\\frac{1}{2} \\left\\|\\vec{w}\\right\\|^2 - \\sum \\alpha_i \\left[ y_i \\left(\\vec{x_i} \\cdot \\vec{w} - \\vec{b} \\right) -1 \\right] \n",
    "\\\\[7pt]\n",
    "&= \\frac{1}{2} \\left(\\sum \\alpha_i y_i \\vec{x_i}\\right)^2-   \\sum \\alpha_i \\left[ y_i \\left(\\vec{x_i} \\left(\\sum \\alpha_j y_j \\vec{x_j}\\right) - b \\right) - 1\\right] \n",
    "\\\\[7pt]\n",
    "&= \\frac{1}{2}\\left(\\sum \\alpha_i y_i \\vec{x_i}\\right)\\left(\\sum \\alpha_j y_j \\vec{x_j}\\right)   -   \\sum \\alpha_i y_i \\vec{x_i} \\left(\\sum \\alpha_j y_j \\vec{x_j}\\right)    - \\left(\\sum \\alpha_i y_i b \\right)    + \\sum \\alpha_i \n",
    "\\\\[7pt]\n",
    "&= \\frac{1}{2}\\sum_{i}\\sum_{j}\\alpha_i\\alpha_j y_i y_j \\left(\\vec{x_i} \\cdot \\vec{x_j} \\right) - \\frac{1}{2}\\sum_{i}\\sum_{j}\\alpha_i\\alpha_j y_i y_j \\left(\\vec{x_i} \\cdot \\vec{x_j} \\right)  + \\sum \\alpha_i \n",
    "\\\\[7pt]\n",
    "=L_d(\\alpha) &= \\sum \\alpha_i - \\frac{1}{2}\\sum_{i}\\sum_{j}\\alpha_i\\alpha_j y_i y_j \\left(\\vec{x_i} \\cdot \\vec{x_j} \\right)\n",
    "\\end{align}\n",
    "$$\n",
    "This reformulation is called the **dual problem**. This is a convex optimization problem which can be solved by quadratic programing.  \n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimization (SMO) ###\n",
    "Lets look at whats actually going on when we are changing our $\\alpha_i's$. Lets Consider $\\vec{w}$:\n",
    "$$ \\vec{w} = \\sum \\alpha_i y_i x_i = \\underbrace{\\sum \\alpha^+_i \\ x_i^+}_{h^+} - \\underbrace{\\sum \\alpha^-_i \\ x_i^-}_{h^-}$$\n",
    "where $x_i^+ \\in + $ and $x_i^- \\in -$. The superscript on $\\alpha_i$ just denotes if it belongs to a positive or negative class. Other than that it doestn have any mathematical meaning.  \n",
    "$h^+$ and $h^-$ are weighted averages of their respective classes, where the $\\alpha_i's$ are the weights. By changing the $\\alpha_i's$ we can move $h^\\pm$ around in the polytope created by its respective class. The SVM will be optimized when $h^+$ and $h^-$ are closest. In this case almost all $\\alpha's$ will be zero except for those that belong to values on the support vectors. Of course $h^+$ and $h^-$ are not allowed to exist outside of their polytopes, so we will need to restrict the weights. With this our objective is:\n",
    "Minimize \n",
    "$$ \\text{Minimize}\\left\\|h^+-h^-\\right\\| \\text{  s.t.  } \\alpha_i \\geq 0 \\text{ and } \\sum\\alpha_i^+=\\sum\\alpha_i^-=1 $$\n",
    "**This is how to do it:**  \n",
    " Repeat until optimized i.e. $h^+$ and $h^-$ cannot move any closer\n",
    "- 1: Pick two alphas from the same class:\n",
    "$$\\text{Without loss of generality: }\\alpha_j \\text{, } \\alpha_k \\in +$$\n",
    "- 2: Change them as follows\n",
    "$$\n",
    "\\begin{align}\n",
    "&c = \\alpha_j + \\alpha_k \\\\\n",
    "&\\alpha_j' = \\frac{(x_j-x_k)\\cdot \\left(h^--h^++cx_k+\\alpha_j x_j + \\alpha_k x_k\\right)}{\\|x_j-x_k\\|^2}\\\\\n",
    "&\\alpha_k' = c - \\alpha_j'\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "**How it works:**\n",
    "\n",
    "$\\left\\|h^+ - h^-\\right\\|$ can be viewed as a (fairly easy) geometric problem\n",
    "$$\n",
    "\\begin{align}\n",
    "\\left\\|h^+ - h^-\\right\\| &= \\Bigg\\| \\sum \\alpha^+_i x_i^+ - \\sum \\alpha^-_i x_i^-\\Bigg\\| \n",
    "\\\\\n",
    "&= \\Bigg\\| \\alpha_jx_j+\\alpha_kx_k+\\sum_{i \\neq j,k} \\alpha^+_i \\ x_i^+ - \\sum \\alpha^-_ix_i^- \\Bigg\\| \n",
    "\\\\\n",
    "&= \\Bigg\\| \\alpha_jx_j+(\\underbrace{c}_{\\alpha_j+\\alpha_k}-\\alpha_j)x_k+\\sum_{i \\neq j,k} \\alpha^+_i \\ x_i^+ - \\sum \\alpha^-_ix_i^- \\Bigg\\| \n",
    "\\\\\n",
    "&= \\Bigg\\| \\underbrace{\\alpha_j}_{a} \\underbrace{(x_j-x_k)}_{\\vec{x}} + \\underbrace{cx_k + \\sum_{i \\neq j,k} \\alpha_i^+ x_i^+}_{\\vec{u}}-\\underbrace{\\sum\\alpha_ix_i}_{\\vec{h}}\\Bigg\\| \n",
    "\\\\\n",
    "&=\\left\\|a\\vec{x}+\\vec{u}-\\vec{h} \\right\\| \\equiv \\left\\| g-\\vec{h} \\right\\|\n",
    "\\end{align}\n",
    "$$\n",
    "<br>\n",
    "where $g:\\vec{u}+a\\vec{x}$ is the line in $\\vec{x}$ direction that passes through $\\vec{u}$.The distance between $g$ and $\\vec{h}$ is minimal at the point on the line where \n",
    "$$\n",
    "\\begin{align}\n",
    "g-\\vec{h} \\perp \\vec{x} &\\implies (g-\\vec{h})\\cdot\\vec{x}=\\vec{x}\\cdot\\vec{u}+a\\vec{x}\\cdot\\vec{x}-h\\vec{x}=0 \\\\\n",
    "&\\implies a = \\frac{\\vec{x}\\cdot(\\vec{h}-\\vec{u})}{\\|x\\|^2}\n",
    "\\end{align}\n",
    "$$\n",
    "<br>\n",
    "Plugging everything we get our formula for $\\alpha_j'$\n",
    "$$\n",
    "\\begin{align}\n",
    "\\alpha_j' &= \\frac{(x_j-x_k)\\cdot \\bigg(\\overbrace{\\sum \\alpha^-_i x_i^-}^{h^-}-\\bigg(cx_k + \\overbrace{\\sum_{i \\neq j,k} \\alpha_i^+ x_i^+}^{h^+-\\alpha_jx_j-\\alpha_kx_k}\\bigg)\\bigg)}{\\|x_j-x_k\\|^2} \n",
    "\\\\[7pt]\n",
    "&=\\frac{(x_j-x_k)\\cdot \\left(h^--\\left(c x_k + h^+-\\alpha_jx_j-\\alpha_kx_k\\right)\\right)}{\\|x_j-x_k\\|^2} \n",
    "\\\\[7pt]\n",
    "&= \\frac{(x_j-x_k)\\cdot \\left(h^--h^++cx_k+\\alpha_j x_j + \\alpha_k x_k\\right)}{\\|x_j-x_k\\|^2}\n",
    "\\end{align}\n",
    "$$\n",
    "<br>\n",
    "Once we have $\\alpha_j'$ its very simple to calculate $\\alpha_k'$ since $\\alpha_j + \\alpha_k \\overset{!}{=} \\alpha_j' + \\alpha_k'$ otherwise $\\sum \\ \\alpha_i y_i \\ne 0$:\n",
    "$$\n",
    "\\begin{align}\n",
    "\\\\\n",
    "\\alpha_k' &= (\\alpha_j + \\alpha_k) - \\alpha_j' \\\\\n",
    "&= c - \\alpha_j'\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iteration 0: h+ & h- distance:  4.024922359499621\n",
      "iteration 1: h+ & h- distance:  4.024922359499621\n",
      "(0, 7, 9.0)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMsAAAD7CAYAAADEgWCeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAZgElEQVR4nO3da1BUZ5oH8H/TDV4C3dDdCEsEEjRuIJJxENwJwywaW5MZE+/xMplYM7lYDpodJRlMdKJTi5gGo1BeqJmxojXBZPCyoqUbo9u6pZuQWqjqiQYhMcQCNCpIdwvNiAG6z37I0oo02N3nPX0u/fw+wYF+z1Mkj/9z3vdcVBzHcSCEPFCY2AUQIhfULIT4iJqFEB9RsxDiI2oWQnxEzUKIjzQsBikvL4fVaoVOp8PWrVsBAAcOHMDp06eh1WoBAEuXLkVGRgaL3REiCibNMnXqVDz77LPYtWvXgO2zZs3C7NmzWeyCENExOQxLS0tDZGQki6EIkSwmyTKUkydP4ty5c0hJScGyZct8aqhr164JWRKMRiPa29sF3QcrVKswhqs1ISFhyM+pWF3u0tbWhuLiYs85y61btzznK/v374fD4UBeXt6gz1ksFlgsFgCA2WxGT08Pi3KGpNFo0NfXJ+g+WKFahTFcrREREUN/TqiCoqOjPV9Pnz4dxcXFXn/PZDLBZDJ5vhf6Xyel/AsoNUqpdbhkEWzq2OFweL6uqalBYmKiULsiJCiYJEtZWRnq6+vhdDqxYsUKLFq0CBcvXkRTUxNUKhViY2OxfPlyFrsiRDRMmmX16tWDtj399NMshiZEMmgFnxAfUbMQ4iNqFiIJFRUVOHHihNhlDIuahYjObrdj06ZNOHLkiNilDIuahYhu9+7d6Orqwpo1a8QuZVjULERUdrsde/bswXPPPYfHH39c7HKGRc1CRCWXVAGoWYiI5JQqADULEZGcUgWgZiEikVuqANQsRCRySxWAmoWIQI6pAlCzEBHIMVUAahYSZHJNFYCahQSZXFMFoGYhQSTnVAGoWUgQyTlVAGoWEiRyTxWAmoUEidxTBaBmIUGghFQBqFlIECghVQBqFiIwpaQKQM1CBKaUVAGoWYiAlJQqADULEZCSUgWgZiECUVqqANQsRCBKSxWAmoUIQImpAlCzEAEoMVUAahbCmFJTBaBmIYwpNVUAahbCkJJTBaBmIQwpOVUAahbCiM1mU3SqAIxek1deXg6r1QqdTud5tXdXVxdKS0tx8+ZNxMbGYs2aNYiMjGSxOyJB27dvV3SqAIySZerUqVi3bt2AbUeOHEF6ejq2b9+O9PR0yb97gwTObrdj165dik4VgFGzpKWlDUqN2tpa5ObmAgByc3NRW1vLYldEgnbv3g2n06noVAEEPGfp6OhATEwMACAmJgadnZ1C7YqIqH8GbMGCBYpOFYDROQsfFosFFosFAGA2m2E0GgXdn0ajEXwfrMih1h07dqCrqwsbNmyQfK39Av27CtYsOp0ODocDMTExcDgc0Gq1Xn/PZDLBZDJ5vm9vbxeqJACA0WgUfB+sSL1Wu92OnTt3es5VpFzrvYb7uyYkJAz5OcEOwzIzM3H27FkAwNmzZ5GVlSXUrohIlL6ucj8myVJWVob6+no4nU6sWLECixYtwty5c1FaWoozZ87AaDQiPz+fxa6IRCh9td4bJs2yevVqr9s3bNjAYngiQaGWKgCt4JMAsE6VzTWbcar5FIPKhEXNQvzGMlUu2i5i1/ldqGuvY1CZsKhZiF9Yp0qZtQzaCC1emfgKg+qERc1C/MI6VT5u+hivTXwNuhE6BtUJi5qF+CyUUwWQwAo+kQ/WM2CmZBNyHs6RRaoA1CzER0KsqyyesJjJOMFCh2HEJyxTpd5Wj91f7kZ3XzeDyoKHmoU8EOtU2Wbdhm3Wbehx9TCoLnjoMCyEqVtaEFVSAvWNG3DFx8NZUABXUtKg32OZKnW2OpxoOoH8jHzZnKv0o2YJUeqWFuiXLEF4c7NnW7jVCntl5YCGEWoG7NWJr/IeK9joMCxERZWUDGgUAAhvbkZUScmAbUKkyqsTX5VdqgDULCFLfeOG9+2trZ6vWaeK2+3G1LFTZZkqAB2GhSxXfLz37XFxnq9Zr6s8GfskPvz5h0zGEgMlS4hyFhSgNzl5wLbe5GQ4CwoAsE+V/V/vx83bN3mPIyZqlhDlSkqCvbISt+fNw/fZ2bg9b96Ak3vW5yr55/JR0VDBeywx0WFYCHMlJeHWzp2Dtof6NWBDoWaRIV/XRwJF6yreUbPIjK/rI4GidZWh0TmLzPi6PhIolqnS5+5DhDoCy9OXyz5VAEoW2fFlfSRQrFNFE6ZB+dPl4DiO91hSQMkiM76sjwSKZapc7riMr+1fAwBUKhXv8aSAmkVmHrQ+EijWqfJuzbuYf3y+7C7DHw4dhslM//pIVEkJ1K2tcMXFMZkNE+Le+vyMfIzSjOI9nlRQs8jQUOsjgaIZMN/QYRgRJFXkemXxcKhZQhzrVGmwNyBudJziUgWgw7CQx/rK4oWPLcTzKc9jhHoEk/GkhJIlhLFOlTpbHTiOU2SjANQsIY1lqpxvPY9nDj+DfV/tY1CZNFGzhCjWqVL0aRG0EVrMTpnNoDppomYJUayvLD566agiZ8DuRc0SgoRYV9GN0ClyBuxe1CwhiGWqdPZ0wtpmxetZrys6VYAgTB2vXLkSI0eORFhYGNRqNcxms9C7JMNgnSraCC0+W/wZ9Ho9ujuVcx2YN0FZZ9m4ceOQr/YmwcUyVWzdNmhHaDFKMwoPRTyEbgTWLC0tapSUROHGDTXi410oKHAiKcnFuz7WaFEyhLBOlbWfrsUV5xV8Mu+TgMdoaVFjyRI9mpsvAUgCEAWrNRyVlXbJNUxQzlmKioqwdu1aWCyWYOyODEGIe+tnJs/kdb9KSUkUmpvDACwCMAsA0NwcjpKSKN41siZ4shQWFkKv16OjowObNm1CQkIC0tLSPD+3WCyeJjKbzTAajYLWo9FoBN8HKyxrtdls2Lt3LxYsWICcnBze45WfLYduhA5rc9ciemR0wLXa7RoABwHUA6i8Z/tIwf47BVqr4M2i1+sBADqdDllZWWhsbBzQLCaTCSaTyfN9e3u7oPUYjUbB98EKy1qLi4vhdDqRl5fHe8z+dZX8jHz0dfWhvas94FpjYqIA/DuANAALPdv1+jtob7/Fq86hDFdrQkLCkJ8T9DDszp076O7u9nx94cIFJDF8ZA/xDetzlYOXDjK7X2Xy5I/wQ6psAKAGACQn96KgwMl7bNYETZaOjg689957AACXy4WcnBxMmjRJyF0SL1hfWbzxJxvxUupLvNdVXC4XKiq24dFHJ+BHP/oF2tq+R1xciM6GxcXFYcuWLULugjwA61Tp7uvGKM0ojI8ez3us48eP49KlSygvL8ecOdJLkvvRCr5ChX97GbGvLMffZsxEV1cXfr/wBd5jXrRdROZHmai+Vs17LJfLhbKyMkyYMAHPPfcc7/GCgdZZFCj828swLl0K+3dXsQPACwB+tn492lNS0DsuJeBxy6xlcHNuPGF4gneN96aKWq3mPV4wULIoUPRmM1TfXYVZOxpOAO8AUH13FdGbA7/UqP/e+tcmvsbkXEVuqQJQsiiSuq0V7QD+1HUHs+JjYHj8YXQ1Xseom4E/tZLlk/DlmCoAJYsiucbEoQxAl9uNP37fA3VPHxxPPoLWMWFw/88pcH19fo3X3Nkc8qkCULMo0uWVv8V2lQovAMh0/ANjqr+C4bID7oSHwX2wE+4/rPCraZK1yTg25xjTVFm9erWsUgWgwzBFKv+vU3ByHN7+6b/Cfacbrtgx+Me6t6BKeRSqOivcx/4G7oOd4P7zAFSzFkH11NNQabz/r8BxHFQqFTLGZPCuS86pAlCzKM696ypxf/4z7n3mvgoA0icjbGIG4GPT/PbMb5HwUAI2/GQD79rkeq7Sjw7DFMaX1XqVSgVV+mSEvb0FYf+2EdBGez08q7PV4djlY3go/CHedck9VQBKFkXxd7VepVINmzSltyuZXQMm91QBqFkUJdBrwLw1Td2hEnwysQFrDHOgVfNLFiWkCkCHYYrB4hqwew/PdvxrDLTucLx8qsnv2bP7yXkG7F6ULArB8spilUqFDTO3ocHWgJhsg1+zZ/dTSqoA1CyKwPrKYgBIjEpEYlQiAPg1e3Y/JZyr9KPDMAVgfW/9iydexBXnFc82X2fP7qekVAEoWWRPiKdLWtus0EYMfnTVcLNntxe/DC59yoCkUVKqANQssifEE1vyM/KHvQbMW9M4y82AYQxUsxbhas4MFGu1OLljB6JSU/HknDm8a5MCahYZE/tdkPc2jbalEbc+/DO4D3ZC9cl/oO2f/hndDQ1AZSVejI1Fpd2OJJf0bhX2B52zyJgQqRLIk/BVKhVGTH4KYW9vwZ/WleCaNho3d+xCQqweWLgQzeHhKImS3nPA/EXNIlOsUyUpKgkFmQW8VutVKhX+O/OnmP/ULHzTdQdR+WuA/z9XaaVzFiIW1k9s0UZo8bsf/473OHE9PUBhIZCWhq9///bd7TI/BAMoWWSJ+Vu7/rcIp5pPMagMmPzRR0B9PbBhgydVknt7UeCU/tNbHoSaRYZYn6uUXyjHl+1f8h7L5XKhYts2PDphAub+4hfI/v57zLt9WxEn9wAdhsmO2DNgwxnwHDAFJMn9KFlkRiozYPdT2mq9N5QsMiLlVDl8+LCiVuu9oWaREdYzYDOSZ+BnD/+MV6q0tKhRXDwap05tRlRUKp58Uhmr9d5Qs8iEEFcWL56wmNfn77616zD636/y4ouxknxrFwt0ziITLFOl3laPv3z5F3T38Xth6t23dt19v4pU39rFAjWLDLBOlW3WbSi1lqLH1cNrnBs31AAO4f73q7S2KvOchZpFBqQ6AxYX1wNvb+2Ki1PeIRhAzSJ5Up4Bk9Nbu1igE3yJE+N+FV/c/9Yuh8MNvf6OZN/axQI1i4TZbDamqcJxHKaNncZ+tX6OE0bjCMFemCoVgjfLF198gb1798LtdmP69OmYO3eu0LtUjO3btzNdV0k3pmPfz/fxHicUVuu9EfScxe124/3338e6detQWlqKzz77DFevXvV7HHVLC6JXrYJh4UJEr1oFdUuLANVKi91ux65du5ilSuXXlWi73cagMuU8B8xfgjZLY2Mj4uPjERcXB41Gg+zsbNTW1vo1hrqlBfolSzC6qgojPv8co6uqoF+yRPENs3v3bjidTmbnKm+cewMVDRW8xwrVVAEEbha73Q6DweD53mAwwG63+zVGVEkJwpubB2wLb25GVEkJkxqlqLOzE3v27MGCBQskNwMWqqkCCHzOwnHcoG0qlWrA9xaLBRaLBQBgNpthNBoH/FwzRHONtNsH/a4vNBpNQJ8LJoPBgEOHDiEpKYl3redbz+NE0wn8IecPGPfwOF5juVwu7NixA6mpqfj1r389oFnk8HftF2itgjaLwWCAzWbzfG+z2RATEzPgd0wmE0wmk+f79vb2AT+P1usx2svYd/R63Lrvd31hNBoH7UOK0tPTmdS68cxGaCO0+GXKL3mPdfToUTQ0NKC8vBwOh2PAz+TydwWGrzUhIWHIzwl6GDZu3Dhcv34dbW1t6OvrQ3V1NTIzM/0aw1lQgN7k5AHbepOT4SwoYFmqIvW5+zBCPQLL05fT/SoMCJosarUaL7/8MoqKiuB2uzFt2jQkJib6NYYrKQn2ykpElZRA3doKV1wcnAUFcCUlCVS1cmjCNNj19C6vh8P+UtrTJQOh4lj8JRm6du2aoOMr5XDhQS53XEaPqweP6/lPELhcLs+hssVi8dosSvm7inYYRsTzbs27WHB8Ae/L8AF+M2AtLWqsWhWNhQsNWLUqGi0t8k0lutxFgS7aLuLjpo/xRsYbGKUZxWssPucqd28OC/dss1rDZXtzGCWLAvWvq4j93vofbg4LH7BNzjeHUbMoTH+qvDbxNdFnwH64OWwwud4cRodhCvOV/SvEj45nmiqBzoDFx3s/1JLrzWE0GyZhgdba4+pBhDqC1759mQG7l7davZ2zJCf3in7OQrNhBF+2fwmO43g3CsDmGrCkJBcqK+2YN+82srO/x7x5t0VvFD7oMEwh6mx1eLbqWZhzzHgp9SVeY7FcrU9KcmHnTmXcFEbNohD9M2CzU2bzHotW672jwzAFoGcW+6Z/gXTmTE1AC6SULAog1JPwlZQqgycbRvu9QErJInOdPZ34e9vfKVUegMUCKSWLzGkjtPh08ad0ZfEDsFggpWaRMVu3DdoRWt7XfwHKThWAzQIpHYbJ2NpP12LWkVlMU0Wp99YXFDiRnNw7YJu/T8+kZpGp/hmwZ5KfGfRcA38pPVWAgQukubnugBZI6TBMpmgGzH/9C6Q/XO7i/0IpJYsM0bqKOChZZOjQpUOUKiKgq44lbKha3ZwblzsuY3z0eF7j+3tl8XCU8HcFhr/qmJJFZrr7ujFKM4p3owCUKv6icxYZqbPVYfKHk1F9rZr3WHSu4j9KFhkps5aBA4cnDE/wHotSxX+ULDJBM2Dio2SRCVpXER8liww0dzZTqkgAJYsMJGuTcXzOcTyqe5T3WJQqgaNmkTg350aYKgw/HvNj3mNRqvBDh2ESl3cmD3/8/I9MxlL6lcVCo2aRsPOt53Hs8jFERfB/3CmlCn90GCZhRZ8W0QyYhFCySFSdrQ5HLx2lGTAJoWSRqO1/3w7dCB2lioRQs0jUhn/ZgKt9VylVJESwZjlw4ABOnz4NrVYLAFi6dCkyMjKE2p3ijI0ai0nGSbwve6dUYUfQZJk1axZmz+b/OFESGEoVtugwTMEoVdgStFlOnjyJc+fOISUlBcuWLUNkZKSQuyP3oFRhj9dtxYWFhbh1a/BTMpYsWYLHHnvMc76yf/9+OBwO5OXlDfpdi8UCi8UCADCbzejp6Qm0HJ9oNBr09fUJug9W+NR68OBB/OpXv8K+ffvwwgsvMK5sMKX8XSMihn63TVDuwW9ra0NxcTG2bt36wN+le/DvCrRWlvfW+0opf1dR3vzlcDg8X9fU1CAxMVGoXZH70DVgwhDsnGXfvn1oamqCSqVCbGwsli9fLtSuyD3oXEU4gjXL66+/LtTQZBg0AyYcujZMQShVhEXrLApCqSIsShaFoFQRHiWLQlCqCI+SRQEoVYKDkkUBKFWCg5JF5ihVgoeSReYoVYKHkkXGKFWCi5JFxihVgouSRaYoVYKPkkWmKFWCj5JFhihVxEHJIkOUKuKgZJEZShXxULLIDKWKeChZZIRSRVyULDJCqSIuShaZoFQRHyWLTFCqiI+SRQYoVaSBkkUGKFWkgZJF4ihVpIOSReIoVaSDkkXCKFWkhZJFwg4fPkypIiGULBLlcrlQVFREqSIhlCwSdfz4cTQ0NFCqSAgliwT1n6ukpqZSqkgIJYsEdXd3IzMzE88//zylioRQs0hQZGQktmzZIqu3aYUCOgwjxEfULIT4iJqFEB/xOmf5/PPPcfDgQXz33XfYvHkzxo0b5/lZVVUVzpw5g7CwMPzmN7/BpEmTeBdLiJh4JUtiYiLefPNNpKamDth+9epVVFdXY9u2bVi/fj3ef/99uN1uXoUSIjZezTJ27Fiv7w2vra1FdnY2wsPDMWbMGMTHx6OxsZHPrggRnSDnLHa7HQaDwfO9Xq+H3W4XYleEBM0Dz1kKCwtx69atQduXLFmCrKwsr5/hOM7nAiwWCywWCwDAbDbDaDT6/NlAaDQawffBCtUqjEBrfWCzvPPOO34PajAYYLPZPN/b7Xbo9Xqvv2symWAymTzfR0RE+L0/fwVjH6xQrcIIpFZBDsMyMzNRXV2N3t5etLW14fr16xg/frwQu/LbW2+9JXYJPqNahRForbymjmtqarBnzx50dnbCbDbjkUcewfr165GYmIinnnoK+fn5CAsLwyuvvIKwMFrSIfLGq1mmTJmCKVOmeP3Z/PnzMX/+fD7DEyIpIffP/b3nR1JHtQoj0FpVnD9TV4SEsJBLFkICFZL3sxw4cACnT5+GVqsFACxduhQZGRkiV3XXF198gb1798LtdmP69OmYO3eu2CUNaeXKlRg5ciTCwsKgVqthNpvFLsmjvLwcVqsVOp0OW7duBQB0dXWhtLQUN2/eRGxsLNasWYPIyEjfBuRC0P79+7mjR4+KXYZXLpeLW7VqFXfjxg2ut7eXe/PNN7krV66IXdaQ8vLyuI6ODrHL8OrixYvct99+y+Xn53u2VVRUcFVVVRzHcVxVVRVXUVHh83h0GCYxjY2NiI+PR1xcHDQaDbKzs1FbWyt2WbKUlpY2KDVqa2uRm5sLAMjNzfXrbxuSh2EAcPLkSZw7dw4pKSlYtmyZ71EssPuvqzMYDPjmm29ErOjBioqKAAAzZsyQ/KxYR0cHYmJiAAAxMTHo7Oz0+bOKbZbhrmmbOXMmFi5cCADYv38/PvjgA+Tl5QW7RK84L5OTKpVKhEp8U1hYCL1ej46ODmzatAkJCQlIS0sTuyxBKLZZfL2mbfr06SguLha4Gt/df12dzWbz/EsoRf3X/Ol0OmRlZaGxsVHSzaLT6eBwOBATEwOHw+GZ5PFFSJ6zOBwOz9c1NTVITEwUsZqBxo0bh+vXr6OtrQ19fX2orq5GZmam2GV5defOHXR3d3u+vnDhApKSkkSuaniZmZk4e/YsAODs2bNDXjnvTUguSu7YsQNNTU1QqVSIjY3F8uXLJfWvt9VqxV//+le43W5MmzZNspcNtba24r333gPww4MBc3JyJFVrWVkZ6uvr4XQ6odPpsGjRImRlZaG0tBTt7e0wGo3Iz8/3+Xw1JJuFkECE5GEYIYGgZiHER9QshPiImoUQH1GzEOIjahZCfETNQoiPqFkI8dH/AbF9HgaL3TPvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# this SVM is a linear hard margin classifier which means it can only deal with linearly separable data\n",
    "# it can only deal with two features \n",
    "class SVM:\n",
    "    def __init__(self):\n",
    "        # initialize with some random values\n",
    "        self.b = 0\n",
    "        self.w = [0,0]\n",
    "        # for visualization later\n",
    "        self.colors = {1:'r',-1:'b'}\n",
    "        self.fig = plt.figure()\n",
    "        self.ax = self.fig.add_subplot(1,1,1)\n",
    "        \n",
    "    def fit(self, x, y, error = 0):\n",
    "        self.xs = x # training data features\n",
    "        self.ys = y # training data classes\n",
    "        \n",
    "        # split data into classes\n",
    "        x_pluses = []\n",
    "        x_minuses = []\n",
    "        for xi, yi in zip(self.xs, self.ys):\n",
    "            if yi==1: x_pluses.append(xi)\n",
    "            elif yi == -1: x_minuses.append(xi)\n",
    "        # convert to numpy arrays for doing calculations\n",
    "        x_pluses = np.array(x_pluses)\n",
    "        x_minuses = np.array(x_minuses)\n",
    "        \n",
    "        # initialize alphas to all have the same values\n",
    "        alpha_pluses = [1 / len(x_pluses) for _ in x_pluses]\n",
    "        alpha_minuses = [1 / len(x_minuses) for _ in x_minuses]\n",
    "        \n",
    "        # initialize h^- and h^+\n",
    "        h_minus = np.dot(alpha_minuses, x_minuses)\n",
    "        h_plus = np.dot(alpha_pluses, x_pluses)\n",
    "        \n",
    "        # this function performs our optimization steps. it optimizes two alphas of the same class per step\n",
    "        def optimize(j, k, alphas, xs, h_opt, h_const):\n",
    "          # parameters:\n",
    "            # j, k: indicies of the alphas to be optimized\n",
    "            # alphas: vector with all alphas (=weights=lagrange multipliers)\n",
    "            # xs: features belogining to the same class as the alphas[j] and alphas[k]\n",
    "            # h_opt: h to be optimized, belongs to same class as xs and alphas[j] and alphas[k]\n",
    "            # h_const: h that stays constant, doesnt belong to same class as xs and alphas[j] and alphas[k]\n",
    "            if j != k: # indicies must be different\n",
    "                c = alphas[j] + alphas[k]\n",
    "                x = xs[j] - xs[k]\n",
    "                h = h_const\n",
    "                u = h_opt - (alphas[j]*xs[j] + alphas[k]*xs[k]) + c*xs[k]\n",
    "                # dont return anything negative or greater than c (otherwise h_opt is outside of its polytope)\n",
    "                a = max(0, min(c, np.dot(x, h-u)/np.dot(x,x))) \n",
    "                return a, c-a\n",
    "            else:\n",
    "                return alphas[j], alphas[k]\n",
    "            \n",
    "        # we need to track the distance between h^+ and h^- so we know when we are done optimizing\n",
    "        last_distance = math.inf\n",
    "        iteration = 0\n",
    "        optimzed = False\n",
    "        while not optimzed:\n",
    "            # optimize all alphas for the plus(+) class\n",
    "            for j in range(len(alpha_pluses)):\n",
    "                for k in range(len(alpha_minuses)):\n",
    "                    alpha_pluses[j], alpha_pluses[k] = optimize(j,k,alpha_pluses, x_pluses, h_plus, h_minus)\n",
    "                    h_plus = np.dot(alpha_pluses, x_pluses)\n",
    "            \n",
    "            # optimize all alphas for the minus(-) class\n",
    "            for j in range(len(alpha_minuses)):\n",
    "                for k in range(len(alpha_pluses)):\n",
    "                    alpha_minuses[j], alpha_minuses[k] = optimize(j,k, alpha_minuses, x_minuses, h_minus, h_plus)\n",
    "                    h_minus = np.dot(alpha_minuses, x_minuses)\n",
    "            \n",
    "            # check if we are done optimizing\n",
    "            new_distance = np.linalg.norm(np.subtract(h_minus, h_plus))\n",
    "            optimzed = last_distance - new_distance == error\n",
    "            last_distance = new_distance\n",
    "            print(f\"iteration {iteration}: h+ & h- distance: \",last_distance)\n",
    "            iteration += 1\n",
    "        \n",
    "        # calculate w & b\n",
    "        self.w = 2*(h_plus-h_minus) / np.linalg.norm(h_plus-h_minus)**2\n",
    "        self.b = -np.dot(self.w,0.5*(h_plus-h_minus)+h_minus)\n",
    "        \n",
    "        # prepare for visualization\n",
    "        self.ax.scatter(h_minus[0], h_minus[1], c = \"cyan\")\n",
    "        self.ax.scatter(h_plus[0], h_plus[1], c = \"pink\")\n",
    "        self.ax.plot([h_plus[0], h_minus[0]], [h_plus[1], h_minus[1]])\n",
    "\n",
    "       \n",
    "\n",
    "    def predict(self, features, visualize = True):\n",
    "        classifications = []\n",
    "        for feature in features:\n",
    "            classifications.append(np.sign(np.dot(np.array(feature), self.w) + self.b))\n",
    "        \n",
    "            if classifications[-1] != 0 and visualize:\n",
    "                self.ax.scatter(feature[0],feature[1],s=200,marker='*', c=self.colors[classifications[-1]])\n",
    "            elif classifications[-1] == 0:\n",
    "                print('featureset',features,'is on the decision boundary')\n",
    "        return classifications\n",
    "    \n",
    "    \n",
    "    \n",
    "    def visualize(self):\n",
    "        self.max_feature_value = 0\n",
    "        self.min_feature_value = math.inf\n",
    "        for feature in self.xs:\n",
    "            for xi in feature:\n",
    "                self.max_feature_value = xi if xi > self.max_feature_value else self.max_feature_value\n",
    "                self.min_feature_value = xi if xi < self.min_feature_value else self.min_feature_value\n",
    "        \n",
    "        [plt.scatter(xi[0], xi[1], color = \"r\" if yi == 1 else \"b\") for xi, yi in zip(self.xs, self.ys)]\n",
    "        \n",
    "        def hyperplane(x,w,b,v):\n",
    "            # v = (w.x+b)\n",
    "            return (-w[0]*x-b+v) / w[1]\n",
    "        \n",
    "        datarange = (self.min_feature_value*0,7,self.max_feature_value*0.9)\n",
    "        print(datarange)\n",
    "        hyp_x_min = datarange[0]\n",
    "        hyp_x_max = datarange[1]\n",
    "        \n",
    "        # w.x + b = 1\n",
    "        # pos sv hyperplane\n",
    "        psv1 = hyperplane(hyp_x_min, self.w, self.b, 1)\n",
    "        psv2 = hyperplane(hyp_x_max, self.w, self.b, 1)\n",
    "        self.ax.plot([hyp_x_min,hyp_x_max], [psv1,psv2], \"k\")\n",
    "        \n",
    "        # w.x + b = -1\n",
    "        # negative sv hyperplane\n",
    "        nsv1 = hyperplane(hyp_x_min, self.w, self.b, -1)\n",
    "        nsv2 = hyperplane(hyp_x_max, self.w, self.b, -1)\n",
    "        self.ax.plot([hyp_x_min,hyp_x_max], [nsv1,nsv2], \"k\")\n",
    "\n",
    "        # w.x + b = 0\n",
    "        # decision\n",
    "        db1 = hyperplane(hyp_x_min, self.w, self.b, 0)\n",
    "        db2 = hyperplane(hyp_x_max, self.w, self.b, 0)\n",
    "        self.ax.plot([hyp_x_min,hyp_x_max], [db1,db2], \"g--\")\n",
    "        \n",
    "        plt.gca().set_aspect('equal', adjustable='box')\n",
    "        plt.show()\n",
    "    \n",
    "    \n",
    "svm = SVM()\n",
    "svm.fit(samples_x, samples_y)\n",
    "\n",
    "predict_samples =  [[0,10],\n",
    "              [1,3],\n",
    "              [3,5],\n",
    "              [5,6],\n",
    "              [6,-5],\n",
    "              [5,3]]\n",
    "svm.predict(predict_samples, False)\n",
    "svm.visualize()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tutorial-env",
   "language": "python",
   "name": "tutorial-env"
  },
  "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}