Type Stability

While Julia is a dynamic language, it can be very important for the compiler to know the type of the variables being used, so that it can specialise on it. We have been alluding to this throughout this course so far, but getting this wrong may have huge performance impacts. Let’s take a concrete example:

julia

function get_fibonacci(n)
    a = 1
    b = 1
    fibonacci_nums = []
    push!(fibonacci_nums, a)
    push!(fibonacci_nums, b)
    for i = 1:n-2
        a, b = b, a+b
        push!(fibonacci_nums, b)
    end
    return fibonacci_nums
end

This calculates the first $n$ Fibonacci numbers and stores the results in an array. Note that [] is shorthand for Vector{Any}() — an array that can hold anything — so even though we only ever push integers into it, the result is a Vector{Any}. Let’s confirm that, and then see what effect it has on performance:

julia

fibonacci_nums = get_fibonacci(50);
@show typeof(fibonacci_nums);

Output

typeof(fibonacci_nums) = Vector{Any}

julia

import BenchmarkTools: @benchmark, @btime, @belapsed
display(@benchmark sum($fibonacci_nums))

Output

BenchmarkTools.Trial: 10000 samples with 59 evaluations per sample.
 Range (min … max):  849.153 ns …  69.536 μs  ┊ GC (min … max): 0.00% … 97.30%
 Time  (median):     902.034 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   978.995 ns ± 984.329 ns  ┊ GC (mean ± σ):  1.91% ±  2.87%

   ▅██▅▅▂▂▃▂▄▆▅▁       ▁         ▁                              ▂
  █████████████████▇██████▇▅▆▆▆▆███▆▇▅▅▆▆▆▇▅▅▆▆▇▅▅▆▃▅▄▄▃▆▅▄▃▁▄▆ █
  849 ns        Histogram: log(frequency) by time       1.77 μs <

 Memory estimate: 608 bytes, allocs estimate: 38.

We can see that this sum had many allocations, and took a significant amount of time. If we simply make another array with the correct type, let us see the difference in performance:

julia

fibonacci_nums_with_type = [f for f in fibonacci_nums];
@show typeof(fibonacci_nums_with_type);

Output

typeof(fibonacci_nums_with_type) = Vector{Int64}

julia

display(@benchmark sum($fibonacci_nums_with_type))

Output

BenchmarkTools.Trial: 10000 samples with 999 evaluations per sample.
 Range (min … max):  7.818 ns … 45.265 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     8.268 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   8.794 ns ±  1.540 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▆█▅   ▆▇▇▄▄▂    ▂▂                                        ▂
  ████▄▄▇██████▆▆▆▆██▅█▅▅▅▆▄▆▅▅▅▄▄▅▅▄▅▄▂▅▅▄▅▇▅█▇▇▇▇▆▇▇▅▆▆▆▅▅ █
  7.82 ns      Histogram: log(frequency) by time     15.6 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Notice that the allocations were reduced to zero. Additionally, the runtime of these algorithms improved dramatically. In this example, when Julia knew the type of the variables, the performance increased by around two orders of magnitude. The reason why performance was so low for the first method, is that Julia had to check the type of each individual element in the array to make sure that it was the right type.

Understanding Type Instability

If the compiler cannot trace the types of the variables used throughout a function (each variable having a known concrete type at compile time), then the function is said to be type unstable. This means that the compiler will have to perform additional checks on variables with unknown types at runtime, which costs performance.

Let’s take a look at an example of a type unstable function:

julia

function example_type_unstable_fn(x)
    s = 0
    for x_i in x
        s += x_i
    end
    s
end

This looks like a normal implementation. We can even run it and benchmark to see if it works:

julia

x = collect(1:100);
display(@benchmark example_type_unstable_fn($x))

Output

BenchmarkTools.Trial: 10000 samples with 999 evaluations per sample.
 Range (min … max):  10.741 ns … 74.434 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     11.271 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   11.804 ns ±  2.403 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▆▇█▇▅▁▂             ▁                                    ▁  ▂
  ███████▇▅▆▇▆████▆▆▇██▇▇▆▇▇▇▇▇▇▇▅▅▅▆▅▆▆▅▅▄▅▆▇█▃▇▅▇▁▄▁▃▃▃▆██▆ █
  10.7 ns      Histogram: log(frequency) by time      22.8 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

There are actually no issues with this code so far that we can see, but let’s show an example where some performance issues come out:

julia

function calc_summary(rows)
    return sum(example_type_unstable_fn, rows)
end

rows = [rand(rand(0:10)) for _ in 1:100];
display(@benchmark calc_summary($rows))

Output

BenchmarkTools.Trial: 10000 samples with 276 evaluations per sample.
 Range (min … max):  283.478 ns …  1.836 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     294.783 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   302.076 ns ± 34.686 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

    ▅█▅▂                                                        
  ▄██████▅▄▃▃▃▃▃▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂ ▃
  283 ns          Histogram: frequency by time          439 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

We can use the @code_warntype macro to identify type instabilities. This macro lives in the InteractiveUtils standard library, which the REPL loads for you automatically, but which you must import explicitly inside a script:

julia

using InteractiveUtils
@code_warntype example_type_unstable_fn(rows[1])

Output

MethodInstance for example_type_unstable_fn(::Vector{Float64})
  from example_type_unstable_fn(x) @ Main ~/dev/website-dev/personal-website/.code-exec/_53f867569c9c8c6a2fbfb8e5cae2618489feecdf3ce20c5bc2d9af61eee5150b.jl:25
Arguments
  #self#::Core.Const(Main.example_type_unstable_fn)
  x::Vector{Float64}
Locals
  @_3::UNION{NOTHING, TUPLE{FLOAT64, INT64}}
  s::UNION{FLOAT64, INT64}
  x_i::Float64
Body::UNION{FLOAT64, INT64}
1 ─       (s = 0)
│   %2  = x::Vector{Float64}
│         (@_3 = Base.iterate(%2))
│   %4  = @_3::UNION{NOTHING, TUPLE{FLOAT64, INT64}}
│   %5  = (%4 === nothing)::Bool
│   %6  = Base.not_int(%5)::Bool
└──       goto #4 if not %6
2 ┄ %8  = @_3::Tuple{Float64, Int64}
│         (x_i = Core.getfield(%8, 1))
│   %10 = Core.getfield(%8, 2)::Int64
│   %11 = Main.:+::Core.Const(+)
│   %12 = s::UNION{FLOAT64, INT64}
│   %13 = x_i::Float64
│         (s = (%11)(%12, %13))
│         (@_3 = Base.iterate(%2, %10))
│   %16 = @_3::UNION{NOTHING, TUPLE{FLOAT64, INT64}}
│   %17 = (%16 === nothing)::Bool
│   %18 = Base.not_int(%17)::Bool
└──       goto #4 if not %18
3 ─       goto #2
4 ┄ %21 = s::UNION{FLOAT64, INT64}
└──       return %21

We can see that the variable s is type unstable (i.e. has a union type). This is because the input array has the possibility of being empty. If this array is empty then the loop gets skipped and the function will return 0, which is an integer. If the array is not empty, then the first assignment to s will promote the type to a floating point number.

Note: This particular type instability is not catastrophic in recent versions of Julia. When the instability covers a small union of concrete types the compiler can union split, generating a branch for each possibility, so you pay a modest cost rather than a full dynamic dispatch. It is still a real cost, as the benchmarks below show. More importantly, if this code is nested within a larger algorithm, the instability can propagate outwards and cause much larger issues.

We can fix this type instability by using generic functions, like zero and eltype:

julia

function example_type_stable_fn(x)
    s = zero(eltype(x))
    for x_i in x
        s += x_i
    end
    s
end

function calc_summary_stable(rows)
    return sum(example_type_stable_fn, rows)
end

Now let’s benchmark both versions:

julia

display(@benchmark calc_summary($rows))

Output

BenchmarkTools.Trial: 10000 samples with 269 evaluations per sample.
 Range (min … max):  283.978 ns …  1.993 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     303.755 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   321.343 ns ± 51.260 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

    ▃█▅       ▁                                                 
  ▂▅████▃▂▂▂▂▄█▆▆▂▁▁▁▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  284 ns          Histogram: frequency by time          521 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia

display(@benchmark calc_summary_stable($rows))

Output

BenchmarkTools.Trial: 10000 samples with 240 evaluations per sample.
 Range (min … max):  311.333 ns …  1.546 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     326.875 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   339.818 ns ± 52.038 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▆███▇▇▆▅▄▃▃▂▂▂▂▁                                             ▂
  ███████████████████▇▆▇▇▇▇▆▆▇▆▆▆▆▇▆▆▅▅▅▇▇█▇█▇▅▅▅▆▅▅▆▅▄▅▄▂▅▅▄▃ █
  311 ns        Histogram: log(frequency) by time       562 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

You can see a solid performance improvement when switching to the type stable version — roughly a factor of two here, with neither version allocating. Making our code type stable also ensures that the code is very generic and reusable across several different datatypes. In general, it is a good idea to avoid hard-coding numeric literals such as 0 and 1 in generic code, and instead build them from the input with the zero, one, eltype and typeof functions, so that they take the same type as the data you are working with.

Barriers

There are some occasions when you cannot predict the type being returned from a function, which leads to type instability. There are still some tools one can use to mitigate the performance hits of this. The antidote to this problem is to break up larger functions into several smaller functions, which act as type barriers so that performance critical parts of your code are not affected.

Let’s take the example of using an array to store parameters which are used within a function:

julia

function params_in_array_test(x)
    a = 0
    b = 1
    for i = 1:100
        c = x[1] + a
        d = x[2] / b
        a, b = c, d
    end
    return a, b
end

x = [2.0, 1]
display(@benchmark params_in_array_test($x))

Output

BenchmarkTools.Trial: 10000 samples with 846 evaluations per sample.
 Range (min … max):  138.570 ns … 399.846 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     146.974 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   151.188 ns ±  14.794 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▁   ▆█▇▇▆▇▇▄▄▄▄▄▆▂▁    ▁             ▁▂▅▂▁▂▁▁                 ▂
  █▇▅▇█████████████████████▇▆▇▅▅▆▆▇▆▆▅██████████▇▆▆▅▆▅▇█▆▇▆▅▆▇▆ █
  139 ns        Histogram: log(frequency) by time        193 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Note that the array constructor has “promoted” the integer 1 to a floating point number, so that every element shares a single concrete type: x here is a Vector{Float64}, not a mixed array. Now, let’s imagine that we add another parameter to this array, which is of a different type:

julia

x = [2.0, 1, "third"]
display(@benchmark params_in_array_test($x))

Output

BenchmarkTools.Trial: 10000 samples with 7 evaluations per sample.
 Range (min … max):  4.453 μs … 92.449 μs  ┊ GC (min … max): 0.00% … 86.79%
 Time  (median):     5.069 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   5.263 μs ±  3.275 μs  ┊ GC (mean ± σ):  2.71% ±  4.09%

              ▂▃▅█▃▃                                          
  ▁▂▂▁▂▂▂▄▆▇▆████████▅▄▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  4.45 μs        Histogram: frequency by time        6.85 μs <

 Memory estimate: 3.16 KiB, allocs estimate: 201.

Now, even though the first two entries of x are unchanged and all we did was add a third parameter, the function has become more than an order of magnitude slower, and it now allocates on every iteration of the loop where before it allocated nothing at all. This is because x is now a Vector{Any}, so the compiler no longer knows the types of x[1] and x[2], and every operation in the loop has to be dispatched dynamically and boxed on the heap.

We can remedy this hit, without changing the input types, by using a barrier function:

julia

function _params_in_array_test_loop(x1, x2, a, b)
    for i = 1:100
        c = x1+a
        d = x2 / b
        a, b = c, d
    end
    return a, b
end

function params_in_array_test_with_barrier(x)
    a = 0
    b = 1
    return _params_in_array_test_loop(x[1], x[2], a, b)
end

display(@benchmark params_in_array_test_with_barrier($x))

Output

BenchmarkTools.Trial: 10000 samples with 575 evaluations per sample.
 Range (min … max):  200.522 ns …  2.791 μs  ┊ GC (min … max): 0.00% … 88.86%
 Time  (median):     212.348 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   215.386 ns ± 53.508 ns  ┊ GC (mean ± σ):  0.52% ±  1.97%

     █       ▁    ▃                                             
  ▁▆▆█▅▃▃▇█▅▅█▅▅▆▇█▆▄▄▄▄▅▄▆▄▄▃▃▄▃▃▃▂▂▃▃▂▂▂▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  201 ns          Histogram: frequency by time          246 ns <

 Memory estimate: 32 bytes, allocs estimate: 1.

This recovers almost all of the lost performance: the loop now runs in a type stable manner and only a single dynamic dispatch remains, at the point where the untyped values are pulled out of x and handed to the barrier function. This is the tool to reach for when it is not possible to remove the type instability entirely.

Using Composite Types

In Julia, we can define our own composite types using struct. By default, structs are immutable. Structs are a common place that people introduce type instability. Take the example of defining one’s own complex number:

julia

struct MyComplex
    real
    imag
end

However, this type is a source of type instability. By not annotating them, we have implicitly declared both the real and imag fields to have type Any, so the compiler cannot know what they hold, and an instance must store pointers to boxed values on the heap. We can fix this by specifying a type:

julia

struct MyComplexFixed
    real::Float64
    imag::Float64
end

However now, this type will only work with type Float64. Instead, we can use a feature of Julia called generics, which lets the field types be a parameter of the struct itself:

julia

struct MyComplexGeneric{T<:Number}
    real::T
    imag::T
end

# Example usage
c1 = MyComplexGeneric(1.0, 2.0)
c2 = MyComplexGeneric(1, 2)
@show typeof(c1) typeof(c2) isbitstype(typeof(c1));

Output

typeof(c1) = MyComplexGeneric{Float64}
typeof(c2) = MyComplexGeneric{Int64}
isbitstype(typeof(c1)) = true

If you have a performance critical struct, make sure that the types are well-defined. If every field of a struct has a concrete, fixed-size type (for example a collection of numeric fields), then the struct itself is an isbits type, and Julia is free to keep it in registers or on the stack rather than allocating it on the heap.

Back Next

High Performance Computing in Julia

1. Introduction

2. Foundations

3. Julia

4. Measuring Performance

5. Optimisation

6. Parallel Programming

7. Multithreading

8. Multiprocessing

9. GPU Programming

Table of Contents

Type Stability

Understanding Type Instability

Barriers

Using Composite Types