# 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
 

### More realistic case: Vandermonde matrix
(modified from [Steve's Julia intro](https://web.mit.edu/18.06/www/Fall17/1806/julia/Julia-intro.pdf))

[Vandermonde matrix:](https://en.wikipedia.org/wiki/Vandermonde_matrix)
\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 a simple Julia implementation


```julia
function vander(x::AbstractVector{T},  n=length(x)) where T
    m = length(x)
    V = Matrix{T}(undef, m, n)
    for j = 1:m
        V[j,1] = one(x[j])
    end
    for i= 2:n
        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.**

When a function is called for the first time, Julia compiles efficient machine code for the given input types.

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) = x^2 + y
```


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

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

**Second call:** running the code


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

If the input types change, Julia compiles a new specialization of the function.


```julia
@time func(1.3,4.8)
@time func(1.3,4.8)
```

We now have two efficient codes, one for all `Int64` inputs and another one for all `Float64` arguments, in the cache.

### *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);
@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 run-time the fun time.

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 run-time 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.

A general strategy is therefore to move parts of the computation to compile-time.

Since Julia specializes on types, at compile-time **only type information is available to the compiler.**


```julia
f1(x::Int) = x + 1
f2(x::Int) = x + 2

function f_slow(x::Int, p::Bool)
    if p                                # check depends on the value of p
        return f1(x)
    else
        return f2(x)
    end
end
```


```julia
@code_llvm debuginfo=:none f_slow(1, true)
```

We can eliminate the if branch by moving the condition check to the type domain. This way, it **will only be evaluated once at compile-time.**


```julia
abstract type Boolean end
struct True <: Boolean end # type domain true
struct False <: Boolean end # type domain false

function f_fast(x::Int, p::Boolean)
    if typeof(p) == True                # check solely based on the type of p
        return f1(x)
    else
        return f2(x)
    end
end
```


```julia
@code_llvm debuginfo=:none f_fast(1, True())
```

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

53783c77e44ad2b688f9bd51ea097e71be93f34a
ipynb
Jupyter Notebook
Day1/3_specialization.ipynb
de-souza/JuliaOulu20
f9685a8b769a66d8f1145d7428bb4894918830b0
MIT
2019-12-05T10:33:42.000Z
2021-07-13T15:33:55.000Z
Day1/3_specialization.ipynb
de-souza/JuliaOulu20
f9685a8b769a66d8f1145d7428bb4894918830b0
MIT
2020-08-03T03:13:27.000Z
2020-08-03T03:13:27.000Z
Day1/3_specialization.ipynb
de-souza/JuliaOulu20
f9685a8b769a66d8f1145d7428bb4894918830b0
MIT
2020-02-11T08:28:07.000Z
2021-06-18T06:21:50.000Z
Qwen/Qwen-72B
1. YES
2. YES