23
$\begingroup$

When applying the classical formula for the angle between two vectors:

$$\alpha = \arccos \frac{\mathbf{v_1} \cdot \mathbf{v_2}}{\|\mathbf{v_1}\| \|\mathbf{v_2}\|}$$

one finds that, for very small/acute angles, there is a loss of precision and the result is not accurate. As explained in this Stack Overflow answer, one solution is to use the arctangent instead:

$$\alpha = \arctan2 \left(\|\mathbf{v_1} \times \mathbf{v_2}\|, \mathbf{v_1} \cdot \mathbf{v_2} \right) $$

And this indeed gives better results. However, I wonder if this would give bad results for angles very close to $\pi / 2$. Is it the case? If so, is there any formula to accurately compute angles without checking for a tolerance inside an if branch?

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ This is going to depend on the implementation of the two parameter inverse tangent function. The slow, stable versions switch conditionally between working with x/y and y/x to maintain precision, while the fast ones just stick things in the right quadrant and thus aren't any more precise than the one parameter version. $\endgroup$ Commented Aug 23, 2017 at 10:51
  • $\begingroup$ You should define "loss of precision": suppose the right answer is $\alpha$ and you get instead $\alpha + \Delta$. Do you need $\Delta \ll \alpha$ or $\Delta \ll \pi$ is sufficient? $\endgroup$ Commented Aug 23, 2017 at 20:54
  • $\begingroup$ In this case, the right answer was $\alpha$ and I got $\alpha \cdot 10^8$, both $\ll 1$. $\endgroup$ Commented Aug 24, 2017 at 7:40
  • $\begingroup$ I wonder whether there are any characterizations of the imprecision scenarios, somewhere between just one example and the cosmos of working up proofs per scenario a la What Every Computer Scientist Should Know About Floating-Point Arithmetic. Looks like a topic filling curricula and lifetimes of research. $\endgroup$ Commented Jul 18, 2023 at 20:57

2 Answers 2

Reset to default
28
$\begingroup$

(I have tested this approach before, and I remember it worked correctly, but I haven't tested it specifically for this question.)

As far as I can tell, both $\|\mathbf{v}_1\times \mathbf{v}_2\|$ and $\mathbf{v}_1\cdot \mathbf{v}_2$ can suffer from catastrophic cancellation if they are almost parallel/perpendicular—atan2 can't give you good accuracy if either input is off.

Start by reformulating the problem as finding the angle of a triangle with side lengths $a=|\mathbf{v}_1|$, $b=|\mathbf{v}_2|$ and $c=|\mathbf{v}_1-\mathbf{v}_2|$ (these are all accurately computed in floating point arithmetic). There is a well-known variant of Heron's formula due to Kahan (Miscalculating Area and Angles of a Needle-like Triangle), which allows you to compute the area and angle (between $a$ and $b$) of a triangle specified by its side lengths, and do so numerically stably. Because the reduction to this subproblem is accurate as well, this approach should work for arbitrary inputs.

Quoting from that paper (see p.3), assuming $a\geq b$, $$ \mu = \begin{cases} c-(a-b),&\text{if }b\geq c\geq 0,\\ b-(a-c),&\text{if }c>b\geq 0,\\ \text{invalid triangle},&\text{otherwise} \end{cases}$$ $$ \mathrm{angle} = 2\arctan\left( \sqrt{\frac{((a-b)+c)\mu}{(a+(b+c))((a-c)+b)}} \right)$$ All the parentheses here are placed carefully, and they matter; if you find yourself taking the square root of a negative number, the input side lengths are not the side lengths of a triangle.

There is an explanation of how this works, including examples of values for which other formulas fail, in Kahan's paper. Your first formula for $\alpha$ is $C''$ on page 4.

The main reason I suggest Kahan's Heron's formula is because it makes a very nice primitive—lots of potentially tricky planar geometry questions can be reduced to finding the area/angle of an arbitrary triangle, so if you can reduce your problem to that, there is a nice stable formula for it, and there's no need to come up with something on your own.

Edit Following Stefano's comment, I made a plot of relative error for $v_1=(1,0)$, $v_2=(\cos\theta, \sin\theta)$ (code). The two lines are the relative errors for $\theta=\epsilon$ and $\theta=\pi/2-\epsilon$, $\epsilon$ going along the horizontal axis. It seems that it works. enter image description here

Improve this answer
$\endgroup$
12
  • $\begingroup$ Thanks for the link and the answer! Unfortunately the second formula I wrote does not appear in the article. On the other hand, this method can get a bit complex, as it requires projection in 2D. $\endgroup$ Commented Aug 23, 2017 at 17:18
  • 2
    $\begingroup$ @astrojuanlu There is no projection to 2d here: whatever the two 3d vectors are, they define a single (planar) triangle between them—you only need to know its side lengths. $\endgroup$ Commented Aug 23, 2017 at 17:22
  • $\begingroup$ You are right, my comment makes no sense. I was thinking in coordinates instead of lengths. Thanks again! $\endgroup$ Commented Aug 23, 2017 at 18:13
  • 3
    $\begingroup$ @astrojuanlu One more thing I want to note: it seems that there is a formal proof that the area formula is accurate in How to Compute the Area of a Triangle: a Formal Revisit, Sylvie Boldo, using Flocq. $\endgroup$ Commented Aug 23, 2017 at 21:51
  • 2
    $\begingroup$ Excellent answer, but one corner case is missing. If the vectors are actually linearly dependent and pointing in opposite directions (angle = pi/2), you will run into div/0 that needs to be handled manually (see also the linked paper). There is also the case a=b=c=0 that one may or may not wish to handle explicitly. $\endgroup$ Commented May 24, 2021 at 11:25
15
$\begingroup$

The efficient answer to this question is, not too surprisingly, in another note by Velvel Kahan:

$$\alpha=2\arctan\left(\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}+\frac{\mathbf v_2}{\|\mathbf v_2\|}\right\|,\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}-\frac{\mathbf v_2}{\|\mathbf v_2\|}\right\|\right)$$

where I use $\arctan(x,y)$ as the angle made by $(x,y)$ with the horizontal axis. (You might have to flip the order of the arguments in some languages.)

(I gave a Mathematica demonstration of Kahan's formula here.)

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ You mean $\arctan2$? $\endgroup$ Commented Sep 1, 2017 at 14:19
  • 1
    $\begingroup$ I'm used to just depicting the two-argument arctangent as $\arctan(x,y)$, yes. In a language like FORTRAN, the equivalent would be ATAN2(Y, X). $\endgroup$ Commented Sep 1, 2017 at 14:42
  • $\begingroup$ Actually, looking at your reference, it looks to me like Kahan goofed: they suggested $2·arctan( || v_1/||v_1|| – v_2/||v_2|| || / || v_1/||v_1|| + v_2/||v_2|| || )$, which will blow up if $v_1$ and $v_2$ point in opposite or nearly opposite directions. Your formula is better. $\endgroup$ Commented Sep 11, 2021 at 20:03
  • $\begingroup$ Would you expect this to be "numerically stable" in other programming languages and across commodity CPU hardware implementations? $\endgroup$ Commented Jul 19, 2023 at 8:02
  • 2
    $\begingroup$ @matanster Yes. arccos/arcsin are basically imprecise (for floating point) by design, as they require arbitrarily high relative precision. atan2 does not. $\endgroup$ Commented Jul 24, 2023 at 11:09