代码之家  ›  专栏  ›  技术社区  ›  nullgeppetto

Python中相似矩阵的高效计算(NumPy)

  •  5
  • nullgeppetto  · 技术社区  · 6 年前

    允许 X 成为一名 Bxn numpy 矩阵,即。,

    import numpy as np
    B = 10
    n = 2
    X = np.random.random((B, n))
    

    现在,我对计算所谓的内核(甚至相似性)矩阵感兴趣 K ,它的形状 BxB ,及其 {i,j} -th元素如下所示:

    K(i,j)=乐趣(x\u i,x\u j)

    哪里 x_t 表示 t -矩阵第行 十、 fun 是的一些函数 x_i , x_j . 例如,该函数可以是所谓的RBF函数,即。,

    K(i,j)=exp(-x\u i-x\u 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)
    

    为了提高效率,我想以矢量化的方式执行上述操作。你能帮忙吗?

    3 回复  |  直到 6 年前
        1
  •  6
  •   MB-F    6 年前

    如果你利用 broadcasting .

    您只需以矢量化的方式对内部距离范数计算进行编码:

    X1 = X[:, np.newaxis, :]
    X2 = X[np.newaxis, :, :]
    K = np.exp(-np.sum((X1 - X2)**2, axis=-1))
    
        2
  •  5
  •   max9111    6 年前

    不要矢量化它,只需编译它

    这几乎每次都更快,代码更容易阅读。 因为一个好的jit编译器 Numba 是可用的,这是一件非常简单的事情。

    在您的情况下:

    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的时间,因为这里的代码已经编译好了,但我想您应该多次调用这个函数。

    如果使用单线程版本,还可以缓存编译结果。多线程代码的缓存可能很快就会实现。

        3
  •  3
  •   Giuseppe Angora    6 年前

    我不确定您是否可以仅使用numpy来完成此任务。我会用这个方法 cdist 从scipy库中,类似这样的内容:

    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.        ]])
    

    希望这能帮助你。干得好

    编辑 对于theano实现,您也只能使用numpy数组:

    dist = (X ** 2).sum(1).reshape((X.shape[0], 1)) + (X ** 2).sum(1).reshape((1, X.shape[0])) - 2 * X.dot(X.T)
    

    这应该是工作!