MIT_6_036_HW05_Colab.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/DanyunHe/MIT6.036/blob/master/MIT_6_036_HW05_Colab.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nU3rlhVsaHpo"
   },
   "source": [
    "#MIT 6.036 Spring 2019: Homework 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "faW-t7llapsY"
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "VkL4viEaJq7u"
   },
   "source": [
    "## Setup\n",
    "First, download the code distribution for this homework that contains test cases and helper functions.\n",
    "\n",
    "Run the next code block to download and import the code for this lab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "vcyC3YnGJybh"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/bin/sh: wget: command not found\n",
      "unzip:  cannot find or open code_and_data_for_hw05.zip, code_and_data_for_hw05.zip.zip or code_and_data_for_hw05.zip.ZIP.\n",
      "mv: rename code_and_data_for_hw05/* to ./*: No such file or directory\n"
     ]
    },
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'code_for_hw5'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-ecd18cc42b9b>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msystem\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'mv code_and_data_for_hw05/* .'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0;32mimport\u001b[0m \u001b[0mcode_for_hw5\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mhw5\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'code_for_hw5'"
     ]
    }
   ],
   "source": [
    "!rm -rf code_and_data_for_hw05*\n",
    "!wget --quiet https://introml.odl.mit.edu/cat-soop/_static/6.036/homework/hw05/code_and_data_for_hw05.zip\n",
    "!unzip code_and_data_for_hw05.zip\n",
    "!mv code_and_data_for_hw05/* .\n",
    "  \n",
    "import code_for_hw5 as hw5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lamb1=np.arange(0,0.1,0.01)\n",
    "lamb2=np.arange(0,200,20)\n",
    "\n",
    "lamb1[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "RMSE1=np.zeros([2,2,10]) \n",
    "for feature in [0,1]:\n",
    "      for order in [0,1]:\n",
    "            for l,lam in enumerate(lamb1):\n",
    "                  RMSE1[feature,order,l]=feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "w1=np.array([[1],[0]])\n",
    "w2=np.array([[-1],[0]])\n",
    "X = np.array([[0, 1, 2],\n",
    "              [0, 1, 2]])\n",
    "Y = np.array([[0, 1, 0]])           "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.5  0.5  1.5]] [[ 1.5  0.5 -0.5]]\n"
     ]
    }
   ],
   "source": [
    "print(w1.T@X-0.5,w2.T@X+1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "0OYsAwoGasfU"
   },
   "source": [
    "## 6) Linear Regression  - going downhill\n",
    "\n",
    "We will now write some general Python code to compute the gradient of the squared-loss objective, following the structure of the expression, and the rules of calculus.  Note that this style of writing the gradient functions maps directly into the chain-rule steps required to compute the gradient, but produces code that is inefficient, because of duplicated computations. It is straightforward to implement more efficient versions if you want to use them for larger problems.\n",
    "\n",
    "### 6.1) Some basic functions\n",
    "\n",
    "We start by defining some basic functions for computing the mean squared loss.  Note that we want these to work for any\n",
    "value of $n$. That is, `x` could be a single feature vector or a full data matrix, and similarly for `y`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "a_WFv1pzd30E"
   },
   "outputs": [],
   "source": [
    "# In all the following definitions:\n",
    "# x is d by n : input data\n",
    "# y is 1 by n : output regression values\n",
    "# th is d by 1 : weights\n",
    "# th0 is 1 by 1 or scalar\n",
    "def lin_reg(x, th, th0):\n",
    "    \"\"\" Returns the predicted y\n",
    "\n",
    "    >>> X = np.array([[ 1.,  2.,  3.,  4.], [ 1.,  1.,  1.,  1.]])\n",
    "    >>> th = np.array([[ 1.  ], [ 0.05]]) ; th0 = np.array([[ 0.]])\n",
    "    >>> lin_reg(X, th, th0).tolist()\n",
    "    [[1.05, 2.05, 3.05, 4.05]]\n",
    "    >>> th = np.array([[ 1.  ], [ 0.05]]) ; th0 = np.array([[ 2.]])\n",
    "    >>> lin_reg(X, th, th0).tolist()\n",
    "    [[3.05, 4.05, 5.05, 6.05]]\n",
    "    \"\"\"\n",
    "    return np.dot(th.T, x) + th0\n",
    "def square_loss(x, y, th, th0):\n",
    "    \"\"\" Returns the squared loss between y_pred and y\n",
    "\n",
    "    >>> X = np.array([[ 1.,  2.,  3.,  4.], [ 1.,  1.,  1.,  1.]])\n",
    "    >>> Y = np.array([[ 1. ,  2.2,  2.8,  4.1]])\n",
    "    >>> th = np.array([[ 1.  ], [ 0.05]]) ; th0 = np.array([[ 2.]])\n",
    "    >>> square_loss(X, Y, th, th0).tolist()\n",
    "    [[4.2025, 3.4224999999999985, 5.0625, 3.8025000000000007]]\n",
    "    \"\"\"\n",
    "    return (y - lin_reg(x, th, th0))**2\n",
    "def mean_square_loss(x, y, th, th0):\n",
    "    \"\"\" Return the mean squared loss between y_pred and y\n",
    "\n",
    "    >>> X = np.array([[ 1.,  2.,  3.,  4.], [ 1.,  1.,  1.,  1.]])\n",
    "    >>> Y = np.array([[ 1. ,  2.2,  2.8,  4.1]])\n",
    "    >>> th = np.array([[ 1.  ], [ 0.05]]) ; th0 = np.array([[ 2.]])\n",
    "    >>> mean_square_loss(X, Y, th, th0).tolist()\n",
    "    [[4.1225]]\n",
    "    \"\"\"\n",
    "    # the axis=1 and keepdims=True are important when x is a full matrix\n",
    "    return np.mean(square_loss(x, y, th, th0), axis = 1, keepdims = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Oe9BjKuwd7uL"
   },
   "source": [
    "### 6.2) Gradients with respect to $\\theta$\n",
    "\n",
    "Now, let's compute the gradients with respect to $\\theta$. Make sure that they work both for data matrices and label vectors.  You can write one function at a time, some of the checks will apply to each function independently."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "VrdkQNFUf3fD"
   },
   "outputs": [],
   "source": [
    "# Write a function that returns the gradient of lin_reg(x, th, th0)\n",
    "# with respect to th\n",
    "def d_lin_reg_th(x, th, th0):\n",
    "    \"\"\" Returns the gradient of lin_reg(x, th, th0) with respect to th\n",
    "\n",
    "    Note that for array (rather than vector) x, we get a d x n \n",
    "    result. That is to say, this function produces the gradient for\n",
    "    each data point i ... n, with respect to each theta, j ... d.\n",
    "\n",
    "    >>> X = np.array([[ 1.,  2.,  3.,  4.], [ 1.,  1.,  1.,  1.]])\n",
    "    >>> th = np.array([[ 1.  ], [ 0.05]]); th0 = np.array([[ 2.]])\n",
    "    >>> d_lin_reg_th(X[:,0:1], th, th0).tolist()\n",
    "    [[1.0], [1.0]]\n",
    "\n",
    "    >>> d_lin_reg_th(X, th, th0).tolist()\n",
    "    [[1.0, 2.0, 3.0, 4.0], [1.0, 1.0, 1.0, 1.0]]\n",
    "    \"\"\"\n",
    "    return None\n",
    "\n",
    "# Write a function that returns the gradient of square_loss(x, y, th, th0) with\n",
    "# respect to th.  It should be a one-line expression that uses lin_reg and\n",
    "# d_lin_reg_th.\n",
    "def d_square_loss_th(x, y, th, th0):\n",
    "    \"\"\"Returns the gradient of square_loss(x, y, th, th0) with respect to\n",
    "       th.\n",
    "\n",
    "       Note: should be a one-line expression that uses lin_reg and\n",
    "       d_lin_reg_th (i.e., uses the chain rule).\n",
    "\n",
    "       Should work with X, Y as vectors, or as arrays. As in the\n",
    "       discussion of d_lin_reg_th, this should give us back an n x d\n",
    "       array -- so we know the sensitivity of square loss for each\n",
    "       data point i ... n, with respect to each element of theta.\n",
    "\n",
    "    >>> X = np.array([[ 1.,  2.,  3.,  4.], [ 1.,  1.,  1.,  1.]])\n",
    "    >>> Y = np.array([[ 1. ,  2.2,  2.8,  4.1]])\n",
    "    >>> th = np.array([[ 1.  ], [ 0.05]]) ; th0 = np.array([[ 2.]])\n",
    "    >>> d_square_loss_th(X[:,0:1], Y[:,0:1], th, th0).tolist()\n",
    "    [[4.1], [4.1]]\n",
    "\n",
    "    >>> d_square_loss_th(X, Y, th, th0).tolist()\n",
    "    [[4.1, 7.399999999999999, 13.5, 15.600000000000001], [4.1, 3.6999999999999993, 4.5, 3.9000000000000004]]\n",
    "\n",
    "    \"\"\"\n",
    "    return None\n",
    "\n",
    "# Write a function that returns the gradient of mean_square_loss(x, y, th, th0) with\n",
    "# respect to th.  It should be a one-line expression that uses d_square_loss_th.\n",
    "def d_mean_square_loss_th(x, y, th, th0):\n",
    "    \"\"\" Returns the gradient of mean_square_loss(x, y, th, th0) with\n",
    "        respect to th.  \n",
    "\n",
    "        Note: It should be a one-line expression that uses d_square_loss_th.\n",
    "\n",
    "    >>> X = np.array([[ 1.,  2.,  3.,  4.], [ 1.,  1.,  1.,  1.]])\n",
    "    >>> Y = np.array([[ 1. ,  2.2,  2.8,  4.1]])\n",
    "    >>> th = np.array([[ 1.  ], [ 0.05]]) ; th0 = np.array([[ 2.]])\n",
    "    >>> d_mean_square_loss_th(X[:,0:1], Y[:,0:1], th, th0).tolist()\n",
    "    [[4.1], [4.1]]\n",
    "\n",
    "    >>> d_mean_square_loss_th(X, Y, th, th0).tolist()\n",
    "    [[10.15], [4.05]]\n",
    "    \"\"\"\n",
    "    # print(\"X =\", repr(X))\n",
    "    # print(\"Y =\", repr(Y))\n",
    "    # print(\"th =\", repr(th), \"th0 =\", repr(th0))\n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UqQVLjZ1gGvd"
   },
   "source": [
    "### 6.3) Gradients with respect to $\\theta_0$\n",
    "\n",
    "Now, let's compute the gradients with respect to $\\theta_0$. Make sure that they work both for data matrices and label vectors.  You can write one function at a time, some of the checks will apply to each function independently."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "RBUSolgxgSgx"
   },
   "outputs": [],
   "source": [
    "# Write a function that returns the gradient of lin_reg(x, th, th0)\n",
    "# with respect to th0. Hint: Think carefully about what the dimensions of the returned value should be!\n",
    "def d_lin_reg_th0(x, th, th0):\n",
    "    \"\"\" Returns the gradient of lin_reg(x, th, th0) with respect to th0.\n",
    "\n",
    "    >>> x = np.array([[ 1.,  2.,  3.,  4.], [ 1.,  1.,  1.,  1.]])\n",
    "    >>> th = np.array([[ 1.  ], [ 0.05]]) ; th0 = np.array([[ 2.]])\n",
    "    >>> d_lin_reg_th0(x, th, th0).tolist()\n",
    "    [[1.0, 1.0, 1.0, 1.0]]\n",
    "    \"\"\"\n",
    "    return None\n",
    "\n",
    "# Write a function that returns the gradient of square_loss(x, y, th, th0) with\n",
    "# respect to th0.  It should be a one-line expression that uses lin_reg and\n",
    "# d_lin_reg_th0.\n",
    "def d_square_loss_th0(x, y, th, th0):\n",
    "    \"\"\" Returns the gradient of square_loss(x, y, th, th0) with\n",
    "        respect to th0.\n",
    "\n",
    "    # Note: uses broadcasting!\n",
    "\n",
    "    >>> X = np.array([[ 1.,  2.,  3.,  4.], [ 1.,  1.,  1.,  1.]])\n",
    "    >>> Y = np.array([[ 1. ,  2.2,  2.8,  4.1]])\n",
    "    >>> th = np.array([[ 1.  ], [ 0.05]]) ; th0 = np.array([[ 2.]])\n",
    "    >>> d_square_loss_th0(X, Y, th, th0).tolist()\n",
    "    [[4.1, 3.6999999999999993, 4.5, 3.9000000000000004]]\n",
    "    \"\"\"\n",
    "    return None\n",
    "\n",
    "# Write a function that returns the gradient of mean_square_loss(x, y, th, th0) with\n",
    "# respect to th0.  It should be a one-line expression that uses d_square_loss_th0.\n",
    "def d_mean_square_loss_th0(x, y, th, th0):\n",
    "    \"\"\" Returns the gradient of mean_square_loss(x, y, th, th0) with\n",
    "    respect to th0.\n",
    "\n",
    "    >>> X = np.array([[ 1.,  2.,  3.,  4.], [ 1.,  1.,  1.,  1.]])\n",
    "    >>> Y = np.array([[ 1. ,  2.2,  2.8,  4.1]])\n",
    "    >>> th = np.array([[ 1.  ], [ 0.05]]) ; th0 = np.array([[ 2.]])\n",
    "    >>> d_mean_square_loss_th0(X, Y, th, th0).tolist()\n",
    "    [[4.05]]\n",
    "    \"\"\"\n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "cb1FsfCTgjMq"
   },
   "source": [
    "## 7) Going down the ridge\n",
    "\n",
    "Now, let's add a regularizer.  The ridge objective can be implemented as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "tIl6pE-ziLFb"
   },
   "outputs": [],
   "source": [
    "# In all the following definitions:\n",
    "# x is d by n : input data\n",
    "# y is 1 by n : output regression values\n",
    "# th is d by 1 : weights\n",
    "# th0 is 1 by 1 or scalar\n",
    "def ridge_obj(x, y, th, th0, lam):\n",
    "    \"\"\" Return the ridge objective value\n",
    "\n",
    "    >>> X = np.array([[ 1.,  2.,  3.,  4.], [ 1.,  1.,  1.,  1.]])\n",
    "    >>> Y = np.array([[ 1. ,  2.2,  2.8,  4.1]])\n",
    "    >>> th = np.array([[ 1.  ], [ 0.05]]) ; th0 = np.array([[ 2.]])\n",
    "    >>> ridge_obj(X, Y, th, th0, 0.0).tolist()\n",
    "    [[4.1225]]\n",
    "    >>> ridge_obj(X, Y, th, th0, 0.5).tolist()\n",
    "    [[4.623749999999999]]\n",
    "    >>> ridge_obj(X, Y, th, th0, 100.).tolist()\n",
    "    [[104.37250000000002]]\n",
    "    \"\"\"\n",
    "    return np.mean(square_loss(x, y, th, th0), axis = 1, keepdims = True) + lam * np.linalg.norm(th)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "R314XINoigFB"
   },
   "source": [
    "Let's extend our previous code for the gradient of the mean square loss to compute the gradient of the ridge objective w\n",
    "ith respect to $\\theta$.  Our previous solutions for the non-ridge case: `d_mean_square_loss_th` and `d_mean_square_loss\n",
    "_th0` will be defined for you in the grader, so feel free to call them!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "XQ4Dy26vijP4"
   },
   "outputs": [],
   "source": [
    "def d_ridge_obj_th(x, y, th, th0, lam):\n",
    "    \"\"\"Return the derivative of tghe ridge objective value with respect\n",
    "    to theta.\n",
    "\n",
    "    Note: uses broadcasting to add d x n to d x 1 array below\n",
    "\n",
    "    >>> X = np.array([[ 1.,  2.,  3.,  4.], [ 1.,  1.,  1.,  1.]])\n",
    "    >>> Y = np.array([[ 1. ,  2.2,  2.8,  4.1]])\n",
    "    >>> th = np.array([[ 1.  ], [ 0.05]]) ; th0 = np.array([[ 2.]])\n",
    "    >>> d_ridge_obj_th(X, Y, th, th0, 0.0).tolist()\n",
    "    [[10.15], [4.05]]\n",
    "    >>> d_ridge_obj_th(X, Y, th, th0, 0.5).tolist()\n",
    "    [[11.15], [4.1]]\n",
    "    >>> d_ridge_obj_th(X, Y, th, th0, 100.).tolist()\n",
    "    [[210.15], [14.05]]\n",
    "    \"\"\"\n",
    "    return None\n",
    "\n",
    "def d_ridge_obj_th0(x, y, th, th0, lam):\n",
    "    \"\"\"Return the derivative of tghe ridge objective value with respect\n",
    "    to theta.\n",
    "\n",
    "    Note: uses broadcasting to add d x n to d x 1 array below\n",
    "\n",
    "    >>> X = np.array([[ 1.,  2.,  3.,  4.], [ 1.,  1.,  1.,  1.]])\n",
    "    >>> Y = np.array([[ 1. ,  2.2,  2.8,  4.1]])\n",
    "    >>> th = np.array([[ 1.  ], [ 0.05]]) ; th0 = np.array([[ 2.]])\n",
    "    >>> d_ridge_obj_th0(X, Y, th, th0, 0.0).tolist()\n",
    "    [[4.05]]\n",
    "    >>> d_ridge_obj_th0(X, Y, th, th0, 0.5).tolist()\n",
    "    [[4.05]]\n",
    "    >>> d_ridge_obj_th0(X, Y, th, th0, 100.).tolist()\n",
    "    [[4.05]]\n",
    "    \"\"\"\n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "vO2Ebhn3i33z"
   },
   "source": [
    "## 8) Stochastic gradient\n",
    "\n",
    "We will now implement stochastic gradient descent in a general way,\n",
    "similar to what we did with gradient descent (`gd`).\n",
    "\n",
    "The calling conventions for `sgd` are similar to those of `gd` except\n",
    "that we need to pass in the data and labels for the problem.\n",
    "\n",
    "(Recall that the *stochastic* part refers to using a randomly selected point and\n",
    "corresponding label from the given dataset to perform an update. Therefore, your objective function for a given step wil\n",
    "l need to take this into account.)\n",
    "\n",
    "* `X`: a standard data array (d by n)\n",
    "* `y`: a standard labels row vector (1 by n)\n",
    "* `J`: a cost function whose input is a data point (a column vector), a\n",
    "  label (1 by 1) and a weight vector `w` (a column vector) (in that order), and which returns a\n",
    "  scalar.\n",
    "* `dJ`: a cost function gradient (corresponding to `J`) whose input is a data point (a column vector), a\n",
    "  label (1 by 1) and a weight vector `w` (a column vector) (also in that order), and which returns a\n",
    "  column vector.\n",
    "* `w0`: an initial value of weight vector $w$, which is a column vector.\n",
    "* `step_size_fn`: a function that is given the (zero-indexed) iteration index (an\n",
    "  integer) and returns a step size.\n",
    "* `max_iter`: the number of iterations to perform\n",
    "\n",
    "It returns a tuple (like `gd`):\n",
    "\n",
    "* `w`: the value of the weight vector at the final step\n",
    "* `fs`: the list of values of $J$ found during all the iterations\n",
    "* `ws`: the list of values of $w$ found during all the iterations\n",
    "\n",
    "**Note:** `w` should be the value one gets after applying stochastic\n",
    "gradient descent to `w0` for `max_iter-1` iterations (we call this the\n",
    "final step). The first element of `fs` should be the value of `J`\n",
    "calculated with `w0`, and `fs` should have length `max_iter`;\n",
    "similarly, the first element of `ws` should be `w0`, and `ws` should\n",
    "have length `max_iter`.\n",
    "\n",
    "You might find the function `np.random.randint(n)` useful in your implementation.\n",
    "\n",
    "**Hint:** This is a short function; our implementation is around 10\n",
    "  lines.\n",
    "  \n",
    "  The main function to implement is below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "oU1DoLD5jYYQ"
   },
   "outputs": [],
   "source": [
    "def sgd(X, y, J, dJ, w0, step_size_fn, max_iter):\n",
    "    \"\"\"Implements stochastic gradient descent\n",
    "\n",
    "    Inputs:\n",
    "    X: a standard data array (d by n)\n",
    "    y: a standard labels row vector (1 by n)\n",
    "\n",
    "    J: a cost function whose input is a data point (a column vector),\n",
    "    a label (1 by 1) and a weight vector w (a column vector) (in that\n",
    "    order), and which returns a scalar.\n",
    "\n",
    "    dJ: a cost function gradient (corresponding to J) whose input is a\n",
    "    data point (a column vector), a label (1 by 1) and a weight vector\n",
    "    w (a column vector) (also in that order), and which returns a\n",
    "    column vector.\n",
    "\n",
    "    w0: an initial value of weight vector www, which is a column\n",
    "    vector.\n",
    "\n",
    "    step_size_fn: a function that is given the (zero-indexed)\n",
    "    iteration index (an integer) and returns a step size.\n",
    "\n",
    "    max_iter: the number of iterations to perform\n",
    "\n",
    "    Returns: a tuple (like gd):\n",
    "    w: the value of the weight vector at the final step\n",
    "    fs: the list of values of JJJ found during all the iterations\n",
    "    ws: the list of values of www found during all the iterations\n",
    "\n",
    "    \"\"\"\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "FMiqXsP3jZNn"
   },
   "source": [
    "The test cases for this problem are provided below (but, as always, you are encouraged to write more if you want to better test your code!). They rely on the function num_grad (taken from the previous week's homework), also provided."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "X9K32gL4jyr4"
   },
   "outputs": [],
   "source": [
    "def num_grad(f):\n",
    "    def df(x):\n",
    "        g = np.zeros(x.shape)\n",
    "        delta = 0.001\n",
    "        for i in range(x.shape[0]):\n",
    "            xi = x[i,0]\n",
    "            x[i,0] = xi - delta\n",
    "            xm = f(x)\n",
    "            x[i,0] = xi + delta\n",
    "            xp = f(x)\n",
    "            x[i,0] = xi\n",
    "            g[i,0] = (xp - xm)/(2*delta)\n",
    "        return g\n",
    "    return df\n",
    "  \n",
    "def downwards_line():\n",
    "    X = np.array([[0.0, 0.1, 0.2, 0.3, 0.42, 0.52, 0.72, 0.78, 0.84, 1.0],\n",
    "                  [1.0, 1.0, 1.0, 1.0, 1.0,  1.0,  1.0,  1.0,  1.0,  1.0]])\n",
    "    y = np.array([[0.4, 0.6, 1.2, 0.1, 0.22, -0.6, -1.5, -0.5, -0.5, 0.0]])\n",
    "    return X, y\n",
    "\n",
    "X, y = downwards_line()\n",
    "\n",
    "def J(Xi, yi, w):\n",
    "    # translate from (1-augmented X, y, theta) to (separated X, y, th, th0) format\n",
    "    return float(ridge_obj(Xi[:-1,:], yi, w[:-1,:], w[-1:,:], 0))\n",
    "\n",
    "def dJ(Xi, yi, w):\n",
    "    def f(w): return J(Xi, yi, w)\n",
    "    return num_grad(f)(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "G9LFTpdkHKNv"
   },
   "source": [
    "##9) Predicting mpg values\n",
    "\n",
    "We will now try to synthesize the functions we have written in order to perform ridge regression on the <a href=\"https://docs.google.com/spreadsheets/d/1NapkzC_jPQjFk5EfIXghiJJu1kNCN8wP6vwGLpG_tYk/edit?usp=sharing\">auto-mpg dataset</a> from <a href=\"https://lms.mitx.mit.edu/courses/course-v1:MITx+6.036+2019_Spring/courseware/Week3/week3_lab/2\">lab03</a>. Unlike in lab03, we will now try to predict the actual mpg values of the cars, instead of whether they are above or below the median mpg!\n",
    "\n",
    "As a reminder, the dataset is as follows:\n",
    "    \n",
    "    1. mpg:           continuous\n",
    "    2. cylinders:     multi-valued discrete\n",
    "    3. displacement:  continuous\n",
    "    4. horsepower:    continuous\n",
    "    5. weight:        continuous\n",
    "    6. acceleration:  continuous\n",
    "    7. model year:    multi-valued discrete\n",
    "    8. origin:        multi-valued discrete\n",
    "    9. car name:      string (many values)\n",
    "\n",
    "For convenience, we will choose to not include `model year` and `car name` as features. For the remaining features, we again have the option to keep the raw values, standardize them, or use a one-hot encoding.\n",
    "\n",
    "With this considered, we decide to standardize or one-hot encode all features in this section (we encourage you, though, to try raw features on your own time to see how their performance matches your expectations!).\n",
    "\n",
    "One additional step we perform is to standardize the output values. Note that we did not have to worry about this in a classification context, as all outputs were $\\pm 1$. In a regression context, standardizing the output values can have practical performance gains, again due to better numerical performance of learning algorithms on data which is smaller in magnitude.\n",
    "\n",
    "The metric we will use to measure the quality of our learned predictors is ** Root Mean Square Error (RMSE). ** RMSE is defined as follows:\n",
    "\n",
    "$$ \\text{RMSE} = \\sqrt{ \\frac{1}{n} \\sum_{i=1}^n \\left( y^{(i)} - f(x^{(i)}) \\right)^2 } $$\n",
    "\n",
    "where $f$ is our learned predictor: in this case, $f(x) = \\theta \\cdot x + \\theta_0$. This gives a measure of how far away the true values are from the predicted values, measured in units of mpg.\n",
    "\n",
    "** Note: ** One very important thing to keep in mind when employing standardization is that we need to reverse the standardization when we want to report results. If we standardize output values in the training set by subtracting $\\mu$ and dividing by $\\sigma$, we need to take care to:\n",
    "\n",
    "1. Perform standardization with the same values of $\\mu$ and $\\sigma$ on the test set (Why?) before predicting outputs using our learned predictor.\n",
    "2. Multiply the RMSE calculated on the test set by a factor of $\\sigma$ to report test error (Why?)\n",
    "\n",
    "Given all of this, we now will try using:\n",
    "\n",
    "* Two choices of feature set:\n",
    "\n",
    "  1. `[cylinders=standard, displacement=standard, horsepower=standard, weight=standard, acceleration=standard, origin=one_hot]`\n",
    "  2. `[cylinders=one_hot, displacement=standard, horsepower=standard, weight=standard, acceleration=standard, origin=one_hot]`\n",
    "\n",
    "* Polynomial features (we will construct the polynomial features after having standardized the input data) of orders 1-3\n",
    "\n",
    "* Different choices of the regularization parameter, $\\lambda$. Although, ideally, you would run a grid search over a large range of $\\lambda$, we will ask you to look at the choices $\\lambda = \\{0.01, 0.02, \\cdots, 0.1\\}$ for polynomial features of orders $1$ and $2$, and the choices $\\lambda = \\{20, 40, \\cdots, 200\\}$ for polynomial features of order $3$ (as this is approximately where we found the optimal $\\lambda$ to lie).\n",
    "\n",
    "We will use $10$-fold cross-validation to try all possible combinations of these feature choices and test which is best.\n",
    "\n",
    "Your functions written above will be called by `ridge_min`, (defined for you below), which takes a dataset $(X, y)$ and a hyperparameter, $\\lambda$ as input and returns $\\theta$ and $\\theta_0$ minimizing the ridge regression objective using SGD (this is the analogue of the `svm_min` function that you wrote for homework last week). The learning rate and number of iterations are fixed in this function, and should not be modified for the purpose of answering the below questions (although you should feel free to experiment with these if you are interested!) This function will then further be called by `xval_learning_alg` (also defined below), which returns the average RMSE across all (here, 10) splits of your data when performing cross-validation.\n",
    "\n",
    "**Note**: Even though these functions are also contained in the code file being imported (`code_for_hw5.py`), you should run the below code block so that they will use the version of the functions you have written above, and not the blank versions in the code file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Uy6R1Vhou1eG"
   },
   "outputs": [],
   "source": [
    "#Concatenates the gradients with respect to theta and theta_0\n",
    "def ridge_obj_grad(x, y, th, th0, lam):\n",
    "    grad_th = d_ridge_obj_th(x, y, th, th0, lam)\n",
    "    grad_th0 = d_ridge_obj_th0(x, y, th, th0, lam)\n",
    "    return np.vstack([grad_th, grad_th0])    \n",
    "  \n",
    "def ridge_min(X, y, lam):\n",
    "    \"\"\" Returns th, th0 that minimize the ridge regression objective\n",
    "\n",
    "    Assumes that X is NOT 1-extended. Interfaces to our sgd by 1-extending\n",
    "    and building corresponding initial weights.\n",
    "    \"\"\"\n",
    "    def svm_min_step_size_fn(i):\n",
    "        return 0.01/(i+1)**0.5\n",
    "\n",
    "    d, n = X.shape\n",
    "    X_extend = np.vstack([X, np.ones((1, n))])\n",
    "    w_init = np.zeros((d+1, 1))\n",
    "\n",
    "    def J(Xj, yj, th):\n",
    "        return float(ridge_obj(Xj[:-1,:], yj, th[:-1,:], th[-1:,:], lam))\n",
    "\n",
    "    def dJ(Xj, yj, th):\n",
    "        return ridge_obj_grad(Xj[:-1,:], yj, th[:-1,:], th[-1:,:], lam)\n",
    "\n",
    "    np.random.seed(0)\n",
    "    w, fs, ws = sgd(X_extend, y, J, dJ, w_init, svm_min_step_size_fn, 1000)\n",
    "    return w[:-1,:], w[-1:,:]\n",
    "  \n",
    "#First finds a predictor on X_train and X_test using the specified value of lam\n",
    "#Then runs on X_test, Y_test to find the RMSE\n",
    "def eval_predictor(X_train, Y_train, X_test, Y_test, lam):\n",
    "    th, th0 = ridge_min(X_train, Y_train, lam)\n",
    "    return np.sqrt(mean_square_loss(X_test, Y_test, th, th0))\n",
    "  \n",
    "#Returns the mean RMSE from cross validation given a dataset (X, y), a value of lam,\n",
    "#and number of folds, k\n",
    "def xval_learning_alg(X, y, lam, k):\n",
    "    _, n = X.shape\n",
    "    idx = list(range(n))\n",
    "    np.random.seed(0)\n",
    "    np.random.shuffle(idx)\n",
    "    X, y = X[:,idx], y[:,idx]\n",
    "\n",
    "    split_X = np.array_split(X, k, axis=1)\n",
    "    split_y = np.array_split(y, k, axis=1)\n",
    "\n",
    "    score_sum = 0\n",
    "    for i in range(k):\n",
    "        X_train = np.concatenate(split_X[:i] + split_X[i+1:], axis=1)\n",
    "        y_train = np.concatenate(split_y[:i] + split_y[i+1:], axis=1)\n",
    "        X_test = np.array(split_X[i])\n",
    "        y_test = np.array(split_y[i])\n",
    "        score_sum += eval_predictor(X_train, y_train, X_test, y_test, lam)\n",
    "    return score_sum/k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "NYEwfuMdu2hV"
   },
   "source": [
    "Below, there is code for creating the two feature sets that you are asked to work with here. Transforming those features further with `make_polynomial_feature_fun` (can be called from the imported file), and running the cross-validation function (defined above), should allow you to answer the questions for this section.\n",
    "\n",
    " **Make sure that you are calling the cross-validation function from above and not the one from the imported file. **\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "w2DtmogNLQcF"
   },
   "outputs": [],
   "source": [
    "# Returns a list of dictionaries.  Keys are the column names, including mpg.\n",
    "auto_data_all = hw5.load_auto_data('auto-mpg-regression.tsv')\n",
    "\n",
    "# The choice of feature processing for each feature, mpg is always raw and\n",
    "# does not need to be specified.  Other choices are hw5.standard and hw5.one_hot.\n",
    "# 'name' is not numeric and would need a different encoding.\n",
    "features1 = [('cylinders', hw5.standard),\n",
    "            ('displacement', hw5.standard),\n",
    "            ('horsepower', hw5.standard),\n",
    "            ('weight', hw5.standard),\n",
    "            ('acceleration', hw5.standard),\n",
    "            ('origin', hw5.one_hot)]\n",
    "\n",
    "features2 = [('cylinders', hw5.one_hot),\n",
    "            ('displacement', hw5.standard),\n",
    "            ('horsepower', hw5.standard),\n",
    "            ('weight', hw5.standard),\n",
    "            ('acceleration', hw5.standard),\n",
    "            ('origin', hw5.one_hot)]\n",
    "\n",
    "# Construct the standard data and label arrays\n",
    "#auto_data[0] has the features for choice features1\n",
    "#auto_data[1] has the features for choice features2\n",
    "#The labels for both are the same, and are in auto_values\n",
    "auto_data = [0, 0]\n",
    "auto_values = 0\n",
    "auto_data[0], auto_values = hw5.auto_data_and_values(auto_data_all, features1)\n",
    "auto_data[1], _ = hw5.auto_data_and_values(auto_data_all, features2)\n",
    "\n",
    "#standardize the y-values\n",
    "auto_values, mu, sigma = hw5.std_y(auto_values)\n",
    "\n",
    "#-------------------------------------------------------------------------------\n",
    "# Analyze auto data\n",
    "#-------------------------------------------------------------------------------     \n",
    "        \n",
    "#Your code for cross-validation goes here\n",
    "#Make sure to scale the RMSE values returned by xval_learning_alg by sigma,\n",
    "#as mentioned in the lab, in order to get accurate RMSE values on the dataset\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "Copy of MIT 6.036 HW05 Colab",
   "provenance": [],
   "version": "0.3.2"
  },
  "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}