用python计算RBF内核的最快方法是什么？

Question

我想计算一个数据矩阵X的RBF或“高斯”核，该数据矩阵有n行和d列。由此得到的平方核矩阵由下式给出：

K[i,j] = var * exp(-gamma * ||X[i] - X[j]||^2)

var和gamma是标量。

在python中，最快的方法是什么？

Answer 1

你在你的^{}中做了很多优化。我想再增加一些（主要是微调）。我将以答案帖子中的获胜者为基础，这似乎是基于numexpr。

调整1

首先，np.sum(X ** 2, axis = -1)可以用np.einsum进行优化。虽然这部分不是最大的开销，但是任何优化都不会有什么坏处。所以，这个总和可以表示为-

X_norm = np.einsum('ij,ij->i',X,X)

调整2

其次，我们可以利用Scipy支持的blas函数，如果允许的话，可以使用单精度的dtype来显著提高其性能。因此，np.dot(X, X.T)可以用^{}来计算，就像-

sgemm(alpha=1.0, a=X, b=X, trans_b=True)

在用gamma重新排列负号方面，再做一些调整，我们就可以为sgemm提供更多信息。此外，我们将把gamma推入alpha项。

调整后的实现

因此，通过这两个优化，我们将拥有numexpr方法的另外两个变体（如果我可以这么说的话），如下所示-

from scipy.linalg.blas import sgemm

def app1(X, gamma, var):
    X_norm = -np.einsum('ij,ij->i',X,X)
    return ne.evaluate('v * exp(g * (A + B + 2 * C))', {\
        'A' : X_norm[:,None],\
        'B' : X_norm[None,:],\
        'C' : np.dot(X, X.T),\
        'g' : gamma,\
        'v' : var\
    })

def app2(X, gamma, var):
    X_norm = -gamma*np.einsum('ij,ij->i',X,X)
    return ne.evaluate('v * exp(A + B + C)', {\
        'A' : X_norm[:,None],\
        'B' : X_norm[None,:],\
        'C' : sgemm(alpha=2.0*gamma, a=X, b=X, trans_b=True),\
        'g' : gamma,\
        'v' : var\
    })

运行时测试

从你的回复帖子中找到一个基于Numexpr的-

def app0(X, gamma, var):
    X_norm = np.sum(X ** 2, axis = -1)
    return ne.evaluate('v * exp(-g * (A + B - 2 * C))', {
            'A' : X_norm[:,None],
            'B' : X_norm[None,:],
            'C' : np.dot(X, X.T),
            'g' : gamma,
            'v' : var
    })

计时和验证-

In [165]: # Setup
     ...: X = np.random.randn(10000, 512)
     ...: gamma = 0.01
     ...: var = 5.0

In [166]: %timeit app0(X, gamma, var)
     ...: %timeit app1(X, gamma, var)
     ...: %timeit app2(X, gamma, var)
1 loop, best of 3: 1.25 s per loop
1 loop, best of 3: 1.24 s per loop
1 loop, best of 3: 973 ms per loop

In [167]: np.allclose(app0(X, gamma, var), app1(X, gamma, var))
Out[167]: True

In [168]: np.allclose(app0(X, gamma, var), app2(X, gamma, var))
Out[168]: True

Answer 2

我将介绍四种不同的方法来计算这样一个内核，然后比较它们的运行时。

使用纯numpy

在这里，我使用的事实是||x-y||^2 = ||x||^2 + ||y||^2 - 2 * x^T * y。

import numpy as np

X_norm = np.sum(X ** 2, axis = -1)
K = var * np.exp(-gamma * (X_norm[:,None] + X_norm[None,:] - 2 * np.dot(X, X.T)))

使用numexpr

^{}是一个python包，它允许对numpy数组执行高效的并行数组操作。我们可以使用它执行与上面相同的计算：

import numpy as np
import numexpr as ne

X_norm = np.sum(X ** 2, axis = -1)
K = ne.evaluate('v * exp(-g * (A + B - 2 * C))', {
        'A' : X_norm[:,None],
        'B' : X_norm[None,:],
        'C' : np.dot(X, X.T),
        'g' : gamma,
        'v' : var
})

使用scipy.spatial.distance.pdist

我们还可以使用^{}来计算成对平方欧氏距离的非冗余数组，计算该数组上的核，然后将其转换为平方矩阵：

import numpy as np
from scipy.spatial.distance import pdist, squareform

K = squareform(var * np.exp(-gamma * pdist(X, 'sqeuclidean')))
K[np.arange(K.shape[0]), np.arange(K.shape[1])] = var

使用sklearn.metrics.pairwise.rbf_kernel

sklearn提供直接计算RBF内核的built-in method：

import numpy as np
from sklearn.metrics.pairwise import rbf_kernel

K = var * rbf_kernel(X, gamma = gamma)

运行时比较

我使用25000个512维随机样本在英特尔酷睿i7-7700HQ（4核，2.8ghz）上进行测试和实验。更准确地说：

X = np.random.randn(25000, 512)
gamma = 0.01
var = 5.0

每个方法运行7次，并报告每次执行时间的平均值和标准偏差。

|               Method                |       Time        |
|-------------------------------------|-------------------|
| numpy                               | 24.2 s ± 1.06 s   |
| numexpr                             | 8.89 s ± 314 ms   |
| scipy.spatial.distance.pdist        | 2min 59s ± 312 ms |
| sklearn.metrics.pairwise.rbf_kernel | 13.9 s ± 757 ms   |

首先，scipy.spatial.distance.pdist是惊人的慢。

numexpr几乎比纯的numpy方法快3倍，但是这个加速因子将随可用cpu的数量而变化。

sklearn.metrics.pairwise.rbf_kernel不是最快的方法，但只比numexpr慢一点。