Cython与Python实现相同卷积函数返回结果不一致原因排查

问题描述 #

我分别在Python和Cython中编写了逻辑完全一致的卷积计算函数，但二者运行后返回的计算结果存在明显差异，其中Python版本输出的结果准确度更高，请问造成该现象的可能原因是什么？

代码实现 #

Python版本 #

def convolution(a,b):
    n = a.size
    rv = np.zeros(n, dtype=np.float64)
    # prods = np.zeros(n,dtype=hsdp.TYPE)
    for i in range(n) :
        prods_sum = 0
        for j in range(n) :
            prods_sum += a[j] * b[(i-j)%n]
        rv[i] = prods_sum#np.sum(prods)

    return rv

Cython版本 #

def cy_convolution(np.float64_t[:] a, np.float64_t[:] b): 
    cdef int n,i,j,ix
    cdef np.float64_t prods_sum
    cdef np.float64_t[:] rv_view

    n = a.size
    rv = np.zeros(n, dtype=np.float64)
    rv_view  = rv
    for i in range(n) :
        prods_sum = 0
        for j in range(n) :
            ix = (i-j) % n
            prods_sum += a[j] * b[ix]
        rv_view[i] = prods_sum

    return rv

运行结果对比 #

In [9]:  cy_convolution(a,b)                                                                                                                                                 
Out[9]: array([ 0.4 , -0.06, -0.18,  0.42,  0.04, -0.31, -0.27, -0.23,  0.19, -0.03])

In [10]: convolution(a,b)                                                                                                                                               
Out[10]: array([ 0.45,  0.11, -0.29,  0.33,  0.02, -0.25, -0.29, -0.24,  0.18, -0.03])

In [11]: a                                                                                                                                                                   
Out[11]: [-0.57  0.4  -0.28 -0.31 -0.07 -0.03  0.53 -0.17 -0.11 -0.04]

In [12]: b                                                                                                                                                                   
Out[12]: [-0.69 -0.38  0.39  0.11  0.1   0.44  0.01 -0.05 -0.01  0.1 ]

In [14]: a.dtype                                                                                                                                                             
Out[14]: dtype('float64')

In [15]: b.dtype                                                                                                                                                             
Out[15]: dtype('float64')

Cython编译配置 #

# cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False,nonecheck=False,cdivision=True
# distutils: define_macros=NPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION

问题根因 #

该差异不是浮点精度误差，是三个编译优化选项叠加导致索引计算错误、越界读取非法内存值：

• 开启cdivision=True后，Cython的%运算符遵循C语言有符号整数取余规则，和Python取模逻辑不一致。Python取模结果符号永远和除数一致，负数取模会返回正的合法索引，例如(-1)%10 = 9；C语言有符号整数取余结果符号和被除数一致，负数取余会返回负值，例如-1%10 = -1。
• 开启wraparound=False后，Cython不再自动将负索引转换为从数组末尾计数的正索引（该行为是Python/Numpy默认的索引逻辑），负索引会直接按内存偏移计算地址。
• 开启boundscheck=False后，数组访问的边界检查被关闭，访问负索引时不会抛出越界错误，会直接读取数组起始内存地址之前的未知值，最终导致计算结果出现大幅偏差。

修复方案 #

将索引计算修改为永远非负的形式，绕开C语言负数取余的问题即可，不需要调整编译优化选项就能得到和Python版本完全一致的结果，同时保留速度优势：

ix = (i - j + n) % n

内容的提问来源于stack exchange，提问作者sten