mapslices performance

[gideonsimpson]

Consider the following piece of code:

n=10^6;
X = randn(1,n);
V(x) = 0.5 * dot(x,x);
f(V, X) = [ V((@view X[:,n])) for n = 1:size(X,2) ]

Timing on this reveals:

julia> @time mapslices(V, X, 1);
  1.177690 seconds (12.00 M allocations: 267.024 MiB, 6.91% gc time)
julia> @time f(V,X);
  0.047745 seconds (1.00 M allocations: 53.406 MiB, 5.58% gc time)

which is substantially different.

Note that in this simple example, my data could be a column vector; however, this is a surrogate for problems where I have time series data and X is d x n with 1< d << n.

[KristofferC]

Better on 0.7, but not parity:

julia> @time mapslices(V, X, 1);
  0.202644 seconds (5.00 M allocations: 114.428 MiB, 5.10% gc time)

julia> @time f(V,X);
  0.032617 seconds (1.00 M allocations: 53.406 MiB, 16.69% gc time)

Edit:

Both idx and ridx at:

julia/base/abstractarray.jl

Lines 1940 to 1958 in 8f46450

	@noinline function inner_mapslices!(safe_for_reuse, indices, nidx, idx, otherdims, ridx, Aslice, A, f, R)
	if safe_for_reuse
	# when f returns an array, R[ridx...] = f(Aslice) line copies elements,
	# so we can reuse Aslice
	for I in indices # skip the first element, we already handled it
	replace_tuples!(nidx, idx, ridx, otherdims, I)
	_unsafe_getindex!(Aslice, A, idx...)
	R[ridx...] = f(Aslice)
	end
	else
	# we can't guarantee safety (#18524), so allocate new storage for each slice
	for I in indices
	replace_tuples!(nidx, idx, ridx, otherdims, I)
	R[ridx...] = f(A[idx...])
	end
	end
	return R
	end

are Vector{Any} and seems to be used in quite tight loops.

[zsoerenm]

JuliennedArrays might help in this case:

julia> @benchmark mapslices(V, X, 1)
BenchmarkTools.Trial: 
  memory estimate:  267.02 MiB
  allocs estimate:  11999527
  --------------
  minimum time:     750.878 ms (2.31% GC)
  median time:      782.824 ms (2.33% GC)
  mean time:        803.105 ms (3.38% GC)
  maximum time:     905.073 ms (2.32% GC)
  --------------
  samples:          7
  evals/sample:     1

julia> @benchmark f(V, X)
BenchmarkTools.Trial: 
  memory estimate:  53.41 MiB
  allocs estimate:  1000006
  --------------
  minimum time:     23.018 ms (5.64% GC)
  median time:      25.374 ms (5.97% GC)
  mean time:        37.715 ms (32.41% GC)
  maximum time:     101.140 ms (67.24% GC)
  --------------
  samples:          133
  evals/sample:     1

julia> using JuliennedArrays
julia> @benchmark map(V, julienne(X, (:,*)))
BenchmarkTools.Trial: 
  memory estimate:  53.41 MiB
  allocs estimate:  1000006
  --------------
  minimum time:     26.012 ms (5.17% GC)
  median time:      28.425 ms (5.71% GC)
  mean time:        41.413 ms (30.04% GC)
  maximum time:     109.614 ms (65.14% GC)
  --------------
  samples:          121
  evals/sample:     1

There are plans to migrate it into Base: #23645 (comment)

[stillyslalom]

Updated timings on latest nightly:

julia> @btime mapslices($V, $X, dims=1);
  377.917 ms (10998504 allocations: 251.75 MiB)

julia> @btime f($V,$X);
  14.712 ms (2 allocations: 7.63 MiB)

julia> @btime map(V, Slices($X, 1));
  14.832 ms (6 allocations: 7.63 MiB)

mapslices is still struggling, but eachslice has improved ~5x since v1.5.2, and is now tied with hand-written & julienned versions:

julia> @btime map(V, eachslice($X, dims=2));
  76.877 ms (4000009 allocations: 129.70 MiB) # v1.5.2

julia> @btime map(V, eachslice($X, dims=2));
  14.937 ms (7 allocations: 7.63 MiB) # v1.6.0-DEV.1208

[stillyslalom]

...and just for fun,

using LoopVectorization, Tullio
g(X) = @tullio y[j] := 0.5*X[i, j]^2

julia> @btime g($X)
  754.101 μs (172 allocations: 7.64 MiB)

[mcabbott]

Timing with #40996:

julia> @btime mapslices($V, $X, dims=1);
  251.758 ms (10998504 allocations: 251.75 MiB) # master
  7.804 ms (22 allocations: 7.63 MiB)           # with PR 40996

julia> @btime f($V,$X);
  5.087 ms (2 allocations: 7.63 MiB)

julia> @btime map(V, eachslice($X, dims=2));
  5.082 ms (7 allocations: 7.63 MiB)

This example is one for which you only need half of mapslices, since f returns a scalar. The second half deals with re-assembling slices returned by f to make another big array.

[mbauman]

Yes, this looks to be close enough to be called fixed now:

julia> using BenchmarkTools, LinearAlgebra

julia> n=10^6;

julia> X = randn(1,n);

julia> V(x) = 0.5 * dot(x,x);

julia> f(V, X) = [ V((@view X[:,n])) for n = 1:size(X,2) ]
f (generic function with 1 method)

julia> @benchmark mapslices($V, $X, dims=1)
BenchmarkTools.Trial: 655 samples with 1 evaluation per sample.
 Range (min … max):  6.955 ms …  12.696 ms  ┊ GC (min … max): 0.00% … 40.19%
 Time  (median):     7.444 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   7.628 ms ± 627.234 μs  ┊ GC (mean ± σ):  5.16% ±  6.22%

   ▁▁██▆▅▂       ▁▂
  ▄███████▆▄▄▅▅▆▆███▇▆▄▄▄▄▄▂▃▂▃▃▂▄▃▄▄▄▄▃▄▄▃▂▂▃▃▃▃▂▁▂▂▂▁▁▁▁▁▁▂ ▄
  6.96 ms         Histogram: frequency by time        9.64 ms <

 Memory estimate: 7.66 MiB, allocs estimate: 18.

julia> @benchmark f($V, $X)
BenchmarkTools.Trial: 949 samples with 1 evaluation per sample.
 Range (min … max):  4.609 ms …   9.122 ms  ┊ GC (min … max): 0.00% … 44.62%
 Time  (median):     5.106 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   5.264 ms ± 608.372 μs  ┊ GC (mean ± σ):  7.77% ±  8.73%

   ▁▄▆██▄          ▂
  ▄███████▇▆▆▃▄▇▅█▇██▇▆▆▅▆▄▄▃▃▄▃▃▂▃▂▃▃▂▁▃▃▃▄▃▃▄▄▄▄▃▄▃▃▂▂▃▂▃▂▃ ▄
  4.61 ms         Histogram: frequency by time        6.96 ms <

julia> @benchmark map($V, eachslice($X, dims=2))
BenchmarkTools.Trial: 985 samples with 1 evaluation per sample.
 Range (min … max):  4.321 ms …   9.597 ms  ┊ GC (min … max): 0.00% … 40.72%
 Time  (median):     4.910 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   5.067 ms ± 693.443 μs  ┊ GC (mean ± σ):  8.66% ±  9.73%

    ▂█▅▂▃        ▂
  ▅▇█████▆▆▆▃▃▄▅▇█▇▇▆▄▄▄▃▃▃▂▂▂▂▂▂▂▁▂▂▂▃▂▂▃▃▃▂▃▂▃▂▃▂▂▂▂▂▂▁▁▂▁▁ ▃
  4.32 ms         Histogram: frequency by time        6.91 ms <

 Memory estimate: 7.66 MiB, allocs estimate: 3.