# Matrix block indexing

**URL:** <https://discourse.julialang.org/t/matrix-block-indexing/22922>\
**Category:** General Usage\
**Tags:** arrays\
**Created:** [April 8, 2019, 2:06pm UTC](https://discourse.julialang.org/t/matrix-block-indexing/22922 "2019-04-08T14:06:51Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Invarianz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/invarianz/32/7825_2.png) [@Invarianz](https://discourse.julialang.org/u/Invarianz)\
**Post date:** [April 8, 2019, 2:06pm UTC](https://discourse.julialang.org/t/matrix-block-indexing/22922/1 "2019-04-08T14:06:51Z")

</div>

I also posted this [question](https://stackoverflow.com/questions/55562897/matrix-block-indexing) on [stackeroverflow.com](http://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):

```julia
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):

```julia
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.

```julia
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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 8, 2019, 3:13pm UTC](https://discourse.julialang.org/t/matrix-block-indexing/22922/2 "2019-04-08T15:13:01Z")

</div>

> [@Invarianz](#):
>
> E.g. I would like to create something like this (this is only a very small example):

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

```julia
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.

```julia
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.)

---

<div class="post-metadata">

**Author:** ![Invarianz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/invarianz/32/7825_2.png) [@Invarianz](https://discourse.julialang.org/u/Invarianz)\
**Post date:** [April 8, 2019, 3:40pm UTC](https://discourse.julialang.org/t/matrix-block-indexing/22922/3 "2019-04-08T15:40:53Z")

</div>

> 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.