用于评估相同长度的1d numpy阵列上的1-d函数数组的并行算法

Question

下面的结果是我有一个令人尴尬的并行for循环,我试图解决.解释这个问题还有一些麻烦,但尽管所有问题都很冗长,但我认为这应该是一个相当简单的问题,多处理模块的设计很容易解决.

我有一个由k个不同函数组成的大长度N数组,以及一个长度为N的abcissa数组.由于@senderle提供的聪明的解决方案在高效算法中描述,用于评估相同长度的1d numpy数组上的1-d函数数组,我有一个快速的基于numpy的算法,我可以用它来评估abcissa的函数返回长度为N的纵坐标数组:

def apply_indexed_fast(abcissa, func_indices, func_table):
    """ Returns the output of an array of functions evaluated at a set of input points 
    if the indices of the table storing the required functions are known. 

    Parameters 
    ----------
    func_table : array_like 
        Length k array of function objects

    abcissa : array_like 
        Length Npts array of points at which to evaluate the functions. 

    func_indices : array_like 
        Length Npts array providing the indices to use to choose which function 
        operates on each abcissa element. Thus func_indices is an array of integers 
        ranging between 0 and k-1. 

    Returns 
    -------
    out : array_like 
        Length Npts array giving the evaluation of the appropriate function on each 
        abcissa element. 
    """
    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(abcissa)

    for i in range(len(func_table)):
        f = func_table[i]
        start = func_ranges[i]
        end = func_ranges[i+1]
        ix = func_argsort[start:end]
        out[ix] = f(abcissa[ix])

    return out

我现在要做的是使用多处理来并行化此函数内的for循环.在描述我的方法之前,为了清楚起见,我将简要描述算法@senderle开发的工作原理.如果您可以阅读上述代码并立即理解,只需跳过下一段文字即可.

首先,我们找到对输入func_indices进行排序的索引数组,我们用它来定义length-k func_ranges整数数组.func_ranges的整数条目控制应用于输入abcissa的相应子数组的函数,其工作方式如下.设f为输入func_table中的第i个函数.那么我们应该应用函数f的输入abcissa的切片是slice(func_ranges [i],func_ranges [i + 1]).因此,一旦计算了func_ranges,我们就可以在输入func_table上运行一个简单的for循环并连续将每个函数对象应用于适当的切片,填充我们的输出数组.请参阅下面的代码,了解该算法的最小示例.

def trivial_functional(i): 
    def f(x):
        return i*x
    return f

k = 250
func_table = np.array([trivial_functional(j) for j in range(k)])

Npts = 1e6
abcissa = np.random.random(Npts)
func_indices = np.random.random_integers(0,len(func_table)-1,Npts)

result = apply_indexed_fast(abcissa, func_indices, func_table)

所以我现在的目标是使用多处理来并行化这个计算.我认为这很简单,使用我常用的技巧来线程尴尬地并行循环.但我在下面的尝试提出了一个我不明白的例外.

from multiprocessing import Pool, cpu_count
def apply_indexed_parallelized(abcissa, 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(abcissa)

    num_cores = cpu_count()
    pool = Pool(num_cores)

    def apply_funci(i):
        f = func_table[i]
        start = func_ranges[i]
        end = func_ranges[i+1]
        ix = func_argsort[start:end]
        out[ix] = f(abcissa[ix])

    pool.map(apply_funci, range(len(func_table)))
    pool.close()

    return out

result = apply_indexed_parallelized(abcissa, func_indices, func_table)
PicklingError: Can't pickle <type 'function'>: attribute lookup __builtin__.function failed

我在其他地方已经看过SO:多处理:如何在类中定义的函数上使用Pool.map？.我一个接一个地尝试了那里提出的每种方法; 在所有情况下,我得到一个"打开太多文件"错误,因为线程从未关闭,或者改编的算法只是挂起.这似乎应该有一个简单的解决方案,因为这只不过是一个令人尴尬的并行循环线程.

Answer 1

警告/警告:

您可能不想申请multiprocessing此问题.你会发现在大型数组上运行相对简单,问题就是内存绑定了numpy.瓶颈是将数据从RAM移动到CPU缓存.CPU缺乏数据,因此在问题上投入更多CPU并没有多大帮助.此外,您当前的方法将为输入序列中的每个项目挑选并制作整个数组的副本,这会增加大量开销.

有很多地方的情况下numpy+ multiprocessing是非常有效的,但你必须确保你有一个CPU绑定的问题的工作.理想情况下,它是一个CPU限制的问题,具有相对较小的输入和输出,以减轻酸洗输入和输出的开销.对于numpy最常用的许多问题,情况并非如此.

你当前的方法有两个问题

回答你的问题:

您的即时错误是由于传入了一个无法从全局范围访问的函数(即函数内定义的函数).

但是,你有另一个问题.您正在将numpy数组视为可以由每个进程修改的共享内存.相反,当使用multiprocessing原始数组时,将被腌制(有效地制作副本)并独立地传递给每个进程.永远不会修改原始数组.

避免 PicklingError

作为重现错误的最小示例,请考虑以下事项:

import multiprocessing

def apply_parallel(input_sequence):
    def func(x):
        pass
    pool = multiprocessing.Pool()
    pool.map(func, input_sequence)
    pool.close()

foo = range(100)
apply_parallel(foo)

这将导致:

PicklingError: Can't pickle <type 'function'>: attribute lookup 
               __builtin__.function failed

当然,在这个简单的例子中,我们可以简单地将函数定义移回__main__命名空间.但是,在你的,你需要它来引用传入的数据.让我们看一个更接近你正在做的事情的例子:

import numpy as np
import multiprocessing

def parallel_rolling_mean(data, window):
    data = np.pad(data, window, mode='edge')
    ind = np.arange(len(data)) + window

    def func(i):
        return data[i-window:i+window+1].mean()

    pool = multiprocessing.Pool()
    result = pool.map(func, ind)
    pool.close()
    return result

foo = np.random.rand(20).cumsum()
result = parallel_rolling_mean(foo, 10)

有多种方法可以解决这个问题,但常见的方法是:

import numpy as np
import multiprocessing

class RollingMean(object):
    def __init__(self, data, window):
        self.data = np.pad(data, window, mode='edge')
        self.window = window

    def __call__(self, i):
        start = i - self.window
        stop = i + self.window + 1
        return self.data[start:stop].mean()

def parallel_rolling_mean(data, window):
    func = RollingMean(data, window)
    ind = np.arange(len(data)) + window

    pool = multiprocessing.Pool()
    result = pool.map(func, ind)
    pool.close()
    return result

foo = np.random.rand(20).cumsum()
result = parallel_rolling_mean(foo, 10)

大!有用!

但是,如果这个规模到大型阵列,你很快就会发现,它要么运行速度很慢(您可以通过增加加快chunksize在pool.map通话),或者你会很快出来的RAM运行(一旦你增加chunksize) .

multiprocessingpickle输入,以便它可以传递给独立和独立的python进程.这意味着您要为每个操作的人制作整个阵列的副本. i

我们稍后会回到这一点......

multiprocessing 不在进程之间共享内存

该multiprocessing模块通过pickling输入并将它们传递给独立进程来工作.这意味着如果您在一个进程中修改某些内容,则其他进程将看不到该修改.

但是,multiprocessing还提供了存在于共享内存中的原语,并且可以由子进程访问和修改.有几种不同的方法可以调整numpy数组以使用共享内存multiprocessing.Array.但是,我建议首先避免使用这些(如果您不熟悉它,请阅读虚假共享).有些情况下它非常有用,但通常是为了节省内存,而不是为了提高性能.

因此,最好在单个进程中对大型数组进行所有修改(对于一般IO,这也是一个非常有用的模式).它不一定是"主要"过程,但最容易思考这种方式.

作为一个例子,假设我们希望让我们的parallel_rolling_mean函数采用输出数组来存储内容.一个有用的模式类似于以下内容.注意使用迭代器并仅在主进程中修改输出:

import numpy as np
import multiprocessing

def parallel_rolling_mean(data, window, output):
    def windows(data, window):
        padded = np.pad(data, window, mode='edge')
        for i in xrange(len(data)):
            yield padded[i:i + 2*window + 1]

    pool = multiprocessing.Pool()
    results = pool.imap(np.mean, windows(data, window))
    for i, result in enumerate(results):
        output[i] = result
    pool.close()

foo = np.random.rand(20000000).cumsum()
output = np.zeros_like(foo)
parallel_rolling_mean(foo, 10, output)
print output

希望这个例子有助于澄清一些事情.

chunksize 和表现

关于性能的一个快速说明:如果我们扩展它,它将非常快地变得非常慢.如果您查看系统监视器(例如top/ htop),您可能会注意到您的核心在大多数时间处于空闲状态.

默认情况下,主进程会为每个进程选择每个输入并立即传入,然后等待它们完成以挑选下一个输入.在许多情况下,这意味着主进程工作,然后在工作进程繁忙时处于空闲状态,然后工作进程处于空闲状态,而主进程正在挑选下一个输入.

关键是增加chunksize参数.这将导致pool.imap为每个过程"预先腌制"指定数量的输入.基本上,主线程可以保持繁忙的酸洗输入,并且工作进程可以保持繁忙的处理.缺点是你正在使用更多的内存.如果每个输入占用大量RAM,这可能是一个坏主意.但是,如果没有,这可以大大加快速度.

作为一个简单的例子:

import numpy as np
import multiprocessing

def parallel_rolling_mean(data, window, output):
    def windows(data, window):
        padded = np.pad(data, window, mode='edge')
        for i in xrange(len(data)):
            yield padded[i:i + 2*window + 1]

    pool = multiprocessing.Pool()
    results = pool.imap(np.mean, windows(data, window), chunksize=1000)
    for i, result in enumerate(results):
        output[i] = result
    pool.close()

foo = np.random.rand(2000000).cumsum()
output = np.zeros_like(foo)
parallel_rolling_mean(foo, 10, output)
print output

使用chunksize=1000,处理200万个元素数组需要21秒:

python ~/parallel_rolling_mean.py  83.53s user 1.12s system 401% cpu 21.087 total

但是chunksize=1(默认值)它需要大约8倍(2分41秒).

python ~/parallel_rolling_mean.py  358.26s user 53.40s system 246% cpu 2:47.09 total

事实上,使用默认的chunksize,它实际上远比同一件事的单进程实现差,只需要45秒:

python ~/sequential_rolling_mean.py  45.11s user 0.06s system 99% cpu 45.187 total