10f. SIMD: Branchless Code

Introduction

SIMD accelerates workloads by applying a single instruction to many elements at once. Its effectiveness critically hinges on all elements following the same control flow. Yet, this pattern can be disrupted when conditionals are introduced, either explicitly in the source code or implicitly by the compiler. And, while the compiler attempts to mitigate the issue by transforming code into SIMD-friendly forms, these adaptations tend to introduce additional work that reduce the net benefit of SIMD.

The current section explores strategies for efficiently applying SIMD in the presence of conditional behavior. In particular, we'll show how to restructure computations to make control-flow branches disappear, allowing the program to preserve a uniform instruction stream.

We'll first examine scenarios in which the compiler introduces branches as a byproduct of its internal computation techniques. In these situations, users may be unaware that conditional control flow has been generated, making it essential to recognize when such hidden branches occur.

After this, we consider explicit conditional statements, written directly by the user. These conditions are intrinsic to program logic, ensuring that certain computations occur only when specific criteria are met. Because they can't simply be eliminated, we revisit common approaches to expressing conditions and analyze their internal implementation. The aim is to compare their relative strengths and limitations, indicating which ones interact more favorably with SIMD execution.

Finally, we show that conditional statements can often be reformulated as algebraic operations. This transformation removes branching logic and restores SIMD’s ability to operate on a uniform instruction stream. In fact, compilers may implement such transformations automatically when it's likely beneficial.

Approaches to Conditional Statements

Certain workflows require conditional statements, making them an integral part of a program's logical flow. Because they can't be avoided, it's important to consider the most effective way to express them.

We've already examined several ways to introduce conditional statements. In the context of SIMD, a relevant distinction among these approaches is whether they evaluate conditional statements eagerly or lazily.

To illustrate, let's consider computing 1 + 1 only when certain condition C is met. A lazy strategy checks C first, performing the addition only if the condition is true. By contrast, an eager approach performs 1 + 1 immediately, regardless of whether C is satisfied. If the condition eventually evaluates to false, the result of 1 + 1 is simply discarded.

When a conditional is evaluated just once, a lazy approach is more performant because it avoids unnecessary work. In contrast, when that conditional appears inside a for-loop and is evaluated repeatedly, SIMD can apply the same operation across many elements simultaneously. In these cases, it's often advantageous to compute all branches upfront, selecting the relevant results only afterward. This eager strategy works particularly well when the branches involve inexpensive computations.

In Julia, whether a conditional expression is evaluated lazily or eagerly depends on how the code is written. Next, we'll examine this aspect in more detail.

ifelse vs if

The ifelse function in Julia follows an eager evaluation strategy, computing all possible outcomes before selecting the one that matches the condition. In contrast, the if construct favors lazy computations, where only the branch corresponding to the true condition is executed.

The following example demonstrates the computational difference between if and ifelse. It uses a reduction operation in which the elements to be summed are contingent on a condition. Because ifelse requires an argument for both the true and false cases, we return zero when the condition isn't met.

x = rand(1_000_000)

function foo(x)
    output = 0.0

    for i in eachindex(x)
        if x[i] > 0.5
            output += x[i]/2
        end
    end

    return output
end

  1.312 ms (0 allocations: 0 bytes)

x = rand(1_000_000)

function foo(x)
    output = 0.0

    @inbounds @simd for i in eachindex(x)
        if x[i] > 0.5
            output += x[i]/2
        end
    end

    return output
end

  1.321 ms (0 allocations: 0 bytes)

x = rand(1_000_000)

function foo(x)
    output = 0.0

    for i in eachindex(x)
        output += ifelse(x[i] > 0.5, x[i]/2, 0)
    end

    return output
end

  1.221 ms (0 allocations: 0 bytes)

x = rand(1_000_000)

function foo(x)
    output = 0.0

    @inbounds @simd for i in eachindex(x)
        output += ifelse(x[i] > 0.5, x[i]/2, 0)
    end

    return output
end

  349.996 μs (0 allocations: 0 bytes)

The example shows that ifelse doesn't automatically trigger SIMD execution, even though it uses eager evaluation. This is precisely why the code required the @inbounds @simd macros.

It's also worth noting that applying SIMD instructions with ifelse doesn't necessarily improve performance. The following example demonstrates this point.

x      = rand(5_000_000)
output = similar(x)

function foo!(output,x)
    for i in eachindex(x)
        output[i] = ifelse(x[i] > 0.5, x[i]/2, 0)
    end
end

  56.122 ms (0 allocations: 0 bytes)

x      = rand(5_000_000)
output = similar(x)

function foo!(output,x)
    @inbounds @simd for i in eachindex(x)
        output[i] = ifelse(x[i] > 0.5, x[i]/2, 0)
    end
end

  50.983 ms (0 allocations: 0 bytes)

Ternary Operator

An alternative approach to expressing conditions is given by the ternary operation, whic follows the syntax <condition> ? <action if true> : <action if false>. Unlike the previous methods, this form doesn't enforce an eager or a lazy approach. Instead, it relies on compiler heuristics to determine which strategy will likely be faster for the specific application.

To illustrate this behavior, we revisit the same reduction example used earlier, now applying @inbounds @simd directly within each construct. In the following, the ternary operator chooses an eager strategy, thus evaluating both expressions before selecting the relevant branch.

x = rand(1_000_000)

function foo(x)
    output = 0.0

    @inbounds @simd for i in eachindex(x)
        if x[i] > 0.5
            output += x[i]/2
        end
    end

    return output
end

  1.318 ms (0 allocations: 0 bytes)

x = rand(1_000_000)

function foo(x)
    output = 0.0

    @inbounds @simd for i in eachindex(x)
        output += ifelse(x[i]>0.5, x[i]/2, 0)
    end

    return output
end

  386.313 μs (0 allocations: 0 bytes)

x = rand(1_000_000)

function foo(x)
    output = 0.0

    @inbounds @simd for i in eachindex(x)
        output += x[i]>0.5 ? x[i]/2 : 0
    end

    return output
end

  392.661 μs (0 allocations: 0 bytes)

On the contrary, the ternary operator opts for a lazy approach in the following example.

x = rand(1_000_000)

function foo(x)
    output = 0.0

    @inbounds @simd for i in eachindex(x)
        if x[i] > 0.99
            output += log(x[i])
        end
    end

    return output
end

  1.267 ms (0 allocations: 0 bytes)

x = rand(1_000_000)

function foo(x)
    output = 0.0

    @inbounds @simd for i in eachindex(x)
        output += ifelse(x[i] > 0.99, log(x[i]), 0)
    end

    return output
end

  11.496 ms (0 allocations: 0 bytes)

x = rand(1_000_000)

function foo(x)
    output = 0.0

    @inbounds @simd for i in eachindex(x)
        output += x[i]>0.99 ? log(x[i]) : 0
    end

    return output
end

  1.358 ms (0 allocations: 0 bytes)

It's worth noting that the method chosen by the ternary operator isn't foolproof, possibly choosing a less performant approach. This is illustrated below.

x      = rand(5_000_000)
output = similar(x)

function foo!(output,x)
    @inbounds @simd for i in eachindex(x)
        if x[i] > 0.5
            output[i] = log(x[i])
        end
    end
end

  79.400 ms (0 allocations: 0 bytes)

x      = rand(5_000_000)
output = similar(x)

function foo!(output,x)
    @inbounds @simd for i in eachindex(x)
        output[i] = ifelse(x[i] > 0.5, log(x[i]), 0)
    end
end

  54.072 ms (0 allocations: 0 bytes)

x      = rand(5_000_000)
output = similar(x)

function foo!(output,x)
    @inbounds @simd for i in eachindex(x)
        output[i] = x[i]>0.5 ? log(x[i]) : 0
    end
end

  80.596 ms (0 allocations: 0 bytes)

Vectors With Conditions

In many situations, you may already have a vector of Boolean conditions defined. This can happen because the condition vector is part of the dataset itself, or because the same conditions will be reused throughout the code, making it worthwhile to compute and store them in advance.

Julia represents such vectors using either Vector{Bool} or BitVector. The latter is the default output for broadcast expressions like x .> 0. While BitVector is memory-efficient and can be faster in some contexts, it can also inhibit SIMD optimizations. When this occurs, converting the BitVector to a Vector{Bool} can improve performance. The example below illustrates this effect: even after accounting for the cost of converting the type, the version using Vector{Bool} achieves a shorter execution time.

x         = rand(1_000_000)
bitvector = x .> 0.5

function foo(x,bitvector)
    output     = similar(x)
        

    @inbounds @simd for i in eachindex(x)
        output[i] = ifelse(bitvector[i], x[i]/i, x[i]*i)
    end

    return output
end

  13.156 ms (3 allocations: 7.629 MiB)

x         = rand(1_000_000)
bitvector = x .> 0.5

function foo(x,bitvector)
    output     = similar(x)
    boolvector = Vector(bitvector)
    
    @inbounds @simd for i in eachindex(x)
        output[i] = ifelse(boolvector[i], x[i]/i, x[i]*i)
    end

    return output
end

  2.588 ms (5 allocations: 8.583 MiB)

No Vector Holding Conditions

When the vector of conditions isn't already available, the conclusions outlined above no longer apply. In those cases, the most efficient strategy to leverage SIMD is to avoid building a condition vector altogether. Instead, you should apply ifelse directly to scalar values inside the for-loop. This approach eliminates unnecessary memory allocations, while still allowing the compiler to generate SIMD instructions effectively. The following example demonstrates this point.

BitVector Vector{Bool} No Vector

x = rand(1_000_000)

function foo(x)
    output     = similar(x)
    bitvector  = x .> 0.5
    
    
    @inbounds @simd for i in eachindex(x)
        output[i] = ifelse(bitvector[i], x[i]/i, x[i]*i)
    end

    return output
end

Output in REPL

julia>

foo($x)

  9.641 ms (7 allocations: 7.749 MiB)

x = rand(1_000_000)

function foo(x)
    output     = similar(x)
    boolvector = Vector{Bool}(undef,length(x))
        boolvector .= x .> 0.5

    @inbounds @simd for i in eachindex(x)
        output[i] = ifelse(boolvector[i], x[i]/i, x[i]*i)
    end

    return output
end

Output in REPL

julia>

foo($x)

  2.314 ms (5 allocations: 8.583 MiB)

x = rand(1_000_000)

function foo(x)
    output     = similar(x)
    
    

    @inbounds @simd for i in eachindex(x)
        output[i] = ifelse(x[i]>0.5, x[i]/i, x[i]*i)
    end

    return output
end

Output in REPL

julia>

foo($x)

  1.436 ms (3 allocations: 7.629 MiB)