Юлия Дифференциал Equations.jl скорость

Новичок в Юлии, пытающейся проверить скорость решения ODE. Я использовал уравнение Лоренца в учебнике

using DifferentialEquations
using Plots
function lorenz(t,u,du)
du[1] = 10.0*(u[2]-u[1])
du[2] = u[1]*(28.0-u[3]) - u[2]
du[3] = u[1]*u[2] - (8/3)*u[3]
end

u0 = [1.0;1.0;1.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob,reltol=1e-8,abstol=1e-8,saveat=collect(0:0.01:100))

Вначале загрузка пакетов заняла около 25 секунд, а код работал в течение 7 секунд на четырехъядерном ноутбуке с Windows 10 в ноутбуке Jupyter. Я понимаю, что Юле нужно предварительно скомпилировать пакеты, и поэтому время загрузки было таким долгим? Я нашел 25 с невыносимым. Кроме того, когда я снова запустил солвер, используя другие начальные значения, это заняло гораздо меньше времени (~1 с), и почему? Это типичная скорость?

1 ответ

Решение

Tl; др:

  1. Пакеты Julia имеют фазу предварительной компиляции. Это помогает сделать все дальше using звонки быстрее, за счет первого, хранящего некоторые данные компиляции. Это срабатывает только при каждом обновлении пакета.
  2. using должен вытащить пакет, который занимает немного (в зависимости от того, сколько можно предварительно скомпилировать).
  3. Прекомпиляция не является "полной", поэтому при первом запуске функции, даже из пакета, она должна будет скомпилироваться.
  4. Разработчики Julia знают об этом, и уже есть планы избавиться от (2) и (3), сделав прекомпиляцию более полной. Также есть планы по сокращению времени компиляции, о которых я не знаю подробностей.
  5. Все функции Julia специализируются на данных типах, и каждая функция является отдельным типом, поэтому внутренние функции DiffEq специализируются на каждой предоставляемой вами функции ODE.
  6. В большинстве случаев с длинными вычислениями (5) на самом деле не имеет значения, так как вы не меняете функции так часто (если это так, рассмотрите возможность изменения параметров).
  7. Но (6) имеет значение при использовании его в интерактивном режиме. Это заставляет это чувствовать себя менее "гладким".
  8. Мы можем избавиться от этой специализации в функции ODE, но она не используется по умолчанию, потому что она вызывает попадание в 2x-4x. Возможно, это будет по умолчанию в будущем.
  9. Наши сроки после предварительной компиляции по этой проблеме все еще лучше, чем такие вещи, как завернутые в SciPy решатели Fortran для подобных задач в 20 раз. Так что это проблема времени компиляции, а не проблема времени выполнения. Время компиляции по существу постоянное (большие проблемы, вызывающие одну и ту же функцию, имеют примерно одинаковую компиляцию), так что это действительно просто проблема интерактивности.
  10. Мы (и Юлия в целом) можем и будем лучше работать с интерактивностью в будущем.

Полное объяснение

Это на самом деле не вещь DifrentialEquations.jl, это просто пакет Julia. 25 с должны включать время предварительной компиляции. При первой загрузке пакета Julia он прекомпилируется. Тогда это не должно повториться до следующего обновления. Вероятно, это самая длинная инициализация, и она довольно длинная для diffrentialEquations.jl, но, опять же, это происходит только каждый раз, когда вы обновляете код пакета. Затем каждый раз, когда есть небольшая стоимость инициализации для using, DiffEq довольно большой, поэтому для его инициализации требуется немного:

@time using DifferentialEquations
5.201393 seconds (4.16 M allocations: 235.883 MiB, 4.09% gc time)

Тогда, как отмечено в комментариях, у вас также есть:

@time using Plots
6.499214 seconds (2.48 M allocations: 140.948 MiB, 0.74% gc time)

Затем в первый раз вы запускаете

function lorenz(t,u,du)
  du[1] = 10.0*(u[2]-u[1])
  du[2] = u[1]*(28.0-u[3]) - u[2]
  du[3] = u[1]*u[2] - (8/3)*u[3]
end

u0 = [1.0;1.0;1.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)
@time sol = solve(prob,reltol=1e-8,abstol=1e-8,saveat=collect(0:0.01:100))

6.993946 seconds (7.93 M allocations: 436.847 MiB, 1.47% gc time)

Но потом второй и третий раз

0.010717 seconds (72.21 k allocations: 6.904 MiB)
0.011703 seconds (72.21 k allocations: 6.904 MiB)

Так что здесь происходит? В первый раз, когда Юлия запускает функцию, она скомпилирует ее. Итак, в первый раз вы запускаете solve, он скомпилирует все свои внутренние функции во время работы. Все продолжающиеся времена будут без компиляции. DiffrentialEquations.jl также специализируется на самой функции, поэтому, если мы изменим функцию:

function lorenz2(t,u,du)
  du[1] = 10.0*(u[2]-u[1])
  du[2] = u[1]*(28.0-u[3]) - u[2]
  du[3] = u[1]*u[2] - (8/3)*u[3]
end

u0 = [1.0;1.0;1.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz2,u0,tspan)

мы возьмем время компиляции снова:

@time sol = 
solve(prob,reltol=1e-8,abstol=1e-8,saveat=collect(0:0.01:100))
3.690755 seconds (4.36 M allocations: 239.806 MiB, 1.47% gc time)

Так вот что, теперь почему. Здесь есть несколько вещей вместе. Прежде всего, пакеты Julia не полностью прекомпилируются. Они не хранят кэшированные скомпилированные версии реальных методов между сессиями. Это то, что нужно сделать в списке выпусков 1.x, и это избавит от этого первого попадания, подобно простому вызову пакета C/Fortran, так как он будет просто скомпилирован с большим опережением времени (AOT) функции. Так что это будет хорошо, но сейчас просто отметьте, что есть время запуска.

Теперь поговорим об изменении функций. Каждая функция в Julia автоматически специализируется на своих аргументах ( подробности см. В этом блоге). Ключевая идея здесь заключается в том, что каждая функция в Юлии является отдельным конкретным типом. Таким образом, поскольку тип проблемы здесь параметризован, изменение функции вызывает компиляцию. Обратите внимание, что это отношение: вы можете изменить параметры функции (если у вас были параметры), вы можете изменить начальные условия и т. Д., Но только изменение типа вызывает перекомпиляцию.

Стоит ли оно того? Хорошо, может быть. Мы хотим специализироваться на быстрых вычислениях, которые сложны. Время компиляции является постоянным (то есть вы можете решить ODE за 6 часов, и это все равно будет несколько секунд), и поэтому вычислительно-дорогостоящие вычисления здесь не выполняются. Моделирование по методу Монте-Карло, когда вы запускаете тысячи параметров и начальных условий, здесь не выполняется, потому что если вы просто изменяете значения начальных условий и параметров, то оно не будет перекомпилировано. Но при интерактивном использовании, когда вы меняете функции, вы получаете около секунды, что нехорошо. Один из ответов от разработчиков Julia на это - потратить время после Julia 1.0 на ускорение компиляции, о чем я не знаю деталей, но я уверен, что здесь есть какой-то низко висящий фрукт.

Можем ли мы избавиться от этого? Да. DiffEq Online не перекомпилируется для каждой функции, потому что он ориентирован на онлайн-использование.

function lorenz3(t,u,du)
  du[1] = 10.0*(u[2]-u[1])
  du[2] = u[1]*(28.0-u[3]) - u[2]
  du[3] = u[1]*u[2] - (8/3)*u[3]
  nothing
end

u0 = [1.0;1.0;1.0]
tspan = (0.0,100.0)
f = NSODEFunction{true}(lorenz3,tspan[1],u0)
prob = ODEProblem{true}(f,u0,tspan)

@time sol = solve(prob,reltol=1e-8,abstol=1e-8,saveat=collect(0:0.01:100))

1.505591 seconds (860.21 k allocations: 38.605 MiB, 0.95% gc time)

И теперь мы можем изменить функцию и не брать на себя стоимость компиляции:

function lorenz4(t,u,du)
  du[1] = 10.0*(u[2]-u[1])
  du[2] = u[1]*(28.0-u[3]) - u[2]
  du[3] = u[1]*u[2] - (8/3)*u[3]
  nothing
end

u0 = [1.0;1.0;1.0]
tspan = (0.0,100.0)
f = NSODEFunction{true}(lorenz4,tspan[1],u0)
prob = ODEProblem{true}(f,u0,tspan)

@time sol = 
solve(prob,reltol=1e-8,abstol=1e-8,saveat=collect(0:0.01
:100))

0.038276 seconds (242.31 k allocations: 10.797 MiB, 22.50% gc time)

И тада, оборачивая функцию в NSODEFunction (который внутренне использует FunctionWrappers.jl), он больше не специализируется для каждой функции, и вы попадаете на время компиляции один раз за сеанс Julia (и затем один раз, когда он кэшируется, один раз за обновление пакета). Но обратите внимание, что это примерно в два-четыре раза дороже, поэтому я не уверен, будет ли он включен по умолчанию. Мы могли бы сделать это по умолчанию внутри конструктора проблемного типа (то есть без дополнительной специализации по умолчанию, но пользователь может выбрать более высокую скорость за счет интерактивности), но я не уверен, что здесь лучше по умолчанию (не стесняйтесь прокомментируйте проблему своими мыслями). Но это определенно будет задокументировано вскоре после того, как Джулия изменит свой аргумент ключевого слова, и поэтому режим "без компиляции" будет стандартным способом его использования, даже если он не используется по умолчанию.

Но просто, чтобы поставить это в перспективе,

import numpy as np
from scipy.integrate import odeint
y0 = [1.0,1.0,1.0]
t = np.linspace(0, 100, 10001)
def f(u,t):
    return [10.0*(u[1]-u[0]),u[0]*(28.0-u[2])-u[1],u[0]*u[1]-(8/3)*u[2]]
%timeit odeint(f,y0,t,atol=1e-8,rtol=1e-8)

1 loop, best of 3: 210 ms per loop

мы смотрим, нужно ли сделать это интерактивное удобство по умолчанию, чтобы оно было в 5 раз быстрее, а не в 20 раз быстрее, чем значение по умолчанию в SciPy (хотя наше значение по умолчанию будет гораздо более точным, чем в SciPy по умолчанию, но это данные для другого времени, которые найти в тестах или просто спросить). С одной стороны, это имеет смысл как простота использования, но с другой стороны, если повторное включение специализации для длинных вычислений и Монте-Карло неизвестно (где вы действительно хотите скорость), тогда многие люди будут там принять 2x-4x удар производительности, который может составить дополнительные дни / недели вычислений. Эхх... тяжелый выбор.

Таким образом, в конце концов, у Джулии отсутствует смесь оптимизирующих вариантов и некоторых функций предварительной компиляции, которые влияют на интерактивность, не влияя на истинную скорость выполнения. Если вы хотите оценить параметры с помощью некоторого большого Монте-Карло, или решить тонну SDE, или решить большой PDE, у нас это не так. Это была наша первая цель, и мы постарались поразить ее как можно лучше. Но у игры в REPL есть 2-3 секунды "глюков", которые мы также не можем игнорировать (лучше, чем играть в C/Fortran, хотя, конечно, но все же не идеально для REPL). Для этого я показал вам, что решения уже разрабатываются и тестируются, и, надеюсь, на этот раз в следующем году у нас будет лучший ответ для этого конкретного случая.

PS

Еще две вещи, на которые стоит обратить внимание. Если вы используете только решатели ODE, вы можете просто сделать using OrdinaryDiffEq продолжать загружать / устанавливать / компилировать / импортировать все DifrentialEquations.jl ( это описано в руководстве). Кроме того, используя saveat похоже, что это не самый быстрый способ решить эту проблему: решение с гораздо меньшим количеством точек и использование плотного вывода при необходимости может быть лучше здесь.

редактировать

Я открыл вопрос о том, как мы можем сократить время компиляции между функциями без потери ускорения, которое дает специализация. Я думаю, что это то, что мы можем сделать краткосрочным приоритетом, так как я согласен, что мы могли бы добиться большего успеха здесь.

Другие вопросы по тегам