Python实现最速下降法(The steepest descent method)详细案例
文章目录
- 最速下降法(The steepest descent method)
- 最速下降法的原理
- Python实现最速下降法实例
- `sympy`包中用到的函数
- 构建符号变量和符号函数
- 对符号函数求导
- 求函数值
- 求解方程的零点
- Python实现最速下降法求解上述算例的完整代码
.[TOC].
最速下降法(The steepest descent method)
本文中的课件来自清华大学深圳国际研究生院，物流与交通学部张灿荣教授《高级运筹学》课程。
在无约束非线性函数的最优化中，最速下降法(The steepest descent method)是一个著名的基础算法。本文就来用一个实例来学习最速下降法。
最速下降法的原理
假设有一个多元非线性函数 f ( x 1 , x 2 , ⋯ , x n ) f(x_1, x_2, \cdots, x_n) f(x1,x2,⋯,xn)，其定义域为 x ∈ R n x \in \mathbb{R}^n x∈Rn, 那么如何求该函数的最大值或者最小值呢？
当然，我们可以用求导的方法。但是有些时候直接求导并不方便，此时我们可以用最速下降法来求解。
最速下降法的基本思想就是：
从一个初始点开始，逐步沿着以当前点为基准，函数值变化最快的方向走，一直走到最优解为止。因此，在有了一个初始点以后，我们就需要决策以下两个事情：
（1）下一步要朝着什么方向走(方向);
（2）沿着该方向走多远(步长)。
具体如何选择方向和步长，见下面的PPT。
下图是一个很直观的例子，在当前点，沿着不同的方向走，函数值的变化速度是不同的，但是在等高线的切线的垂直方向，是变化最快的。
也就是说，我们得去找平行于该点梯度的方向，沿着这个方向(当为max问题)或者沿着这个方向的反方向(当为min问题)去更新当前位置。
就像上图一样，一步一步，最终走到最优解对应的点。
Python实现最速下降法实例
我们看下面的函数
f ( x 1 , x 2 ) = 2 x 1 x 2 + 2 x 2 − x 1 2 − 2 x 2 2 , x 1 , x 2 ∈ R 1 \begin{aligned} f(x_1,x_2)=&2x_1x_2 +2x_2-x_1^2-2x_2^2 , \\ &x_1, x_2 \in \mathbb{R^1} \end{aligned} f(x1,x2)=2x1x2+2x2−x12−2x22,x1,x2∈R1
给定初始点 ( 0.5 , 0.5 ) (0.5, 0.5) (0.5,0.5)，用最速下降法找到该函数的最大值。
我们用python中的sympy包来实现最速下降法求解上面的问题。
首先，我们来介绍我们将要用到的sympy包中的几个函数。结合jupyter notebook来给大家做个展示。ps.这个包还是挺好用的，结合jupyter notebook之后，可视化非常不错，所见即所得。
sympy包中用到的函数
构建符号变量和符号函数
from sympy import *
x_1 = symbols('x_1')
x_2 = symbols('x_2') fun = 2 * x_1 * x_2 + 2 * x_2 - x_1**2 - 2 * x_2**2
fun
这个是用来构造两个符号变量 x 1 , x 2 x_1, x_2 x1,x2，就像代数中用字母代替变量一样。然后可以定义出我们的函数
f ( x 1 , x 2 ) = 2 x 1 x 2 + 2 x 2 − x 1 2 − 2 x 2 2 , \begin{aligned} f(x_1,x_2)=&2x_1x_2 +2x_2-x_1^2-2x_2^2 , \end{aligned} f(x1,x2)=2x1x2+2x2−x12−2x22,
jupyter notebook中的显示效果是这样的
可以看到jupyter notebook中直接就显示出了数学公式格式的形式，这是因为jupyter notebook中内嵌了LaTeX相关支持包的缘故。总之这样可视化就非常不错。
对符号函数求导
下面我们来计算 f ( x 1 , x 2 ) f(x_1,x_2) f(x1,x2)的偏导数， ∂ f ∂ x 1 \frac{\partial f}{\partial x_1} ∂x1∂f，代码很简单，用函数diff(函数, 变量)
grad_1 = diff(fun, x_1)
grad_1
如下图
求函数值
有了符号函数，我们怎么知道自变量 x 1 , x 2 x_1, x_2 x1,x2取具体值得时候，符号函数的取值呢？我们用函数函数.subs({变量1:变量1的取值, ....}).evalf()，具体代码为
fun_value = fun.subs({x_1:4, x_2: 2}).evalf()
fun_value
如果说只有一个变量已知，那就是下面的情况
fun_value = fun.subs({x_1:4, x_2: 2}).evalf()
fun_value
还是非常智能的。
求解方程的零点
为了寻找下降速度最快的方向，我们需要利用之前PPT中的方法去求解方程组。这里我们用到函数solve(函数，变量)
我们举一个例子，假设有函数
g ( x ) = 4 x + 5 \begin{aligned} g(x) = 4x + 5 \end{aligned} g(x)=4x+5
我们求解
g ( x ) = 4 x + 5 = 0 \begin{aligned} g(x) = 4x + 5 = 0 \end{aligned} g(x)=4x+5=0
OK，准备好了上述函数，我们就可以开心的把最速下降法写出来了。具体代码见接下来的部分。
Python实现最速下降法求解上述算例的完整代码
import math
from sympy import *# define symbol variable
x_1 = symbols('x_1')
x_2 = symbols('x_2')# define objective function
fun = 2 * x_1 * x_2 + 2 * x_2 - x_1 ** 2 - 2 * x_2 ** 2
fun# take derivative of x_1 and x_2
grad_1 = diff(fun, x_1)
grad_2 = diff(fun, x_2)# define parameters
MaxIter = 100
epsilon = 0.0001# define initial point
x_1_value = 0.5
x_2_value = 0.5iter_cnt = 0
current_step_size = 10000 grad_1_value = (float)(grad_1.subs({x_1:x_1_value, x_2: x_2_value}).evalf())
grad_2_value = (float)(grad_2.subs({x_1:x_1_value, x_2: x_2_value}).evalf()) current_obj = fun.subs({x_1: x_1_value, x_2: x_2_value}).evalf()print('itCnt: %2d cur_point (%3.2f, %3.2f) cur_Obj: %5.4f grad_1: %5.4f grad_2 : %5.4f step_size : %5.4f' % (iter_cnt, x_1_value, x_2_value, current_obj, grad_1_value, grad_2_value, current_step_size)) # while (iter_cnt <= MaxIter and abs(grad_1_value) + abs(grad_2_value) >= epsilon):
while(abs(grad_1_value) + abs(grad_2_value) >= epsilon): iter_cnt += 1# find the step sizet = symbols('t')x_1_updated = x_1_value + grad_1_value * tx_2_updated = x_2_value + grad_2_value * tFun_updated = fun.subs({x_1: x_1_updated, x_2: x_2_updated})grad_t = diff(Fun_updated, t)t_value = solve(grad_t, t)[0] # solve grad_t == 0# update x_1_value and x_2_valuegrad_1_value = (float)(grad_1.subs({x_1: x_1_value, x_2: x_2_value}).evalf()) grad_2_value = (float)(grad_2.subs({x_1: x_1_value, x_2: x_2_value}).evalf()) x_1_value = (float)(x_1_value + t_value * grad_1_value)x_2_value = (float)(x_2_value + t_value * grad_2_value) current_obj = fun.subs({x_1: x_1_value, x_2: x_2_value}).evalf()current_step_size = t_valueprint('itCnt: %2d cur_point (%3.2f, %3.2f) cur_Obj: %5.4f grad_1: %5.4f grad_2 : %5.4f step_size : %5.4f' % (iter_cnt, x_1_value, x_2_value, current_obj, grad_1_value, grad_2_value, current_step_size))
运行结果如下
itCnt: 0 cur_point (0.50, 0.50) cur_Obj: 0.7500 grad_1: 0.0000 grad_2 : 1.0000 step_size : 10000.0000
itCnt: 1 cur_point (0.50, 0.75) cur_Obj: 0.8750 grad_1: 0.0000 grad_2 : 1.0000 step_size : 0.2500
itCnt: 2 cur_point (0.50, 0.75) cur_Obj: 0.8750 grad_1: 0.5000 grad_2 : 0.0000 step_size : 0.0000
itCnt: 3 cur_point (0.75, 0.75) cur_Obj: 0.9375 grad_1: 0.5000 grad_2 : 0.0000 step_size : 0.5000
itCnt: 4 cur_point (0.75, 0.75) cur_Obj: 0.9375 grad_1: 0.0000 grad_2 : 0.5000 step_size : 0.0000
itCnt: 5 cur_point (0.75, 0.88) cur_Obj: 0.9688 grad_1: 0.0000 grad_2 : 0.5000 step_size : 0.2500
itCnt: 6 cur_point (0.75, 0.88) cur_Obj: 0.9688 grad_1: 0.2500 grad_2 : 0.0000 step_size : 0.0000
itCnt: 7 cur_point (0.88, 0.88) cur_Obj: 0.9844 grad_1: 0.2500 grad_2 : 0.0000 step_size : 0.5000
itCnt: 8 cur_point (0.88, 0.88) cur_Obj: 0.9844 grad_1: 0.0000 grad_2 : 0.2500 step_size : 0.0000
itCnt: 9 cur_point (0.88, 0.94) cur_Obj: 0.9922 grad_1: 0.0000 grad_2 : 0.2500 step_size : 0.2500
itCnt: 10 cur_point (0.88, 0.94) cur_Obj: 0.9922 grad_1: 0.1250 grad_2 : 0.0000 step_size : 0.0000
itCnt: 11 cur_point (0.94, 0.94) cur_Obj: 0.9961 grad_1: 0.1250 grad_2 : 0.0000 step_size : 0.5000
itCnt: 12 cur_point (0.94, 0.94) cur_Obj: 0.9961 grad_1: 0.0000 grad_2 : 0.1250 step_size : 0.0000
itCnt: 13 cur_point (0.94, 0.97) cur_Obj: 0.9980 grad_1: 0.0000 grad_2 : 0.1250 step_size : 0.2500
itCnt: 14 cur_point (0.94, 0.97) cur_Obj: 0.9980 grad_1: 0.0625 grad_2 : 0.0000 step_size : 0.0000
itCnt: 15 cur_point (0.97, 0.97) cur_Obj: 0.9990 grad_1: 0.0625 grad_2 : 0.0000 step_size : 0.5000
itCnt: 16 cur_point (0.97, 0.97) cur_Obj: 0.9990 grad_1: 0.0000 grad_2 : 0.0625 step_size : 0.0000
itCnt: 17 cur_point (0.97, 0.98) cur_Obj: 0.9995 grad_1: 0.0000 grad_2 : 0.0625 step_size : 0.2500
itCnt: 18 cur_point (0.97, 0.98) cur_Obj: 0.9995 grad_1: 0.0312 grad_2 : 0.0000 step_size : 0.0000
itCnt: 19 cur_point (0.98, 0.98) cur_Obj: 0.9998 grad_1: 0.0312 grad_2 : 0.0000 step_size : 0.5000
itCnt: 20 cur_point (0.98, 0.98) cur_Obj: 0.9998 grad_1: 0.0000 grad_2 : 0.0312 step_size : 0.0000
itCnt: 21 cur_point (0.98, 0.99) cur_Obj: 0.9999 grad_1: 0.0000 grad_2 : 0.0312 step_size : 0.2500
itCnt: 22 cur_point (0.98, 0.99) cur_Obj: 0.9999 grad_1: 0.0156 grad_2 : 0.0000 step_size : 0.0000
itCnt: 23 cur_point (0.99, 0.99) cur_Obj: 0.9999 grad_1: 0.0156 grad_2 : 0.0000 step_size : 0.5000
itCnt: 24 cur_point (0.99, 0.99) cur_Obj: 0.9999 grad_1: 0.0000 grad_2 : 0.0156 step_size : 0.0000
itCnt: 25 cur_point (0.99, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0156 step_size : 0.2500
itCnt: 26 cur_point (0.99, 1.00) cur_Obj: 1.0000 grad_1: 0.0078 grad_2 : 0.0000 step_size : 0.0000
itCnt: 27 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0078 grad_2 : 0.0000 step_size : 0.5000
itCnt: 28 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0078 step_size : 0.0000
itCnt: 29 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0078 step_size : 0.2500
itCnt: 30 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0039 grad_2 : 0.0000 step_size : 0.0000
itCnt: 31 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0039 grad_2 : 0.0000 step_size : 0.5000
itCnt: 32 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0039 step_size : 0.0000
itCnt: 33 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0039 step_size : 0.2500
itCnt: 34 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0020 grad_2 : 0.0000 step_size : 0.0000
itCnt: 35 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0020 grad_2 : 0.0000 step_size : 0.5000
itCnt: 36 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0020 step_size : 0.0000
itCnt: 37 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0020 step_size : 0.2500
itCnt: 38 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0010 grad_2 : 0.0000 step_size : 0.0000
itCnt: 39 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0010 grad_2 : 0.0000 step_size : 0.5000
itCnt: 40 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0010 step_size : 0.0000
itCnt: 41 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0010 step_size : 0.2500
itCnt: 42 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0005 grad_2 : 0.0000 step_size : 0.0000
itCnt: 43 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0005 grad_2 : 0.0000 step_size : 0.5000
itCnt: 44 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0005 step_size : 0.0000
itCnt: 45 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0005 step_size : 0.2500
itCnt: 46 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0002 grad_2 : 0.0000 step_size : 0.0000
itCnt: 47 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0002 grad_2 : 0.0000 step_size : 0.5000
itCnt: 48 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0002 step_size : 0.0000
itCnt: 49 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0002 step_size : 0.2500
itCnt: 50 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0001 grad_2 : 0.0000 step_size : 0.0000
itCnt: 51 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0001 grad_2 : 0.0000 step_size : 0.5000
itCnt: 52 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0001 step_size : 0.0000
itCnt: 53 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0000 grad_2 : 0.0001 step_size : 0.2500
itCnt: 54 cur_point (1.00, 1.00) cur_Obj: 1.0000 grad_1: 0.0001 grad_2 : 0.0000 step_size : 0.0000
这样就完成了。
当然了，在这个例子中，从第3步左右的迭代开始，后续的点就非常近了，因此，步长就需要动态的调整。具体文献之后补充。
作者：刘兴禄，清华大学，清华伯克利深圳学院 (博士在读)
邮箱：hsinglul@163.com
Python实现最速下降法(The steepest descent method)详细案例相关推荐
- 最优化学习 最速下降法（steepest Descent)
最速下降法(steepest Descent) 最速(陡)下降法(steepest Descent) v正则化为2范数 v为1范数和v为无穷范数 Steepest Gradient的变种 坐标轴交替下 ...
- 最速下降法steepest descent详细解释
----------------------先闲聊几句------------------------------------ [1]首次提出了梯度下降法和最速下降法,既然柯西写出来了,所以这两个算法 ...
- python求最值_用Python实现最速下降法求极值的方法
对于一个多元函数 ,用最速下降法(又称梯度下降法)求其极小值的迭代格式为 其中 为负梯度方向,即最速下降方向,αkαk为搜索步长. 一般情况下,最优步长αkαk的确定要用到线性搜索技术,比如精确线性搜 ...
- Fletcher-Reevers Conjugate Descent和Steepest Descent两种算法中伪代码的区别
本文主要用来比较两个算法到底差别在哪里 step Fletcher-Reevers Conjugate Descent Steepest Descent 1st1st1st 选择初始点x(1)选择初始 ...
- Supervised Descent Method(人脸对齐之SDM论文解析）
Supervised Descent Method(人脸对齐之SDM论文解析) 标签: SDM NLS Jacobian Hessian FaceAlignment 作者:贾金让 本人博客链接:htt ...
- python迭代法求极值_用Python实现最速下降法求极值的方法
对于一个多元函数 ,用最速下降法(又称梯度下降法)求其极小值的迭代格式为 其中 为负梯度方向,即最速下降方向,αkαk为搜索步长. 一般情况下,最优步长αkαk的确定要用到线性搜索技术,比如精确线性搜 ...
- LASSO坐标下降法Coordinate Descent Method公式推导及代码
文章目录 LASSO by Coordinate Descent Method Coordinate Descent Method Framework Coordinate Descent Metho ...
- 【 ML 】Steepest Descent Iteration Procedure of TOA - Based Positioning Simulation
steepest descent algorithms are, 用函数表示为: function g = grad_ml(X,x,r,sigma2) % ML gradient computatio ...
- 【 NLS 】Steepest Descent Algorithm Iteration Procedure of TOA - Based Positioning
这种类型的仿真,前面已经有两篇了: Newton – Raphson Iteration Procedure of TOA - Based Positioning Gauss-Netwon algor ...
最新文章
- pointnet 结果可视化_PointNet论文复现及代码详解
- 为什么“ cd”在shell脚本中不起作用？
- 【Codeforces Round #767 (Div. 2)】 C. Meximum Array 题解
- Servlet—07—Cookie; Seesion;
- php获取上海时间代码,PHP获取星期的方法及代码
- react中实现异步请求的方法一，react-thunk
- Android 系统（254）---Android libphonenumber Demo 手机号码归属地
- 重庆计算机教师招聘 专业技能测试什么,教师招聘考试面试，专业技能测试考什么？全在这了...
- 代写php代码作业,代写phpmyadmin留学生作业、代做SQL语言作业、SQL程序设计作业调试、代做PHP script作业...
- idea+spring boot+jrebel7.0.14热启动
- j2me模拟器linux,J2ME HELLOWORLD 小试牛刀（转）
- excel的图表里如何添加上下标
- 七夕情人节送什么礼物给女朋友？音质好的蓝牙耳机推荐
- SQLite主键自动增长
- 谷歌浏览器打开金格在线编辑插件
- 关于人际关系的一些问答
- MMFNet: A Multi-modality MRI Fusion Network for Segmentation of Nasopharyngeal Carcinoma
- python发送邮件功能
- 【室内温度+树莓派性能监控】树莓派+DS18B20温度传感器+0.96寸OLED显示屏使用及安装经验分享
- get请求400神坑