如何提高Python中odeint的速度？

Question

我正在使用Python和scipy包中的odeint来解决大量（约10e6）耦合ODE。方程组可以表示为一些矩阵乘法的总和，我使用numpy和blas支持进行计算。我的问题是这需要很长时间。当我对代码进行分析时，我发现大部分时间都花在odeint上，而不是评估rhs。以下是分析器中消耗时间最多的五个调用：

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
5       1547.915  309.583 1588.170  317.634 {scipy.integrate._odepack.odeint}
60597   11.535    0.000   23.751    0.000   terms3D.py:5(two_body_evolution)
121194  11.242    0.000   11.242    0.000   {numpy.core._dotblas.dot}
60597   10.145    0.000   15.460    0.000   generator.py:13(Gs2)
121203   3.615    0.000   3.615     0.000   {method 'repeat' of 'numpy.ndarray' objects}

rhs基本上由two_body_evolution和Gs2组成。这个配置文件包含大约7000个耦合的ODE，而这是大约4000个ODE的相同情况：

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
5        259.427  51.885   273.316  54.663 {scipy.integrate._odepack.odeint}
30832    3.809    0.000    7.864    0.000  terms3D.py:5(two_body_evolution)
61664    3.650    0.000    3.650    0.000  {numpy.core._dotblas.dot}
30832    3.464    0.000    5.637    0.000  generator.py:13(Gs2)
61673    1.280    0.000    1.280    0.000  {method 'repeat' of 'numpy.ndarray' objects}

我的主要问题是，odeint中的“隐藏”时间随着方程数量的增加而急剧增加。你有什么想法，为什么会这样，并且如何提高性能？

感谢您的时间。

Oscar Åkerlund

- Oscar Åkerlund

Answer 1

以下是至少一种计算时间的可能来源：

如果您没有向odeint（即LSODA）提供Jacobian，它会尝试通过有限差分来计算。此外，如果它认为问题很棘手，它可能会尝试反转Jacobian，它的缩放比例为O(m^3)。当变量数量较大时，这两个步骤都是昂贵的。

您可以通过强制odeint使用带状Jacobian来尝试减少这些操作所需的时间，通过向该程序传递适当的ml和mu参数值。您不需要提供Dfun，这些参数也适用于通过微分计算的Jacobian。

- pv.

---

在我有限的经验中，如果没有提供雅可比矩阵，请不要考虑使用ODEPACK积分器进行非平凡计算。 - timday

Answer 2

一千万个方程是一个不小的数字。为什么说它“扩展性很差”？矩阵加法对于m x m矩阵是O(m^2)，乘法是O(m^3)。您引用了时间和方程数/自由度的两个点，但这只能描述一条直线。我会选择4K到10M之间的几个中间点，看看大O符号如何显示它的扩展性。将结果拟合为壁钟时间与DOF的三次多项式；这将告诉您事物如何扩展。 您的方程是线性还是非线性？静态还是瞬态？根据问题的类型，您可能可以玩一些其他参数，比如时间步长、收敛标准、积分方案选择等。

- duffymo

---

我已经将时间适应于m（mxm矩阵），但它的规模略差于m^6，而我原本期望的是类似于m^3或者m^4的结果，因为这涉及到矩阵乘法。 - Oscar Åkerlund

这些方程是线性和瞬态的，我也尝试降低收敛容限，但它并没有显著影响。 - Oscar Åkerlund