【问题标题】:Dot product of two numpy arrays with 3D Vectors具有 3D 向量的两个 numpy 数组的点积
【发布时间】:2020-11-27 17:53:25
【问题描述】:

我的目标是找到离单个点最近的段(在段数组中)。 获取 2D 坐标数组之间的点积是可行的,但使用 3D 坐标会出现以下错误:

*ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 2 is different from 3)*


A = np.array([[1,1,1],[2,2,2]])
B = np.array([[3,3,3], [4,4,4]])

dp = np.dot(A,B)

dp 应该返回 2 个值, [1,1,1]@[3,3,3][2,2,2]@[4,4,4] 的点积

// 谢谢大家。

这是找到离单点最近的线段的最终解决方案。
欢迎任何优化。

import numpy as np
import time

#find closest segment to single point

then = time.time()

#random line segment
l1 = np.random.rand(1000000, 3)*10   
l2 = np.random.rand(1000000, 3)*10

#single point
p = np.array([5,5,5]) #only single point

#set to origin
line = l2-l1
pv = p-l1  

#length of line squared
len_sq = np.sum(line**2, axis = 1) #len_sq = numpy.einsum("ij,ij->i", line, line)

#dot product of 3D vectors with einsum
dot = np.einsum('ij,ij->i',line,pv) #np.sum(line*pv,axis=1)


#percentage of line the pv vector travels in
param = np.array([dot/len_sq])

#param<0 projected point=l1, param>1 pp=l2
clamped_param = np.clip(param,0,1)

#add line fraction to l1 to get projected point
pp = l1+(clamped_param.T*line)

##distance vector between single point and projected point
pp_p = pp-p

#sort by smallest distance between projected point and l1
index_of_mininum_dist = np.sum(pp_p**2, axis = 1).argmin()

print(index_of_mininum_dist)
print("FINISHED IN: ", time.time()-then)

【问题讨论】:

    标签: python arrays numpy vectorization dot-product


    【解决方案1】:

    点积是 numpy 显然不是为与数组一起使用而设计的。围绕它编写一些包装器非常容易。比如这样:

    def array_dot(A, B):
        return [A[i]@B[i] for i in range(A.shape[0])]
    

    【讨论】:

    • 谢谢,我没有指定我想要矢量化版本。在我的情况下,np.einsum 快了约 15 倍。
    【解决方案2】:

    np.dot 仅适用于向量,不适用于矩阵。当传递矩阵时,它期望进行矩阵乘法,由于传递的维度,这将失败。

    在矢量上它会像你预期的那样工作:

    np.dot(A[0,:],B[0,:])
    np.dot(A[1,:],B[1,:])
    

    一次性完成:

    np.sum(A*B,axis=1)
    

    【讨论】:

      【解决方案3】:

      你是这个意思吗:

      np.einsum('ij,ij->i',A,B)
      

      输出:

      [ 9 24]
      

      但是,如果您想要 A 中的每一行与 B 中的每一行的点积,您应该这样做:

      A@B.T
      

      输出:

      [[ 9 12]
       [18 24]]
      

      【讨论】:

      • 谢谢,np.einsum 比 np.sum 快一点,但对于来自纯 python 的人来说不是很清楚。
      • @Tortenrandband 欢迎您。代码的选择绝对是您认为更具可读性的主观因素。希望对其他用户有所帮助。
      【解决方案4】:
      In [265]: A = np.array([[1,1,1],[2,2,2]]) 
           ...: B = np.array([[3,3,3], [4,4,4]]) 
      

      元素明智的乘法,然后是 sum 工作正常:

      In [266]: np.sum(A*B, axis=1)                                                                        
      Out[266]: array([ 9, 24])
      

      einsum 也使表达变得容易:

      In [267]: np.einsum('ij,ij->i',A,B)                                                                  
      Out[267]: array([ 9, 24])
      

      dot 带有 2d 数组(此处为 (2,3) 形状),执行矩阵乘法,这是经典的跨行、向下列。在 einsum 表示法中,这是 'ij,jk->ik'。

      In [268]: np.dot(A,B)                                                                                
      ---------------------------------------------------------------------------
      ValueError                                Traceback (most recent call last)
      <ipython-input-268-189f80e2c351> in <module>
      ----> 1 np.dot(A,B)
      
      <__array_function__ internals> in dot(*args, **kwargs)
      
      ValueError: shapes (2,3) and (2,3) not aligned: 3 (dim 1) != 2 (dim 0)
      

      使用转置,尺寸匹配 (2,3) 和 (3,2),但结果是 (2,2):

      In [269]: np.dot(A,B.T)                                                                              
      Out[269]: 
      array([[ 9, 12],
             [18, 24]])
      

      所需的值在对角线上。

      思考问题的一种方式是我们想做一批一维产品。添加了matmul/@ 以执行批量矩阵乘法(dot 不能这样做)。但是数组必须扩展为 3d,因此批处理维度是前导维度(并且 3 在各自的最后一个维度和第二个到最后一个维度上):

      In [270]: A[:,None,:]@B[:,:,None]       # (2,1,3) with (2,3,1)                                                              
      Out[270]: 
      array([[[ 9]],
      
             [[24]]])
      

      但结果是 (2,1,1) 形的。正确的数字在那里,但我们必须挤出额外的维度。

      总体而言,前 2 个解决方案是最简单的 - sum 或 product 或 einsum 等效项。

      【讨论】:

        猜你喜欢
        • 1970-01-01
        • 1970-01-01
        • 2021-04-18
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 2019-11-08
        • 1970-01-01
        • 2016-05-15
        相关资源
        最近更新 更多