如何使用numpy在for循环中对矩阵和进行矢量化？答案

【问题标题】：how to vectorize a matrix sum in a for loop using numpy?如何使用numpy在for循环中对矩阵和进行矢量化？
【发布时间】：2014-08-21 22:13:45
【问题描述】：

基本上我有一个行 = 3600 和列 = 5 的矩阵，并希望将其下采样为 60 行的地块：

import numpy as np

X = np.random.rand(3600,5)

down_sample = 60
ds_rng = range(0,X.shape[0],down_sample)
X_ds = np.zeros((ds_rng.__len__(),X.shape[1]))

i = 0
for j in ds_rng:
    X_ds[i,:] = np.sum( X[j:j+down_sample,:], axis=0 )
    i += 1

【问题讨论】：

最后有 40 个“孤儿”行，这是个问题吗？
你说得对，那 40 行 'orhpan' 不应该存在，我已将其更改为倍数。

标签： python numpy scipy vectorization broadcast

【解决方案1】：

另一种方法可能是：

def blockwise_sum(X, down_sample=60):
    n, m = X.shape

    ds_n = n / down_sample
    N = ds_n * down_sample

    if N == n:
        return np.sum(X.reshape(-1, down_sample, m), axis=1)

    X_ds = np.zeros((ds_n + 1, m))
    X_ds[:ds_n] = np.sum(X[:N].reshape(-1, down_sample, m), axis=1)
    X_ds[-1] = np.sum(X[N:], axis=0)

    return X_ds

我不知道它是否更快。

【讨论】：

谢谢！这正是我所期望的，N == n case 做得很好，而且比 strides 更快。
@pvstrln 很酷。鉴于您编辑的问题，函数体可以简化为：return np.sum(X.reshape(-1, down_sample, X.shape[1]), axis=1)

【解决方案2】：

至少在这种情况下，einsum 比 sum 快。

np.einsum('ijk->ik',x.reshape(-1,down_sample,x.shape[1]))

比 blockwise_sum 快 2 倍。

我的时间安排：

OP iterative  - 1.59 ms
with strided  -   198 us
blockwise_sum -   179 us
einsum        -    76 us

【讨论】：

【解决方案3】：

看起来您可以使用一些大步小技巧来完成工作。

这是我们需要的设置代码：

import numpy as np
X = np.random.rand(1000,5)
down_sample = 60

现在我们欺骗 numpy 认为 X 被分割成包裹：

num_parcels = int(np.ceil(X.shape[0] / float(down_sample)))
X_view = np.lib.stride_tricks.as_strided(X, shape=(num_parcels,down_sample,X.shape[1]))

X_ds = X_view.sum(axis=1)  # sum over the down_sample axis

最后，如果您的下采样间隔没有完全均匀地划分您的行，您需要修复 X_ds 中的最后一行，因为我们采用的步幅技巧使其回绕。

rem = X.shape[0] % down_sample
if rem != 0:
  X_ds[-1] = X[-rem:].sum(axis=0)

【讨论】：

谢谢！您的方法大约花费不到一半的时间：0.000556945800781 vs 0.000196933746338。我想知道这是否可以完全使用矢量化/广播来完成。
结束行的不匹配已得到纠正（OP 选择了 3600 而不是 1000）。
当求和窗口中有重叠时，stride_tricks 更有价值。没有重叠，基本上和reshaping一样。