Matrix block indexing

General Usage
1

I also posted this question on stackeroverflow.com, so if anyone wants to also answer it there I would happily accept it.

I am using Julia version 1.1. I am working a lot with matrices that can be constructed from smaller matrices, e.g., the Pauli matrices. It is not clear to me how to efficiently construct big matrices by using a set of smaller matrices in Julia, i.e., directly write the smaller matrix into a certain index position.

Julias kron is not satisfactory since I would need to generate several “big matrices” to get my final result. E.g. I would like to create something like this (this is only a very small example):

sy = [[0 -im]; [im 0]]
M = [[0 sy adjoint(sy)]; 
     [adjoint(sy) 0 sy];
     [sy adjoint(sy) 0]]

It would be possible to do so by doing two kronecker products adding the two results. However this would be a huge waste, especially if the matrices got bigger.

I also already tried to work with the package BlockArrays.jl but realised that it does not fullfill my need, since it only allows to divide the matrix into 2x2 blocks of arbitrary size. However, I need an arbitrary amount of blocks depending on the size of my submatrix and “big” matrix.

In the end I want to be able to address “matrix blocks” of my big matrix so that I can directly assign the construction matrices to the right position, for the above example this would look like the following (I did not use a loop here to make my point clearer):

M[1, 2] = sy
M[1, 3] = adjoint(sy)
M[2, 1] = adjoint(sy)
M[2, 3] = sy
M[3, 1] = sy
M[3, 2] = adjoint(sy)

I realise that this means reducing my original big array indices to something like array “block indices”.

I thought about doing this with views, where I create a matrix of SubArrays that I can then address with the matrix block index notation e.g.

S0 = view(M, 1:2, 1:2)
S1 = view(M, 1:2, 2:4)
S2 = view(M, 1:2, 4:6)
...

Viewmatrix = [[S0 S1 S2]; [S3 S4 S5]; [S6 S7 S8]]
Viewmatrix[1, 2] .= sy
Viewmatrix[1, 3] .= adjoint(sy)
...

Now it is unclear to me how one would actually go about this and write such a view matrix in general or if this is even a feasable way to address the problem. If there is a better way to approach this problem I would like to know it.

2

Given sy = [0 -im; im 0], your example can be done easily with:

julia> M = [0I sy sy'; sy' 0I sy; sy sy' 0I]
6×6 Array{Complex{Int64},2}:
 0+0im  0+0im  0+0im  0-1im  0+0im  0-1im
 0+0im  0+0im  0+1im  0+0im  0+1im  0+0im
 0+0im  0-1im  0+0im  0+0im  0+0im  0-1im
 0+1im  0+0im  0+0im  0+0im  0+1im  0+0im
 0+0im  0-1im  0+0im  0-1im  0+0im  0+0im
 0+1im  0+0im  0+1im  0+0im  0+0im  0+0im

Note that 0I automatically expands to a zero block of the right size. You can also do e.g.

M[1:2, 3:4] = sy

to write into a block. That might be the easiest thing to do if you are constructing M programmatically.

However this would be a huge waste, especially if the matrices got bigger.

How big are your matrices and how often do you need to construct them? I would do the simplest thing first, and then only optimize if it is a problem.

(If the matrices are really big, of course, then you probably want to exploit sparsity or other structure.)

3

How big are your matrices and how often do you need to construct them? I would do the simplest thing first, and then only optimize if it is a problem.

This very much depends on the exact physical implementation at hand. I would say the matrices I work with usually have sizes of around (1000x1000) to (10000x10000) entries. You are right in the respect that I overexaggerated my wording with “huge” waste. Mostly I create one reference matrix where I can exchange the entries blockwise with other matrices.

In the end it is much more about; how to efficiently and intuitively index matrix blocks. I realise that your suggestion M[1:2, 3:4] = sy works perfectly well, but would lead to for loops with awkward ranges, if one would like to assign many matrices blockwise.

After I also got an answer on stackoverflow, I realise that my use case is covered by BlockArrays.jl. The documentation was just a bit sparse with that information. With that package I can use PseudoBlockArrays to directly achieve what I want and still use “simple” indexing.