【问题标题】:Broadcasting N-dim array to (N+1)-dim array and summing on all but 1 dim将 N-dim 数组广播到 (N+1)-dim 数组并在除 1dim 之外的所有数组上求和
【发布时间】:2019-09-18 00:02:18
【问题描述】:

假设您有一个形状为 (a,b,c) 的 numpy 数组和一个形状为 (a,b,c,d) 的布尔掩码。 我想将掩码应用于遍历最后一个轴的数组,沿前三个轴对掩码数组求和,并获得长度/形状(d,)的列表(或数组)。 我尝试通过列表理解来做到这一点:

Result = [np.sum(Array[Mask[:,:,:,i]], axis=(0,1,2)) for i in range(d)]

它可以工作,但它看起来不是很pythonic,而且它也有点慢。 我也尝试过类似的东西

Array = Array[:,:,:,np.newaxis]
Result = np.sum(Array[Mask], axis=(0,1,2))

但这当然行不通,因为 Mask 沿最后一个轴的维度 d 大于数组最后一个轴的维度 1。 此外,考虑到每个轴的维度可能为 100 或 200,因此使用 np.repeat 沿新的最后一个轴重复 Array d 次确实会占用大量内存,我想避免这种情况。 列表理解还有其他更快、更 Python 的替代方案吗?

【问题讨论】:

  • 一种廉价的repeat 方法是:np.broadcast_to(arr[...,None], mask.shape)[mask]。但结果是 1d,丢失了所有“按行”信息。通常,每行(最后一个维度)的 True 值的数量会有所不同。

标签: python arrays numpy array-broadcasting


【解决方案1】:

将 N 维数组广播到匹配的 (N+1) 维数组的最直接方法是使用np.broadcast_to()

import numpy as np


arr = np.random.randint(0, 100, (2, 3))
mask = np.random.randint(0, 2, (2, 3, 4), dtype=bool)
b_arr = np.broadcast_to(arr[..., None], mask.shape)
print(mask.shape == b_arr.shape)
# True

但是,正如@hpaulj 已经指出的那样,您不能使用maskb_arr 进行切片而不丢失尺寸。


鉴于您只想将元素加在一起并将零相加“不会造成伤害”,您可以简单地将数组和掩码按元素相乘,以保持正确的维度,但 False 中的元素掩码与相应数组元素的后续sum 无关:

result = np.sum(b_arr * mask, axis=tuple(range(mask.ndim - 1)))

或者,因为* 会自动进行广播:

result = np.sum(arr[..., None] * mask, axis=tuple(range(mask.ndim - 1)))

首先不需要使用np.broadcast_to()(但您仍然需要匹配维度的数量,即使用arr[..., None] 而不仅仅是arr)。


作为@PaulPanzer already pointed out,由于您想要总结除一维之外的所有维度,因此可以使用np.matmul()/@ 进一步简化:

result2 = arr.ravel() @ mask.reshape(-1, mask.shape[-1])
print(np.all(result == result2))
# True

对于涉及求和的更高级的操作,请查看np.einsum()


编辑

广播的问题是它会在评估表达式期间创建临时数组。

对于您似乎正在处理的数字,当我遇到MemoryError 时,我根本无法使用广播数组,但从时间上看,元素乘法可能仍然是比您最初建议的更好的方法。

或者,如果您追求速度,您可以通过 Cython 或 Numba 中的显式循环在较低级别上执行此操作。

您可以在下面找到几个基于 Numba 的解决方案(处理 ravel()-ed 数据):

  • _vector_matrix_product(): 不使用任何临时数组
  • _vector_matrix_product_mp(): 部分同上,但使用并行执行
  • _vector_matrix_product_sum():使用np.sum() 和并行执行
import numpy as np
import numba as nb


@nb.jit(nopython=True)
def _vector_matrix_product(
        vect_arr,
        mat_arr,
        result_arr):
    rows, cols = mat_arr.shape
    if vect_arr.shape == result_arr.shape:
        for i in range(rows):
            for j in range(cols):
                result_arr[i] += vect_arr[j] * mat_arr[i, j]
    else:
        for i in range(rows):
            for j in range(cols):            
                result_arr[j] += vect_arr[i] * mat_arr[i, j]


@nb.jit(nopython=True, parallel=True)
def _vector_matrix_product_mp(
        vect_arr,
        mat_arr,
        result_arr):
    rows, cols = mat_arr.shape
    if vect_arr.shape == result_arr.shape:
        for i in nb.prange(rows):
            for j in nb.prange(cols):
                result_arr[i] += vect_arr[j] * mat_arr[i, j]
    else:
        for i in nb.prange(rows):
            for j in nb.prange(cols):        
                result_arr[j] += vect_arr[i] * mat_arr[i, j]


@nb.jit(nopython=True, parallel=True)
def _vector_matrix_product_sum(
        vect_arr,
        mat_arr,
        result_arr):
    rows, cols = mat_arr.shape
    if vect_arr.shape == result_arr.shape:
        for i in nb.prange(rows):
            result_arr[i] = np.sum(vect_arr * mat_arr[i, :])
    else:
        for j in nb.prange(cols):
            result_arr[j] = np.sum(vect_arr * mat_arr[:, j])


def vector_matrix_product(
        vect_arr,
        mat_arr,
        swap=False,
        dtype=None,
        mode=None):
    rows, cols = mat_arr.shape
    if not dtype:
        dtype = (vect_arr[0] * mat_arr[0, 0]).dtype
    if not swap:
        result_arr = np.zeros(cols, dtype=dtype)
    else:
        result_arr = np.zeros(rows, dtype=dtype)
    if mode == 'sum':
        _vector_matrix_product_sum(vect_arr, mat_arr, result_arr)
    elif mode == 'mp':
        _vector_matrix_product_mp(vect_arr, mat_arr, result_arr)
    else:
        _vector_matrix_product(vect_arr, mat_arr, result_arr)
    return result_arr


np.random.seed(0)
arr = np.random.randint(0, 100, (2, 3, 4))
mask = np.random.randint(0, 2, (2, 3, 4, 5), dtype=bool)
target = arr.ravel() @ mask.reshape(-1, mask.shape[-1])
print(target)
# [820 723 861 486 408]
result1 = vector_matrix_product(arr.ravel(), mask.reshape(-1, mask.shape[-1]))
print(result1)
# [820 723 861 486 408]
result2 = vector_matrix_product(arr.ravel(), mask.reshape(-1, mask.shape[-1]), mode='mp')
print(result2)
# [820 723 861 486 408]
result3 = vector_matrix_product(arr.ravel(), mask.reshape(-1, mask.shape[-1]), mode='sum')
print(result3)
# [820 723 861 486 408]

与任何基于 list-comprehension 的解决方案相比,时间有所改进:

arr = np.random.randint(0, 100, (256, 256, 256))
mask = np.random.randint(0, 2, (256, 256, 256, 128), dtype=bool)


%timeit np.sum(arr[..., None] * mask, axis=tuple(range(mask.ndim - 1)))
# MemoryError

%timeit arr.ravel() @ mask.reshape(-1, mask.shape[-1])
# MemoryError

%timeit np.array([np.sum(arr * mask[..., i], axis=tuple(range(mask.ndim - 1))) for i in range(mask.shape[-1])])
# 24.1 s ± 105 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit np.array([np.sum(arr[mask[..., i]]) for i in range(mask.shape[-1])])
# 46 s ± 119 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit vector_matrix_product(arr.ravel(), mask.reshape(-1, mask.shape[-1]))
# 408 ms ± 2.12 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit vector_matrix_product(arr.ravel(), mask.reshape(-1, mask.shape[-1]), mode='mp')
# 1.63 s ± 3.58 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit vector_matrix_product(arr.ravel(), mask.reshape(-1, mask.shape[-1]), mode='sum')
# 7.17 s ± 258 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

正如预期的那样,JIT 加速版本是最快的,并且对代码强制执行并行性并不会提高速度。 另请注意,逐元素乘法的方法比切片更快(这些基准的速度大约是切片的两倍)。


编辑 2

按照@max9111 的建议,首先按行循环,然后按列循环会导致最耗时的循环在连续数据上运行,从而显着提高速度。 如果没有这个技巧,_vector_matrix_product_sum()_vector_matrix_product_mp() 将以基本相同的速度运行。

【讨论】:

  • 这些都是非常棒的想法,而且都很有效;但是,对于我正在使用的数组的大小(通常我有 a=b=c=256 和 d =128),我要么得到大量内存消耗,要么当我解决这个问题时,这些解决方案似乎比列表理解。这是正常的吗?还是我错过了什么?
  • @Quasark 查看编辑。基本上,是的,广播是内存效率低下的,您需要变通方法来获得一些速度。一般来说,这个操作需要时间处理您拥有的数字(仅mask 就有 20 亿个条目)。
  • 循环排序在您的代码中不是最理想的。如果您交换循环顺序,您可以在 _vector_matrix_product 中获得很大的加速。就我而言,它是从 15.3 秒到 193 毫秒。并行版本不是必需的,因为它比单线程版本运行得慢
  • @max9111 不错!我更新了答案以反映您的建议。
  • @Quasark 查看最新更新以显着提高内存效率。
【解决方案2】:

怎么样

Array.reshape(-1)@Mask.reshape(-1,d)

既然你是在前三个轴上求和,你也可以合并它们,之后很容易看出该操作可以写成矩阵向量积

例子:

a,b,c,d = 4,5,6,7
Mask = np.random.randint(0,2,(a,b,c,d),bool)
Array = np.random.randint(0,10,(a,b,c))
[np.sum(Array[Mask[:,:,:,i]]) for i in range(d)]
# [310, 237, 253, 261, 229, 268, 184]    
Array.reshape(-1)@Mask.reshape(-1,d)
# array([310, 237, 253, 261, 229, 268, 184])

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 2022-12-11
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2021-06-18
    • 1970-01-01
    相关资源
    最近更新 更多