是否有可能在Python中加速这个循环？

Question

映射函数的常规方法是numpy.narray类似np.array[map(some_func,x)]或vectorize(f)(x)不能提供索引.以下代码只是许多应用程序中常见的简单示例.

dis_mat = np.zeros([feature_mat.shape[0], feature_mat.shape[0]])

for i in range(feature_mat.shape[0]):
    for j in range(i, feature_mat.shape[0]):
        dis_mat[i, j] = np.linalg.norm(
            feature_mat[i, :] - feature_mat[j, :]
        )
        dis_mat[j, i] = dis_mat[i, j]

有没有办法加快速度？

---

谢谢您的帮助!加速此代码的最快方法是使用@ user2357112评论的函数:

    from scipy.spatial.distance import pdist,squareform
    dis_mat = squareform(pdist(feature_mat))

@Julien的方法也很好,如果feature_mat很小,但是当feature_mat1000是2000,那么它需要近40 GB的内存.

Answer 1

SciPy具有专门用于计算您正在计算的成对距离类型的功能.它是scipy.spatial.distance.pdist,并且它以精简格式生成距离,基本上只存储距离矩阵的上三角形,但您可以将结果转换为方形scipy.spatial.distance.squareform:

from scipy.spatial.distance import pdist, squareform

distance_matrix = squareform(pdist(feature_mat))

这样可以避免使用直接矢量化解决方案所需的巨大中间阵列,因此速度更快,适用于更大的输入.然而,它失去了使用代数操作来dot处理繁重工作的方法的时机.

pdist 如果您决定要使用欧几里德距离以外的其他东西,也会支持各种各样的备用距离指标.

# Manhattan distance!
distance_matrix = squareform(pdist(feature_mat, 'cityblock'))

# Cosine distance!
distance_matrix = squareform(pdist(feature_mat, 'cosine'))

# Correlation distance!
distance_matrix = squareform(pdist(feature_mat, 'correlation'))

# And more! Check out the docs.

Answer 2

您可以创建新轴并广播:

dis_mat = np.linalg.norm(feature_mat[:,None] - feature_mat, axis=-1)

定时:

feature_mat = np.random.rand(100,200)

def a():
    dis_mat = np.zeros([feature_mat.shape[0], feature_mat.shape[0]])
    for i in range(feature_mat.shape[0]):
        for j in range(i, feature_mat.shape[0]):
            dis_mat[i, j] = np.linalg.norm(
                feature_mat[i, :] - feature_mat[j, :]
            )
            dis_mat[j, i] = dis_mat[i, j]

def b():
    dis_mat = np.linalg.norm(feature_mat[:,None] - feature_mat, axis=-1)



%timeit a()
100 loops, best of 3: 20.5 ms per loop

%timeit b()
100 loops, best of 3: 11.8 ms per loop

Answer 3

考虑可以做什么,并np.dot在k x k矩阵中使用优化,在小内存位置(kxk):

def c(m): 
    xy=np.dot(m,m.T) # O(k^3)
    x2=y2=(m*m).sum(1) #O(k^2)
    d2=np.add.outer(x2,y2)-2*xy  #O(k^2)
    d2.flat[::len(m)+1]=0 # Rounding issues
    return np.sqrt(d2)  # O (k^2)

为了比较:

def d(m):
   return  squareform(pdist(m))

以下是ak*k初始矩阵的'时间(it)':

这两种算法是O(k ^ 3),但是c(m)通过np.dot线性代数的关键节点,使得O(k ^ 3)部分工作,这有利于所有优化(多核等).pdist只是在源中看到的循环.

这解释了大数组的15倍因子,即使pdist通过仅计算一半项来利用矩阵的对称性.