如何提高Python中odeint的速度？

Question

我正在使用scipy包中的Python和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是相同的:

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中的"隐藏"时间与方程的数量可怕地成比例.您对此有何想法以及如何提高性能？

感谢您的时间

奥斯卡Åkerlund

Answer 1

这至少是一个可能的时间来源:

如果你不提供雅可比行列式odeint(即LSODA),它将尝试通过有限差分来计算它.此外,如果它认为问题是僵硬的,它可能会试图将雅可比变换为O(m ^ 3).当变量数量很大时,这两个步骤都很昂贵.

您可以尝试odeint通过将适当的值ml和mu参数传递给例程来强制使用带状雅可比行列式来减少这些操作所花费的时间.您不需要提供Dfun,这些参数也适用于通过微分计算的雅可比.

Answer 2

一千万个方程是一个不小的数字。

为什么说它“规模可怕”？矩阵加法适用O(m^2)于 mxm 矩阵，乘法适用于O(m^3). 您引用了墙时间与方程数/自由度的两个点，但这只能描述一条直线。我会在 4K 和 10M 之间选择几个中间点，看看 Big-Oh 表示法是否可以显示它的缩放方式。将结果拟合为壁时间与自由度的三阶多项式；这会告诉你事情是如何扩展的。

你的方程是线性的还是非线性的？静态还是瞬态？根据问题的类型，您可能可以使用一些其他参数，例如时间步大小、收敛标准、集成方案选择等。