CODEWITHSUNDEEP
arraysmatrixmultidimensional-arrayjulia
Invarianz
Invarianz

Reputation: 461

Matrix block indexing

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.

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.

Upvotes: 2

Views: 731

Answers (1)

phipsgabler
phipsgabler

Reputation: 20950

BlockArrays.jl does not only support 2×2 blocked arrays, although they used to use only them in their documentation. You can easily create create the 3×3-blocked 6×6 array you want as follows:

M = BlockArray(fill(0im, 6, 6), [2, 2, 2], [2, 2, 2])
M[Block(1, 2)] = sy
M[Block(1, 3)] = adjoint(sy)
M[Block(2, 1)] = adjoint(sy)
M[Block(2, 3)] = sy
M[Block(3, 1)] = sy
M[Block(3, 2)] = adjoint(sy)

julia> M
3×3-blocked 6×6 BlockArray{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

But be careful: the blocks are stored by reference. So if you modify sy afterwards, all blocks containing it will be changed as well, and vice versa. If you want to avoid that, use broadcast assignment (.= instead of =).

If your problem is actually as simple as the example, and more on the dense side, it might be simpler to use the mortar function to "stick together" the available blocks:

julia> mortar(reshape([z, sy, sy', sy', z, sy, sy, sy', z], (3, 3)))
3×3-blocked 6×6 BlockArray{Complex{Int64},2,Array{AbstractArray{Complex{Int64},2},2},BlockArrays.BlockSizes{2,Array{Int64,1}}}:
 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

Although that uses an abstract type internally, instead of promoting the assigned arrays.

Upvotes: 2

Related Questions