Sparse block matrices all wrong

Question

I need to represent $M-lambdacdottextrm{Identity}$ where $M$ is an operator on $V_1opluscdotsoplus V_n$ made from
$$
V_1xrightarrow{M_1}V_2xrightarrow{M_2}cdotsxrightarrow{M_{n-1}}V_nxrightarrow{M_n}V_1.
$$
With, say, $n=3$, $dim V_k=k+1$, I am trying

sa =
 With[{dim = {2, 3, 4}},   
  With[{n = Length[dim]},    
   SparseArray[
    {
     Band[{1, 1}] -> Table[λ IdentityMatrix[dim[[k]]], {k, n}],
     Band[{2, 1}] -> 
      Table[
       Table[Indexed[M[k], {i, j}], {i, dim[[k + 1]]}, {j, dim[[k]]}], 
       {k, n - 1}
      ],
     Band[{1, n}] -> Table[Indexed[M[n], {i, j}], {i, dim[[1]]}, {j, dim[[n]]}]
    }
   ]
  ]
 ]

The result is very strange: both sa//TableForm and ArrayFlatten[sa]//TableForm give

The same happens with Band[{1, n}] -> {Table[...]}.

What am I doing wrong here?

Henrik Schumacher · Accepted Answer

SparseArray reacts different when lists of matrices appear on the right hand side of Rule for Band. In this case, it assembles the ArrayFlattened block matrix instead of the matrix of blocks and Band[{i,j}] refers to the positions in the final, assembled matrix, not to the position within the matrix of blocks. I have to admit that this is really counterintuitive.
sa[dim_] := With[{n = Length[dim]},
  SparseArray[{
    Band[{1, 1}] -> λ,
    Band[{dim[[1]] + 1, 1}] -> 
     Table[Table[ Indexed[M[k], {i, j}], {i, dim[[k + 1]]}, {j, dim[[k]]}], {k, n - 1}], 
    Band[{1, Total[Most[dim]] + 1}] -> 
     Table[Indexed[M[n], {i, j}], {i, dim[[1]]}, {j, dim[[n]]}]},
   {Total[dim], Total[dim]}
   ]
  ]
sa[{2, 3, 4}]

There is also the undocumented function SparseArray`SparseBlockMatrix in which you can use {i,j}-> ... for the block in the i-row and j-th column. But SparseArray`SparseBlockMatrix does not go well with Band. And it is undocumented, so one has to go mostly by trial and error.

kglr · Answer

An alternative way without using Band:

sA = Module[{dim = #, n = Length @ #, l = Array[[FormalX], Length @ #], 
     lim = λ IdentityMatrix[Total @ #], mats},
 mats = Table[Indexed[M[k], {i, j}], {k, n}, {i, dim[[Mod[k + 1, n, 1]]]}, {j, dim[[k]]}];
 SparseArray[ArrayFlatten[RotateRight[DiagonalMatrix[l]] /. Thread[l -> mats]] + lim]] &;

Examples:

MatrixForm @ sA @ {2, 3, 4}

MatrixForm @ sA @ {2, 3, 4, 3}

Sparse block matrices all wrong

2 Answers

Add your own answers!

Ask a Question