Writing Fast Julia

Writing Fast Julia#

Everything in this book has been building toward this chapter.

Julia’s promise is that you can write high-level code and get low-level speed. That promise is real, but it isn’t automatic. The gap between slow Julia and fast Julia is often 10x or 100x, and it usually comes down to a handful of specific, learnable mistakes. The good news is that Julia gives you excellent tools for finding them — so this chapter is less “memorize these rules” and more “here’s how to ask the compiler what it thinks.”

The one rule underneath all the others#

The compiler must be able to predict the type of every value.

That’s it. When Julia can infer that a variable holds a Float64, it emits a machine instruction that adds two floats. When it can only infer Any, it emits code that inspects the value at runtime, looks up which + method applies, allocates a box for the result, and moves on. The second is easily fifty times slower and allocates on every operation.

Code where every type is predictable is called type stable. Everything below is a specific way to achieve it, or a specific way people accidentally break it.

Rule 1: Put your work inside functions#

I planted this back in the control flow chapter. Here’s the payoff.

# In a script, at top level — SLOW
total = 0.0
for i in 1:10_000_000
    global total += sqrt(i)
end

# Inside a function — FAST
function accumulate_sqrt(n)
    total = 0.0
    for i in 1:n
        total += sqrt(i)
    end
    return total
end

The difference is often an order of magnitude, and the reason is simple: total at the top level is a global variable. Any other code, at any time, could reassign it to a string. The compiler must assume nothing about its type. Inside a function, total is a local — the compiler can see every assignment to it and knows it’s always a Float64.

Wrap your work in functions. If you take one piece of advice from this chapter, take that one.

If you genuinely need a global, make it const:

const TAX_RATE = 0.0825      # type is now known and fixed

Rule 2: Watch out for type instability#

A function is type unstable when its return type depends on the values of its arguments rather than their types.

# UNSTABLE: returns Int or Float64 depending on the value
function bad_sign(x)
    if x > 0
        return 1          # Int
    else
        return -1.0       # Float64
    end
end

Julia can only infer Union{Int64, Float64} here, and every caller pays for that uncertainty.

# STABLE
function good_sign(x)
    return x > 0 ? 1.0 : -1.0
end

A subtler and much more common version:

# UNSTABLE — accumulator starts as Int, becomes Float64
function sum_scaled(xs)
    total = 0                # Int64!
    for x in xs
        total += x * 1.5     # now it needs to be Float64
    end
    return total
end

# STABLE — start with the right type
function sum_scaled(xs)
    total = zero(eltype(xs)) * 1.5      # or just 0.0 if you know
    for x in xs
        total += x * 1.5
    end
    return total
end

Initializing an accumulator with the wrong type is, I’d guess, the single most common type instability in beginner Julia. Watch for it.

@code_warntype — your best diagnostic#

Julia will tell you exactly where inference failed:

julia> @code_warntype bad_sign(3)

The output lists every variable with its inferred type. Anything printed in red (or shown as Any or a big Union) is a problem. Anything in blue or plain is fine.

You don’t need to understand the whole output. Scan for red. That’s the skill.

For a more automated version, JET.jl analyses a whole call graph and reports instabilities without you having to inspect each function:

using JET
@report_opt my_function(data)

Rule 3: Give containers concrete element types#

bad = []                      # Vector{Any} — every element is a pointer
good = Float64[]              # Vector{Float64} — contiguous, dense

bad_dict = Dict()             # Dict{Any,Any}
good_dict = Dict{String,Int}()

This is the “abstract in signatures, concrete in storage” rule from the type system chapter, applied. The same principle covers struct fields:

struct Slow            # every field is Any
    a
    b
end

struct Fast
    a::Float64
    b::Int
end

And it covers a case that catches people out:

struct AlsoSlow
    values::Vector{Real}      # abstract element type — pointers again
end

struct AlsoFast{T<:Real}
    values::Vector{T}         # parametric — concrete once constructed
end

If you want a struct field to hold different numeric types in different instances, make the struct parametric. Don’t reach for an abstract field type.

Rule 4: Don’t allocate in hot loops#

We covered the mechanics in the last chapter; here’s the discipline.

# Allocates a temporary array on every iteration
function slow_norms(rows)
    out = Float64[]
    for i in axes(rows, 1)
        push!(out, sqrt(sum(rows[i, :] .^ 2)))     # rows[i,:] copies!
    end
    return out
end

# Preallocated output, views instead of copies
function fast_norms(rows)
    out = Vector{Float64}(undef, size(rows, 1))
    for i in axes(rows, 1)
        r = @view rows[i, :]
        out[i] = sqrt(sum(abs2, r))
    end
    return out
end

Three changes: preallocate out, use @view instead of slicing, and use sum(abs2, r) which fuses the squaring into the sum instead of materializing a squared array.

Measuring properly#

Do not trust @time for anything small. It measures a single run, includes compilation on the first call, and is noisy.

using BenchmarkTools

@btime sum($xs)              # runs many times, reports the minimum
@benchmark sum($xs)          # full statistics

Note the $. In a benchmark, $xs interpolates the variable so you’re not accidentally measuring the cost of looking up a global. Forgetting the dollar sign is the classic benchmarking mistake in Julia and will make everything look slower than it is.

@time remains useful for one thing: its allocation count. If @time says a function that should do no allocation allocated 4 million times, you have your lead.

Profiling: find the hot spot before optimizing it#

using Profile

@profile my_analysis(data)
Profile.print()

Better, use ProfileView.jl or the VS Code extension’s profiler for a flame graph. Look for the wide bars — those are where the time goes.

I want to state the obvious rule that everybody knows and nobody follows: profile first. Programmer intuition about where time goes is unreliable, mine very much included. I have spent afternoons carefully optimizing a function that turned out to account for 2% of the runtime.

The sharper tools#

Once your code is type stable and not allocating, there are a few more levers. Use them last, and measure each one.

@inbounds for i in eachindex(xs)     # skip bounds checks
    ...
end

@simd for i in eachindex(xs)         # permit vectorization
    ...
end

@fastmath a * b + c                  # relax floating-point strictness

Warning

These are genuinely dangerous @inbounds disables the check that catches out-of-bounds indexing. Get the index wrong and instead of a clean error you read arbitrary memory, or corrupt it. Only use it on loops where you can prove the indices are valid — and eachindex is how you prove it.

@fastmath permits reassociating floating-point operations, which changes results. Fine for a graphics kernel; not fine for a financial calculation or a numerically delicate algorithm.

I’d add: check whether you need these at all. The compiler often eliminates bounds checks on its own when it can see the loop is safe, and modern LLVM vectorizes many loops without @simd.

A worked optimization#

Let’s take a realistic function through the whole process. The task: for each account in a ledger, compute how far its balance deviates from the mean, in standard deviations.

Version 1 — the way you’d naturally write it first:

function zscores_v1(balances)
    results = []
    m = sum(balances) / length(balances)
    s = sqrt(sum((balances .- m) .^ 2) / length(balances))
    for b in balances
        push!(results, (b - m) / s)
    end
    return results
end

@btime on a million elements: slow, and lots of allocations. @code_warntype shows results is a Vector{Any}.

Version 2 — fix the container type:

function zscores_v2(balances)
    results = Float64[]
    ...
end

Better. Still allocating as the vector grows.

Version 3 — preallocate:

function zscores_v3(balances)
    n = length(balances)
    m = sum(balances) / n
    s = sqrt(sum(abs2, balances .- m) / n)
    results = Vector{Float64}(undef, n)
    for i in eachindex(balances)
        results[i] = (balances[i] - m) / s
    end
    return results
end

Now it’s one allocation for the output plus one temporary in the standard deviation.

Version 4 — remove the last temporary and let broadcasting fuse:

using Statistics

function zscores_v4(balances)
    m = mean(balances)
    s = std(balances; corrected = false)
    return @. (balances - m) / s
end

Shorter than version 1, and fast. That’s the arc I want you to notice: the readable version and the fast version converged. The optimization didn’t make the code uglier — the intermediate versions were uglier than both ends.

That is the thing I like most about optimizing Julia. In many languages, making code fast means making it worse. Here, it usually means making it clearer about what it actually does.

A checklist#

When something is slower than you expect, in order:

  1. Is the work inside a function?

  2. Are you benchmarking with @btime and $-interpolated variables?

  3. Does @code_warntype show red?

  4. Are your containers and struct fields concretely typed?

  5. Is @time reporting allocations you didn’t expect?

  6. Have you profiled to confirm you’re optimizing the right function?

  7. Only then: views, preallocation, @inbounds, @simd.

Try it yourself#

  1. Write the top-level loop version and the function version of summing sqrt(i) for a million values. Time both. Note the ratio.

  2. Write a type-unstable function and run @code_warntype on it. Find the red.

  3. Take zscores_v1 above and work it up to v4 yourself, benchmarking at each step.

  4. Benchmark sum(xs) with and without the $ in @btime. Explain the difference.

You can now write Julia that’s genuinely fast on one core. Next: using all of them.