Compute $\sqrt{a}$ with Newton's method.

2.6k Views Asked by At

It's possible to compute $\sqrt{a}$, $a>0$, using only addition and division and to compute $1/b$, $b>0$ through addition and multiplication by solving the equations

i) $x^2-a=0$

ii) $\frac{1}{x}-b=0$

using Newton's method. Write down the iteration formulas. What is a good starting value $x_0?$

The iteration formula is $$x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}, \quad k=1,2,\ldots$$

How do I a choose a good starting value in the case of $f(x)=x^2-a$?

6

There are 6 best solutions below

4
On

You can just use $a$ as a starting value. Newton's algorithm works fine from there. A few iterations, and you will be right on track to finding the square root of $a$.

If you want a "smart" starting valuse, then you can use the fact that square roots basically halve the number of digits. So write $a$ on the form $x\cdot 100^n$, where $1\leq x<100$ and $n$ is an integer. Then you can have as a starting value anything close to $\sqrt x\cdot 10^n$. For instance, $\lfloor\sqrt x\rfloor\cdot10^n$. This can of course be done in any base, like binary, in place of base ten if you so wish.

Picking a good initial value can be important. For instance, the fast inverse square root was invented to quickly get a good initial value for the function $f(x) = \frac1{\sqrt x}$. It exploited the fact that this was to be done on a computer, on numbers stored as 32-bit floating points, to exploit the exact way comnputers store and interpret numbers in order to get a final result about four times quicker than previous methods (directly computing the square root and inverting it). It consists of basically what I've done here: halving (and flipping the sign of) the exponent, and get an approximation to $1/\sqrt x$ (which they did by subtracting a magic number).

2
On

You'll find the point at which $a=x^2$ by root approximation in this method. It doesn't necessarily matter what your $a$ is. $$x_{k+1}=x_k-\frac{(x_k)^2-a}{2x_k}$$ $$\rightarrow x_{k+1}=\frac{(x_k)^2+a}{2x_k}$$ Inserting a value for $a$ of your choice will cause this iteration to approximate its square root, thereby solving $x^2-a=0$.

I've done loads of work on this iteration,so I can present this to you: https://www.desmos.com/calculator/qdbskhvmjg (the first six iterations graphed)

They also solve the following sets of equations. $$\sqrt{\mu_k}+\sqrt{a+\mu_k}=x_k$$ $$\sqrt{a+\mu_k}=x_{k+1}$$ $$\sqrt{\mu_{k+1}}+\sqrt{a+\mu_{k+1}}=x_{k+1}$$


Applying a similar thing to (ii) yields $x_{k+1}=2x_k-b(x_k)^2$, which upon testing with some $b$ values, rapidly accelerates to $-\infty$ for $b>1$, but for $b<1$ it iterates very well.

2
On

You don't tell us the data type of $a$.

If the floating-point representation is available, it suffices to extract the exponent, let $e$, and use the starting value $2^{\lfloor e/2\rfloor}$. If the exponent is odd, $\sqrt2\,2^{\lfloor e/2\rfloor}$ is even better.

If $a$ is an integer, you can perform successive multiplications of a variable by $4$ and simultaneously of another by $2$, both starting from $1$, until the value $a$ is exceeded. Then $4^k\approx a$ implies $2^k\approx\sqrt a$.

If $a$ is a real number with unspecified representation and is smaller than $1$, perform the above procedure with the powers of $\dfrac14$ and $\dfrac12$.


For the inverse, the initial value can be $2^{-e}$. You can also find a $2^{\pm k}$ approximation and turn it to $2^{\mp k}$.


It is possible to be a little faster, by using exponential search rather than linear (hence successive squarings). This is useful only when large/small orders of magnitude are expected.

12
On

In the case of finding $\sqrt{a}$, it turns out that any positive $x_0$ will work. However, very large $x_0$ will spend a long time basically just dividing by $2$ over and over again until getting close to the solution, while very small $x_0$ will immediately be sent to very large values and then go through the same issue. You can see this by looking at the Newton formula carefully: if $x_n \ll \sqrt{a}$ then $x_{n+1} \approx \frac{a}{2x_n} \gg \sqrt{a}$. Similarly if $x_n \gg \sqrt{a}$ then $x_{n+1} \approx \frac{x_n}{2} \gg \sqrt{a}$. As a result, Newton's method with a bad initial guess is basically the same as bisection using the bracket $\{ x_0,a/x_0 \}$ (in fact it is slightly worse than that), until a certain level of accuracy is achieved. At that point convergence accelerates.

$x_0=a$ and $x_0=1$ are both decent compromises that take no real work to compute up front, but based on the thinking above, you will still spend around $|\log_2(a)|/2$ time just dividing by $2$ a whole bunch of times with those, which could be a while if $a$ is large or small.

If $a$ is given to you in the form $a=x \cdot 2^n$ with $x \in [1,2)$ and $n$ an integer, then you can "cheat" by taking something like $x_0=2^{n/2}$, which in some sense "does all the divisions by $2$ at once". (To deal with the case where $n$ is odd, you would need to either precompute $\sqrt{2}$ or just pick one way or the other to round.) You could do this in any other base $b$, e.g. $10$, as well. But if you have to compute this $n$ (which is basically the first digit of $\log_b(a)$), then you're spending more time doing that then it would've taken to just do the problem with $x_0=1$.

For $1/b$ you need to be more careful, because not any initial guess will work.

2
On

Newton's Method $$x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}, \quad k=1,2,....$$

works fine with an arbitrary initial value $ x_0$ as long as the $f'(x_k)$ stays away from zero.

For example if we want to solve $\sin x =0$ with Newton's method and we are looking for $x=\pi$ we better stay away from $x=\pi/2$ for an initial value. Otherwise we get into a chaotic pattern.

In your case of $ f(x)= x^2-a$ , the initial value could be $a$ itself.

0
On

The goal is to compute $\sqrt{a}$ with a small relative error using a small number of iterations.

I will now analyse Newton's method for the equation $$ f(x) = x^2 - a = 0$$ and explain how to construct an excellent initial guess from scratch. I conclude with some general observations.

Let $$e_n = \sqrt{a} - x_n$$ denote the error and let $$r_n = \frac{e_n}{\sqrt{a}}$$ denote the relative error.

Newton's method takes the form $$ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} = x_n - \frac{x_n^2 - a}{2x_n} = \frac{x_n^2 + a}{2x_n} = g(x_n) $$ where $g : (0, \infty) \rightarrow \mathbb{R}$ is given by $$ g(x) = \frac{1}{2} \left( x + \frac{a}{x} \right). $$ Any sensible choice of $x_0$ must satisfy $x_0 > 0$ and it is clear that $x_0 > 0$ implies that $x_n > 0$ for all $n$. We have $$ e_{n+1} = \sqrt{a} - x_{n+1} = - \frac{x_n^2 - 2 a x_n + a}{2x_n} = - \frac{(x_n - \sqrt{a})^2}{2x_n} = - \frac{e_n^2}{2 x_n} $$ which implies $$ r_{n+1} = - \frac{\sqrt{a}}{2 x_n} \frac{r_n^2}{2}. $$ An elementary analysis reveals that $$ g(x) \ge \sqrt{a}, \quad x > 0. $$ It follows that $$ x_n \ge \sqrt{a}, \quad n \ge 1, $$ regardless of the choice of $x_0 > 0$. We can therefore assume that $x_n \ge \sqrt{a}$. If necessary, we simply discard $x_0$ and renumber, i.e., $x_n' = x_{n+1}$ for $n \ge 0$. It follows that $$ |r_{n+1}| \leq \frac{1}{2} r_n^2. $$ Let $y_n = \frac{1}{2} |r_n|$, then $$ y_{n+1} \leq y_n^2 $$ By induction on $n$ we discover, that $$ y_n \leq \left(y_0\right)^{2^n} $$ or equivalently $$ |r_n| \leq 2 \left(\frac{|r_0|}{2} \right)^{2^n}. $$ We conclude that Newton's method converges quadratically, if the initial relative error is less than $2$.

We now turn to the question of constructing a good initial guess. This is not easy for a general $a > 0$. However, the general case can be reduced to the special case of $a \in [1,4)$. To this end, we consider the binary representation of $a > 0$. We have $$ a = f \cdot 2^m,$$ where $f \in [1,2)$ and $m$ is an integer. Let $m = 2k + r$, where $r \in \{0,1\}$. Then, $$ \sqrt{a} = (\sqrt{f \cdot 2^r}) 2^k.$$ We conclude that any square root can be computed provided that we have the ability to compute the square root of any number in the interval $[1,4]$.

An initial guess can be constructed form the best uniform approximation of the square root on this interval. It is easy to verify that the function $h$ given by $$ h(s) = \sqrt{s} - \left(\frac{1}{3} s + \frac{17}{24}\right) $$ satisfies $$ \forall \: s \in [1,4] \: : \quad -\frac{1}{24} \leq h(s) \leq \frac{1}{24}. $$ The initial guess $$ x_0(a) = \frac{1}{3} a + \frac{17}{24} $$ will ensure that the initial residual satisfies $$ |r_0| \leq \frac{1}{24}. $$ In this case, we find that $n=4$ is the smallest integer such that $$ |r_n| \leq u, $$ where $u = 2^{-53} \approx 10^{-16}$ is the unit roundoff in IEEE double precision floating point arithmetic.


Be wary of the difference between exact and computer arithmetic. In exact arithmetic, Newton's method applied to the equation $$ f(x) = x^2 - a = 0$$ will converge for any initial guess $x_0 > 0$. This is not true in floating point arithmetic. The product $x \cdot x$ needed for the computation of $f(x)$ will overflow if $$x > \sqrt{\Omega},$$ where $\Omega$ is the largest representable number. The initial guess must therefore satisfy $$x_0 \leq \sqrt{\Omega}$$ or your implementation will fail immediately. This significantly reduces the choice of $x_0$. In particular, $x_0 = a$ is not a good choice for large values of $a$.

It is not irrelevant how you implement Newton's method. The expressions $$ x_{n+1} = x_n - \frac{x_n^2 - a}{2x_n}$$ and $$ x_{n+1} = \frac{x_n^2 + a}{2x_n}$$ and $$ x_{n+1} = \frac{1}{2} \left( x_n + \frac{a}{x_n} \right)$$ are equivalent in exact arithmetic, but not in floating point arithmetic. Here the first expression is of the preferred form, i.e., a good approximation plus a small correction. It does not matter that $x_n^2 - a$ suffers from catastrophic cancellation when $x_n \approx \sqrt{a}$. If $x_n$ is a good approximation, then you do not need to compute the correction with a tiny relative error. The first expression also has the virtue of utilizing the residual $x_n^2 - a$. This number is often used to determine if the iteration has converged.


Even in exact arithmetic it is not true that avoiding points where $f'(x_0) \not = 0$ is sufficient to ensure convergence. An example is the equation $$f(x) = \sin(x) = 0, \quad |x| < \frac{1}{2} \pi.$$ On this interval $f$ has one zero namely $x = 0$ and $f'(x) = \cos(x)$ has no zeros. Newton's method takes the form $$x_{n+1} = x_n - \tan(x_n).$$ Now if $x_0 = \nu$, where $\nu$ solves the equation $$ \tan(x) = 2x $$ then $$ x_1 = - \nu, \quad x_2 = \nu$$ and we are back where we started. It is easy to verify that $\nu \approx 1.1656$ is one such choice.