21

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?

J. M.
  • 3,155
  • 28
  • 37
astrojuanlu
  • 1,187
  • 2
  • 12
  • 25
  • 1
    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. – origimbo Aug 23 '17 at 10:51
  • 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? – Stefano M Aug 23 '17 at 20:54
  • In this case, the right answer was $\alpha$ and I got $\alpha \cdot 10^8$, both $\ll 1$. – astrojuanlu Aug 24 '17 at 07:40
  • 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. – matanox Jul 18 '23 at 20:57

2 Answers2

25

(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

Kirill
  • 11,438
  • 2
  • 27
  • 51
  • 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. – astrojuanlu Aug 23 '17 at 17:18
  • 2
    @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. – Kirill Aug 23 '17 at 17:22
  • You are right, my comment makes no sense. I was thinking in coordinates instead of lengths. Thanks again! – astrojuanlu Aug 23 '17 at 18:13
  • 2
    @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. – Kirill Aug 23 '17 at 21:51
  • Excellent answer, but I dispute that you can always accurately compute $c$ in floating point arithmetic. In fact if $c < \epsilon \cdot \min(a,b)$ then catastrophic cancellations occurs in computing the components of $(\mathbf{v}_1-\mathbf{v}_2)$. – Stefano M Aug 23 '17 at 21:51
  • @StefanoM I accept that I didn't think carefully enough about that to be certain, but it seems to me your example doesn't quite work: if $v_1$ and $v_2$ are so close, then component-wise subtraction would be exact (in fact: $x-y$ is computed exactly when $x$ and $y$ are at most a factor of two from each other). This only works since $v_{1,2}$ are input data and therefore exact. I tried doing the usual relative error analysis on $|v_1-v_2|^2$ with epsilons, and I got the relative error as $O(\epsilon)$. Do you have a specific counterexample in mind? So far, I think my claim was right. – Kirill Aug 23 '17 at 22:07
  • @StefanoM See edit, I still think it's okay. What do you think? – Kirill Aug 23 '17 at 22:40
  • My wording was imprecise, sorry. In fact I was not assuming that $v_{1,2}$ are exact, but thinking about the ill-conditioning problem, with respect to a change of base. E.g: $v_1 = (1,0,0)$ and $v_2 = (1, \epsilon, 0)$ is not problematic, but if we rotate by $\pi/2$ then we have $v_1 = (\sqrt{2}/2,\sqrt{2}/2,0)$, and $v_2 = ((1-\epsilon)\sqrt{2}/2, (1+\epsilon)\sqrt{2}/2,0)$. But this is a different problem, so I accept your claim. – Stefano M Aug 23 '17 at 22:41
  • @StefanoM Then I agree with what you said. – Kirill Aug 23 '17 at 22:55
  • In my comment above, the rotation is $\pi/4$, of course... (but it is not possible to edit comments, sorry again for all this fuss.) – Stefano M Aug 24 '17 at 06:37
  • 1
    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. – FirefoxMetzger May 24 '21 at 11:25
  • For the rest of us whom haven't taken a related math course on this topic, what does it rigorously mean to say that ... something is done numerically stably, as mentioned in this answer with regard to the adaptation of Heron's Formula? I don't think we have that in the referenced paper either, at least not explicitly. – matanox Jul 18 '23 at 21:07
14

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.)

J. M.
  • 3,155
  • 28
  • 37
  • 1
    You mean $\arctan2$? – astrojuanlu Sep 01 '17 at 14:19
  • 1
    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). – J. M. Sep 01 '17 at 14:42
  • 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. – Don Hatch Sep 11 '21 at 20:03
  • Would you expect this to be "numerically stable" in other programming languages and across commodity CPU hardware implementations? – matanox Jul 19 '23 at 08:02
  • Numpy adaptation added here – matanox Jul 19 '23 at 14:24
  • @matanster Yes. arccos/arcsin are basically imprecise (for floating point) by design, as they require arbitrarily high relative precision. atan2 does not. – Sneftel Jul 24 '23 at 11:09