Newton's method and trig functions on a computer

I'm trying to use Newton's method to find roots for the function $A \cos(\Theta_2 - \Theta_1) + B \sin(\Theta_1)$. (That is, iterate $x_{i+1} = x_i - f(x_i) / f'(x_i)$). I've got a working implementation and everything's great, except that I'm not getting all the precision I'd expect.

I'm using 64 bit floating point numbers (doubles) throughout, and I'm expecting effective precision to be at least $1e^{-10}$, but I'm getting only maybe $1e^{-9}$. I traced the problem to the fact that $\cos(\frac{\Pi}{2}) \approx 1e^{-9}$ on my machine. Which means when I evaluate my function it returns 0 at a place that isn't exactly a root of my function.

This isn't entirely unexpected. I'm well aware of numeric issues on machines. I think most modern CPUs use CORDIC for trig evaluations, which interpolate between table lookup values. Some numeric noise on the least significant bits is sort of par for the course.

I'm just wondering if there's some known method to get around the limitations of the machine implementation of sine and cosine for root finding with Newton's Method? Is there some way to maybe transform my initial equation in some way that will give me a function with the same roots in some given neighborhood but better numerical behavior? Taylor series approximations come to mind but they don't always have great numerical behavior (although I don't have a lot of experience in that regard so I might be wrong).

UPDATE

I was slightly wrong. The real problem is that: $\cos(3.141592650897981) = -1$ on my machine. That angle value is roughly $2.69e^{-9}$ off from actual $\Pi$. You can reproduce this with google, actually. Try googling "cos(3.141592650897981)" and "PI - 3.141592650897981" to see what I mean.

-
Do you mean $\cos(\pi/2) \approx 1e-9$? $\sin(\pi/2)$ should be $1$. – Robert Israel Jan 3 '12 at 22:26
Er, indeed I do :P That would be quite a numerical error! I've fixed my post above – Jay Lemmon Jan 3 '12 at 22:44
@JayLemmon : I might write $x_{i+1}=x_i-f(x_i)/f'(x_i)$ or $t_{i+1}=t_i-f(t_i)/f'(t_i)$ or maybe something like $s=t-f(t)/f'(t)$ or $w=x-f(x)/f'(x)$, but the way you wrote it, "$x_{i+1}=x_i-f(t)/f'(t)$" leaves some things to be desired. – Michael Hardy Jan 3 '12 at 22:55
@Michael Hardy: I see your point. I've updated my post. – Jay Lemmon Jan 3 '12 at 23:04

Somewhere your program seems to be using single precision rather than double precision. Any reasonable implementation of double-precision arithmetic should give you $\cos(\pi/2)$ with an error of $10^{-16}$ or so. What do you get for $\pi$?

-
Ah, you were right to be suspicious. The actual problem is that I'm finding the value of sine/cosine from $\Theta = \frac{-\Pi}{2} + \frac{\Pi}{2} * 2.9999999983$. The calculated cosine of that is apparently -1 exactly, but the sine of that is 0.0000000026918122290084216. I'll update my original post – Jay Lemmon Jan 3 '12 at 22:56

Okay, after quite a bit of research and work I think I've found a robust fix for the problem (fixes all the test cases that I had failing before). It's non trivial enough that I wrote up the explanation and hopefully it'll be useful to someone someday.

Basically, cosine and sine lose digits of precision as you get near to $1$ (because $1 - 1e^{-20} = 1$, even in double precision, because the mantissa can't contain all the digits). So you need to separate out the $1$ from the digits of precision and treat them separately. That is, get the trig call in to a form like $1 - \epsilon$, and then calculate $\epsilon$ directly. So to fix, it, I transform the trig calls thusly:

• $\cos(\Theta) = 1 - 2 \sin^2\frac{\Theta}{2}$ when $\cos(\Theta)$ is near $1$
• $\cos(\Theta) = -1 + 2 \cos^2\frac{\Theta}{2}$ when $\cos(\Theta)$ is near $-1$

For sine, I combined the two identities above with the identity $\sin(\Theta) = \cos(\frac{\Pi}{2} - \Theta)$. So basically just plugged in $\frac{\Pi}{2} - \Theta$ for $\Theta$ in the above equations. These transformations are exact for any values of $\Theta$, but they lose digits when they're away from $1$ or $-1$ (same as the normal trig calls do when they're near $1$ or $-1$).

The next problem is that, since the aggregate function evaluates to approximately $0$ in the neighborhood of its root, there's a lot of cancellation issues at work when you actually evaluate the functions. So instead of just naively summing up all the terms, I use something like Knuth's summation scheme (algorithm 4.1 (using algorithm 3.1) in this paper: Accurate Sum and Dot Product) (basically keeping a running sum and a running underflow value, and taking the final sum as the sum of the running sum and the underflow at the end).

The final function evaluation seems accurate to something like at least $1e^{-20}$ or so near the root (that is, basically the full double precision), which is sufficient for Newton's method to actually converge to an answer with a good number of significant digits (I stop the iteration at a step size of $1e^{-12}$, but it could probably converge even further).

Anyway, hope that's useful to someone in the future.

-
"Compensated summation" is indeed a good thing to do when you're confronted with summing terms of different signs and magnitudes. As for cosine behaving badly when it's near $1$ or $-1$, this was precisely the reason why the versine was used a fair bit in the olden days. – J. M. Jan 7 '12 at 2:28
I actually came across versine during my research. Never heard of it before, but various versine identities helped me figure some of the above out. – Jay Lemmon Jan 7 '12 at 6:13

With your update at least the problem is clear : $\cos(\pi-\epsilon)=-\cos(\epsilon) \approx -1+\epsilon^2/2$.
$\epsilon^2/2 \approx 3.6\ 10^{-18}$ (that is the 18 digits you may hope using double precision) but $\epsilon \approx 2.69\ 10^{-9}$.

So that your program is probably correct.
Perhaps that some rewriting like $x_{i+1} = (x_i f'(x_i) - f(x_i)) / f'(x_i)$ could help... (quad. precision surely with the adapted mathematical libraries!).
EDIT : This rewriting isn't useful in this case so let's forget that and resume the situation.

The problem is the evaluation of $\cos$ or $\sin$ function at points were the derivative is $0$ so that no first order term will appear after evaluation but only terms of order two or more. You can't 'cure' this problem inside the cos or sin function because of the constant term : $-1 + \frac{\epsilon^2}{2}$ is considered $-1$ for $\epsilon$ of order $10^{-8}$ or less (IEEE 64 bits double).

If you need more precision you'll have to combine two or more terms : in your specific example you may rewrite $1+\cos(u)$ as $\ 2\cos(u/2)^2$ (or directly as $(u-\pi)^2/2!-(u-\pi)^4/4!$) when $|u-\pi|\lt 10^{-6}$ (for example). Note that there should be no loss of significance for $u$ near $0$.

In the more general case $A \cos(\Theta2-\Theta1) + B \sin(\Theta1)$ you could perhaps proceed at some rewriting when : ( this sum is very small ) and ( if $\Theta2-\Theta1$ is 'near' $k \pi$ or if $\Theta1$ is 'near' $(k+1/2) \pi$ ). Another technical solution could be to use higher precision (slower!) arithmetic in the few cases where the sum is small!

The real problem is to use the Newton(–Raphson) method outside of its applicability limits. Quadratic convergence is not warranted when $f'(s)=0$ for the solution $s$. In your case $f(s)$ with $s=\pi+\epsilon$ is of order $\epsilon^2$ and $f'(s)$ is of order $\epsilon$ but $\epsilon^2 \lt 10^{-16}$ was replaced by $0$ so that instead of adding a term of order $\epsilon=\frac{\epsilon^2}{\epsilon}$ at each iteration we add merely $0$ and the Newton algorithm will not converge. Another solution could be to try an alternative algorithm when $f(s)$ and $f'(s)$ become too small...

-
I haven't seen that formulation of Newton's method before... It would at first glance seem to handle the cases where $f(x_i)$ converges before its derivative does. But you're also changing the order of adds and divides, and I'm wondering if that affects the convergence/numerical properties at all? Any links you can point me to? – Jay Lemmon Jan 4 '12 at 0:10
I tried it out (writing some C code for this) but the error remained of the same order ($10^{-16}$ with the constants I used) so... Could you provide the $A, B, \theta2$ producing your worst result? – Raymond Manzoni Jan 4 '12 at 0:30
$A$,$B$ in my test case are both 1. $\Theta_1$ is 0, and $\Theta_2(t) = -\frac{\Pi}{2} + \frac{\Pi}{2} * t$, and there's a constant term I didn't include in my original problem. So the actual equation I'm trying to find roots for in this specific failing test is $1 + \cos(-\frac{\Pi}{2} + \frac{\Pi}{2} * t)$ (specifically the root at 3.0) – Jay Lemmon Jan 4 '12 at 0:55
I reproduced the problem (I got an error of 2.2 E-10 or 1.8 E-10 for the root at 3.0 on a PowerPC processor, not much better...). Well the cosinus near $k \pi$ can't be made more precise so that some rewriting may be necessary for specific cases like $1+\cos(u)=2\cos(u)^2/2$ (by the way $\sin(\theta1)=0$). I fear that for $A=B$ some solutions will remain less precise as you observed... – Raymond Manzoni Jan 4 '12 at 1:45
Thanks for your explanation. I at least understand the problem pretty well now. I'll play around with various hacks and see if I can get reasonable behavior. – Jay Lemmon Jan 4 '12 at 18:57