Newton's square root

Question

This is the code I wrote to compute the newton's square root of a number.

I aimed for extreme clarity over efficiency and I have put an abundant docstring at the start of the newton_sqrt function where I describe the main formula and explain all the parameters.

This function is not meant to be fast (standard library is 4x times faster and more accurate), I wrote just to learn some good math coding techniques in Python.

• This was written in Python 3 but is fully compatible with Python 2.
• I created constants for the error messages. Is that a good practice? Or should I just write the errors down verbatim when need occurs?
• I chose n / 2 as my initial guess. Do you know a better initial guess?
• Are all the option parameters that I am using necessary or should I delete some of them?

from __future__ import division
import time
def newton_sqrt(n,x0=None,tolerance=10 ** (-3),max_iter=15,epsilon=10 ** (-5)):
 """
 This function computes the sqare root of a number
 using Newton's method.
 After chosing a starting approximate value for x runs:
 x = x - f(x) / derivative_of_f(x)
 After a given number of iterations or if the result is good enough
 the programme will stop and return x.
 Parameters:
 Compulsory:
 @ n: the number of which the square root must be taken.
 Optional:
 @ x0: the initial approximation of n.
 @ tolerance: the maximum error permitted.
 @ max_iter: the maximum number of iterations.
 @ epsilon: a number so small that dividing by it would break
 @ the programme, giving vastly inapproximate results.
 """
 DIVISION_BY_TINY_DERIVATIVE_ERROR = """The derivative of f is {} and therefore under the {} limit chosen.
 Dividing by it would give a result with a vast error."""
 NOT_ENOUGH_ACCURACY_IN_ITER_ERROR = """Solution of {} accuracy for {} was not found in {} iterations"""
 if x0 is None:
 x0 = n / 2
 f = lambda x: x ** 2 - n
 fprime = lambda x: 2 * x # The derivative of f(x)
 for i in range(max_iter):
 y = f(x0)
 yprime = fprime(x0)
 if(abs(yprime) < epsilon):
 raise Exception(DIVISION_BY_TINY_DERIVATIVE_ERROR.format(
 yprime, epsilon))
 x1 = x0 - y / yprime
 if(abs(x1 - x0) / abs(x1) < tolerance):
 return x1
 x0 = x1
 raise Exception(NOT_ENOUGH_ACCURACY_IN_ITER_ERROR.format(
 tolerance,n,max_iter))
def test_and_timing(max_):
 start = time.time()
 n_sqrt = [newton_sqrt(i) for i in range(1, max_)]
 print("Newton's sqrt took {} seconds to execute for i 1:{}".format(
 time.time() - start, max_))
 start = time.time()
 b_sqrt = [i ** 0.5 for i in range(1, max_)]
 print("Built-in sqrt took {} seconds to execute for i 1:{}".format(
 time.time() - start, max_))
 print("""Sum of all newton sqrts is: {}
 Sum of all built-in sqrts is: {}""".format(sum(n_sqrt), sum(b_sqrt)))
if __name__ == "__main__":
 test_and_timing(100000)

Answer 1

The spacing in the function definition is nonstandard; it should be one of:

def newton_sqrt(n, x0=None, tolerance=10 ** -3, max_iter=15, epsilon=10 ** -5):
def newton_sqrt(n, x0=None, tolerance=10**-3, max_iter=15, epsilon=10**-5):
def newton_sqrt(n, x0=None, tolerance=10e-3, max_iter=15, epsilon=10e-5):

With my preference for the 3rd. (FWIW, the second technically breaks PEP 8.)

For some reason, all of your line-wrapped strings are indented 3 spaces. It should be 4, like the rest of the code.

DIVISION_BY_TINY_DERIVATIVE_ERROR would probably be better without a line break and indent; you can avoid this as so:

NOT_ENOUGH_ACCURACY_IN_ITER_ERROR doesn't need tripple quotes.

IMO, both of the above shouldn't be called _ERROR if they are strings. I would also drop _IN_ITER.

Don't write

f = lambda x: x ** 2 - n
fprime = lambda x: 2 * x # The derivative of f(x)

Instead, write

def f(x): return x ** 2 - n
def fprime(x): return 2 * x # The derivative of f(x)

Your

if(abs(yprime) < epsilon):

doesn't need brackets:

if abs(yprime) < epsilon:

You should almost never raise Exception; instead raise appropriate specific exceptions. In this case I would use ValueError.

Don't line break onto the indentation; add a hanging indent.

This code is overcomplicated:

x1 = x0 - y / yprime
if(abs(x1 - x0) / abs(x1) < tolerance):
 return x1
x0 = x1

You can just do

diff = y / yprime
x0 -= diff
if abs(diff / x0) < tolerance:
 return x0

Your test_and_timing could utilize better naming and the last print should be split into two. I would also reduce the number of decimal places printed:

def test_and_timing(up_to):
 start = time.time()
 newton_roots = [newton_sqrt(i) for i in range(1, up_to)]
 took = time.time() - start
 print("Newton's sqrt took {:.3f} seconds to execute for i in range(1, {})".format(took, up_to))
 start = time.time()
 builtin_roots = [i ** 0.5 for i in range(1, up_to)]
 took = time.time() - start
 print("Built-in sqrt took {:.3f} seconds to execute for i range(1, {})".format(took, up_to))
 print("Sum of all newton sqrts is: {:.2f}".format(sum(newton_roots)))
 print("Sum of all built-in sqrts is: {:.2f}".format(sum(builtin_roots)))

Your timing is also slightly worsened by the cost of a list comprehension. Removing that would be more honest:

def test_and_timing(up_to):
 start = time.time()
 for i in range(1, up_to):
 newton_sqrt(i)
 took = time.time() - start
 print("Newton's sqrt took {:.3f} seconds to execute for i in range(1, {})".format(took, up_to))
 start = time.time()
 for i in range(1, up_to):
 i ** 0.5
 took = time.time() - start
 print("Built-in sqrt took {:.3f} seconds to execute for i range(1, {})".format(took, up_to))
 newton_roots = (newton_sqrt(i) for i in range(1, up_to))
 builtin_roots = (i ** 0.5 for i in range(1, up_to))
 print("Sum of all newton sqrts is: {:.2f}".format(sum(newton_roots)))
 print("Sum of all built-in sqrts is: {:.2f}".format(sum(builtin_roots)))

Note that f(x) / fprime(x) is actually cheaper to calculate at once, simplifying to just 0.5 * (x - n / x). Further, note that we can just compare x to epsilon / 2 instead of fprime to epsilon; given the approximation of epsilon in the input this is basically the same as comparing x to epsilon.

Now, this might be problematic if you want to generalize, but there's actually little reason to have a DIVISION_BY_TINY_DERIVATIVE_ERROR. If the number doesn't stabalize, you'll just get a NOT_ENOUGH_ACCURACY_IN_ITER_ERROR.

That said, I'd actually avoid generalizing too much; a good square-root algorithm is not in this case a good general algorithm. It makes sense to specialize this use-case because you get both speed and accuracy improvements from doing so. A general algorithm should also support several cases you currently don't; look at how scipy.optimize.newton works for useful ideas along that line.

Answer 2

Looks good overall. Some suggestions:

• The newton_sqrt is written in a quite generic manner. Unfortunately it may only calculate square roots - because lambdas are defined inside the function. I would rewrite it along the lines of

  def newton_raphson(x, x0, func, derivative, governing_options):
   your code verbatim, except the lambda definitions 
  def newton_sqrt(x, governing_options):
   return newton_raphson(x, x/2, lambda y: y*y - x, lambda y: 2*y, governing_options)
  

• In any case, you most likely want to test for n > 0 and raise InvalidArgiment accordingly.

• I am not sure about the epsilon. If an abs(prime) < epsilon condition triggers, the iteration would not converge anyway. I don't see an extra value here, but maybe it is just me.

• Caridorc: What other functions could the newton_raphson method calculate too?
• vnp: Any, as long as a derivative is correct.
• Caridorc: Mathematics always menages to amaze me: I followed an algorithm that I though would only calculate square roots and I discover through a comment that it can calculate anything... it is like building something like a car thinking that it can only go on the ground and then discovering that it can fly, swim, and dig too... The beauty of a formula depends on how much general it is, given how general this algorithm is, it is very beatiful. Thank you for giving me this perl of knowledge I will investigate about it soon.