AW-SGD Paper Recreation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Imports of all packages used trough out this tutorial\n",
    "%matplotlib inline\n",
    "import os\n",
    "import copy\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.stats import rv_discrete"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recreating the AW-SGD Paper\n",
    "\n",
    "In this notebook I try to replicate the results of the paper \"Accerlerating Stochastic Gradient Descent by Online learning to Sample\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Least Squares Matrix Factorisation\n",
    "\n",
    "The first experiment in the paper that I will attempt to replicate is low rank matrix factorisation by minimisation of the sum square loss between the elements of the reconstructed matrix and the original matrix.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Define Some Small and Useful Functions\n",
    "def matrix_square_loss(U,V,Y):\n",
    "    return np.sum((np.dot(U,V.T) - Y)**2)\n",
    "    \n",
    "def grad_u(U,V,Y,i,j):\n",
    "    \"\"\" Gradient of loss with respect to an element of U\n",
    "    \"\"\"\n",
    "    return 2*(np.dot(U[i,:],V[j,:].T) - Y[i,j])*V[j,:]\n",
    "\n",
    "def grad_v(U,V,Y,i,j):\n",
    "    \"\"\" Gradient of loss with respect to an element of V\n",
    "    \"\"\"\n",
    "    return 2*(np.dot(U[i,:],V[j,:].T) - Y[i,j])*U[i,:]\n",
    "\n",
    "def one_hot(i,N):\n",
    "    \"\"\" Return a vector of length N thats all 0's apart from at place i\n",
    "    \"\"\"\n",
    "    x = np.zeros(N)\n",
    "    x[i] = 1\n",
    "    return x\n",
    "\n",
    "def softmax(tau):\n",
    "    return np.exp(tau)/np.sum(np.exp(tau))\n",
    "\n",
    "\n",
    "# Some functions for sampling from discrete distributions\n",
    "class distribution(rv_discrete):\n",
    "    \n",
    "    def plot(self,N):\n",
    "        xk = np.arange(N)\n",
    "        fig, ax = plt.subplots(1, 1)\n",
    "        ax.plot(xk, self.pmf(xk), 'ro', ms=12, mec='r')\n",
    "        ax.vlines(xk, 0, self.pmf(xk), colors='r', lw=4)\n",
    "        \n",
    "def create_dist(tau):\n",
    "    xk = np.arange(np.size(tau))\n",
    "    pk = softmax(tau)\n",
    "    custm = distribution(name='custm', values=(xk, pk))\n",
    "    return custm\n",
    "\n",
    "# Define Updates for SGD and AWSGD\n",
    "def sgd_update(U,V,Y,i,j,mu):\n",
    "    \"\"\" One Update of SGD\n",
    "    \"\"\"\n",
    "    u_new = U[i,:] - mu*grad_u(U,V,Y,i,j)\n",
    "    v_new = V[j,:] - mu*grad_v(U,V,Y,i,j)\n",
    "    return u_new,v_new\n",
    "\n",
    "def AW_sgd_update(U,V,Y,i,j,mu,dist_u,dist_v,rho,tau1,tau2):\n",
    "    \"\"\" One update of Awsgd\n",
    "    \"\"\"\n",
    "    du = grad_u(U,V,Y,i,j)/dist_u._pmf(i)\n",
    "    dv = grad_v(U,V,Y,i,j)/dist_v._pmf(j)\n",
    "    u_new = U[i,:] - mu*du\n",
    "    v_new = V[j,:] - mu*dv\n",
    "    dtau1 = one_hot(i,N) - softmax(tau1)\n",
    "    dtau2 = one_hot(j,M) - softmax(tau2)\n",
    "    tau1 = tau1 + rho*((np.linalg.norm(du))**2)*dtau1\n",
    "    tau2 = tau2 + rho*((np.linalg.norm(dv))**2)*dtau2\n",
    "    return u_new,v_new,tau1,tau2\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define functions to perform SGD and AWSGD\n",
    "\n",
    "def factor_sgd(M,k,num_iter,initial_mu):\n",
    "    \"\"\" Takes in a Matrix M and returns two matrices A,B such that M ~= AB'\n",
    "        A and B are found by minimising the square loss\n",
    "        \n",
    "        Inputs\n",
    "        ------\n",
    "            M : type - np.array\n",
    "                Matrix to be factorised\n",
    "            k:  type - int\n",
    "                Rank of approximated matrix\n",
    "      num_iter: type - int\n",
    "     Inital_mu: type - float\n",
    "                Initial learning rate\n",
    "        \n",
    "       OutPuts\n",
    "       -------\n",
    "           A,B: np.array\n",
    "                Factorised matrices\n",
    "        errors: np.array\n",
    "                Square error as a function of iteration number\n",
    "    \"\"\"\n",
    "    \n",
    "    # Initialse the factored matrices randomly    \n",
    "    num_rows = np.shape(M)[0]\n",
    "    num_cols = np.shape(M)[1]\n",
    "    A = np.random.randn(num_rows,k)\n",
    "    B = np.random.randn(num_cols,k)\n",
    "\n",
    "    # Storage for the error\n",
    "    errors = np.zeros(num_iter)\n",
    "                \n",
    "    for step in range(num_iter):\n",
    "        # Define Learning Rate\n",
    "        rate = initial_mu/(num_rows/2.0 + 1.0*step)\n",
    "        # Sample a row and a column\n",
    "        i,j = np.random.randint(num_rows), np.random.randint(num_cols)\n",
    "        # Perform one step of sgd\n",
    "        A[i,:],B[j,:] = sgd_update(A,B,M,i,j,rate)\n",
    "        # Calculate and store the error\n",
    "        errors[step] = matrix_square_loss(A,B,M)\n",
    "    \n",
    "    return A,B,errors\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def factor_aswsgd(M,k,num_iter,initial_mu,initial_tau = None):\n",
    "    \"\"\" Takes in a Matrix M and returns two matrices A,B such that M ~= AB'\n",
    "        A and B are found by minimising the square loss\n",
    "        \n",
    "        Inputs\n",
    "        ------\n",
    "            M : type - np.array\n",
    "                Matrix to be factorised\n",
    "            k:  type - int\n",
    "                Rank of approximated matrix\n",
    "      num_iter: type - int\n",
    "     Inital_mu: type - float\n",
    "                Initial learning rate\n",
    "        \n",
    "       OutPuts\n",
    "       -------\n",
    "           A,B: np.array\n",
    "                Factorised matrices\n",
    "        errors: np.array\n",
    "                Square error as a function of iteration number\n",
    "    \"\"\"\n",
    "    \n",
    "    # Initialse the factored matrices randomly    \n",
    "    num_rows = np.shape(M)[0]\n",
    "    num_cols = np.shape(M)[1]\n",
    "    A = np.random.randn(num_rows,k)\n",
    "    B = np.random.randn(num_cols,k)\n",
    "    \n",
    "    # Storage for the error\n",
    "    errors = np.zeros(num_iter)\n",
    "    \n",
    "    #Initialise Tau if not already initialised\n",
    "    if initial_tau is None:\n",
    "        tau_rows = np.random.randn(num_rows)\n",
    "        tau_cols = np.random.randn(num_cols)\n",
    "    else:\n",
    "        tau_rows = initial_tau[0:num_rows]\n",
    "        tau_cols = initial_tau[num_rows + 1:]\n",
    "    \n",
    "    for step in range(num_iter):\n",
    "        # Define the leanring rates\n",
    "        rate = initial_mu/(num_rows/2.0 + 1.0*step)\n",
    "        rho = initial_mu/100000.0\n",
    "        # Sample a row and a column from a soft-max parameterised by tau\n",
    "        row_dist = create_dist(tau_rows)\n",
    "        col_dist = create_dist(tau_cols)\n",
    "        i,j = row_dist.rvs(),col_dist.rvs()\n",
    "        # Perform one step of awsgd\n",
    "        A[i,:],B[j,:],tau_rows,tau_cols = AW_sgd_update(A,B,M,i,j,rate,row_dist,col_dist,rho,tau_rows,tau_cols)\n",
    "        # Calculate and return the Error\n",
    "        errors[step] = matrix_square_loss(A,B,M)\n",
    "    \n",
    "    return A,B,errors,row_dist,col_dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "source": [
    "## Experiments with matrix factorisation SGD vs AW-SGD\n",
    "\n",
    "In the original paper, the authors provide two demonstrations of awsgd for matrix factorisation. First on a randomly generated low rank matrix of dimension 100x100 with rank 10 and second on mnist.\n",
    "\n",
    "For the first experiment the model learning rate was set to $\\frac{\\rho_0}{0.5N +t}$ and the sampling learning rate was kept constant. The results of each experiment were averaged over 10 runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Generate a Low Rank Matrix Y = UV'\n",
    "N = 100 # Dimension 1 of M\n",
    "M = 100 # Dimension 2 0f M \n",
    "K = 10 # Rank of the matrix\n",
    "\n",
    "U = np.random.randn(N,K)\n",
    "V = np.random.randn(M,K)\n",
    "Y = np.dot(U,V.T)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/site-packages/ipykernel/__main__.py:23: RuntimeWarning: overflow encountered in exp\n",
      "/usr/local/lib/python2.7/site-packages/ipykernel/__main__.py:23: RuntimeWarning: invalid value encountered in divide\n",
      "/usr/local/lib/python2.7/site-packages/ipykernel/__main__.py:37: RuntimeWarning: divide by zero encountered in divide\n",
      "/usr/local/lib/python2.7/site-packages/ipykernel/__main__.py:43: RuntimeWarning: invalid value encountered in multiply\n"
     ]
    },
    {
     "ename": "KeyError",
     "evalue": "nan",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-162-60dc41732b3f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     13\u001b[0m     \u001b[0;31m#average_sgd_error = average_sgd_error + esgd\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m     \u001b[0;31m#awsgd\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 15\u001b[0;31m     \u001b[0m_\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0m_\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0meawsgd\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0m_\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0m_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfactor_aswsgd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mY\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnum_iter\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0minitial_mu_awsgd\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     16\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     17\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-126-95c1f3284fa9>\u001b[0m in \u001b[0;36mfactor_aswsgd\u001b[0;34m(M, k, num_iter, initial_mu, initial_tau)\u001b[0m\n\u001b[1;32m     45\u001b[0m         \u001b[0mrow_dist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcreate_dist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtau_rows\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     46\u001b[0m         \u001b[0mcol_dist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcreate_dist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtau_cols\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 47\u001b[0;31m         \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrow_dist\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrvs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcol_dist\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrvs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     48\u001b[0m         \u001b[0;31m# Perform one step of awsgd\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     49\u001b[0m         \u001b[0mA\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mB\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mtau_rows\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mtau_cols\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mAW_sgd_update\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mB\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mM\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mrate\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mrow_dist\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcol_dist\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mrho\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mtau_rows\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mtau_cols\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.pyc\u001b[0m in \u001b[0;36mrvs\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m   2810\u001b[0m         \"\"\"\n\u001b[1;32m   2811\u001b[0m         \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'discrete'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2812\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrv_discrete\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrvs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2813\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2814\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mpmf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.pyc\u001b[0m in \u001b[0;36mrvs\u001b[0;34m(self, *args, **kwds)\u001b[0m\n\u001b[1;32m    869\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_random_state\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcheck_random_state\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrndm\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    870\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 871\u001b[0;31m         \u001b[0mvals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_rvs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    872\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_size\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    873\u001b[0m             \u001b[0mvals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvals\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msize\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.pyc\u001b[0m in \u001b[0;36m_rvs\u001b[0;34m(self, *args)\u001b[0m\n\u001b[1;32m    803\u001b[0m         \u001b[0;31m## Use basic inverse cdf algorithm for RV generation as default.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    804\u001b[0m         \u001b[0mU\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_random_state\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom_sample\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_size\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 805\u001b[0;31m         \u001b[0mY\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_ppf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mU\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    806\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mY\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    807\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python2.7/site-packages/numpy/lib/function_base.pyc\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1809\u001b[0m             \u001b[0mvargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mextend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0m_n\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0m_n\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mnames\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1810\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1811\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_vectorize_call\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mvargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1812\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1813\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_get_ufunc_and_otypes\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python2.7/site-packages/numpy/lib/function_base.pyc\u001b[0m in \u001b[0;36m_vectorize_call\u001b[0;34m(self, func, args)\u001b[0m\n\u001b[1;32m   1878\u001b[0m                       for _a in args]\n\u001b[1;32m   1879\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1880\u001b[0;31m             \u001b[0moutputs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mufunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1881\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1882\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mufunc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnout\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.pyc\u001b[0m in \u001b[0;36m_drv_ppf\u001b[0;34m(self, q, *args)\u001b[0m\n\u001b[1;32m   2349\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_drv_ppf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2350\u001b[0m     \u001b[0mindx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0margmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mqvals\u001b[0m \u001b[0;34m>=\u001b[0m \u001b[0mq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2351\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mFinv\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mqvals\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mindx\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2352\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2353\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: nan"
     ]
    }
   ],
   "source": [
    "# Perform 10 rounds of sgd and awsgd\n",
    "\n",
    "num_iter = 1000\n",
    "average_sgd_error = np.zeros(num_iter)\n",
    "average_awsgd_error = np.zeros(num_iter)\n",
    "\n",
    "initial_mu_sgd = 5.0\n",
    "initial_mu_awsgd = 0.001\n",
    "\n",
    "for iter in range(10):\n",
    "    # sgd\n",
    "    #_,_,esgd = factor_sgd(Y,10,num_iter,initial_mu_sgd)\n",
    "    #average_sgd_error = average_sgd_error + esgd\n",
    "    #awsgd\n",
    "    _,_,eawsgd,_,_ = factor_aswsgd(Y,10,num_iter,initial_mu_awsgd,)\n",
    "    \n",
    "\n",
    "\n",
    "plt.plot(eawsgd/10,c = 'blue')\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plot and compare the average error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}