Is Julia fast?¶

Julia isn't fast per se.

One can write terribly slow code in any language, including Julia.

So let's ask a different question.

Can Julia be fast?¶

Microbenchmarks¶

More realistic case: Vandermonde matrix¶

Vandermonde matrix: \begin{align}V=\begin{bmatrix}1&\alpha _{1}&\alpha _{1}^{2}&\dots &\alpha _{1}^{n-1}\\1&\alpha _{2}&\alpha _{2}^{2}&\dots &\alpha _{2}^{n-1}\\1&\alpha _{3}&\alpha _{3}^{2}&\dots &\alpha _{3}^{n-1}\\\vdots &\vdots &\vdots &\ddots &\vdots \\1&\alpha _{m}&\alpha _{m}^{2}&\dots &\alpha _{m}^{n-1}\end{bmatrix}\end{align}

In [ ]:

using PyCall

In [ ]:

np = pyimport("numpy")

In [ ]:

np.vander(1:5, increasing=true)

The source code for this function is here. It calls np.multiply.accumulate which is implemented in C here. However, this code doesn't actually perform the computation, it basically only checks types and stuff. The actual kernel that gets called is here. This isn't even C code but a template for C code which is used to generate type specific kernels.

Overall, this setup only supports a limited set of types, like Float64, Float32, and so forth.

Here is a simple Julia implementation

In [ ]:

function vander(x::AbstractVector{T},  n=length(x)) where T
    m = length(x)
    V = Matrix{T}(undef, m, n)
    for j = 1:m
        V[j,1] = one(x[j])
    end
    for i= 2:n
        for j = 1:m
            V[j,i] = x[j] * V[j,i-1]
            end
        end
    return V
end

In [ ]:

vander(1:5)

A quick speed comparison¶

Show Code

using BenchmarkTools, Plots
ns = exp10.(range(1, 4, length=30));

tnp = Float64[]
tjl = Float64[]
for n in ns
    x = 1:n |> collect
    push!(tnp, @belapsed np.vander(\$x) samples=3 evals=1)
    push!(tjl, @belapsed vander(\$x) samples=3 evals=1)
end
plot(ns, tnp./tjl, m=:circle, xscale=:log10, xlab="matrix size", ylab="NumPy time / Julia time", legend=:false)

Note that the clean and concise Julia implementation is beating numpy's C implementation for small matrices and is on-par for large matrix sizes.

At the same time, the Julia code is generic and works for arbitrary types!

In [ ]:

vander(Int32[4, 8, 16, 32])

It even works for non-numerical types. The only requirement is that the type has a one (identity element) and a multiplication operation defined.

In [ ]:

vander(["this", "is", "a", "test"])

Here, one(String) == "" since the empty string is the identity under multiplication (string concatenation).

How can Julia be fast?¶

drawing

AST = Abstract Syntax Tree

SSA = Static Single Assignment

LLVM = Low Level Virtual Machine

Specialization and code inspection¶

Julia specializes on the types of function arguments.

When a function is called for the first time, Julia compiles efficient machine code for the given input types.

If it is called again, the already existing machine code is reused, until we call the function with different input types.

In [ ]:

func(x,y) = x^2 + y

In [ ]:

@time func(1,2)
@time func(1,2)

First call: compilation + running the code

Second call: running the code

In [ ]:

@time func(1,2)

If the input types change, Julia compiles a new specialization of the function.

In [ ]:

@time func(1.3,4.8)
@time func(1.3,4.8)

We now have two efficient codes, one for all Int64 inputs and another one for all Float64 arguments, in the cache.

But I really want to see what happens!¶

We can inspect the code at all transformation stages with a bunch of macros:

The AST after parsing (@macroexpand)
The AST after lowering (@code_typed, @code_warntype)
The AST after type inference and optimization (@code_lowered)
The LLVM IR (@code_llvm)
The assembly machine code (@code_native)

In [ ]:

@code_typed func(1,2)

In [ ]:

@code_lowered func(1,2)

In [ ]:

@code_llvm func(1,2)

We can remove the comments (lines starting with ; using debuginfo=:none).

In [ ]:

@code_llvm debuginfo=:none func(1,2)

In [ ]:

@code_native func(1,2)

Let's compare this to Float64 input.

In [ ]:

@code_native func(1.2,2.9)

How important is code specialization?¶

Let's try to estimate the performance gain by specialization.

We wrap our numbers into a custom type which internally stores them as Any to prevent specialization.

(This is qualitatively comparable to what Python does.)

In [ ]:

struct Anything
    value::Any
end

operation(x::Number) = x^2 + sqrt(x)
operation(x::Anything) = x.value^2 + sqrt(x.value)

In [ ]:

using BenchmarkTools

@btime operation(2);
@btime operation(2.0);

x = Anything(2.0)
@btime operation($x);

That's about an 40 times slowdown!

In [ ]:

@code_native operation(2.0)

In [ ]:

@code_native operation(x)

Make run-time the fun time.¶

In scientific computations, we typically run a piece of code many times over and over again. Think of a Monte Carlo simulation, for example, where we perform the update and the Metropolis check millions of times.

Therefore, we want our run-time to be as short as possible.

On the other hand, for a given set of input arguments, Julia compiles the piece of code only once, as we have seen above. The time it takes to compile our code is almost always negligible compared to the duration of the full computation.

A general strategy is therefore to move parts of the computation to compile-time.

Since Julia specializes on types, at compile-time only type information is available to the compiler.

In [ ]:

f1(x::Int) = x + 1
f2(x::Int) = x + 2

function f_slow(x::Int, p::Bool)
    if p                                # check depends on the value of p
        return f1(x)
    else
        return f2(x)
    end
end

In [ ]:

@code_llvm debuginfo=:none f_slow(1, true)

We can eliminate the if branch by moving the condition check to the type domain. This way, it will only be evaluated once at compile-time.

In [ ]:

abstract type Boolean end
struct True <: Boolean end # type domain true
struct False <: Boolean end # type domain false

function f_fast(x::Int, p::Boolean)
    if typeof(p) == True                # check solely based on the type of p
        return f1(x)
    else
        return f2(x)
    end
end

In [ ]:

@code_llvm debuginfo=:none f_fast(1, True())

Are explicit type annotations necessary? (like in C or Fortran)¶

Note that Julia's type inference is powerful. Specifying types is not necessary for best performance!

In [ ]:

function my_function(x)
    y = rand()
    z = rand()
    x+y+z
end

function my_function_typed(x::Int)::Float64
    y::Float64 = rand()
    z::Float64 = rand()
    x+y+z
end

In [ ]:

@btime my_function(10);
@btime my_function_typed(10);

However, annotating types explicitly can serve a purpose.

Define a user interface/type filter (will throw error if incompatible type is given)
Enforce conversions
Rarely, help the compiler infer types in tricky situations

Core messages of this Notebook¶

Julia can be fast.
A function is compiled when called for the first time with a given set of argument types.
The are multiple compilation steps all of which can be inspected through macros like @code_warntype.
Code specialization based on the types of all of the input arguments is important for speed.
Calculations can be moved to compile-time to make run-time faster.
In virtually all cases, explicit type annotations are irrelevant for performance.
Type annotations in function signatures define a type filter/user interface.