不同数组和索引的函数答案

【问题标题】：Function over different arrays and indexes不同数组和索引的函数
【发布时间】：2022-01-21 05:00:03
【问题描述】：

我想对两个数组 x = [1, 2, 3, ...] 和 y = [4, NaN, NaN, ...] 执行操作。目标是使用标量和前一个索引来计算数组 y 的剩余值，例如，y[1] = (y[0] * 7 + x[1]) / 8、y[2] = (y[1] * 7 + x[2]) / 8、......如何在 numpy 中创建这种操作？ for 循环效率太低了？

【问题讨论】：

由于y 的每个元素都依赖于前一个值，所以你必须一个一个地计算它们……向量化的操作在这里没有意义；矢量化需要将 something 加载到内存中，然后在两个数值数据数组上执行一条指令（如“add”）。换句话说，您需要所有数据先验才能利用矢量化。
从更数学的角度来看，您所描述的是通常所说的“递归关系”。它是一个纯粹以迭代术语定义的数学对象。换句话说，没有其他方法可以生成递归关系除了迭代。
@ddejohn 谢谢你的回复！我会保留我当前的 for 循环。
这是一个严肃的问题。一般来说，这不能在numpy 中优化，除非你使用一些铅笔和纸来获得A * y[i+1] + B * y[i] == x[i+1] 形式的循环关系的解析表达式。我已经在我的解决方案中展示了如何做到这一点，但它仍然是O(n²)，它比简单的循环慢得多。但是，如果您的数组x 像您的情况一样非常简单，并且您知道表达式(n-1) + (n-2)*7/8 + (n-3)*(7/8)**2+ ... 1*(7/8)**(n-1) 的一般公式，则可以改进它。

标签： python numpy

【解决方案1】：

你正试图这样做：

for i in range(N):
    y[i+1] = (y[i] * 7 + x[i+1]) / 8

所以我很确定将其称为“Python 级别的循环”是值得的。来自Introduction of numpy vectorized operations：

NumPy 提供了高度优化的函数，用于对数字数组执行数学运算。在 Python 中执行大量迭代（例如，通过“for-loops”）以执行重复的数学计算几乎总是应该由在数组上使用向量化函数来代替。这为 NumPy 的整个设计范式提供了信息。

在您的情况下，没有直接的矢量化方法。至少，您可以尝试数学方法，因为这种关系广为人知

A * y[i+1] + B * y[i] == f(i) 形式的关系称为线性一阶非齐次差分递推方程

为了更进一步，您可以将其重新排列为：

y[n+1] = (7/8) * y[n] + x[n+1]/8

根据这个公式，解决方案的剩余部分如下所示：

y[1] = (7/8) * y[0] + x[1]/8
y[2] = 7/8 * [(7/8) * y[0] + x[1]/8] + x[2]/8 = (7/8)**2*y[0] + 7/8*x[1]/8 + x[2]/8
...

所以你可以得出结论

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

因此不再需要查找上一个项目 y[i-1] 以查找 y[i]。现在让我们根据general numpy broadcasting rules的公式，一口气找到[y[1], y[2], ... y[n]]。

您需要添加Y =

y[0]*(7/8)
y[0]*(7/8)**2
y[0]*(7/8)**3
       ...
y[0]*(7/8)**(n

带有产品V =的结果

1   0        0        0  ...     0
1 7/8        0        0  ...     0
1 7/8 (7/8)**2        0  ...     0
1 7/8 (7/8)**2 (7/8)**3  ...     0
.  .      .        .             .
·  ·      ·        ·             ·
'  '      '        '             '
1 7/8 (7/8)**2 (7/8)**3  ... (7/8)**(n-1)

和X =

 1  0  0  0  0  0  0  0  0  0 ...
 2  1  0  0  0  0  0  0  0  0 ...
 3  2  1  0  0  0  0  0  0  0 ...
 4  3  2  1  0  0  0  0  0  0 ...
 5  4  3  2  1  0  0  0  0  0 ...
 6  5  4  3  2  1  0  0  0  0 ...
...

我会为此提供一个numpy 解决方案，不是最好看的，但你可以重构它：

n = 10
y = np.empty(n+1)
y[0] = 4
R = np.repeat([1, 7/8], [1, n])
Y = np.cumprod(R[1:])
V = np.tril(np.vander(R, increasing=True))[:-1,:-1]
X = np.lib.stride_tricks.sliding_window_view(np.r_[np.zeros(n-1), np.arange(1, n+1)], n)[:, ::-1]
T1, T2 = (Y*y[0], np.sum(V*X/8, axis=1))
y[1:] = T1 + T2
>>> y
[4.         3.625      3.421875   3.36914062 3.44799805 3.64199829
 3.9367485  4.31965494 4.77969807 5.30723581 5.89383134]

一般来说它是O(n²)，而初始方法是O(n)。如果您知道np.sum(V*X, axis=1) 的术语公式，则可以达到O(n) 并优于第一种方式。

【讨论】：