# Specialization

To be fast, Julia needs to **specialize** code, that is compile specific native versions of the code. **The better the specialization the faster the code!** In the following we will investigate how Julia achieves good code specialization while retaining the power of generic programming.

## 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
 

### Vandermonde matrix (once again)
(modified from [Steven's Julia intro](https://web.mit.edu/18.06/www/Fall17/1806/julia/Julia-intro.pdf))

\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}


```julia
using PyCall
```


```julia
np = pyimport("numpy")
```


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

The source code for this function is [here](https://github.com/numpy/numpy/blob/v1.16.1/numpy/lib/twodim_base.py#L475-L563). It calls `np.multiply.accumulate` which is implemented in C [here](https://github.com/numpy/numpy/blob/deea4983aedfa96905bbaee64e3d1de84144303f/numpy/core/src/umath/ufunc_object.c#L3678). However, this code doesn't actually perform the computation, it basically only checks types and stuff. The actual kernel that gets called is [here](https://github.com/numpy/numpy/blob/deea4983aedfa96905bbaee64e3d1de84144303f/numpy/core/src/umath/loops.c.src#L1742). 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 our simple generic Julia implementation


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


```julia
vander(1:5)
```

#### A quick speed comparison

<details>
  <summary>Show Code</summary>
<br>
    
```julia
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)
```
</details>

 

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!


```julia
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.


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

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

# How can Julia be fast?

<p></p>
 
**AST = Abstract Syntax Tree**

**SSA = Static Single Assignment**

**[LLVM](https://de.wikipedia.org/wiki/LLVM) = Low Level Virtual Machine**

### Specialization and code inspection

**Julia specializes on the types of function arguments**, i.e. Julia compiles efficient machine code for the given input types, **when a function is called for the first time**.

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


```julia
func(x,y) = 2x + y
```


```julia
x = [1.2, 3.4, 5.6] # Vector{Float64}
y = [0.4, 0.7, 0.9] # Vector{Float64}

@time func(x,y);
@time func(x,y);
```

**First call:** compilation + running the code

**Second call:** running the code


```julia
@time func(x,y);
```

If one of the input types changes, Julia compiles a new specialization of the function!


```julia
typeof(x)
```


```julia
x = [1, 3, 5]
```


```julia
typeof(x)
```


```julia
@time func(x,y); # Vector{Int64}, Vector{Float64}
@time func(x,y);
```

We now have two efficient native codes in the cache: one for all `Vector{Float64}` inputs and another one for `Vector{Int64}` as the first and `Vector{Float64}` as the second argument type.

### *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`**)


```julia
@code_typed func(1,2)
```


```julia
@code_lowered func(1,2)
```


```julia
@code_llvm func(1,2)
```

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


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


```julia
@code_native debuginfo=:none func(1,2)
```

Let's compare this to `Float64` input.


```julia
@code_native debuginfo=:none 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.)


```julia
struct Anything
    value::Any
end

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


```julia
using BenchmarkTools

@btime operation(2.0);

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

**That's about an 40 times slowdown!**


```julia
@code_native debuginfo=:none operation(2.0)
```


```julia
@code_native debuginfo=:none operation(x)
```

# Make runtime the fun time.

Julia specializes on input types (not values). As a rule of thumb: **only type information is available to the compiler when specializing code.**

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 runtime 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.

### Example: Determinant of a 2x2 matrix

Let's say your task would be to write a function computing the determinant of a 2x2 matrix. How would you implement it?

Probably you'd say, well I know the formula for computing the determinant of a 2x2 matrix! Let's just implement it.


```julia
det_2x2(X::AbstractMatrix) = X[1,1] * X[2,2] - X[1,2] * X[2,1]
```


```julia
M = [1 2; 3 4]
```


```julia
det_2x2(M)
```


```julia
@btime det_2x2(M);
```

Let's see how Julia's built-in `det` function compares to our algorithm:


```julia
using LinearAlgebra

det(M)
```


```julia
@btime det(M);
```

It's much slower!! But why?

The reason is that, as we've discussed above, the compiler only has the type information available for producing a specialization of the `det` function, i.e. in this case:


```julia
typeof(M)
```

Note that, in particular, the **size of the matrix is not encoded in the type and therefore not available to the compiler!**


```julia
size(M)
```


```julia
size(typeof(M))
```

So, in spirit, when calling `det(M)` for the first time, we ask the compiler to create a specialization which can handle matrices of all sizes! That's obviously much harder than the 2x2 case. Hence, the produced code is more general - it will also work for 3x3 matrices etc. - but much less efficient.


```julia
@code_typed debuginfo=:none det(M)
```

Let's now move the size information to the type domain and see how things change.


```julia
using StaticArrays
```


```julia
S = SMatrix{2,2}(1, 2, 3, 4)
```

or equivalently


```julia
S = @SMatrix [1 2; 3 4]
```

For static arrays, we can extract the matrix size solely from the type:


```julia
size(typeof(S)) # doesn't error
```

Hence, the compiler can utilize the fact that `S` is a 2x2 matrix - it can create different specializations for different matrix sizes.


```julia
@btime det(S);
```

In fact, let us check what `det(S)` is actually doing under the hood.


```julia
@code_typed debuginfo=:none det(S)
```

Let's translate this into something more readable:
```julia
%2 = %3 = S[1,1] = 1
%4 = %5 = S[2,1] = 3
%6 = %7 = S[1,2] = 2
%8 = %9 = S[2,2] = 4

%10 = %3 * %9 = S[1,1] * S[2,2]
%11 = %7 * %5 = S[1,2] * S[2,1]
%12 = %10 - %11 = S[1,1] * S[2,2] - S[1,2] * S[2,1]
```

**Overall it just corresponds to the explicit formula that we hand-coded in `det_2x2`!**
```julia
det(S) = S[1,1] * S[2,2] - S[1,2] * S[2,1]
```

(Side remark: of course, here it is not really the compiler that does the job but it gets help from multiple dispatch: StaticArrays.jl just implements a specific method of the function `det` for the 2x2 case, see `@edit det(S)`. But you hopefully get the point.)


```julia
@edit det(S)
```

**Other operations sped up by using StaticArrays.jl**

```
============================================
    Benchmarks for 3×3 Float64 matrices
============================================
Matrix multiplication               -> 5.9x speedup
Matrix multiplication (mutating)    -> 1.8x speedup
Matrix addition                     -> 33.1x speedup
Matrix addition (mutating)          -> 2.5x speedup
Matrix determinant                  -> 112.9x speedup
Matrix inverse                      -> 67.8x speedup
Matrix symmetric eigendecomposition -> 25.0x speedup
Matrix Cholesky decomposition       -> 8.8x speedup
Matrix LU decomposition             -> 6.1x speedup
Matrix QR decomposition             -> 65.0x speedup
```

Of course, by putting more information in the type you are putting more stress on the compiler to optimize things. If the static arrays are too big compile time might explode or the compiler might just give up and fall back to a slow default version.

Hence, static arrays are only useful as small fixed-size arrays.


```julia
# # might take longer to compile and the speedup is gone
# N = 20
# M = rand(N,N);
# m = SMatrix{N,N}(M);

# println("Inversion")
# @btime inv($m);
# @btime inv($M);
```

# 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!


```julia
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
```


```julia
@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**.

b3f54d238d534650d4f9be951258a388b611b1c8
ipynb
Jupyter Notebook
Day1/4_specialization.ipynb
crstnbr/JuliaNRWSS21
52ce594ea2986d4cb670ab2b704e231404a63617
MIT
2021-06-10T18:53:50.000Z
2021-06-22T10:15:56.000Z
Day1/4_specialization.ipynb
carstenbauer/JuliaNRWSS21
52ce594ea2986d4cb670ab2b704e231404a63617
MIT
2021-06-16T14:50:16.000Z
2021-06-16T14:50:16.000Z
Day1/4_specialization.ipynb
carstenbauer/JuliaNRWSS21
52ce594ea2986d4cb670ab2b704e231404a63617
MIT
2021-06-16T09:34:17.000Z
2021-06-16T09:34:17.000Z
Qwen/Qwen-72B
1. YES
2. YES