CODEWITHSUNDEEP
numpymatrixscipy
drerD
drerD

Reputation: 689

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

enter image description here

A is NxN, B is NxM. They are the matrices for a dynamical system x=Ax+Bu. I could propagate the matrix using np.block(), but I hope there's a way of forming this matrix can scale based on the size of 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, thanks.

EDIT: I can do the following to generate a hardcoded solution.

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]
         ] )

input:

A: [[0.35227098 0.98853832 0.42468451]
 [0.49288462 0.34457356 0.79954007]
 [0.52171255 0.63167711 0.27997043]]
B: [[0.53736681]
 [0.05086935]
 [0.42558792]]

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 with 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]])

Upvotes: 0

Views: 320

Answers (1)

a_guest
a_guest

Reputation: 36289

You can use itertools.accumulate to compute the block matrices which reoccur throughout the array:

import itertools as it

AB = list(it.accumulate(it.chain([B], it.repeat(A, N-1)), lambda x, y: y @ x))
zero = np.zeros_like(B)
result = np.block([AB[i::-1] + [zero]*(N-1-i) for i in range(N)])

Upvotes: 3

Related Questions