计算多输入多输出数据序列汉克尔矩阵的矢量化方法

Question

我正在构造一个Hankel矩阵，并想知道是否有一种方法可以进一步矢量化以下计算（即，没有for循环或列表理解）

# Imagine this is some time-series data
q = 2  # Number of inputs
p = 2  # Number of outputs
nt = 6  # Number of timesteps
y = np.array(range(p*q*nt)).reshape([nt, p, q]).transpose()
assert y.shape == (q, p, nt)
print(y.shape)

(2, 2, 6)

print(y[:,:,0])

[[0 2]
 [1 3]]

print(y[:,:,1])

[[4 6]
 [5 7]]

print(y[:,:,2])

[[ 8 10]
 [ 9 11]]

预期结果

# Desired Hankel matrix
m = n = 3  # dimensions
assert nt >= m + n
h = np.zeros((q*m, p*n), dtype=int)
for i in range(m):
    for j in range(n):
        h[q*i:q*(i+1), p*j:p*(j+1)] = y[:, :, i+j]
print(h)

[[ 0  2  4  6  8 10]
 [ 1  3  5  7  9 11]
 [ 4  6  8 10 12 14]
 [ 5  7  9 11 13 15]
 [ 8 10 12 14 16 18]
 [ 9 11 13 15 17 19]]

（注意2x2块的堆叠方式）

# Alternative method using stacking
b = np.hstack([y[:,:,i] for i in range(y.shape[2])])
assert b.shape == (q, p*nt)
print(b)

[[ 0  2  4  6  8 10 12 14 16 18 20 22]
 [ 1  3  5  7  9 11 13 15 17 19 21 23]]

h = np.vstack([b[:, i*p:i*p+n*q] for i in range(m)])
print(h)

[[ 0  2  4  6  8 10]
 [ 1  3  5  7  9 11]
 [ 4  6  8 10 12 14]
 [ 5  7  9 11 13 15]
 [ 8 10 12 14 16 18]
 [ 9 11 13 15 17 19]]

Answer 1

您可以使用stride_tricks：

>>> from numpy.lib.stride_tricks import as_strided
>>> 
>>> a = np.arange(20).reshape(5,2,2)
>>> s0,s1,s2 = a.strides
>>> as_strided(a,(3,2,3,2),(s0,s2,s0,s1)).reshape(6,6)
array([[ 0,  2,  4,  6,  8, 10],
       [ 1,  3,  5,  7,  9, 11],
       [ 4,  6,  8, 10, 12, 14],
       [ 5,  7,  9, 11, 13, 15],
       [ 8, 10, 12, 14, 16, 18],
       [ 9, 11, 13, 15, 17, 19]])