Slicing a matrix of matrices

General Usage
,
1

I would like to better understand how to slice matrices of matrices.

For example, I initialise a matrix of matrices this way:

M = [[1 0; 0 2] for i = 1:2, j=1:2]

If I slice it this way: M[1,1][:,:] I obtain the result [1 0; 0 2] , which is what I would expect. It correspond to the entry (1,1) of M.

However of I slice is that way M[:,:][1,1] I obtain a similar result. I would expect to obtain the result [1 1; 1 1] because each matrix of M has 1 in position (1,1).

Is it possible to obtain this results using slicing? Or does one need to iterate explicitly over M to extract this information ?

thanks in advance & best regards !

2

What’s going on here is simply that M[:, :] == M, so that M[:, :][1, 1] == M[1, 1] == M[1, 1][:, :].

I doubt you can get [M[i, j][1, 1] for i = axes(M, 1), j = axes(M, 2)] via a direct slice, because this operation just does not really make sense for a Matrix{Matrix}. Note that there is nothing stopping us from the inner matrices having different dimensions.

julia> M[1, 1] = rand(1:10, 3, 4);

julia> M
2×2 Matrix{Matrix{Int64}}:
 [7 2 7 4; 1 4 10 5; 8 4 10 10]  [1 0; 0 2]
 [1 0; 0 2]                      [1 0; 0 2]

Of course, if you stack everything into a fourth order tensor, you can slice just fine:

julia> M = [rand(2, 2) for i = 1:2, j = 1:2]
2×2 Matrix{Matrix{Float64}}:
 [0.19189 0.794161; 0.0113205 0.159683]  [0.838033 0.61838; 0.733269 0.217573]
 [0.468851 0.748258; 0.563997 0.146174]  [0.154156 0.509213; 0.913721 0.309242]

julia> S = stack(M)  # Note that M[i, j] == S[:, :, i, j]
2×2×2×2 Array{Float64, 4}:
[:, :, 1, 1] =
 0.19189    0.794161
 0.0113205  0.159683

[:, :, 2, 1] =
 0.468851  0.748258
 0.563997  0.146174

[:, :, 1, 2] =
 0.838033  0.61838
 0.733269  0.217573

[:, :, 2, 2] =
 0.154156  0.509213
 0.913721  0.309242

julia> S[1, 1, :, :]
2×2 Matrix{Float64}:
 0.19189   0.838033
 0.468851  0.154156

julia> [M[i, j][1, 1] for i = axes(M, 1), j = axes(M, 2)]
2×2 Matrix{Float64}:
 0.19189   0.838033
 0.468851  0.154156
3 Likes
3

I think you can do:

getindex.(M, 1, 1)
5 Likes
4

This really sounds like you should just use a 4D array. That supports then all the slicing operations you can think of.

5

For a matrix, slicing with [:, :] isn’t very meaningful, since it just creates a copy of the matrix, which is not what you seem to want.

Instead of M[:, :], it’s better to be explicit and write copy(M), since that is easier to understand.

6

thanks a lot for the suggestions!
In practice, I’m working with small (2x2 or 3x3) collections of small (2x2 or 3x3) static matrices.
The getindex.() works smoothly. I will also try the rank 4 tensor version using a MArray but this require some more coding.

here’s a MWE:

using StaticArrays, LinearAlgebra

let 
    # Overall data structure - a potentially large size (nx x ny) matrix os SMatrices 
    nx, ny = 11, 11
    D = [@MMatrix(zeros(4,4)) for _ in 1:nx, _ in 1:ny]

    # Fill information into the matrices contained in D (not done here)
    for 𝐼 in CartesianIndices(D)
        D[𝐼] .= I(4)
    end

    # Loop over the data set access information
    for 𝐼 in CartesianIndices(D)

        i, j = 𝐼[1], 𝐼[2]

        # Avoid boundaries
        if i>1 && i<nx && j>1 && j<ny

            ########### CONSTRUCT LOCAL OBJECTS ###########

            # Collect information in 3x3 matrices 
            D_loc_v1 = SMatrix{3,3}( D[ii,jj] for ii in i-1:i+1, jj in j-1:j+1)

            # Collect information in a tensor of size {3, 3, 4, 4} 
            D_loc_v2 = @MArray(zeros(3,3,4,4))
            for ii=1:3, jj=1:3, kk=1:4, ll=1:4
                D_loc_v2[ii,jj,kk,ll] = D[i-2+ii,jj-1+jj][kk,ll] 
            end

            ########### ACCESS LOCAL OBJECTS ###########
 
            # Version 1
            D_loc_11_v1 =  getindex.(D_loc_v1, 1, 1)
            D_loc_12_v1 =  getindex.(D_loc_v1, 1, 2)

            # Version 2
            D_loc_11_v2 = D_loc_v2[:,:,1,1]
            D_loc_12_v2 = D_loc_v2[:,:,1,2]

        end
    end
end