Performance comparison of Julia, MATLAB and Python

Introduction

When choosing a tool for engineering calculations, not only convenience and the number of libraries are important, but also the speed of calculations.

MATLAB is considered the standard for scientific computing, Python has become widespread thanks to the NumPy and SciPy ecosystem, and Julia was created specifically for high-performance computing.

In this example, we will compare the performance of Julia, MATLAB, and Python on several typical engineering tasks.:

Vector calculations;
Matrix operations;
solving a system of linear equations;
Fast Fourier Transform (FFT) calculation;
numerical integration;
The Monte Carlo method for finding the number π;
Solving the equation of thermal conductivity.

In all the tests, we will measure the execution time and compare the results.

For the sake of objectivity, we will measure time using the built-in tools.:

@belapsedin Julia;

tic … toc in MATLAB;

time.perf_counter()in Python.

Connecting libraries

To run the tests, we will need the following libraries:

MATLAB — calling MATLAB from Julia;
PyCall — interaction with Python;
LinearAlgebra— for matrix operations;
FFTW— for fast Fourier transform;
BenchmarkTools — Accurate time measurement.

using MATLAB, PyCall, BenchmarkTools, LinearAlgebra, FFTW

Preparation of Python modules

We will connect the necessary Python libraries.

np = pyimport("numpy")
scipy_integrate = pyimport("scipy.integrate")

PyObject <module 'scipy.integrate' from '/opt/python3.11/lib/python3.11/site-packages/scipy/integrate/__init__.py'>

Test 1. Vector calculations

Calculate the expression

for an array of 10 million elements.

This expression is specially chosen as an example of intensive piecemeal calculations.

N = 10^7
x = rand(N)

10000000-element Vector{Float64}:
 0.3209199928094635
 0.1865690245097057
 0.027519003346574533
 0.1035696000771914
 0.7617551461014622
 0.3285472660096638
 0.8569073684267615
 0.4177394290269377
 0.403488227552495
 0.6147429411864976
 0.8956454977439492
 0.9337120189173626
 0.9543570461710522
 ⋮
 0.6089946931248793
 0.9624634504691789
 0.6232579150388496
 0.49009958177103763
 0.774540432270273
 0.7570462764844246
 0.8760435688206902
 0.3420754138897466
 0.6462666791429025
 0.6322249655087943
 0.21683558167577033
 0.49816246221767047

Julia

julia_time = @belapsed sin.($x).^2 .+ cos.($x).^2
println("Julia: ", julia_time, " seconds")

Julia: 0.213971722 seconds

MATLAB

Let's transfer the data to MATLAB and measure the execution time.

@mput x

mat"""
tic;
y = sin(x).^2 + cos(x).^2;
matlab_time = toc;
"""
@mget matlab_time
println("MATLAB: ", matlab_time, " seconds")

MATLAB: 0.279133

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, " seconds")

Python: 0.3749271952547133 seconds

Julia demonstrates the lowest execution time due to JIT compilation, which generates machine code without the overhead of interpretation.

Test 2. The product of large matrices

We will generate two random 5000×5000 matrices and perform their multiplication.

This is one of the most common engineering tests.

A = rand(5000, 5000)
B = rand(5000, 5000)

5000×5000 Matrix{Float64}:
 0.247848   0.275075   0.103041   …  0.955899   0.442469  0.142161
 0.829696   0.570888   0.779525      0.803322   0.261279  0.416469
 0.302144   0.470861   0.929515      0.879647   0.463827  0.659836
 0.622374   0.715513   0.749483      0.790127   0.532119  0.0226104
 0.22009    0.807178   0.496316      0.0163167  0.143018  0.68864
 0.0202647  0.932929   0.196466   …  0.137595   0.750025  0.830217
 0.737694   0.466092   0.52982       0.103175   0.552223  0.872671
 0.182882   0.327355   0.803583      0.274724   0.901884  0.520132
 0.684589   0.0492455  0.0692971     0.276476   0.475649  0.246488
 0.161581   0.674723   0.786201      0.141705   0.594002  0.243299
 0.834808   0.484064   0.874741   …  0.366673   0.909212  0.195076
 0.78583    0.849583   0.549752      0.937795   0.488873  0.396699
 0.482646   0.689754   0.191593      0.287839   0.90025   0.0900362
 ⋮                                ⋱                       
 0.0490212  0.362037   0.0800845     0.814074   0.688598  0.0435109
 0.646176   0.398827   0.342848      0.520501   0.293243  0.948626
 0.0814196  0.786904   0.694019   …  0.717057   0.894803  0.214561
 0.205205   0.599291   0.604951      0.284944   0.892169  0.670361
 0.746486   0.425566   0.838346      0.281594   0.641695  0.319097
 0.160965   0.64698    0.371449      0.730004   0.395501  0.110665
 0.64245    0.0561828  0.924955      0.713472   0.702081  0.368934
 0.331956   0.273443   0.749144   …  0.622444   0.8223    0.736889
 0.494684   0.795712   0.231628      0.986563   0.707547  0.639998
 0.0958006  0.190157   0.610949      0.561864   0.527178  0.922572
 0.188368   0.56108    0.954863      0.0105769  0.689971  0.779974
 0.277134   0.151909   0.227931      0.534004   0.353346  0.575454

Julia

julia_time = @belapsed $A * $B
println("Julia: ", julia_time, " seconds")

Julia: 2.176358714 seconds

Matlab

@mput A B

mat"""
tic;
C = A * B;
matlab_time = toc;
"""
@mget matlab_time

println("Matlab: ", matlab_time, " seconds")

Matlab: 2.319368 seconds

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, " seconds")

Python: 7.0526929018087685 seconds

The results of Julia and MATLAB are close, while Python is noticeably inferior in the speed of organizing calls.

Test 3. Solving a system of linear equations

Let's solve the system of equations:

for the matrix the size is 8000×8000.

A = rand(8000, 8000)
b = rand(8000)

Julia

julia_time = @belapsed A\b
println("Julia: ", julia_time, " seconds")

Julia: 5.758521299 seconds

MATLAB

@mput A b

mat"""
tic;
x = A\\b;
matlab_time = toc;
"""
@mget matlab_time
println("Matlab: ", matlab_time, " seconds")

Matlab: 6.065439 seconds

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, " seconds")

Python: 11.172098383773118 seconds

The performance of Julia and MATLAB is identical; Python is slightly inferior as well as in matrix multiplication.

Test 4. Fast Fourier Transform

Let's calculate the FFT of a signal with a length of 10 million samples.

signal = rand(10^7)

10000000-element Vector{Float64}:
 0.19964052302019852
 0.5334663195022673
 0.9121321253899346
 0.04129394526409169
 0.6589395955777433
 0.5388352658418127
 0.26895775008531053
 0.5829179391433303
 0.6515824751637956
 0.836612931381444
 0.6559235223843917
 0.6026254909782042
 0.5890734109652789
 ⋮
 0.6905007831051632
 0.2881370605087836
 0.5270581989024726
 0.9312542593606681
 0.5850647418223617
 0.11650802263000237
 0.9507322240086514
 0.10666018251459808
 0.21469656971547735
 0.4365831131657778
 0.5078941848918891
 0.657329997060914

Julia

fft(signal)
julia_time = @belapsed fft(signal)
println("Julia: ", julia_time, " seconds")

Julia: 1.01158086 seconds

Matlab

@mput signal
mat"""
tic;
Y = fft(signal);
matlab_time = toc;
"""
@mget matlab_time
println("Matlab: ", matlab_time, " seconds")

Matlab: 0.486747 seconds

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, " seconds")

Python: 2.5993731832131743 seconds

In this test, MATLAB shows the best result; Julia is slightly inferior, Python demonstrates acceptable, but significantly lower speed.

Test 5. Numerical integration

Setting the task

Calculate the integral

using the trapezoid method on a grid of 100 million points.

N = 10^8

100000000

Julia

Let's define the integration function.

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

integrate_julia (generic function with 1 method)

julia_time = @belapsed integrate_julia(N)
println("Julia: ", julia_time, " seconds")

Julia: 2.705959406 seconds

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, " seconds")

MATLAB: 5.584325 seconds

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, " seconds")

Python: 25.88057485362515 seconds

Julia is almost twice as good as MATLAB and an order of magnitude better than Python, which clearly demonstrates the efficiency of compiling sequential loops.

Test 6. Monte Carlo method for calculating the number pi

Let's generate 100 million random points inside the square and determine how many of them fall within the unit circle.

This is a classic example where a lot of simple calculations are required in a loop.

Julia

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, " seconds")

Julia: 0.328965659 seconds

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, " seconds")

MATLAB: 9.227858 seconds

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, " seconds")

Python: 15.37143736006692 seconds

On a loop with intensive calls to the random number generator, Julia runs an order of magnitude faster than MATLAB and Python, eliminating the overhead of iteration.

Test 7. Solving the heat equation

Let's perform 1000 time steps for the heat equation

on a grid of 100 thousand nodes.

This is one of the typical engineering tasks.

u = rand(10^5)
steps = 10^3

1000

Julia

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, " seconds")

Julia: 0.071473353 seconds

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, " seconds")

MATLAB: 0.632743 seconds

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, " seconds")

Python: 76.2125775879249 seconds

Julia outperforms MATLAB by almost an order of magnitude, and Python by more than three orders of magnitude, revealing the fundamental advantage of nested loop compilation.

Results

Let's make a comparative table with the measurement results of the calculation time.

Test	Julia time, c	Matlab time, c	Python time, c
Vector calculations	0.2139	0.2791	0.3749
The product of large matrices	2.1764	2.3193	7.0526
Solving a system of linear equations	5.7585	6.0654	11.1721
Fast Fourier transform	1.0116	0.4867	2.5993
Numerical integration	2.7059	5.5843	25.8805
The Monte Carlo method for calculating the number π	0.3290	9.2279	15.3714
Solving the heat equation	0.0715	0.6327	76.2126

According to the set of tests, Julia demonstrates the best weighted average performance, showing parity with MATLAB on some operations and multiple superiority in tasks with intensive cycles. Python, despite its highly optimized libraries, is systematically inferior where non-trivial logic is required.

Conclusion

Testing has demonstrated Julia's high efficiency for engineering computing. In tasks that come down to calling optimized libraries (matrix multiplication, SLOUGH solving), Julia is at the MATLAB level, ahead of Python. The key advantage is revealed in scenarios that cannot be reduced to pure vectorization, such as numerical integration, the Monte Carlo method, and evolutionary schemes, where the gap with MATLAB reaches an order of magnitude, and with Python — two or three orders of magnitude.

From an engineering point of view, Julia is particularly attractive in Engee: unlike proprietary MATLAB, Julia is an open technology that, in combination with Engee, provides reproducibility of calculations and flexibility of scaling; in addition, Python's fundamental limitation is overcome — the low speed of cyclic code execution, which is critical for simulation. Thus, for a wide range of computationally intensive engineering tasks, Julia provides an optimal balance of performance, cost-effectiveness, and scientific reproducibility.