python - 如何在蒙特卡罗积分中实现多处理

Question

我创建了一个 Python 程序，该程序使用蒙特卡罗模拟在给定时间间隔内对给定函数进行积分。它工作得很好，除了当你想要更高的准确度(更大的 N 值)时它运行得非常慢。我想我应该尝试多处理以加快速度，但后来我意识到我不知道如何实现它。这是我现在所拥有的:

from scipy import random
import numpy as np
import matplotlib.pyplot as plt
from multiprocessing import Process
import os

# GOAL: Approximate the integral of a function f(x) from lower bound a to upper bound b using Monte Carlo simulation

# bounds of integration
a = 0
b = np.pi

# function to integrate
def f(x):
    return np.sin(x)

N = 10000
areas = []


def mcIntegrate():
    
    for i in range(N):
        
        # array filled with random numbers between limits
        xrand = random.uniform(a, b, N)
        
        # sum the return values of the function of each random number
        integral = 0.0
        for i in range(N):
            integral += f(xrand[i])
        
        # scale integral by difference of bounds divided by amount of random values
        ans = integral * ((b - a) / float(N))
        
        # add approximation to list of other approximations
        areas.append(ans)


if __name__ == "__main__":

    processes = []
    numProcesses = os.cpu_count()

    for i in range(numProcesses):
        process = Process(target=mcIntegrate)
        processes.append(process)

    for process in processes:
        process.start()

    for process in processes:
        process.start()

    # graph approximation distribution
    plt.title("Distribution of Approximated Integrals")
    plt.hist(areas, bins=30, ec='black')
    plt.xlabel("Areas")
    plt.show()

我可以获得有关此实现的帮助吗？

---

听取了评论的建议并使用了 multiprocessor.Pool，并且还通过使用 NumPy 减少了一些操作。从大约 5 分钟运行到现在大约 6 秒(对于 N = 10000)。这是我的实现:

import scipy
import numpy as np
import matplotlib.pyplot as plt
import multiprocessing
import os

# GOAL: Approximate the integral of function f from lower bound a to upper bound b using Monte Carlo simulation

a = 0       # lower bound of integration
b = np.pi   # upper bound of integration
f = np.sin  # function to integrate
N = 10000   # sample size

def mcIntegrate(p):
    xrand = scipy.random.uniform(a, b, N)       # create array filled with random numbers within bounds
    integral = np.sum(f(xrand))                 # sum return values of function of each random number
    approx = integral * ((b - a) / float(N))    # scale integral by difference of bounds divided by sample size
    return approx


if __name__ == "__main__":
    
    # run simulation N times in parallel and store results in array
    with multiprocessing.Pool(os.cpu_count()) as pool:
        areas = pool.map(mcIntegrate, range(N))

    # graph approximation distribution
    plt.title("Distribution of Approximated Integrals")
    plt.hist(areas, bins=30, ec='black')
    plt.xlabel("Areas")
    plt.show()

Answer 1

当我开始优化它时，事实证明这是一个比我想象的更有趣的问题。基本方法非常简单:

from multiprocessing import pool

def f(x):
    return x

results = pool.map(f, range(100))

这是适合多处理的 mcIntegerate:

from tqdm import tqdm

def mcIntegrate(steps):
    tasks = []

    print("Setting up simulations")

    # linear
    for _ in tqdm(range(steps)):
        xrand = random.uniform(a, b, steps)
        for i in range(steps):
            tasks.append(xrand[i])

    pool = Pool(cpu_count())

    print("Simulating (no progress)")
    results = pool.map(f, tasks)
    pool.close()

    print("summing")
    areas = []
    for chunk in tqdm(range(steps)):
        vals = results[chunk * steps : (chunk + 1) * steps]
        integral = sum(vals)
        ans = integral * ((b - a) / float(steps))
        areas.append(ans)

    return areas

tqdm 只是用来显示进度条。

这是多处理的基本工作流程:将问题分解为任务，解决所有任务，然后将它们重新添加到一起。确实，给定的代码有效。 (请注意，我已将您的 N 更改为 steps)。

为了完整起见，脚本现在开始:

from scipy import random
import numpy as np
import matplotlib.pyplot as plt
from multiprocessing import Pool, cpu_count
from tqdm import tqdm

# function to integrate
def f(x):
    return np.sin(x)

结束

areas = mcIntegrate(3_000)

a = 0
b = np.pi

plt.title("Distribution of Approximated Integrals")
plt.hist(areas, bins=30, ec="black")
plt.xlabel("Areas")
plt.show()

优化

我故意将问题分解到尽可能小的层面。这是个好主意吗？为了回答这个问题，请考虑:我们如何优化生成任务的线性过程？目前这确实需要相当长的时间。我们可以并行化它:

def _prepare(steps):
    xrand = random.uniform(a, b, steps)
    return [xrand[i] for i in range(steps)]

def mcIntegrate(steps):
    ...
    tasks = []
    for res in tqdm(pool.imap(_prepare, (steps for _ in range(steps))), total=steps):
        tasks += res  # slower except for very large steps

这里我使用了pool.imap，它返回一个迭代器，我们可以在结果可用后立即对其进行迭代，从而允许我们构建进度条。如果您这样做并进行比较，您会发现它的运行速度比线性解决方案慢。删除进度条(在我的机器上)并替换为:

    import time

    start = time.perf_counter()
    results = pool.map(_prepare, (steps for _ in range(steps)))
    tasks = []
    for res in results:
        tasks += res
    print(time.perf_counter() - start)

只是稍微快一点:它仍然比线性运行慢。将数据序列化到进程然后反序列化会产生开销。如果你试图获得整个事情的进度条，它会变得极其缓慢:

    results = []
    for result in tqdm(pool.imap(f, tasks), total=len(tasks)):
        results.append(result)

那么在更高层次上迭代又如何呢？这是 mcIterate 的另一个改编:

a = 0
b = np.pi

def _mcIntegrate(steps):
    xrand = random.uniform(a, b, steps)
    integral = 0.0
    for i in range(steps):
        integral += f(xrand[i])
    ans = integral * ((b - a) / float(steps))

    return ans


def mcIntegrate(steps):
    areas = []
    p = Pool(cpu_count())
    for ans in tqdm(p.imap(_mcIntegrate, ((steps) for _ in range(steps))), total=steps):
        areas.append(ans)

    return areas

这在我的机器上要快得多。它也简单得多。我期待有差异，但没想到差异这么大。

要点

多重处理不是免费的。像 np.sin() 这样简单的东西对于多进程来说太便宜了:我们花钱序列化、反序列化、追加等等，所有这些都是为了一个 sin( )计算。但如果您进行太多计算，就会因为失去粒度而浪费时间。这里的效果比我预期的更显着。了解特定问题的正确粒度级别的唯一方法...就是分析并尝试。