python - 矢量化 Python 代码

Question

我编写了一个克里金算法，但我发现它很慢。特别是，您是否知道我如何将下面 cons 函数中的代码片段矢量化:

import time
import numpy as np

B = np.zeros((200, 6))
P = np.zeros((len(B), len(B)))

def cons():
  time1=time.time()
  for i in range(len(B)):
    for j in range(len(B)):
      P[i,j] = corr(B[i], B[j])
  time2=time.time()
  return time2-time1

def corr(x,x_i):
  return np.exp(-np.sum(np.abs(np.array(x) - np.array(x_i))))    

time_av = 0.
for i in range(30):
  time_av+=cons()
print "Average=", time_av/100.

---

编辑:奖励问题

1. 如果我希望 corr(B[i], C[j]) 的 C 维度与 B 相同，广播解决​​方案会怎样

2. 如果我的 p-范数订单是一个数组，scipy 解决方案会发生什么:

   p=np.array([1.,2.,1.,2.,1.,2.])
   def corr(x, x_i):
     return np.exp(-np.sum(np.abs(np.array(x) - np.array(x_i))**p))  
   

   对于 2.，我尝试了 P = np.exp(-cdist(B, C,'minkowski', p)) 但 scipy 需要一个标量。

Answer 1

您的问题似乎很容易向量化。对于您要计算的 B 的每一对行

P[i,j] = np.exp(-np.sum(np.abs(B[i,:] - B[j,:])))

您可以利用数组广播并引入第三个维度，对最后一个维度求和:

P2 = np.exp(-np.sum(np.abs(B[:,None,:] - B),axis=-1))

想法是将第一次出现的 B reshape 为 (N,1,M) 而第二个 B 则保留形状 (N,M)。对于数组广播，后者等价于(1,N,M)，所以

B[:,None,:] - B

形状为(N,N,M)。沿着最后一个索引求和将得到您正在寻找的 (N,N)-shape 相关数组。

---

请注意，如果您使用的是 scipy，则可以使用 scipy.spatial.distance.cdist 执行此操作(或者，等效地，scipy.spatial.distance.pdist 和 scipy.spatial.distance.squareform 的组合)，无需计算此对称矩阵的下三角一半.使用@Divakar以这种方式对最简单的解决方案的评论中的建议:

from scipy.spatial.distance import cdist
P3 = 1/np.exp(cdist(B, B, 'minkowski',1))

cdist 将计算 1-范数的 Minkowski 距离，它恰好是坐标差的绝对值之和。