Python中高效的逐元素函数计算答案

【问题标题】：Efficient element-wise function computation in PythonPython中高效的逐元素函数计算
【发布时间】：2015-07-17 03:42:54
【问题描述】：

我有以下优化问题。给定两个 np.arrays X,Y 和一个函数 K 我想尽可能快地计算矩阵关联 gram_matrix，其中 (i,j)-th 元素被计算为 K(X[i],Y[j])。

这里有一个使用嵌套 for 循环的实现，它被认为是解决这类问题最慢的。

def proxy_kernel(X,Y,K):
    gram_matrix = np.zeros((X.shape[0], Y.shape[0]))
    for i, x in enumerate(X):
        for j, y in enumerate(Y):
            gram_matrix[i, j] = K(x, y)
    return gram_matrix

非常感谢任何帮助。

【问题讨论】：

最好在 Code Review 上问这个问题。
如果不知道 K 是什么，你不会比嵌套循环做得更好。
不幸的是 K 作为参数传递。很遗憾，因为这个实现真的很慢，但到目前为止是唯一有效的。
这是一个基本的numpy 操作方法问题；它会在这里得到比 Code Review 更多的答案。
这已经被问过很多次了。只要K 是一个只接受标量的黑盒，您就无法提高速度。但首先，将典型 K 的时间与最简单的 K 的时间进行比较。迭代花费了多少时间，以及评估K 无数次的时间？

标签： python numpy scikit-learn vectorization

【解决方案1】：

你肯定至少可以矢量化内循环：

def proxy_kernel_vect(X, Y, K):
    K_vect = np.vectorize(K)
    gram_matrix = np.zeros((X.shape[0], Y.shape[0]))
    for i, x in enumerate(X):
        gram_matrix[i] = K_vect(x, Y)
    return gram_matrix

这对相对较长的数组产生了很好的改进：

In [15]: a = np.array(range(1000))
    ...: b = np.array(range(1000))
    ...: 

In [16]: %timeit proxy_kernel(a, b, k)
1 loops, best of 3: 665 ms per loop

In [17]: %timeit proxy_kernel_vect(a, b, k)
1 loops, best of 3: 266 ms per loop

其中k 就是lambda x, y: x+y。

【讨论】：

这里的k 是什么？我不认为以这种方式应用vectorize 实际上保留了原始语义。
@user2357112 它在答案的最后一行中定义。它只是一个标量函数（因为这是 OP 似乎使用的）。
我不认为它应该是一个标量函数。在“格拉姆矩阵”的通常定义中，函数是内积。
@Bakuriu 谢谢。不幸的是，这在我的代码中不起作用gram_matrix[i] = vect(x, Y) ValueError: could not broadcast input array from shape (600,2) into shape (600)

【解决方案2】：

np.vectorize 确实使速度有所提高 - 大约 2 倍（这里我使用 math.atan2 作为一个黑盒函数，它需要 2 个标量参数）。

In [486]: X=np.linspace(0,1,100)
In [487]: K_vect=np.vectorize(math.atan2)

In [488]: timeit proxy_kernel(X,X,math.atan2)
100 loops, best of 3: 7.84 ms per loop

In [489]: timeit K_vect(X[:,None],X)
100 loops, best of 3: 3.49 ms per loop

In [502]: timeit np.arctan2(X[:,None],X)  # numpy compiled loop
1000 loops, best of 3: 871 µs per loop

在哪里

def proxy_kernel(X,Y,K):
    gram_matrix = np.zeros((X.shape[0], Y.shape[0]))
    for i, x in enumerate(X):
            for j, y in enumerate(Y):
                    gram_matrix[i, j] = K(x, y)
    return gram_matrix

只要K 是一个黑匣子，你就会受到调用K 和X.shape[0]*Y.shape[0] 次数的限制。您可以尝试最小化迭代时间，但您仍然受到所有这些函数调用的限制。

https://stackoverflow.com/a/29733040/901925 通过利用 np.linalg.norm 函数的 axis 参数来加快计算速度。

【讨论】：

鉴于名称 gram_matrix，我很确定 K 应该是两个向量的函数，而不是两个标量。
在 CR 上他引用了机器学习内核函数；但他仍然不清楚尺寸； codereview.stackexchange.com/questions/90005/…
为什么使用矢量化可以提高速度？我认为矢量化只是让代码更漂亮，但在内部仍然只是一个 for 循环？
在这种情况下我不知道。 vectorize 使用 frompyfunc。 frompyfunc 的测试速度通常比显式循环快 2 倍，但 vectorize 的开销更大，而且速度通常有点慢。但无论如何，我都不会过多关注 2 倍的差异。 np.arctan2 提供的 10 倍差异是我们想要的。

【解决方案3】：

您也可以尝试numba 模块中的vectorize 装饰器。

使用vectorize 和 numpy broadcasting 可以轻松解决您的特定问题：

from numba import vectorize

@vectorize(['float64(float64, float64)'])
def K(x,y):
    return x + y

a = np.arange(1000)
b = np.arange(1000)

gram_array = K(a[None,:], b[:, None])

【讨论】：