Why does b *= A[r] allocate memory?

General Usage Performance
,
1

The following code does a whole bunch of allocations:

N = 10000000
A = randn(Float64, N)
b = 1.0
@show typeof(b)
@time for r = 1:N
    b *= A[r]
 end

Why is that? How can I write this in a way that doesn’t allocate ~ a gigabyte of memory? (In this particular case the right answer is prod, but this is a MWE we’ve run into in a more complicated context, where it’s not completely clear how to use prod.)

In particular:

julia> N = 10000000; A = randn(Float64, N); b = 1.0;@show typeof(b); @time for r = 1:N; b = b*A[r] end
typeof(b) = Float64
 0.636152 seconds (50.00 M allocations: 915.513 MiB, 3.44% gc time, 0.94% compilation time)

julia> N = 10000000; A = randn(Float64, N); b = 1.0;@show typeof(b); @time for r = 1:N; b = b*A[r] end
typeof(b) = Float64
 0.612025 seconds (50.00 M allocations: 915.512 MiB, 2.22% gc time)
2

don’t benchmark in the global scope. You’re seeing type instability since A is a non-type asserted global variable.

julia> N = 1000;
julia> A = randn(Float64, N);
julia> b = 1.0;
julia> @time for r = 1:N
           b *= A[r]
        end
  0.007243 seconds (4.08 k allocations: 82.250 KiB, 98.37% compilation time)
julia> function myprod(A)
           b = 1.0
           for r = 1:length(A)
               b *= A[r]
           end
           b
       end
julia> @time myprod(A);
  0.006203 seconds (3.66 k allocations: 203.109 KiB, 99.70% compilation time)

julia> @time myprod(A);
  0.000006 seconds (1 allocation: 16 bytes)
8 Likes
3

Hmm. Yes, right, this is a very good point about my MWE. Our actual problem isn’t in global scope (I don’t think), so it’s not completely obvious how that’s going to play out, but I suppose there could be some similar type-instability effect.

Probably we should be wrapping the inside of the inner loop in a function anyway, which will help.

In any case thank you!

2 Likes
4

Hello,

I ran the same codes as yours but didn’t get your result in the last line(i.e. 1 allocation: 16 bytes). Could you tell me why?

julia> N = 100000000;

julia> A = randn(Float16,N)
100000000-element Vector{Float16}:
0.1344
-1.894
0.5503
-0.282
0.829
-1.03
-0.1345
0.0896
0.1259
2.053
⋮
-0.173
2.73
0.1725
0.2515
-0.5327
-0.8926
0.195
-0.3242
1.868

julia> function myprod(A)
b = 1.0
for r = 1:N
b *=A[r]
end
b
end
myprod (generic function with 1 method)

julia> @time myprod(A);
5.137363 seconds (500.01 M allocations: 8.941 GiB, 2.59% gc time, 0.19% compilation time)

julia> @time myprod(A);
5.201536 seconds (500.00 M allocations: 8.941 GiB, 2.53% gc time)

Thanks!

5

That would because I had a typo. Specifically, before, I had 1:N which is slow because N is a global variable instead of 1:length(A)

6

I see the difference! Thank you very much

7

The 16 bytes are also a benchmarking artefact, btw. It should be zero, which is what you get when you interpolate the input argument:

julia> @btime myprod($A);
  142.350 ms (0 allocations: 0 bytes)
5 Likes