Why enforcing order onto Julia matrix vector multiplications has speed up?

user_3187369 0
Viewed 226

I've tried the following.

a = randn(100,100);
b = randn(100,100);
c = randn(100,1);

@time a*b*c
@time a*(b*c)

Results

julia> @time a*b*c;
  0.000591 seconds (7 allocations: 79.234 KiB)

julia> @time a*(b*c);
  0.000101 seconds (6 allocations: 1.906 KiB)

Results are quite consistent with the above. While it does make intuitive sense why the the second one is better (A matrix-vector multiplication twice instead of a large matrix-matrix multiplication).

I was wondering, shouldn't Julia have optimized this knowing that the dimensions of the matrices and it could've reordered the operation to have optimized it? Or am I just being too lazy to wish for that or are there other technical issues that I don't see.

So, this is what I get when I use dump() on (a*b*c)

  head: Symbol call
  args: Array{Any}((4,))
    1: Symbol *
    2: Array{Float64}((100, 100)) [0.290788 -0.0601455 … -0.408164 1.16261; -0.539274 -1.56979 … 2.56233 0.806247; … ; 1.30981 -1.31929 … 1.38655 -1.89169; -1.58483 0.318804 … -0.0500151 2.13105]
    3: Array{Float64}((100, 100)) [-0.464882 1.60371 … -0.390234 0.605401; -1.06837 0.296049 … 0.759708 0.0124688; … ; -0.149613 -1.38653 … 0.284494 1.47524; 0.34351 0.420449 … 0.544973 1.85736]
    4: Array{Float64}((100, 1)) [1.64066; 0.593296; … ; 0.908361; 0.486164]
  typ: Any

dump(a*(b*c))

Expr
  head: Symbol call
  args: Array{Any}((3,))
    1: Symbol *
    2: Array{Float64}((100, 100)) [0.290788 -0.0601455 … -0.408164 1.16261; -0.539274 -1.56979 … 2.56233 0.806247; … ; 1.30981 -1.31929 … 1.38655 -1.89169; -1.58483 0.318804 … -0.0500151 2.13105]
    3: Expr
      head: Symbol call
      args: Array{Any}((3,))
        1: Symbol *
        2: Array{Float64}((100, 100)) [-0.464882 1.60371 … -0.390234 0.605401; -1.06837 0.296049 … 0.759708 0.0124688; … ; -0.149613 -1.38653 … 0.284494 1.47524; 0.34351 0.420449 … 0.544973 1.85736]
        3: Array{Float64}((100, 1)) [1.64066; 0.593296; … ; 0.908361; 0.486164]
      typ: Any
  typ: Any

So, is there an issue trying to optimize on this assuming operators are associative? Does it become very complex or intractable otherwise?

julia
1 Answers

Edit: This whole answer assumes that * is a binary operator in Julia. This is incorrect and as such the OP question is perfectly valid. See https://discourse.julialang.org/t/is-there-any-way-to-make-custom-binary-infix-operators-right-associative/3202/5 for some discussion (specifically mentioning that * is parsed as n-arry, which can be seen here). I will leave the original answer below.

Matrix-matrix multiplication takes n^3 steps (for smallish matrices). Matrix-vector multiplication takes n^2 steps.

a*b*c is the same as (a*b)*c which takes a total of n^3+n^2 operations, while a*(b*c) takes 2*n^2 operations because it contains only matrix-vector multiplications.

If you want something to automatically optimize your matrix algebra expression, you can try a library like https://github.com/Jutho/TensorOperations.jl (the first hit on google for that type of question). Another option mentioned in the comments is https://github.com/AustinPrivett/MatrixChainMultiply.jl

As to why Julia does not do this automatically: That would have to be part of the parser (precedence of operators), however you do not want the parser to be able to peek inside the runtime in order to learn what the object is. On the other hand, this is certainly achievable with a macro applied before compilation as done by the aforementioned packages. Maybe it is also feasible if you use some compile-time type information, however this will be a rather ugly special-case for only a couple of types in the otherwise elegant multiple-dispatch mechanism of Julia.

user_6693290
Related