在相同长度的一维 numpy 数组上评估一维函数数组的有效算法

小能豆

在相同长度的一维 numpy 数组上评估一维函数数组的有效算法

我有一个长度为 N 的 (大) 数组，其中包含 k 个不同的函数，以及一个长度为 N 的横坐标数组。我想在横坐标处计算函数以返回长度为 N 的纵坐标数组，而且至关重要的是，我需要非常快地完成此操作。

我尝试通过调用 np.where 进行以下循环，但速度太慢：

创建一些虚假数据来说明问题：

def trivial_functional(i): return lambda x : i*x
k = 250
func_table = [trivial_functional(j) for j in range(k)]
func_table = np.array(func_table) # possibly unnecessary

我们有一张包含 250 个不同函数的表格。现在我创建一个大数组，其中包含这些函数的许多重复条目，以及一组长度相同的点，这些函数应在这些点处进行求值。

Npts = 1e6
abcissa_array = np.random.random(Npts)
function_indices = np.random.random_integers(0,len(func_table)-1,Npts)
func_array = func_table[function_indices]

最后，循环遍历数据使用的每个函数并在相关点集上对其进行评估：

desired_output = np.zeros(Npts)
for func_index in set(function_indices):
    idx = np.where(function_indices==func_index)[0]
    desired_output[idx] = func_table[func_index](abcissa_array[idx])

这个循环在我的笔记本电脑上花费约 0.35 秒，这是我的代码中最大的瓶颈。

阅读 18

2024-11-17

共1个答案

小能豆

在我的计算机上，它似乎也略快一些——根据粗略测试，大约 30 毫秒。

def apply_indexed_fast(array, func_indices, func_table):
    func_argsort = func_indices.argsort()
    func_ranges = list(np.searchsorted(func_indices[func_argsort], range(len(func_table))))
    func_ranges.append(None)
    out = np.zeros_like(array)
    for f, start, end in zip(func_table, func_ranges, func_ranges[1:]):
        ix = func_argsort[start:end]
        out[ix] = f(array[ix])
    return out

像您的一样，这会将索引序列拆分argsort为块，每个块对应于中的一个函数func_table。然后，它使用每个块为其对应的函数选择输入和输出索引。为了确定块边界，它使用np.searchsorted而不是np.unique– 其中searchsorted(a, b)可以被认为是二分搜索算法，它返回中第一个a等于或大于中给定值的值的索引b。

然后 zip 函数简单地并行迭代其参数，从每个参数中返回一个项目，收集到一个元组中，然后将它们串在一起形成一个列表。（因此zip([1, 2, 3], ['a', 'b', 'c'], ['b', 'c', 'd'])返回[(1, 'a', 'b'), (2, 'b', 'c'), (3, 'c', 'd')]。）这与语句内置的“解包”这些元组的功能一起for，允许以简洁但富有表现力的方式并行迭代多个序列。

在本例中，我使用它来迭代中的函数func_tables以及的两个不同步副本。这确保了变量中func_ranges的项始终比变量中的项领先一步。通过附加到，我确保最终块得到妥善处理——当其任何一个参数用尽项时停止，从而切断序列中的最终值。方便的是，该值还可以用作开放式切片索引！func_ranges``end``start``None``func_ranges``zip``None

另一个完成相同事情的技巧需要多几行，但内存开销较低，尤其是与,itertools等效一起使用时：zipizip

range_iter_a = iter(func_ranges)   # create generators that iterate over the 
range_iter_b = iter(func_ranges)   # values in `func_ranges` without making copies
next(range_iter_b, None)           # advance the second generator by one
for f, start, end in itertools.izip(func_table, range_iter_a, range_iter_b):
    ...

但是，这些基于生成器的低开销方法有时会比原始列表慢一点。另外，请注意，在 Python 3 中，zip行为更像izip。

2024-11-17