Newton-Raphson算法

簡介

牛頓法又叫做牛頓-拉裴森（Newton-Raphson）方法，是一維求根方法中最著名的一種。其特點是在計算時需要同時計算函數值與其一階導數值，從幾何上解釋，牛頓法是將當前點處的切線延長，使之與橫軸相交，然后把交點處值作為下一估值點。

 

圖1

從數學上解釋，牛頓法可以從函數的泰勒展開得到。\(f(x)\)的泰勒展開可以表示為：

\(f(x+\delta)=f(x)+f’(x)\delta+\frac{f’’(x)}{2}\delta^2+O(\delta^3)\)

 

對於足夠小的\(\delta\)，可以將只保留上式右端關於的一階項，得到：

 

\(\delta=-\frac{f(x)}{f’(x)}\)

於是得到由到的遞推公式：

 

\(x_{i+1}=x_{i}+\delta=x_i-\frac{f(x_i)}{f’(x_i)}\)

 

可見牛頓法是讓\(x\)沿着\(f(x)\)梯度的方向下降，類似於最優化方法中的梯度下降法。牛頓法也可以作為最優化算法，只不過那時需要求函數的二階導數。

牛頓法相比二分法、截弦法的優點是收斂速度可以達到二階，在根附近沒迭代一次，結果的有效數字幾乎可以翻倍。當然牛頓法也可能可能失敗，比如收斂到一個局部極值，其切線方向與橫軸水平，從而無法計算下一個迭代值。

另外，牛頓法的實現需要用戶提供一個函數用於計算函數值\(f(x)\)與其一階導數值\(f'(x)\)，因此比較適合函數的導數可以解析求出的情況，如果需要求數值導數，則牛頓法的收斂速度和精度都會受影響。

我們可以將牛頓法和二分法綜合起來形成一個混合算法，一旦牛頓法在運行過程中出現解跳出給定區間或者猜測值遠離實際根導致收斂速度較慢時，就采取一步二分法。

實現一：利用預先求出的一階導函數

import numpy as np
import matplotlib.pyplot as plt


def f(FV, PMT, r, n):

    return PMT * (1 + r) * (((1 + r)**n - 1)) / r + FV


def df(FV, PMT, r, n):
    r_plus_1_power_n = (1 + r)**n

    p1 = N * PMT * r_plus_1_power_n / r
    p2 = -PMT * (r + 1) * (r_plus_1_power_n - 1) / r / r
    p3 = PMT * (r_plus_1_power_n - 1) / r
    return p1 + p2 + p3


def newtonRaphson2(FV,
                   PMT,
                   n,
                   f,
                   df,
                   xmin,
                   xmax,
                   maxit,
                   shift=0.0001,
                   tol=1.0e-9):
    '''
    函數作用說明：計算組合收益率
    FV：目標金額
    PMT：每期投資金額
    n:定投期數
    f：函數值（根據要求的方程自定義）
    df：導數值（根據要求的方程自定義）
    xmin：根的下限
    xmax：根的上限
    maxit：最大迭代次數
    tol：計算精度
    '''
    import math
    fxmin = f(FV, PMT, xmin, n)
    if fxmin == 0.0:
        return xmin
    fxmax = f(FV, PMT, xmax, n)
    if fxmax == 0.0:
        return xmax
    if fxmin * fxmax > 0.0:
        print(
            'Root is not bracketed')  # 在[xmin, xmax]內函數不變號（沒根），或者是變了偶數次號（多個根）
        return 1
    if fxmin < 0.0:  # 確定搜索方向使f(xl)<0
        xl = xmin
        xh = xmax
    else:
        xl = xmax
        xh = xmin
    x = 0.5 * (xmin + xmax)  # 根的預測值
    if x == 0:
        x += shift
    fx, dfx = f(FV, PMT, x, n), df(FV, PMT, x, n)  # 求f(x)和其一階導數
    dxold = math.fabs(xmax - xmin)  # 儲存步長
    dx = dxold
    for ii in range(maxit):
        # 牛頓法的解跳出解區間或者收斂速度太慢，則下一步改用二分法
        if ((x - xh) * dfx - fx) * ((x - xl) * dfx - fx) > 0.0 or (
                math.fabs(2 * fx) > math.fabs(dxold * dfx)):
            # 二分法
            dxold = dx
            dx = 0.5 * (xh - xl)
            x = xl + dx
        else:
            # 牛頓法
            dxold = dx
            dx = fx / dfx
            temp = x
            x -= dx
            if temp == x:
                print("total iterate time:%s " % ii)
                return x
        if math.fabs(dx) < tol:  # 達到要求精度，返回找到的根
            print("total iterate time:%s " % ii)
            return x
        if x == 0:
            x += shift
        fx, dfx = f(FV, PMT, x, n), df(FV, PMT, x, n)  # 否則繼續迭代，求f(x)和其一階導數
        if fx < 0.0:  # 使根保持在解區間內
            xl = x
        else:
            xh = x
    print('Maximum number of iterations exceeded')
    return 1


### 測試用例：首先給定PMT，n，r_analytical，計算FV，然后利用PMT，n，FV計算r_numerical，兩者應該相等
##給定r_analytical計算FV
R=0.1
r_analytical = R / 12
PMT = -4e3
N = 30
n = N * 12
FV = -PMT * (1 + r_analytical) * (((1 + r_analytical)**n - 1)) / r_analytical

##給定FV反解r_numerical
r_numerical = newtonRaphson2(FV, PMT, n, f, df, -1, 1, 100, tol=1.0e-8)

print('\nr_analytical=%s,\nr_numerical=%s\n' % (r_analytical, r_numerical))

 

實現二：利用TensorFlow提供的自動微分計算導函數

import numpy as np
import math
import pandas as pd
import tensorflow as tf
import matplotlib.pyplot as plt

##一個利用tensorflow的自動微分功能實現牛頓法解方程的小程序


class NewtonRaphson:
    def __init__(self, y, x, session):
        self.y = y
        self.x = x
        self.grad = tf.gradients(y, x)
        self.sess = session
        sess.run(tf.global_variables_initializer())

    def _fx(self, x_value):
        # 盡量避免出現f(x)不能計算的情況，比如函數試圖計算a/0，log(0)等，如果計算結果為inf則x+0.0001再進行計算

        temp = self.sess.run(y, feed_dict={x: [x_value]})[0]
        if np.isinf(temp):
            return self.sess.run(y, feed_dict={x: [x_value + 0.0001]})[0]
        else:
            return temp

    def _dfx(self, x_value):
        return self.sess.run(self.grad, feed_dict={x: [x_value]})[0][0]

    def solve(self, xmin, xmax, maxiter, tol):
        fmin = self._fx(xmin)
        fmax = self._fx(xmax)
        if fmin == 0:
            return xmin
        if fmax == 0:
            return xmax
        if fmin * fmax > 0.0:
            raise ValueError('Root is not brackted!!')
        if fmin < 0:
            xl = xmin
            xh = xmax
        else:
            xl = xmax
            xh = xmin

        x = (xmin + xmax) / 2
        fx, dfx = self._fx(x), self._dfx(x)

        dxold = math.fabs(xmax - xmin)
        dx = dxold

        for ii in range(maxiter):
            if ((x - xh) * dfx - fx) * ((x - xl) * dfx - fx) > 0.0 or (
                    math.fabs(2 * fx) > math.fabs(dxold * dfx)):
                dxold = dx
                dx = 0.5 * (xh - xl)
                x = xl + dx
            else:
                dxold = dx
                dx = fx / dfx
                temp = x
                x -= dx
                # newton
                if temp == x:
                    print("total iterate time:%s " % ii)
                    return x
            fx, dfx = self._fx(x), self._dfx(x)
            if fx < 0.0:
                xl = x
            else:
                xh = x
        print('Maximum number of iterations exceeded')
        return 1



PV = 1e4
FV = 3e6
N = 20
cpi = 0.018
RATE = 0.15
r = RATE / 12
PMT = 10000


x = tf.placeholder(shape=[1], dtype=tf.float32)
y=r * (FV * (1 + cpi)**(N) - PV * (r + 1)**x) / ((r + 1)**x - 1 - r) - PMT

sess = tf.InteractiveSession()
solver=NewtonRaphson(y,x,sess)


nmin = 2
nmax = 300

solver.solve(nmin,nmax,100,1e-9)