首页 文章

广播视图不规则地numpy

提问于
浏览
9

假设我想要一个大小为 (n,m) 的numpy数组,其中 n 非常大,但是有很多重复,即 . 0:n1 是相同的, n1:n2 是相同的等等( n2%n1!=0 ,即不是规则的间隔) . 有没有办法只为每个重复项存储一组值,同时拥有整个数组的视图?

例:

unique_values = np.array([[1,1,1], [2,2,2] ,[3,3,3]]) #these are the values i want to store in memory
index_mapping = np.array([0,0,1,1,1,2,2]) # a mapping between index of array above, with array below

unique_values_view = np.array([[1,1,1],[1,1,1],[2,2,2],[2,2,2],[2,2,2], [3,3,3],[3,3,3]]) #this is how I want the view to look like for broadcasting reasons

我计划将数组(视图)乘以一些大小为 (1,m) 的其他数组,并取这个产品的点积:

other_array1 = np.arange(unique_values.shape[1]).reshape(1,-1) # (1,m)
other_array2 = 2*np.ones((unique_values.shape[1],1)) # (m,1)
output = np.dot(unique_values_view * other_array1, other_array2).squeeze()

输出是长度为 n 的1D数组 .

2 回答

  • 6

    您的表达式允许两个重要的优化:

    • 做最后的索引

    • 首先将 other_array1other_array2 相乘,然后使用 unique_values

    让我们应用这些优化:

    >>> output_pp = (unique_values @ (other_array1.ravel() * other_array2.ravel()))[index_mapping]
    
    # check for correctness
    >>> (output == output_pp).all()
    True
    
    # and compare it to @Yakym Pirozhenko's approach
    >>> from timeit import timeit
    >>> print("yp:", timeit("np.dot(unique_values * other_array1, other_array2).squeeze()[index_mapping]", globals=globals()))
    yp: 3.9105667411349714
    >>> print("pp:", timeit("(unique_values @ (other_array1.ravel() * other_array2.ravel()))[index_mapping]", globals=globals()))
    pp: 2.2684884609188884
    

    如果我们观察到两件事情,很容易发现这些优化:

    (1)如果 Amxn -matrix而 bn -vector那么

    A * b == A @ diag(b)
    A.T * b[:, None] == diag(b) @ A.T
    

    (2)如果 A 是一个 mxn -matrix而且 Ikk 整数,则 range(m) 然后

    A[I] == onehot(I) @ A
    

    onehot 可以定义为

    def onehot(I, m, dtype=int):
        out = np.zeros((I.size, m), dtype=dtype)
        out[np.arange(I.size), I] = 1
        return out
    

    使用这些事实并缩写 uvimoa1oa2 我们可以写

    uv[im] * oa1 @ oa2 == onehot(im) @ uv @ diag(oa1) @ oa2
    

    现在,上述优化只是选择这些矩阵乘法的最佳阶数

    onehot(im) @ (uv @ (diag(oa1) @ oa2))
    

    使用(1)和(2)向后,我们从这篇文章的开头获得优化的表达式 .

  • 1

    根据您的示例,您可以简单地将通过计算的索引映射计算到最后:

    output2 = np.dot(unique_values * other_array1, other_array2).squeeze()[index_mapping]
    
    assert (output == output2).all()
    

相关问题