Toeplitz矩阵的toeplitz矩阵

Question

我要创建toeplitz矩阵的toeplitz矩阵。H1、H2和H3已经是toeplitz矩阵了。我的结果应该是这样的：H1 0 0 H2 H1 0 H3 H2 H1 0 H3 H2 0 0 H3

现有的toeplitz函数只接受向量，所以我不能用它作为矩阵。目前，我使用vstack创建第一列，然后创建第二列等等，然后使用hstack合并所有列。这需要很大的努力，因为我必须在某些地方专门添加np.zeros矩阵。我想不出更好的方法来连接numpy数组，因为只有少数几个函数，它们都不适合我的问题。

Answer 1

与对vstack和hstack的嵌套调用不同，更有效的方法是预先分配最终数组，然后使用嵌套循环填充数组。您最初可以使用高维数组来保持代码的整洁。

例如，这个脚本

import numpy as np

H1 = np.array([[11, 11], [11, 11]])
H2 = np.array([[22, 22], [22, 22]])
H3 = np.array([[33, 33], [33, 33]])

inputs = (H1, H2, H3)

# This assumes all the arrays in `inputs` have the same shape,
# and that the data type of all the arrays is the same as H1.dtype.
nh = len(inputs)
nrows = 2*nh - 1
m, n = H1.shape
# T is a 4D array.  For a given i and j, T[i, :, j, :] is a 2D array
# with shape (m, n).  T can be intepreted as a 2D array of 2D arrays. 
T = np.zeros((nrows, m, nh, n), dtype=H1.dtype)
for i, H in enumerate(inputs):
    for j in range(nh):
        T[i + j, :, j, :] = H

# Partially flatten the 4D array to a 2D array that has the desired
# block structure.
T.shape = (nrows*m, nh*n)

print(T)

版画

[[11 11  0  0  0  0]
 [11 11  0  0  0  0]
 [22 22 11 11  0  0]
 [22 22 11 11  0  0]
 [33 33 22 22 11 11]
 [33 33 22 22 11 11]
 [ 0  0 33 33 22 22]
 [ 0  0 33 33 22 22]
 [ 0  0  0  0 33 33]
 [ 0  0  0  0 33 33]]

(注意，结果不是Toeplitz矩阵；它是块Toeplitz矩阵。)

Answer 2

对于任何对这个问题感兴趣的人来说，这里有一个替代的方法

from pylab import *
import scipy.linalg

H1 = array([[11, 11], [11, 11]])
H2 = array([[22, 22], [22, 22]])
H3 = array([[33, 33], [33, 33]])

# Setup blocks 
t = array([zeros_like(H1), H1, H2, H3])
# Create index to each block, using toeplitz
idx = scipy.linalg.toeplitz(r_[1, 2, 3, zeros(2)], r_[1, zeros(2)]).astype(int)
# Index into blocks, transpose and reshape to get re-ordered array
#  copy is used to ensure memory is nicely ordered
T = t[idx].transpose(0, 2, 1, 3).reshape(10, 6).copy()

大部分时间都花在scipy.linalg.toeplitz上，这比在这里使用的小矩阵在数组中填充内存要慢，所以在使用这种方法之前，我建议进行分析。

Answer 3

下面是一个使用Kronecker产品的示例，它为所需的行和列数(nr, nc)提供了一定的灵活性：

import numpy as np

H1 = np.array([[11, 11], [11, 11]])
H2 = np.array([[22, 22], [22, 22]])
H3 = np.array([[33, 33], [33, 33]])

inputs = (H1, H2, H3)
nr = 5
nc = 3
m, n = H1.shape

T = np.zeros([m*nr,n*nc])
for ii, H in enumerate(inputs):
    T = T + np.kron(np.eye(nr,nc,-ii),H)

print(T)

版画

[[11. 11.  0.  0.  0.  0.]
 [11. 11.  0.  0.  0.  0.]
 [22. 22. 11. 11.  0.  0.]
 [22. 22. 11. 11.  0.  0.]
 [33. 33. 22. 22. 11. 11.]
 [33. 33. 22. 22. 11. 11.]
 [ 0.  0. 33. 33. 22. 22.]
 [ 0.  0. 33. 33. 22. 22.]
 [ 0.  0.  0.  0. 33. 33.]
 [ 0.  0.  0.  0. 33. 33.]]