3-factorization-machine-simple-NN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Factorization Machines (FM)\n",
    "\n",
    "Factorization machines (FM) model all nested interactions up to order $d$ between the $p$ input variables in $x$ using factorized interaction parameters. The factorization machine (FM) model is defined as\n",
    "\n",
    "$$\\hat{y}(\\text{x})= w_0 + \\sum_{j=1}^d w_j x_j + \\sum_{j=1}^d \\sum_{j'=j+1}^d x_j x_{j'} < v_j, v_{j'} > $$\n",
    "\n",
    "The first part corresponds to linear model, the second part contains all pairwise interactions between input variables. The important difference to standard polynomial regression is that the effect of the interaction is not modeled by an independent parameter $w_{j, j'}$ but with a factorized parametrization $w_{j, j'} = <v_j, v_{j'}>=\\sum_{f=1}^k v_{j,f}v_{j',f}$, which corresponds to the assumption that the effect of pairwise interactions has a low rank.\n",
    "\n",
    "Three learning methods (for minimization of given loss function $\\mathcal{L}$) are widely used for FMs:\n",
    "- stochastic gradient descent (SGD) and modifications\n",
    "- alternating least-squares (ALS)\n",
    "- and Markov Chain Monte Carlo (MCMC) inference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy\n",
    "import astropy\n",
    "from astropy.io import ascii\n",
    "from astropy.table import Table\n",
    "from sklearn.metrics import roc_auc_score\n",
    "\n",
    "# fastFM is python-friendly implementation of factorization machines\n",
    "import fastFM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read data for flight delay challenge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = ascii.read(\"data/training.csv\", delimiter=',')  \n",
    "test_kaggle = ascii.read(\"data/test.csv\", delimiter=',')  \n",
    "\n",
    "binary_target = (data['ARRIVAL_DELAY'] > 10) * 1\n",
    "data.remove_column('ARRIVAL_DELAY')\n",
    "\n",
    "# take small part of data\n",
    "small_data = data[::20]\n",
    "small_binary_target = binary_target[::20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### helping functions from the first notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def train_test_split(*arrays, **kargs):\n",
    "    '''modification of sklearn's train_test_split to support astropy. See sklearn documentation for parameters '''\n",
    "    from sklearn.cross_validation import train_test_split\n",
    "    arrays2 = map(lambda x: numpy.array(x) if isinstance(x, Table) else x, arrays)\n",
    "    results = list(train_test_split(*arrays2, **kargs))\n",
    "    \n",
    "    for i in range(len(results) // 2):\n",
    "        if isinstance(arrays[i], Table):\n",
    "            results[2 * i] = Table(results[2 * i])\n",
    "            results[2 * i + 1] = Table(results[2 * i + 1])\n",
    "    return results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import FileLink\n",
    "\n",
    "def create_solution(predictions, filename='flight-delay-predictions.csv'):\n",
    "    result = astropy.table.Table({'ID': numpy.arange(len(predictions)), 'ARRIVAL_DELAY': predictions})\n",
    "    result.write('data/{}'.format(filename), format='csv', delimiter=',', overwrite=True)\n",
    "    return FileLink('data/{}'.format(filename))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "non_categorical_features = ['SCHEDULED_DEPARTURE', 'DISTANCE', 'SCHEDULED_ARRIVAL']\n",
    "categorical_features = list(set(data.columns) - set(non_categorical_features))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/antares/.virtualenvs/rep/lib/python2.7/site-packages/sklearn/cross_validation.py:44: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.\n",
      "  \"This module will be removed in 0.20.\", DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "# prepare training and test samples\n",
    "trainX, testX, trainY, testY = train_test_split(small_data, small_binary_target, random_state=42, train_size=0.5)\n",
    "\n",
    "# normalize non categorical features for optimization process stability\n",
    "for column in non_categorical_features:\n",
    "    coeff = numpy.std(trainX[column])\n",
    "    testX[column]  = testX[column] / coeff\n",
    "    trainX[column] = trainX[column] / coeff\n",
    "    # prepare test kaggle samples\n",
    "    test_kaggle[column] = test_kaggle[column] / coeff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare data for factorization machines\n",
    "one-hot encoding is used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OneHotEncoder(categorical_features=array([ True,  True,  True,  True,  True,  True,  True,  True, False,\n",
       "       False, False], dtype=bool),\n",
       "       dtype=<type 'numpy.float64'>, handle_unknown='ignore',\n",
       "       n_values='auto', sparse=True)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.preprocessing import OneHotEncoder\n",
    "coder = OneHotEncoder(sparse=True, categorical_features=numpy.in1d(trainX.colnames, categorical_features),\n",
    "                      handle_unknown='ignore')\n",
    "coder.fit(trainX.to_pandas())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# transform categorical features into vectors with zeros and ones\n",
    "trainX_sparse = coder.transform(trainX.to_pandas())\n",
    "testX_sparse = coder.transform(testX.to_pandas())\n",
    "\n",
    "test_kaggle_sparse = coder.transform(test_kaggle.to_pandas())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ALS \n",
    "\n",
    "Alternating least-squares (ALS) is a coordinate descent method, frequently used in factorization models.\n",
    "\n",
    "<!--\n",
    "The optimization approach of SGD is based on iterating over cases (rows) of the training data and performing small steps in the direction of a smaller loss. Coordinate descent or alternating least-squares (ALS) takes another approach by minimizing the loss per model parameter: we iterate over parameters, fix all parameters except the current parameter $\\theta^*$, find $\\hat{\\theta}^*=\\arg\\min_{\\theta^*} \\mathcal L$ and use it as estimation of parameter $\\theta^*$.\n",
    "-->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# import model, it follows sklearn interface\n",
    "from fastFM.als import FMClassification as FMClassificationALS\n",
    "\n",
    "fm = FMClassificationALS(n_iter=100)\n",
    "# target should be coded as {-1, 1}\n",
    "fm.fit(trainX_sparse, trainY * 2 - 1)\n",
    "proba_als = fm.predict_proba(testX_sparse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.54313173055361774"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "roc_auc_score(testY, proba_als)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regularization and initialization of parameters (they are initialized with normal distribution) play significant role in FM training:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.60078186815056023"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fm = FMClassificationALS(n_iter=50, l2_reg_V=5, l2_reg_w=5, init_stdev=0.01)\n",
    "# target should be coded as {-1,1}\n",
    "fm.fit(trainX_sparse, trainY * 2 - 1)\n",
    "roc_auc_score(testY, fm.predict_proba(testX_sparse))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Averaging for different regularizations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.555581409726\n",
      "0.575374847872\n",
      "0.584520700631\n",
      "0.591473818172\n",
      "0.595013148803\n",
      "0.59923225186\n",
      "0.603069565927\n",
      "0.605306563354\n",
      "0.606520811103\n",
      "0.607179476674\n"
     ]
    }
   ],
   "source": [
    "proba = numpy.zeros(len(testY))\n",
    "for reg_V in [0.1, 0.5, 1, 2, 3, 4, 5, 6, 7, 8]:\n",
    "    fm = FMClassificationALS(n_iter=50, l2_reg_V=reg_V, l2_reg_w=5, init_stdev=0.01)\n",
    "    # target should be coded as {-1,1}\n",
    "    fm.fit(trainX_sparse, trainY * 2 - 1)\n",
    "    proba += fm.predict_proba(testX_sparse)\n",
    "    print roc_auc_score(testY, proba)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MCMC\n",
    "\n",
    "Both ALS and SGD learn the best parameters which are used for a point estimate of $\\hat{y}$. MCMC is a Bayesian\n",
    "inference technique that generates the distribution of $\\hat{y}$ by sampling. There is a [demo](http://arogozhnikov.github.io/2016/12/19/markov_chain_monte_carlo.html) with playground which explains the MCMC approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from fastFM.mcmc import FMClassification as FMClassificationMCMC\n",
    "fm = FMClassificationMCMC(n_iter=100)\n",
    "# keeping all models takes too much space, so predictions are accumulated in the train-time\n",
    "proba_mcmc = fm.fit_predict_proba(trainX_sparse, trainY * 2 - 1, testX_sparse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.63711272853109091"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "roc_auc_score(testY, proba_mcmc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<a href='data/flight-delay-mcmc.csv' target='_blank'>data/flight-delay-mcmc.csv</a><br>"
      ],
      "text/plain": [
       "/Users/antares/Yandex.Disk.localized/2017-04-Heidelberg/practice/data/flight-delay-mcmc.csv"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "proba_mcmc_kaggle = fm.fit_predict_proba(trainX_sparse, trainY * 2 - 1, test_kaggle_sparse)\n",
    "create_solution(proba_mcmc_kaggle, 'flight-delay-mcmc.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple mixing with GB\n",
    "\n",
    "As soon as FM and GB are very different model, mixing of two models is expected to be better than separate models. Let's try:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.64552997411209456"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from xgboost.sklearn import XGBClassifier\n",
    "\n",
    "xgb_clf = XGBClassifier(n_estimators=100, learning_rate=0.1, max_depth=4)\n",
    "xgb_clf.fit(trainX.to_pandas(), trainY)\n",
    "proba_gb = xgb_clf.predict_proba(testX.to_pandas())[:, 1]\n",
    "roc_auc_score(testY, proba_gb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simple averaging of two models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.65215540498275426"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "roc_auc_score(testY, (proba_mcmc + proba_gb) / 2.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- ## Stacking (stacked generalization)\n",
    "\n",
    "<img src='http://arogozhnikov.github.io/images/etc/rogozhnikov_stacking.png' />\n",
    "-->\n",
    "\n",
    "###  You can do better!\n",
    "- Increase number of samples\n",
    "- Play with regularization in FM model: regularization is the most important in FM training\n",
    "- Feature engineering matters for FMs too!\n",
    "- Mix up different models using also machine learning?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neural Networks & pytorch\n",
    "\n",
    "There are many different libraries for deep learning, we'll use [pytorch](http://pytorch.org), a fastest to dive library that benefits from torch infrastructure and [chainer](http://chainer.org) approach 'define-by-run'.\n",
    "\n",
    "> Warning: `pytorch` is unstable, not even beta.\n",
    "\n",
    "### Autograd 101\n",
    "\n",
    "The main point of NN training is gradient computation with respect to parameters (weights). \n",
    "\n",
    "Consider the expression \n",
    "$$ f = <w, x> = \\sum_j x_j w_j$$\n",
    "and the gradient with respect to $w$\n",
    "$$ \\frac{\\partial{f}}{\\partial{w_i}} = x_i$$\n",
    "or in vector form\n",
    "$$ \\frac{\\partial{f}}{\\partial{w}} = x$$\n",
    "\n",
    "Let see how this is done in `pytorch`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import torch\n",
    "from torch.autograd import Variable\n",
    "dtype = torch.FloatTensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create sample with shape [1, 10] and set that gradients are not needed for this variable\n",
    "x = Variable(torch.randn(1, 10).type(dtype), requires_grad=False)\n",
    "# create weights vector with shape [10, 1] and set that we want to compute gradients with respect to this varibale\n",
    "w = Variable(torch.randn(10, 1).type(dtype) / 2, requires_grad=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "-0.6530 -1.6536 -0.9376  1.0105 -1.0674  1.4680  0.5759  0.4326 -0.0700 -0.4823\n",
      "[torch.FloatTensor of size 1x10]\n",
      "\n",
      "\n",
      " 0.1229\n",
      " 0.4837\n",
      " 0.5096\n",
      "-0.4697\n",
      " 0.2409\n",
      " 1.3351\n",
      " 0.2483\n",
      " 0.8510\n",
      " 0.1855\n",
      " 0.5166\n",
      "[torch.FloatTensor of size 10x1]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# print values of x and w\n",
    "print x.data\n",
    "print w.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.65301806 -1.6536448  -0.93764013  1.01047695 -1.06735981  1.46797395\n",
      "   0.57594097  0.43258277 -0.0699688  -0.4823184 ]] (1, 10)\n"
     ]
    }
   ],
   "source": [
    "# convert to numpy\n",
    "print x.data.numpy(), x.data.numpy().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "-0.6530\n",
       "-1.6536\n",
       "-0.9376\n",
       " 1.0105\n",
       "-1.0674\n",
       " 1.4680\n",
       " 0.5759\n",
       " 0.4326\n",
       "-0.0700\n",
       "-0.4823\n",
       "[torch.FloatTensor of size 10x1]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute expression (as in numpy), where mm is matrix multiplication\n",
    "f = x.mm(w)\n",
    "# Use autograd to compute the backward pass. This call will compute the\n",
    "# gradient of f with respect to all Variables with requires_grad=True.\n",
    "f.backward()\n",
    "# call w.grad will be Variables holding the gradient\n",
    "# of the loss with respect to w.\n",
    "w.grad.data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, gradient is equal to $x$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A bit harder\n",
    "\n",
    "Compute gradients with respect to $w_1$:\n",
    "\n",
    "$$\\mathcal{L} = \\sum_i \\dfrac{\\sin{(e^{(2 w_{1, i} - 1) x_i})}}{3}$$\n",
    "\n",
    "for $x = (1, .., 1)\\in \\mathbb{R}^{10}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = Variable(torch.ones(10).type(dtype), requires_grad=False)\n",
    "w1 = Variable(torch.FloatTensor([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]), requires_grad=True)\n",
    "\n",
    "w2 = 2 * w1 - 1\n",
    "w3 = torch.exp(w2 * x)\n",
    "w4 = torch.sin(w3) / 3\n",
    "loss = w4.sum()\n",
    "loss.backward()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Variable containing:\n",
       "-1.6522e+00\n",
       " 4.4000e+00\n",
       "-7.1832e+01\n",
       "-7.1389e+02\n",
       "-3.2876e+03\n",
       "-3.8672e+03\n",
       "-5.6167e+04\n",
       "-6.0108e+04\n",
       " 1.6065e+07\n",
       " 1.0951e+08\n",
       "[torch.FloatTensor of size 10]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w1.grad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question: how was gradient computed?\n",
    "\n",
    "<details>\n",
    " <summary>Answer</summary>\n",
    " <p>\n",
    " $$ \\mathcal{L} = \\mathcal{L}(w_5) = \\mathcal{L}(w_5 (w_4)) = \\dots = \\mathcal{L}(w_5 (w_4 (w_3 (w_2 (w_1)))))\n",
    " $$\n",
    " </p>\n",
    " <p>\n",
    "     Computations (forward propagation) are done in the following order:\n",
    "     $$\n",
    "         w_1 \\rightarrow w_2 \\rightarrow w_3 \\rightarrow w_4 \\rightarrow w_5\n",
    "     $$\n",
    " </p>\n",
    " <p>\n",
    "     Computation of gradient:\n",
    " </p>\n",
    " <p>\n",
    "     \\begin{aligned}\n",
    "         \\dfrac{ \\partial \\mathcal{L}}{\\partial w_5} & \\text{is straightforward computation} \\\\\n",
    "         \\dfrac{ \\partial \\mathcal{L}}{\\partial w_4} &= \\dfrac{ \\partial \\mathcal{L}}{\\partial w_5} \\dfrac{ \\partial \\mathcal{w_5}}{\\partial w_4} \\\\\n",
    "         \\dfrac{ \\partial \\mathcal{L}}{\\partial w_3} &= \\dfrac{ \\partial \\mathcal{L}}{\\partial w_4} \\dfrac{ \\partial \\mathcal{w_4}}{\\partial w_3} \\\\\n",
    "         \\dfrac{ \\partial \\mathcal{L}}{\\partial w_2} &= \\dfrac{ \\partial \\mathcal{L}}{\\partial w_3} \\dfrac{ \\partial \\mathcal{w_3}}{\\partial w_2} \\\\\n",
    "         \\dfrac{ \\partial \\mathcal{L}}{\\partial w_1} &= \\dfrac{ \\partial \\mathcal{L}}{\\partial w_2} \\dfrac{ \\partial \\mathcal{w_2}}{\\partial w_1} \\\\\n",
    "     \\end{aligned}\n",
    " </p>\n",
    " <p>\n",
    "     These computations are done in reverse order \n",
    "     $$ \\dfrac{ \\partial \\mathcal{L}}{\\partial w_1}  \\leftarrow \n",
    "         \\dfrac{ \\partial \\mathcal{L}}{\\partial w_2} \\leftarrow \n",
    "         \\dfrac{ \\partial \\mathcal{L}}{\\partial w_3} \\leftarrow \n",
    "         \\dfrac{ \\partial \\mathcal{L}}{\\partial w_4} \\leftarrow \n",
    "         \\dfrac{ \\partial \\mathcal{L}}{\\partial w_5}  \n",
    "     $$ \n",
    "     this process is called `backward propagation` or just `backpropagation`.\n",
    " </p>\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "w4 Variable containing:\n",
      " 1\n",
      " 1\n",
      " 1\n",
      " 1\n",
      " 1\n",
      " 1\n",
      " 1\n",
      " 1\n",
      " 1\n",
      " 1\n",
      "[torch.FloatTensor of size 10]\n",
      "\n",
      "w3 Variable containing:\n",
      "-0.3039\n",
      " 0.1095\n",
      "-0.2420\n",
      "-0.3255\n",
      "-0.2029\n",
      "-0.0323\n",
      "-0.0635\n",
      "-0.0092\n",
      " 0.3325\n",
      " 0.3068\n",
      "[torch.FloatTensor of size 10]\n",
      "\n",
      "w2 Variable containing:\n",
      "-8.2612e-01\n",
      " 2.2000e+00\n",
      "-3.5916e+01\n",
      "-3.5694e+02\n",
      "-1.6438e+03\n",
      "-1.9336e+03\n",
      "-2.8084e+04\n",
      "-3.0054e+04\n",
      " 8.0326e+06\n",
      " 5.4753e+07\n",
      "[torch.FloatTensor of size 10]\n",
      "\n",
      "w1 Variable containing:\n",
      "-1.6522e+00\n",
      " 4.4000e+00\n",
      "-7.1832e+01\n",
      "-7.1389e+02\n",
      "-3.2876e+03\n",
      "-3.8672e+03\n",
      "-5.6167e+04\n",
      "-6.0108e+04\n",
      " 1.6065e+07\n",
      " 1.0951e+08\n",
      "[torch.FloatTensor of size 10]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Let see the backpropagation\n",
    "def printgrad(name):\n",
    "    def print_hook(grad):\n",
    "        print name, grad \n",
    "    return print_hook\n",
    "    \n",
    "x = Variable(torch.ones(10).type(dtype), requires_grad=False)\n",
    "w1 = Variable(torch.FloatTensor([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]), requires_grad=True)\n",
    "w1.register_hook(printgrad('w1'))\n",
    "w2 = 2 * w1 - 1\n",
    "w2.register_hook(printgrad('w2'))\n",
    "w3 = torch.exp(w2 * x)\n",
    "w3.register_hook(printgrad('w3'))\n",
    "w4 = torch.sin(w3) / 3\n",
    "w4.register_hook(printgrad('w4'))\n",
    "loss = w4.sum()\n",
    "\n",
    "\n",
    "# first the gradient over w4 is printed, then over w3, then over w2, and then over w1\n",
    "loss.backward()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple NN\n",
    "\n",
    "Let code the simple NN with one hidden layer with `pytorch`\n",
    "\n",
    "Parameters:\n",
    "- $W$, $v$\n",
    "\n",
    "Calculations:\n",
    "- hidden activations: $h_{ik} = \\sigma(\\sum_j X_{ij} W_{jk} )$\n",
    "- $p_i$ = $\\sigma(\\sum_k h_{ik} v_{k} )   $\n",
    "- loss function can be written as\n",
    "  $$\\mathcal{L}=-\\sum_i y_i \\log{p_i} + (1-y_i)\\log{(1 - p_i)}\\qquad,$$ where $y \\in \\{0, 1\\}$\n",
    "- compute loss function gradient with respect to parameters\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare data for NNs with a `StandardScaler`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "scaler = StandardScaler().fit(trainX.to_pandas())\n",
    "trainX_scaled = scaler.transform(trainX.to_pandas())\n",
    "testX_scaled = scaler.transform(testX.to_pandas())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from torch.utils.data import TensorDataset, DataLoader\n",
    "\n",
    "train_dataset = TensorDataset(torch.FloatTensor(trainX_scaled), torch.FloatTensor(trainY.reshape((len(trainY), 1))))\n",
    "loader_train = DataLoader(dataset=train_dataset, batch_size=64, shuffle=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 \t0.706803307415\n",
      "5 \t0.58766642908\n",
      "10 \t0.556104045263\n",
      "15 \t0.541402690094\n",
      "20 \t0.534092096746\n",
      "25 \t0.530062030943\n",
      "30 \t0.527607891954\n",
      "35 \t0.525937170661\n",
      "40 \t0.524670135692\n",
      "45 \t0.523656870161\n",
      "50 \t0.522795454201\n",
      "55 \t0.522036619202\n",
      "60 \t0.52134101868\n",
      "65 \t0.520706337014\n",
      "70 \t0.520127066595\n",
      "75 \t0.519583558773\n",
      "80 \t0.519081023461\n",
      "85 \t0.518599956505\n",
      "90 \t0.518144798859\n",
      "95 \t0.517732902518\n"
     ]
    }
   ],
   "source": [
    "dtype = torch.FloatTensor\n",
    "\n",
    "# batch_size is batch size; dim_input is input dimension;\n",
    "# dim_hidden is hidden dimension; dim_output is output dimension.\n",
    "batch_size, dim_input, dim_hidden, dim_output = 64, len(trainX.colnames), 15, 1\n",
    "\n",
    "# Create random Tensors for weights, and wrap them in Variables.\n",
    "# Setting requires_grad=True indicates that we want to compute gradients with\n",
    "# respect to these Variables during the backward pass.\n",
    "w = Variable(torch.randn(dim_input, dim_hidden).type(dtype), requires_grad=True)\n",
    "v = Variable(torch.randn(dim_hidden, dim_output).type(dtype), requires_grad=True)\n",
    "\n",
    "learning_rate = 1e-4\n",
    "for iteration in range(100):\n",
    "    # use SGD optimization and iterate over batch\n",
    "    average_loss = 0\n",
    "    for nbatch, (batch_X, batch_y) in enumerate(loader_train):\n",
    "        batch_X, batch_y = Variable(batch_X), Variable(batch_y)\n",
    "        \n",
    "        # Forward pass: compute predicted y using operations on Variables;\n",
    "        hidden_activations = batch_X.mm(w).sigmoid()\n",
    "        y_pred = hidden_activations.mm(v).sigmoid()\n",
    "        # Compute and print loss using operations on Variables.\n",
    "        # Now loss is a Variable of shape (1,) and loss.data is a Tensor of shape\n",
    "        # (1,); loss.data[0] is a scalar value holding the loss.\n",
    "        loss = - (batch_y * torch.log(y_pred) + (1 - batch_y) * torch.log(1 - y_pred)).sum()\n",
    "        average_loss += loss.data[0]\n",
    "        \n",
    "        # Use autograd to compute the backward pass. This call will compute the\n",
    "        # gradient of loss with respect to all Variables with requires_grad=True.\n",
    "        # After this call w.grad and v.grad will be Variables holding the gradient\n",
    "        # of the loss with respect to w and v respectively.\n",
    "        loss.backward()\n",
    "\n",
    "        # Update weights using gradient descent; w1.data and w2.data are Tensors,\n",
    "        # w1.grad and w2.grad are Variables and w1.grad.data and w2.grad.data are\n",
    "        # Tensors.\n",
    "        w.data -= learning_rate * w.grad.data\n",
    "        v.data -= learning_rate * v.grad.data\n",
    "\n",
    "        # Manually zero the gradients after updating weights\n",
    "        w.grad.data.zero_()\n",
    "        v.grad.data.zero_()\n",
    "    if iteration % 5 == 0:\n",
    "        print iteration, '\\t', average_loss / len(trainY)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Predict test samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.58307875856574343"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "testX_var = Variable(torch.FloatTensor(testX_scaled))\n",
    "proba_var = testX_var.mm(w).sigmoid().mm(v).sigmoid()\n",
    "roc_auc_score(testY, proba_var.data.numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** how to modify previous code to have SGD with momentum?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple NN with layers and optimization methods\n",
    "\n",
    "All the needed components are written already for us, so let's just use them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.66745367231\n",
      "1 0.523617102622\n",
      "2 0.511197137234\n",
      "3 0.509856051213\n",
      "4 0.509322984625\n",
      "5 0.509123538793\n",
      "6 0.508960434119\n",
      "7 0.508864008313\n",
      "8 0.508786980206\n",
      "9 0.508743424029\n",
      "10 0.508654685386\n",
      "11 0.508629862152\n",
      "12 0.508581621509\n",
      "13 0.508713238715\n",
      "14 0.50860274833\n",
      "15 0.508476358781\n",
      "16 0.508499256466\n",
      "17 0.508510646376\n",
      "18 0.508560741483\n",
      "19 0.50847081916\n"
     ]
    }
   ],
   "source": [
    "batch_size, dim_input, dim_hidden, dim_output = 64, len(trainX.colnames), 15, 1\n",
    "\n",
    "# Use the nn package to define our model and loss function.\n",
    "# sequence of layers\n",
    "model = torch.nn.Sequential(\n",
    "    torch.nn.Linear(dim_input, dim_hidden),\n",
    "    torch.nn.Sigmoid(),\n",
    "    torch.nn.Linear(dim_hidden, dim_output),\n",
    ")\n",
    "# define logistic loss (take NN output, apply sigmoid and compute logistic loss)\n",
    "loss_fn = torch.nn.SoftMarginLoss(size_average=False)\n",
    "\n",
    "# Use the optim package to define an Optimizer that will update the weights of\n",
    "# the model for us. Here we will use Adam optimization methos; the optim package contains many other\n",
    "# optimization algoriths. The first argument to the Adam constructor tells the\n",
    "# optimizer which Variables it should update.\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)\n",
    "\n",
    "for iteration in range(20):\n",
    "    average_loss = 0\n",
    "    # use optimization over batch\n",
    "    for nbatch, (batch_X, batch_y) in enumerate(loader_train):\n",
    "        batch_X, batch_y = Variable(batch_X), Variable(batch_y)\n",
    "        # loss function takes {-1, 1} as target\n",
    "        batch_y = 2 * batch_y - 1\n",
    "        \n",
    "        # Forward pass: compute predicted y by passing x to the model.\n",
    "        pred = model(batch_X)\n",
    "\n",
    "        # Compute and print loss.\n",
    "        loss = loss_fn(pred, batch_y)\n",
    "        average_loss += loss.data[0]\n",
    "\n",
    "        # Before the backward pass, use the optimizer object to zero all of the\n",
    "        # gradients for the variables it will update \n",
    "        # (which are the learnable weights of the model)\n",
    "        optimizer.zero_grad()\n",
    "\n",
    "        # Backward pass: compute gradient of the loss with respect to model parameters\n",
    "        loss.backward()\n",
    "\n",
    "        # Calling the step function on an Optimizer makes an update of parameters\n",
    "        optimizer.step()\n",
    "    print iteration, average_loss / len(trainX)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.60418980220645291"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from torch.nn.functional import sigmoid\n",
    "\n",
    "testX_var = Variable(torch.FloatTensor(testX_scaled))\n",
    "proba_var = sigmoid(model(testX_var))\n",
    "roc_auc_score(testY, proba_var.data.numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# References\n",
    "- [Paper](https://www.csie.ntu.edu.tw/~b97053/paper/Factorization%20Machines%20with%20libFM.pdf) on factorization machines\n",
    "- [Pytorch documentation](http://pytorch.org/tutorials/index.html)"
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}