实现 numpy.sin(x) / x 的最佳方法，其中 x 可能包含 0

Question

我现在做的是：

import numpy as np

eps = np.finfo(float).eps

def sindiv(x):
    x = np.abs(x)
    return np.maximum(eps, np.sin(x)) / np.maximum(eps, x)

但是还有很多额外的数组操作。有没有更好的办法？

Answer 1

您可以使用numpy.sinc，它计算 sin(pi x)/(pi x)：

In [20]: x = 2.4

In [21]: np.sin(x)/x
Out[21]: 0.28144299189631289

In [22]: x_over_pi = x / np.pi

In [23]: np.sinc(x_over_pi)
Out[23]: 0.28144299189631289

In [24]: np.sinc(0)
Out[24]: 1.0

Answer 2

在 numpy 数组表示法中（所以你得到一个 np 数组）：

def sindiv(x):
    return np.where(np.abs(x) < 0.01, 1.0 - x*x/6.0, np.sin(x)/x)

在这里，我将“epsilon”设置得相当大以进行测试，并使用泰勒级数的前两项进行近似。在实践中，我会将 0.01 更改为您的 eps（机器 epsilon）的一些小倍数。

xx = np.arange(-0.1, 0.1, 0.001)
yy = sinxdiv(xx)
type(yy)

输出 numpy.ndarray 并且值在原点附近是连续的（并且可微分，如果这很重要）。

如果您不想要双重评估（即在上面评估两个分支），那么我认为您必须使用循环，因为我不相信有任何“懒惰的地方”选项。

def sindiv(x):
    sox = np.zeros(x.size)
    for i in xrange(x.size):
        xv = x[i]
        if np.abs(xv) < 0.001: # For testing, use a small multiple of machine epsilon
            sox[i] = 1.0 - xv * xv / 6.0
        else:
            sox[i] = np.sin(xv) / xv
    return sox

要使这个真正pythonic，最好检查x 的类型，如果它不是数组，则只执行非数组版本。

Answer 3

正如其他人所说，numpy.sinc() 是最简单的。

我想在 NumPy 1.21.2 (link) 中包含其当前实现的副本，以表明没有特殊技巧：

y = pi * where(x == 0, 1.0e-20, x)
return sin(y)/y

基本上就是sin(x)/x。请注意，在创建y 时：乘以pi、where() 和x == 0 将创建至少2 个中间数组加上y 的最终数组。然后sin(y)/y 再创建两个数组。 numpy.sinc() 总共创建了至少 5 个数组；根据我的统计，您的 sindiv() 还创建了至少 5 个数组，因此实际上并没有那么浪费。

---

这是另一个实现：

TINY = np.finfo(float).tiny  # ≈ 2e-308 (smallest 'normal' float)
def mysinc(x):
    y = np.abs(np.pi*x) + TINY
    return np.sin(y)/y

我很确定这会将 相同 值返回到 numpy.sinc()。原因是sin(x) == x 相对“大”值x：

x = np.ldexp(1, -26, dtype=np.double)  # x = 2**-26 ≈ 1.5e-8
print(np.sin(x) == x)  # True

x = np.ldexp(1, -32, dtype=np.longdouble)  # x = 2**-32 ≈ 2.3e-10
print(np.sin(x) == x)  # True

所以对于足够小的x（忽略 pi 因子），mysinc(x) = (x+TINY)/(x+TINY) = x/x = np.sinc(x)。发生这种情况的确切阈值并不重要，只要TINY <np.spacing(x) 在这种情况下发生，以便x + TINY = x 在这种情况下。

（从泰勒级数sin(x) = x - x**3/6 + ... = x(1-x**2/6) + ... 可以理解，截止值在机器ε 的平方根附近。所以TINY 总是足够小而无所谓。）

时间

import numpy as np

eps = np.finfo(float).eps
tiny = np.finfo(float).tiny

def npsinc(x):
    y = np.pi * np.where(x == 0, 1.0e-20, x)
    return np.sin(y)/y

def sindiv(x):
    x = np.pi * np.abs(x)
    return np.maximum(eps, np.sin(x)) / np.maximum(eps, x)

def mysinc(x):
    y = np.abs(np.pi*x) + tiny
    return np.sin(y)/y

def mysinc2(x):
    y = np.abs(np.pi*x)
    y += tiny  # in-place addition
    return np.sin(y)/y

# Test data
x = np.random.rand(100)
x[np.random.randint(100, size=10)] = 0

%timeit npsinc(x)
# 10.9 µs ± 18.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit sindiv(x)
# 9.4 µs ± 12.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit mysinc(x)
# 7.38 µs ± 15.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit mysinc2(x)
# 8.64 µs ± 20.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

奇怪地使用mysinc2() 和就地添加似乎更慢，使用就地numpy.abs() 和就地numpy.sin() 甚至更慢。不完全确定原因，但请参阅 this related question。

不管怎样，如果你真的需要性能，你可以尝试使用Cython生成C代码并正确地做事，而不是玩弄NumPy：

%%cython

from libc.math cimport M_PI, sin
cimport cython
cimport numpy as np
import numpy as np

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cdef _cysinc(double[:] x, double[:] out):
    cdef size_t i
    for i in range(x.shape[0]):
        if x[i] == 0:
            out[i] = 1
        else:
            out[i] = sin(M_PI*x[i])/(M_PI*x[i])

def cysinc(np.ndarray x):
    out = np.empty_like(x)
    _cysinc(x.ravel(), out.ravel())
    return out

%timeit cysinc(x)
# 4.38 µs ± 11.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

一如既往，不要过早优化，直接使用numpy.sinc()即可。

---

旁注

有一个问题Is boost::math::sinc_pi unnecessarily complicated? 询问使用关于 x=0 的泰勒展开的好处。总之，几乎没有，但也许他们这样做是出于其他原因。

需要强调的是，浮点除法或将一个小数除以一个小数并没有什么不稳定的地方，因为您只是将有效数字相除并减去指数。

如果您将sinc(x) 计算为sin(x)/x，而不是直接泰勒级数或其他总和超出机器εnp.spacing(sinc(x)) 的方法，则您最多会偏离来自回合的np.spacing(sinc(x)) -除法/ 中的关闭错误，就像乘法* 一样。 （假设没有subnormal 业务，即使在sin(x)/x 的处理中也无所谓。）

Answer 4

允许 div 为零并稍后替换 NaN 怎么样？

import numpy as np

def sindiv(x):
    a = np.sin(x)/x
    a = np.nan_to_num(a)
    return a

如果您不想要警告，请通过seterr 禁止它们

当然，使用a 可以消除：

def sindiv(x):
    return np.nan_to_num(np.sin(x)/x)