Processes within a multiprocessing.Pool exhibit very different performance when using map()

Viewed 48

I'm looking to solve a differential equation N times, each time with a different set of parameters. multiprocessing therefore sounds like the right tool for the job:let's define a solve function that does the actual work of computing the solution and use a multiprocessing.Pool with map to distribute the work among several processes.

Here the code with N=8 and 2 processes:

from multiprocessing import Pool
import time
import numpy as np

from bench import init_model, get_initial_solution

model = init_model()
sol_init = get_initial_solution(model, np.linspace(0, 1, 2), {"Current": 0.67})

Nsteps = 10
step_solver = model.default_solver

def solve(ind):
    st = time.time()
    step_solution = sol_init
    for step in range(0, Nsteps):
        step_solution = step_solver.step(
            step_solution, 
            model, 
            dt=1, 
            npts=2, 
            inputs={"Current": 2.0}, 
            save=False
    )
    return f"Task {ind} took {time.time() - st:.2f}s"


if __name__ == "__main__":
    with Pool(processes=2) as p:
        times = p.map(solve, np.arange(1, 9))
        print("\n".join(times))

For debugging purposes, solve does not return the solution, but instead the time spent by the process inside the function.

Executing the above (my computer has 4 cores), I get:

Task 1 took 4.41s
Task 2 took 5.59s
Task 3 took 1.67s
Task 4 took 0.62s
Task 5 took 0.61s
Task 6 took 0.72s
Task 7 took 0.68s
Task 8 took 0.53s

As you can see the time spent in the function solve varies widely across the process pool, and across a large range of values. Note that these result's aren't deterministic. i.e if I executed the above script again, a very different set of times would be observed. There is no reason for this randomness though, the amount of work to be done is the same across processes and executions.

Let's profile the execution of the processes to get a bit more information about what's going on in there.

import cProfile
import pstats

def profile(ind):
    cProfile.runctx("solve(ind)", globals(), locals(), "report_"+str(ind)+".txt")

with Pool(processes=2) as p:
    times = p.map(profile, np.arange(1,9))

for ind in range(1,9):
    stats = pstats.Stats("report_"+str(ind)+".txt").strip_dirs()
    stats.sort_stats("cumulative")
    stats.print_stats(11)

For instance if we look at the reports for tasks 1 and 7:

Thu Oct 29 18:07:18 2020    report_1.txt

         75858 function calls (75426 primitive calls) in 0.895 seconds

   Ordered by: cumulative time
   List reduced from 249 to 11 due to restriction <11>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.895    0.895 {built-in method builtins.exec}
        1    0.000    0.000    0.895    0.895 <string>:1(<module>)
        1    0.000    0.000    0.895    0.895 python-1mtcbY:15(solve)
       10    0.001    0.000    0.895    0.090 base_solver.py:712(step)
       10    0.001    0.000    0.892    0.089 scipy_solver.py:35(_integrate)
       10    0.003    0.000    0.889    0.089 ivp.py:156(solve_ivp)
       88    0.001    0.000    0.738    0.008 base.py:159(step)
       88    0.010    0.000    0.738    0.008 bdf.py:296(_step_impl)
       45    0.000    0.000    0.668    0.015 bdf.py:216(lu)
       45    0.666    0.015    0.668    0.015 decomp_lu.py:15(lu_factor)
     1284    0.003    0.000    0.139    0.000 base_solver.py:906(__call__)


Thu Oct 29 18:07:27 2020    report_7.txt

         75831 function calls (75399 primitive calls) in 6.773 seconds

   Ordered by: cumulative time
   List reduced from 244 to 11 due to restriction <11>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    6.773    6.773 {built-in method builtins.exec}
        1    0.000    0.000    6.773    6.773 <string>:1(<module>)
        1    0.000    0.000    6.773    6.773 python-1mtcbY:15(solve)
       10    0.001    0.000    6.773    0.677 base_solver.py:712(step)
       10    0.000    0.000    6.770    0.677 scipy_solver.py:35(_integrate)
       10    0.002    0.000    6.769    0.677 ivp.py:156(solve_ivp)
       88    0.001    0.000    6.612    0.075 base.py:159(step)
       88    0.011    0.000    6.612    0.075 bdf.py:296(_step_impl)
       45    0.000    0.000    6.520    0.145 bdf.py:216(lu)
       45    6.519    0.145    6.520    0.145 decomp_lu.py:15(lu_factor)
     1284    0.003    0.000    0.146    0.000 base_solver.py:906(__call__)

The above tells us that in the two cases (it is the case for all tasks), the process spends the most of its lifetime inside Scipy's lu_factor. The suprising bit however (at least surprising to me) is that the time spent in lu_factor varies - widely - across tasks. This is what I would like to understand.

2 Answers

presumably this is due to global mutable state; the first two tasks are "slow" which corresponds to the number of processes in the pool.

I'd suggest making things more deterministic by moving

model = init_model()
sol_init = get_initial_solution(model, np.linspace(0, 1, 2), {"Current": 0.67})

inside solve

lu_factor is a wrapper around LAPACK's getrf routines, which itself relies on the BLAS. In my case, the underlying BLAS implementation is OpenBlas, provided inside the scipy wheel when installing scipy from PyPI with pip.
As Sam Mason suggested, running lu_factor utilizes the four cores on my machine, because the underlying OpenBLAS is multithreaded. Disabling multithreading for OpenBLAS with export OPENBLAS_NUM_THREADS=1 largely solves the he performance issues described above. Indeed, running the same script with 2 processes again:

Task 1 took 0.19s
Task 2 took 0.19s
Task 3 took 0.18s
Task 4 took 0.18s
Task 5 took 0.18s
Task 6 took 0.18s
Task 7 took 0.18s
Task 8 took 0.18s

I note that the first two tasks are consistently slightly (around 10%) slower than the rest, which may be linked to Sam Mason's original answer about the impact of global mutable state, although I don't fully understand why at this point.

Related