Turn vector of matrices into 3D array

New to Julia
, , , ,
1

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[1])..., 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

2

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
3 
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

4

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
5

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.

6

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[1] 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[3][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.

5 Likes
7

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

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

is more efficient.

1 Like