10f. SIMD: Branchless Code

Branches

SIMD relies on the ability to implement the same set of instructions on multiple data. This parallel processing allows operations like scalar additions to be performed on several elements at once, significantly boosting performance. However, conditional statements pose a particular challenge for SIMD execution. When code branches into different paths and each requires distinct instructions, it becomes difficult to maintain the "same instruction" requirement that SIMD depends on. And, although compilers can sometimes transform such code to make it SIMD-compatible, these adaptations often come at the cost of reduced efficiency.

This section will focus on the handle of conditional operations. Firstly, some conditional statements could be inherent to the nature of our program, and hence unavoidable. This occurs when different operations are executed based on the state of the program. In these cases, we'll show that the approaches to implementing conditionals are dealt different internally. Depending on the workflow at hand, this results in advantages and disadvantages for each approach, with some being more SIMD-friendly than others.

Secondly, we'll consider cases where conditional statements can be eliminated entirely. This reflects scenarios where our code introduces conditions implicitly, which can be bypassed by writing code in an alternative way. In other scenarios, we can avoid conditional statements by transforming a conditional statement into an algebraic operation.

Approaches to Conditional Statements

Conditional statements can be evaluated eagerly or lazily. To illustrate this distinction, let's consider a simple example where we need to compute the result 1 + 1, but only if some condition C is met. A lazy approach first evaluates whether C holds true, before proceeding with the computation of 1 + 1. Thus, the operation is deferred until it's confirmed that C is satisfied. In contrast, an eager approach computes 1 + 1, regardless of whether C is satisfied or not. After this, if C turns out to be false, the computation remains unused.

When the conditional statement is applied only once, a lazy approach is always more performant as avoids unnecessary computations. However, neither approach consistently outperforms the other when the conditional statement is part of a for-loop. This is largely due the way modern CPUs handle simultaneous computations through SIMD. If the CPU can compute multiple operations simultaneously, it's possible that computing all conditions and branches upfront, selecting the relevant branches only afterwards, could actually result in performance gains. This is especially true when the branches consist of simple operations.

In Julia, whether a conditional statement is evaluated lazily or eagerly depends on the method used to compute it. Next, we explore this topic in more detail.

Ternary Operators

Ternary operators are an alternative approach for conditional statements, consisting of the form <condition> ? <action if true> : <action if false>. Unlike the previous methods, this form relies on heuristics to determine whether an eager or lazy approach should be used, depending on which would more likely be faster in the application considered.

For the illustrations, we'll consider examples where we directly add @inbounds and @simd to each approach.

Different Choices

Starting with the same example as above, we show that the ternary operator could choose an eager approach.

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

  422.480 μs (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

  85.895 μ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

  87.881 μs (0 allocations: 0 bytes)

Instead, the ternary operator implements 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

  405.304 μs (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

  3.470 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

  421.493 μs (0 allocations: 0 bytes)

Ternary Operator Could Choose A Less Performant Approach

It's worth remarking that the method chosen by the ternary operator isn't foolproof. In the following scenario, it actually chooses the slowest approach.

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

  26.620 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

  16.864 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

  26.517 ms (0 allocations: 0 bytes)

Vector of Conditions

Next, we consider scenarios where you already have a vector holding conditions. This could occur either because the vector is already part of your dataset, or because the conditions will be reused multiple times over you code, in which case storing the conditions is worthy.

Storing conditions in a vector could be done through an object with type Vector{Bool} or BitVector. The latter is the default type returned by Julia, as when you define objects like x .> 0. Although this type offer certain performance advantages, it can also hinder the application of SIMD. In cases like this, transforming BitVector to Vector{Bool} could speed up computations.

The following example demonstrates this, where the execution time is faster even considering the vector transformation.

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

  3.393 ms (2 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

  862.798 μs (4 allocations: 8.583 MiB)

No Vector of Conditions

The conclusions stated here assumes that you already have the vector holding the conditions. If this isn't the case and you want to apply SIMD instructions, you should implement ifelse without a vector of conditions. This allows you to avoid memory allocations, while still applying SIMD effectively. The following example illustrates this point. [note] Note that the approach for Vector{Bool} is somewhat different to the examples we considered above. As we don't have a vector of conditions already defined, it's optimal to create Vector{Bool} directly, rather than defining it as a transformation of the BitVector. In this way, we avoid unnecessary memory allocations too.

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)

  3.628 ms (6 allocations: 7.753 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)

  774.952 μs (4 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)

  501.114 μs (2 allocations: 7.629 MiB)