Notebook10_NumpyMatrix_part2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-467487ba59ea183e",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "# Part 2: Dense matrix storage \n",
    "This part of the lab is a brief introduction to efficient storage of matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 0** (ungraded). Import Numpy!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-4263c0d16078cf0a",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dense matrix storage: Column-major versus row-major layouts\n",
    "\n",
    "For linear algebra, we will be especially interested in 2-D arrays, which we will use to store matrices. For this common case, there is a subtle performance issue related to how matrices are stored in memory.\n",
    "\n",
    "By way of background, physical storage---whether it be memory or disk---is basically one big array. And because of how physical storage is implemented, it turns out that it is much faster to access consecutive elements in memory than, say, to jump around randomly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A matrix is a two-dimensional object. Thus, when it is stored in memory, it must be mapped in some way to the one-dimensional physical array. There are many possible mappings, but the two most common conventions are known as the _column-major_ and _row-major_ layouts:\n",
    "\n",
    "<img src=\"matrix-layout.png\" alt=\"Exercise: Extract these slices\" width=\"640\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1** (2 points). Let $A$ be an $m \\times n$ matrix stored in column-major format. Let $B$ be an $m \\times n$ matrix stored in row-major format.\n",
    "\n",
    "Based on the preceding discussion, recall that these objects will be mapped to 1-D arrays of length $mn$, behind the scenes. Let's call the 1-D array representations $\\hat{A}$ and $\\hat{B}$. Thus, the $(i, j)$ element of $a$, $a_{ij}$, will map to some element $\\hat{a}_u$ of $\\hat{A}$; similarly, $b_{ij}$ will map to some element $\\hat{b}_v$ of $\\hat{B}$.\n",
    "\n",
    "Determine formulae to compute the 1-D index values, $u$ and $v$, in terms of $\\{i, j, m, n\\}$. Assume that all indices are 0-based, i.e., $0 \\leq i \\leq m-1$, $0 \\leq j \\leq n-1$, and $0 \\leq u, v \\leq mn-1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "e628bb9dc9f0e8a68ad52ba1d43caca4",
     "grade": false,
     "grade_id": "calc_u",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def linearize_colmajor(i, j, m, n): # calculate `u`\n",
    "    \"\"\"\n",
    "    Returns the linear index for the `(i, j)` entry of\n",
    "    an `m`-by-`n` matrix stored in column-major order.\n",
    "    \"\"\"\n",
    "    # YOUR CODE HERE\n",
    "    # colmajor_linear_index \n",
    "    u = i + m * j\n",
    "    return u\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "a7ca53a5a658c1c36cbdf1a3ef1a03d9",
     "grade": false,
     "grade_id": "calc_v",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def linearize_rowmajor(i, j, m, n): # calculate `v`\n",
    "    \"\"\"\n",
    "    Returns the linear index for the `(i, j)` entry of\n",
    "    an `m`-by-`n` matrix stored in row-major order.\n",
    "    \"\"\"\n",
    "    # YOUR CODE HERE\n",
    "    # rowmajor_linear_index\n",
    "    v = i * n + j\n",
    "    return v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "917a5dd8ab9f91bc5fbae92435acd42d",
     "grade": true,
     "grade_id": "calc_uv_test",
     "locked": true,
     "points": 2,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(Passed.)\n"
     ]
    }
   ],
   "source": [
    "# Test cell: `calc_uv_test`\n",
    "\n",
    "# Quick check (not exhaustive):\n",
    "assert linearize_colmajor(7, 4, 10, 20) == 47\n",
    "assert linearize_rowmajor(7, 4, 10, 20) == 144\n",
    "\n",
    "assert linearize_colmajor(10, 8, 86, 26) == 698\n",
    "assert linearize_rowmajor(10, 8, 86, 26) == 268\n",
    "\n",
    "assert linearize_colmajor(8, 34, 17, 40) == 586\n",
    "assert linearize_rowmajor(8, 34, 17, 40) == 354\n",
    "\n",
    "assert linearize_colmajor(32, 48, 37, 55) == 1808\n",
    "assert linearize_rowmajor(32, 48, 37, 55) == 1808\n",
    "\n",
    "assert linearize_colmajor(24, 33, 57, 87) == 1905\n",
    "assert linearize_rowmajor(24, 33, 57, 87) == 2121\n",
    "\n",
    "assert linearize_colmajor(10, 3, 19, 74) == 67\n",
    "assert linearize_rowmajor(10, 3, 19, 74) == 743\n",
    "\n",
    "print (\"(Passed.)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Requesting a layout in Numpy\n",
    "\n",
    "In Numpy, you can ask for either layout. The default in Numpy is row-major.\n",
    "\n",
    "Historically numerical linear algebra libraries were developed assuming column-major layout. This layout happens to be the default when you declare a 2-D array in the Fortran programming language. By contrast, in the C and C++ programming languages, the default convention for a 2-D array is row-major layout. So the Numpy default is the C/C++ convention.\n",
    "\n",
    "In your programs, you can request either order of Numpy using the `order` parameter. For linear algebra operations (common), we recommend using the column-major convention.\n",
    "\n",
    "In either case, here is how you would create column- and row-major matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 5000\n",
    "A_colmaj = np.ones((n, n), order='F') # column-major (Fortran convention)\n",
    "A_rowmaj = np.ones((n, n), order='C') # row-major (C/C++ convention)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2** (1 point). Given a matrix $A$, write a function that scales each column, $A(:, j)$ by $j$. Then compare the speed of applying that function to matrices in row and column major order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "8abc750df11036d09bd1e787452d2682",
     "grade": false,
     "grade_id": "scale_colwise",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def scale_colwise(A):\n",
    "    \"\"\"Given a Numpy matrix `A`, visits each column `A[:, j]`\n",
    "    and scales it by `j`.\"\"\"\n",
    "    assert type(A) is np.ndarray\n",
    "    \n",
    "    n_cols = A.shape[1] # number of columns\n",
    "    # YOUR CODE HERE\n",
    "    # A = n_cols * A   # my code...not matrix notation\n",
    "    \n",
    "    # their answer\n",
    "    for j in range(n_cols):\n",
    "        A[:,j] *= j\n",
    "    \n",
    "    return A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "scale_colwise_test",
     "locked": true,
     "points": 1,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "120 ms ± 109 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "31.9 ms ± 25.9 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "# Test (timing) cell: `scale_colwise_test`\n",
    "\n",
    "# Measure time to scale a row-major input column-wise\n",
    "%timeit scale_colwise(A_rowmaj)\n",
    "\n",
    "# Measure time to scale a column-major input column-wise\n",
    "%timeit scale_colwise(A_colmaj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Python vs. Numpy example: Matrix-vector multiply\n",
    "\n",
    "Look at the definition of matrix-vector multiplication from [Da Kuang's linear algebra notes](https://www.dropbox.com/s/f410k9fgd7iesdv/kuang-linalg-notes.pdf?dl=0). Let's benchmark a matrix-vector multiply in native Python, and compare that to doing the same operation in Numpy.\n",
    "\n",
    "First, some setup. (What does this code do?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Dimensions; you might shrink this value for debugging\n",
    "n = 2500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Generate random values, for use in populating the matrix and vector\n",
    "from random import gauss\n",
    "\n",
    "# Native Python, using lists\n",
    "A_py = [gauss(0, 1) for i in range(n*n)] # Assume: Column-major\n",
    "x_py = [gauss(0, 1) for i in range(n)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Convert values into Numpy arrays in column-major order\n",
    "A_np = np.reshape(A_py, (n, n), order='F')\n",
    "x_np = np.reshape(x_py, (n, 1), order='F')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  13.89207985],\n",
       "       [ 104.68760455],\n",
       "       [  76.12002576],\n",
       "       ..., \n",
       "       [   0.41580127],\n",
       "       [ -19.92075803],\n",
       "       [  41.50614829]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_np.dot(x_np)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.49 ms ± 182 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "# Here is how you do a \"matvec\" in Numpy:\n",
    "%timeit A_np.dot(x_np)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3** (3 points). Implement a matrix-vector product that operates on native Python lists. Assume the 1-D **column-major** storage of the matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "nbgrader": {
     "checksum": "407dbfe44ba36ac12065b96d9142e3a4",
     "grade": false,
     "grade_id": "matvec_py",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def matvec_py(m, n, A, x):\n",
    "    \"\"\"\n",
    "    Native Python-based matrix-vector multiply, using lists.\n",
    "    The dimensions of the matrix A are m-by-n, and x is a\n",
    "    vector of length n.\n",
    "    \"\"\"\n",
    "    assert type(A) is list and all([type(aij) is float for aij in A])\n",
    "    assert type(x) is list\n",
    "    assert len(x) >= n\n",
    "    assert len(A) >= (m*n)\n",
    "\n",
    "    y = []\n",
    "    \n",
    "    # YOUR CODE HERE\n",
    "    A_np = np.reshape(A, (m, n), order='F')\n",
    "    \n",
    "    for row in A_np:\n",
    "        row_sum = []\n",
    "        for i in range(m):\n",
    "            row_sum.append(row[i]*x[i])\n",
    "        y.append(sum(row_sum))\n",
    "#         print(row_sum)\n",
    "#     print(y)\n",
    "    return y\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "caf39b3909e91acb9536e726b45bc4cf",
     "grade": true,
     "grade_id": "matvec_py_test",
     "locked": true,
     "points": 3,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==> Error bound estimate:\n",
      "         C*n*eps\n",
      "         == 10*2500*2.22045e-16\n",
      "         == 5.55112e-12\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$$||y_{\\textrm{np}} - y_{\\textrm{py}}||_{\\infty} = \\textrm{5.11591e-13} \\leq \\textrm{5.55112e-12}\\ (\\textrm{estimated bound})$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "(Passed!)\n"
     ]
    }
   ],
   "source": [
    "# Test cell: `matvec_py_test`\n",
    "\n",
    "# Estimate a bound on the difference between these two\n",
    "EPS = np.finfo (float).eps # \"machine epsilon\"\n",
    "CONST = 10.0 # Some constant for the error bound\n",
    "dy_max = CONST * n * EPS\n",
    "\n",
    "print (\"\"\"==> Error bound estimate:\n",
    "         C*n*eps\n",
    "         == %g*%g*%g\n",
    "         == %g\n",
    "\"\"\" % (CONST, n, EPS, dy_max))\n",
    "\n",
    "# Run the Numpy version and your code\n",
    "y_np = A_np.dot (x_np)\n",
    "y_py = matvec_py (n, n, A_py, x_py)\n",
    "\n",
    "# Compute the difference between these\n",
    "dy = y_np - np.reshape (y_py, (n, 1), order='F')\n",
    "dy_norm = np.linalg.norm (dy, ord=np.inf)\n",
    "\n",
    "# Summarize the results\n",
    "from IPython.display import display, Math\n",
    "\n",
    "comparison = \"\\leq\" if dy_norm <= dy_max else \"\\gt\"\n",
    "display (Math (\n",
    "        r'||y_{\\textrm{np}} - y_{\\textrm{py}}||_{\\infty}'\n",
    "        r' = \\textrm{%g} %s \\textrm{%g}\\ (\\textrm{estimated bound})'\n",
    "        % (dy_norm, comparison, dy_max)\n",
    "    ))\n",
    "\n",
    "if n <= 4: # Debug: Print all data for small inputs\n",
    "    print (\"@A_np:\\n\", A_np)\n",
    "    print (\"@x_np:\\n\", x_np)\n",
    "    print (\"@y_np:\\n\", y_np)\n",
    "    print (\"@A_py:\\n\", A_py)\n",
    "    print (\"@x_py:\\n\", x_np)\n",
    "    print (\"@y_py:\\n\", y_py)\n",
    "    print (\"@dy:\\n\", dy)\n",
    "\n",
    "# Trigger an error on likely failure\n",
    "assert dy_norm <= dy_max\n",
    "print(\"\\n(Passed!)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-f0155950b35ebcf2",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.98 s ± 11.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit matvec_py (n, n, A_py, x_py)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-3c70eea14727218a",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Fin!** If you've reached this point and everything executed without error, you can submit this part and move on to the next one."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Create Assignment",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}