美文网首页
牛顿法和最速下降法的Python实现

牛顿法和最速下降法的Python实现

作者: 一个扫地的垃圾 | 来源:发表于2019-04-18 22:07 被阅读0次

1 牛顿法

1.1 牛顿法的Python程序

from sympy import *
import numpy as np
# 假设多元函数是二维形式
# x_init为二维向量(x1, x2)
def newton_dou(step, x_init, obj):
    i = 1 # 记录迭代次数的变量
    while i <= step:
        if i == 1:
            grandient_obj = np.array([diff(obj, x1).subs(x1, x_init[0]).subs(x2, x_init[1]), diff(obj, x2).subs(x1, x_init[0]).subs(x2, x_init[1])], dtype=float) # 初始点的梯度值
            hessian_obj = np.array([[diff(obj, x1, 2), diff(diff(obj, x1), x2)], [diff(diff(obj, x2), x1), diff(obj, x2, 2)]], dtype=float) # 初始点的hessian矩阵
            inverse = np.linalg.inv(hessian_obj) # hessian矩阵求逆
            x_new = x_init - np.matmul(inverse, grandient_obj) # 第一次迭代公式
            print(x_new)
            # print('迭代第%d次:%.5f' %(i, x_new))
            i = i + 1
        else:
            grandient_obj = np.array([diff(obj, x1).subs(x1, x_new[0]).subs(x2, x_new[1]), diff(obj, x2).subs(x1, x_new[0]).subs(x2, x_new[1])], dtype=float) # 当前点的梯度值
            hessian_obj = np.array([[diff(obj, x1, 2), diff(diff(obj, x1), x2)], [diff(diff(obj, x2), x1), diff(obj, x2, 2)]], dtype=float) # 当前点的hessian矩阵
            inverse = np.linalg.inv(hessian_obj) # hessian矩阵求逆
            x_new = x_new - np.matmul(inverse, grandient_obj) # 迭代公式
            print(x_new)
            # print('迭代第%d次:%.5f' % (i, x_new))
            i = i + 1
    return x_new

x0 = np.array([0, 0], dtype=float)
x1 = symbols("x1")
x2 = symbols("x2")
newton_dou(5, x0, x1**2+2*x2**2-2*x1*x2-2*x2)

1.2 牛顿法的结果分析

    程序执行的结果如下:

[1. 1.]
[1. 1.]
[1. 1.]
[1. 1.]
[1. 1.]
Process finished with exit code 0

    经过实际计算函数f(x)=x_1^2+2x_2^2-2x_1x_2-2x_2的极值点为(1,1),只需一次迭代就能收敛到极值点。

2 最速下降法

2.1 最速下降法的Python程序

# !/usr/bin/python
# -*- coding:utf-8 -*-
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
def sinvar(fun):
    s_p = solve(diff(fun)) #stationary point
    return s_p
def value_enter(fun, x, i):
    value = fun[i].subs(x1, x[0]).subs(x2, x[1])
    return value
def grandient_l2(grand, x_now):
    grand_l2 = sqrt(pow(value_enter(grand, x_now, 0), 2)+pow(value_enter(grand, x_now, 1), 2))
    return grand_l2
def msd(x_init, obj):
    i = 1
    grandient_obj = np.array([diff(obj, x1), diff(obj, x2)])
    error = grandient_l2(grandient_obj, x_init)
    plt.plot(x_init[0], x_init[1], 'ro')
    while(error>0.001):
        if i == 1:
            grandient_obj_x = np.array([value_enter(grandient_obj, x_init, 0), value_enter(grandient_obj, x_init, 1)])
            t = symbols('t')
            x_eta = x_init - t * grandient_obj_x
            eta = sinvar(obj.subs(x1, x_eta[0]).subs(x2, x_eta[1]))
            x_new = x_init - eta*grandient_obj_x
            print(x_new)
            i = i + 1
            error = grandient_l2(grandient_obj, x_new)
            plt.plot(x_new[0], x_new[1], 'ro')
        else:
            grandient_obj_x = np.array([value_enter(grandient_obj, x_new, 0), value_enter(grandient_obj, x_new, 1)])
            t = symbols('t')
            x_eta = x_new - t * grandient_obj_x
            eta = sinvar(obj.subs(x1, x_eta[0]).subs(x2, x_eta[1]))
            x_new = x_new - eta*grandient_obj_x
            print(x_new)
            i = i + 1
            error = grandient_l2(grandient_obj, x_new)
            plt.plot(x_new[0], x_new[1], 'ro')
    plt.show()
x_0 = np.array([0, 0], dtype=float)
x1 = symbols("x1")
x2 = symbols("x2")
result = msd(x_0, x1**2+2*x2**2-2*x1*x2-2*x2)
print(result)

2.2 最速下降法结果分析

    程序执行的结果如下:

[0 0.500000000000000]
[0.500000000000000 0.500000000000000]
[0.500000000000000 0.750000000000000]
[0.750000000000000 0.750000000000000]
[0.750000000000000 0.875000000000000]
[0.875000000000000 0.875000000000000]
[0.875000000000000 0.937500000000000]
[0.937500000000000 0.937500000000000]
[0.937500000000000 0.968750000000000]
[0.968750000000000 0.968750000000000]
[0.968750000000000 0.984375000000000]
[0.984375000000000 0.984375000000000]
[0.984375000000000 0.992187500000000]
[0.992187500000000 0.992187500000000]
[0.992187500000000 0.996093750000000]
[0.996093750000000 0.996093750000000]
[0.996093750000000 0.998046875000000]
[0.998046875000000 0.998046875000000]
[0.998046875000000 0.999023437500000]
[0.999023437500000 0.999023437500000]
[0.999023437500000 0.999511718750000]
None

Process finished with exit code 0

    最速下降法选择的停止条件为迭代的梯度值如果小于0.001,则停止迭代,通过结果可看出设定停止条件的情况下,极值点并没有完全收敛到(1,1)点,在Python中画出了迭代的极值点(如图1所示)。如果继续缩小停止条件,则程序会延长时间增加迭代次数,最后极值点达到(1,1)

图1 最速下降法迭代的极值点

相关文章

  • 牛顿法和最速下降法的Python实现

    1 牛顿法 1.1 牛顿法的Python程序 1.2 牛顿法的结果分析     程序执行的结果如下:     经过...

  • 第十章 共轭方向法

    10.1引言 从效率上,共轭方向法位于最速下降法和牛顿法之间,具有以下特性:1、对于n维二次型问题,能够在n步之内...

  • 无约束最优化(二) 共轭方向法与共轭梯度法

    基本思想   之前文章最速下降法、Newton法、修正Newton法介绍的最速下降法存在锯齿现象,Newton法需...

  • 牛顿法python3实现

    本文参考链接:牛顿法实现偏导数求解 输出结果 从输出结果可以看出最终收敛至(0,0),即f在(0,0)处取得极小值。

  • 第九章 牛顿法

    9.1 引言 在确定搜索方向时,最速下降法只利用了一阶导数(梯度)。这并不是最高效的。因此引入牛顿法,它同时使用一...

  • 最速下降法的python实现

    代码一: 分析1:这个采用的是backtracking line search来寻找alpha。

  • Newton's method and Quasi Ne

    Welcome To My Blog 牛顿法和拟牛顿法是求解无约束最优化问题的常用方法,优点是收敛速度快.牛顿法...

  • Python梯度下降法

    本文主要讲解梯度下降算法,以及Python的实现一个简单的例子 梯度下降法又称为最速下降法,是 1847 年有数学...

  • 【转】常见的几种最优化方法

    转自Poll 的笔记 阅读目录 梯度下降法(Gradient Descent) 牛顿法和拟牛顿法(Newton's...

  • 牛顿线性逼近求开方

    牛顿法线性逼近计算开方 python Code: 第一步:开方问题转化为求根问题: 第二步:牛顿逼近法原理如下:

网友评论

      本文标题:牛顿法和最速下降法的Python实现

      本文链接:https://www.haomeiwen.com/subject/ojvxgqtx.html