在 numpy 中构造块 Hankel 矩阵的矢量化方法(或 scipy)
Vectorized way to construct a block Hankel matrix in numpy (or scipy)
我想构造以下矩阵:
[v0 v1 v2 v3 .... v(M-d+1)
v1 .
v2 . .
. .
.
vd . . v(M) ]
其中每个 v(k) 是一个 (ndarray) 向量,例如来自矩阵
X = np.random.randn(100, 8)
M = 7
d = 3
v0 = X[:, 0]
v1 = X[:, 1]
...
使用 for 循环,我可以这样做,例如:
v1 = np.array([1, 2, 3]).reshape((-1, 1))
v2 = np.array([10, 20, 30]).reshape((-1, 1))
v3 = np.array([100, 200, 300]).reshape((-1, 1))
v4 = np.array([100.1, 200.1, 300.1]).reshape((-1, 1))
v5 = np.array([1.1, 2.2, 3.3]).reshape((-1, 1))
X = np.hstack((v1, v2, v3, v4, v5))
d = 2
X_ = np.zeros((d * X.shape[0], X.shape[1]+1-d))
for i in range (d):
X_[i*X.shape[0]:(i+1) * X.shape[0], :] = X[:X.shape[0], i:i+(X.shape[1]+1-d)]
我得到:
X_ = array([[ 1. , 10. , 100. , 100.1],
[ 2. , 20. , 200. , 200.1],
[ 3. , 30. , 300. , 300.1],
[ 10. , 100. , 100.1, 1.1],
[ 20. , 200. , 200.1, 2.2],
[ 30. , 300. , 300.1, 3.3]]) #Which is the wanted matrix
有没有什么方法可以用矢量化的方式构造这个矩阵(我想当涉及到大矩阵时,这会比 for 循环更快?)。
谢谢。
这看起来是最佳的;你已经很好地矢量化了它。我能做的唯一改进是用 np.empty
替换 np.zeros
,这会跳过初始化数组。我尝试使用 np.vstack
和 np.lib.stride_tricks.sliding_window_view
(在 之后)并获得了与使用 np.empty
.
的 for 循环相同的性能
# sliding window:
X_ = np.lib.stride_tricks.sliding_window_view(X, (X.shape[0], X.shape[1]+1-d)).reshape(d*X.shape[0], -1)
# np.vstack:
X_ = np.vstack([X[:, i:i+(X.shape[1]+1-d)] for i in range(d)])
我想构造以下矩阵:
[v0 v1 v2 v3 .... v(M-d+1)
v1 .
v2 . .
. .
.
vd . . v(M) ]
其中每个 v(k) 是一个 (ndarray) 向量,例如来自矩阵
X = np.random.randn(100, 8)
M = 7
d = 3
v0 = X[:, 0]
v1 = X[:, 1]
...
使用 for 循环,我可以这样做,例如:
v1 = np.array([1, 2, 3]).reshape((-1, 1))
v2 = np.array([10, 20, 30]).reshape((-1, 1))
v3 = np.array([100, 200, 300]).reshape((-1, 1))
v4 = np.array([100.1, 200.1, 300.1]).reshape((-1, 1))
v5 = np.array([1.1, 2.2, 3.3]).reshape((-1, 1))
X = np.hstack((v1, v2, v3, v4, v5))
d = 2
X_ = np.zeros((d * X.shape[0], X.shape[1]+1-d))
for i in range (d):
X_[i*X.shape[0]:(i+1) * X.shape[0], :] = X[:X.shape[0], i:i+(X.shape[1]+1-d)]
我得到:
X_ = array([[ 1. , 10. , 100. , 100.1],
[ 2. , 20. , 200. , 200.1],
[ 3. , 30. , 300. , 300.1],
[ 10. , 100. , 100.1, 1.1],
[ 20. , 200. , 200.1, 2.2],
[ 30. , 300. , 300.1, 3.3]]) #Which is the wanted matrix
有没有什么方法可以用矢量化的方式构造这个矩阵(我想当涉及到大矩阵时,这会比 for 循环更快?)。
谢谢。
这看起来是最佳的;你已经很好地矢量化了它。我能做的唯一改进是用 np.empty
替换 np.zeros
,这会跳过初始化数组。我尝试使用 np.vstack
和 np.lib.stride_tricks.sliding_window_view
(在 np.empty
.
# sliding window:
X_ = np.lib.stride_tricks.sliding_window_view(X, (X.shape[0], X.shape[1]+1-d)).reshape(d*X.shape[0], -1)
# np.vstack:
X_ = np.vstack([X[:, i:i+(X.shape[1]+1-d)] for i in range(d)])