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

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) 实际上并不代表三角形的边长。然而，考虑到我是如何构造三角形的，这应该——事实上——是这种情况（参见代码中的 cmets）。

从这里开始，我有点迷失在哪里寻找错误的根源。有人可以向我解释发生了什么吗？

---

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

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

我刚刚尝试了将 vectorB 设置为 vectorA 的倍数的情况 - 有趣的是 - 它有时会产生 nan，有时会产生 0，有时会失败并产生 1e-8 大小的小角度......有什么想法吗？

是的，我认为这就是您的问题归结为的原因。这是您一直在使用的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 不为零但小于错误时，相同的参数有效。