尝试使用numpy向量化迭代计算
问题内容：
我试图通过使用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循环？
问题答案：
这样的线性递归可以使用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.]])
