计算我的 for 循环并希望使用某些函数在没有 for 循环的情况下进行计算

Question

dec = 0.1
data = np.array([100,200,300,400,500])

我有一个这样的 for 循环

y = np.zeros(len(data))
for i in range(len(data)):
    if i == 0:
        y[i] = (1.0 - dec) * data[i]
    else:
        y[i] = (1.0 - dec) * data[i] + (dec * y[i - 1])

输出 y 是：

array([ 90.   , 189.   , 288.9  , 388.89 , 488.889])

现在我想在没有循环的情况下进行上述计算，所以如果我破坏代码并这样做

data[0] = (1.0 - dec) * data[0]
data[1:] = (1.0 - dec) * data[1:] + (dec * data[0])

输出数据为：

array([ 90, 189, 279, 369, 459])

当您比较 y 和 data 输出时，前两个值是正确的，因为它与 data[0] 相乘，这是有道理的，但稍后它应该像循环代码中的循环一样继续，那么我们如何实现呢？有没有可以处理这个的函数？这主要是为了优化我的代码，使其对数千个数据运行得更快。

预期输出与 y 输出相同。

Answer 1

我们可以使用scipy.linalg.toeplitz 来制作数据移位矩阵，然后将其乘以dec 的幂并求和列：

import numpy as np
from scipy.linalg import toeplitz

dec = 0.1
data = np.array([100,200,300,400,500])

decs = np.power(dec, np.arange(len(data)))

r = np.zeros_like(data)
r[0] = data[0]
toep = toeplitz(r, data)

output = (1 - dec) * np.sum(toep * decs.reshape(-1, 1), axis=0)

第一个decs 是dec 的幂的向量：

print(decs) 
#[1.e+00 1.e-01 1.e-02 1.e-03 1.e-04]

接下来，我们使用toeplitz 制作data 的移位矩阵：

print(toep)
#[[100 200 300 400 500]
# [  0 100 200 300 400]
# [  0   0 100 200 300]
# [  0   0   0 100 200]
# [  0   0   0   0 100]])

最后，我们将decs 整形为一列，将其乘以toep 并沿列求和。这个结果需要用1 - dec缩放。

这一切都有效，因为我们可以将 data[i] 的定义重写为求和而不是递归：

y[i] = (1.0 - dec) * data[i] + (dec * y[i - 1])
y[i] = (1.0 - dec) * data[i] + (dec * ((1.0 - dec) * data[i - 1] + (dec * y[i - 2])))
...
y[i] = (1.0 - dec) * (data[i] + dec * data[i - 1] + dec ** 2 * data[i - 2] + ... dec ** i * data[0])
y[i] = (1.0 - dec) * sum(dec ** j * data[i - j] for j in range(i + 1))

这可以通过归纳来证明。

从那里开始，将这些总和重写为矩阵的列，并将该矩阵转换为 numpy/scipy 中的计算。

Answer 2

我知道你说过不要使用 python for 循环

但 np.vectorize 也不是真正的矢量化（它不会移动你的代码 C）它只是方便

既然你说的是 1000 条数据，所以你应该尝试 numba，因为它会将 for 循环移动到机器代码

到目前为止，我无法想到仅使用 numpy ufuncs(np.add、np.dot 等) 获得相同的 correct 输出，因为已知 ufuncs 可以矢量化（真正的 simd 矢量化），具体取决于是否编译器可以做到，所以也许现在你可以试试 numba

import numba as nb
import numpy as np

dec = 0.1
data = np.array([100,200,300,400,500])
@nb.jit((nb.int64[:],))
def f(data):
    y = np.zeros(data.shape[0])
    for i in range(y.shape[0]):
        if i == 0:
            y[i] = (1.0 - dec) * data[i]
        else:
            y[i] = (1.0 - dec) * data[i] + (dec * y[i - 1])
    return y
print(f(data))

也可以并行化（here），但我不确定在这种情况下如何正确执行。事实上，我怀疑在不添加更复杂代码的情况下是否可以解决这个问题的向量化和并行性