最优化算法python实现篇（4）——无约束多维极值（梯度下降法）
摘要
本文介绍了多维无约束极值优化算法中的梯度下降法，通过python进行实现，并可视化展示了算法过程。
算法简介
给定初始点，沿着负梯度方向（函数值下降最快的方向）按一定步长（机器学习中也叫学习率）进行搜索，直到满足算法终止条件，则停止搜索。
注意事项
学习率不能太小，也不能太大，可以多尝试一些值。当然每次沿着负梯度方向搜索时，总会存在一个步长使得该次搜索的函数值最低，也就是一个一维无约束极值问题，可调用黄金分割法的一维无约束优化方法求取最佳步长（学习率）。
算法适用性
1、有可能会陷入局部小值。
2、适用于凸函数，由于线性回归的损失函数（Loss Function）是凸函数，所以该算法的应用之一就是解决线性回归问题。
python实现
基本参数：
func:优化的目标函数
x0：初始化变量值
alpha：学习率，一般指定为（0-1），若不指定，则调取一维极值搜索法（黄金分割法）进行求取最优学习率值。黄金分割法代码可参考我的博客：黄金分割法.
黄金分割法内部嵌套了进退法求取一个凸区间。进退法代码参考我的博客：进退法.
epoch：最大迭代次数，若不指定默认为1000
eps:精度，默认为：1e-6
from sympy import *
import numpy as np
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
class CyrusGradientDescent(object):
"""
func:优化的目标函数
x0：初始化变量值
alpha：学习率，一般指定为（0-1），若不指定，则调取一维极值搜索法（黄金分割法）进行求取最优学习率值
黄金分割法代码可参考我的博客：https://blog.csdn.net/Cyrus_May/article/details/105877363
黄金分割法内部嵌套了进退法求取一个凸区间。
进退法代码参考我的博客：https://blog.csdn.net/Cyrus_May/article/details/105821131
epoch：最大迭代次数，若不指定默认为1000
eps:精度，默认为：1e-6
"""
# 1、初始化输入参数
def __init__(self,func,x0,**kargs):
self.var = [Symbol("x"+str(i+1)) for i in range(int(len(x0)))]
func_input = "func(("
for i in range(int(len(x0))):
if i != int(len(x0))-1:
func_input += "self.var[" + str(i) + "]" + ","
else:
func_input += "self.var[" + str(i) + "]" + "))"
self.func = eval(func_input)
self.x = np.array(x0).reshape(-1,1)
if "alpha" in kargs.keys():
self.alpha = kargs["alpha"]
else:
self.alpha = None
if "epoch" in kargs.keys():
self.epoch = kargs["epoch"]
else:
self.epoch = 1e3
if "eps" in kargs.keys():
self.eps = kargs["eps"]
else:
self.eps = 1e-6
self.process = []
self.process.append(self.x)
# 2、定义计算函数值函数
def cal_func_value(self,x):
func = self.func
for i in range(x.shape[0]):
func = func.subs(self.var[i],x[i,0])
return func
# 3、定义计算雅克比矩阵，即梯度的函数
def cal_gradient(self):
f = Matrix([self.func])
v = Matrix(self.var)
gradient = f.jacobian(v)
gradient_value = []
for diff_func in list(gradient):
for i in range(len(self.var)):
diff_func = diff_func.subs(self.var[i],self.x[i,0])
gradient_value.append(diff_func)
return np.array(gradient_value).reshape(-1,1)
# 4、定义 若未指定学习率α时，计算最优学习率的函数
def cal_alpha(self,gradient_value):
if self.alpha != None:
return self.alpha
else:
def alpha_func(alpha):
x = self.x - alpha*gradient_value
return self.cal_func_value(x)
from minimize_golden import Minimize_Golden
return Minimize_Golden(func = alpha_func).run()[0]
# 5、定义更新变量值的函数
def update_x(self,alpha,gradient_value):
self.x = self.x - alpha*gradient_value
self.process.append(self.x)
# 6、定义可视化函数(当目标函数只有两个自变量时才使用)
def visual(self,x1,x2):
X1,X2 = np.meshgrid(x1,x2)
Z = np.ones(X1.shape)
for i in range(X1.shape[0]):
for j in range(X1.shape[1]):
Z[i,j] = self.cal_func_value(np.array([X1[i,j],X2[i,j]]).reshape(-1,1))
fig = plt.figure(figsize=(16,8))
z = []
x = []
y = []
for i in range(len(self.process)):
z.append(self.cal_func_value(self.process[i]))
x.append(self.process[i][0,0])
y.append(self.process[i][1,0])
ax = fig.add_subplot(1,1,1,projection = "3d")
ax.plot_wireframe(X1,X2,Z,rcount = 20,ccount = 20)
ax.plot(x,y,z,color = "r",marker = "*")
# 7、统筹运行
def run(self):
for i in range(int(self.epoch)):
# 1、计算梯度
gradient_value = self.cal_gradient()
if (gradient_value == 0).all():
return self.x,self.cal_func_value(self.x)
# 2、计算学习率α
alpha = self.cal_alpha(gradient_value)
# 3、更新变量值
x_old = self.x
self.update_x(alpha,gradient_value)
if np.abs(self.cal_func_value(x_old)-self.cal_func_value(self.x)) < self.eps:
return self.x,self.cal_func_value(self.x)
return self.x,self.cal_func_value(self.x)
if __name__ == "__main__":
def func(x):
return x[0]**2+x[1]**2+100
gd_model = CyrusGradientDescent(func = func,x0 = (-5,5),alpha = 0.1)
x,y_min = gd_model.run()
print("*"*10,"Gradient Descent Algorithm","*"*10)
print("x:",x)
print("y_min:",y_min)
x1 = np.linspace(-5,5,100)
x2 = np.linspace(-5,5,100)
gd_model.visual(x1,x2)
实例运行结果
********** Gradient Descent Algorithm **********
x: [[-0.000830767497365573]
[0.000830767497365573]]
y_min: 100.000001380349
算法过程可视化
by CyrusMay 2020 05 08
直到文明又毁灭
一千世纪后的第一天
伊甸园里肩并肩
我们笑看太阳也熄灭
——五月天（一千个世纪）——