Concatenate matrix of matrices in Julia

Question

I have a matrix that I am building up block-wise, like so:

parts = Matrix{SparseMatrixCSC{Float64}}(undef, r, r)
for i in 1:r
    for k in 1:r
        parts[i, k] = make_part(i, k)
    end
end

Now I'd like to concatenate all these parts into a larger matrix, with each block in its place. Is there a way to do this using the concatenation functions in Base? I can't figure out how to use cat to do it. For example, reduce((a, acc) -> cat(acc, a, dims=(1, 2), parts) will create a block-diagonal matrix with r^2 blocks on the diagonal.

user_9957710 · Answer

Assuming that you have something like:

parts = [sprand(2,2,0.7) for i in 1:2, j in 1:3]
2×3 Matrix{SparseMatrixCSC{Float64, Int64}}:
…

than using BlockArrays you can do:

julia> mortar(parts)
2×3-blocked 4×6 BlockMatrix{Float64, Matrix{SparseMatrixCSC{Float64, Int64}}, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
  ⋅        0.543703  │   ⋅       0.773536  │  0.390049  0.0287028
 0.113945  0.422259  │   ⋅       0.764168  │  0.140334   ⋅
 ────────────────────┼─────────────────────┼─────────────────────
 0.295535  0.286277  │  0.81173   ⋅        │  0.646139  0.0616408
 0.641121  0.033793  │   ⋅       0.754504  │  0.995975   ⋅

user_1319631 · Answer

I was able to do this by a nesting vcat and hcat like so:

reduce(vcat, [
    reduce(hcat, [parts[i, k] for k in 1:r])
    for i in 1:r
])

However, this doesn't nicely generalize to higher dimensions (although I guess you could nest further loops and use plain cat). It will also allocate for each block row, rather than only allocating once for the entire result.

user_13706735 · Answer

I think the official command for this is hvcat, which is what's used to digest literal matrix syntax. Although whether it's the best solution may depend e.g. on how big theproblem is.

julia> parts = [fill(i + 10j, 2,2) for i in 1:2, j in 1:3]
2×3 Matrix{Matrix{Int64}}:
 [11 11; 11 11]  [21 21; 21 21]  [31 31; 31 31]
 [12 12; 12 12]  [22 22; 22 22]  [32 32; 32 32]

julia> parts[end] = [91 92; 93 94]  # secretly calls hvcat
2×2 Matrix{Int64}:
 91  92
 93  94

julia> hvcat(3, permutedims(parts)...)  # note that it works along the rows!
4×6 Matrix{Int64}:
 11  11  21  21  31  31
 11  11  21  21  31  31
 12  12  22  22  91  92
 12  12  22  22  93  94

julia> parts = [sprand(2,2,0.7) for i in 1:2, j in 1:3];

julia> hvcat(3, permutedims(parts)...) |> summary  # preserves sparseness
"4×6 SparseMatrixCSC{Float64, Int64} with 17 stored entries"

Concatenate matrix of matrices in Julia

3 Answers