Initializing a matrix of matrices

New to Julia
1

Hi everyone,

I have found these posts:

Neither seems to do exactly what I want. I am trying to create a matrix of matrices. I need the structure to be in this way for what I am working on. I will provide an example. What I am trying to do is the following:


for j = 1:2
      M[1,j] = ones(2,2)
end

for j = 1:2
     M[2,j] = [2 2 ; 2 2]
end

I am doing a larger example, but the error message is the same. The error I am getting is MethodError: Cannot ‘convert’ object of type Matrix{Int64} to an object of type Float64.

I know exactly why the error is occurring. When I initialize M with the zeros command, it creates a matrix of floats. Therefore, when I try to replace an entry, it expects floats, but I am trying to give it matrices. To try on a smaller example, I then tried initializing M as follows:

M = [ [0 0 ; 0 0] [0 0 ; 0 0] ; [0 0 ; 0 0] [0 0 ; 0 0]

The above is what I want, a matrix within a matrix. However, when I find M[1,1], for example, it gives me the single entry 0. I want

M[1,1] = [0 0 ; 0 0]

as it is the first component of the big matrix M. Is there a way I can do this?
Thanks!

2

Why do you want a matrix of matrices, when you can simply have one big matrix?

julia> M = [ [0 0 ; 0 0] [0 0 ; 0 0] ; [0 0 ; 0 0] [0 0 ; 0 0] ]
4×4 Matrix{Int64}:
 0  0  0  0
 0  0  0  0
 0  0  0  0
 0  0  0  0

julia> M[1:2, 1:2]
2×2 Matrix{Int64}:
 0  0
 0  0

If you insist, though…

julia> M = [[0 0; 0 0] for i=1:2, j=1:2]
2×2 Matrix{Matrix{Int64}}:
 [0 0; 0 0]  [0 0; 0 0]
 [0 0; 0 0]  [0 0; 0 0]

julia> M[1, 1]
2×2 Matrix{Int64}:
 0  0
 0  0
1 Like
3 
julia> M = [zeros(2,2) for i = 1:2, j=1:2]
2×2 Matrix{Matrix{Float64}}:
 [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]
 [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]

julia> M[1,1]
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0

Note that if you are working with lots of tiny matrices like this you may want to consider StaticArrays.jl, which can be much more efficient for small matrices and vectors of fixed sizes:

julia> using StaticArrays

julia> M = [@SMatrix(zeros(2,2)) for i = 1:2, j=1:2]
2×2 Matrix{SMatrix{2, 2, Float64, 4}}:
 [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]
 [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]

julia> M = fill(@SMatrix(zeros(2,2)), 2,2)
2×2 Matrix{SMatrix{2, 2, Float64, 4}}:
 [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]
 [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]

julia> M = zeros(SMatrix{2,2,Float64,4}, 2,2)
2×2 Matrix{SMatrix{2, 2, Float64, 4}}:
 [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]
 [0.0 0.0; 0.0 0.0]  [0.0 0.0; 0.0 0.0]
2 Likes
4

I would need to provide context to the problem, but I do have a reason for preferring the matrix of matrices. However, it seems I’ve run into an issue I initially didn’t expect. I need to take the inverse of the matrix I’m working with (it’s not all zeros in my case), and the inverse exists globally, but not over each block. Therefore, what I need to do is convert the matrix into the full 4 x 4, take the inverse, and then go back to a matrix of matrices to do the operation I need. Is there a way to do that? Sorry, I didn’t expect this issue initially. Since it’s a different question, I can accept this one and make a new thread if that’s better.

5

I will keep this in mind! At the moment I don’t have a pressing need to do this all the time, but these types of problems I will be solving a lot in the future. I may need this!

Does the package have a way to convert the matrix of matrices into a full matrix (so a 4 x 4 in this case) and then back to the 2 x 2? I found out I need to take the inverse, and the inverse exist globally, but not in each block. Therefore, I need to get the full matrix, and then convert back.

6

Can you just leave everything in a single large matrix, and operate on the submatrices through views?

struct View{A} a::A end  # indexable `View` object
Base.size(v::View) = size(v.a)
Base.getindex(v::View, args...) = view(v.a, args...)

M = [i+j for i=1:4, j=0:4:12]

# construct sub-matrices 
vm = View(M)
A, B, C, D = vm[1:2, 1:2], vm[1:2, 3:4], vm[3:4, 1:2], vm[3:4, 3:4]

# do stuff with sub-matrices
A .+= C
D .*= B
M
8

Sure, this is possible. But this is sounding more and more like an XY problem … what underlying problem are you actually trying to solve?

1 Like
9

In this small example I could, but in general I’m doing this with large matrices. I see you’re defining A,B,C, and D all separately. I don’t want to have to do that 25+ times! Converting would be easiest, but I don’t know if that’s possible

10

I mean I could take multiple paragraphs to explain if you’d like, but I’d rather avoid that. Could you explain the way I’d do it? I think this is the last step in my problem!

11

Julia has lots of looping constructs, so you would not have to do this one at a time. But again, can you give us a bit of background on what problem you are trying to solve?

13

It would take too long. This is nearly the ending step of a multiple step problem. It would take me probably 3-4 paragraphs to explain and show everything. I don’t want to take that much of your time. I believe this is the last step, there’s no more additional functions after this.

14

Like

julia> n = 100
100

julia> M = rand(n,n); # 100x100 matrix

julia> M2 = [@view(M[2i-1:2i,2j-1:2j]) for i=1:n÷2, j=1:n÷2]; # 50x50 matrix of 2x2 blocks

julia> M2[1,1]
2×2 view(::Matrix{Float64}, 1:2, 1:2) with eltype Float64:
 0.663167   0.712037
 0.0853913  0.887108

julia> M[1:2,1:2]
2×2 Matrix{Float64}:
 0.663167   0.712037
 0.0853913  0.887108
1 Like
16

Nearly there! Now I see how to get from the full matrix to the small 2 x 2 ones. The last thing I need is how to get from the small 2 x 2 ones to the full matrix. I need to do that first since I’m storing the 2 x 2 form with my loop. Unless you believe there’s an easier way to refine my earlier loop in my first post to begin with the full matrix? I’m open to that as well! Then I’ll be all set

17

Well, if B is a matrix of matrices, e.g. B = [rand(2,2) for i = 1:50, j=1:50], then you could convert it to a single 100x100 matrix M with e.g.:

M = zeros(100,100)
for i=1:50, j=1:50
    M[2i-1:2i,2j-1:2j] = B[i,j]
end

Another option is to use the BlockArrays.jl package. Then if you have an array of arrays, like B above, you can treat it as a single 100x100 matrix (and do linear algebra directly) with mortar(B), or convert it to an ordinary matrix with Matrix(mortar(B)).

2 Likes
18

I would use @views here.

@views A, B, C, D = M[1:2, 1:2], M[1:2, 3:4], M[3:4, 1:2], M[3:4, 3:4]

produces the same result.

1 Like
19

It can also be done this way:

julia> M = [i+j for i=1:4, j=0:4:12];

julia> B = [M[i:i+1, j:j+1] for i in 1:2:size(M,1), j in 1:2:size(M,2)]
2×2 Matrix{Matrix{Int64}}:
 [1 5; 2 6]  [9 13; 10 14]
 [3 7; 4 8]  [11 15; 12 16]

julia> reduce(hcat, reduce.(vcat, eachcol(B)))
4×4 Matrix{Int64}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

Although I echo the sentiments above that this is probably an XY problem :man_shrugging:

1 Like
20

Yeah, the @views macro is definitely more ergonomic here. I wanted to take this opportunity to toy around more with an indexable view object though :sweat_smile:.

21

This worked perfectly! Thank you so much! I will keep BlockArrays in mind as well.