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

Question

我想为X具有n行和d列的数据矩阵计算RBF或"高斯"内核.得到的方形核矩阵由下式给出:

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

var并且gamma是标量.

在python中执行此操作的最快方法是什么？

Answer 1

我将介绍四种不同的计算这种内核的方法,然后比较它们的运行时间.

使用纯粹的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

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

我们还可以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内核的内置方法:

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

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

运行时比较

我使用了25,000个512维随机样本进行测试,并在英特尔酷睿i7-7700HQ(4核@ 2.8 GHz)上进行实验.更确切地说:

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.

Answer 2

好吧，您在中做了很多优化answer post。我想添加更多（主要是调整）。我会以似乎是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)可以这样计算SciPy's sgemm：

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