6
$\begingroup$

When compared to Python the main selling point of Julia is its speed as it is often argued. However, from my own personal experience I never noticed any significant difference in speed between Julia and Python. If there is a trivial few seconds difference that would not matter practically. This could be problem for models using big-data but outside machine learning are there any (even large) DSGE models where this makes difference?

I am looking for some examples of significant differences between the two.

PS: Despite reference request tag I am happy to accept answer that showcases some examples directly.

$\endgroup$
2
  • 1
    $\begingroup$ Random effects can heavily influence any measurements in elapsed time if the typical task takes a few seconds. You would have to run a large task, say one that takes 2 minutes or 2 hours or 2 days. Then implement the same in the other language in order to compare. If your task runs for 2 days consider the support for parallel processing. The language that makes parallel processing easier might then be preferred. If it takes 2 hours, consider how often you need to run it. If it is seldom then it might not matter if the ratio of elapsed times between the two languages is say 5 to 1. $\endgroup$
    – H2ONaCl
    Nov 18, 2021 at 4:59
  • 2
    $\begingroup$ “I never noticed any significant difference in speed between Julia and Python” – you're probably talking about NumPy's (or even e.g. PyTorch's) speed, rather than Python's. IOW, the speed of C / C++, not of Python. –Now, in Julia something similar could be probably said too... but at least, you can actually write general-purpose code in Julia that runs at reasonable speed. In Python, this is only possible with e.g. Numba, which is a very different experience from standard Python. $\endgroup$ Nov 18, 2021 at 17:42

3 Answers 3

5
$\begingroup$

Julia is actually a lot faster than Python, also when running DSGE models.

The NY FED moved their DSGE model to Julia because it allows them to:

  • Estimate models 10x faster
  • Complete 'solve' test 11x faster
  • Reduce number of lines of code in half, saving time, increasing readability and reducing errors

The original results can be found here.

These speedups are compared to Matlab which the FED used before. However, with regards to speed, Matlab should be quicker than Python for DSGE models (any matrix statistics, iteration and recursion and so forth). That is also what the first suggested paper in the accepted answer finds (I only looked at the revised 2018 version). By the way, the conclusion in the paper reads as follows:

... C++ is the fastest alternative, Julia offers great balance of speed and ease of use, and Python is too slow.

The screenshot below is directly from the paper. enter image description here

The slowest Julia implementation takes only 1.6% of the fastest Python version.

The second paper (2020) essentially tests a model that runs only a few seconds (less than 1 in Julia). enter image description here

Nonetheless, as soon as it's more than a fraction of a second (where compile time latency of Julia will almost certainly be an issue), Julia is significantly quicker. It may not be a difference you really feel or that matters for you, because waiting for 10 or 15 seconds is essentially not a deal breaker.

Unlike Python, Julia is significantly more complex when it comes to potential code implementation. Everything that runs fast in Python is not Python, but C, Fortran and co. Julia can also turn slow when type inference fails or the JIT compiler has insufficient information to optimize effectively. If you notice no speedup to Python, it is either because you use examples (or benchmark that way) where compile time latency matters more than computation, or it is inefficient code. I did some work a while ago on Julia's speed and comparison to Python. I will copy paste some sections in case anyone is interested in some details why Julia is actually fast.

In dynamic languages like Python, classes could be subclassed. This is also possible in some statically typed languages that are object oriented such as Java. For this reason, a single integer in Python 3.x actually contains four pieces:

  • ob_refcnt, a reference count that helps Python silently handle memory allocation and deallocation
  • ob_type, which encodes the type of the variable
  • ob_size, which specifies the size of the following data members
  • ob_digit, which contains the actual integer value that we expect the Python variable to represent.

This means that there is some overhead in storing an integer in Python as compared to an integer in say Julia or C. Python uses 28 bytes to store an integer, Julia only 8 bytes.

A great feature of Julia is that types are extremely lightweight, and user-defined types are exactly as performant as anything built-in: this is a big difference to languages like Python where constructing a new instance of a class is very expensive. The amount of data required to store a Struct (Class) called Foo is the same as the amount of data required to store a pointer to a string. At compilation, Julia (just like Fortran and C) doesn't have additional overhead contrary to how an object (in say Python) is represented in memory (includes type information, reference counter etc).

Julia:

enter image description here

Python:

enter image description here

You can see the massive overhead in Python. It's probably even better to say that Julia isn't slow (not that it is fast), meaning it doesn't get in the way of the hardware and lets the hardware do its job efficiently. Mathematically, the integer 5 is exactly identical to the floating point 5.0, but it has a different bitwise representation. This is why Types are important in computers. Slow programming languages have a very hard time to know in advance what types things are. So they test for all sorts of combinations to figure out what routine to use to call the correct CPU instructions, because at the lowest level, the microprocessor has to perform different instructions when dealing with for example integers or floats. A microprocessor does not know types. A CPU needs special instructions for all different types; e.g. signed and unsigned integers, floats with single precision or double precision, or combinations of floating points and integers, for which you will need conversion (promotion) of types. Usually there will be different special instructions for every possible type.

As long as your task is computational, Julia will always beat Python (provided both are optimized). Where it starts to matter is when you really work with lots of data and / or complex models.

For example look at the question Why is this task faster in Python than Julia? on stackoverflow. It is almost always code written in non-performant ways if you get slow Julia code. Without any optimization, just getting rid of compile time latency, resulted in Julia taking 8s vs 30s in Python. Regarding the stackoverflow question, Quantinsti has some good tests for Multiple operations on large datasets. For 100 groups of ~10,000,000 rows, Python (pandas package) and R (dplyr package) resulted in an internal error and out of memory error, respectively while Julia took 2.4 seconds the first time and 1.8 seconds the second time.

Julia is fast because of its design decisions that allow a compiler to make efficient code. This is achieved with type-stability through specialization via multiple-dispatch as a core design decision. To illustrate this, you can try to compute 2^3.

In Julia, you can run @code_llvm as well as @code_native macros to print the LLVM bitcodes as well as the native assembly instructions generated for running the method matching the given generic function and type.

If you use @code_llvm exponential(2,3) as well as @code_llvm exponential(2,3) and compared it to the 2^3 counterpart, you will be surprised how messy the “manual” exponential solution is compared to Julia’s built in.

Manual (ignoring almost half the output):

enter image description here

Built-in (entire output):

enter image description here

The built in solution is something you can decipher fairly quickly with some basic assembly.

  • The lines beginning with ; are comments, and explain which section of the code the following instructions come from. Calling it with debuginfo=:none removes the comments in the assembly code
  • pushq , movq, subq and addq are to save stack frames, move content to or from registers (the suffix q is short for quad, which means 64-bit integer), allocate to registers and dealocate. Neither is important for our purpose. For short, a CPU operates only on data present in registers. These are small, fixed size slots (e.g. 8 bytes in size) inside the CPU itself. A register is meant to hold one single piece of data, like an integer or a floating point number.
  • movabsq is the so called Exponentiation by squaring. For example if you compute $x^{15}$, it will look like this
x^15 = (x^7)*(x^7)*x 
x^7 = (x^3)*(x^3)*x 
x^3 = x*x*x

Or only the last line in our simple example above, which is $x^3$ only.

The reason this works is that Julia's built in function made some clever design decisions. For instance, x^-n works for an integer x and literal n but not for expressions because Julia developers defined a special syntactic treatment of literals.

  • The benefit is quick and reliable exponentiation for something like $x^-3$. The advantage is that literal powers automatically become type-stable, because they turn into a fixed number of * calls. This is why Julia now lowers x^-3 to inv(x)^3.
  • The downside is that $p = -3$ and $x^p$ does not work and throws an error which seems potentially quite confusing. However, the x^literal has a different meaning than x^expression. In essence, referential transparency was sacrificed, and type stability "extended": That is why ^ to a literal integer power is different than raising to a variable with the same integer value.

Note: The crux here is that type-stability doesn’t mean that the type returned by a function is the same as the input. It means that the type can be inferred at every step of the way. If inv of an integer is always a floating point number, then that’s type stable.

The problem however, is that most casual Python users will expect x^p to work for p = -3.

Similarly, scope matters a lot. You can look at this simple function.

b = 1.0
function g(a)
    global b
    tmp = a
    for i in 1:1_000_000
        tmp = tmp + a + b
    end
    return tmp
end

b is global here, which is like poison for performance because the type of a global variable is never certain.

If you instead write

function g(a, b)
    tmp = a
    for i in 1:1_000_000
        tmp = tmp + a + b
    end
    return tmp
end

you are eliminating the global variable which

• reduces the number of allocations to zero (from 3000001 allocations: 45.78 MiB)
• speeds up execution and
• produces clean and fairly simple machine code.

For short, with some care, Julia can be on par with the fastest languages.

These allocations matter massively. On my private laptop for example, finding a single file (including opening and closing the file) takes about 1.4 milliseconds. Accessing 1,000,000 integers from memory takes 75.7 milliseconds. So my RAM is almost 20,000 times faster than my disk. I ran this with a reasonably fast SSD but even the newer Optane technology disks are usually thousands of times slower than RAM.

I ran a speed test between C, Python and Julia, which is based on a notebook from an introductory MIT math course.

It is a speed comparisons between C, Python and Julia by implementing a sum function sum(a), which computes

$$ \mathrm{sum}(a) = \sum_{i=1}^n a_i $$

for an array a with n elements. I compared built in sum functions in Julia and Python along with hand-coded implementations in C, Python, and Julia. If you are a Windows user and want to replicate this on your own, you will need to install the MinGW (GCC/G++) compiler.

The C code is:

C_code = """
#include <stddef.h>
double c_sum(size_t n, double *X) {
    double s = 0.0;
    for (size_t i = 0; i < n; ++i) {
        s += X[i];
    }
    return s;
}
"""

const Clib = tempname()   # make a temporary file

# compile to a shared library by piping C_code to gcc
open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code) 
end

Julia's

function my_sum1(x)
    result = zero(eltype(x))
    for element in x
        result += element
    end
    return result
end

Not only is the Julia code significantly easier, but also, unlike C, it works in Julia for any type of array (any iterable that provides an eltype method).

Python (within a Julia Kernel using PyCall)

py"""
def mysum(a):
    s = 0.0
    for x in a:
        s = s + x
    return s
"""
mysum_py = py"mysum"

The results are

enter image description here

There are a few interesting aspects.

  • Starting with C, testing the function on an array of $10^7$ random numbers in [0,1] gives a result of ≈ 10ms for summing these numbers, which means that we have about 1 billion additions per second. This sounds a lot, but is well below the clock-speed of 2.5GHZ my processor has. This is because each Floating-Point addition needs to perform several additional calculations to load the next element of the array from memory, plus the time needed to access memory itself.
  • Probably most interesting is that Python's built in sum is fairly slow, even though the Python sum function is written in C (Julia's sum is built using only Julia). It takes almost 4x longer than the equivalent C code and allocates memory. Python code pays a price for being generic and being able to handle arbitrary iterable data structures. See above for memory allocation differences between Julia and Python. Therefore, there are not only the computations involving addition, but also the overhead from fetching each item from memory.
  • Julia equivalent of a Python list is a Vector{Any}: internally, this is an array of pointers to "boxes" that can hold any type (Any). This makes things much slower: each + computation on an Any value must dynamically look up the type of object, figure out what + function to call, and allocate a new "box" to store the result. If you run this logic in Julia, it is even slower than Python. This is the same problem as with the manual Numpy array implementation. Python is faster (better optimized for dealing) with untyped Any values. This is because Julia expects you to use "concretely" typed arrays in all performance-critical cases and why you should always avoid abstract types in Julia.
  • On the other hand, Numpy arrays can take advantage of the fact that all of the elements are of the same type. This is actually faster than C! The reason is that NumPy gets and extra turbo boost from exploiting SIMD instructions.

You can run SIMD for C as well as Julia as well.

Julia

function mysum(A)
    s = zero(eltype(A)) 
    @simd for a in A
        s += a
    end
    return s
end

C

const Clib_fastmath = tempname()   # make a temporary file

# The same as above but with a -ffast-math flag added
open(`gcc -fPIC -O3 -msse3 -xc -shared -ffast-math -o $(Clib_fastmath * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code) 
end

# define a Julia function that calls the C function:
c_sum_fastmath(X::Array{Float64}) = ccall(("c_sum", Clib_fastmath), Float64, (Csize_t, Ptr{Float64}), length(X), X)

enter image description here

However, once you work with SIMD, you really need to know what you do when you implement and use this. Incorrect use of the @simd macro may cause unexpected results.

While this may sound like esoteric computer science stuff, you can look at the following question to see what can happen with floating point math results in a reasonable normal situation. The most intriguing example demonstrating the dangers of floating point math I have seen so far comes from Stefan Karpinksi. You can use this to illustrate the difference between left-to-right summation and left-to-right summation with a @simd annotation.

Essentially, if you use this code snippet in Julia,

function sumsto(x::Float64)
    0 <= x < exp2(970) || throw(ArgumentError("sum must be in [0,2^970)"))
    n, p₀ = Base.decompose(x) # integers such that `n*exp2(p₀) == x`
    [floatmax(); [exp2(p) for p in reverse(-1074:969) if iseven(n >> (p-p₀))]

    -floatmax(); [exp2(p) for p in reverse(-1074:969) if isodd(n >> (p-p₀))]]
end    

s = rand()*10^6
v = sumsto(s)

you order the random vector (10^6) from highest to lowest value. None of the standard algorithms in Julia or Python provide the true sum, although interestingly, in this example SIMD comes closest in this case. The reason is that in the case of summation, strict left-to-right summation is actually one of the least accurate algorithms, so putting @simd on the loop makes it faster and more accurate at computing the true sum because SIMD largely allows for reassociation of mathematically associative operations which are not actually associative in floating point. This includes + and *. Generally, (a + b) + c can be evaluated as a + (b + c) instead, but this does not compute the same thing for floating-point values.

The correct answer can be obtained with the Kahan-Babuska-Neumaier (KBN) summation algorithm for computing sums (and cumulative sums). This comes at a speed difference of about 20x, which is why none of the conventional implementations uses it.

Bottom line is, in my opinion, that Julia is serving a niche market. While the paper states that Julia comes with a great mix of speed and usability, the main problem is that it is still quite a task to reach C (C++ or Fortran speed) in Julia. If you require real speed, firms will prefer C++ (at least in my field, quant finance, C++ is the by far most frequently used language, with some niche products like OCAML used by Lexifi or Bloomberg's DLIB, which focus on complex derivs). Most users don't need that speed because the models are to small, or the data is not big enough. E.g. an everyday Bloomberg user could never download as much data (BBG has a limit on daily data usage) that it really makes a noticeable difference, just like ~150 seconds (compared to 2 in Julia) don't make a huge difference if you run a DSGE model because there is no need to rerun it every 3 minutes or at even shorter periods.

That said, packages are really developing in Julia, and charting functionalities and the DataFrames package are really well done. You can look here to see how little code is needed to make some kinda decent interactive charts (dynamic option PnL charts).

Also, since the release of 1.0, you don't usually face issues with breaking code due to deprecation anymore. If everyday datasets grow, and computational power becomes more important for casual programmers, Julia could very well become a big player (subjective opinion).

$\endgroup$
5
$\begingroup$

There are certainly cases where the performance differences between Python and Julia matter, but solving a simple DSGE model is not one of them.

There have been some formal comparison exercises (e.g., Aruoba and Fernández-Villaverde 2015 and 2018 update). Coleman et al. 2020 also consider a NK model (see Table 4). The real gains are only found when considering higher-order polynomials/perturbations.

In practice, I would say that it is mostly with Bayesian DSGE models and HANK models that you could realize improvements.

$\endgroup$
4
  • 1
    $\begingroup$ A common criticism with some of these language comparison papers is that the authors don't use the most efficient code for each language. E.g., in some cases the authors could have avoided slowdowns associated with for-loops by vectorizing. Just something to watch out for. Regarding Python v Julia case, notice in the first paper that Python with Numba is faster than Julia. It's strange to consider Numba as a hybrid language or special case when it usually only requires adding a line or two more of code. $\endgroup$
    – jmbejara
    Nov 18, 2021 at 17:40
  • 4
    $\begingroup$ @jmbejara “usually only requires adding a line or two more” is a gross overstatement. Yes, sometimes that's possible – but only for algorithms that are very array-focused and monomorphic to begin with. IOW the kind of thing that could probably have been done with pre-existing libraries as well. But for actual nontrivial code, Numba can be a huge pain. $\endgroup$ Nov 18, 2021 at 17:46
  • $\begingroup$ @leftaroundabout Fair enough! $\endgroup$
    – jmbejara
    Nov 18, 2021 at 17:56
  • $\begingroup$ @jmbejara +1 That's a valid point, those comparison papers are usually flawed in that sense... $\endgroup$
    – Joe
    Nov 19, 2021 at 9:37
2
$\begingroup$

The comparison between julia and python is apples to oranges. They have completely different architectures and performance characteristics. Julia uses JIT compilation to native machine code, which means that there is a big performance hit at startup, but after that it's quite efficient. Python is an interpreted language that was never designed for speed, and its performance is extremely poor compared to other interpreted languages such as perl and javascript. However, python's startup time is negligible for most purposes.

None of this actually matters much in many of the numerical analysis tasks that people use julia and numpy for, because operations like matrix inversion or Fourier analysis are being done by library routines that are not written native code anyway.

In my experience, the big difference between python and julia is the maturity and quality of their implementations. When I tried julia in 2015 it was extremely buggy, to the point of not being usable in my experience. It may have improved since then.

The huge win with julia is if you already have some matlab code, and you don't want to pay for matlab or have a political affinity for open source. Converting matlab to julia is very straightforward. Julia also has a very friendly user community and excellent documentation.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.