How to dynamically define matrix with blocks in subdiagonals?

Question

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?

@HenrikSchumacher Thanks for pointing out. Sigma are pauli matrices as shown in the code now. — Bob Lin, Jul 20 '21 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! — CA Trevillian, Jul 20 '21 at 21:45

Henrik Schumacher · Answer 1 · 2021-07-20T23:37:55.040

2

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.

edited Jul 20 '21 at 23:37

answered Jul 20 '21 at 23:23

Henrik Schumacher

106,770
7
179
309

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. – Bob Lin Jul 20 '21 at 23:34
See this post https://mathematica.stackexchange.com/questions/761/how-to-enter-matrices-in-block-matrix-format – Bob Lin Jul 20 '21 at 23:39
1

Please see updated code. The point is to use Band[{3, 1}] instead of Band[{2, 1}]. – Henrik Schumacher Jul 20 '21 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. – Bob Lin Jul 20 '21 at 23:49
1

Yes, I agree. This is somewhat counterintuitive. – Henrik Schumacher Jul 20 '21 at 23:50
Also, the size of the sparse matrix is depending on the size and number of Block matrices, which should be calculated accordingly. – Bob Lin Jul 21 '21 at 00:04

LouisB · Accepted Answer · 2021-07-21T02:36:58.657

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? — Bob Lin, Jul 21 '21 at 00:50
Thanks for the comment. I am editing the second version of the function to make each step clearer. — LouisB, Jul 21 '21 at 02:32

How to dynamically define matrix with blocks in subdiagonals?

2 Answers2