Recommendation System .ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import *\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SVD++"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradAscent(data,K,alpha,beta,maxCycles):\n",
    "    '''\n",
    "    data: user-item rating matrix\n",
    "    K: the number of features (assume smaller than the no. of user and items)\n",
    "    alpha: learning rate\n",
    "    beta: regularization parameter\n",
    "    maxCycles: no. of iterations\n",
    "    '''\n",
    "    R=mat(data)\n",
    "    m,n=shape(R)\n",
    "    \n",
    "#   initialize user and item latent feature matrix\n",
    "    p=mat(random.random((m,K)))\n",
    "    q=mat(random.random((K,n)))\n",
    "    \n",
    "#   initialize the bisases\n",
    "    bu=zeros(m)\n",
    "    bd=zeros(n)\n",
    "    b=mean(R[where(R!=0)])\n",
    "    \n",
    "#   gradient descent\n",
    "    ls=[]\n",
    "    for step in range(maxCycles):\n",
    "        for i in range(m):\n",
    "            for j in range(n):\n",
    "                if R[i,j] >0:\n",
    "                    r=R[i,j]\n",
    "                    pq=0.0\n",
    "                    for k in range(K):\n",
    "                        pq+=p[i,k]*q[k,j]\n",
    "                    e = r-b-bu[i]-bd[j]-pq\n",
    "#                   Update user and item latent feature matrix\n",
    "                    bu[i] = bu[i] + alpha*(2*e-beta*bu[i])\n",
    "                    bd[j] = bd[j] + alpha*(2*e-beta*bd[j])\n",
    "#                   Update biases\n",
    "                    for k in range(K):\n",
    "                        p[i,k] += alpha*(2*e*q[k,j]-beta*p[i,k])\n",
    "                        q[k,j] += alpha*(2*e*p[i,k]-beta*q[k,j])    \n",
    "\n",
    "#       calculate loss         \n",
    "        loss=0.0\n",
    "        for i in range(m):\n",
    "            for j in range(n):\n",
    "                bb=0.0\n",
    "                dd=0.0\n",
    "                pq2=0.0\n",
    "                pq=0\n",
    "                if R[i,j] > 0:\n",
    "                    r=R[i,j]\n",
    "                    bb +=bu[i]**2\n",
    "                    dd +=bd[j]**2\n",
    "                    for k in range(K):\n",
    "                        pq2 +=p[i,k]**2+q[k,j]**2\n",
    "                        pq +=p[i,k]*q[k,j]\n",
    "                    \n",
    "                    loss=(r-b-bu[i]-bd[j]-pq)**2+beta*(bb+dd+pq2)/2\n",
    "\n",
    "        ls.append(loss)\n",
    "\n",
    "        if loss<0.0001:\n",
    "            break\n",
    "        if step %1000==0:\n",
    "            print(loss)\n",
    "    \n",
    "    return p,q,b,bu,bd,ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# test on a small sample set\n",
    "R = mat([\n",
    "    [5, 3, 0, 1],\n",
    "    [4, 0, 0, 1],\n",
    "    [1, 1, 0, 5],\n",
    "    [1, 0, 0, 4],\n",
    "    [0, 1, 5, 4]])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.856523276163\n",
      "0.0875302625331\n",
      "0.0524067385546\n",
      "0.0530254037105\n",
      "0.052268654279\n",
      "0.0521082435059\n",
      "0.0520507254201\n",
      "0.0519239481276\n",
      "0.0517095163844\n",
      "0.0514249544425\n",
      "[[ 4.98595827  2.99234488  3.62115856  1.0156667 ]\n",
      " [ 3.99325266  2.18754982  3.23347954  1.01275902]\n",
      " [ 1.02102858  0.98361577  5.41321644  4.98486837]\n",
      " [ 1.00130492  0.71038106  4.68801997  3.98955391]\n",
      " [ 1.5303263   1.0352638   4.98624638  3.9956678 ]]\n"
     ]
    }
   ],
   "source": [
    "if __name__ == \"__main__\":\n",
    "    p, q, b, bu, bd, ls = gradAscent(R,3,0.0002,0.02,10000)\n",
    "    \n",
    "    m,n=shape(R)\n",
    "    \n",
    "    bu=list(bu.reshape(m,1))\n",
    "    bu=broadcast_to(bu,(m,n))\n",
    "    \n",
    "    bd=broadcast_to(bd,(m,n))\n",
    "    \n",
    "    result = b+bu+bd+p*q\n",
    "    \n",
    "    print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x11a5e22e8>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHL5JREFUeJzt3X2UXXV97/H3JxmSkAcCMSFAEhIeQgGLgkZMlCIWytOy\noFerwWfl4sMVrReuLRQv11rbZaXLYtfFB8TqrRUQVGhUaLSKWlsDSRQiSUgJkJAJiQkQSEgUEvK9\nf+zfnLNzmJmcmTl79pmzP6+1zprf/u3fOfPdZ4f58Nv77H0UEZiZmQGMKrsAMzNrHw4FMzOrcSiY\nmVmNQ8HMzGocCmZmVuNQMDOzGoeCmZnVOBRs2Eh6q6Rlkp6RtEnSnZJOK7suM6tzKNiwkHQZcC3w\nN8B04Ejg88CFZdaVJ6mr7BqGk6TRZddg7cehYIWTNBn4JPChiPhOROyMiN0R8d2I+FgaM1bStZIe\nS49rJY1N686Q1C3pcklb0izjPWndfEmb83/gJL1B0orUHiXpCkkPSXpC0i2SpqR1cySFpIslPQr8\nOPW/U9L6NP5/S1on6awBvN67JD0q6XFJV+XqGi3pL9Jzd0haLmlWWne8pB9KelLSGklv7uf9nCLp\nq+l92ibp9tT/bkk/bxgbko5N7a9J+oKkOyTtBK4c7HtnncuhYMNhATAOuK2fMVcB84GTgZcCpwIf\nz60/DJgMzAAuBq6TdEhELAF2An+YG/tW4MbU/gjweuA1wBHANuC6ht/9GuAE4BxJJ5LNYN4GHJ77\nnT2aeb3TgN8DzgSulnRC6r8MuAg4HzgIeC+wS9IE4Iep5kPTmM9LenEf79XXgfHAi9P4v+9jXG/e\nCvw1MAn4O4b+3lmniQg//Cj0QfYHdvN+xjwEnJ9bPgdYl9pnAL8FunLrtwDzU/tTwD+m9iSyP3Sz\n0/Jq4Mzc8w4HdgNdwBwggKNz668GbsotjweeA84awOvNzK2/B1iY2muAC3vZ9rcA/97Q9yXg//Qy\n9nBgL3BIL+veDfy8oS+AY1P7a8A/Nawf1HtX9r8pP4p7VOoYqpXmCWCqpK6I2NPHmCOA9bnl9amv\n9hoNz90FTEztG4H/lPRB4L8Bv4yInteaDdwmaW/uuc+TndfosaGhjtpyROyS9ERufTOvt7mPOmeR\nhV+j2cArJT2V6+simxE0mgU8GRHbelnXjA0Ny4N97zYO8vdbm/PhIxsOvwB+R3Yooi+Pkf0R6nFk\n6tuviFhFFiLnse/hD8j+CJ4XEQfnHuMiIv9HLX+r4E3AzJ4FSQcCLxrg6/VlA3BMH/0/bXjNiRHx\nwT7GTpF0cC/rdpLNbHpqP6yXMfvcFrkF7511GIeCFS4iniY7LHOdpNdLGi/pAEnnSfpMGnYT8HFJ\n0yRNTeP/eQC/5kayY+CnA7fm+r8I/LWk2QDp9fv7xNO3gD+W9CpJY4C/BDSE18u7AfgrSXOVeYmk\nFwHfA46T9I70vhwg6RW5cxE1EbEJuJPsnMMhaezpafV9wIslnSxpHPCJJutq1XtnHcChYMMiIj5L\ndqL148BWsv8LvRS4PQ35FLAMWAH8Gvhl6mvWTWTnHn4cEY/n+j8HLAJ+IGkHsAR4ZT91rgQ+DNxM\nNmvYQXb+4tnBvF6DzwK3AD8AtgNfAQ6MiB3A2cBCstnRZuBvgbF9vM47yI7tP5Bq+2iq/b/IPuX1\nb8CDwM/7eH6jlrx31hkU4S/ZMeuLpInAU8DciHik7HrMiuaZglkDSX+cDnFNIPvY5q+BdeVWZTY8\nHApmL3Qh2WGcx4C5ZB8p9ZTaKsGHj8zMrMYzBTMzqxlxF69NnTo15syZU3YZZmYjyvLlyx+PiGn7\nGzfiQmHOnDksW7as7DLMzEYUSev3P8qHj8zMLMehYGZmNQ4FMzOrcSiYmVmNQ8HMzGocCmZmVlO9\nUIiAe++FLVvKrsTMrO1UKxQi4C1vgVNOgaOOgsWLy67IzKytVCsUbr8dbk3fIbJrF1xyCezeXW5N\nZmZtpFqh8M1v7ru8YQPceWc5tZiZtaHqhMLevb0fLvrud4e/FjOzNlWdUFi3Dp566oX9P/zhsJdi\nZtauqhMK991Xb59+Ohx4YNZevx42bSqnJjOzNlOdUFixot6eNy979FiyZPjrMTNrQ9UJhYcfrreP\nPx7mz68vL106/PWYmbWhEfd9CoO2bl29PWcOTJhQX77//uGuxsysLVUnFNbnvl9i9myYPr2+7FAw\nMwOqEgp79kB3d335yCNBgtGj4fnn4ZFHYOfOfWcPZmYVVI1zChs3Zn/8IZshjBsHY8fCccfVx6xa\nVU5tZmZtpBqh0HjoqMeLX1xvr1w5fPWYmbWpaoRCVxeccUZ2gvnYY+v9+ZnCQw8Nd1VmZm2nGucU\nXvUquOuurB1R788HxNq1w1uTmVkbqsZMIU+qt/Oh4JmCmVkFQyHvmGPqbc8UzMwqHgqHH16/B9K2\nbfDkk+XWY2ZWsmqHguRDSGZmOdUOBfAhJDOzHIdCPhTyN80zM6sgh8KcOfX2o4+WVoaZWTsoNBQk\nnStpjaS1kq7oZf2Rku6S9CtJKySdX2Q9vTryyHo7f+WzmVkFFRYKkkYD1wHnAScCF0k6sWHYx4Fb\nIuIUYCHw+aLq6VP+theeKZhZxRU5UzgVWBsRD0fEc8DNwIUNYwI4KLUnA48VWE/vGmcK+Suezcwq\npshQmAFsyC13p768TwBvl9QN3AF8uLcXkvQ+ScskLdu6dWtrqzz4YJg4MWvv2uVrFcys0ooMBfXS\n1/i/4RcBX4uImcD5wNclvaCmiLg+IuZFxLxp06a1uErtewjJ5xXMrMKKDIVuYFZueSYvPDx0MXAL\nQET8AhgHTC2wpt7lDyH5vIKZVViRobAUmCvpKEljyE4kL2oY8yhwJoCkE8hCocXHh5rgmYKZGVBg\nKETEHuBSYDGwmuxTRislfVLSBWnY5cAlku4DbgLeHVHCmV7PFMzMgIK/TyEi7iA7gZzvuzrXXgW8\nusgamuJrFczMAF/RnPG1CmZmgEMhMyt3PnzDhr7HmZl1OIcCZN+r0PONbL/5DezeXW49ZmYlcSgA\njBkDhx6atSNg06Zy6zEzK4lDoceM3MXWGzeWV4eZWYkcCj0cCmZmDoUah4KZmUOhxqFgZuZQqHEo\nmJk5FGocCmZmDoUah4KZmUOhZubMenvjRn8Dm5lVkkOhx0EHwYQJWfu3v4Wnniq3HjOzEjgUekg+\nhGRmledQyHMomFnFORTyHApmVnEOhbx8KHR3l1eHmVlJHAp5nimYWcU5FPIcCmZWcQ6FPIeCmVWc\nQyEvHwqPPVZeHWZmJXEo5B12GIxKb8mWLfDcc+XWY2Y2zBwKeV1dMH16fXnz5vJqMTMrgUOh0RFH\n1Ns+r2BmFeNQaOTzCmZWYQ6FRp4pmFmFORQa5UPBMwUzqxiHQiNfq2BmFeZQaOSZgplVmEOhkWcK\nZlZhDoVGnimYWYU5FBpNmQJjx2btHTuyh5lZRTgUGkmeLZhZZTkUeuML2MysohwKvfEFbGZWUQ6F\n3nimYGYV5VDojWcKZlZRhYaCpHMlrZG0VtIVfYx5s6RVklZKurHIeprmmYKZVVRXUS8saTRwHfBH\nQDewVNKiiFiVGzMXuBJ4dURsk3RoUfUMiGcKZlZRRc4UTgXWRsTDEfEccDNwYcOYS4DrImIbQERs\nKbCe5nmmYGYVVWQozAA25Ja7U1/eccBxkv5D0hJJ5xZYT/Mar1PYu7e8WszMhlGRoaBe+qJhuQuY\nC5wBXATcIOngF7yQ9D5JyyQt27p1a8sLfYHx4+HgVMbu3fDEE8X/TjOzNlBkKHQDs3LLM4HGYzHd\nwL9ExO6IeARYQxYS+4iI6yNiXkTMmzZtWmEF78PnFcysgooMhaXAXElHSRoDLAQWNYy5HXgtgKSp\nZIeTHi6wpub5vIKZVVBhoRARe4BLgcXAauCWiFgp6ZOSLkjDFgNPSFoF3AV8LCLa41iNZwpmVkGF\nfSQVICLuAO5o6Ls61w7gsvRoL74pnplVkK9o7osPH5lZBTkU+uLDR2ZWQQ6FvnimYGYV5FDoi2cK\nZlZBDoW+HHZY9i1sAFu2ZBexmZl1OIdCX7q6YPr0+vKmTeXVYmY2TBwK/fF5BTOrGIdCf3xewcwq\nxqHQH88UzKxiHAr98UzBzCrGodAfzxTMrGIcCv3xTMHMKsah0B/PFMysYhwK/fFMwcwqxqHQnylT\nYOzYrL1jR/YwM+tgDoX+SP5eBTOrFIfC/vi8gplVSFOhIOlPJR2kzFck/VLS2UUX1xY8UzCzCml2\npvDeiNgOnA1MA94DfLqwqtqJTzabWYU0GwrpHtKcD3w1Iu7L9XU2Hz4yswppNhSWS/oBWSgsljQJ\n2FtcWW3EMwUzq5CuJsddDJwMPBwRuyRNITuE1Pk8UzCzCml2prAAWBMRT0l6O/Bx4OniymojnimY\nWYU0GwpfAHZJeinwZ8B64J8Kq6qdNH76KKK8WszMCtZsKOyJiAAuBD4XEZ8DJhVXVhuZMAEmT87a\nu3fD44+XW4+ZWYGaDYUdkq4E3gF8X9Jo4IDiymozPq9gZhXRbCi8BXiW7HqFzcAM4JrCqmo3Pq9g\nZhXRVCikIPgGMFnS64DfRUQ1zimAZwpmVhnN3ubizcA9wJ8AbwbulvSmIgtrK54pmFlFNHudwlXA\nKyJiC4CkacC/Ad8qqrC24pmCmVVEs+cURvUEQvLEAJ478nmmYGYV0exM4V8lLQZuSstvAe4opqQ2\n5JmCmVVEU6EQER+T9Ebg1WQ3wrs+Im4rtLJ24pmCmVVEszMFIuLbwLcLrKV9HXZY9i1sEbBlS3YR\n2wHVuUzDzKqj3/MCknZI2t7LY4ek7cNVZOm6umD69Pry5s3l1WJmVqB+ZwoRUY1bWTRjxox6GGzc\nCLNmlVuPmVkBqvMJoqHy13KaWQU4FJrlk81mVgGFhoKkcyWtkbRW0hX9jHuTpJA0r8h6hsQfSzWz\nCigsFNKdVK8DzgNOBC6SdGIv4yYBHwHuLqqWlvBMwcwqoMiZwqnA2oh4OCKeA24m+z6GRn8FfAb4\nXYG1DJ1nCmZWAUWGwgxgQ265O/XVSDoFmBUR3+vvhSS9T9IyScu2bt3a+kqb4ZmCmVVAkaGgXvpq\n32UpaRTw98Dl+3uhiLg+IuZFxLxp06a1sMQB8EzBzCqgyFDoBvIf5p8J5P+aTgJ+H/iJpHXAfGBR\n255snjIFxo7N2tu3w44d5dZjZlaAIkNhKTBX0lGSxgALgUU9KyPi6YiYGhFzImIOsAS4ICKWFVjT\n4En7XrD26KPl1WJmVpDCQiEi9gCXAouB1cAtEbFS0iclXVDU7y3U7Nn19vr15dVhZlaQpm+INxgR\ncQcNt9iOiKv7GHtGkbW0hEPBzDqcr2geiHworFtXWhlmZkVxKAyEZwpm1uEcCgPhUDCzDudQGAiH\ngpl1OIfCQMycCaPSW7ZpEzz7bLn1mJm1mENhIA44YN8rmzds6HusmdkI5FAYKB9CMrMO5lAYKH8s\n1cw6mENhoDxTMLMO5lAYKIeCmXUwh8JAORTMrIM5FAbKoWBmHcyhMFBHHllvd3fDnj3l1WJm1mIO\nhYEaPx4OPTRr79njr+Y0s47iUBiMY4+tt9euLa8OM7MWcygMxjHH1NsPPVReHWZmLeZQGAyHgpl1\nKIfCYPjwkZl1KIfCYHimYGYdyqEwGI2hEFFeLWZmLeRQGIypU+Ggg7L2M8/Ali3l1mNm1iIOhcGQ\nfAjJzDqSQ2Gw8iebHQpm1iEcCoOVnyn4E0hm1iEcCoPlw0dm1oEcCoPlw0dm1oEcCoOVD4U1a/yx\nVDPrCA6FwZoxAyZOzNrbtsHWreXWY2bWAg6FwZLg+OPry6tXl1eLmVmLOBSG4oQT6m2Hgpl1AIfC\nUDgUzKzDOBSGwqFgZh3GoTAUPqdgZh3GoTAUxxwDXV1Zu7sbduwotx4zsyFyKAzFAQfA3Ln15Qce\nKK8WM7MWcCgMVf68gkPBzEY4h8JQ5UNh5cry6jAza4FCQ0HSuZLWSFor6Ype1l8maZWkFZJ+JGl2\nkfUU4iUvqbfvu6+8OszMWqCwUJA0GrgOOA84EbhI0okNw34FzIuIlwDfAj5TVD2FeelL6+177y2v\nDjOzFihypnAqsDYiHo6I54CbgQvzAyLirojYlRaXADMLrKcYxx4L48dn7c2b4Te/KbceM7MhKDIU\nZgAbcsvdqa8vFwN39rZC0vskLZO0bGu73Xhu9Gg46aT6sg8hmdkIVmQoqJe+Xu8vLentwDzgmt7W\nR8T1ETEvIuZNmzathSW2yMkn19s+hGRmI1hXga/dDczKLc8EHmscJOks4CrgNRHxbIH1FCcfCp4p\nmNkIVuRMYSkwV9JRksYAC4FF+QGSTgG+BFwQEVsKrKVYPtlsZh2isFCIiD3ApcBiYDVwS0SslPRJ\nSRekYdcAE4FbJd0raVEfL9feTjop+34FyC5g27mz3HrMzAapyMNHRMQdwB0NfVfn2mcV+fuHzcSJ\n2UVsq1bB3r2wfDmcfnrZVZmZDZivaG6V+fPr7SVLyqvDzGwIHAqtsmBBvf2LX5RXh5nZEDgUWqVx\nphC9fvrWzKytORRa5YQTYNKkrL15Mzz6aLn1mJkNgkOhVUaPhle+sr7sQ0hmNgI5FFopf17hpz8t\nrw4zs0FyKLTSa19bb//oR+XVYWY2SA6FVlqwAMaNy9oPPggbNvQ/3syszTgUWmncODjttPqyZwtm\nNsI4FFrtzDPrbYeCmY0wDoVWy4fC4sXw/PPl1WJmNkAOhVZ7+cth+vSsvXWrP5pqZiOKQ6HVRo2C\nC3PfOnrbbeXVYmY2QA6FIrzhDfX27bf7lhdmNmI4FIrw2tfWb3nx8MOwbFm59ZiZNcmhUISxY/ed\nLXz1q+XVYmY2AA6Forz3vfX2TTfBb39bXi1mZk1yKBTl9NPh6KOz9lNPwa23lluPmVkTHApFkeDi\ni+vL11zjE85m1vYcCkX6wAdgwoSsff/98P3vl1uPmdl+OBSKNGUKvP/99eWrrvIVzmbW1hwKRbv8\nchg/PmuvWAFf/nK59ZiZ9cOhULQjjoArr6wv//mfZ9cumJm1IYfCcLj8cjj22Ky9fTssXAg7d5Zb\nk5lZLxwKw+HAA+HGG6GrK1teuhTe+EYHg5m1HYfCcHnFK+Daa+vLixdnX8izcmV5NZmZNegqu4BK\n+dCH4PHH4ROfyJbvvRdOPhne+c7smoZTT63PJpqxdy888ww8/fTAHs88kz33+eezn721I7JrLUaN\nqj/yy82u62lL5balF75/zfYNZGzVXjP/s7e+Zn8O5bll//7hrv2cc+CggyiKYoRdUDVv3rxYNtJv\nMPflL2cBsXv3vv2TJsEJJ8CRR2bXN4wdC88+C7t2ZY/GP+7bt/uCOLOqWb0ajj9+wE+TtDwi5u1v\nnGcKZbjkEpg/Hy69FH72s3r/jh1wzz3Zw8ysBA6Fspx0EvzkJ7BkCdxwQ3aOYePGgb/O+PEwefLA\nHpMmZYepRo2C0aPrh3rybSmbhUTUDyvt3bvvcn/r8n09/WW2GzXbN5CxVXvN/M/e+pr9OZTnlv37\ny3huz235C+JQKJMECxZkjwh47DFYuxY2bcruqvq738G4cdkf/vHjs+OIkyfv+/OAA8reCjPrIA6F\ndiHBjBnZw8ysJP5IqpmZ1TgUzMysxqFgZmY1DgUzM6txKJiZWY1DwczMahwKZmZWM+LufSRpK7B+\nkE+fCjzewnJGAm9zNXibq2Eo2zw7Iqbtb9CIC4WhkLSsmRtCdRJvczV4m6thOLbZh4/MzKzGoWBm\nZjVVC4Xryy6gBN7mavA2V0Ph21ypcwpmZta/qs0UzMysHw4FMzOrqUwoSDpX0hpJayVdUXY9gyVp\nlqS7JK2WtFLSn6b+KZJ+KOnB9POQ1C9J/5C2e4Wkl+Ve611p/IOS3lXWNjVL0mhJv5L0vbR8lKS7\nU/3flDQm9Y9Ny2vT+jm517gy9a+RdE45W9IcSQdL+pakB9L+XtDp+1nS/0z/ru+XdJOkcZ22nyX9\no6Qtku7P9bVsv0p6uaRfp+f8gyQNqMCI6PgHMBp4CDgaGAPcB5xYdl2D3JbDgZel9iTgv4ATgc8A\nV6T+K4C/Te3zgTsBAfOBu1P/FODh9POQ1D6k7O3bz7ZfBtwIfC8t3wIsTO0vAh9M7f8BfDG1FwLf\nTO0T074fCxyV/k2MLnu7+tne/wf899QeAxzcyfsZmAE8AhyY27/v7rT9DJwOvAy4P9fXsv0K3AMs\nSM+5EzhvQPWV/QYN005YACzOLV8JXFl2XS3atn8B/ghYAxye+g4H1qT2l4CLcuPXpPUXAV/K9e8z\nrt0ewEzgR8AfAt9L/+AfB7oa9zGwGFiQ2l1pnBr3e35cuz2Ag9IfSDX0d+x+TqGwIf2h60r7+ZxO\n3M/AnIZQaMl+TeseyPXvM66ZR1UOH/X8Y+vRnfpGtDRdPgW4G5geEZsA0s9D07C+tn2kvSfXAn8G\n7E3LLwKeiog9aTlff23b0vqn0/iRtM1HA1uBr6ZDZjdImkAH7+eI2Aj8HfAosIlsvy2ns/dzj1bt\n1xmp3djftKqEQm/H1Eb0Z3ElTQS+DXw0Irb3N7SXvuinv+1Ieh2wJSKW57t7GRr7WTditpns/3xf\nBnwhIk4BdpIdVujLiN/mdBz9QrJDPkcAE4DzehnaSft5fwa6jUPe9qqEQjcwK7c8E3ispFqGTNIB\nZIHwjYj4Tur+jaTD0/rDgS2pv69tH0nvyauBCyStA24mO4R0LXCwpK40Jl9/bdvS+snAk4ysbe4G\nuiPi7rT8LbKQ6OT9fBbwSERsjYjdwHeAV9HZ+7lHq/Zrd2o39jetKqGwFJibPsUwhuyk1KKSaxqU\n9EmCrwCrI+KzuVWLgJ5PILyL7FxDT/8706cY5gNPp+npYuBsSYek/0M7O/W1nYi4MiJmRsQcsn33\n44h4G3AX8KY0rHGbe96LN6XxkfoXpk+tHAXMJTsp13YiYjOwQdLvpa4zgVV08H4mO2w0X9L49O+8\nZ5s7dj/ntGS/pnU7JM1P7+E7c6/VnLJPuAzjiZ3zyT6p8xBwVdn1DGE7TiObDq4A7k2P88mOpf4I\neDD9nJLGC7gubfevgXm513ovsDY93lP2tjW5/WdQ//TR0WT/sa8FbgXGpv5xaXltWn907vlXpfdi\nDQP8VEYJ23oysCzt69vJPmXS0fsZ+EvgAeB+4OtknyDqqP0M3ER2zmQ32f/ZX9zK/QrMS+/fQ8D/\npeHDCvt7+DYXZmZWU5XDR2Zm1gSHgpmZ1TgUzMysxqFgZmY1DgUzM6txKFjlSfrP9HOOpLe2+LX/\nopWvZ1Y0fyTVLJF0BvC/IuJ1A3jO6Ih4vp/1z0TExFbUZzYcPFOwypP0TGp+GvgDSfem+/qPlnSN\npKXpXvbvT+PPUPadFjeSXVCEpNslLU/fBfC+1Pdp4MD0et9IfZel7wq4X9JHU98cZd+X8OX0/B9I\nOnC43wcz8EzBrPZ/840zhfTH/dCI+JSkscB/AH8CzAa+D/x+RDySxk6JiCfTH/OlwGsi4on8TEHS\ny4Gvkd0XX2R3t307sI3sqtR5EXGvpFuARRHxz8P1Hpj16Nr/ELPKOht4iaSe++5MJruPznPAPT2B\nkHxE0htSe1Ya90TD650G3BYROwEkfQf4A7L72zwSEfemccvJ7rdvNuwcCmZ9E/DhiNjnBnJpRrGz\nYfkssi9y2SXpJ2T35ent9frybK79PODDR1YKn1Mwq9tB9hWnPRYDH0y3KkfScemLbhpNBralQDie\n7PBQj909zwd+Brw+3QV0AvAG4N9bvhVmQ+CZglndCmCPpPvIjv1/juwwzi/TbYi3Aq/v5Xn/CnxA\n0gqyu3Iuya27Hlgh6ZcR8TZJX6N+G+cbIuJXyn3hvFnZfKLZzMxqfPjIzMxqHApmZlbjUDAzsxqH\ngpmZ1TgUzMysxqFgZmY1DgUzM6v5/8K0g/wwfI0/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111572a58>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Check Convergence \n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "plt.plot(range(len(ls)), ls, color='r',linewidth=3)\n",
    "plt.title('Convergence curve')\n",
    "plt.xlabel('iteraton')\n",
    "plt.ylabel('loss')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Moive Recommder System based on SVD++"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This data set contains 100,000 ratings (1-5) from 943 users on 1682 movies. The data was collected through the MovieLens web site (movielens.umn.edu) during the seven-month period from September 19th, 1997 through April 22nd, 1998. Source URL: https://grouplens.org/datasets/movielens/100k/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from surprise import SVDpp\n",
    "from surprise import Dataset\n",
    "from surprise.model_selection import cross_validate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dataset ml-100k could not be found. Do you want to download it? [Y/n] y\n",
      "Trying to download dataset from http://files.grouplens.org/datasets/movielens/ml-100k.zip...\n",
      "Done! Dataset ml-100k has been saved to /Users/chenxian/.surprise_data/ml-100k\n"
     ]
    }
   ],
   "source": [
    "data = Dataset.load_builtin('ml-100k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "algo = SVDpp()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evaluating RMSE, MAE of algorithm SVDpp on 5 split(s).\n",
      "\n",
      "                  Fold 1  Fold 2  Fold 3  Fold 4  Fold 5  Mean    Std     \n",
      "RMSE (testset)    0.9238  0.9154  0.9257  0.9207  0.9130  0.9197  0.0048  \n",
      "MAE (testset)     0.7268  0.7222  0.7258  0.7200  0.7179  0.7225  0.0033  \n",
      "Fit time          239.69  232.94  224.31  231.68  258.98  237.52  11.79   \n",
      "Test time         3.85    3.88    3.46    3.69    3.79    3.73    0.15    \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'fit_time': (239.69034695625305,\n",
       "  232.93804383277893,\n",
       "  224.30923509597778,\n",
       "  231.67732906341553,\n",
       "  258.9827971458435),\n",
       " 'test_mae': array([ 0.72678327,  0.72218598,  0.72575526,  0.72003387,  0.71792656]),\n",
       " 'test_rmse': array([ 0.92379001,  0.91542607,  0.92566785,  0.92073923,  0.9130041 ]),\n",
       " 'test_time': (3.846879005432129,\n",
       "  3.8829410076141357,\n",
       "  3.461089849472046,\n",
       "  3.692612886428833,\n",
       "  3.7882609367370605)}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cross_validate(algo,data,measures=['RMSE','MAE'], cv = 5, verbose = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from surprise import accuracy\n",
    "from surprise.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE: 0.9234\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.92335218625386495"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#split traning and testing set (25%)\n",
    "trainset, testset = train_test_split(data, test_size = 0.25)\n",
    "\n",
    "#train SVD++ on the training set and predict ratings for the testing set\n",
    "predictions = algo.fit(trainset).test(testset)\n",
    "\n",
    "#compute RMSE\n",
    "accuracy.rmse(predictions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradAscent2(data, K,alpha,beta,maxCycles):\n",
    "    '''\n",
    "    data: user-item rating matrix\n",
    "    K: the number of latent dimensions \n",
    "    alpha: learning rate\n",
    "    beta: regularization parameter\n",
    "    maxCycles: no. of iterations\n",
    "    '''\n",
    "    dataMat = mat(data)\n",
    "    m, n = shape(dataMat)\n",
    "    p = mat(random.random((m, K)))\n",
    "    q = mat(random.random((K, n)))\n",
    "\n",
    "    ls=[]\n",
    "    for step in range(maxCycles):\n",
    "        for i in range(m):\n",
    "            for j in range(n):\n",
    "                if dataMat[i,j] > 0:\n",
    "                    error = dataMat[i,j]\n",
    "                    for k in range(K):\n",
    "                        error = error - p[i,k]*q[k,j]\n",
    "                    for k in range(K):\n",
    "                        p[i,k] = p[i,k] + alpha * (2 * error * q[k,j] - beta * p[i,k])\n",
    "                        q[k,j] = q[k,j] + alpha * (2 * error * p[i,k] - beta * q[k,j])\n",
    "\n",
    "        loss = 0.0\n",
    "        for i in range(m):\n",
    "            for j in range(n):\n",
    "                if dataMat[i,j] > 0:\n",
    "                    error = 0.0\n",
    "                    for k in range(K):\n",
    "                        error = error + p[i,k]*q[k,j]\n",
    "                    loss = (dataMat[i,j] - error) * (dataMat[i,j] - error)\n",
    "                    for k in range(K):\n",
    "                        loss = loss + beta * (p[i,k] * p[i,k] + q[k,j] * q[k,j]) / 2\n",
    "        \n",
    "        ls.append(loss)\n",
    "        \n",
    "        if loss < 0.001:\n",
    "            break\n",
    "        \n",
    "        if step % 1000 == 0:\n",
    "            print (loss)\n",
    "\n",
    "    return p, q, ls\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.2481523434\n",
      "0.113763832946\n",
      "0.10484625791\n",
      "0.101993644909\n",
      "0.100610463125\n",
      "0.0998412596373\n",
      "0.0993424134791\n",
      "0.0989453915735\n",
      "0.0985898714034\n",
      "0.098256509064\n",
      "[[ 4.97781815  2.98142358  4.61777884  1.00281913]\n",
      " [ 3.98185777  2.4891703   3.58797606  1.00091279]\n",
      " [ 1.00505184  0.98795238  3.82224815  4.97036883]\n",
      " [ 0.9976844   0.66753871  3.82432313  3.98328448]\n",
      " [ 1.8845386   1.01259227  4.9836503   3.98857821]]\n"
     ]
    }
   ],
   "source": [
    "if __name__ == \"__main__\":\n",
    "    '''\n",
    "    test the model on a small sample set\n",
    "    '''\n",
    "    p, q, ls= gradAscent2(R,3,alpha = 0.0002,beta = 0.02,10000) #alpha: learning rate, beta: regularization parameter\n",
    "    result = p*q\n",
    "    print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x11b20f898>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGjlJREFUeJzt3X2UXVWd5vHvkxdC3iAJKTBAIGBnEOxWxFJABRlRFLoV\ne8ZWVBCUXlnLmW6l7dZFBrodZ7odu8el0GsQzYBit4KtiIoMCKjYSKtoJbwFAoYAgfBiCoEQkjRJ\nzG/+OLuqTt3cW3WT1Lmn6u7ns9Zdtc/LvWefe+A+2XufF0UEZmaWr0l1V8DMzOrlIDAzy5yDwMws\ncw4CM7PMOQjMzDLnIDAzy5yDwMwscw4Cq5Sk90nqk/SCpCcl3SDpDXXXy8yGOAisMpI+BlwEfBo4\nADgE+AJwep31KpM0pe46dJKkyXXXwcYfB4FVQtK+wP8A/mtEXBMRmyJiW0R8PyI+ntaZJukiSU+k\n10WSpqVlJ0laJ+kvJa1PrYkPpmXHSXqq/KMm6Y8l3Z3KkySdL2mNpN9K+qakeWnZIkkh6VxJjwI/\nTvM/IGltWv+vJT0i6c278HlnS3pU0tOSLijVa7Kk/5beu1HSckkL07KXSbpZ0jOSHpD07hG+z3mS\nvpK+p2clfTfNP0fSbQ3rhqTfS+UrJF0q6XpJm4Clu/vdWfdyEFhVjgf2Br4zwjoXAMcBRwOvBF4L\nXFha/hJgX+Ag4FzgEklzI+IXwCbgTaV13wdcmcofAd4JvBE4EHgWuKRh228EjgTeKukoipbK+4EF\npW0OaOfz3gAcAZwM/I2kI9P8jwHvBU4D9gE+BGyWNBO4OdV5/7TOFyS9vMV39c/ADODlaf3Pt1iv\nmfcBfwfMBj7Lnn931m0iwi+/xvxF8aP61CjrrAFOK02/FXgklU8CtgBTSsvXA8el8t8CX07l2RQ/\nboem6VXAyaX3LQC2AVOARUAAh5eW/w1wVWl6BrAVePMufN7BpeW/BM5I5QeA05vs+3uAnzbM+xLw\nySbrLgB2AHObLDsHuK1hXgC/l8pXAP/UsHy3vru6/5vyq7pXVv2j1lG/BeZLmhIR21uscyCwtjS9\nNs0b/IyG924GZqXylcDPJH0Y+E/AiogY+KxDge9I2lF67+8oxikGPNZQj8HpiNgs6bel5e183lMt\n6rmQIvAaHQocK+m50rwpFP/yb7QQeCYinm2yrB2PNUzv7nf3+G5u38Y5dw1ZVX4O/DtFN0MrT1D8\n8Aw4JM0bVUTcRxEcpzK8awOKH75TI2JO6bV3RJR/yMq33X0SOHhgQtJ0YL9d/LxWHgNe2mL+vzZ8\n5qyI+HCLdedJmtNk2SaKFsxA3V/SZJ1htxgeg+/OuoyDwCoRERsoulwukfROSTMkTZV0qqR/SKtd\nBVwoqUfS/LT+13ZhM1dS9GmfCHyrNP+LwN9JOhQgff5IZypdDbxd0usk7QV8CtAefF7ZZcD/lLRY\nhVdI2g+4DvgPks5K38tUSa8pjS0MiogngRsoxhDmpnVPTIvvAl4u6WhJewP/vc16jdV3Z13AQWCV\niYjPUQyWXgj0U/xr88+A76ZV/hboA+4G7gFWpHntuopiLOHHEfF0af7FwLXATZI2Ar8Ajh2hnvcC\nfw58g6J1sJFiPOLF3fm8Bp8DvgncBDwPXA5Mj4iNwCnAGRStoKeAvwemtficsyj66u9PdTsv1f3X\nFGdn/RBYDdzW4v2NxuS7s+6gCD+YxqxM0izgOWBxRDxcd33MquYWgRkg6e2p+2omxSmW9wCP1Fsr\ns85wEJgVTqfoonkCWExx+qeby5YFdw2ZmWXOLQIzs8xNiAvK5s+fH4sWLaq7GmZmE8ry5cufjoie\n0dabEEGwaNEi+vr66q6GmdmEImnt6Gu5a8jMLHsOAjOzzDkIzMwy5yAwM8ucg8DMLHMOAjOzzE2I\n00d3y+OPw8MPwwsvwGGHwRFH1F0jM7NxqXtbBJdfDiecAKeeCl/blVvcm5nlpXuDYObMofKmTfXV\nw8xsnMsjCF54ob56mJmNc5UFgaQvS1ovaWVp3jxJN0tanf7OrWr7zJo1VHaLwMyspSpbBFcAb2uY\ndz7wo4hYDPwoTVfDXUNmZm2pLAgi4lbgmYbZpwNfTeWvAu+savsOAjOz9nR6jOCAiHgSIP3dv9WK\nkpZI6pPU19/fv+tbchCYmbVl3A4WR8SyiOiNiN6enlFvp70zDxabmbWl00HwG0kLANLf9ZVtyYPF\nZmZt6XQQXAucncpnA9+rbEvuGjIza0uVp49eBfwcOELSOknnAp8B3iJpNfCWNF0NB4GZWVsqu9dQ\nRLy3xaKTq9rmMI1BEAFSRzZtZjaRjNvB4j02eTLsvXdRjoAtW+qtj5nZONW9QQDuHjIza4ODwMws\ncw4CM7PM5RMEvqjMzKyp7g4CX1RmZjaq7g4Cdw2ZmY3KQWBmlrl8gsBjBGZmTXV3EHiMwMxsVN0d\nBO4aMjMblYPAzCxz+QSBxwjMzJrKJwjcIjAza6q7g8CDxWZmo+ruIHCLwMxsVA4CM7PM5RMEHiw2\nM2uqu4PAYwRmZqPq7iBw15CZ2agcBGZmmcsnCDxGYGbWVHcHwYwZQ+UtW2DHjvrqYmY2TnV3EEya\nNDwMNm+ury5mZuNUdwcBeJzAzGwUDgIzs8zlFQQeMDYz20n3B4EvKjMzG1H3B4G7hszMRuQgMDPL\nXC1BIOkvJN0raaWkqyTtXdnGPEZgZjaijgeBpIOAjwC9EfH7wGTgjMo26DECM7MR1dU1NAWYLmkK\nMAN4orItuWvIzGxEHQ+CiHgc+CzwKPAksCEibmpcT9ISSX2S+vr7+3d/gw4CM7MR1dE1NBc4HTgM\nOBCYKenMxvUiYllE9EZEb09Pz+5v0EFgZjaiOrqG3gw8HBH9EbENuAZ4XWVb82CxmdmI6giCR4Hj\nJM2QJOBkYFVlW/NgsZnZiOoYI7gduBpYAdyT6rCssg26RWBmNqIpdWw0Ij4JfLIjG/MYgZnZiLr/\nymJ3DZmZjaj7g8AtAjOzEXV/EJRbBB4jMDPbSfcHgVsEZmYjyisI3CIwM9tJ9weBB4vNzEbU/UEw\nbRpMSru5dSts21ZvfczMxpnuDwLJrQIzsxF0fxCAB4zNzEaQXxB4wNjMbJg8gsBdQ2ZmLeURBO4a\nMjNrKY8g8NXFZmYt5REEbhGYmbWUXxC4RWBmNkweQeDBYjOzlvIIAncNmZm1lEcQeLDYzKylPILA\nLQIzs5byCAK3CMzMWsojCNwiMDNryUFgZpa5PILAXUNmZi3lEQRuEZiZtZRHELhFYGbWUh5B4BaB\nmVlL+QWBWwRmZsPkEQS+15CZWUt5BMGMGUPlzZthx4766mJmNs7kEQSTJ8P06UU5ArZsqbc+Zmbj\nSC1BIGmOpKsl3S9plaTjK9+oB4zNzJqqq0VwMfCDiHgZ8EpgVeVb9ICxmVlTUzq9QUn7ACcC5wBE\nxFZga+Ub9oCxmVlTdbQIDgf6ga9IukPSZZJmNq4kaYmkPkl9/f39e75Vdw2ZmTVVRxBMAY4BLo2I\nVwGbgPMbV4qIZRHRGxG9PT09e75VX11sZtZUHUGwDlgXEben6aspgqFabhGYmTXV8SCIiKeAxyQd\nkWadDNxX+YY9WGxm1lTHB4uTPwe+Lmkv4CHgg5Vv0YPFZmZN1RIEEXEn0NvRjbpryMysqTyuLAYP\nFpuZtZBPELhFYGbWVJ5B4BaBmdmgfILAg8VmZk21FQSSPippHxUul7RC0ilVV25MuWvIzKypdlsE\nH4qI54FTgB6K0z0/U1mtquDBYjOzptoNAqW/pwFfiYi7SvMmBrcIzMyaajcIlku6iSIIbpQ0G5hY\nj/nyYLGZWVPtXlB2LnA08FBEbJY0j05cDTyWPFhsZtZUuy2C44EHIuI5SWcCFwIbqqtWBdw1ZGbW\nVLtBcCmwWdIrgU8Aa4F/qqxWVfBgsZlZU+0GwfaICOB04OKIuBiYXV21KuAWgZlZU+0GwUZJS4Gz\ngP8naTIwtbpqVWCvvWDy5KK8bRtsrf7pmGZmE0G7QfAe4EWK6wmeAg4C/ndltaqC5AFjM7Mm2gqC\n9OP/dWBfSX8E/HtETKwxAvAppGZmTbR7i4l3A78E/gR4N3C7pHdVWbFKuEVgZraTdq8juAB4TUSs\nB5DUA/yQ4nnDE4cHjM3MdtLuGMGkgRBIfrsL7x0/3DVkZraTdlsEP5B0I3BVmn4PcH01VaqQu4bM\nzHbSVhBExMcl/Wfg9RQ3m1sWEd+ptGZVKAfBxo311cPMbBxp++H1EfFt4NsV1qV6s0vXwDkIzMyA\nUYJA0kYgmi0CIiL2qaRWVdmnVF0HgZkZMEoQRMTEuo3EaMotguefr68eZmbjyMQ782dPuEVgZraT\nfIPALQIzMyC3IHDXkJnZTvIKAncNmZntJK8gcIvAzGwneQWBWwRmZjvJNwjcIjAzA2oMAkmTJd0h\n6bqObdRXFpuZ7aTOFsFHgVUd3aJbBGZmO6klCCQdDPwhcFlHNzxtGkxJF1Nv3QovvtjRzZuZjUd1\ntQguAj4B7Gi1gqQlkvok9fX394/NViUPGJuZNeh4EKRnHq+PiOUjrRcRyyKiNyJ6e3p6xq4C7h4y\nMxumjhbB64F3SHoE+AbwJklf69jWPWBsZjZMx4MgIpZGxMERsQg4A/hxRJzZsQq4RWBmNkxe1xGA\nWwRmZg3afkJZFSLiJ8BPOrpRtwjMzIZxi8DMLHP5BYFbBGZmwzgIzMwyl18QuGvIzGyY/ILALQIz\ns2HyCwK3CMzMhskvCMotgg0b6quHmdk4kV8QzJkzVHYQmJllHgTPPltfPczMxon8gmDu3KHyc8/V\nVw8zs3EivyBobBFE1FcXM7NxIL8gmD69eFIZFE8p27Kl3vqYmdUsvyCA4a0Cdw+ZWebyDILyOIEH\njM0sc3kGgVsEZmaD8gwCtwjMzAY5CNwiMLPM5RkEvqjMzGxQnkHgriEzs0F5BoEHi83MBuUZBG4R\nmJkNyjMI3CIwMxuUZxC4RWBmNijPIHCLwMxsUJ5BUG4RPPNMffUwMxsH8gyC+fOHyk8/XV89zMzG\ngTyDYNasoVtRb9kCmzbVWx8zsxrlGQSSWwVmZkmeQQDQ0zNU7u+vrx5mZjVzEICDwMyy1vEgkLRQ\n0i2SVkm6V9JHO10HYHgQuGvIzDI2pYZtbgf+MiJWSJoNLJd0c0Tc19FalMcI3CIws4x1vEUQEU9G\nxIpU3gisAg7qdD3cNWRmVqh1jEDSIuBVwO1Nli2R1Cepr7+KH2p3DZmZATUGgaRZwLeB8yLi+cbl\nEbEsInojoren/KM9Vtw1ZGYG1BQEkqZShMDXI+KaOurgriEzs0IdZw0JuBxYFRGf6/T2B7lryMwM\nqKdF8HrgLOBNku5Mr9M6Xov99x8q/+Y3Hd+8mdl40fHTRyPiNkCd3u5O5s2DvfaCrVthw4bifkMz\nZ9ZdKzOzjsv3ymIJDjxwaPqJJ+qri5lZjfINAnAQmJnhIBgqP/54ffUwM6tR3kFwUOmCZrcIzCxT\neQeBWwRmZg6CQW4RmFmm8g4Cdw2ZmWUeBOUWwbp19dXDzKxGeQfBwoVD5cceg+3b66uLmVlN8g6C\nGTNgwYKi/LvfwaOP1lsfM7Ma5B0EAC996VD5oYfqq4eZWU0cBIcfPlR2EJhZhhwE5SBYs6a+epiZ\n1cRB4K4hM8ucg6DcIli9ur56mJnVxEFw5JFD5VWrfAqpmWXHQTB3Lhx8cFHeuhV+/et662Nm1mEO\nAoBXvGKofM899dXDzKwGDgKAP/iDobKDwMwy4yCA4UGwYkV99TAzq4GDAODYY4fKP/tZcbsJM7NM\nOAiguJbgJS8pyhs2wMqV9dbHzKyDHAQAEpxwwtD0rbfWVxczsw5zEAx44xuHyt//fn31MDPrMAfB\ngHe8Y6h8yy3wzDP11cXMrIMcBAMWLoTXvrYob98OV15Zb33MzDrEQVB21llD5Ysu8u0mzCwLDoKy\nc86BOXOK8po18PnP11odM7NOmFJ3BcaVWbPg/POLF8DSpXDAAUVLQRr5vRHw4ouwaVNxHcKOHcP/\nDpCGPmusy1V97kDZzLqSIqLzG5XeBlwMTAYui4jPjLR+b29v9PX1daRubNsGxx8Py5cPzVu8GF79\napg9G7ZsgeefL643aPy7bVtn6li3TgVPp7bRWN7d76SO99a5be9zZ967337wwx/u5ia1PCJ6R1uv\n4y0CSZOBS4C3AOuAX0m6NiLu63Rdmpo6tTh99KSThu5Eunq1n1VQFlG8zKx6BxxQ+Sbq6Bp6LfBg\nRDwEIOkbwOnA+AgCgAULoK8PLrwQLr+86O5px9SpRffS5MnFa9Kkob/S8B/QKstj/blm1tXqCIKD\ngMdK0+uAY1usW5/Zs+Hii+HTny66idauLQJh+nTYZ5/ite++xWugPG1a9/alVx02nQi0qkKvrvfW\nuW3vc+feO6X6n+k6gqDZL+VO35KkJcASgEMOOaTqOrU2cyaceGJ92x8vPGhs1rXqOH10HbCwNH0w\n8ETjShGxLCJ6I6K3p6enY5UzM8tNHUHwK2CxpMMk7QWcAVxbQz3MzIwauoYiYrukPwNupDh99MsR\ncW+n62FmZoVaLiiLiOuB6+vYtpmZDedbTJiZZc5BYGaWOQeBmVnmarnX0K6S1A+s3c23zweeHsPq\nTATe5zx4n7vfnu7voREx6vn3EyII9oSkvnZuutRNvM958D53v07tr7uGzMwy5yAwM8tcDkGwrO4K\n1MD7nAfvc/fryP52/RiBmZmNLIcWgZmZjcBBYGaWua4OAklvk/SApAclnV93fXaXpIWSbpG0StK9\nkj6a5s+TdLOk1env3DRfkv4x7ffdko4pfdbZaf3Vks6ua5/aJWmypDskXZemD5N0e6r/v6Q72CJp\nWpp+MC1fVPqMpWn+A5LeWs+etEfSHElXS7o/He/ju/04S/qL9N/1SklXSdq7246zpC9LWi9pZWne\nmB1XSa+WdE96zz9Ku/jgkIjoyhfFnU3XAIcDewF3AUfVXa/d3JcFwDGpPBv4NXAU8A/A+Wn++cDf\np/JpwA0UDwE6Drg9zZ8HPJT+zk3luXXv3yj7/jHgSuC6NP1N4IxU/iLw4VT+L8AXU/kM4F9S+ah0\n7KcBh6X/JibXvV8j7O9XgT9N5b2AOd18nCmeWPgwML10fM/ptuMMnAgcA6wszRuz4wr8Ejg+vecG\n4NRdql/dX1CFX/zxwI2l6aXA0rrrNUb79j3gLcADwII0bwHwQCp/CXhvaf0H0vL3Al8qzR+23nh7\nUTy06EfAm4Dr0n/kTwNTGo8xxW3Nj0/lKWk9NR738nrj7QXsk34U1TC/a48zQ4+unZeO23XAW7vx\nOAOLGoJgTI5rWnZ/af6w9dp5dXPXULNnIx9UU13GTGoKvwq4HTggIp4ESH/3T6u12veJ9p1cBHwC\n2JGm9wOei4jtabpc/8F9S8s3pPUn0j4fDvQDX0ndYZdJmkkXH+eIeBz4LPAo8CTFcVtOdx/nAWN1\nXA9K5cb5bevmIGjr2cgTiaRZwLeB8yLi+ZFWbTIvRpg/7kj6I2B9RCwvz26yaoyybMLsM8W/cI8B\nLo2IVwGbKLoMWpnw+5z6xU+n6M45EJgJnNpk1W46zqPZ1X3c433v5iBo69nIE4WkqRQh8PWIuCbN\n/o2kBWn5AmB9mt9q3yfSd/J64B2SHgG+QdE9dBEwR9LAA5XK9R/ct7R8X+AZJtY+rwPWRcTtafpq\nimDo5uP8ZuDhiOiPiG3ANcDr6O7jPGCsjuu6VG6c37ZuDoKueTZyOgPgcmBVRHyutOhaYODMgbMp\nxg4G5n8gnX1wHLAhNT1vBE6RNDf9S+yUNG/ciYilEXFwRCyiOHY/joj3A7cA70qrNe7zwHfxrrR+\npPlnpLNNDgMWUwysjTsR8RTwmKQj0qyTgfvo4uNM0SV0nKQZ6b/zgX3u2uNcMibHNS3bKOm49B1+\noPRZ7al7AKXiwZnTKM6wWQNcUHd99mA/3kDR1LsbuDO9TqPoG/0RsDr9nZfWF3BJ2u97gN7SZ30I\neDC9Plj3vrW5/ycxdNbQ4RT/gz8IfAuYlubvnaYfTMsPL73/gvRdPMAunk1Rw74eDfSlY/1dirND\nuvo4A58C7gdWAv9MceZPVx1n4CqKMZBtFP+CP3csjyvQm76/NcD/oeGEg9FevsWEmVnmurlryMzM\n2uAgMDPLnIPAzCxzDgIzs8w5CMzMMucgMKuYpPMkzShNXy9pTp11Mivz6aNmeyhdxKOI2NFi+SMU\n54I/3dGKmbXJLQLrWpL+WsV9/W9O97n/K0kvlfQDScsl/VTSy9K6V6T7uP9M0kOS3lX6nI9L+lW6\nN/yn0rxFKp4X8AVgBbBQ0qWS+lTcW39gvY9Q3EPnFkm3pHmPSJqfyh9TcR/+lZLOa/js/5s+6yZJ\n0zv53Vlm6r7izi+/qnhRXGl5JzCd4hkOq4G/oriCc3Fa51iKWxQAXEFxxeokinvbP5jmn0LxAHGl\nZddR3Ft+EcVdUY8rbXPgytDJwE+AV6TpR4D5pfUeAeYDr6a4cnQmMAu4l+LOsouA7cDRaf1vAmfW\n/Z361b2vgZs6mXWbNwDfi4gtAJK+T3F7gtcB3yo9wGla6T3fjaJ75z5JB6R5p6TXHWl6FsV9bB4F\n1kbEL0rvf7ekJRR3EV1AESh3j1LH70TEplTHa4ATKO4183BE3JnWW04RDmaVcBBYt2p2a95JFPe5\nP7rFe15s8n4B/ysivjTsw4vnQmwqTR9G0eJ4TUQ8K+kKiuDZ1To2q8vvKFo2ZpXwGIF1q9uAt6t4\n/u0s4A+BzcDDkv4EBp8N+8pRPudG4EPpM5B0kKT9m6y3D0UwbEitifI99TdSdE81uhV4Z7rz5kzg\nj4Gftr+LZmPDLQLrShHxK0nXUjzHdi3FHT03AO8HLpV0ITCV4lkHd43wOTdJOhL4eepOegE4k+Jf\n6eX17pJ0B0U//0PAv5UWLwNukPRkRPzH0ntWpJbDwO2SL4uIO1R6ILtZJ/j0UetakmZFxAvpHP5b\ngSURsaLuepmNN24RWDdbJukoir76rzoEzJpzi8DMLHMeLDYzy5yDwMwscw4CM7PMOQjMzDLnIDAz\ny9z/B1kjdUYXQgiWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11a66a2e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# check convergence\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "\n",
    "plt.plot(range(len(ls)), ls, color='r',linewidth=3)\n",
    "plt.title('Convergence curve')\n",
    "plt.xlabel('generation')\n",
    "plt.ylabel('loss')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}