Speed up resampling from array

user_2063292 0
Viewed 175

I'm doing a lot of resampling from large arrays, by row, and I was wondering if there's a way to speed it up. For example

n = 10^4
a = rand(Float64, (n,n))
@time r = a[sample(1:n,n),:]

I usually get about 0.8 seconds on my machine. sample() itself is quite fast. Indeed, r = a[1:n, :] is just about as slow as the above. I know these are large arrays, but am I missing something obvious? An order of magnitude speedup would be wonderful...

EDIT: I have selected Przemyslaw Szufel's answer as it is nice and comprehensive, and definitely much faster if you're not going to further manipulate the array. Unfortunately in my case, Fredrik Bagge's caution proved true: it was not overall faster to use views, because subsequent operations to the array become slower--it was basically a wash in my testing. Also, Oscar Smith made a good point about column major. In my case, I later do something across the other dimension that's even more expensive the resampling so it made sense to leave the array as it is.

julia
2 Answers

Use views for a 2000x speedup!

julia> const n = 10^4; const a = rand(Float64, (n,n));

julia> using BenchmarkTools

julia> @btime a[sample(1:n,n),:];
  336.622 ms (6 allocations: 763.02 MiB)

julia> @btime a[:, sample(1:n,n)];
  230.512 ms (6 allocations: 763.02 MiB)

julia> @btime view(a,sample(1:n,n),:);
  164.601 μs (5 allocations: 78.31 KiB)

julia> @btime view(a,:, sample(1:n,n));
  165.601 μs (5 allocations: 78.31 KiB)

Note that when creating the view column or row selection does not matter. It will matter however when the data is going to be read from the view.

EDIT

@Fredrik Bagge made a very important comment that using later the view will be much slower. While naturally, when copying data happens there is cost to be incurred, there are the following issues to remember:

  1. In practice very often not all elements of a view might be used in data processing - in all those cases one gets an immediate speedup.

  2. The memory footprint of a view approach will be much lower.

  3. One can always materialize the view. Of course if you materialize the view you need to incur the cost of copying the data in memory. However it will not be slower than materializing directly. Let's benchmark:

julia> @btime a[sample(1:n,n),:];
  351.572 ms (6 allocations: 763.02 MiB)

julia> const myview = view(a, sample(1:n,n), :);

julia> @btime collect(myview);
  297.866 ms (2 allocations: 762.94 MiB)

These are results from my machine. Actually, for row-oriented querying of a matrix, creating a view and materializing it later seems to be faster than materializing the Array up-front.

Let's have a look at the column oriented matrix:

julia> @btime a[:, sample(1:n,n)];
  255.276 ms (6 allocations: 763.02 MiB)

julia> const myview2 = view(a, :, sample(1:n,n));

julia> @btime collect(myview2);
  297.952 ms (2 allocations: 762.94 MiB)

At least on my machine materializing the view has the same speed regardless of using column or row major. So for the example presented in the question is it actually faster, while for the optimized column-oriented code it is slower. However note again point (1.) - not always all data is subsequently used.

EDIT 2: Additional explanation following the second comment by Frederik. Benchmarking the collect code we have skipped the call to sample when measuring the times. However, note that the time to run sample is orders of magnitude smaller than the variance of subsequent runs of @btime it does not affect the results.

user_99577100

Slicing into an array always creates a copy of the data, so this can get quite expensive for large slices. In your case, it would probably make sense to use a view, which doesn't create a copy, but references the underlying data in a instead. You can create a view either with the @view macro as r = @view a[sample(1:n, n), :] or by using the view function directly: r = view(a, sample(1:n, n), :). Hope that helps!

user_32631640
Related