julia - Julia 微分方程.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 s难以忍受。另外，当我使用不同的初始值再次运行求解器时，它花费的时间要少得多（〜1秒），为什么？这是典型的速度吗？

Answer 1

Tl; dr：

Julia程序包有一个预编译阶段。这有助于使所有进一步的using调用更快，而第一个调用会存储一些编译数据。这仅在每次软件包更新时触发。
using必须提取需要一点点的软件包（取决于可以预编译的数量）。
预编译不是“完整的”，因此，即使您是从程序包中首次运行函数，也必须对其进行编译。
朱莉娅（Julia）开发人员知道这一点，并且已经计划通过使预编译更加完整来摆脱（2）和（3）。还计划减少编译时间，我不知道这方面的细节。
Julia的所有函数都专门处理给定的类型，每个函数都是单独的类型，因此DiffEq的内部函数专门处理您提供的每个ODE函数。
在大多数情况下，计算时间长，（5）实际上并不重要，因为您不经常更改函数（如果需要的话，请考虑更改参数）。
但是（6）在交互使用时确实很重要。它使它不那么“平滑”。
我们可以摆脱ODE函数的这种专长，但这不是默认设置，因为它会导致2x-4x的性能冲击。也许它将是将来的默认设置。
我们在此问题上进行预编译后的时机仍然比SciPy封装的Fortran解算器等问题好20倍。因此，这全是编译时间问题，而不是运行时问题。编译时间本质上是恒定的（调用同一函数的较大问题大约具有相同的编译），因此这实际上只是一个交互性问题。
我们（以及一般来说，朱莉娅）将来会并且将来会在互动方面做得更好。

完整说明
这确实不是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中的每个函数都会自动专注于其参数（ see this blog post for details）。这里的关键思想是Julia中的每个函数都是单独的具体类型。因此，由于此处的问题类型已参数化，因此更改函数会触发编译。请注意就是这种关系：您可以更改函数的参数（如果有参数），可以更改初始条件等，但这仅更改触发重新编译的类型。
这值得么？也许会。我们想专门处理快速的困难的事情。编译时间是恒定的（即，您可以求解6小时的ODE，并且仍然需要几秒钟的时间），因此此处不会进行计算量大的计算。在此处运行数千个参数和初始条件的蒙特卡洛模拟不会受到影响，因为如果您仅更改初始条件和参数的值，那么它将不会重新编译。但是，在要更改功能的地方进行交互使用确实能获得一秒钟左右的效果，但这并不好。 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会话中一次遇到编译时间（然后缓存一次，每次打包更新一次）。但是请注意，这大约是2倍至4倍的成本 so I am not sure if it will be enabled by default。我们可以默认在问题类型构造函数内部实现此目的（即默认情况下无需额外的专业化，但用户可以选择以交互性为代价提高速度），但是我不确定这里有更好的默认值（可以用您的想法对问题发表评论）。但是它一定会在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所使用的默认设置准确得多，但这是另一时间的数据，可以在基准中找到或直接问）。一方面，这很容易使用，但另一方面，如果重新启用长时间计算的专业化功能，而蒙特卡洛不为人所知（这是您真正想要的速度），那么那里的很多人都会造成2到4倍的性能下降，这可能需要额外的天数/周的计算。嗯...艰难的选择。
因此，最终，Julia缺少了优化选择和一些预编译功能的混合，这些功能会影响交互性而不会影响真正的运行速度。如果您要使用一些较大的Monte Carlo估计参数，或者要解决大量的SDE，或者要解决较大的PDE，那么我们不建议这样做。那是我们的第一个目标，我们确保做到最好。但是在REPL中玩游戏确实有2-3秒的“故障”，我们也不能忽略（当然，比在C / Fortran中玩游戏更好，但是对于REPL来说仍然不理想）。为此，我向您展示了已经开发和测试了解决方案，因此希望明年明年这个时候我们可以针对该特定情况提供更好的答案。
聚苯乙烯
还有两件事要注意。如果仅使用ODE求解器，则只需执行 using OrdinaryDiffEq即可继续下载/安装/编译/导入所有DifferentialEquations.jl（ this is described in the manual）。同样，使用 saveat可能不是解决此问题的最快方法：用更少的点来解决它，并在需要时使用密集输出可能更好。
编辑
I opened an issue detailing how we can reduce the "between function" compilation time without losing the speedup that specializing gives。我认为这是我们可以短期优先考虑的事情，因为我同意我们可以在这里做得更好。