unit03_mult_lin_reg/demo2_python_broadcasting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numpy Array Operations:  Axes and Broadcasting \n",
    "\n",
    "There is an excellent introduction to `numpy` multi-dimensional arrays on the [scipy](https://docs.scipy.org/doc/numpy-dev/user/quickstart.html) website.  In this note, we cover two concepts in a little more detail:\n",
    "* Using the `axis` feature \n",
    "* Python broadcasting\n",
    "\n",
    "We will need both of these for performing many of the numerical operations for the ML class.\n",
    "\n",
    "As usual, we begin by loading the `numpy` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Axis Parameter\n",
    "\n",
    "Many operations in the `numpy` package can take an optional `axis` parameter to specify which dimensions the operation is to be applied.  This is extremely useful for multi-dimensional data.  To illustrate the `axis` parameter, consider a matrix the `(3,2)` array `X` defined as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 1]\n",
      " [2 3]\n",
      " [4 5]]\n"
     ]
    }
   ],
   "source": [
    "X = np.arange(6).reshape(3,2)\n",
    "print(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An operation like `np.mean` or `np.sum` takes the mean or sum of *all* elements in the array. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.5\n",
      "15\n"
     ]
    }
   ],
   "source": [
    "print(np.mean(X))\n",
    "print(np.sum(X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To take only the `sum` along each column, we can use the `axis` parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6 9]\n"
     ]
    }
   ],
   "source": [
    "print(np.sum(X,axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since `X` has shape `(3,2)`, the output `np.sum(X,axis=0)` is of shape `(2,)`.  Similarly, we can take the `sum` along each row:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 5 9]\n"
     ]
    }
   ],
   "source": [
    "print(np.sum(X,axis=1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can apply this to higher-order arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Y1 = \n",
      "[[12 14 16 18]\n",
      " [20 22 24 26]\n",
      " [28 30 32 34]]\n",
      "Y2 = \n",
      "[[12 15 18 21]\n",
      " [48 51 54 57]]\n"
     ]
    }
   ],
   "source": [
    "X = np.arange(24).reshape(2,3,4)  # shape = (2,3,4)\n",
    "Y1 = np.sum(X,axis=0)             # shape = (3,4)\n",
    "Y2 = np.sum(X,axis=1)             # shape = (2,4)\n",
    "print('Y1 = ')\n",
    "print(Y1)\n",
    "print('Y2 = ')\n",
    "print(Y2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Broadcasting\n",
    "\n",
    "**Broadcasting** is a powerful tool in Python for performing operations on matrices that we will use throughout the ML class.  A good tutorial on broadcasting can be found on the [scipy broadcasting page](https://docs.scipy.org/doc/numpy-1.13.0/user/basics.broadcasting.html).  Here we provide some examples.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 1:  Mean Removal\n",
    "\n",
    "Suppose that `X` is a data matrix of shape `(n,p)`.  That is, there are `n` data points and `p` features per point.  Often, we have to remove the mean from each feature.  That is, we want to compute the mean for each feature and then remove the mean from each column.  We could do this with a for-loop as:\n",
    "   \n",
    "    Xm = np.zeros(p)      # Mean for each feature\n",
    "    X1_demean = np.zeros((n,p))  # Transformed features with the means removed\n",
    "    for j in range(p):\n",
    "       Xm[j] = np.mean(X[:,j])\n",
    "       for i in range(n):\n",
    "           X_demean[i,j] = X[i,j] - Xm[j]\n",
    "           \n",
    "The code below does this without a for loop using the `axis` parameter and broadcasting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Generate some random data\n",
    "n = 100\n",
    "p = 5\n",
    "X = np.random.rand(n,p)\n",
    "\n",
    "# Compute the mean per column using the axis command\n",
    "Xm = np.mean(X,axis=0)  # This is a p-dim matrix\n",
    "\n",
    "# Subtract the mean\n",
    "X_demean = X - Xm[None,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The command `Xm = np.mean(X,axis=0)` computes the mean of each column which is a `p` dimensional array.  Then, `Xm[None,:]` converts this to a `(1,p)` shape array.  Using Python broadcasting, we can then subtract the `Xm[None,:]` from `X`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2:  Standardizing variables\n",
    "\n",
    "A variant of the above example is to *standardize* the features, where we compute the transform variables,\n",
    "\n",
    "    Z[i,j] = (X[i,j] - Xm[j])/ Xstd[j]\n",
    "    \n",
    "where `Xstd[j]` is the standard deviation per feature.  This can be done as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Xstd = np.std(X,axis=0)\n",
    "Z = (X-Xm[None,:])/Xstd[None,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 3:  Distances\n",
    "\n",
    "Here is a more complicated example.  Suppose we have a data matrix `X` of shape `(nx,p)` and a second set of points, `Y` of shape `(ny,p)`. For each `i` and `j`, we want to compute the distances, \n",
    "\n",
    "     d[i,j] = np.sum((X[i,:] - Y[j,:])**2)\n",
    "     \n",
    "This represents the distances between the vectors `X[i,:]` and `Y[j,:]`.  This sort of computation is used for clustering and nearest neighbors.  We can do this without a for loop as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Some random data\n",
    "nx = 100\n",
    "ny = 10\n",
    "p = 2\n",
    "X = np.random.rand(nx,p)\n",
    "Y = np.random.rand(ny,p)\n",
    "\n",
    "# Computing the distances in two lines.  No for loop!\n",
    "DXY = X[:,None,:]-Y[None,:,:]\n",
    "d = np.sum(DXY**2,axis=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "How does this work? First, we use `None` keyword to reshape the matrices `X` and `Y` to compatible sizes\n",
    "\n",
    "     X[:,None,:]    # Shape nx,  1, p\n",
    "     Y[None,:,:]    # Shape 1,  ny, p\n",
    "     \n",
    "The two matrices can be subtracted so that\n",
    "\n",
    "     DXY[i,j,k]  = X[i,k] - Y[j,k]\n",
    "     \n",
    "Then, `d[i,j] = sum_k (X[i,k] - Y[j,k])**2`, which is the norm squared of the vector differences."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 4:  Outer product\n",
    "\n",
    "The *outer product* of vectors `x` and `y` is the matrix `Z[i,j] = x[i]y[j]`.  This can be performed in one line as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Some random data\n",
    "nx = 100\n",
    "ny = 10\n",
    "x = np.random.rand(nx)\n",
    "y = np.random.rand(ny)\n",
    "\n",
    "# Compute the outer product in one line\n",
    "Z = x[:,None]*y[None,:]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here:\n",
    "\n",
    "     x[:,None] # Has shape (nx,  1)\n",
    "     y[None,:] # Has shape ( 1, ny)\n",
    "     \n",
    "So, with python broadcasting:\n",
    "\n",
    "     Z = x[:,None]*y[None,:] # has shape (nx,  ny)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Exercise 1:**  Given a matrix `X`, compute the matrix `Y`, where the rows of `X` are normaized to one.  That is:\n",
    "\n",
    "     Y[i,j] = X[i,j] / sum_j X[i,j]   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = np.random.rand(4,3)\n",
    "# Y = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2:** Diagonal multiplication.  Given a matrix `X` and a vector `d`, compute `Y = diag(d)*X`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X = np.random.rand(5,3)\n",
    "d = np.random.rand(5)\n",
    "# Y = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}