Question

我们想加快速度。 我们将使用 Numba，一个直接将代码编译为 C 的 Python 库。

Numba 是一个编译器。

资源

Jake VanderPlas 的 这个教程 是一个很好的介绍。 在这里，Jake 使用 Numba 实现了一个 非平凡的算法 （非均匀快速傅里叶变换）。

Cython 是另一种选择。 我发现 Cython 主要比 Numba 更多的知识（它更接近 C），但提供类似 Numba 的加速。

这里是预先编译（AOT）编译器，即时编译（JIT）编译器和解释器之间差异的 全面回答 。

使用向量化和原生代码进行实验

让我们先了解一下 Numba，然后我们将回到我们的糖尿病数据集回归的多项式特征问题。

%matplotlib inline

import math, numpy as np, matplotlib.pyplot as plt
from pandas_summary import DataFrameSummary
from scipy import ndimage

from numba import jit, vectorize, guvectorize, cuda, float32, void, float64

我们将展示以下方面的影响：

• 避免内存分配和副本（比 CPU 计算慢）
• 更好的局部性
• 向量化

如果我们一次在整个数组上使用 numpy，它会创建大量的临时值，并且不能使用缓存。 如果我们一次使用 numba 循环遍历数组项，那么我们就不必分配大型临时数组，并且可以复用缓存数据，因为我们正在对每个数组项进行多次计算。

# 无类型和没有向量化
def proc_python(xx,yy):
    zz = np.zeros(nobs, dtype='float32')
    for j in range(nobs):   
        x, y = xx[j], yy[j] 
        x = x*2 - ( y * 55 )
        y = x + y*2         
        z = x + y + 99      
        z = z * ( z - .88 ) 
        zz[j] = z           
    return zz

    nobs = 10000
x = np.random.randn(nobs).astype('float32')
y = np.random.randn(nobs).astype('float32')

%timeit proc_python(x,y) 

# 49.8 ms ± 1.19 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

NumPy

Numpy 让我们对其向量化：

# 有类型和向量化
def proc_numpy(x,y):
    z = np.zeros(nobs, dtype='float32')
    x = x*2 - ( y * 55 )
    y = x + y*2         
    z = x + y + 99      
    z = z * ( z - .88 ) 
    return z

np.allclose( proc_numpy(x,y), proc_python(x,y), atol=1e-4 )

# True

%timeit proc_numpy(x,y)    # Typed and vectorized

# 35.9 µs ± 166 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Numba

Numba 提供几种不同的装饰器。 我们将尝试两种不同的方法：

• @jit ：非常一般
• @vectorize ：不需要编写 for 循环。操作相同大小的向量时很有用

首先，我们将使用 Numba 的 jit （即时）编译器装饰器，而无需显式向量化。 这避免了大量内存分配，因此我们有更好的局部性：

@jit()
def proc_numba(xx,yy,zz):
    for j in range(nobs):   
        x, y = xx[j], yy[j] 
        x = x*2 - ( y * 55 )
        y = x + y*2         
        z = x + y + 99      
        z = z * ( z - .88 ) 
        zz[j] = z           
    return zz

z = np.zeros(nobs).astype('float32')
np.allclose( proc_numpy(x,y), proc_numba(x,y,z), atol=1e-4 )

# True

%timeit proc_numba(x,y,z)

# 6.4 µs ± 17.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

现在我们将使用 Numba 的 vectorize 装饰器。 Numba 的编译器以比普通 Python 和 Numpy 更聪明的方式优化它。 它为你写了一个 Numpy ufunc ，传统上它涉及编写 C 并且不那么简单。

@vectorize
def vec_numba(x,y):
    x = x*2 - ( y * 55 )
    y = x + y*2         
    z = x + y + 99      
    return z * ( z - .88 ) 

np.allclose(vec_numba(x,y), proc_numba(x,y,z), atol=1e-4 )

# True

%timeit vec_numba(x,y)

# 5.82 µs ± 14.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Numba 很棒。 看看这有多快！

Numba 多项式特征

@jit(nopython=True)
def vec_poly(x, res):
    m,n=x.shape
    feat_idx=0
    for i in range(n):
        v1=x[:,i]
        for k in range(m): res[k,feat_idx] = v1[k]
        feat_idx+=1
        for j in range(i,n):
            for k in range(m): res[k,feat_idx] = v1[k]*x[k,j]
            feat_idx+=1

行序和列序存储

来自 Eli Bendersky 的博客文章 ：

“矩阵的行序布局将第一行放在连续的内存中，然后是第二行放在它后面，然后是第三行，依此类推。列序布局将第一列放在连续内存中，然后放入第二列，等等....虽然知道特定数据集使用哪种布局对于良好的性能至关重要，但对于哪种布局“更好”的问题，没有单一的答案。”

“事实证明，匹配算法与数据布局的工作方式，可以决定应用程序的性能。”

“简短的说法是：始终按照布局顺序遍历数据。”

列序布局：Fortran，Matlab，R 和 Julia

行序布局：C，C ++，Python，Pascal，Mathematica

trn = np.asfortranarray(trn)
test = np.asfortranarray(test)

m,n=trn.shape
n_feat = n*(n+1)//2 + n
trn_feat = np.zeros((m,n_feat), order='F')
test_feat = np.zeros((len(y_test), n_feat), order='F')

vec_poly(trn, trn_feat)
vec_poly(test, test_feat)

regr.fit(trn_feat, y_trn)

# LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

regr_metrics(y_test, regr.predict(test_feat))

# (55.74734592292935, 42.836164292252306)

%timeit vec_poly(trn, trn_feat)

# 7.33 µs ± 19.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

回想一下，这是 sklearn PolynomialFeatures 实现的时间，它是由专家创建的：

%timeit poly.fit_transform(trn)

# 635 µs ± 9.25 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

605/7.7

# 78.57142857142857

这是一个大问题！ Numba 太神奇了！ 只需一行代码，我们就可以获得比 scikit 学习快 78 倍的速度（由专家优化）。

正则化和噪声

正则化是一种减少过拟合，并创建更好地泛化到新数据的模型的方法。

正则化

Lasso 回归使用 L1 惩罚，产生稀疏系数。 参数 α 用于加权惩罚项。 Scikit Learn 的 LassoCV 使用许多不同的 α 值进行交叉验证。

观看 Lasso 回归的 Coursera 视频 ，了解更多信息。

reg_regr = linear_model.LassoCV(n_alphas=10)

reg_regr.fit(trn_feat, y_trn)

'''
/home/jhoward/anaconda3/lib/python3.6/site-packages/sklearn/linear_model/coordinate_descent.py:484: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Fitting data with very small alpha may cause precision problems.
  ConvergenceWarning)

LassoCV(alphas=None, copy_X=True, cv=None, eps=0.001, fit_intercept=True,
    max_iter=1000, n_alphas=10, n_jobs=1, normalize=False, positive=False,
    precompute='auto', random_state=None, selection='cyclic', tol=0.0001,
    verbose=False)
'''

reg_regr.alpha_

# 0.0098199431661591518

regr_metrics(y_test, reg_regr.predict(test_feat))

# (50.0982471642817, 40.065199085003101)

噪声

现在我们将为数据添加一些噪音。

idxs = np.random.randint(0, len(trn), 10)

y_trn2 = np.copy(y_trn)
y_trn2[idxs] *= 10 # label noise

regr = linear_model.LinearRegression()
regr.fit(trn, y_trn)
regr_metrics(y_test, regr.predict(test))

# (51.1766253181518, 41.415992803872754)

regr.fit(trn, y_trn2)
regr_metrics(y_test, regr.predict(test))

# (62.66110319520415, 53.21914420254862)

Huber 损失是一种损失函数，对异常值的敏感度低于平方误差损失。 对于小的误差值，它是二次的，对于大的值，它是线性的。

hregr = linear_model.HuberRegressor()
hregr.fit(trn, y_trn2)
regr_metrics(y_test, hregr.predict(test))

# (51.24055602541746, 41.670840571376822)