3
\$\begingroup\$

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)
\$\endgroup\$
2
\$\begingroup\$

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.

\$\endgroup\$
3
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ What other functions could the newton_raphson method calculate too? \$\endgroup\$ – Caridorc Dec 30 '14 at 22:03
  • 1
    \$\begingroup\$ Any, as long as a derivative is correct. \$\endgroup\$ – vnp Dec 30 '14 at 22:15
  • \$\begingroup\$ 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. \$\endgroup\$ – Caridorc Dec 30 '14 at 22:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.