如何在numpy中对不规则数组形状进行向量化/张量化操作

Question

我想执行这个操作： 如果 有一个规则的形状，那么我可以使用np.einsum，我相信语法应该是：

np.einsum('ijp,ipk->ijk',X, alpha)

很不幸，我的数据X在第一（如果我们从零索引）轴上具有非常规结构。

为了更好地说明， 指的是第i组中第j个成员的第p个特征。由于不同组大小不同，因此实际上它是一个长度不同的列表的列表，每个列表的长度相同。 具有规则结构，因此可以保存为标准的numpy数组（它首先是1维的，然后我使用alpha.reshape（a，b，c），其中a，b，c是问题特定的整数）。

我想避免将X存储为列表的列表或不同维度的np.array列表，并编写类似以下内容的代码：

A = []
for i in range(num_groups):
    temp = np.empty(group_sizes[i], dtype=float)
    for j in range(group_sizes[i]):
        temp[i] = np.einsum('p,pk->k',X[i][j], alpha[i,:,:])
    A.append(temp)

有没有一些不错的numpy函数或数据结构可以用来完成这个操作？或者我必须妥协，采用部分向量化的实现方式？

- gazza89

---

你能添加一个样例吗？ - Divakar

同意；特别是，我想知道涉及的不同轴的典型大小是多少；以及group_sizes中的整数是如何分布的。 group_sizes是否偶然具有相对较少的唯一元素？ - Eelco Hoogendoorn

组大小从几十到数千不等。总共有82个组。有人建议我可以通过添加零来使结构规律化。这在数学上是可行的，但由于最大的组比最小的组大得多（约2个数量级），所以将包含比必要更多的操作。 - gazza89

@gazza89 你可能希望使用一些分区策略，例如将所有0-20组和所有20-50、50-500等等的组聚类到它们自己的矩阵中，然后你可以在这些矩阵上运行向量化操作，这将避免许多不必要的计算。不过找到最佳分区可能会很困难（我认为是NP难问题）。 - Tom Wyllie

Answer 1

我知道这听起来很明显，但是如果你有足够的内存，我建议你首先检查通过将数据填充为统一大小（即仅添加零并执行操作）所获得的性能。有时候，一个更简单的解决方案比一个更优化的解决方案更快，尽管后者在Python/C往返中具有更多的优势。

如果这种方法不起作用，那么你最好的选择就像Tom Wyllie建议的那样采用桶策略。假设X是你的列表列表列表，alpha是一个数组，你可以从收集第二个索引的大小开始（也许你已经有了这个）。

X_sizes = np.array([len(x_i) for x_i in X])

并将它们排序：

idx_sort = np.argsort(X_sizes)
X_sizes_sorted = X_sizes[idx_sort]

然后你需要选择一定数量的桶，这就是你工作分割的数量。假设你选择了BUCKETS = 4。你只需要将数据分成相同大小的块：

sizes_cumsum = np.cumsum(X_sizes_sorted)
total = sizes_cumsum[-1]
bucket_idx = []
for i in range(BUCKETS):
    low = np.round(i * total / float(BUCKETS))
    high = np.round((i + 1) * total / float(BUCKETS))
    m = sizes_cumsum >= low & sizes_cumsum < high
    idx = np.where(m),
    # Make relative to X, not idx_sort
    idx = idx_sort[idx]
    bucket_idx.append(idx)

然后你需要针对每个桶进行计算：

bucket_results = []
for idx in bucket_idx:
    # The last index in the bucket will be the biggest
    bucket_size = X_sizes[idx[-1]]
    # Fill bucket array
    X_bucket = np.zeros((len(X), bucket_size, len(X[0][0])), dtype=X.dtype)
    for i, X_i in enumerate(idx):
        X_bucket[i, :X_sizes[X_i]] = X[X_i]
    # Compute
    res = np.einsum('ijp,ipk->ijk',X, alpha[:, :bucket_size, :])
    bucket_results.append(res)

在这部分中，填充数组X_bucket可能会很慢。如果您可以承受内存，将X放在单个填充的数组中，然后只需切片X[idx,:bucket_size,:]将更有效率。

最后，您可以将结果放回到列表中：

result = [None] * len(X)
for res, idx in zip(bucket_results, bucket_idx):
    for r, X_i in zip(res, idx):
        result[X_i] = res[:X_sizes[X_i]]

抱歉，我没有提供一个完整的函数，但是我不确定您的输入或期望输出的确切情况，因此我只是提供了一些代码片段，您可以根据自己的需要使用它们。

- jdehesa