朱莉娅 DifferentialEquations.jl 速度

8

我正在尝试测试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秒太长了。而且,当我使用不同的初始值再次运行求解器时,它所需的时间要少得多(约1秒),这是为什么?这是典型的速度吗?


DifferentialEquations是预编译的,因此25秒绝对不对。可能是Plots,但即使如此,25秒的时间对我来说听起来也太长了(但我在Ubuntu上,不确定在Windows 10上绘图是什么样子)。关于第二次更短的运行时间,那是因为没有编译,这是可以预料的。 - Colin T Bowers
可能是花了25秒钟来加载库(这仍然不应该需要那么长时间,但你正在使用Windows)?如果您更改微分方程而不仅仅是初始条件,求解器的第二次运行需要多长时间? - Wrzlprmft
是的,Plots需要很长时间才能加载。 - Michael K. Borregaard
“Plots”在我的工作机器上需要> 10秒钟才能加载,而在旧机器上需要更长时间。 - Alexander Morley
1个回答

23

简而言之:

  1. Julia的包有预编译阶段。这可以使所有后续的using调用更快,代价是第一个调用需要存储一些编译数据。这仅在每个包更新时触发。
  2. using必须拉取包,这需要一点时间(取决于可以预编译多少内容)。
  3. 预编译不是“完整”的,所以第一次运行函数时,即使是来自包的函数,也需要进行编译。
  4. Julia开发人员知道这一点,已经计划通过使预编译更完整来消除(2)和(3)。还计划减少编译时间,但我不知道具体细节。
  5. 所有Julia函数都专门针对给定的类型,并且每个函数都是单独的类型,因此DiffEq的内部函数是针对您提供的每个ODE函数进行专门化的。
  6. 在大多数长时间计算的情况下,(5)实际上并不重要,因为您并不经常更改函数(如果您确实更改,请考虑更改参数而不是函数)。
  7. 但是,在交互使用时,(6)很重要。它会使操作感觉不那么“流畅”。
  8. 我们可以消除对ODE函数的专业化,但这不是默认设置,因为会导致2倍至4倍的性能损失。也许将来会成为默认设置。
  9. 我们在此问题上进行预编译后的时间仍比像SciPy的封装Fortran求解器等问题快20倍。因此,这只是一个编译时间问题,而不是运行时问题。编译时间本质上是恒定的(调用相同函数的更大问题具有大约相同的编译),因此这实际上只是一个交互性问题。
  10. 我们(以及Julia总体)未来可以做得更好。

完整解释

这不是 DifferentialEquations.jl 的问题,这只是 Julia 包的问题。25 秒可能包括预编译时间。第一次加载 Julia 包时会进行预编译。然后在下一次更新之前就不需要再次进行预编译。这可能是最长的初始化时间,并且对于 DifferentialEquations.jl 来说相当长,但是这只会在每次更新包代码时发生。然后,每次使用时都会有一个小的初始化成本。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)

所以这就是“what”,现在来说“why”。这里有几个因素。首先,Julia包没有完全预编译。它们不会在会话之间保留实际方法的缓存编译版本。这是1.x版本发布清单上要做的事情,这将消除第一次调用时的开销,类似于调用C/Fortran包,因为它只会命中许多提前编译(AOT)的函数。所以这很好,但现在请注意存在启动时间。
现在让我们谈谈如何更改这些函数。Julia中的每个函数都会自动专门处理其参数(详见此博客文章)。关键思想在于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倍的成本因此我不确定是否默认启用。我们可以在问题类型构造函数内部默认使其生效(即默认情况下不进行额外的专门化,但用户可以选择更快的速度以换取互动性),但我不确定这里哪种更好的默认值(欢迎在问题上发表您的想法)。但是,一旦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倍的性能损失,这可能会导致额外的计算天数/周数。嗯...抉择很艰难。
因此,最终存在一些优化选择和一些预编译功能缺失的混合因素会影响交互性,而不会影响真正的运行时速度。如果您想使用一些大型蒙特卡罗来估计参数,或解决大量SDE或大型PDE问题,我们已经掌握了这个技巧。这是我们的第一个目标,我们确保尽可能好地完成了它。但是,在REPL中玩耍确实有2-3秒钟的“故障”,我们也不能忽视它(当然比在C / Fortran中玩耍要好,但对于REPL来说仍然不是理想的)。为此,我向您展示了已经正在开发和测试的解决方案,希望明年的这个时候我们可以在这个特定情况下有一个更好的答案。
附言
请注意两件事。如果您仅使用ODE求解器,可以只执行using OrdinaryDiffEq,而不必下载/安装/编译/导入DifferentialEquations.jl(在手册中有描述)。另外,像那样使用saveat可能不是解决此问题的最快方法:使用更少的点进行求解,并根据需要使用密集输出可能更好。

编辑

我提出了一个问题,详细说明了我们如何在不失去专业化加速的情况下减少“函数之间”的编译时间。我认为这是我们可以在短期内优先考虑的事情,因为我同意我们可以做得更好。

谢谢!这真的很有帮助。所以基本上总运行时间是编译时间和计算时间的总和,而编译时间应该始终为几秒钟?如果是这样的话,我不介意花时间在这上面。在这种特殊情况下,问题非常简单,因此编译时间比计算时间长得多,我是对的吗? - JianghuiDu
是的,这就是情况。虽然只有几秒钟,但如果您只是快速地做一次什么事情,您会感觉到它。 - Chris Rackauckas
我的机器上预编译大约需要5到10分钟,之后每次加载库需要大约1分钟。我计时发现,分配的数量比你在上面的回答中报告的要多得多。
@time using DifferentialEquations 54.007511秒(49.86 M分配:2.870 GiB,3.24%gc时间)
我的Arch Linux设置有问题吗?
- lattitude
自此文写作以来,编译已经发生了三年半的变化。而 Julia v1.6 又是一个重大变革。到了 v1.6,预编译时间将再次降低近一个数量级,但您会注意到在 v1.3 中,所有内容都进行了更多的预缓存,因此第一次绘图的时间也大大减少了。 - Chris Rackauckas

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接