vectorized operations vs loops

[true pebble]

coming from numpy, I thought that like A .+ B would be faster than A + B because it is vectorized. But I just read that in julia loops are faster. Then what is the point of A .+ B?

[true pebble]

why would you add an extra symbol to make your code slower

[uneven maple]

if A and B are arrays A .+ B is equivalent to A + B

[true pebble]

whats the point of A .+ B

then

[uneven maple]

in your case nothing

julia> x = rand(10);

julia> x .= sin.(x)
10-element Vector{Float64}:
 0.14254193236205576
 0.3701762327696455
 0.3623631091392756
 0.8036946167310883
 0.2740764904500096
 0.5152639962240843
 0.5124770784935754
 0.2804176660697805
 0.5793032650627981
 0.6758465850947173

[pseudo egret]

the . turns any operation that would normally work on individual "things" into an operation that works on Arrays/Matrices/Tensors of things

[true pebble]

so its the same as a for loop?

[pliant eagle]

Unlike python loop, julia loops are fast so writing a loop isn't something that you must avoid for performance

[true pebble]

wait

so if I define a function

and then can I add . in front of it and it automatically works?

[pseudo egret]

it's not explicitly forbidden/disabled for something like A .+ B, even though the + operator already has an overload for arrays of the same size

[pliant eagle]

Yes, that's beautiful, isn't it?

[pseudo egret]

exactly

[pseudo egret]

a bit, but it can go a bit further, too. For multidimensional arrays, it loops over the non-shared dimensions, so you can do stuff like

rand(5,3) .+ rand(5)
5×3 Matrix{Float64}:
 1.84255  1.12342   1.16323
 1.44004  1.49876   1.02867
 1.10553  0.667457  0.583658
 1.1654   0.90351   0.424651
 1.72491  1.13988   1.48848

^ adding a 5-vector to each of the columns of a 5x3 matrix

and since operators in julia are just functions, the . works for essentially any function

[true pebble]

would rand(3,5) .+ rand(5) work (thats what I'd expect in numpy, and the code above would fail)

[pseudo egret]

no but you can do rand(3,5) .+ rand(1,5)

or rand(3,5) .+ rand(5)'

[true pebble]

interesting thats like having last dimensions as batch dimensions

[pliant eagle]

That's indeed the case. It's not a rule but that's the case most of the times. This happends for multiple reasons one of which is the fact that Julia is column major so having the batch dimension last will result in faster code/loops.

[true pebble]

is this true for something like solve or factorize, if I do solve(rand(100,100,100, rand(100,100) what do I expect, the docs don't seem to mention this

[pliant eagle]

There is also something called "fused vectorized operations"(see this paragraph of the julia manual).

Sometimes you can't write code that contains loops, for example when doing array programming with CUDA.jl

You can "add BenchmarkTools", "@btime" and try to benchmark the different styles of code and see the performance difference and then comprare it to standard python or numpy.

Also sometimes they recommend you to write the function that works on scalars first and then leverage the dot and multiple dispatch to "vectorize" it. I think that in the end the dot operation will be lowered to a loop. That's also the same for numpy functions, if I'm not mistaken numpy calls are C calls and loops in C are fast.

[pliant eagle]

I tried it and it doesn't work as is. This is what I tried: randn(100,100,100) \ randn(100,100)
The solve operation is done using ''\', for example to solve Ax=b you can do A \ b
It's about getting used to Julia syntax and multiple dispatch

[true pebble]

interesting so after googling linear algebra doesn't support batched matrices but I don't yet know if that will matter or if for loop is fine

[pseudo egret]

if I'm not mistaken numpy calls are C calls and loops in C are fast
If applicable, numpy calls out to some BLAS implementation that is faster than "just loops". Julia Arrays do the same

e.g. multiplying big matrices can be made faster by going through the matrices is chunks, rather than just applying the definition of the matrix product directly

[pseudo egret]

generally, you don't have to try to avoid loops as hard as with python. In python, it can be desirable to be kind of wasteful with computation, if that means you can "squeeze" your problem into big matrices/tensors. In Julia, if the more "natural" way to solve a problem is with a lot of small matrix operations in a loop, that's totally fine

(for working with small matrices/vectors, you might also want to look at StaticArrays.jl)