Python(NumPy)中相似度矩阵的高效计算
Efficient computation of similarity matrix in Python (NumPy)
设X
为Bxn
numpy
矩阵,即
import numpy as np
B = 10
n = 2
X = np.random.random((B, n))
现在,我感兴趣的是计算所谓的核(或相似度)矩阵 K
,它的形状是 BxB
,它的第 {i,j}
个元素给出如下:
K(i,j) = 乐趣(x_i, x_j)
其中 x_t
表示矩阵 X
的第 t
行,fun
是 x_i
、x_j
的某个函数。例如,这个函数可以是所谓的 RBF 函数,即
K(i,j) = exp(-|x_i - x_j|^2).
为此,一种天真的方法如下:
K = np.zeros((B, B))
for i in range(X.shape[0]):
x_i = X[i, :]
for j in range(X.shape[0]):
x_j = X[j, :]
K[i, j] = np.exp(-np.linalg.norm(x_i - x_j, 2) ** 2)
我想要的是把上面的操作向量化,为了效率。你能帮忙吗?
我不确定您能否仅使用 numpy 来解决这个问题。我会使用 scipy 库中的方法 cdist,像这样:
import numpy as np
from scipy.spatial.distance import cdist
B=5
X=np.random.rand(B*B).reshape((B,B))
dist = cdist(X, X, metric='euclidean')
K = np.exp(dist)
dist
array([[ 0. , 1.2659804 , 0.98231231, 0.80089176, 1.19326493],
[ 1.2659804 , 0. , 0.72658078, 0.80618767, 0.3776364 ],
[ 0.98231231, 0.72658078, 0. , 0.70205336, 0.81352455],
[ 0.80089176, 0.80618767, 0.70205336, 0. , 0.60025858],
[ 1.19326493, 0.3776364 , 0.81352455, 0.60025858, 0. ]])
K
array([[ 1. , 3.5465681 , 2.67062441, 2.22752646, 3.29783084],
[ 3.5465681 , 1. , 2.06799756, 2.23935453, 1.45883242],
[ 2.67062441, 2.06799756, 1. , 2.01789192, 2.25584482],
[ 2.22752646, 2.23935453, 2.01789192, 1. , 1.82259002],
[ 3.29783084, 1.45883242, 2.25584482, 1.82259002, 1. ]])
希望对您有所帮助。干得好
编辑
您也可以只使用 numpy 数组,用于 theano 实现:
dist = (X ** 2).sum(1).reshape((X.shape[0], 1)) + (X ** 2).sum(1).reshape((1, X.shape[0])) - 2 * X.dot(X.T)
应该可以了!
如果您利用 broadcasting.
的力量,这在 numpy 中当然是可能的
您只需以矢量化方式编写内部 distance-norm 计算代码:
X1 = X[:, np.newaxis, :]
X2 = X[np.newaxis, :, :]
K = np.exp(-np.sum((X1 - X2)**2, axis=-1))
不要向量化,直接编译
这几乎每次都更快并且代码更易于阅读。
由于像 Numba 这样的好 jit 编译器可用,所以这是一件非常简单的事情。
你的情况:
import numpy as np
import numba as nb
@nb.njit(fastmath=True)
def Test_1(X):
K = np.zeros((B, B))
for i in range(X.shape[0]):
x_i = X[i, :]
for j in range(X.shape[0]):
x_j = X[j, :]
K[i, j] = np.exp(-np.linalg.norm(x_i - x_j, 2) ** 2)
return K
函数并行化也很容易:
import numpy as np
import numba as nb
@nb.njit(fastmath=True,parallel=True)
def Test_1(X):
K = np.zeros((B, B))
for i in nb.prange(X.shape[0]):
x_i = X[i, :]
for j in range(X.shape[0]):
x_j = X[j, :]
K[i, j] = np.exp(-np.linalg.norm(x_i - x_j, 2) ** 2)
return K
这很容易胜过目前提供的所有其他解决方案。第一个函数调用需要大约 0.5s 的时间,因为这里你的代码是编译的,但我猜你想多次调用这个函数。
如果使用single-threaded版本,还可以缓存编译结果。多线程代码的缓存可能很快就会实现。
设X
为Bxn
numpy
矩阵,即
import numpy as np
B = 10
n = 2
X = np.random.random((B, n))
现在,我感兴趣的是计算所谓的核(或相似度)矩阵 K
,它的形状是 BxB
,它的第 {i,j}
个元素给出如下:
K(i,j) = 乐趣(x_i, x_j)
其中 x_t
表示矩阵 X
的第 t
行,fun
是 x_i
、x_j
的某个函数。例如,这个函数可以是所谓的 RBF 函数,即
K(i,j) = exp(-|x_i - x_j|^2).
为此,一种天真的方法如下:
K = np.zeros((B, B))
for i in range(X.shape[0]):
x_i = X[i, :]
for j in range(X.shape[0]):
x_j = X[j, :]
K[i, j] = np.exp(-np.linalg.norm(x_i - x_j, 2) ** 2)
我想要的是把上面的操作向量化,为了效率。你能帮忙吗?
我不确定您能否仅使用 numpy 来解决这个问题。我会使用 scipy 库中的方法 cdist,像这样:
import numpy as np
from scipy.spatial.distance import cdist
B=5
X=np.random.rand(B*B).reshape((B,B))
dist = cdist(X, X, metric='euclidean')
K = np.exp(dist)
dist
array([[ 0. , 1.2659804 , 0.98231231, 0.80089176, 1.19326493],
[ 1.2659804 , 0. , 0.72658078, 0.80618767, 0.3776364 ],
[ 0.98231231, 0.72658078, 0. , 0.70205336, 0.81352455],
[ 0.80089176, 0.80618767, 0.70205336, 0. , 0.60025858],
[ 1.19326493, 0.3776364 , 0.81352455, 0.60025858, 0. ]])
K
array([[ 1. , 3.5465681 , 2.67062441, 2.22752646, 3.29783084],
[ 3.5465681 , 1. , 2.06799756, 2.23935453, 1.45883242],
[ 2.67062441, 2.06799756, 1. , 2.01789192, 2.25584482],
[ 2.22752646, 2.23935453, 2.01789192, 1. , 1.82259002],
[ 3.29783084, 1.45883242, 2.25584482, 1.82259002, 1. ]])
希望对您有所帮助。干得好
编辑 您也可以只使用 numpy 数组,用于 theano 实现:
dist = (X ** 2).sum(1).reshape((X.shape[0], 1)) + (X ** 2).sum(1).reshape((1, X.shape[0])) - 2 * X.dot(X.T)
应该可以了!
如果您利用 broadcasting.
的力量,这在 numpy 中当然是可能的您只需以矢量化方式编写内部 distance-norm 计算代码:
X1 = X[:, np.newaxis, :]
X2 = X[np.newaxis, :, :]
K = np.exp(-np.sum((X1 - X2)**2, axis=-1))
不要向量化,直接编译
这几乎每次都更快并且代码更易于阅读。 由于像 Numba 这样的好 jit 编译器可用,所以这是一件非常简单的事情。
你的情况:
import numpy as np
import numba as nb
@nb.njit(fastmath=True)
def Test_1(X):
K = np.zeros((B, B))
for i in range(X.shape[0]):
x_i = X[i, :]
for j in range(X.shape[0]):
x_j = X[j, :]
K[i, j] = np.exp(-np.linalg.norm(x_i - x_j, 2) ** 2)
return K
函数并行化也很容易:
import numpy as np
import numba as nb
@nb.njit(fastmath=True,parallel=True)
def Test_1(X):
K = np.zeros((B, B))
for i in nb.prange(X.shape[0]):
x_i = X[i, :]
for j in range(X.shape[0]):
x_j = X[j, :]
K[i, j] = np.exp(-np.linalg.norm(x_i - x_j, 2) ** 2)
return K
这很容易胜过目前提供的所有其他解决方案。第一个函数调用需要大约 0.5s 的时间,因为这里你的代码是编译的,但我猜你想多次调用这个函数。
如果使用single-threaded版本,还可以缓存编译结果。多线程代码的缓存可能很快就会实现。