Best way to subtract Matrix from Vector of Matrix

user_5828964 0
Viewed 78

I want to subtract a Matrix from a Vector of Matrix. For Example, In below defined variable i want to subtract y for each element of xl

xl = [Matrix{Float64}(undef, 2, 3) for i in 1:1000]
y = Matrix{Float64}(undef, 2, 3)

I know this can be done with map function

map(x -> x-y, xl)

But is there any better way to do this or can we reduce runtime by multithreading? I tried benchmarking simple map and ThreadsX.map and found ThreadsX.map is slower any suggestion to reduce runtime.

@btime map(x -> x-y, xl);          # 55.558 μs (1010 allocations: 117.59 KiB)
@btime ThreadsX.map(x -> x-y, xl); # 99.764 μs (5414 allocations: 285.47 KiB)

Note: Here Matrix is just an example I want to do this on large vector of custom struct that support subtraction.

julia
2 Answers

First of all, map is not the fastest thing, loops will almost always beat it. This array comprehension is 70% faster than map.

@btime [[x[i]-$y[i] for i in eachindex($y)] for x in $xl];
  33.300 μs (1001 allocations: 117.31 KiB)

Now, for why multithreading is slower, the answer might be that multithreading setup for this simple calculation has non-negligible overhead. You want a heavier calculation for multithreading to have a noticeable effect. Also, don't forget to interpolate ($) when timing with @btime.

@btime map(x -> x-$y, $xl);
  55.600 μs (1001 allocations: 117.31 KiB)

@btime $xl .- ($y,);     # broadcasting
  56.100 μs (1001 allocations: 117.31 KiB)

@btime ThreadsX.map(x -> x-$y, $xl);
  53.600 μs (1602 allocations: 242.05 KiB)

See this heavier example, you get about 6X speedup with multithreading:

julia> @btime ThreadsX.map(x -> exp.(sin.(x))-exp.(sin.($y)), $xl);
  303.700 μs (3605 allocations: 460.89 KiB)

julia> @btime map(x -> exp.(sin.(x))-exp.(sin.($y)), $xl);
  1.673 ms (3001 allocations: 336.06 KiB)
user_39431700

By default, Julia start with just one thread. To start Julia with more than one thread you can use --threads=n option:

PS C:\Users\Elias> julia --threads=4
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.7.2 (2022-02-06)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> xl = [rand(2, 3) for i in 1:1000];

julia> y = rand(2, 3);

julia> using ThreadsX

julia> using BenchmarkTools

julia> @btime ThreadsX.map(x -> x-$y, $xl);
  42.700 μs (1353 allocations: 206.77 KiB)

For more details, see: https://docs.julialang.org/en/v1/manual/multi-threading/

user_192013700
Related