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.
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).
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 deallocationob_type, which encodes the type of the variableob_size, which specifies the size of the following data membersob_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:
Python:
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):
Built-in (entire output):
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
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)
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).
