Optimized version of cython code at http://nbviewer.ipython.org/github/rasbt/One-Python-benchmark-per-day/blob/master/ipython_nbs/day10_fortran_lstsqr.ipynb

cython_code.pyx

import numpy as np
cimport numpy as np
cimport cython
@cython.boundscheck(False) 
@cython.wraparound(False)
cpdef cython_lstsqr(x_ary, y_ary):
    """ Computes the least-squares solution to a linear matrix equation. """
    cdef double x_avg, y_avg, var_x, cov_xy,\
         slope, y_interc
    cdef double[:] x = x_ary # memory view
    cdef double[:] y = y_ary
    cdef long N
    
    N = x.shape[0]
    x_avg = np.sum(x)/N
    y_avg = np.sum(y)/N
    var_x = 0
    cov_xy = 0
    for i in range(N):
        temp = (x[i] - x_avg)
        var_x += temp**2
        cov_xy += temp*(y[i] - y_avg)
    slope = cov_xy / var_x
    y_interc = y_avg - slope*x_avg
    return (slope, y_interc)

@cython.boundscheck(False) 
@cython.wraparound(False)
@cython.cdivision(True)
cdef double cysum(double[:] x):
    cdef: 
        unsigned int i
        double s
        int N
    
    N = x.shape[0]
    
    for i in xrange(N):
        s += x[i]
        
    return s / N

@cython.boundscheck(False) 
@cython.wraparound(False)
@cython.cdivision(True)
cpdef cython_lstsqr2(x_ary, y_ary):
    """ Computes the least-squares solution to a linear matrix equation. """
    cdef double x_avg, y_avg, var_x, cov_xy,\
         slope, y_interc, temp
    cdef double[:] x = x_ary # memory view
    cdef double[:] y = y_ary
    cdef long N
    cdef int i
    
    N = x.shape[0]
    x_avg = cysum(x) #np.sum(x)/N
    y_avg = cysum(y) #np.sum(y)/N
    var_x = 0
    cov_xy = 0
    for i in range(N):
        temp = (x[i] - x_avg)
        var_x += temp**2
        cov_xy += temp*(y[i] - y_avg)
    slope = cov_xy / var_x
    y_interc = y_avg - slope*x_avg
    return (slope, y_interc)