尝试使用numpy对迭代计算进行矢量化

Question

我试图通过在numpy中使用矢量化形式来使一些代码更有效.让我举个例子让你知道我的意思.

给出以下代码:

a = np.zeros([4,4])
a[0] = [1., 2., 3., 4.]
for i in range(len(a)-1):
    a[i+1] = 2*a[i]
print a

它输出

[[  1.   2.   3.   4.]
 [  2.   4.   6.   8.]
 [  4.   8.  12.  16.]
 [  8.  16.  24.  32.]]

当我现在尝试像这样向量化代码:

a = np.zeros([4,4])
a[0] = [1., 2., 3., 4.]
a[1:] = 2*a[0:-1]
print a

我只是让第一次迭代正确:

[[ 1.  2.  3.  4.]
 [ 2.  4.  6.  8.]
 [ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]]

是否有可能以矢量化形式有效地编写上面的代码(下一次迭代总是访问上一次迭代)或者我是否必须保持for循环？

Answer 1

可以使用scipy.signal.lfilter以下公式计算此类线性递归:

In [19]: from scipy.signal import lfilter

In [20]: num = np.array([1.0])

In [21]: alpha = 2.0

In [22]: den = np.array([1.0, -alpha])

In [23]: a = np.zeros((4,4))

In [24]: a[0,:] = [1,2,3,4]

In [25]: lfilter(num, den, a, axis=0)
Out[25]: 
array([[  1.,   2.,   3.,   4.],
       [  2.,   4.,   6.,   8.],
       [  4.,   8.,  12.,  16.],
       [  8.,  16.,  24.,  32.]])

详情请参见以下内容:蟒蛇递归矢量与时间序列,在熊猫递归定义

请注意,lfilter如果您正在解决非均匀问题(例如x[i+1] = alpha*x[i] + u[i],u给定输入数组的位置),那么使用才真正有意义.对于简单的重复a[i+1] = alpha*a[i],您可以使用确切的解决方案a[i] = a[0]*alpha**i.可以使用广播对多个初始值的解决方案进行矢量化.例如,

In [271]: alpha = 2.0

In [272]: a0 = np.array([1, 2, 3, 4])

In [273]: n = 5

In [274]: a0 * (alpha**np.arange(n).reshape(-1, 1))
Out[274]: 
array([[  1.,   2.,   3.,   4.],
       [  2.,   4.,   6.,   8.],
       [  4.,   8.,  12.,  16.],
       [  8.,  16.,  24.,  32.],
       [ 16.,  32.,  48.,  64.]])

Answer 2

Numpy 的向量计算作用于向量，而不是一系列步骤，因此您必须对整个表达式进行向量化。例如：

np.multiply(np.arange(1,5), 2**np.arange(0,4)[np.newaxis].T)

要解决“最终”问题，是的，for如果要进行顺序计算，则必须保持循环。您可能会通过map或[... for ...]但以这种方式进行优化需要大量反复试验才能提高效率。用向量术语思考并使用 Numpy 来实现的美妙之处在于，您无需反复试验即可有效地获得结果。

该cumsum和cumprod功能可以做类似于你在问什么东西。取而代之的是2**np.arange(...)，您可以从

np.multiply(np.arange(1,5), np.cumprod([1,2,2,2,])[np.newaxis].T)