numba and cython benchmarks

pairwise_benchmark.ipynb

{
 "metadata": {
  "name": "pairwise_benchmark"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Benchmarking\n",
      "\n",
      "We want to show how different tools can be used with python to get great performance.\n",
      "\n",
      "To do this analysis we will compute the pairwise distance between each point in a set of points.\n",
      "\n",
      "Below we have included 5 implementations of this algorithm:\n",
      "\n",
      "1. python + numpy\n",
      "2. python + numpy + numba\n",
      "3. cython with array buffers for accessing slices\n",
      "4. cython with C memory views for accessing slices\n",
      "5. MatLab\n",
      "\n",
      "For each of these implementations we have done some preliminary benchmarking and we report the findings at the bottom of this page.\n",
      "\n",
      "Note that the idea for these benchmarks came from http://jakevdp.github.com\n",
      "\n",
      "Also, I recognize that the innermost loop (with k as the iterator) in each of these pieces of code could easily be vectorized, but explicitly writing it out does a better job at showing the performance differences among the cases."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext cythonmagic\n",
      "\n",
      "import pyximport as pyx\n",
      "import numpy as np\n",
      "pyx.install(setup_args={\"include_dirs\": np.get_include()})  # points cython to numpy bindings\n",
      "X = np.random.randn(800, 3)  # 800 3 dimensional vectors\n",
      "D = np.empty((800, 800))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Python version\n",
      "def pairs_py(X, dist):\n",
      "    \"\"\"\n",
      "    Compute the pairwise distance between all points in the matrix X.\n",
      "    Each row of X is treated as an n-tuple of coordinates specifying\n",
      "    an n-dimensional vector.\n",
      "\n",
      "    The distance is the Euclidean (L2) norm.\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    X: array-like, dtype=float, shape=(m x n)\n",
      "        The array X for which the pairwise distance is to be computed\n",
      "\n",
      "    dist: array-like, dtype=float\n",
      "        The empty (m x m) matrix to be filled with the distances.\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    dist: array-like, dtype=float\n",
      "        The (m x m) matrix representing the pairwise distance of each\n",
      "        row in X.\n",
      "    \"\"\"\n",
      "    M = X.shape[0]\n",
      "    N = X.shape[1]\n",
      "    for i in range(M):\n",
      "        for j in range(M):\n",
      "            d = 0.0\n",
      "            for k in range(N):\n",
      "                tmp = X[i, k] - X[j, k]\n",
      "                d += tmp * tmp\n",
      "            dist[i, j] = np.sqrt(d)\n",
      "\n",
      "    return dist"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "timeit pairs_py(X, D)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 7.24 s per loop\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Numba version\n",
      "from numba import autojit\n",
      "\n",
      "@autojit()\n",
      "def pairs_numba(X ,dist):\n",
      "    M = X.shape[0]\n",
      "    N = X.shape[1]\n",
      "    for i in range(M):\n",
      "        for j in range(M):\n",
      "            d = 0.0\n",
      "            for k in range(N):\n",
      "                tmp = X[i, k] - X[j, k]\n",
      "                d += tmp * tmp\n",
      "            dist[i, j] = np.sqrt(d)\n",
      "\n",
      "    return dist"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "timeit pairs_numba(X, D)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 7.4 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "\n",
      "# Cython using array buffers\n",
      "\n",
      "cimport cython\n",
      "cimport numpy as np\n",
      "from libc.math cimport sqrt\n",
      "\n",
      "@cython.boundscheck(False)\n",
      "@cython.wraparound(False)\n",
      "def pairs_cy1(np.ndarray[double, ndim=2, mode='c'] X,\n",
      "              np.ndarray[double, ndim=2, mode='c'] dist):\n",
      "    cdef int M = X.shape[0]\n",
      "    cdef int N = X.shape[1]\n",
      "    cdef double tmp, d\n",
      "    cdef int i, j, k\n",
      "    for i in range(M):\n",
      "        for j in range(M):\n",
      "            d = 0.0\n",
      "            for k in range(N):\n",
      "                tmp = X[i, k] - X[j, k]\n",
      "                d += tmp * tmp\n",
      "            dist[i, j] = sqrt(d)\n",
      "\n",
      "    return dist"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "timeit pairs_cy1(X, D)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 7.04 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "\n",
      "# Cython using C memory views on arrays\n",
      "\n",
      "cimport cython\n",
      "from libc.math cimport sqrt\n",
      "\n",
      "@cython.boundscheck(False)\n",
      "@cython.wraparound(False)\n",
      "def pairs_cy2(double[:, ::1] X, double[:, ::1] dist):\n",
      "    cdef int M = X.shape[0]\n",
      "    cdef int N = X.shape[1]\n",
      "    cdef double tmp, d\n",
      "    cdef int i, j, k\n",
      "    for i in range(M):\n",
      "        for j in range(M):\n",
      "            d = 0.0\n",
      "            for k in range(N):\n",
      "                tmp = X[i, k] - X[j, k]\n",
      "                d += tmp * tmp\n",
      "            dist[i, j] = sqrt(d)\n",
      "\n",
      "    return dist"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "timeit pairs_cy2(X, D)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 7.02 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### MatLab testing\n",
      "\n",
      "This is the function comparable to the pairs functions above:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "    function dist = Pairs(X, dist)\n",
      "    \n",
      "    %Compute the pairwise distance between all points in the matrix X\n",
      "    %Each row of X is treated as an n-tuple of coordinates specifying\n",
      "    %an n dimensional vector.\n",
      "    %----------------------------------------------------------------\n",
      "    %Parameters:\n",
      "    %X: matrix, shape = (m x n)\n",
      "    %\tThe array X for which the pairwise distance is to be computed\n",
      "    %dist: matrix, shape=(m x m)\n",
      "    %   The array to store difference data\n",
      "    \n",
      "    %Returns:\n",
      "    %dist: matrix, shape = (m x m)\n",
      "    %\t The m x m matrix representing the pairwise distance of each\n",
      "    %\t row in X.\n",
      "    \n",
      "    \n",
      "    %First we get the dimensions of X\n",
      "    [m , n] = size(X);\n",
      "    \n",
      "    for i = 1:m\n",
      "        for j = 1:m\n",
      "            d = 0.0;\n",
      "            for k = 1:n\n",
      "                tmp = X(i, k) - X(j, k);\n",
      "                d = d + tmp * tmp;\n",
      "            dist(i, j) = sqrt(d);\n",
      "            end\n",
      "        end\n",
      "    end"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "These are the commands used to do the benchmark:\n",
      "\n",
      "    X = randn(800, 3);\n",
      "    D = zeros(800, 800);\n",
      "    x = @() Pairs(X, D);\n",
      "    \n",
      "    for i=1:15\n",
      "        y(i) = timeit(x);\n",
      "    end\n",
      "    \n",
      "    mean(y)\n",
      "\n",
      "And this is the result (in seconds):\n",
      "\n",
      "    ans =\n",
      "    \n",
      "        0.2118\n",
      "\n",
      "    "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Results\n",
      "\n",
      "Just taking a ratio of computation time. Because cython1 was the fastest, I will divide all the other ones by that to get relative times:\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pandas as pd\n",
      "py = 7.24\n",
      "numba = 0.00739\n",
      "cy1 = 0.00704\n",
      "cy2 = 0.00702\n",
      "matlab = 0.2118\n",
      "data = pd.DataFrame([py, numba , cy1, cy2, matlab], \n",
      "            index=['python', 'numba', 'cython 1', 'cython 2', 'MatLab'], columns=['Raw times (s)'])\n",
      "data[\"Relative to cython 2\"] = data['Raw times (s)'] / cy2\n",
      "data"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>Raw times (s)</th>\n",
        "      <th>Relative to cython 2</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <td><strong>python</strong></td>\n",
        "      <td> 7.24000</td>\n",
        "      <td> 1031.339031</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <td><strong>numba</strong></td>\n",
        "      <td> 0.00739</td>\n",
        "      <td>    1.052707</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <td><strong>cython 1</strong></td>\n",
        "      <td> 0.00704</td>\n",
        "      <td>    1.002849</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <td><strong>cython 2</strong></td>\n",
        "      <td> 0.00702</td>\n",
        "      <td>    1.000000</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <td><strong>MatLab</strong></td>\n",
        "      <td> 0.21180</td>\n",
        "      <td>   30.170940</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "</div>"
       ],
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "          Raw times (s)  Relative to cython 2\n",
        "python          7.24000           1031.339031\n",
        "numba           0.00739              1.052707\n",
        "cython 1        0.00704              1.002849\n",
        "cython 2        0.00702              1.000000\n",
        "MatLab          0.21180             30.170940"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}