# Is there an efficient way to form this block matrix with numpy or scipy?

Is there an efficient way to form this block matrix with numpy or scipy?

$$\left[ \begin{array}{cccc} \mathbf{B} & \mathbf{0} & \cdots & \mathbf{0}\\ \mathbf{AB} & \mathbf{B} & \cdots & \mathbf{0}\\ \vdots & \vdots & \ddots & \vdots\\ \mathbf{A}^{N-1}\mathbf{B} & \mathbf{A}^{N-2}\mathbf{B} & \cdots & \mathbf{B}\\ \end{array} \right]$$

Here, $$\mathbf{A}$$ is $$N\times N$$, $$\mathbf{B}$$ is $$N\times M$$. They are the matrices for a dynamical system $$x=\mathbf{A}x+\mathbf{B}u$$.

I could propagate the matrix using np.block(), but I hope there's a way of forming this matrix that can scale based on $$N$$. I was thinking maybe Kronecker product np.kron() can help, but I can't think of a way. Any advice would help, and if anyone knows what this matrix is called please let me know.

I can do the following to generate a hardcoded solution for $$N=3$$.

N = 3
A = np.random.random((N, N))
B = np.random.random((N, 1))
print("A:", A)
print("B:", B)
np.block([ [B, np.zeros((B.shape[0], B.shape[1]*(N-1)))],
[A@B, B, np.zeros((B.shape[0], B.shape[1]*(N-1-1)))],
[np.linalg.matrix_power(A,N-1)@B, np.linalg.matrix_power(A,N-1-1)@B, B]
] )


output:

array([[0.53736681, 0.        , 0.        ],
[0.05086935, 0.        , 0.        ],
[0.42558792, 0.        , 0.        ],
[0.42032564, 0.53736681, 0.        ],
[0.62266266, 0.05086935, 0.        ],
[0.43163605, 0.42558792, 0.        ],
[0.94690357, 0.42032564, 0.53736681],
[0.76683544, 0.62266266, 0.05086935],
[0.73345624, 0.43163605, 0.42558792]])


What I am hoping to achieve is an efficient way to generate such a matrix that scales with $$N$$, and I am able to do it with np.block() and list but it just seems to be not efficient for me.

mat_list = []
for i in range(N): # generate row
tmp = []
for j in range(N): # loop through A^j*B
if j <= i:
tmp.append(np.linalg.matrix_power(A,i-j)@B)
if i < N-1:
tmp.append(np.zeros((B.shape[0], B.shape[1]*(N-1-i))))
mat_list.append(tmp)

np.block(mat_list)


output:

array([[0.53736681, 0.        , 0.        ],
[0.05086935, 0.        , 0.        ],
[0.42558792, 0.        , 0.        ],
[0.42032564, 0.53736681, 0.        ],
[0.62266266, 0.05086935, 0.        ],
[0.43163605, 0.42558792, 0.        ],
[0.94690357, 0.42032564, 0.53736681],
[0.76683544, 0.62266266, 0.05086935],
[0.73345624, 0.43163605, 0.42558792]])

• Two questions: (1) By scales with N do you mean $O(N)$ in big-O notation? If so, that's impossible since you want to construct a full matrix of size $N^2 \times NM$. Simply touching every entry is $O(N^3)$. (2) Do you really need to explicitly construct this matrix? If you just need to, for example, perform a matrix-vector-product there might be ways to do that more efficiently. – LedHead Aug 15 at 4:19
• @LedHead (1) no I want a function that takes in N to construct the matrix. (2) that'll be what I prefer, and that was why I mentioned Kronecker product which might be able to do it, but I haven't figured how. – drerD Aug 15 at 5:14

The code proposed by the OP can indeed made be more efficient, mainly by noting the fact that to form the sequence $$A^i B$$, with $$i=0\,\dots,N$$ you do not have to compute $$A^i$$ at each step, but you can exploit the fact that $$A^i B = A\,(A^{i-1}B)$$, reusing the result of the previous step.

My proposed implementation is

import numpy as np

N = 3
M = 1

A = np.random.random((N, N))
B = np.random.random((N, M))

X = np.zeros((N ** 2, N * M))
rsl = slice(0, N)
X[rsl, :M] = B
for i in range(1, N):
rsl_p, rsl = rsl, slice(i * N, (i + 1) * N)
X[rsl, :M] = A @ X[rsl_p, :M]
X[rsl, M : (i + 1) * M] = X[rsl_p, : i * M]


This question is more related to programming that computational science, so I add some general remarks, useful for a novice programmer.

1. never blindly translate mathematical expressions into code: think algorithmically.

X = np.linalg.matrix_power(A,i)@B


is much slower than

X=B;
for j in range(i): X = A*X


if $$M \ll N$$

(exercise: compute operation count for both code fragments).

2. whenever possible reuse previous calculations.

3. pre-allocate output matrices, and do not store intermediate results in temp. matrices

• thanks, based on my benchmark this is about 10 times as fast as what I had. I am guessing this is already optimal solution for constructing such as matrix. – drerD Aug 15 at 5:15
• @drerD glad that this was useful to you. The proposed implementation should also have better asymptotics with respect to the old one, so I expect that the speed-up should increase with $N$. – Stefano M Aug 15 at 21:09

If you need to explicitly construct the entire matrix, then Stefano M's answer is your best bet. If, however, you don't really need the whole matrix, but just need to be able to perform a matrix-vector-product (MVP), then you might want to consider the following approach. Based on your comments it's unclear to me if this helps you, and I realize this doesn't technically answer your question, but I thought I'd post it just in case.

Suppose you want to calculate the MVP:

$$\left( \begin{array}{c}y_1 \\ y_2 \\ \vdots \\ y_N \end{array} \right) = \left( \begin{array}{cccc} B & 0 & ... & 0 \\ A B & B & ... & 0 \\ \vdots & \vdots & \ddots & \vdots \\ A^{N-1} B & A^{N-2} B & ... & B \end{array} \right) \left( \begin{array}{c}x_1 \\ x_2 \\ \vdots \\ x_N \end{array} \right)$$

where $$x_i$$'s are size $$Mx1$$ and known and we want to calculate $$y_i$$'s which are size $$Nx1$$ and unknown. This can be expressed block-by-block as

$$\begin{array}{lll} y_i & = & \sum\limits_{j=1}^{i} A^{i-j} B x_j\\ & = & B x_i + \sum\limits_{j=1}^{i-1} A^{i-j} B x_j \\ & = & B x_i + A \sum\limits_{j=1}^{i-1} A^{(i-1)-j} B x_j \\ & = & B x_i + A y_{i-1}\\ \end{array}$$

for $$i=1,2,...,N$$. (The last step works if we assume $$y_0 \equiv \mathbf{0}^{[N \times 1]}$$.)

This should be much faster than using the full dense matrix.

Example Python code:

import numpy as np

def mvp(A,B,x):

n=A.shape[0]
m=B.shape[1]

# assume that x is 1D array of size [n*m]

y=np.zeros( (n*n) )

y[:n] = B @ x[:m]
for i in range(1,n):
y[i*n:(i+1)*n] = B @ x[i*m:(i+1)*m] + A @ y[(i-1)*n:i*n]

return y

N = 3
M = 1

A = np.random.random((N, N))
B = np.random.random((N, M))
x=np.random.random( (N*M) )

y=mvp(A,B,x)


• thanks, that's exactly how the original dynamical equations get turned into the matrix. When you mentioned MVP, I was thinking you might be thinking of doing the following, for example N=3, M=1 case: – drerD Aug 20 at 18:25
• $\left( \begin{array}{cccc} B & 0 & 0 \\ A B & B & 0 \\ A^2 B & AB & B \end{array} \right) = \left( \left(\begin{array}{c}0 \\ 0 \\ A \end{array} \right) \left(\begin{array}{c}A \ I \ 0 \end{array} \right) + \left(\begin{array}{c}0 \\ A \\ 0 \end{array} \right) \left(\begin{array}{c}I \ 0 \ 0 \end{array} \right) + \left( \begin{array}{cccc} I & 0 & 0 \\ 0 & I & 0 \\ 0 & 0 & I \end{array} \right) \right) \left( \begin{array}{cccc} B & 0 & 0 \\ 0 & B & 0 \\ 0 & 0 & B \end{array} \right)$ – drerD Aug 20 at 18:30
• This is something I just derived, but I think there might be efficient way of utilizing matrix vector products to create such a matrix, especially without for loop. – drerD Aug 20 at 18:34
• One thing I just tried is scipy.linalg.toeplitz([range(N)], np.hstack([[1], np.zeros(N-1) ])), which produces a matrix of the powers for A, then if I can just map it, then multiply by the diagonal matrix of Bs, it'll be quite efficient. – drerD Aug 20 at 18:39
• A question I probably should have asked earlier: how large of $N$ and $M$ are you typically dealing with? – LedHead Aug 20 at 18:39