如何在py3中向量化这个函数答案

【问题标题】：How to vectorize this function in py3如何在py3中向量化这个函数
【发布时间】：2016-08-18 07:30:27
【问题描述】：

我似乎无法找出如何向量化这个 py3 循环

import numpy as np

a = np.array([-72, -10, -70, 37, 68, 9, 1, -3, 2, 3, -6, -4, ], np.int16)
result = np.array([-72, -10, -111, -23, 1, -2, 1, -3, 1, 2, -5, -5, ], np.int16)

b = np.copy(a)
for i in range(2, len(b)):
    b[i] += int( (b[i-1] + b[i-2]) / 2)

assert (b == result).all()

我尝试使用np.convolve 和pandas.rolling_apply，但无法正常工作。也许现在是学习 c-extensions 的时候了？

对于大约 500k 个元素的输入数组，将时间缩短到 50..100 毫秒会很棒。

@hpaulj 在他的回答中要求b[k] 的封闭表达式为a[:k]。我不认为它存在，但我对其进行了一些研究，并且确实发现封闭形式包含一堆 Jacobsthal 数字，正如@Divakar 指出的那样。

这是一个封闭的形式：

J_n 这里是Jacobsthal number，当这样展开时：

J_n = (2^n - (-1)^n) / 3

最终得到一个表达式，我可以想象它使用矢量化实现......

【问题讨论】：

但是b 和result 是不同的，对吧？哪个是正确的预期输出？
@Divakar 整数除法问题可能吗？它们在 Python 3.5 中是相同的。
@ayhan 没错，这应该在 py3.5 上运行。 @Divakar 正确的预期输出是 result
@ayhan 是的，因为当我们使用2.0 进行除法时它可以工作。
好吧，似乎封闭形式将遵循 Jacobsthal 数。所以，也许调查一下。

标签： python-3.x numpy

【解决方案1】：

大多数numpy 代码一次对整个数组进行操作。好的，它在 C 代码中迭代，但以与先使用哪个元素无关的方式进行缓冲。

此处对b[2] 的更改会影响为b[3] 计算的值，并且会持续下去。

add.at 和其他类似的ufunc 进行无缓冲计算。这允许您向一个元素重复添加一些值。在那种情况下我玩了一下，但到目前为止没有运气。

cumsum 和 cumprod 对于值依赖于早期问题的问题也很方便。

是否可以概括计算，以便根据所有a[:i]定义b[i]。我们知道b[2] 是a[:2] 的函数，但是b[3] 呢？

即使我们为浮点数工作，它也可能在进行整数除法时关闭。

【讨论】：

【解决方案2】：

我认为您已经有了理智的解决方案。任何其他矢量化都将依赖浮点计算，并且很难跟踪误差累积。例如，假设你想要一个矩阵向量乘法：对于前七个项，矩阵看起来像

array([[ 1.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ],
       [ 0.     ,  1.     ,  0.     ,  0.     ,  0.     ,  0.     ,  0.     ],
       [ 0.5    ,  0.5    ,  1.     ,  0.     ,  0.     ,  0.     ,  0.     ],
       [ 0.25   ,  0.75   ,  0.5    ,  1.     ,  0.     ,  0.     ,  0.     ],
       [ 0.375  ,  0.625  ,  0.75   ,  0.5    ,  1.     ,  0.     ,  0.     ],
       [ 0.3125 ,  0.6875 ,  0.625  ,  0.75   ,  0.5    ,  1.     ,  0.     ],
       [ 0.34375,  0.65625,  0.6875 ,  0.625  ,  0.75   ,  0.5    ,  1.     ]])

关系可以用迭代公式来描述

                       [ a[i-2] ]
b[i] = [0.5 , 0.5 , 1] [ a[i-1] ]
                       [ a[i]   ]

这定义了一系列具有单位矩阵形式的基本矩阵

[0 ... 0.5 0.5 1 0 ... 0]

在第 i 行。连续乘法给出了前七项的上述矩阵。确实有一个次对角结构，但是这些项很快就变得太小了。正如您所展示的，2 次方 500k 并不好玩。

为了跟踪浮点噪声，需要一个迭代解决方案，这就是您所拥有的。

【讨论】：