numpy 如何加速卷积函数?

q3aa0525  于 9个月前  发布在  其他
关注(0)|答案(1)|浏览(163)

我写了一个这样的卷积函数:

import numpy as np 
import numba as nb

# Generate sample input data
num_chans = 111
num_bins = 47998 
num_rad = 8
num_col = 1000

rng = np.random.default_rng()

wvl_sensor = rng.uniform(low=1000, high=11000, size=(num_chans, num_col))
fwhm_sensor = rng.uniform(low=0.01, high=2.0, size=num_chans)

wvl_lut = rng.uniform(low=1000, high=11000, size=num_bins) 
rad_lut = rng.uniform(low=0, high=1, size=(num_rad, num_bins))

# Original convolution implementation
def original_convolve(wvl_sensor, fwhm_sensor, wvl_lut, rad_lut):

    sigma = fwhm_sensor / (2.0 * np.sqrt(2.0 * np.log(2.0)))  
    var = sigma ** 2
    denom = (2 * np.pi * var) ** 0.5
    
    numer = np.exp(-(wvl_lut[:, None] - wvl_sensor[None, :])**2 / (2*var)) 
    response = numer / denom
    
    response /= response.sum(axis=0)
    resampled = np.dot(rad_lut, response)
    
    return resampled

字符串
numpy版本运行时间约为45秒:

# numpy version
num_chans, num_col = wvl_sensor.shape
num_bins = wvl_lut.shape[0]
num_rad = rad_lut.shape[0]

original_res = np.empty((num_col, num_rad, num_chans), dtype=np.float64)

for x in range(wvl_sensor.shape[1]):
    original_res[x, :, :] = original_convolve(wvl_sensor[:, x], fwhm_sensor, wvl_lut, rad_lut)


我试着用numba加速它:

@nb.jit(nopython=True)
def numba_convolve(wvl_sensor, fwhm_sensor, wvl_lut, rad_lut):
    num_chans, num_col = wvl_sensor.shape
    num_bins = wvl_lut.shape[0]
    num_rad = rad_lut.shape[0]

    output = np.empty((num_col, num_rad, num_chans), dtype=np.float64)

    sigma = fwhm_sensor / (2.0 * np.sqrt(2.0 * np.log(2.0)))  
    var = sigma ** 2
    denom = (2 * np.pi * var) ** 0.5

    for x in nb.prange(num_col):
        numer = np.exp(-(wvl_lut[:, None] - wvl_sensor[None, :, x])**2 / (2*var))
        response = numer / denom

        response /= response.sum(axis=0)
        resampled = np.dot(rad_lut, response)
        output[x, :, :] = resampled

    return output


它仍然花费大约32秒。请注意,如果我使用@nb.jit(nopython=True, parallel=True),输出是全零值。
有什么想法可以正确应用numba?或者改进卷积函数?

t5fffqht

t5fffqht1#

您可以将Numba实现中类似Numpy的代码替换为普通循环。Numba喜欢普通循环。这种策略有助于避免创建许多小的临时数组或一些大的临时数组。这两者都会导致可扩展性问题。然后您可以并行计算外部循环(就像您所做的那样)。
这种策略还有其他好处:

  • np.dot使用一个 *BLAS库 *,它通常已经是并行的,所以在并行代码中使用它可能会导致性能问题(尽管它可以在应用程序级别进行调优);
  • 并行Numba代码中的高级Numpy特性往往更 * 虚假 *(由于Numba本身),因此使用不太高级的特性可以避免此类问题(我遇到过几次这个问题)。

下面是生成的代码:

@nb.jit(nopython=True, parallel=True)
def numba_convolve_faster(wvl_sensor, fwhm_sensor, wvl_lut, rad_lut):
    num_chans, num_col = wvl_sensor.shape
    num_bins = wvl_lut.shape[0]
    num_rad = rad_lut.shape[0]

    original_res = np.empty((num_col, num_rad, num_chans), dtype=np.float64)
    sigma = fwhm_sensor / (2.0 * np.sqrt(2.0 * np.log(2.0)))
    var = sigma ** 2
    denom = (2 * np.pi * var) ** 0.5
    inv_denom = 1.0 / denom
    factor = -1 / (2*var)

    for x in nb.prange(wvl_sensor.shape[1]):
        wvl_sensor_col = wvl_sensor[:, x].copy()
        response = np.empty(num_bins)
        for j in range(num_chans):
            response_sum = 0.0
            for i in range(num_bins):
                diff = wvl_lut[i] - wvl_sensor_col[j]
                response[i] = np.exp(diff * diff * factor[j]) * inv_denom[j]
                response_sum += response[i]
            inv_response_sum = 1.0 / response_sum
            for i in range(num_bins):
                response[i] *= inv_response_sum
            for k in range(num_rad):
                s = 0.0
                for i in range(num_bins):
                    s += rad_lut[k, i] * response[i]
                original_res[x, k, j] = s

    return original_res

res = numba_convolve_faster(wvl_sensor, fwhm_sensor, wvl_lut, rad_lut)

字符串
输出结果在我的机器上是正确的(有一个~ 1 e-16的错误)。注意,初始实现产生了像这样的NaN值。

结果

以下是我的i5- 9600 KF CPU(6核)的性能结果:

OP Numpy code:    4min 37s
OP Numba code:    4min 24s
With Numexpr:     1min 13s
This Numba code:       46s     <-----


因此,这个新的并行Numba实现比最初的实现快5.74倍。它也比使用Numexpr快。这在我的CPU上接近最佳,使用默认的指数函数。

注意事项

在x86-64 CPU上使用英特尔SVML库的指数函数,代码肯定会更快一些。
如果你有一个服务器端的GPU,那么代码可以更有效地计算(因为这种计算是GPU友好的)。主流PC GPU可能不会比CPU更快地计算,因为主流PC GPU是为简单精度而设计的。
有了simple-precision,代码在CPU和GPU上都可以明显更快,特别是如果你不需要非常高的精度。这是因为:

  • 简单精度的指数近似可以用一种非常SIMD友好的方式计算,而精确的双精度指数函数却不能(通常不值得);
  • 一些基本的简单精度运算比双精度运算更快(消耗更少的能量);
  • SIMD寄存器可以比双精度寄存器多保存两倍的简单精度数(通常每个SIMD操作的成本相同,因此吞吐量是简单精度数的两倍)。

相关问题