cython memoryview not faster than ndarray

Question

我有一个用常规numpy ndarray和另一个用typed memoryview写的函数。然而，我无法让memoryview版本比常规版本更快地工作(不像许多博客，如内存视图基准测试)。

任何提高内存视图代码相对于numpy替代方案的速度的指针/建议将非常感谢!…还是……如果有人能指出为什么memoryview版本没有比常规numpy版本快得多

在下面的代码中有两个函数，它们都接受两个向量bi和xi并返回一个矩阵。第一个函数shrink_correl是常规的numpy版本，第二个函数shrink_correl2是memoryview的备选方案(让文件为sh_cor.pyx)。

# cython: boundscheck=False
# cython: wraparound=False
# cython: cdivision=True
cimport cython
cimport numpy as np
import numpy as np
from numpy cimport ndarray as ar
# -- ***this is the Regular Cython version*** -
cpdef ar[double, ndim=2, mode='c'] shrink_correl(ar[double, ndim=1, mode='c'] bi, ar[double, ndim=1, mode='c'] xi):
    cdef:
        int n_ = xi.shape[0]
        int n__ = int(n_*(n_-1)/2)
        ar[double, ndim=2, mode='c'] f = np.zeros([n__, n_+1])
        int x__ = 0
        ar[double, ndim=2, mode='c'] f1 = np.zeros([n_, n_+1])
        ar[double, ndim=2, mode='c'] f2 = np.zeros([n__, n_+1])
        ar[double, ndim=1, mode='c'] g = np.zeros(n_+1)
        ar[double, ndim=1, mode='c'] s = np.zeros(n__)
        ar[double, ndim=2, mode='c'] cori_ = np.zeros([n_, n_])
        Py_ssize_t j, k
    with nogil:
        for j in range(0, n_-1):
            for k in range(j+1, n_):
                x__ += 1
                f[x__-1, j] = bi[k]*xi[k]*1000
                f[x__-1, k] = bi[j]*xi[j]*1000
    f1 = np.dot(np.transpose(f), f)      
    with nogil:
        for j in range(0, n_):
            f1[n_, j] = xi[j]*1000
            f1[j, n_] = f1[n_, j]
    f2 = np.dot(f, np.linalg.inv(f1))
    with nogil:
        for j in range(0, n_):
            g[j] = -bi[j]*xi[j]*1000
    s = np.dot(f2, g)
    with nogil:
        for j in range(0, n_):
            cori_[j, j] = 1.0
    x__ = 0
    with nogil:
        for j in range(0, n_-1):
            for k in range(j+1, n_):
                x__ += 1
                cori_[j, k] = s[x__-1]
                cori_[k, j] = cori_[j, k]
    return cori_
# -- ***this is the MemoryView Cython version*** -    
cpdef ar[double, ndim=2, mode='c'] shrink_correl2(double[:] bi, double[:] xi):
    cdef:
        int n_ = xi.shape[0]
        int n__ = int(n_*(n_-1)/2)
        double[:, ::1] f = np.zeros([n__, n_+1])
        int x__ = 0
        double[:, ::1] f1 = np.zeros([n_, n_+1])
        double[:, ::1] f2 = np.zeros([n__, n_+1])
        double[:] g = np.zeros(n_+1)
        double[:] s = np.zeros(n__)
        double[:, ::1] cori_ = np.zeros([n_, n_])
        ar[double, ndim=2, mode='c'] cori__ = np.zeros([n_, n_])
        Py_ssize_t j, k
    with nogil:
        for j in range(0, n_-1):
            for k in range(j+1, n_):
                x__ += 1
                f[x__-1, j] = bi[k]*xi[k]*1000
                f[x__-1, k] = bi[j]*xi[j]*1000
    f1 = np.dot(np.transpose(f), f)      
    with nogil:
        for j in range(0, n_):
            f1[n_, j] = xi[j]*1000
            f1[j, n_] = f1[n_, j]
    f2 = np.dot(f, np.linalg.inv(f1))
    with nogil:
        for j in range(0, n_):
            g[j] = -bi[j]*xi[j]*1000
    s = np.dot(f2, g)
    with nogil:
        for j in range(0, n_):
            cori_[j, j] = 1.0
    x__ = 0
    with nogil:
        for j in range(0, n_-1):
            for k in range(j+1, n_):
                x__ += 1
                cori_[j, k] = s[x__-1]
                cori_[k, j] = cori_[j, k]
    cori__[:, :] = cori_
    return cori__

使用以下setup.py代码编译

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy as np
import os
ext_modules = [Extension('sh_cor', ['sh_cor.pyx'], include_dirs=[np.get_include(),
                                                                 os.path.join(np.get_include(), 'numpy')],
                         define_macros=[('NPY_NO_DEPRECATED_API', None)],
                         extra_compile_args=['-O3', '-march=native', '-ffast-math', '-flto'],
                         libraries=['m']
                         )]
setup(
    name="Sh Cor",
    cmdclass={'build_ext': build_ext},
    ext_modules=ext_modules
)

用来测试速度的代码是

import numpy as np
import sh_cor  # this the library created by the setup.py file
import time
b = np.random.random(400)
b = b/np.sum(b)
x = np.random.random(400)-0.5
n = 10 
t0 = time.time()
for i in range(n):
    v1 = sh_cor.shrink_correl(b, x)
t1 = time.time()
print((t1-t0)/n)
t0 = time.time()
for i in range(n):
    v2 = sh_cor.shrink_correl2(b, x)
t1 = time.time()
print((t1-t0)/n)

我的PC上的输出是:

0.7070999860763549   # regular numpy
0.6726999998092651   # memoryview

使用memoryview(在上面的代码中)只给了我5%的速度提升(不像博客中巨大的速度提升)。

Answer 1

@uday给我大约一个星期，因为我的电脑少了，但这里是加快速度让你开始的地方:1)而不是使用np.transpose创建一个与你想要在任何循环之前转置的相同的内存视图(即你将变量f声明为不需要gil的内存视图，只是在f_t上创建一个视图，即cdef double[:, ::1] f_T = np.transpose(f)或=f.T。

2)这一步有点棘手，因为你需要一个C/c++风格的np.dot包装器版本(所以在这种情况下，确保对dgemm函数的调用是with nogil:在它上面&缩进函数的下一行以释放4个空格缩进(SO要求):https://gist.github.com/pv/5437087。这个例子看起来很有效(尽管你必须保存包含f2pyptr.h文件，并把它放在你的项目正在构建的地方;我还怀疑你应该添加cimport numpy as np);如果不是，它需要mod，你可以看到我在另一篇文章中所做的:调用BLAS/LAPACK直接使用SciPy接口和Cython(指针问题?)/-也如何添加MKL然后你需要在顶部添加from cython.parallel cimport prange，并将所有循环从range更改为prange，并确保所有prange部分都是nogil，所有变量都是cdef在操作之前声明。此外，您必须在编译器参数中将-openmp添加到setup.py中，并添加到其包含库的链接。如果你需要澄清，问更多的问题。这并不像它应该的那么容易，但有了一点指导就变得相当简单。基本上，一旦你的setup.py被修改为包含所有内容，它将继续工作。

3)虽然可能是最容易修复的-摆脱那个列表。如果需要文本和数据，将其设置为numpy数组或pandas数据框架。每当我使用列表来存储数据时，速度的减缓都是令人难以置信的。