数字集成Python与Matlab

Question

我的python代码大约需要6.2秒才能运行。 Matlab代码在0.05秒内运行。为什么这样做以及我可以做些什么来加速Python代码？ Cython是解决方案吗？

Matlab的：

function X=Test

nIter=1000000;
Step=.001;
X0=1;

X=zeros(1,nIter+1); X(1)=X0;

tic
for i=1:nIter
    X(i+1)=X(i)+Step*(X(i)^2*cos(i*Step+X(i)));
end
toc

figure(1) plot(0:nIter,X)

蟒蛇：

nIter = 1000000
Step = .001
x = np.zeros(1+nIter)
x[0] = 1
start = time.time()
for i in range(1,1+nIter):
      x[i] = x[i-1] + Step*x[i-1]**2*np.cos(Step*(i-1)+x[i-1])
end = time.time()
print(end - start)

Answer 1

How to speed up your Python code

你最大的时间接收器是np.cos，它对输入的格式进行了几次检查。这些是相关的，通常可以忽略不计的高维输入，但对于你的一维输入，这成为瓶颈。对此的解决方案是使用math.cos，它只接受一维数字作为输入，因此更快（虽然不太灵活）。

另一个时间汇是多次索引x。您可以通过更新一个状态变量来加快速度，并且每次迭代只写一次x。

有了这一切，你可以将事情加速大约十倍：

import numpy as np
from math import cos

nIter = 1000000
Step = .001
x = np.zeros(1+nIter)
state = x[0] = 1
for i in range(nIter):
    state += Step*state**2*cos(Step*i+state)
    x[i+1] = state

现在，你的主要问题是你真正的内心循环完全发生在Python中，即你有很多包装操作耗费时间。您可以通过使用uFuncs（例如，使用SymPy的ufuncify创建）和使用NumPy的accumulate来避免这种情况：

import numpy as np
from sympy.utilities.autowrap import ufuncify
from sympy.abc import t,y
from sympy import cos

nIter = 1000000
Step = 0.001
state = x[0] = 1
f = ufuncify([y,t],y+Step*y**2*cos(t+y))

times = np.arange(0,nIter*Step,Step)
times[0] = 1
x = f.accumulate(times)

这实际上在一瞬间运行。

… and why that’s not what you should worry about

如果您确切的代码（并且只是那个）是您关心的，那么您不应该担心运行时，因为它无论如何都非常短。另一方面，如果使用它来衡量运行时相当大的问题的效率，那么您的示例将失败，因为它只考虑一个初始条件并且是一个非常简单的动态。

此外，您使用的是Euler方法，根据您的步长，它不是非常有效或强大。后者（Step）在你的情况下是荒谬的低，产生的数据比你可能需要的多得多：步长为1，你可以看到发生了什么。

如果你想在这种情况下进行强大的集成，那么使用现代自适应积分器几乎总是最好的，它可以自己调整步长，例如，这里是使用本机Python集成器解决问题的方法：

from math import cos
import numpy as np
from scipy.integrate import solve_ivp

T = 1000
dt = 0.001

x = solve_ivp(
        lambda t,state: state**2*cos(t+state),
        t_span = (0,T),
        t_eval = np.arange(0,T,dt),
        y0 = [1],
        rtol = 1e-5
    ).y

这会自动将步长调整为更高的值，具体取决于容错rtol。它仍然返回相同数量的输出数据，但这是通过解决方案的插值。它对我来说以0.3秒的速度运行。

How to speed up things in a scalable manner

如果你仍然需要加速这样的事情，很可能你的衍生物（f）比你的例子复杂得多，因此它是瓶颈。根据您的问题，您可以对其计算进行矢量化（使用NumPy或类似方法）。

如果你不能矢量化，我写了一个module专门关注这个，通过硬编码你的衍生物。以下是采样步骤为1的示例。

import numpy as np
from jitcode import jitcode,y,t
from symengine import cos

T = 1000
dt = 1

ODE = jitcode([y(0)**2*cos(t+y(0))])
ODE.set_initial_value([1])
ODE.set_integrator("dop853")
x = np.hstack([ODE.integrate(t) for t in np.arange(0,T,dt)])

这在一瞬间再次运行。虽然这可能不是相关的速度提升，但这可以扩展到庞大的系统。

Answer 2

区别在于jit-compilation，Matlab默认使用它。让我们用Numba（Python jit-compiler）尝试你的例子

码

import numba as nb
import numpy as np
import time

nIter = 1000000
Step = .001

@nb.njit()
def integrate(nIter,Step):
  x = np.zeros(1+nIter)
  x[0] = 1
  for i in range(1,1+nIter):
    x[i] = x[i-1] + Step*x[i-1]**2*np.cos(Step*(i-1)+x[i-1])
  return x

#Avoid measuring the compilation time,
#this would be also recommendable for Matlab to have a fair comparison
res=integrate(nIter,Step)

start = time.time()
for i in range(100):
  res=integrate(nIter,Step)

end=time.time()
print((end - start)/100)

这导致每次调用0.022s运行时。