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

Question

I implemented a function (angle_between) to calculate the angle between two vectors. It makes use of needle-like triangles and is based on Miscalculating Area and Angles of a Needle-like Triangle and this related question.

The function appears to work fine most of the time, except for one weird case where I don't understand what is happening:

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

Digging into my function, the np.nan is produced by taking the square root of a negative number, and the negative number seems to be the result of the increased accuracy of the method:

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

According to Kahan's paper, this means that the triplet (foo, bar, baz) actually doesn't represent the side lengths of a triangle. However, this should - in fact - be the case given how I constructed the triangle (see the comments in the code).

From here, I feel a bit lost as to where to look for the source of the error. Could somebody explain to me what is happening?

---

For completeness, here is the full code of my function:

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]

Edit: The problem is definitely related to float64 and disappears when performing the computation with larger floats:

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?

Yea and I think that's what your question boils down to. Here is the formula from the berkeley paper due to Kahan that you've been using. Assuming that a≥b, a≥c (only then is the formula valid) and b+c≈a. If we ignore mu for a second and look at everything else under the square root it must all be positive since a is the longest side. And mu is c-(a-b) which is 0 ± a small error. If that error is zero you get zero which is btw. the correct result. If the error is negative the square root gives you nan and if the error is positive you get a small angle.

Notice that the same argument works when b+c-a is non zero but smaller than the error.