Numpy：根据前一个元素计算？答案

【问题标题】：Numpy: calculate based on previous element?Numpy：根据前一个元素计算？
【发布时间】：2015-07-16 01:48:30
【问题描述】：

假设我有数组x 和y：

x = numpy.array([1,2,3,4,5,6,7,8,9,10])  # actual content is the a result of another calculation step

y有一个公式，每个元素都是基于前一个元素，让i表示y的索引，每个元素是：

y[i] = y[i-1] * 2 + x[i]

计算第一个元素时，设y[i-1] = 50。换句话说，y 应该是：

[101, 204, 411, 826, 1657, 3320, 6647, 13302, 26613, 53236]

如何使用 numpy 计算 y？

【问题讨论】：

y 的起始值是多少，即y[0]？
嗯，有一个人为的起始值 50。
你写的代码有什么问题？
@moose，这只是伪代码，不是 numpy。
这是一个非齐次线性递推方程。您是否考虑过先找到解析解，然后看看是否可以有效编码？

标签： python numpy

【解决方案1】：

让我们在你的序列中构建一些项目：

y[0] = 2*y[-1] + x[0]
y[1] = 2*y[0] + x[1] = 4*y[-1] + 2*x[0] + x[1]
y[2] = 2*y[1] + x[2] = 8*y[-1] + 4*x[0] + 2*x[1] + x[2]
...
y[n] = 2**(n+1)*y[-1] + 2**n*x[0] + 2**(n-1)*x[1] + ... + x[n]

这可能不是很明显，但您可以使用 numpy 构建上述序列，执行以下操作：

n = len(x)
y_1 = 50
pot = 2**np.arange(n-1, -1, -1)
y = np.cumsum(pot * x) / pot + y_1 * 2**np.arange(1, n+1)
>>> y
array([  101,   204,   411,   826,  1657,  3320,  6647, 13302, 26613, 53236])

这种解决方案的缺点是它们不是很通用：问题的微小变化可能会使整个方法无用。但只要你能用一点代数解决一个问题，它几乎肯定会远远超过任何算法方法。

【讨论】：

numpy.cumsum() 比使用 for 循环快得多，对于我的测试应用程序来说快了大约 8 倍。

【解决方案2】：

如果您需要递归计算，如果您的 y[i] 应该依赖于同一运行中计算出的 y[i-1]，那么 numpy 中似乎没有内置解决方案，您需要使用简单的for循环：

y = np.empty(x.size)
last = 50
for i in range(x.size):
    y[i] = last = last * 2 + x[i]

看到这个问题：Is a "for" loop necessary if elements of the a numpy vector are dependant upon the previous element?

否则，您可以使用 numpy 在一行中实现您的公式：

y = np.concatenate(([50], y[:-1])) * 2 + x

解释：

y[:-1]

创建一个N-1 大小的数组：y_0, y_1, ... y_N-1。

np.concatenate(([50], y[:-1]))

创建一个N 大小的数组，第一个元素的起始值为 50。所以这个表达式基本上就是你的y[i-1]。

然后您可以使用 numpy 数组算法逐元素进行数学运算。

【讨论】：

不错！然后更改y的第一个元素。
输出不应该是[101, 204, 411, 826, 1657, 3320, 6647, 13302, 26613, 53236]吗？我猜与第二个元素不匹配？应该是 204 而不是 102
这似乎无法正确处理重复发生。看起来它的设计目的是从旧的y 和x 构造一个全新的y 数组，而所需的行为似乎是让给定的关系保持在单个y 数组中；期望的结果将是这个解决方案的一个固定点。
好的，我想我当时误解了要求。如果 OP 想要一个递归公式，这个答案对他不起作用。（我写这个答案时没有示例结果）
确实我想要的是递归，我没有意识到这一点，我希望 numpy 能够通过一些魔法切片来实现这一点。谢谢你的解释。我将继续使用 for 循环。

【解决方案3】：

也许最快最简洁的方法是使用scipy.signal.lfilter，它完全实现了你描述的那种递归关系：

from scipy.signal import lfilter
import numpy as np

x = np.array([1,2,3,4,5,6,7,8,9,10])

b = [1., 0.]
a = [1., -2.]
zi = np.array([2*50])  # initial condition
y, _ = lfilter(b, a, x, zi=zi)

结果将为np.float64，但您可以转换为例如np.int32 如果你需要的话：

>>> y.astype(np.int32)
array([  101,   204,   411,   826,  1657,  3320,  6647, 13302, 26613, 53236])

【讨论】：

【解决方案4】：

以下是使用 Python 的方法：

x = [1, 2, 3, 4, 5, 6, 7, 8 ,9, 10]
y = [50]
for i in range(len(x)):
    y.append(y[-1] * 2 + x[i])
y = y[1:]

您可能希望从最后一个元素计算它，以防止将新值用于下一个 i。

【讨论】：

我会假设“我如何用 numpy 计算 *”意味着一个有效的矢量化版本。
@AmiTavory 显然不是。我问过他在 cmets 中的代码有什么问题，答案是“它是伪代码”。但我也认为 Ferdinand Beyers 的回答应该是被接受（并赞成）的答案。
像y = [50*2 + x[0]]; for v in x[1:]: y.append(y[-1]*2 + v) 这样更简单，但您的解决方案还不正确。
嗯...我的意思是这是伪代码，要转换成numpy。我以为我很清楚我想要一个 numpy 解决方案。
那么您可以简单地将我在这里的代码块复制到您的问题中并说您想要相同的结果，但使用矢量化/高效的 numpy 解决方案完成。

【解决方案5】：

这是使用 numpy 的方法：

import numpy as np
x = np.array([ 1, 2, 3, 4, 5, 6, 7, 8 ,9, 10 ])
y = np.array([ 50 ])
for i in np.arange(len(x)):
    y = np.append(
                  y,
                  ( y[-1] * 2 + x[i] )
                  )
y = y[1:]

print(y)

【讨论】：