我在 Numpy 中无法解决这个问题。我需要模拟一个模拟最大跟踪器（电阻二极管电容）。我有一些很长的一维数组 X，我想从中计算输出数组 Y，这样
Y[0] = X[0]
Y[i] = max(0.99 * Y[i - 1], X[i])
我通过用Y^30 = ExpDecayFunc * X^30 近似我的上述规则来伪造它，其中星号是卷积。当然，我缺少一些更直接的东西吗？非常感谢！
【问题讨论】：
标签： python arrays numpy iteration filtering
您是否正在尝试模拟非对称信号滤波器（电阻器、二极管、电容器）？这是一个令人讨厌的非线性运算，无法并行计算。所以，这对于 NumPy 来说确实不是一件好事。
简单的解决方案是：
import numpy as np
# just do something random
X = np.random.random(1000000)
def my_filter(X):
Y = np.empty(len(X))
Y[0] = X[0]
for i in range(1, len(X)):
Y[i] = max(.99*Y[i-1], X[i])
return Y
这需要时间，我的机器需要 1.36 秒（物品需要 1.36 秒）。不大好。 （编辑：np.arange 的愚蠢用法改为range。）
可以通过重新排列算法来加快算法速度以避免查找：
def my_filter_2(X):
Y = np.empty(len(X))
Y[0] = X[0]
a = .99 * Y[0]
for i in range(1, len(X)):
a = max(a, X[i])
Y[i] = a
a *= .99
return Y
现在我们有 1.16 ms（每个元素 1.16 us）。一种改进，但毕竟不是很快。
但是我们有cython。这是通过 IPython 的 %%cython 完成的（不是我的解决方案，Andrew Jaffe 在他的出色答案中显示了这一点）：
%%cython
import numpy as np
cimport numpy as np
# just do something random
cdef np.ndarray cX = np.random.random(1000000)
def cy_filter(np.ndarray[np.double_t] X):
cdef int i
cdef np.ndarray[np.double_t] Y = np.empty(len(X))
Y[0] = X[0]
for i in range(1, len(X)):
Y[i] = max(.99*Y[i-1], X[i])
return Y
这很快！我的电脑声称 6.43 毫秒（6.43 ns/元素）。
另一个几乎是 Pythonic 的解决方案是 numba，正如 DSM 在他们的回答中所建议的那样：
from numba import autojit
import numpy as np
@autojit
def my_filter_nb(X, Y):
Y[0] = X[0]
for i in range(1, len(X)):
Y[i] = max(.99*Y[i-1], X[i])
return Y
def my_filter_fast(X):
Y = np.empty(len(X))
my_filter_nb(X, Y)
return Y
这给出了 4.18 毫秒（4.18 ns/元素）。
但如果我们仍然需要速度，让我们 C:
import numpy as np
import scipy.weave
X = np.random.random(1000000)
def my_filter_c(X):
x_len = len(X)
Y = np.empty(x_len)
c_source = """
#include <math.h>
int i;
double a, x;
Y(0) = X(0);
a = .99 * Y(0);
for (i = 1; i < x_len; i++)
{
x = X(i);
if (x > a)
a = x;
Y(i) = a;
a *= .99;
}
"""
scipy.weave.inline(c_source, ["X","Y","x_len"],
compiler="gcc",
headers=["<math.h>"],
type_converters=scipy.weave.converters.blitz)
return Y
这个给出了 3.72 毫秒（3.72 纳秒/轮）。 （顺便说一句，我的大脑不是多线程的，将内联 C 写入 Python 需要两个线程——用 C 编写一个简单的程序时会漏掉多少分号，这真是令人惊讶。）改进并没有那么大，麻烦的是。
看看这与普通 C 相比有多差或多好：
#include <stdio.h>
#include <stdlib.h>
#include <sys/resource.h>
#include <time.h>
#define NUMITER 100000000
int main(void)
{
double *x, *y;
double a, b, time_delta;
int i;
struct rusage ru0, ru1;
x = (double *)malloc(NUMITER * sizeof(double));
y = (double *)malloc(NUMITER * sizeof(double));
for (i = 0; i < NUMITER; i++)
x[i] = rand() / (double)(RAND_MAX - 1);
getrusage(RUSAGE_SELF, &ru0);
y[0] = x[0];
a = .99 * y[0];
for (i = 0; i < NUMITER; i++)
{
b = x[i];
if (b > a)
a = b;
y[i] = a;
a *= .99;
}
getrusage(RUSAGE_SELF, &ru1);
time_delta = ru1.ru_utime.tv_sec + ru1.ru_utime.tv_usec * 1e-6
- ru0.ru_utime.tv_sec - ru0.ru_utime.tv_usec * 1e-6;
printf("Took %.6lf seconds, %.2lf nanoseconds per element", time_delta, 1e9 * time_delta / NUMITER);
return (int)y[1234] % 2; // just to make sure the optimizer is not too clever
}
使用gcc -Ofast 编译需要 318 毫秒或 3.18 ns/元素（注意元素数量较多），因此是赢家。
所有 Python 计时都是使用 IPython 的 %timeit 执行的，它们包括来自 np.empty 的一些开销，但这非常微不足道。然而，可能由于内存管理问题，每次运行的结果都会有所不同，因此无论如何都需要小心谨慎。
我还尝试了包含 5 亿个元素的更快解决方案以避免调用开销：
%cython: 7.5 ns/元素numba: 7.3 ns/元素我也尝试了一些纯 C 的手动优化技巧，但至少在不查看编译结果的情况下，gcc 似乎至少和我一样聪明。
在这个堆栈中，我可能会选择numba 或纯 C，这取决于我的匆忙。与这个具体问题scipy.weave.inline相比优势太大了。
另外——取决于数据——并行处理可能会稍微快一点，但最坏的情况会更糟，而且整个事情可能会受到内存带宽的限制。
【讨论】：
您也可以使用numba，但需要进行一些更改：
from numba import autojit
import numpy as np
@autojit
def my_filter_nb(X, Y):
Y[0] = X[0]
for i in range(1, len(X)):
Y[i] = max(.99*Y[i-1], X[i])
return Y
def my_filter_fast(X):
Y = np.empty(len(X))
my_filter_nb(X, Y)
return Y
def my_filter(X):
Y = np.empty(len(X))
Y[0] = X[0]
for i in np.arange(1, len(X)):
Y[i] = max(.99*Y[i-1], X[i])
return Y
这给了我：
>>> X = np.random.random(1000000)
>>> %timeit my_filter(X)
1 loops, best of 3: 936 ms per loop
>>> %timeit my_filter_fast(X)
100 loops, best of 3: 3.83 ms per loop
>>> (my_filter(X) == my_filter_fast(X)).all()
True
【讨论】：
Cython 非常快。我在 iPython 中使用cython magic 运行了这个。
%% cython
import numpy as np
cimport numpy as np
# just do something random
cdef np.ndarray cX = np.random.random(1000000)
def cy_filter(np.ndarray[np.double_t] X):
cdef int i
cdef np.ndarray[np.double_t] Y = np.empty(len(X))
Y[0] = X[0]
for i in range(1, len(X)):
Y[i] = max(.99*Y[i-1], X[i])
return Y
使用%timeit，我从
1 loops, best of 3: 1.52 s per loop
到
100 loops, best of 3: 4.67 ms per loop
（对于它的价值，当我错过 cdef int i 时，它只是大约 3 倍的加速，而不是 300 倍！）
【讨论】：