python - 迭代定义的 Numpy 数组创建

Question

我在 Numpy 中表达这个问题时遇到了问题。我需要模拟一个模拟最大值跟踪器(电阻二极管电容器)。我有一些很长的一维数组 X，我想从中计算输出数组 Y，这样

Y[0] = X[0]
Y[i] = max(0.99 * Y[i - 1], X[i])

我通过使用 Y^30 = ExpDecayFunc * X^30 近似我的上述规则来伪造它，其中星号是卷积。当然，我还缺少一些更直接的东西吗？非常感谢!

Answer 1

您是否正在尝试模拟非对称信号滤波器(电阻器、二极管、电容器)？这是一个讨厌的非线性操作，不能并行计算。所以，这对 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 毫秒(每个元素 1.16 微秒)。有进步，但毕竟不是很快。

但是我们有了 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 纳秒/元素)。

另一种几乎是 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 纳秒/元素)。

但如果我们仍然需要速度，让我们 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 纳秒/元素
• numba:7.3 ns/元素
• 内联 C(编织):5.7 ns/元素
• 纯 C:3.2 ns/元素

我还尝试了一些纯 C 的手动优化技巧，但至少在不查看编译结果的情况下，似乎 gcc 至少和我一样聪明。

在这个堆栈中，我可能会选择 numba 或纯 C，具体取决于我的匆忙。有了这个具体的问题，scipy.weave.inline 相比优势就太麻烦了。

另外——根据数据——这可能会通过并行处理变得稍微快一些，但最坏的情况会更糟，而且整个事情可能无论如何都受内存带宽限制。