Numpy quirk：将函数应用于两个 1D 数组的所有对，以获得一个 2D 数组

Posted 2023-02-23

【中文标题】Numpy quirk：将函数应用于两个 1D 数组的所有对，以获得一个 2D 数组

【英文标题】：Numpy quirk: Apply function to all pairs of two 1D arrays, to get one 2D array

【发布时间】：2014-02-09 04:28:39

【问题描述】：

假设我有 2 个一维 (1D) numpy 数组，a 和 b，长度分别为 n1 和 n2。我还有一个函数F(x,y)，它接受两个值。现在我想将该函数应用于我的两个一维数组中的每一对值，因此结果将是一个形状为n1, n2 的二维 numpy 数组。二维数组的i, j 元素将是F(a[i], b[j])。

如果没有大量的 for 循环，我无法找到一种方法，而且我确信在 numpy 中有一种更简单（而且更快！）的方法。

提前致谢！

【问题讨论】：

您是否可能正在寻找带有标量的外部产品类型函数？

【参考方案1】：

您可以使用列表推导来创建一个数组数组：

import numpy as np

# Arrays
a = np.array([1, 2, 3]) # n1 = 3
b = np.array([4, 5]) # n2 = 2

# Your function (just an example)
def f(i, j):
    return i + j

result = np.array([[f(i, j)for j in b ]for i in a])
print result

输出：

[[5 6]
 [6 7]
 [7 8]]

【讨论】：

列表推导比 NumPy 代码的 for 循环好一点。

那是我的。我已经提供了反馈，但我会详细说明。列表推导可以缩短代码，但它们与基于循环的解决方案一样慢。使用 NumPy 时，始终首先使用 NumPy 的矢量化操作寻找解决方案是一个好习惯。

这比其他选项慢得多。

【参考方案2】：

您可以使用numpy broadcasting对这两个数组进行计算，使用newaxis将a变成一个垂直二维数组：

In [11]: a = np.array([1, 2, 3]) # n1 = 3
    ...: b = np.array([4, 5]) # n2 = 2
    ...: #if function is c(i, j) = a(i) + b(j)*2:
    ...: c = a[:, None] + b*2

In [12]: c
Out[12]: 
array([[ 9, 11],
       [10, 12],
       [11, 13]])

基准测试：

In [28]: a = arange(100)

In [29]: b = arange(222)

In [30]: timeit r = np.array([[f(i, j) for j in b] for i in a])
10 loops, best of 3: 29.9 ms per loop

In [31]: timeit c = a[:, None] + b*2
10000 loops, best of 3: 71.6 us per loop

【讨论】：

newaxis 可能是向 ndarray 添加新轴的正式方式，即使它们相同。

是的，它们在内部代码中是相同的。但是newaxis 更容易理解，更像是语法糖。

这看起来像我要找的，我马上试试！

这很好用！这可以扩展到比较两个二维数组吗？所以结果将是两个二维数组之间每对列的函数

假设我有一个 numpy 数组，我想在行对中应用一个函数并累积结果，例如row0+row1 ， this + row3 的输出等，我怎样才能聪明地做到这一点？

【参考方案3】：

我可以建议，如果您的用例更局限于产品，that you use the outer-product?

例如：

import numpy

a = array([0, 1, 2])
b = array([0, 1, 2, 3])

numpy.outer(a,b)

返回

array([[0, 0, 0, 0],
       [0, 1, 2, 3],
       [0, 2, 4, 6]])

然后您可以应用其他转换：

numpy.outer(a,b) + 1

返回

array([[1, 1, 1, 1],
       [1, 2, 3, 4],
       [1, 3, 5, 7]])

这要快得多：

>>> import timeit
>>> timeit.timeit('numpy.array([[i*j for i in a] for j in b])', 'import numpy; a=numpy.arange(3); b=numpy.arange(4)')
31.79583477973938

>>> timeit.timeit('numpy.outer(a,b)', 'import numpy; a=numpy.arange(3); b=numpy.arange(4)')
9.351550102233887
>>> timeit.timeit('numpy.outer(a,b)+1', 'import numpy; a=numpy.arange(3); b=numpy.arange(4)')
12.308301210403442

【讨论】：

很有趣，但这不是假设函数总是取两个输入值的乘积吗？如果你想要一个更复杂的函数怎么办？

这很公平，但相对于您可能想要的其他操作而言，这是一个非常常见的操作。如果你想要一个更复杂的功能，我可以在这里建议我的另一个答案：***.com/questions/21226610/…如果我有点冒昧，你愿意分享你的神秘功能吗？

【参考方案4】：

作为比点积更具可扩展性的另一种替代方法，在嵌套列表推导不到 1/5 到 1/9 的时间内，使用 numpy.newaxis (took a bit more digging to find)：

>>> import numpy
>>> a = numpy.array([0,1,2])
>>> b = numpy.array([0,1,2,3])

这一次，使用幂函数：

>>> pow(a[:,numpy.newaxis], b)
array([[1, 0, 0, 0],
       [1, 1, 1, 1],
       [1, 2, 4, 8]])

与替代品相比：

>>> numpy.array([[pow(i,j) for j in b] for i in a])
array([[1, 0, 0, 0],
       [1, 1, 1, 1],
       [1, 2, 4, 8]])

并比较时机：

>>> import timeit
>>> timeit.timeit('numpy.array([[pow(i,j) for i in a] for j in b])', 'import numpy; a=numpy.arange(3); b=numpy.arange(4)')
31.943181037902832
>>> timeit.timeit('pow(a[:, numpy.newaxis], b)', 'import numpy; a=numpy.arange(3); b=numpy.arange(4)')
5.985810041427612

>>> timeit.timeit('numpy.array([[pow(i,j) for i in a] for j in b])', 'import numpy; a=numpy.arange(10); b=numpy.arange(10)')
109.74687385559082
>>> timeit.timeit('pow(a[:, numpy.newaxis], b)', 'import numpy; a=numpy.arange(10); b=numpy.arange(10)')
11.989138126373291

【参考方案5】：

如果 F() 与广播参数一起使用，请务必使用它，正如其他人所描述的那样。 另一种方法是使用 np.fromfunction （function_on_an_int_grid 会是一个更好的名字。） 以下只是将 int 网格映射到您的 a-b 网格，然后映射到 F()：

import numpy as np

def func_allpairs( F, a, b ):
    """ -> array len(a) x len(b):
        [[ F( a0 b0 )  F( a0 b1 ) ... ]
         [ F( a1 b0 )  F( a1 b1 ) ... ]
         ...
        ]
    """
    def fab( i, j ):
        return F( a[i], b[j] )  # F scalar or vec, e.g. gradient

    return np.fromfunction( fab, (len(a), len(b)), dtype=int )  # -> fab( all pairs )


#...............................................................................
def F( x, y ):
    return x + 10*y

a = np.arange( 100 )
b = np.arange( 222 )
A = func_allpairs( F, a, b )
# %timeit: 1000 loops, best of 3: 241 µs per loop -- imac i5, np 1.9.3

【讨论】：

我真的很喜欢这种方法对于具有两个一维向量的所有对矩阵的泛化性。

【参考方案6】：

如果F 超出您的控制范围，您可以使用numpy.vectorize 自动将其包装为“矢量感知”。我在下面提供了一个工作示例，为了完整起见，我定义了自己的F。这种方法具有简单的优势，但如果您可以控制F，则稍微小心地重写它以正确矢量化可以带来巨大的速度优势

import numpy

n1 = 100
n2 = 200

a = numpy.arange(n1)
b = numpy.arange(n2)

def F(x, y):
    return x + y

# Everything above this is setup, the answer to your question lies here:
fv = numpy.vectorize(F)
r = fv(a[:, numpy.newaxis], b)

在我的电脑上，找到了以下时间，显示了您为“自动”矢量化支付的价格：

%timeit fv(a[:, numpy.newaxis], b)
100 loops, best of 3: 3.58 ms per loop

%timeit F(a[:, numpy.newaxis], b)
10000 loops, best of 3: 38.3 µs per loop

【讨论】：

我建议这个In [5]: %timeit a[:, numpy.newaxis] + b 49.9 µs ± 337 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

漂亮！非常笼统的答案，本身就值得一个明确的问题。其中很多是关于语言的——广播、外积、reduce、具有任意函数的相关器，......在任何情况下你的方法会失败吗？ np.vectorize 可以在任意函数/lambda上工作吗？如果没有？

@jtlz2 numpy.vectorise 适用于所有函数，因为它有效地用作 for 循环（请参阅 the docs 中的注释）。因此，vectorise 是避免自己编写 for 循环的一种快速方法，但正如我所展示的那样，它实际上并没有使代码运行得更快。为此，您需要使用您提到的其他技术之一。