Julia、MATLAB和Python的性能对比
导言
在选择用于工程计算的工具时，不仅方便性和库的数量很重要，而且计算的速度也很重要。
MATLAB被认为是科学计算的标准，Python由于NumPy和SciPy生态系统而变得广泛，而Julia是专门为高性能计算而创建的。
在这个例子中，我们将比较Julia，MATLAB和Python在几个典型工程任务上的性能。:
*矢量计算;
*矩阵运算;
*求解线性方程组;
*快速傅立叶变换（FFT）计算;
*数值积分;
*求数π的蒙特卡洛方法;
*求解热导率方程。
在所有测试中，我们将测量执行时间并比较结果。
为了客观起见，我们将使用内置工具测量时间。:
@belapsed在朱莉娅;
tic … toc 在MATLAB中;
time.perf_counter()在Python中。
连接图书馆
要运行测试，我们需要以下库:
MATLAB-从Julia调用MATLAB;PyCall-与Python的交互;LinearAlgebra-用于矩阵运算;FFTW-用于快速傅立叶变换;BenchmarkTools-精确的时间测量。
neededLibs = ["MATLAB", "PyCall", "BenchmarkTools", "LinearAlgebra", "FFTW"]
for lib in neededLibs
try
eval(Meta.parse("using $lib"))
catch ex
Pkg.add(lib)
eval(Meta.parse("using $lib"))
end
end
Python模块的准备
我们将连接必要的Python库。
np = pyimport("numpy")
scipy_integrate = pyimport("scipy.integrate")
测试1。 矢量计算
计算表达式
为1000万个元素的数组。
这个表达式被特别选择作为密集的零碎计算的一个例子.
N = 10^7
x = rand(N)
茱莉亚
julia_time = @belapsed sin.($x).^2 .+ cos.($x).^2
println("Julia: ", julia_time, " 秒数")
MATLAB
让我们将数据传输到MATLAB并测量执行时间。
@mput x
mat"""
tic;
y = sin(x).^2 + cos(x).^2;
matlab_time = toc;
"""
@mget matlab_time
println("MATLAB: ", matlab_time, " 秒数")
Python
py_x = PyObject(x)
py"""
import time
import numpy as np
start = time.perf_counter()
y = np.sin($py_x)**2 + np.cos($py_x)**2
python_time = time.perf_counter() - start
"""
python_time = py"python_time"
println("Python: ", python_time, " 秒数")
Julia演示了由于JIT编译而产生的最低执行时间，该编译生成的机器代码没有解释的开销。
测试2。 大矩阵的乘积
我们将生成两个随机的5000×5000矩阵并执行它们的乘法。
这是最常见的工程测试之一。
A = rand(5000, 5000)
B = rand(5000, 5000)
茱莉亚
julia_time = @belapsed $A * $B
println("Julia: ", julia_time, " 秒数")
Matlab
@mput A B
mat"""
tic;
C = A * B;
matlab_time = toc;
"""
@mget matlab_time
println("Matlab: ", matlab_time, " 秒数")
Python
pyA = PyObject(A)
pyB = PyObject(B)
py"""
tmp = $pyA @ $pyB
"""
py"""
import timeit
python_time = timeit.timeit( lambda: $pyA @ $pyB, number=1)
"""
python_time = py"python_time"
println("Python: ", python_time, " 秒数")
Julia和MATLAB的结果接近，而Python在组织调用的速度上明显逊色。
测试3。 求解线性方程组
让我们解决方程组:
对于矩阵 尺寸为8000×8000。
A = rand(8000, 8000)
b = rand(8000)
茱莉亚
julia_time = @belapsed A\b
println("Julia: ", julia_time, " 秒数")
MATLAB
@mput A b
mat"""
tic;
x = A\\b;
matlab_time = toc;
"""
@mget matlab_time
println("Matlab: ", matlab_time, " 秒数")
Python
pyA = PyObject(A)
pyb = PyObject(b)
py"""
import time
import numpy as np
start = time.perf_counter()
x = np.linalg.solve($pyA, $pyb)
python_time = time.perf_counter() - start
"""
python_time = py"python_time"
println("Python: ", python_time, " 秒数")
Julia和MATLAB的性能相同;Python在矩阵乘法方面略逊一筹。
测试4。 快速傅立叶变换
让我们计算长度为1000万个样本的信号的FFT。
signal = rand(10^7)
茱莉亚
fft(signal)
julia_time = @belapsed fft(signal)
println("Julia: ", julia_time, " 秒数")
Matlab
@mput signal
mat"""
tic;
Y = fft(signal);
matlab_time = toc;
"""
@mget matlab_time
println("Matlab: ", matlab_time, " 秒数")
Python
py_signal = PyObject(signal)
py"""
import time
import numpy as np
start = time.perf_counter()
Y = np.fft.fft($py_signal)
python_time = time.perf_counter() - start
"""
python_time = py"python_time"
println("Python: ", python_time, " 秒数")
在这个测试中，MATLAB显示了最好的结果;Julia略逊一筹，Python显示可以接受，但速度显着降低。
测试5。 数值积分
设置任务
计算积分
在1亿点的网格上使用梯形方法。
N = 10^8
茱莉亚
让我们定义积分函数。
function integrate_julia(N)
h = 1000 / N
s = 0.0
@inbounds for i in 1:N
x = i*h
s += sin(x)*cos(x)*exp(-x/1000)
end
return s*h
end
julia_time = @belapsed integrate_julia(N)
println("Julia: ", julia_time, " 秒数")
MATLAB
@mput N
mat"""
N = double(N);
tic;
h = 1000 / N;
s = 0.0;
for i = 1:N
x = i * h;
s = s + sin(x) * cos(x) * exp(-x / 1000);
end
result = s * h;
matlab_time = toc;
"""
@mget matlab_time
println("MATLAB: ", matlab_time, " 秒数")
Python
py"""
import math
import time
def integrate_python(N):
h = 1000/N
s = 0.0
for i in range(N):
x = i*h
s += math.sin(x)*math.cos(x)*math.exp(-x/1000)
return s*h
start = time.perf_counter()
integrate_python(100000000)
python_time = time.perf_counter() - start
"""
python_time = py"python_time"
println("Python: ", python_time, " 秒数")
Julia几乎是MATLAB的两倍，比Python好一个数量级，这清楚地表明了编译顺序循环的效率。
测试6。 计算数pi的蒙特卡洛方法
让我们在正方形内生成1亿个随机点，并确定它们中有多少属于单位圆内。
这是一个典型的例子，在一个循环中需要很多简单的计算。
茱莉亚
function montecarlo_pi(N)
inside = 0
@inbounds for i in 1:N
x = rand()
y = rand()
if x*x + y*y <= 1.0
inside += 1
end
end
return 4 * inside / N
end
julia_time = @belapsed montecarlo_pi(100_000_000)
println("Julia: ", julia_time, " 秒数")
MATLAB
mat"""
tic
inside = 0;
for i = 1:100000000
x = rand();
y = rand();
if x*x + y*y <= 1
inside = inside + 1;
end
end
pi_est = 4 * inside / 100000000;
matlab_time = toc;
"""
@mget matlab_time
println("MATLAB: ", matlab_time, " 秒数")
Python
py"""
import random
import time
def montecarlo_pi(N):
inside = 0
for _ in range(N):
x = random.random()
y = random.random()
if x*x + y*y <= 1.0:
inside += 1
return 4.0 * inside / N
start = time.perf_counter()
montecarlo_pi(100_000_000)
python_time = time.perf_counter() - start
"""
python_time = py"python_time"
println("Python: ", python_time, " 秒数")
在对随机数生成器进行密集调用的循环中，Julia的运行速度比MATLAB和Python快一个数量级，从而消除了迭代的开销。
测试7。 求解热方程
让我们为热方程执行1000个时间步长
在10万个节点的网格上。
这是典型的工程任务之一。
u = rand(10^5)
steps = 10^3
茱莉亚
function heat1d(u, α, steps)
tmp = similar(u)
for n in 1:steps
@inbounds for i in 2:length(u)-1
tmp[i] =
u[i] +
α*(u[i+1] - 2u[i] + u[i-1])
end
u, tmp = tmp, u
end
return u
end
julia_time = @belapsed heat1d(copy($u), 0.1, steps)
println("Julia: ", julia_time, " 秒数")
Matlab
@mput u
@mput steps
mat"""
u = double(u);
tmp = zeros(size(u));
tic
for n = 1:steps
for i = 2:length(u)-1
tmp(i) = u(i) + 0.1*(u(i+1) - 2*u(i) + u(i-1));
end
t = u;
u = tmp;
tmp = t;
end
matlab_time = toc;
"""
@mget matlab_time
println("MATLAB: ", matlab_time, " 秒数")
Python
pyu = PyObject(u)
psteps = PyObject(steps)
py"""
import numpy as np
import time
def heat1d(u, alpha, steps):
tmp = np.empty_like(u)
for _ in range(steps):
for i in range(1, len(u)-1):
tmp[i] = (
u[i]
+ alpha*(u[i+1] - 2*u[i] + u[i-1])
)
u, tmp = tmp, u
return u
start = time.perf_counter()
heat1d($pyu, 0.1, $psteps)
python_time = time.perf_counter() - start
"""
python_time = py"python_time"
println("Python: ", python_time, " 秒数")
Julia的性能优于MATLAB几乎一个数量级，Python超过三个数量级，揭示了嵌套循环编译的根本优势。
结果
让我们与计算时间的测量结果做一个比较表。
|
测试 |
朱莉娅时间，c |
Matlab时间，c |
Python时间，c |
|
矢量计算 |
0.2139 |
0.2791 |
0.3749 |
|
大矩阵的乘积 |
2.1764 |
2.3193 |
7.0526 |
|
求解线性方程组 |
5.7585 |
6.0654 |
11.1721 |
|
快速傅立叶变换 |
1.0116 |
0.4867 |
2.5993 |
|
数值积分 |
2.7059 |
5.5843 |
25.8805 |
|
计算数π的蒙特卡洛方法 |
0.3290 |
9.2279 |
15.3714 |
|
求解热方程 |
0.0715 |
0.6327 |
76.2126 |
根据该组测试，Julia演示了最佳的加权平均性能，在某些操作上显示出与MATLAB的奇偶校验以及在具有密集周期的任务中的多重优越性。 尽管Python具有高度优化的库，但在需要非平凡逻辑的地方，它在系统上是逊色的。
结论
测试证明了Julia在工程计算方面的高效率。 在归结为调用优化库（矩阵乘法，SLOUGH求解）的任务中，Julia处于MATLAB级别，领先于Python。 关键优势体现在无法简化为纯矢量化的场景中，例如数值积分，蒙特卡洛方法和进化方案，其中与MATLAB的差距达到一个数量级，与Python的差距—两个或三个数量级。
从工程的角度来看，Julia在AnyMath中特别有吸引力：与专有MATLAB不同，Julia是一种开放技术，与AnyMath结合使用，提供计算的可重复性和缩放的灵活性;此外，python的基本限制被 因此，对于广泛的计算密集型工程任务，Julia提供了性能，成本效益和科学再现性的最佳平衡。