# Vector of matrices vs. multidimensional arrays

**URL:** <https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602>\
**Category:** Performance\
**Created:** [March 8, 2018, 5:50pm UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602 "2018-03-08T17:50:05Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![plapplop](https://avatars.discourse-cdn.com/v4/letter/p/bbce88/32.png) [@plapplop](https://discourse.julialang.org/u/plapplop)\
**Post date:** [March 8, 2018, 5:50pm UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/1 "2018-03-08T17:50:05Z")

</div>

I have to deal with 3 dimensional structures, I was hesitating between vectors of vectors of vectors, vectors of matrices or tridimensional arrays. For performance comparison, I wrote this small test:

```julia
julia> L = M = N = Int(5e2);

julia> vec_vec_vec = Array{Array{Vector}}(L);

julia> for i = 1:L vec_vec_vec[i] = [zeros(N) for j = 1:M] end;

julia> vec_mat = Array{Matrix}(L);

julia> fill!(vec_mat, rand(M, N));

julia> arr = rand(L, M, N);

julia> tic(); for i = 1:L for j = 1:M for k = 1:N vec_vec_vec[i][j][k] += 1; end; end; end; t = toc(); println(t)
elapsed time: 28.747643437 seconds
28.747643437

julia> tic(); for i = 1:L for j = 1:M for k = 1:N vec_mat[i][j,k] += 1; end; end; end; t = toc(); println(t)
elapsed time: 25.837854348 seconds
25.837854348

julia> tic(); for i = 1:L for j = 1:M for k = 1:N arr[i,j,k] += 1; end; end; end; t = toc(); println(t)
elapsed time: 100.412620007 seconds
100.412620007

```

Can someone explain to me why it is so much more expensive to use a 3d array rather than a vector of matrices, and why the vec\_vec\_vec and vec\_mat solutions are similar? I understand that this example just looks at loops over the whole structures, and even in my case I need to perform matrix-related operations so vec\_vec\_vec is not very suitable. The ideal solution for me would have been to use the 3d array, but this wouldn’t be reasonable with such performances…

Many thanks

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [March 8, 2018, 5:58pm UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/2 "2018-03-08T17:58:53Z")

</div>

One reason your `arr` version is slow is that you’re accessing the indices in exactly the least cache-friendly way possible (this may be different than what you’re used to because Julia/Fortran/MATLAB use a column-major convention while C/NumPy use row major). Try reversing the order of your loops so that the _first_ index changes in the innermost loop. See also [https://docs.julialang.org/en/stable/manual/performance-tips/#Access-arrays-in-memory-order,-along-columns-1](https://docs.julialang.org/en/stable/manual/performance-tips/#Access-arrays-in-memory-order,-along-columns-1)

Edit: this is relevant, but not nearly as important as the reminder not to use non-constant global variables (see below).

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [March 8, 2018, 6:07pm UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/3 "2018-03-08T18:07:41Z")

</div>

Before anything else, read the [Performance Tips](https://docs.julialang.org/en/latest/manual/performance-tips/)–you shouldn’t work with non-constant variables in global scope.

```julia
using BenchmarkTools

# declare global-scope variables as constants, otherwise the compiler can't
# stably infer the variable's type
const L = 500 # integer literals in Julia don't require conversion
const v3 = [[rand(L) for i = 1:L] for j = 1:L]
const vec_mat = Array{Matrix}(L);
const arr = rand(L, L, L);

fill!(vec_mat, rand(L, L));

@btime v3 .+= 1.0
@btime vec_mat .+= 1.0
@btime arr .+= 1.0

```

```julia
  5.767 s (250500 allocations: 993.80 MiB)
  5.411 s (1000 allocations: 953.71 MiB)
  86.659 ms (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [March 8, 2018, 6:08pm UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/4 "2018-03-08T18:08:47Z")

</div>

```julia
function test()
L = M = N = Int(5e2);

vec_vec_vec = Array{Array{Vector}}(L);

for i = 1:L vec_vec_vec[i] = [zeros(N) for j = 1:M] end;

vec_mat = Array{Matrix}(L);

fill!(vec_mat, rand(M, N));

arr = rand(L, M, N);

@time for i = 1:L for j = 1:M for k = 1:N vec_vec_vec[i][j][k] += 1; end; end; end; 

@time for i = 1:L for j = 1:M for k = 1:N vec_mat[i][j,k] += 1; end; end; end; 

@time for k = 1:N for j = 1:M for i = 1:L arr[i,j,k] += 1; end; end; end; 

end
test()

```

results in

```julia
julia> include("test/playground.jl")
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(uninitialized, m)` instead.
│ caller = test() at playground.jl:4
└ @ Main playground.jl:4
┌ Warning: `Array{T}(m::Int) where T` is deprecated, use `Array{T}(uninitialized, m)` instead.
│ caller = test() at playground.jl:8
└ @ Main playground.jl:8
 13.718289 seconds (250.00 M allocations: 3.725 GiB, 3.44% gc time)
 11.875066 seconds (250.00 M allocations: 3.725 GiB, 0.58% gc time)
  0.193893 seconds

```

How’s that?

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [March 8, 2018, 6:10pm UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/5 "2018-03-08T18:10:37Z")

</div>

And another reason is that you are timing your function in the global scope, so _every single access_ to `arr`, `vec_mat`, or `vec_vec_vec` or the sizes L, M, or N is essentially performing a global lookup. These global accesses make it harder to generate efficient code, and the difference can be significant. Performing global lookup inside your innermost loop is essentially the worst possible case.

For example, with no changes at all, here’s your `arr` version just wrapped in a function:

```julia
         L = size(arr, 1)
         M = size(arr, 2)
         N = size(arr, 3)
         for i = 1:L for j = 1:M for k = 1:N arr[i,j,k] += 1; end; end; end
       end
f (generic function with 1 method)

julia> @time f(arr)
  2.066944 seconds (3.50 k allocations: 175.565 KiB)

julia> @time f(arr)
  2.058089 seconds (4 allocations: 160 bytes)

```

Note that the second call is faster because the first call included Julia’s compilation overhead, but both versions are _50 times faster_ than your original code.

We can do even better by accessing the elements in the memory order:

```julia
julia> function f2(arr)
         L = size(arr, 1)
         M = size(arr, 2)
         N = size(arr, 3)
         for k = 1:N for j = 1:M for i=1:L arr[i,j,k] += 1; end; end; end
       end
f2 (generic function with 1 method)

julia> @time f2(arr)
  0.161044 seconds (3.50 k allocations: 175.644 KiB)

julia> @time f2(arr)
  0.144406 seconds (4 allocations: 160 bytes)

```

that’s almost _1000 times faster_ than the original version.

Remember to do these things when profiling is hard, which is why we have [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl). Here’s a much more representative benchmark:

```julia
julia> @benchmark for k = 1:$N for j = 1:$M for i = 1:$L; $arr[i,j,k] += 1; end; end; end
BenchmarkTools.Trial: 
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 115.313 ms (0.00% GC)
  median time: 115.966 ms (0.00% GC)
  mean time: 116.025 ms (0.00% GC)
  maximum time: 118.212 ms (0.00% GC)
  --------------
  samples: 44
  evals/sample: 1

```

Note that this is about the same time we measured from the hand-written `f2` function, but BenchmarkTools actually performs multiple trials to get a statistically meaningful estimate of the function’s performance.

Try using BenchmarkTools to measure the speed of each of your versions. You should find that they’re all quite fast in Julia 😄

---

<div class="post-metadata">

**Author:** ![plapplop](https://avatars.discourse-cdn.com/v4/letter/p/bbce88/32.png) [@plapplop](https://discourse.julialang.org/u/plapplop)\
**Post date:** [March 8, 2018, 6:20pm UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/6 "2018-03-08T18:20:52Z")

</div>

Thank you very much to the three of you, I understood many things. So I should be fine with using 3d arrays, but then if I have to perform matrix operations on such an array, I should ensure that I slice it according to the last two dimensions, i.e. operating on `arr[i,:,:]` should be much more efficient than operating on `arr[:,:,k]`…? (And the same goes for `arr[i,j,:]` and `arr[:,j,k]`.)

I will time the different versions in a proper manner anyway.

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [March 8, 2018, 6:28pm UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/7 "2018-03-08T18:28:38Z")

</div>

Close (but exactly backwards 😉). In a normal Julia array, the values `arr[1, j, k]` and `arr[2, j, k]` are stored next to each other in memory, so it’s most efficient to slice along the _first_ dimension, since all those elements live next to each other. Likewise, slicing along the _first N_ dimensions will be more efficient than the _last N_ dimensions.

By the way, if you find it hard to remember which order to use when looping (I do), then I have good news: for looping over all the elements of an array, there’s an easier, better way. Just do:

```julia
for i in eachindex(arr)
  arr[i] += 1
end

```

`eachindex()` is a Julia function that will produce the “right” kind and order of indices for a particular array-like object of any number of dimensions. For example, for a normal Array it will produce indices in memory order, but it will also do the most efficient thing for a sparse array, some kind of weird row-major array, or lots of other array-like types. For example:

```julia
julia> @benchmark for i in eachindex($arr); $arr[i] += 1; end
BenchmarkTools.Trial: 
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 120.798 ms (0.00% GC)
  median time: 121.303 ms (0.00% GC)
  mean time: 121.723 ms (0.00% GC)
  maximum time: 128.112 ms (0.00% GC)
  --------------
  samples: 42
  evals/sample: 1

```

just as fast as the correctly-ordered nested loops, but with no need to write them out.

---

<div class="post-metadata">

**Author:** ![plapplop](https://avatars.discourse-cdn.com/v4/letter/p/bbce88/32.png) [@plapplop](https://discourse.julialang.org/u/plapplop)\
**Post date:** [March 8, 2018, 6:32pm UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/8 "2018-03-08T18:32:57Z")

</div>

Haha, thanks for your tolerance, I obviously did not think long enough before typing. Great additional information. Thank you for everything.

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [January 12, 2019, 6:17pm UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/9 "2019-01-12T18:17:54Z")

</div>

Sorry to dig up an old thread, but since this was referenced in [another post](https://discourse.julialang.org/t/help-in-understanding-efficiently-handling-multidimensional-arrays/19512/11) I figured it was worth correcting the record.

The entire difference you’re observing between the vector-of-matrix approach and the 3D array approach here is due to your vector having a non-concrete element type. Your `vec_mat` is of type `Array{Matrix}`, but `Matrix` is not a concrete type. You need `Array{Matrix{Float64}}`. Likewise for your vector-of-vectors-of-vectors approach, you need an `Array{Vector{Vector{Float64}}}` not `Array{Vector{Vector}}`.

Switching to concretely typed structures makes the 3D array and vector-of-matrix approaches perform almost exactly the same:

```julia
julia> function f(x::Array{T, 3}) where T
         x .+= 1
       end
f (generic function with 2 methods)

julia> function f(x::Array{Matrix{T}}) where T
         for matrix in x
           matrix .+= 1
         end
       end
f (generic function with 2 methods)

julia> x_array = zeros(N, N, N);

julia> x_vec_mat = [zeros(N, N) for i in 1:N];

julia> using BenchmarkTools

julia> @btime f($x_array);
  609.838 μs (0 allocations: 0 bytes)

julia> @btime f($x_vec_mat);
  505.012 μs (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [January 13, 2019, 2:46am UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/10 "2019-01-13T02:46:45Z")

</div>

True:

```julia
function test()
	L = M = N = Int(5e2);
	vec_vec_vec = Array{Vector{Vector{Float64}}}(undef, L);
	for i = 1:L vec_vec_vec[i] = [zeros(N) for j = 1:M] end;
	vec_mat = Array{Matrix{Float64}}(undef, L);
	fill!(vec_mat, rand(M, N));
	arr = rand(L, M, N);
	@time for i = 1:L for j = 1:M for k = 1:N vec_vec_vec[i][j][k] += 1; end; end; end; 
	@time for i = 1:L for j = 1:M for k = 1:N vec_mat[i][j,k] += 1; end; end; end; 
	@time for k = 1:N for j = 1:M for i = 1:L arr[i,j,k] += 1; end; end; end; 
end
test()

```

gives

```julia
julia> include("playground.jl");                                                        
  0.203685 seconds                                                                      
  0.224881 seconds                                                                      
  0.147060 seconds    

```

with

```julia
julia> versioninfo()                                                                    
Julia Version 1.2.0-DEV.158                                                             
Commit 7772486a69 (2019-01-12 19:21 UTC)                                                
Platform Info:                                                                          
  OS: Windows (x86_64-w64-mingw32)                                                      
  CPU: Intel(R) Core(TM) i7-6650U CPU @ 2.20GHz                                         
  WORD_SIZE: 64                                                                         
  LIBM: libopenlibm                                                                     
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)      

```

---

<div class="post-metadata">

**Author:** ![spaceLem](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/spacelem/32/217628_2.png) [@spaceLem](https://discourse.julialang.org/u/spaceLem)\
**Post date:** [January 17, 2019, 8:35am UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/11 "2019-01-17T08:35:55Z")

</div>

Sorry to add just one more small thing to this, but instead of  
`for i=1:L for j=1:M for k=1:N … end; end; end`  
you can write  
`for i=1:L, j=1:M, k=1:N … end`  
which is shorter and easier to read.

---

<div class="post-metadata">

**Author:** ![bennedich](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bennedich/32/4894_2.png) [@bennedich](https://discourse.julialang.org/u/bennedich)\
**Post date:** [January 17, 2019, 9:58am UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/12 "2019-01-17T09:58:27Z")

</div>

You can rearrange the indexing to optimal memory order for the vector-of-matrices version to speed it up. Also, `@inbounds` helps a lot here:

```julia
function test()
	L = M = N = Int(5e2);
	vec_vec_vec = Array{Vector{Vector{Float64}}}(undef, L);
	for i = 1:L vec_vec_vec[i] = [zeros(N) for j = 1:M] end;
	vec_mat = Array{Matrix{Float64}}(undef, L);
	fill!(vec_mat, rand(M, N));
	arr = rand(L, M, N);
	@time @inbounds for i = 1:L, j = 1:M, k = 1:N vec_vec_vec[i][j][k] += 1 end
	@time @inbounds for i = 1:L, k = 1:N, j = 1:M vec_mat[i][j,k] += 1 end
	@time @inbounds for k = 1:N, j = 1:M, i = 1:L arr[i,j,k] += 1 end
end
test()

```

Gives me:

```julia
  0.124387 seconds
  0.059382 seconds
  0.097464 seconds

```

There’s quite a difference still in performance between the three versions IMO (the fastest being twice as fast as the slowest), I haven’t looked into what’s causing that. I’d guess that either some versions are being vectorized (SIMD) more efficiently, or it has to do with memory accesses and caching.

* * *

**Edit:** BenchmarkTools gives similar timings:

```julia
julia> @btime @inbounds for i = 1:$L, j = 1:$M, k = 1:$N; $vec_vec_vec[i][j][k] += 1; end
  111.676 ms (0 allocations: 0 bytes)

julia> @btime @inbounds for i = 1:$L, k = 1:$N, j = 1:$M; $vec_mat[i][j,k] += 1; end
  58.472 ms (0 allocations: 0 bytes)

julia> @btime @inbounds for k = 1:$N, j = 1:$M, i = 1:$L; $arr[i,j,k] += 1; end
  76.587 ms (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [January 17, 2019, 12:00pm UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/13 "2019-01-17T12:00:27Z")

</div>

> [@bennedich](#):
>
> ```julia
> fill!(vec_mat, rand(M, N))
> 
> ```

Will give you the same matrix in all the slots for `vec_mat`. Probably explains the timings. There is no reason why the last one should not be the fastest.

---

<div class="post-metadata">

**Author:** ![bennedich](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bennedich/32/4894_2.png) [@bennedich](https://discourse.julialang.org/u/bennedich)\
**Post date:** [January 17, 2019, 12:05pm UTC](https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/14 "2019-01-17T12:05:24Z")

</div>

Oh, good point, I missed that! Indeed, fixing that, the last one is now fastest:

```julia
  113.182 ms (0 allocations: 0 bytes)
  94.457 ms (0 allocations: 0 bytes)
  79.743 ms (0 allocations: 0 bytes)

```