Running scipy.integrate.ode in multiprocessing Pool results in huge performance hit
<p>I'm using python's <code>scipy.integrate</code> to simulate a 29-dimensional linear system of differential equations. Since I need to solve several problem instances, I thought I could speed it up by doing computations in parallel using <code>multiprocessing.Pool</code>. Since there is no shared data or synchronization necessary between threads (the problem is embarrassingly parallel), I thought this should obviously work. After I wrote the code to do this, however, I got very strange performance measurements:</p>

<ul>
<li>Single-threaded, without jacobian: 20-30 ms per call</li>
<li>Single-threaded, with jacobian: 10-20 ms per call</li>
<li>Multi-threaded, without jacobian: 20-30 ms per call</li>
<li><strong>Multi-threaded, with jacobian: 10-5000 ms per call</strong></li>
</ul>

<p>What's shocking is that what I thought should be the fastest setup, was actually the slowest, and the variability was <em>two orders of magnitude</em>. It's a deterministic computation; computers aren't supposed to work this way. What could possibly be causing this?</p>

<h3>Effect seems system-dependent</h3>

<p>I tried the same code on another computer and I didn't see this effect. </p>

<p>Both machines were using Ubuntu 64 bit, Python 2.7.6, scipy version 0.18.0, and numpy version 1.8.2. I didn't see the variability with an Intel(R) Core(TM) i5-5300U CPU @ 2.30GHz processor. I did see the issue with an <a href="http://www.cpu-world.com/CPUs/Core_i7/Intel-Core%20i7%20Mobile%20i7-2670QM.html">Intel(R) Core(TM) i7-2670QM CPU @ 2.20GHz</a>.</p>

<h3>Theories</h3>

<p>One thought was that there might be a shared cache among processors, and by running it in parallel I can't fit two instances of the jacobian matrix in the cache, so they constantly battle each other for the cache slowing each other down compared with if they are run serially or without the jacobian. But it's not a million variable system. The jacobian is a 29x29 matrix, which takes up 6728 bytes. The level 1 cache on the processor is <a href="http://www.cpu-world.com/CPUs/Core_i7/Intel-Core%20i7%20Mobile%20i7-2670QM.html">4 x 32 KB</a>, much larger. Are there any other shared resources between processors that might be to blame? How can we test this?</p>

<p>Another thing I noticed is that each python process seems to take several hundred percent of the CPU as it's running. This seems to mean that the code is already parallelized at some point (perhaps in the low-level library). This could mean that further parallelization wouldn't help, but I wouldn't expect such a dramatic slowdown.</p>

<h3>Code</h3>

<p>It would be good to try out the on more machines to see if (1) other people can experience the slowdown at all and (2) what are the common features of systems where the slowdown occurs. The code does 10 trials of two parallel computations using a multiprocessing pool of size two, printing out the time per scipy.ode.integrate call for each of the 10 trials. </p>

<pre><code>'odeint with multiprocessing variable execution time demonsrtation'

from numpy import dot as npdot
from numpy import add as npadd
from numpy import matrix as npmatrix
from scipy.integrate import ode
from multiprocessing import Pool
import time

def main():
    "main function"

    pool = Pool(2) # try Pool(1)
    params = [0] * 2

    for trial in xrange(10):
        res = pool.map(run_one, params)
        print "{}. times: {}ms, {}ms".format(trial, int(1000 * res[0]), int(1000 * res[1]))

def run_one(_):
    "perform one simulation"

    final_time = 2.0
    init_state = [0.1 if d &lt; 7 else 0.0 for d in xrange(29)]
    (a_matrix, b_vector) = get_dynamics()

    derivative = lambda dummy_t, state: npadd(npdot(a_matrix, state), b_vector)
    jacobian = lambda dummy_t, dummy_state: a_matrix
    #jacobian = None # try without the jacobian

    #print "jacobian bytes:", jacobian(0, 0).nbytes

    solver = ode(derivative, jacobian)
    solver.set_integrator('vode')
    solver.set_initial_value(init_state, 0)

    start = time.time()
    solver.integrate(final_time)
    dif = time.time() - start

    return dif

def get_dynamics():
    "return a tuple (A, b), which are the system dynamics x' = Ax + b"

    return \
    (
        npmatrix([
        [0, 0, 0, 0.99857378006, 0.053384274244, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ],
        [0, 0, 1, -0.003182219341, 0.059524655342, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ],
        [0, 0, -11.570495605469, -2.544637680054, -0.063602626324, 0.106780529022, -0.09491866827, 0.007107574493, -5.20817921341, -23.125876742495, -4.246931301528, -0.710743697134, -1.486697327603, -0.044548215175, 0.03436637817, 0.022990248611, 0.580153205353, 1.047552018229, 11.265023544535, 2.622275290571, 0.382949404795, 0.453076470454, 0.022651889536, 0.012533628369, 0.108399390974, -0.160139432044, -6.115359574845, -0.038972389136, 0, ],
        [0, 0, 0.439356565475, -1.998182296753, 0, 0.016651883721, 0.018462046981, -0.001187470742, -10.778778281386, 0.343052863546, -0.034949331535, -3.466737362551, 0.013415853489, -0.006501746896, -0.007248032248, -0.004835912875, -0.152495086764, 2.03915052839, -0.169614300211, -0.279125393264, -0.003678218266, -0.001679708185, 0.050812027754, 0.043273505033, -0.062305315646, 0.979162836629, 0.040401368402, 0.010697028656, 0, ],
        [0, 0, -2.040895462036, -0.458999156952, -0.73502779007, 0.019255757332, -0.00459562242, 0.002120360732, -1.06432932386, -3.659159530947, -0.493546966858, -0.059561101143, -1.953512259413, -0.010939065041, -0.000271004496, 0.050563886711, 1.58833954495, 0.219923768171, 1.821923233098, 2.69319056633, 0.068619628466, 0.086310028398, 0.002415425662, 0.000727041422, 0.640963888079, -0.023016712545, -1.069845542887, -0.596675149197, 0, ],
        [-32.103607177734, 0, -0.503355026245, 2.297859191895, 0, -0.021215811372, -0.02116791904, 0.01581159234, 12.45916782984, -0.353636907076, 0.064136531117, 4.035326800046, -0.272152744884, 0.000999589868, 0.002529691904, 0.111632959213, 2.736421830861, -2.354540136198, 0.175216915979, 0.86308171287, 0.004401276193, 0.004373406589, -0.059795009475, -0.051005479746, 0.609531777761, -1.1157829788, -0.026305051933, -0.033738880627, 0, ],
        [0.102161169052, 32.057830810547, -2.347217559814, -0.503611564636, 0.83494758606, 0.02122657001, -0.037879735231, 0.00035400386, -0.761479736492, -5.12933410588, -1.131382179292, -0.148788337148, 1.380741054924, -0.012931029503, 0.007645723855, 0.073796656681, 1.361745395486, 0.150700793731, 2.452437244444, -1.44883919298, 0.076516270282, 0.087122640348, 0.004623192159, 0.002635233443, -0.079401941141, -0.031023369979, -1.225533436977, 0.657926151362, 0, ],
        [-1.910972595215, 1.713829040527, -0.004005432129, -0.057411193848, 0, 0.013989634812, -0.000906753354, -0.290513515472, -2.060635522957, -0.774845915178, -0.471751979387, -1.213891560083, 5.030515136324, 0.126407660877, 0.113188603433, -2.078420624662, -50.18523312358, 0.340665548784, 0.375863242926, -10.641168797333, -0.003634153255, -0.047962774317, 0.030509705209, 0.027584169642, -10.542357589006, -0.126840767097, -0.391839285172, 0.420788121692, 0, ],
        [0.126296110212, -0.002898250629, -0.319316070797, 0.785201711657, 0.001772374259, 0.00000584372, 0.000005233812, -0.000097899495, -0.072611454126, 0.001666291957, 0.195701043078, 0.517339177294, 0.05236528267, -0.000003359731, -0.000003009077, 0.000056285381, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ],
        [-0.018114066432, 0.077615035084, 0.710897211118, 2.454275059389, -0.012792968774, 0.000040510624, 0.000036282541, -0.000678672106, 0.010414324729, -0.044623231468, 0.564308412696, -1.507321670112, 0.066879720068, -0.000023290783, -0.00002085993, 0.000390189123, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ],
        [-0.019957254425, 0.007108972111, 122.639137999354, 1.791704310155, 0.138329792976, 0.000000726169, 0.000000650379, -0.000012165459, -8.481152717711, -37.713895394132, -93.658221074435, -4.801972165378, -2.567389718833, 0.034138340146, -0.038880106034, 0.044603217363, 0.946016722396, 1.708172458034, 18.369114490772, 4.275967542224, 0.624449778826, 0.738801257357, 0.036936909247, 0.020437742859, 0.176759579388, -0.261128576436, -9.971904607075, -0.063549647738, 0, ],
        [0.007852964982, 0.003925745426, 0.287856349997, 58.053471054491, 0.030698062827, -0.000006837601, -0.000006123962, 0.000114549925, -17.580742026275, 0.55713614874, 0.205946900184, -43.230778067404, 0.004227082975, 0.006053854501, 0.006646690253, -0.009138926083, -0.248663457912, 3.325105302428, -0.276578605231, -0.455150962257, -0.005997822569, -0.002738986905, 0.082855748293, 0.070563187482, -0.101597078067, 1.596654829885, 0.065879787896, 0.017442923517, 0, ],
        [0.011497315687, -0.012583019909, 13.848373855148, 22.28881517216, 0.042287331657, 0.000197558695, 0.000176939544, -0.003309689199, -1.742140233901, -5.959510415282, -11.333020298294, -14.216479234895, -3.944800806497, 0.001304578929, -0.005139259078, 0.08647432259, 2.589998222025, 0.358614863147, 2.970887395829, 4.39160430183, 0.111893402319, 0.140739944934, 0.003938671797, 0.001185537435, 1.045176603318, -0.037531801533, -1.744525005833, -0.972957942438, 0, ],
        [-16.939142002537, 0.618053512295, 107.92089190414, 204.524147386814, 0.204407545189, 0.004742101706, 0.004247169746, -0.079444150933, -2.048456967261, -0.931989524708, -66.540858220883, -116.470289129818, -0.561301215495, -0.022312225275, -0.019484747345, 0.243518778973, 4.462098610572, -3.839389874682, 0.285714413078, 1.40736916669, 0.007176864388, 0.007131419303, -0.097503691021, -0.083171197416, 0.993922379938, -1.819432085819, -0.042893874898, -0.055015718216, 0, ],
        [-0.542809857455, 7.081822285872, -135.012404429101, 460.929268260027, 0.036498617908, 0.006937238413, 0.006213200589, -0.116219147061, -0.827454697348, 19.622217613195, 78.553728334274, -283.23862765888, 3.065444785639, -0.003847616297, -0.028984525722, 0.187507140282, 2.220506417769, 0.245737625222, 3.99902408961, -2.362524402134, 0.124769923797, 0.142065016461, 0.007538727793, 0.004297097528, -0.129475392736, -0.050587718062, -1.998394759416, 1.072835822585, 0, ],
        [-1.286456393795, 0.142279456389, -1.265748910581, 65.74306027738, -1.320702989799, -0.061855995532, -0.055400100872, 1.036269854556, -4.531489334771, 0.368539277612, 0.002487097952, -42.326462719738, 8.96223401238, 0.255676968878, 0.215513465742, -4.275436802385, -81.833676543035, 0.555500345288, 0.612894852362, -17.351836610113, -0.005925968725, -0.078209662789, 0.049750119549, 0.044979645917, -17.190711833803, -0.206830688253, -0.638945907467, 0.686150823668, 0, ],
        [0, 0, 0, 0, 0, -0.009702263896, -0.008689641059, 0.162541456323, 0, 0, 0, 0, 0, 0, 0, 0, -0.012, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ],
        [-8.153162937544, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.005, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ],
        [0, -3.261265175018, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.005, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ],
        [0, 0, 0, 0.17441246156, -3.261265175018, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.01, 0, 0, 0, 0, 0, 0, 0, 0, 0, ],
        [0, 0, -3.261265175018, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -8.5, -18, 0, 0, 0, 0, 0, 0, 0, ],
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, ],
        [0, 0, 0, -8.153162937544, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -8.5, -18, 0, 0, 0, 0, 0, ],
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, ],
        [0, 0, 0, 0, 0, 0, 0, 0, 0.699960862226, 0.262038222227, 0.159589891262, 0.41155156501, -1.701619176699, -0.0427567124, -0.038285155304, 0.703045934017, 16.975651534025, -0.115788018654, -0.127109026104, 3.599544290134, 0.001229743857, 0.016223661959, -0.01033400498, -0.00934235613, -6.433934989563, 0.042639567847, 0.132540852847, -0.142338323726, 0, ],
        [0, 0, 0, 0, 0, 0, 0, 0, -37.001496211974, 0.783588795613, -0.183854784348, -11.869599790688, -0.106084318011, -0.026306590251, -0.027118088888, 0.036744952758, 0.76460150301, 7.002366574508, -0.390318898363, -0.642631203146, -0.005701671024, 0.003522251111, 0.173867535377, 0.147911422248, 0.056092715216, -6.641979472328, 0.039602243105, 0.026181724138, 0, ],
        [0, 0, 0, 0, 0, 0, 0, 0, 1.991401999957, 13.760045912368, 2.53041689113, 0.082528789604, 0.728264862053, 0.023902766734, -0.022896554363, 0.015327568208, 0.370476566397, -0.412566245022, -6.70094564846, -1.327038338854, -0.227019235965, -0.267482033427, -0.008650986307, -0.003394359441, 0.098792645471, 0.197714179668, -6.369398456151, -0.011976840769, 0, ],
        [0, 0, 0, 0, 0, 0, 0, 0, 1.965859332057, -3.743127938662, -1.962645156793, 0.018929412474, 11.145046656101, -0.03600197464, -0.001222148117, 0.602488409354, 11.639787952728, -0.407672972316, 1.507740702165, -12.799953897143, 0.005393102236, -0.014208764492, -0.000915158115, -0.000640326416, -0.03653528842, 0.012458973237, -0.083125038259, -5.472831842357, 0, ],
        [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ],
        ])
    , 
        npmatrix([1.0 if d == 28 else 0.0 for d in xrange(29)])
    )



if __name__ == "__main__":
    main()
</code></pre>

<h3>Example Output</h3>

<p>Here's an example of the output that demonstrates the problem (each run is slightly different). Notice the large variability in execution times (over two orders of magnitude!). Again, this all goes away if I either use a pool of size 1 (or run the code without a pool), or if I don't use an explicit jacobian in the call to <code>integrate</code>.</p>

<blockquote>
  <ol>
  <li>times: 5847ms, 5760ms</li>
  <li>times: 4177ms, 3991ms</li>
  <li>times: 229ms, 36ms</li>
  <li>times: 1317ms, 1544ms</li>
  <li>times: 87ms, 100ms</li>
  <li>times: 113ms, 102ms</li>
  <li>times: 4747ms, 5077ms</li>
  <li>times: 597ms, 48ms</li>
  <li>times: 9ms, 49ms</li>
  <li>times: 135ms, 109ms</li>
  </ol>
</blockquote>

<p>This is intended as a formatted comment regarding the mathematical background raised in <a href="http://stackoverflow.com/questions/39350557/running-scipy-integrate-ode-in-multiprocessing-pool-results-in-huge-performance#comment66277221_39350557">a comment by @Dietrich</a>. As it doesn't address the programming question, I intend to delete this answer in a little while until the bounty blows over.</p>

<p>As @Dietrich noted, you can solve your ODE exactly, since if</p>

<pre><code>x' = A*x,
</code></pre>

<p>then the exact solution is</p>

<pre><code>x(t) = exp(A*t)*x0
</code></pre>

<p>Already I'd say that an exact solution is always superior than a numerical approximation, but this can indeed be faster than a numerical integration. As you noted in a comment, you're worried about efficiency. So don't compute the matrix exponential for each <code>t</code>: compute the eigensystem of <code>A</code> only once:</p>

<pre><code>A*v_i = L_i*v_i
</code></pre>

<p>then</p>

<pre><code>x(t) = sum_i c_i*v_i*exp(L_i*t),
</code></pre>

<p>and the coefficients <code>c_i</code> can be determined from the linear equations</p>

<pre><code>x0 = sum_i c_i*v_i.
</code></pre>

<p>Now, having an inhomogeneous term doesn't change much, as long as your matrix is not singular:</p>

<pre><code>x' = A*x + b
(x - A^(-1)*b)' = A*(x - A^(-1)*b)
</code></pre>

<p>so we can solve the homogeneous equation for <code>y = x - A^(-1)*b</code> and in a final step recover <code>x = y + A^(-1)*b</code>.</p>

<p>This all works nicely while the matrix is regular, but in your specific case it's singular. But it turns out that this is due to your final dimension:</p>

<pre><code>&gt;&gt;&gt; np.linalg.det(A)
0.0
&gt;&gt;&gt; np.linalg.det(A[:-1,:-1])
1920987.0461154305
</code></pre>

<p>And also note that the final row of <code>A</code> is all zeros (this is the reason for the singularity of <code>A</code>). So the last dimension of <code>x</code> is constant (or changes linearly due to <code>b</code>).</p>

<p>I suggest eliminating this variable, rewriting your equation for the rest of the variables, and solving the non-singular inhomogeneous linear system of ODEs using the above procedure, exactly. It should be faster and precise.</p>

<hr>

<p>The following will be a bit speculative, see also the caveat at the end.</p>

<p>In case of user-input <code>A</code> and <code>b</code>, things might get trickier. Finding a zero row/column in your matrix would be easy, but <code>A</code> can be singular even though none of its rows/columns are fully zero. I'm not an expert in the subject, but I think your best bet is using something akin to <a href="https://en.wikipedia.org/wiki/Principal_component_analysis" rel="nofollow">principal component analysis</a>: transforming your system of equations according to the eigensystem of <code>A</code>. My following thoughts will still assume that <code>A</code> is diagonalizable, but mostly because I'm unfamiliar with <a href="https://en.wikipedia.org/wiki/Singular_value_decomposition" rel="nofollow">singular value decomposition</a>. In realistic cases I'd expect your matrices to be diagonalizable, even if singular.</p>

<p>So I'll assume that the matrix <code>A</code> can be decomposed as</p>

<pre><code>A = V * D * V^(-1),
</code></pre>

<p>where <code>D</code> is a diagonal matrix containing the eigenvalues of <code>A</code>, and columns of <code>V</code> are the eigenvectors of <code>A</code> corresponding to each respective eigenvalue. The very same decomposition can be obtained in numpy using</p>

<pre><code>DD,V = np.linalg.eig(A)
D = np.asmatrix(np.diag(DD))
</code></pre>

<p>I usually prefer using <code>ndarray</code>s instead of matrices, but this way <code>V*D*np.linalg.inv(V)</code> would really correspond to the matrix product of the three matrices, rather than calling <code>np.dot</code> twice.</p>

<p>Now, rewrite your equation again:</p>

<pre><code>x' = A*x + b
x' = V*D*V^(-1)*x + b
V^(-1)*x' = D*V^(-1)*x + V^(-1)*b
</code></pre>

<p>By defining the auxiliary variables</p>

<pre><code>X = V^(-1)*x
B = V^(-1)*b
</code></pre>

<p>we obtain</p>

<pre><code>X' = D*X + B
</code></pre>

<p>i.e. the usual inhomogeneous form, but now <code>D</code> is a diagonal matrix containing the eigenvalues of <code>A</code> in the diagonal.</p>

<p>Since <code>A</code> is singular, some of the eigenvalues are zero. Look for zero elements in <code>D</code> (well, you can do that already with <code>DD</code> from <code>eig()</code>), and you'll know that they behave trivially during time-evolution. The remaining variables behave well, although at this point we see that the equations for <code>X</code> are decoupled due to <code>D</code> being diagonal, so you could integrate each independently and analytically. For this you need to first go from your initial condition <code>x0</code> to <code>X0 = np.linalg.inv(V)*x0</code>, then after solving the equations, back to <code>x = V*X</code>.</p>

<p>Caveat: as I said, I'm not an expert in this subject. I can easily imagine that the inversions involved in the diagonalization can be a numerical issue in practical applications. So I'd first test if the matrix is singular, and only carry on with this procedure if it is (or nearly is). It's possible that the above carries a lot of error, in which case numerical integration might be better (I really can't tell).</p>

python numpy scipy python-multiprocessing