使用针状三角形计算两个向量之间的角度

Question

我实现了一个函数 (angle_between) 来计算两个向量之间的角度。它使用针状三角形并基于 Miscalculating Area and Angles of a Needle-like Triangle and this related question.

该函数在大多数情况下似乎工作正常，除了一个我不明白发生了什么的奇怪情况：

import numpy as np
vectorA = np.array([0.008741225033460295, 1.1102230246251565e-16], dtype=np.float64)
vectorB = np.array([1, 0], dtype=np.float64)
angle_between(vectorA, vectorB)  # is np.nan

深入研究我的函数，np.nan是通过取负数的平方根产生的，负数似乎是该方法精度提高的结果：

foo = 1.0                  # np.linalg.norm(vectorA)
bar = 0.008741225033460295 # np.linalg.norm(vectorB)
baz = 0.9912587749665397   # np.linalg.norm(vectorA- vectorB)

# algebraically equivalent ... numerically not so much
order1 = baz - (foo - bar)
order2 = bar - (foo - baz)

assert order1 == 0
assert order2 == -1.3877787807814457e-17

根据 Kahan 的论文，这意味着三元组 (foo, bar, baz) 实际上并不表示三角形的边长。但是，考虑到我构建三角形的方式（请参阅代码中的注释），这应该 - 事实上 - 是这种情况。

从这里开始，我对去哪里寻找错误源感到有点迷茫。有人可以向我解释发生了什么吗？

为了完整性，这里是我的函数的完整代码：

import numpy as np
from numpy.typing import ArrayLike

def angle_between(
    vec_a: ArrayLike, vec_b: ArrayLike, *, axis: int = -1, eps=1e-10
) -> np.ndarray:
    """Computes the angle from a to b

    Notes
    -----
    Implementation is based on this post:
    https://scicomp.stackexchange.com/a/27694
    """

    vec_a = np.asarray(vec_a)[None, :]
    vec_b = np.asarray(vec_b)[None, :]

    if axis >= 0:
        axis += 1

    len_c = np.linalg.norm(vec_a - vec_b, axis=axis)
    len_a = np.linalg.norm(vec_a, axis=axis)
    len_b = np.linalg.norm(vec_b, axis=axis)

    mask = len_a >= len_b
    tmp = np.where(mask, len_a, len_b)
    np.putmask(len_b, ~mask, len_a)
    len_a = tmp

    mask = len_c > len_b
    mu = np.where(mask, len_b - (len_a - len_c), len_c - (len_a - len_b))

    numerator = ((len_a - len_b) + len_c) * mu
    denominator = (len_a + (len_b + len_c)) * ((len_a - len_c) + len_b)

    mask = denominator > eps
    angle = np.divide(numerator, denominator, where=mask)
    np.sqrt(angle, out=angle)
    np.arctan(angle, out=angle)
    angle *= 2
    np.putmask(angle, ~mask, np.pi)
    return angle[0]

编辑： 该问题肯定与 float64 有关，并且在使用较大的浮点数执行计算时消失：

import numpy as np

vectorA = np.array([0.008741225033460295, 1.1102230246251565e-16], dtype=np.float128)
vectorB = np.array([1, 0], dtype=np.float128)
assert angle_between(vectorA, vectorB) == 0

Answer 1

I just tried the case of setting vectorB as a multiple of vectorA and - interestingly - it sometimes produces nan, sometimes 0 and sometimes it fails and produces a small angle of magnitude 1e-8 ... any ideas why?

是的，我认为这就是您的问题的归结所在。这是您一直在使用的来自 the berkeley paper due to Kahan 的公式。假设a≥b、a≥c（只有这样公式才有效）和b+c≈a。如果我们暂时忽略 mu 并查看平方根下的所有其他内容，它必须都是正的，因为 a 是最长的边。 mu 是 c-(a-b) 即 0 ± a small error。如果该错误为零，您将得到零，顺便说一句。正确的结果。如果误差为负，则平方根为您提供 nan，如果误差为正，您将得到一个小角度。

请注意，当 b+c-a 非零但小于误差时，相同的参数有效。

使用针状三角形计算两个向量之间的角度

Calculating the angle between two vectors using a needle-like triangle

python

numpy

geometry

numeric

floating-accuracy