Julia DifferentialEquations.jl速度

Question

Julia的新手,试图测试ODE求解器的速度.我在教程中使用了Lorenz方程

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秒,代码在Jupyter笔记本电脑的Windows 10四核笔记本电脑上运行了7秒.我知道Julia需要先预编译包,这就是加载时间这么长的原因吗？我发现25秒难以忍受.此外,当我使用不同的初始值再次运行求解器时,运行时间(~1s)要少得多,为什么会这样？这是典型的速度吗？

Answer 1

TL;博士:

1. Julia包具有预编译阶段.这有助于using更快地进行所有进一步的调用,代价是第一个存储一些编译数据.这仅触发每个包更新.
2. using 必须拉入需要一点点的包(取决于可以预编译多少).
3. 预编译不是"完整",因此第一次运行函数时,即使从包中,也必须进行编译.
4. Julia devs知道这一点,并且已经计划通过使预编译更完整来摆脱(2)和(3).还有计划减少编译时间,我不知道有关详细信息.
5. 所有Julia函数都专注于给定的类型,每个函数都是一个单独的类型,因此DiffEq的内部函数专注于您提供的每个ODE函数.
6. 在大多数使用长计算的情况下,(5)实际上并不重要,因为您不经常更改函数(如果是,请考虑更改参数).
7. 但是(6)在交互使用它时很重要.它让它感觉不那么"平滑".
8. 我们可以在ODE函数上摆脱这种专业化,但它不是默认值,因为它会导致2x-4x高性能命中.也许这将是未来的默认.
9. 关于这个问题的预编译发布后的时间仍然比SciPy包装的Fortran解决方案更好,比如20x这样的问题.所以这都是编译时问题,而不是运行时问题.编译时间基本上是不变的(调用相同函数的大问题具有大致相同的编译),因此这实际上只是一个交互问题.
10. 我们(和Julia一般)可以并且将来会在交互性方面做得更好.

完整解释

这真的不是一个DifferentialEquations.jl的东西,这只是一个Julia包的东西.25s必须包括预编译时间.第一次加载Julia包时,它会预编译.然后在下次更新之前不需要再次发生.这可能是最长的初始化,对于DifferentialEquations.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)

那么这里发生了什么？Julia第一次运行一个函数,它将编译它.因此,第一次运行时solve,它将在运行时编译其所有内部函数.所有进行的时间都没有编译.DifferentialEquations.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中的每个函数都会自动专注于其参数(有关详细信息,请参阅此博客文章).这里的关键思想是Julia中的每个函数都是一个单独的具体类型.因此,由于此处的问题类型是参数化的,因此更改函数会触发编译.注意它是关系:您可以更改函数的参数(如果您有参数),您可以更改初始条件等,但它只是更改触发重新编译的类型.

这值得么？也许.我们希望专注于快速进行难以计算的事情.编译时间是恒定的(即你可以解决6小时的ODE并且它仍然是几秒钟),因此计算成本高昂的计算在这里不受影响.蒙特卡罗模拟在这里运行数以千计的参数和初始条件不受影响,因为如果您只是更改初始条件和参数的值,那么它将不会重新编译.但是你正在改变功能的交互式使用确实会在那里获得第二次左右,这并不好.朱莉娅开发人员对此的一个答案就是花费朱莉娅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会话都会达到一次编译时间(然后一旦缓存,每个包更新一次).但请注意,这大约是2x-4x的成本,所以我不确定它是否会默认启用.我们可以在问题类型构造函数中默认发生这种情况(即默认情况下没有额外的特殊化,但是用户可以以交互性为代价选择更快的速度)但是我不确定这里更好的默认值是什么(请随意用你的想法评论这个问题).但是,在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

我们正在研究这种交互式便利是否应该默认为比SciPy的默认速度快5倍,而不是20倍速(尽管我们的默认值通常比默认的SciPy使用的更准确,但这是另一个时间的数据,可以在基准测试中找到或只是问).一方面,它易于使用,但另一方面,如果重新启用长时间计算的专业化和蒙特卡洛不知道(这是你真正想要的速度),那么很多人将会采取2x-4x的性能命中,这可能相当于额外的几天/几周的计算.呃......艰难的选择.

因此,最终会有优化选择和Julia缺少的一些预编译功能,这些功能会影响交互性,而不会影响真正的运行时速度.如果您正在寻找使用一些大型蒙特卡罗估算参数,或者解决大量的SDE,或者解决一个大的PDE,那么我们就是这样.这是我们的第一个目标,我们确保尽可能地达到目标.但是在REPL中玩游戏确实有2-3秒的"gliches",我们也不能忽视它(当然比在C/Fortran中游玩更好,但对于REPL来说仍然不理想).为此,我已经向您展示了已经开发和测试的解决方案,所以希望明年的这个时候我们可以为这个具体案例找到更好的答案.

PS

还有两点需要注意.如果您只使用ODE解算器,您可以using OrdinaryDiffEq继续下载/安装/编译/导入所有DifferentialEquations.jl(这在手册中有描述).此外,使用saveat这样可能不是解决此问题的最快方法:用更少的点来解决它并在必要时使用密集输出可能更好.

编辑

我打开了一个问题,详细说明了如何减少"功能间"编译时间而不会失去专业化的加速.我认为这是我们可以作为短期优先事项的事情,因为我同意我们可以在这里做得更好.