Creating large matrices from sub blocks

Question

I am trying to generate a matrix from square blocks. Effectively, I have three blocks, call them $A_{diag}$, $A_{up}$, $A_{down}$, and a block of zeroes, which I'll just call $0$. I have some integer $L$ such that if $L=1$ then I am generating the matrix

$$M = \begin{pmatrix} A_{diag} \end{pmatrix}$$

If $L=2$ then I should generate

$$M = \begin{pmatrix} A_{diag} & A_{up} \\ A_{down} & A_{diag} \end{pmatrix}$$

If $L=3$ then

$$M = \begin{pmatrix} A_{diag} & A_{up} & 0 \\ A_{down} & A_{diag} & A_{up} \\ 0 & A_{down} & A_{diag} \end{pmatrix}$$

If $L=4$ then

$$M = \begin{pmatrix} A_{diag} & A_{up} & 0 & 0\\ A_{down} & A_{diag} & A_{up} & 0 \\ 0 & A_{down} & A_{diag} & A_{up} \\ 0 & 0 & A_{down} & A_{diag}\end{pmatrix}$$

and so forth. I suspect this is relatively simple, I found other questions which were similar but more complicated and specific. Any help appreciated, thank you.

Something like SparseArray with Band? (See 1st example in the Basic Examples section.) — corey979, Jan 31 '18 at 18:03

score 4 · Answer 1 · answered Jan 31 '18 at 18:13

You can use ArrayFlatten and DiagonalMatrix:

mat[diag_, up_, down_, n_] := ArrayFlatten @ Block[{Identity},
    DiagonalMatrix[ConstantArray[Identity[diag], n]]+
    DiagonalMatrix[ConstantArray[Identity[up], n-1], 1]+
    DiagonalMatrix[ConstantArray[Identity[down], n-1], -1]
]

I used a blocked Identity head so that DiagonalMatrix didn't complain.

Example:

mat[Array[a, {2,2}], Array[b, {2,2}], Array[c, {2,2}], 3] //TeXForm

$\left( \begin{array}{cccccc} a(1,1) & a(1,2) & b(1,1) & b(1,2) & 0 & 0 \\ a(2,1) & a(2,2) & b(2,1) & b(2,2) & 0 & 0 \\ c(1,1) & c(1,2) & a(1,1) & a(1,2) & b(1,1) & b(1,2) \\ c(2,1) & c(2,2) & a(2,1) & a(2,2) & b(2,1) & b(2,2) \\ 0 & 0 & c(1,1) & c(1,2) & a(1,1) & a(1,2) \\ 0 & 0 & c(2,1) & c(2,2) & a(2,1) & a(2,2) \\ \end{array} \right)$

Thanks, the only downside to this is it doesn't work when n=1 — Kai, Feb 01 '18 at 01:16
I simply added If[n>1, (your code) , Adiag] which works perfectly it seems — Kai, Feb 01 '18 at 01:31

Kai · Accepted Answer · 2018-02-25T21:42:26.623

4

While the other posted methods work, this is the one I wound up using, based on corey979's suggestion:

If[L == 1,
  A[r_] := Adiag[r],
  mup   = SparseArray[Band[{1, 2}] -> 1, {L, L}];
  mdown = SparseArray[Band[{2, 1}] -> 1, {L, L}];
  mdiag = SparseArray[Band[{1, 1}] -> 1, {L, L}];
  A[r_] :=
    KroneckerProduct[mup  ,   Aup[r]] + 
    KroneckerProduct[mdown, Adown[r]] + 
    KroneckerProduct[mdiag, Adiag[r]];
  ];

It has the advantage of using SparseArray's, which is desirable if L is large. One can always use Normal to convert the SparseArray's into regular matrices for small L. For my application for example I found that using normal matrices was faster for L < 3 but using SparseArray's was faster for L > 3, with both being equally fast for L = 3.

edited Feb 25 '18 at 21:42

answered Feb 24 '18 at 04:37

Kai

2,099
8
17

Wow, this one is really fast. Respect. – Henrik Schumacher Mar 11 '18 at 22:46
1

For your mdiag, I would instead use IdentityMatrix[L, SparseArray]. – J. M.'s missing motivation Mar 17 '18 at 11:34
Thanks, I don't notice much change in speed but it certainly makes more sense. – Kai Mar 18 '18 at 05:38

Henrik Schumacher · Answer 3 · 2018-03-11T22:54:39.617

The "SparseArray`" context contains many hidden gems, for example this one:

SparseArray`SparseBlockMatrix[{
   Band[{1, 1}] -> SparseArray@Array[a, {2, 2}],
   Band[{1, 2}] -> SparseArray@Array[b, {2, 2}],
   Band[{2, 1}] -> SparseArray@Array[c, {2, 2}]
   }, {3, 3}] // MatrixForm

$$\left( \begin{array}{cccccc} a(1,1) & a(1,2) & b(1,1) & b(1,2) & 0 & 0 \\ a(2,1) & a(2,2) & b(2,1) & b(2,2) & 0 & 0 \\ c(1,1) & c(1,2) & a(1,1) & a(1,2) & b(1,1) & b(1,2) \\ c(2,1) & c(2,2) & a(2,1) & a(2,2) & b(2,1) & b(2,2) \\ 0 & 0 & c(1,1) & c(1,2) & a(1,1) & a(1,2) \\ 0 & 0 & c(2,1) & c(2,2) & a(2,1) & a(2,2) \\ \end{array} \right)$$

It's not very efficient, though...

score 2 · Answer 4 · answered Jan 31 '18 at 18:14

There it is,

M[n_] := Adiag IdentityMatrix[n] + 
Table[Aup KroneckerDelta[i + 1, j], {i, 1, n}, {j, 1, n}] + 
Table[Adown KroneckerDelta[i - 1, j], {i, 1, n}, {j, 1, n}];

MatrixForm[M[3]]

to elaborate more, you can have Aup[i, j] and Adown[i, j] and Adiag[i, i].

Creating large matrices from sub blocks

4 Answers4

Linked

Related