 # Turn vector of matrices into 3D array

Background: My simulation program returns a 4x4 matrix. Running many replicates using `map()` or a vector comprehension, I get a long vector of small matrices. Similar to `[rand(4, 4) for _ in 1:10000]`.

Question: Is there a one-line way to turn a vector of 10,000 small (4x4) matrices into a 3D array (size 4x4x1000)? If the vector was short, one could do: `cat(vector_of_vectors..., dims=3)`. But that fails for longer vectors as pointed out by others before.

The analogous problem of turning a vector of vectors into a matrix is easy: `reduce(hcat, vector_of_vectors)`.

For my matrices, I tried `reduce(cat, vector_of_matrices)` and `reduce(x -> cat(x, dims=3), vector_of_matrices)` without success. At the moment, I use this function, which works and is fast (`mv` being the vector of matrices):

``````function matVec2Arr(mv)
a = zeros(size(mv)..., length(mv))
for i in 1:length(mv)
a[:,:,i] = mv[i]
end
return a
end
``````

I am just wondering is there a more ‘elegant’ way?

Thanks,
Hannes

If the matrices are always 4x4 (thus small) use static arrays and store a vector of these, the memory layout will be the same of a 3D array and a “reinterpret” should do the trick.

1 Like
``````julia> @benchmark reduce(hcat, b) setup=b=[rand(4, 4) for _ in 1:10000]
BenchmarkTools.Trial: 3898 samples with 1 evaluation.
Range (min … max):  158.440 μs …   5.706 ms  ┊ GC (min … max):  0.00% … 93.17%
Time  (median):     242.868 μs               ┊ GC (median):     0.00%
Time  (mean ± σ):   327.088 μs ± 232.752 μs  ┊ GC (mean ± σ):  16.53% ± 18.96%

▁▇█▅▅▃▃▃▃▂▁▁ ▁                            ▁                 ▁
▅▅███████████████▇▆▆▆▅▃▂▂▂▂▃▄▄▅▃▂▄▂▄▅▆█▄▆▆▆██▆▇▆▆▅▆▅▆▃▄▅▄▄▅▄▅ █
158 μs        Histogram: log(frequency) by time       1.22 ms <

Memory estimate: 1.22 MiB, allocs estimate: 2.
``````

this is fine (non-stipid), see the allocation

`TensorCast` makes such operations easy one-liners:

``````using TensorCast
m = [rand(4, 4) for _ in 1:10_000]
@cast M[i,j,k] := m[k][i,j]   # := returns a view
``````

It returns `M` which is `4x4x10_000`

NB:
I wish this was the default language.

3 Likes

Consider RecursiveArrayTools:

``````julia> using RecursiveArrayTools

julia> A = [rand(2,3) for _=1:4]
4-element Vector{Matrix{Float64}}:
[0.35766183647954497 0.00945510515311243 0.18756680509300394; 0.9738104612662519 0.10918176238450261 0.7844365367234503]
[0.8572574873275078 0.606234830814707 0.8911587693644292; 0.29714178204328334 0.34785399506946235 0.4215022402770998]
[0.31381806075696206 0.09938494049668578 0.33507384871214696; 0.8384919930451804 0.9773745902770696 0.5043700409936012]
[0.7990079554473333 0.39911052428069627 0.49464323002688526; 0.057411978409010134 0.39387907129167 0.28244091355499146]

julia> B = VectorOfArray(A)
VectorOfArray{Float64,3}:
4-element Vector{Matrix{Float64}}:
[0.35766183647954497 0.00945510515311243 0.18756680509300394; 0.9738104612662519 0.10918176238450261 0.7844365367234503]
[0.8572574873275078 0.606234830814707 0.8911587693644292; 0.29714178204328334 0.34785399506946235 0.4215022402770998]
[0.31381806075696206 0.09938494049668578 0.33507384871214696; 0.8384919930451804 0.9773745902770696 0.5043700409936012]
[0.7990079554473333 0.39911052428069627 0.49464323002688526; 0.057411978409010134 0.39387907129167 0.28244091355499146]

julia> B[1,2,3]
0.09938494049668578
``````

There is also ArraysOfArrays.jl offering `VectorOfArrays` type, but I have not explored it.

I think ArraysOfArrays solves the opposite problem, i.e. it’s like `eachslice` not like `cat`. There is a small zoo of packages for doing `cat`-like things, my summary of all the ones I found is here. RecursiveArrayTools has the slightly weird feature that linear indexing doesn’t do what you’d expect for a 3-array: `B isa Matrix`.

You can use that here, with a reshape:

``````julia> A = [rand(4, 4) for _ in 1:10000];

julia> @time B = reduce(hcat, A); summary(B)  # done efficiently
0.000543 seconds (2 allocations: 1.221 MiB)
"4×40000 Matrix{Float64}"

julia> C = reshape(B, 4, 4, :); summary(C)  # almost free
"4×4×10000 Array{Float64, 3}"

julia> C[1,2,3] == A[1,2]
true

julia> @time C == reduce((x,y) -> cat(x, y, dims=3), A)  # done pairwise, slow!
0.193491 seconds (262.96 k allocations: 396.856 MiB, 9.68% gc time, 28.20% compilation time)
true
``````

Although this reshape doesn’t work beyond matrices, and it would be convenient to have a specific function. In julia#21672 it sounds like someone just has to get around to writing it.

4 Likes

This is what I was looking for.
Great to know that

``````B = reduce(hcat, A);
C = reshape(B, 4, 4, :);
``````

is more efficient.