Homework 09

Write code to solve all problems. The grading rubric includes the following criteria:

• Correctness
• Readability
• Efficiency

For this exercise, the most important grading criteria is how much your optimizations improved run-times.

Please do not copy answers found on the web or elsewhere as it will not benefit your learning. Searching the web for general references etc. is OK. Some discussion with friends is fine too - but again, do not just copy their answer.

Honor Code: By submitting this assignment, you certify that this is your original work.

Exercise 1 (100 points)

The code given below performs a stochastic gradient descent to fit a quadratic polynomila to \(n\) data points. Maake the code run faster by:

1. Using numba JIT (20 points)
2. using Cython (30 poits)
3. Rewrite the sgd function in C or C++ and wrap for use in Python (50 points)

Replace the code stubs below with your otpimized code. Reprot the ratio optimized_time/original_time for each of the four parts.

In [45]:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

In [46]:

def sgd(b, x, y, max_iter, alpha):
    n = x.shape[0]
    for i in range(max_iter):
        for j in range(n):
            b[0] -= alpha * (2*(b[0] + b[1]*x[j] + b[2]*x[j]**2 - y[j]))
            b[1] =- alpha * (2*x[j] * (b[0] + b[1]*x[j] + b[2]*x[j]**2 - y[j]))
            b[2] -= alpha * (2*x[j]**2 * (b[0] + b[1]*x[j] + b[2]*x[j]**2 - y[j]))
    return b

In [47]:

np.random.seed(12345)
n = 10000
x = np.linspace(0, 10, n)
y = 2*x**2 + 6*x + 3 + np.random.normal(0, 5, n)
k = 100
alpha = 0.00001
b0 = np.random.random(3)

In [48]:

%%time
np.random.seed(123)
b = sgd(b0, x, y, k, alpha)
print(b)

[  9.67934503e+00   7.71388025e-04   2.53519494e+00]
CPU times: user 5.92 s, sys: 0 ns, total: 5.92 s
Wall time: 5.94 s

In [49]:

yhat = b[0] + b[1]*x+ b[2]*x**2
idx = sorted(np.random.choice(n, 100))
plt.scatter(x[idx], y[idx])
plt.plot(x[idx], yhat[idx], c='red')
pass

/opt/conda/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):

Using numba JIT

In [17]:

from numba import jit

In [50]:

@jit
def sgd_numba(b, x, y, max_iter, alpha):
    n = x.shape[0]
    for i in range(max_iter):
        for j in range(n):
            b[0] -= alpha * (2*(b[0] + b[1]*x[j] + b[2]*x[j]**2 - y[j]))
            b[1] =- alpha * (2*x[j] * (b[0] + b[1]*x[j] + b[2]*x[j]**2 - y[j]))
            b[2] -= alpha * (2*x[j]**2 * (b[0] + b[1]*x[j] + b[2]*x[j]**2 - y[j]))
    return b

In [51]:

%%time
np.random.seed(123)
b = sgd_numba(b0, x, y, k, alpha)
print(b)

[  9.82084322e+00   7.71416420e-04   2.53377882e+00]
CPU times: user 223 ms, sys: 3.77 ms, total: 227 ms
Wall time: 228 ms

In [52]:

yhat = b[0] + b[1]*x+ b[2]*x**2
idx = sorted(np.random.choice(n, 100))
plt.scatter(x[idx], y[idx])
plt.plot(x[idx], yhat[idx], c='red')
pass

/opt/conda/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):

Speed-up ratio

In [ ]:

Using Cython

In [53]:

%load_ext cython

The cython extension is already loaded. To reload it, use:
  %reload_ext cython

In [56]:

%%cython -a

import numpy as np
cimport numpy as np
import cython
from libc.math cimport pow

@cython.boundscheck(False)
@cython.wraparound(False)
def sgd_cython(double[:] b, double[:] x, double[:] y, int max_iter, double alpha):
    cdef int i, j
    cdef int n = x.shape[0]
    for i in range(max_iter):
        for j in range(n):
            b[0] -= alpha * (2*(b[0] + b[1]*x[j] + b[2]*pow(x[j], 2) - y[j]))
            b[1] =- alpha * (2*x[j] * (b[0] + b[1]*x[j] + b[2]*pow(x[j], 2) - y[j]))
            b[2] -= alpha * (2*pow(x[j], 2) * (b[0] + b[1]*x[j] + b[2]*pow(x[j], 2) - y[j]))
    return np.array(b)

/opt/conda/lib/python3.4/site-packages/IPython/utils/path.py:264: UserWarning: get_ipython_cache_dir has moved to the IPython.paths module
  warn("get_ipython_cache_dir has moved to the IPython.paths module")

Out[56]:

Generated by Cython 0.22.1

Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.

 01: 

+02: import numpy as np

  __pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, -1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 2; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 2; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
/* … */
  __pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 2; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 2; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

 03: cimport numpy as np

 04: import cython

 05: from libc.math cimport pow

 06: 

 07: @cython.boundscheck(False)

 08: @cython.wraparound(False)

+09: def sgd_cython(double[:] b, double[:] x, double[:] y, int max_iter, double alpha):

/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_e085510efc29f001ac6c143701661c75_1sgd_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static PyMethodDef __pyx_mdef_46_cython_magic_e085510efc29f001ac6c143701661c75_1sgd_cython = {"sgd_cython", (PyCFunction)__pyx_pw_46_cython_magic_e085510efc29f001ac6c143701661c75_1sgd_cython, METH_VARARGS|METH_KEYWORDS, 0};
static PyObject *__pyx_pw_46_cython_magic_e085510efc29f001ac6c143701661c75_1sgd_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
  __Pyx_memviewslice __pyx_v_b = { 0, 0, { 0 }, { 0 }, { 0 } };
  __Pyx_memviewslice __pyx_v_x = { 0, 0, { 0 }, { 0 }, { 0 } };
  __Pyx_memviewslice __pyx_v_y = { 0, 0, { 0 }, { 0 }, { 0 } };
  int __pyx_v_max_iter;
  double __pyx_v_alpha;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("sgd_cython (wrapper)", 0);
  {
    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_b,&__pyx_n_s_x,&__pyx_n_s_y,&__pyx_n_s_max_iter,&__pyx_n_s_alpha,0};
    PyObject* values[5] = {0,0,0,0,0};
    if (unlikely(__pyx_kwds)) {
      Py_ssize_t kw_args;
      const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
      switch (pos_args) {
        case  5: values[4] = PyTuple_GET_ITEM(__pyx_args, 4);
        case  4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = PyDict_Size(__pyx_kwds);
      switch (pos_args) {
        case  0:
        if (likely((values[0] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_b)) != 0)) kw_args--;
        else goto __pyx_L5_argtuple_error;
        case  1:
        if (likely((values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_x)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("sgd_cython", 1, 5, 5, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
        }
        case  2:
        if (likely((values[2] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_y)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("sgd_cython", 1, 5, 5, 2); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
        }
        case  3:
        if (likely((values[3] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_max_iter)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("sgd_cython", 1, 5, 5, 3); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
        }
        case  4:
        if (likely((values[4] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_alpha)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("sgd_cython", 1, 5, 5, 4); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
        }
      }
      if (unlikely(kw_args > 0)) {
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "sgd_cython") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
      }
    } else if (PyTuple_GET_SIZE(__pyx_args) != 5) {
      goto __pyx_L5_argtuple_error;
    } else {
      values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
      values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
      values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
      values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
      values[4] = PyTuple_GET_ITEM(__pyx_args, 4);
    }
    __pyx_v_b = __Pyx_PyObject_to_MemoryviewSlice_ds_double(values[0]); if (unlikely(!__pyx_v_b.memview)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
    __pyx_v_x = __Pyx_PyObject_to_MemoryviewSlice_ds_double(values[1]); if (unlikely(!__pyx_v_x.memview)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
    __pyx_v_y = __Pyx_PyObject_to_MemoryviewSlice_ds_double(values[2]); if (unlikely(!__pyx_v_y.memview)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
    __pyx_v_max_iter = __Pyx_PyInt_As_int(values[3]); if (unlikely((__pyx_v_max_iter == (int)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
    __pyx_v_alpha = __pyx_PyFloat_AsDouble(values[4]); if (unlikely((__pyx_v_alpha == (double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("sgd_cython", 1, 5, 5, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_e085510efc29f001ac6c143701661c75.sgd_cython", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_46_cython_magic_e085510efc29f001ac6c143701661c75_sgd_cython(__pyx_self, __pyx_v_b, __pyx_v_x, __pyx_v_y, __pyx_v_max_iter, __pyx_v_alpha);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_e085510efc29f001ac6c143701661c75_sgd_cython(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_b, __Pyx_memviewslice __pyx_v_x, __Pyx_memviewslice __pyx_v_y, int __pyx_v_max_iter, double __pyx_v_alpha) {
  CYTHON_UNUSED int __pyx_v_i;
  int __pyx_v_j;
  int __pyx_v_n;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("sgd_cython", 0);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_28);
  __Pyx_XDECREF(__pyx_t_29);
  __Pyx_XDECREF(__pyx_t_30);
  __Pyx_XDECREF(__pyx_t_31);
  __Pyx_XDECREF(__pyx_t_32);
  __Pyx_AddTraceback("_cython_magic_e085510efc29f001ac6c143701661c75.sgd_cython", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __PYX_XDEC_MEMVIEW(&__pyx_v_b, 1);
  __PYX_XDEC_MEMVIEW(&__pyx_v_x, 1);
  __PYX_XDEC_MEMVIEW(&__pyx_v_y, 1);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple__19 = PyTuple_Pack(8, __pyx_n_s_b, __pyx_n_s_x, __pyx_n_s_y, __pyx_n_s_max_iter, __pyx_n_s_alpha, __pyx_n_s_i, __pyx_n_s_j, __pyx_n_s_n); if (unlikely(!__pyx_tuple__19)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_tuple__19);
  __Pyx_GIVEREF(__pyx_tuple__19);
/* … */
  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_e085510efc29f001ac6c143701661c75_1sgd_cython, NULL, __pyx_n_s_cython_magic_e085510efc29f001ac); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_sgd_cython, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_codeobj__20 = (PyObject*)__Pyx_PyCode_New(5, 0, 8, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__19, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_home_jovyan_cache_ipython_cytho, __pyx_n_s_sgd_cython, 9, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__20)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 9; __pyx_clineno = __LINE__; goto __pyx_L1_error;}

 10:     cdef int i, j

+11:     cdef int n = x.shape[0]

  __pyx_v_n = (__pyx_v_x.shape[0]);

+12:     for i in range(max_iter):

  __pyx_t_1 = __pyx_v_max_iter;
  for (__pyx_t_2 = 0; __pyx_t_2 < __pyx_t_1; __pyx_t_2+=1) {
    __pyx_v_i = __pyx_t_2;

+13:         for j in range(n):

    __pyx_t_3 = __pyx_v_n;
    for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) {
      __pyx_v_j = __pyx_t_4;

+14:             b[0] -= alpha * (2*(b[0] + b[1]*x[j] + b[2]*pow(x[j], 2) - y[j]))

      __pyx_t_5 = 0;
      __pyx_t_6 = 1;
      __pyx_t_7 = __pyx_v_j;
      __pyx_t_8 = 2;
      __pyx_t_9 = __pyx_v_j;
      __pyx_t_10 = __pyx_v_j;
      __pyx_t_11 = 0;
      *((double *) ( /* dim=0 */ (__pyx_v_b.data + __pyx_t_11 * __pyx_v_b.strides[0]) )) -= (__pyx_v_alpha * (2.0 * ((((*((double *) ( /* dim=0 */ (__pyx_v_b.data + __pyx_t_5 * __pyx_v_b.strides[0]) ))) + ((*((double *) ( /* dim=0 */ (__pyx_v_b.data + __pyx_t_6 * __pyx_v_b.strides[0]) ))) * (*((double *) ( /* dim=0 */ (__pyx_v_x.data + __pyx_t_7 * __pyx_v_x.strides[0]) ))))) + ((*((double *) ( /* dim=0 */ (__pyx_v_b.data + __pyx_t_8 * __pyx_v_b.strides[0]) ))) * pow((*((double *) ( /* dim=0 */ (__pyx_v_x.data + __pyx_t_9 * __pyx_v_x.strides[0]) ))), 2.0))) - (*((double *) ( /* dim=0 */ (__pyx_v_y.data + __pyx_t_10 * __pyx_v_y.strides[0]) ))))));

+15:             b[1] =- alpha * (2*x[j] * (b[0] + b[1]*x[j] + b[2]*pow(x[j], 2) - y[j]))

      __pyx_t_12 = __pyx_v_j;
      __pyx_t_13 = 0;
      __pyx_t_14 = 1;
      __pyx_t_15 = __pyx_v_j;
      __pyx_t_16 = 2;
      __pyx_t_17 = __pyx_v_j;
      __pyx_t_18 = __pyx_v_j;
      __pyx_t_19 = 1;
      *((double *) ( /* dim=0 */ (__pyx_v_b.data + __pyx_t_19 * __pyx_v_b.strides[0]) )) = ((-__pyx_v_alpha) * ((2.0 * (*((double *) ( /* dim=0 */ (__pyx_v_x.data + __pyx_t_12 * __pyx_v_x.strides[0]) )))) * ((((*((double *) ( /* dim=0 */ (__pyx_v_b.data + __pyx_t_13 * __pyx_v_b.strides[0]) ))) + ((*((double *) ( /* dim=0 */ (__pyx_v_b.data + __pyx_t_14 * __pyx_v_b.strides[0]) ))) * (*((double *) ( /* dim=0 */ (__pyx_v_x.data + __pyx_t_15 * __pyx_v_x.strides[0]) ))))) + ((*((double *) ( /* dim=0 */ (__pyx_v_b.data + __pyx_t_16 * __pyx_v_b.strides[0]) ))) * pow((*((double *) ( /* dim=0 */ (__pyx_v_x.data + __pyx_t_17 * __pyx_v_x.strides[0]) ))), 2.0))) - (*((double *) ( /* dim=0 */ (__pyx_v_y.data + __pyx_t_18 * __pyx_v_y.strides[0]) ))))));

+16:             b[2] -= alpha * (2*pow(x[j], 2) * (b[0] + b[1]*x[j] + b[2]*pow(x[j], 2) - y[j]))

      __pyx_t_20 = __pyx_v_j;
      __pyx_t_21 = 0;
      __pyx_t_22 = 1;
      __pyx_t_23 = __pyx_v_j;
      __pyx_t_24 = 2;
      __pyx_t_25 = __pyx_v_j;
      __pyx_t_26 = __pyx_v_j;
      __pyx_t_27 = 2;
      *((double *) ( /* dim=0 */ (__pyx_v_b.data + __pyx_t_27 * __pyx_v_b.strides[0]) )) -= (__pyx_v_alpha * ((2.0 * pow((*((double *) ( /* dim=0 */ (__pyx_v_x.data + __pyx_t_20 * __pyx_v_x.strides[0]) ))), 2.0)) * ((((*((double *) ( /* dim=0 */ (__pyx_v_b.data + __pyx_t_21 * __pyx_v_b.strides[0]) ))) + ((*((double *) ( /* dim=0 */ (__pyx_v_b.data + __pyx_t_22 * __pyx_v_b.strides[0]) ))) * (*((double *) ( /* dim=0 */ (__pyx_v_x.data + __pyx_t_23 * __pyx_v_x.strides[0]) ))))) + ((*((double *) ( /* dim=0 */ (__pyx_v_b.data + __pyx_t_24 * __pyx_v_b.strides[0]) ))) * pow((*((double *) ( /* dim=0 */ (__pyx_v_x.data + __pyx_t_25 * __pyx_v_x.strides[0]) ))), 2.0))) - (*((double *) ( /* dim=0 */ (__pyx_v_y.data + __pyx_t_26 * __pyx_v_y.strides[0]) ))))));
    }
  }

+17:     return np.array(b)

  __Pyx_XDECREF(__pyx_r);
  __pyx_t_29 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_29)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_29);
  __pyx_t_30 = __Pyx_PyObject_GetAttrStr(__pyx_t_29, __pyx_n_s_array); if (unlikely(!__pyx_t_30)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_30);
  __Pyx_DECREF(__pyx_t_29); __pyx_t_29 = 0;
  __pyx_t_29 = __pyx_memoryview_fromslice(__pyx_v_b, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_29)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_29);
  __pyx_t_31 = NULL;
  if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_30))) {
    __pyx_t_31 = PyMethod_GET_SELF(__pyx_t_30);
    if (likely(__pyx_t_31)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_30);
      __Pyx_INCREF(__pyx_t_31);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_30, function);
    }
  }
  if (!__pyx_t_31) {
    __pyx_t_28 = __Pyx_PyObject_CallOneArg(__pyx_t_30, __pyx_t_29); if (unlikely(!__pyx_t_28)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    __Pyx_DECREF(__pyx_t_29); __pyx_t_29 = 0;
    __Pyx_GOTREF(__pyx_t_28);
  } else {
    __pyx_t_32 = PyTuple_New(1+1); if (unlikely(!__pyx_t_32)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    __Pyx_GOTREF(__pyx_t_32);
    __Pyx_GIVEREF(__pyx_t_31); PyTuple_SET_ITEM(__pyx_t_32, 0, __pyx_t_31); __pyx_t_31 = NULL;
    __Pyx_GIVEREF(__pyx_t_29);
    PyTuple_SET_ITEM(__pyx_t_32, 0+1, __pyx_t_29);
    __pyx_t_29 = 0;
    __pyx_t_28 = __Pyx_PyObject_Call(__pyx_t_30, __pyx_t_32, NULL); if (unlikely(!__pyx_t_28)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
    __Pyx_GOTREF(__pyx_t_28);
    __Pyx_DECREF(__pyx_t_32); __pyx_t_32 = 0;
  }
  __Pyx_DECREF(__pyx_t_30); __pyx_t_30 = 0;
  __pyx_r = __pyx_t_28;
  __pyx_t_28 = 0;
  goto __pyx_L0;

In [57]:

%%time
np.random.seed(123)
b = sgd_cython(b0, x, y, k, alpha)
print(b)

[  9.82306353e+00   7.71416866e-04   2.53375660e+00]
CPU times: user 28.7 ms, sys: 166 µs, total: 28.9 ms
Wall time: 28.5 ms

In [58]:

yhat = b[0] + b[1]*x+ b[2]*x**2
idx = sorted(np.random.choice(n, 100))
plt.scatter(x[idx], y[idx])
plt.plot(x[idx], yhat[idx], c='red')
pass

/opt/conda/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):

Speed-up ratio

In [ ]:

Using C or C++

In [59]:

%%file sgd.h

#pragma once
void sgd(double *b, double *x, double *y, int max_iter, double alpha, int n);

Overwriting sgd.h

In [60]:

%%file sgd.c
#include "sgd.h"
#include <math.h>

void sgd(double *b, double *x, double *y, int max_iter, double alpha, int n) {
    for (int i=0; i<max_iter; i++) {
        for (int j=0; j<n; j++) {
            b[0] -= alpha * (2*(b[0] + b[1]*x[j] + b[2]*pow(x[j], 2) - y[j]));
            b[1] =- alpha * (2*x[j] * (b[0] + b[1]*x[j] + b[2]*pow(x[j], 2) - y[j]));
            b[2] -= alpha * (2*pow(x[j], 2) * (b[0] + b[1]*x[j] + b[2]*pow(x[j], 2) - y[j]));
        }
    }
}

Overwriting sgd.c

In [61]:

%%file sgd_wrap.pyx

cdef extern from 'sgd.h':
    void sgd(double *b, double *x, double *y, int max_iter, double alpha, int n)

def sgd_wrap(double[::1] b0, double[::1] x, double[::1] y, int max_iter, double alpha):
    sgd(&b0[0], &x[0], &y[0], max_iter, alpha, len(y))
    return b0

Overwriting sgd_wrap.pyx

In [62]:

%%file setup.py
from distutils.core import setup, Extension
from Cython.Build import cythonize

ext = Extension("sgd_wrap",
                sources=["sgd_wrap.pyx", "sgd.c"],
                libraries=["m"],
                extra_compile_args=["-std=c99"])

setup(
    ext_modules = cythonize(
            ext
    ))

Overwriting setup.py

In [63]:

! python setup.py build_ext --inplace

Compiling sgd_wrap.pyx because it changed.
Cythonizing sgd_wrap.pyx
running build_ext
building 'sgd_wrap' extension
gcc -pthread -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/opt/conda/include/python3.4m -c sgd_wrap.c -o build/temp.linux-x86_64-3.4/sgd_wrap.o -std=c99
gcc -pthread -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/opt/conda/include/python3.4m -c sgd.c -o build/temp.linux-x86_64-3.4/sgd.o -std=c99
gcc -pthread -shared -L/opt/conda/lib -Wl,-rpath=/opt/conda/lib,--no-as-needed build/temp.linux-x86_64-3.4/sgd_wrap.o build/temp.linux-x86_64-3.4/sgd.o -L/opt/conda/lib -lm -lpython3.4m -o /home/jovyan/work/sta-663-2016/homework/sgd_wrap.cpython-34m.so

In [64]:

from sgd_wrap import sgd_wrap

In [65]:

%%time
np.random.seed(123)
b = sgd_wrap(b0, x, y, k, alpha)
print(np.array(b))

[  9.82309837e+00   7.71416873e-04   2.53375625e+00]
CPU times: user 32.2 ms, sys: 34 µs, total: 32.2 ms
Wall time: 35.1 ms

In [66]:

yhat = b[0] + b[1]*x+ b[2]*x**2
idx = sorted(np.random.choice(n, 100))
plt.scatter(x[idx], y[idx])
plt.plot(x[idx], yhat[idx], c='red')
pass

/opt/conda/lib/python3.4/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if self._edgecolors == str('face'):

Speed-up ratio

In [ ]: