在numpy中计算两个矩阵行之间的夹角

Question

我有两个由3D向量（numpy 1D数组）组成的矩阵，我需要逐行计算向量之间的角度，并将结果返回到一个1D数组中。我知道如何计算两个1D向量之间的角度。请问应该如何正确实现？

*** 计算结果的单位是度而不是弧度。

目前我已经有以下代码：

import numpy as np

A = np.array([[1,0,0],
              [0,1,0],
              [0,0,1]])

B = np.array([[1,0,1],
              [1,1,0],
              [0,1,0]])

def angle(V1,V2):
    """
    angle between vectors V1 and V2 in degrees using
    angle = arccos ( V1 dot V2 / norm(V1) * norm(V2) ) *180/np.pi
    """

    cos_of_angle = V1.dot(V2) / (np.linalg.norm(V1) * np.linalg.norm(V2)) 
    return np.arccos(np.clip(cos_of_angle,-1,1))  * 180/np.pi

注意缩放项180/np.pi用于将弧度转换为角度。

我希望有一个数组：

C = [ angle(A[0],B[0]) , angle(A[1],B[1])...... and so on]

如果有人能帮忙，将不胜感激。

- RRS

---

缩进和括号看起来不对。请查看并进行进一步编辑（如果需要）。 - Divakar

已更正。感谢@Divakar。 - RRS

Answer 1

我们可以使用 einsum 来替代点积计算，使用 axis 参数来进行 norm 运算，从而得到向量化的解决方案，代码如下 -

def angle_rowwise(A, B):
    p1 = np.einsum('ij,ij->i',A,B)
    p2 = np.linalg.norm(A,axis=1)
    p3 = np.linalg.norm(B,axis=1)
    p4 = p1 / (p2*p3)
    return np.arccos(np.clip(p4,-1.0,1.0))

我们可以进一步优化并引入更多的 einsum，特别是用它计算norms。因此，我们可以像这样使用它 -

def angle_rowwise_v2(A, B):
    p1 = np.einsum('ij,ij->i',A,B)
    p2 = np.einsum('ij,ij->i',A,A)
    p3 = np.einsum('ij,ij->i',B,B)
    p4 = p1 / np.sqrt(p2*p3)
    return np.arccos(np.clip(p4,-1.0,1.0))

因此，为了解决我们的问题并获得以度为单位的输出 -

out = angle_rowwise(A, B)*180/np.pi

- Divakar

---

这些角度是相对于原点（0, 0, 0）的吗？我正在尝试将这个想法转化为使用3D点。 - NaN

由于您得到的结果与此处实现不同，因此它必须是不同的a = np.array([[1, 0, 0], [1, 0, 0], [1, 0, 0]])和b = np.array([[1, 0, 0], [0, 1, 0], [0, -1, 0]])表示相同、相似和相反的方向向量配对。 - NaN

@NaN，这个实现看起来和OP在这里发布的不同。我正在遵循这个。 - Divakar

亲爱的朋友们，感谢你们的热心帮助。我编辑了原始帖子，因为我注意到在角度函数的返回语句中忘记包括np.arccos了。正如@Nan指出的那样，我已经包括了所建议的更正，网址为(https://dev59.com/g3E85IYBdhLWcg3wVR-g)，以防两个向量平行时返回Nan值。 - RRS

@RaphaelSantos 如果需要最终输出以度为单位，则可以在问题中添加说明，并使用“180/np.pi”进行缩放。除此之外看起来很不错。 - Divakar

Answer 2

如果您正在使用3D向量，可以使用工具包vg来简洁地完成此操作。它是建立在numpy之上的轻量级层，可与单个向量和向量堆栈一样有效地工作。

import numpy as np
import vg

A = np.array([[1,0,0],
              [0,1,0],
              [0,0,1]])

B = np.array([[1,0,1],
              [1,1,0],
              [0,1,0]])

vg.angle(A, B)

在我的上一个创业公司中，我创建了这个库，它的动机是为了像这样的用途：将在NumPy中冗长或不透明的简单想法变得更加简洁。

- paulmelnikow