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

Question

我实现了一个函数 (angle_between) 来计算两个向量之间的角度。它利用针状三角形并基于 Miscalculating Area and Angles of a Needle-like Triangle和 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 有关，并且在使用较大的 float 执行计算时消失:

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 ± 一个小错误。如果该错误为零，您将得到零，顺便说一句。正确的结果。如果误差为负，则平方根给出 nan，如果误差为正，则得到一个小角度。

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