Creating a block matrix from arrays of blocks

Question

I am trying to generate a matrix from square blocks. Effectively, I have a $n×n$ matrix polynomial $P(l)$, the $qth$ derivative of $P(l)$ with respect to $l$, which is denoted by $P^{(q)} (l)$, and a block of zeroes, which I’ll just call $0$. I have some integer $k$ such that if $k=1$ then I am generating the matrix

$$ R= \begin{pmatrix} P(l) \end{pmatrix} $$

If $k=2$ then I should generate

$$ R = \begin{pmatrix} P(l) & 0 \\ \frac{1}{1!} P^{(1)}(l) & P(l) \end{pmatrix} $$

If $k=3$ then

$$ R = \begin{pmatrix} P(l) & 0 & 0 \\ \frac{1}{1!} P^{(1)}(l) & P(l) & 0 \\ \frac{1}{2!} P^{(2)}(l) & \frac{1}{1!} P^{(1)}(l) & P(l) \end{pmatrix} $$

and so forth. Generally,

$$ R = \begin{pmatrix} P(l) & 0 & \cdots & 0 & 0 \\ \frac{1}{1!} P^{(1)}(l) & P(l) & \cdots & 0 & 0 \\ \frac{1}{2!} P^{(2)}(l) & \frac{1}{1!} P^{(1)}(l) & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \frac{1}{(k-1)!} P^{(k-1)}(l) & \frac{1}{(k-2)!} P^{(k-2)}(l) & \cdots & \frac{1}{1!} P^{(1)}(l) & P(l) \end{pmatrix} $$

is an $nk×nk$ matrix.

I prefer a simple and understandable way and for that I thought to start with a zero matrix $R$ of dimensions $nk×nk$ and then with two "for" loops to full the initial zero matrix, putting the corresponding derivative which is needed. I’m not sure in what should go as my statement in "for" loops. I found other questions which were similar but more complicated and specific. Any help appreciated, thank you.

Answer 1

Because of the banded structure of your matrix, you may use the undocumented function SparseArray`SparseBlockMatrix for that as follows (using a dummy matrix-valued function P here):

P[i_, l_, n_] := ConstantArray[l, {n, n}];
R[k_, l_, n_] := SparseArray`SparseBlockMatrix[
  Table[Band[{i, 1}] -> 1/i! P[i, l, n], {i, 1, k}], 
  {k, k}
  ]

E.g., R[4, ell, 3] produces a matrix of size {12, 12} with block entries specified by the function P.

Answer 2

First we need same data. I denote your "P(l)" bl p[l] and `P^(1)" by der[n,p,x]. Then I define a function "makeblmat" that assembles the block matrix:

m = Table[RandomInteger[{-10, 10}], 2, 2];
p[x_] = Sum[MatrixPower[m, i] x^i, {i, 0, 2}];
derp[n_, p_, x_] := D[p, {x, n}];
makeblmat[k_] := 
 Table[If[i >= j, derp[i - j, p[x], x]/! (i - j), 0], {i, k}, {j, k}]//ArrayFlatten;
MatrixForm[makeblmat[2]]