3

I'm trying to implement Newton's method for fun in Python but I'm having a problem conceptually understanding the placement of the check.

A refresher on Newton's method as a way to approximate roots via repeated linear approximation through differentiation:

Newton's Method

I have the following code:

# x_1 = x_0 - (f(x_0)/f'(x_0))
# x_n+1 - x_n = precision

def newton_method(f, f_p, prec=0.01):
    x=1
    x_p=1
    tmp=0

    while(True):
        tmp = x
        x = x_p - (f(x_p)/float(f_p(x_p)))
        if (abs(x-x_p) < prec):
            break;
        x_p = tmp

    return x

This works, however if I move the if statement in the loop to after the x_p = tmp line, the function ceases to work as expected. Like so:

# x_1 = x_0 - (f(x_0)/f'(x_0))
# x_n+1 - x_n = precision

def newton_method(f, f_p, prec=0.01):
    x=1
    x_p=1
    tmp=0

    while(True):
        tmp = x
        x = x_p - (f(x_p)/float(f_p(x_p)))
        x_p = tmp
        if (abs(x-x_p) < prec):
            break;

    return x

To clarify, function v1 (the first piece of code) works as expected, function v2 (the second) does not.

Why is this the case?
Isn't the original version essentially checking the current x versus the x from 2 assignments back, rather than the immediately previous x ?

Here is the test code I am using:

def f(x):
    return x*x - 5

def f_p(x):
    return 2*x

newton_method(f,f_p)

EDIT

I ended up using this version of the code, which forgoes the tmp variable and is much clearer for me, conceptually:

# x_1 = x_0 - (f(x_0)/f'(x_0))
# x_n+1 - x_n = precision

def newton_method(f, f_p, prec=0.01):
    x=1
    x_p=1
    tmp=0

    while(True):
        x = x_p - (f(x_p)/float(f_p(x_p)))
        if (abs(x-x_p) < prec):
            break;
        x_p = x

    return x
Aristides
  • 3,805
  • 7
  • 34
  • 41

2 Answers2

2

Let x[i] be the new value to be computed in an iteration.

What is happening in version 1:

The statement x = x_p - (f(x_p)/float(f_p(x_p))) translates to:

x[i] = x[i-2] - f(x[i-2])/f'(x[i-2]) - 1

But according to the actual mathematical formula, it should have been this:

x[i] = x[i-1] - f(x[i-1])/f'(x[i-1])

Similarly, x[i-1] = x[i-2] - f(x[i-2])/f'(x[i-2]) - 2

Comparing 1 and 2, we can see that the x[i] in 1 is actually x[i-1] according to the math formula.

The main point to note here is that x and x_p are always one iteration apart. That is, x is the actual successor to x_p, unlike what it might seem by just looking at the code.

Hence, it is working correctly as expected.

What is happening in version 2:

Just like the above case, the same thing happens at the statement x = x_p - (f(x_p)/float(f_p(x_p))).
But by the time we reach if (abs(x-x_p) < prec), x_p had changed its value to temp = x = x[i-1].

But as deduced in the case of version 1, x too is x[i-1] rather than x[i].

So, abs(x - x_p) translates to abs(x[i-1] - x[i-1]), which will turn out to be 0, hence terminating the iteration.

The main point to note here is that x and x_p are actually the same values numerically, which always results in the algorithm terminating after just 1 iteration itself.

Anmol Singh Jaggi
  • 8,376
  • 4
  • 36
  • 77
  • Sorry, I should have clarified. The code I posted, which as I understand it compares to x from two iterations ago, seems to work as I would expect, while with the if statement after, seems to work incorrectly. Why is this? I will update the question to help explain – Aristides May 16 '16 at 23:06
  • What confuses me is that the x[i-2] version seems to be the one that works while the x[i-1] version does not. – Aristides May 16 '16 at 23:08
  • @Aristides Updated the answer! Hope it's clear this time! – Anmol Singh Jaggi May 17 '16 at 00:03
  • Much much better, however something still isn't clicking for me. Why is `x[i]` actually `x[i-1]` in the first version? – Aristides May 17 '16 at 02:13
  • Actually `x[i]` is `x[i-1]` in both the versions. If it was `x = x - (f(x)/float(f_p(x)))`, it would have indeed been `x[i]`. Anyway, don't take the absolute values of the indices too seriously. I just introduced them to establish an ordering between the variables. Their difference is all that matters. – Anmol Singh Jaggi May 17 '16 at 09:50
  • See [this](https://github.com/Anmol-Singh-Jaggi/Machine-Learning/blob/master/newton_raphson/newton.py) for an easier implementation of the algorithm. – Anmol Singh Jaggi May 17 '16 at 09:51
0

Update

This saves the current value of x

tmp = x

In this statement, the next value of x is created from the current value x_p

x = x_p - (f(x_p)/float(f_p(x_p)))

If convergence (eg next value - current value < threshold), break

if (abs(x-x_p) < prec):
    break;

Set x_p to the next iteration's value

x_p = tmp

If you pull x_p = tmp above the if statement, you are actually checking x vs x` from 2 iterations ago, which is not what you want to do. This actually causes weird behavior where the correctness of the outcome is dependent on the starting values. If you start x at you get the correct response whereas if you start with 1, you will not.

To test it out and see why, it can be helpful to add in a print statement as below.

def newton_method(f, f_p, prec=0.01):
    x=7
    x_p=1
    tmp=0

    while(True):
        tmp = x
        x = x_p - (f(x_p)/float(f_p(x_p)))
        print (x,x_p,_tmp)
        if (abs(x-x_p) < prec):
            break;
        x_p = tmp

Are you trying to check X vs X from 2 iterations ago? Or X from the previous iteration of the loop?

If you have x_p=tmp before the if statement, if (abs(x-x_p) < prec): will check the current value of x versus the previous version of x, instead of x from 2 assignments ago

David
  • 11,245
  • 3
  • 41
  • 46
  • okay. so it is indeed checking x from two assignments ago. However why does checking the x from one assignment ago cause the algorithm to fail? That's essentially my question. – Aristides May 16 '16 at 20:09
  • updated question to clarify the nature of my question a little more – Aristides May 16 '16 at 23:10