Как улучшить скорость odeint в Python?

Question

Я использую Python и odeint из пакета scipy для решения большого числа (~ 10e6) связанных ODE. Систему уравнений можно сформулировать как сумму некоторых матричных умножений, и я использую numpy с поддержкой blas для этого. Моя проблема в том, что это занимает очень много времени. Когда я профилирую код, я вижу, что большую часть времени уходит в odeint, делая что-то еще, кроме оценки правых. Это пять самых трудоемких звонков от профилировщика:

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 ужасно масштабируется с количеством уравнений. У вас есть идеи, почему это так и как улучшить производительность?

Спасибо, что уделили время

Оскар Акерлунд

Answer 1

Это как минимум один из возможных источников количества времени:

Если вы не предоставите якобиан для odeint (т.е. LSODA), он попытается вычислить его с помощью конечных разностей. Более того, он может попытаться инвертировать якобиан, который масштабируется как O (m ^ 3), если он считает, что проблема жесткая. Оба этих шага дороги, когда число переменных велико.

Вы можете попытаться сократить время, затрачиваемое этими операциями, заставив odeint использовать полосатый якобиан, передав в подпрограмму подходящие значения параметров ml и mu. Вам не нужно указывать Dfun, эти параметры также применяются к якобиану, вычисленному дифференцированием.

Answer 2

Десять миллионов уравнений - это нетривиальное число.

Почему вы говорите, что это "ужасно масштабируется"? Матричное сложение составляет O(m^2) для m x m матриц, а умножение - O(m^3). Вы указываете две точки настенного времени в зависимости от числа уравнений / степеней свободы, но это может описать только прямую линию. Я бы выбрал пару промежуточных точек между 4K и 10M и посмотрел, показывает ли обозначение Big-Oh, как оно масштабируется. Подгоните результаты к полиному 3-го порядка времени стены по сравнению с DOF; это скажет вам, как все масштабируется.

Вы уравнения линейные или нелинейные? Статический или переходный? В зависимости от типа проблемы вы можете играть с некоторыми другими параметрами, такими как размер временного шага, критерии сходимости, выбор схемы интеграции и т. Д.