使用类似指针的三角形计算两个向量之间的夹角

Question

使用类似指针的三角形计算两个向量之间的夹角

pythonnumpygeometrynumericfloating-accuracy

3

我实现了一个函数(angle_between)来计算两个向量之间的夹角。它利用类似于针尖的三角形，基于“针尖三角形的面积和角度计算错误”以及“计算向量之间角度的数值稳定方法”的相关问题。

这个函数在大多数情况下运行良好，但有一个奇怪的情况我不理解：

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

- FirefoxMetzger

请注意，如果三边满足强三角不等式，即两个较短的边之和必须严格大于较长的边，则它们组成一个三角形。但是对于您来说，情况并非如此，因为 bar + baz == 1 == foo。 - Lukas S

@user2640045 我猜 bar + baz == 1 == foo 是由于浮点数不精确造成的？三个向量 vectorA、vectorB 和 vectorA - vectorB 的长度应该始终形成一个有效的三角形，对吧？除此之外，函数应正确处理两种退化情况：vectorA == vectorB 和 vectorA == -vectorB。前者通过 len_c 为0来处理，后者通过 np.putmask(angle, ~mask, np.pi) 来处理。 - FirefoxMetzger

不，还有一种情况，即当vectorA和vectorB是彼此的倍数时。这里几乎就是这种情况。如果我用零替换1.1102230246251565e-16，它们就会是倍数关系。我想1.1102230246251565e-16与零之间的差异不足以避免问题。 - Lukas S

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

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Lukas S · Accepted Answer

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

是的，我认为这就是你的问题的核心。这是由于Kahan所写的伯克利论文中使用的公式，你一直在使用它。

假设a≥b，a≥c（只有这样公式才有效），并且b+c≈a。如果我们暂时忽略mu，并查看根号下的所有其他内容，它必须都是正数，因为a是最长的边。而mu是c-(a-b)，即0±一个小误差。如果该误差为零，则得到零，这是正确的结果。如果误差为负，则平方根会给出nan，如果误差为正，则会得到一个小角度。

请注意，当b+c-a不为零但小于误差时，相同的论证也适用。