# How to dynamically define matrix with blocks in subdiagonals?

Clear[V0, V1, V];
Nm = 10;
\[Sigma] = {({
{0, 1.},
{1., 0}
}), ({
{0, -I},
{I, 0}
}), ({
{1., 0},
{0, -1.}
})};(*Pauli-spin matrices*)
V0 = Normal[
SparseArray[{Band[{2, 1}] -> Table[V[n], {n, 1, Nm - 1}],
Band[{1, 2}] -> Table[ConjugateTranspose[V[n]], {n, 1, Nm - 1}]},
Nm]];
V[n_] := IdentityMatrix[2] + I ({n, n, n}. \[Sigma]);
V1 = ArrayFlatten[V0];


Above is the code with fix size. In this way I can use arrayFlatten to define the normal matrix. See example below:

However, the thing will change if I use a variable input (See below).

Clear[V0, V1, V]
\[Sigma] = {({
{0, 1.},
{1., 0}
}), ({
{0, -I},
{I, 0}
}), ({
{1., 0},
{0, -1.}
})};(*Pauli-spin matrices*)
V0[Nm_] :=
Normal[SparseArray[{Band[{2, 1}] -> Table[V[n], {n, 1, Nm - 1}],
Band[{1, 2}] -> Table[ConjugateTranspose[V[n]], {n, 1, Nm - 1}]},
Nm]];
V[n_] := IdentityMatrix[2] +  I ({n, n, n}. \[Sigma]);
V1[Nm_] := ArrayFlatten[V0[Nm]];


If I evaluate above, then I got error as follows:

How to solve this problem?

• What is \[Sigma]? Please post complete code. Jul 20 at 18:41
• @HenrikSchumacher Thanks for pointing out. Sigma are pauli matrices as shown in the code now. Jul 20 at 20:37
• What is α? If it is not necessary to show the problem, it is better to remove it, but if it is needed, can you, please, provide such detail of its definition? Thanks! Jul 20 at 21:45
• @CATrevillian Sorry, α doesn't matter so I removed them. Jul 20 at 22:06

Here are two versions of the same solution. First, using a general 2x2 matrix along the subdiagonal

ClearAll[v]
v[Nm_] := Block[{array, band},
array = With[{mat =
Normal@SparseArray[Band[{2, 1}] -> Array[band, Nm - 1], Nm]},
ArrayFlatten[mat /. {0 -> ConstantArray[0, {2, 2}], band[n_] :>
n{{b11, b12}, {b21, b22}}}]];
(array + ConjugateTranspose[array])
]

v[4] /. Conjugate[b_] :> Superscript[b, "*"] // MatrixForm


The code creates a variable array with the 2x2 matrices along the subdiagonal. The code takes advantage of the Band notation in SparseArray. band is an indexed variable that stands for 2x2 matrices. band variables are replaced by actual 2x2 matrices. Zeroes in the matrix are also replaced by 2x2 matrices. Then the complex transpose is added to produce the superdiagonal.

Here is a version that uses the Pauli matrices. This should reproduce the results of an earlier version of the question. This version uses ComplexExpand, which assumes that $$\alpha$$ is a real number.

ClearAll[v]
v[Nm_] := Block[{array, subdiag, b},
subdiag = Array[b, Nm - 1];
array = Normal@SparseArray[Band[{2, 1}] -> subdiag, Nm];
array = array /. {0 -> ConstantArray[0, {2, 2}]};
array = array /. b[n_] :> Cos[α] IdentityMatrix[2] +
I n Sin[α] Total[Array[PauliMatrix, 3]];
array = ArrayFlatten[array];
(array + ConjugateTranspose[array]) // ComplexExpand]

v[4] // MatrixForm


• Thanks for the solution. Do you think the ArrayFlatten should automatically do all the things above to make it more intuitive? Jul 21 at 0:50
• Thanks for the comment. I am editing the second version of the function to make each step clearer. Jul 21 at 2:32

I don't know why the first case does not through an error, but SparseArray expectes the correct dimensions of the array as second argument. So V0 should read as follows:

V0[Nm_] := SparseArray[{
Band[{3, 1}] -> Table[V[n], {n, 1, Nm - 1}],
Band[{1, 3}] -> Table[ConjugateTranspose[V[n]], {n, 1, Nm - 1}]
},
{2 Nm, 2 Nm}
];


The result is a matrix of dimensions{2 Nm, 2 Nm} and you do not have to apply ArrayFlatten afterwards.

• This is not correct since what I want to do is really arrayFlatten the Block Matrices V[n] (each of them is 2x2 matrix). By using arrayFlatten, I can convert the Block matrices into the normal matrix. See the example in my post. Jul 20 at 23:34
• See this post mathematica.stackexchange.com/questions/761/… Jul 20 at 23:39
• Please see updated code. The point is to use Band[{3, 1}] instead of Band[{2, 1}]. Jul 20 at 23:39
• Ok, now this works. Thank! But from my point of view, the indices should not depend on the size of the block and should be automatically handled by Mathematica. Jul 20 at 23:49
• Yes, I agree. This is somewhat counterintuitive. Jul 20 at 23:50