作业笔记:基于二次插值的Wolfe-Powell非精确线搜索算法及Python代码实现
写在前面的说明:
- 查阅了很多资料,发现资料里对于非精确线搜索的Wolfe-Powell准则求步长都讲得很粗糙,中英文的教材里(黄红选.数学规划 4.2.1节、Jorge Nocedal.Numerical Optimization 3.1节等)只介绍了Wolfe-Powell准则(Wolfe Conditions)的原理,没有给出算法步骤;而能查到的W-P代码实现,都只给了代码,也没有给算法步骤。
- 关于Wolfe Conditions的数学原理可以看Numerical Optimization一书,不再具体介绍。这篇笔记主要介绍了怎样把原理转化为算法步骤,再转化为Python代码。
- 关于二次插值法,大多数教材里面只讲了三点二次插值法(陈宝林.最优化理论与算法 9.3.3节),这篇笔记里二点二次插值的内插和外插公式是自己推的,如果有问题也可以指出~
1 准备知识:二次插值法
1.1 概述
二次插值法(抛物线法)基本思路:在极小点附近,用二次三项式逼近目标函数
分为三点二次插值法和二点二次插值法:
- 三点二次插值法:已知三点的函数值,求极小点
- 二点二次插值法:
- 已知两点的函数值,和其中一点的导数值,求极小点(内插)
- 已知两点的导数值,求极小点(外插)
1.2 二点二次插值法
(1)内插:找两点之间的极小点
令二次三项式
已知,,,其中
为保证二次插值函数有极小点,要求或者
得到以下方程组:
解方程组,得,,。
为求的极小点,则令 ,,将代入得极小点
(2)外插:找两点之外的极小点
令二次三项式
已知,
为保证二次插值函数有极小点,要求
得到以下方程组:
解方程组,得,。
为求的极小点,则令 ,,将代入得极小点
2 准备知识:Wolfe Conditions
其中
3 基于二点二次插值的Wolfe-Powell法
3.1 基本思路
- 满足,且满足,直接输出步长
- 满足,但不满足,用外插公式增大步长(见Step 3)
- 不满足,用内插公式缩小步长(见Step 2)
3.2 计算步骤
Step 1.
选取初始数据,给定初始搜索区间,给出参数,令,取,计算,。
def WolfePowell(f,d,x,alpha_,rho,sigma):
a1 = 0
a2 = alpha_
alpha = (a1+a2)/2
输入的初始数据有,f
:目标函数,d
:初始方向,x
:初始点,alpha_
:,rho
:,sigma
:。另外,a1
:,a2
:,alpha
:,算法中取。
phi1 = fun(x)
phi1_ = np.dot(grad,d)
phi1
:,phi1_
:,求梯度的方法会在后面说明。
Step 2.
计算。
若(条件1)成立,转到Step 3;否则,由二点二次插值法的内插公式计算:
令,,转到Step 2
if phi <= phi1 + rho * alpha * phi1_:
……
else:
alpha_new = a1 - 0.5*(alpha-a1)**2*phi1_/((phi-phi1)-(alpha-a1)*phi1_)
a2 = alpha
alpha = alpha_new
Step 3.
计算。
若(条件2)成立,则令,输出,停止;否则,由二点二次插值法的外插公式计算:
if phi <= phi1 + rho * alpha * phi1_:
……
if phi_ >= sigma*phi0_:
break
else:
alpha_new = alpha + (alpha-a1) *phi_ / (phi1_-phi_)
a1 = alpha
alpha = alpha_new
phi1 = phi
phi1_ = phi_
else:
……
4 对Wolfe-Powell法的改进
4.1 计算步长方法的改进
- 网上大多数的W-P程序中,采用如下的方法计算步长α
本文的算法中,二点二次插值法的内插和外插,分别对应上面第(3)步和第(4)步。但不同之处在于对步长放缩的方式:
-
转到第(3)步时,不满足,所以要缩小步长。此时步长太大,取和步长下界的中值为新步长,就能使步长缩小。因为,所以步长最多缩小到原来的。但是,这种方法得到的新步长,并不一定能使目标函数最小。采用二次插值法的内插法,可以找到之间的极小点,能使目标函数近似最小,提高了算法效率。
-
转到第(4)步时,满足,但不满足,所以要增大步长。此时步长太小,取和步长上界的中值为新步长,就能使步长放大。因为,步长最多放大到原来的倍。但是,这种方法得到的新步长,并不一定能使目标函数最小。采用二次插值法的外插法,可以找到之外(大于一侧)的极小点,能使目标函数近似最小,提高了算法效率。
4.2 求梯度方法的改进
-
在有些程序中,需要手动计算梯度并以数组形式输入
def fun(x): return x[0]**2+2*x[1]**2-2*x[0]*x[1]+2*x[1]+2 def gfun(x): return np.array([2*x[0]-2*x[1], 4*x[1]-2*x[0]+2])
-
在另一些程序中,采用了数值法求梯度,数值法在参数比较少的情况下,效果较好;但参数极多的情况下,计算量非常大、运行效率极差;常用于检验所写梯度的正确性。
def num_gradient(f, x, delta = DELTA):
""" 计算数值梯度 """
# 中心差分
try:
params = x.tolist()
gradient_vector = np.array([(f(*(params[:i] + [params[i] + delta] + params[i + 1:]))
-f(*(params[:i] + [params[i] - delta] + params[i + 1:])))
/(2 * delta) for i in range(len(params))])
return gradient_vector
# 若中心差分报错(例如根号x在0处求梯度,定义域不可行),改为前向差分或后向差分
except:
try:
params = x.tolist()
gradient_vector = np.array([(f(*(params[:i] + [params[i] + delta] + params[i + 1:]))
-f(*(params[:])))
/delta for i in range(len(params))])
return gradient_vector
except:
params = x.tolist()
gradient_vector = np.array([(-f(*(params[:i] + [params[i] - delta] + params[i + 1:]))
+f(*(params[:])))
/delta for i in range(len(params))])
return gradient_vector
因此,我们对求梯度方法进行了改进,采用了解析法,也就是先计算出偏导表达式,再形成梯度向量。用解析法求得的梯度更为精确,效果更好。
dx = []
grad = []
for a in range(dimensions):
dx.append(diff(f, symbols('x'+str(a),real=True))) #偏导表达式,梯度
item={}
for i in range(dimensions):
item[x_[i]] = x[i]
grad.append(dx[a].subs(item)) #梯度值
若输入为:
f = 100*(x_[1]-x_[0]**2)**2+(1-x_[0])**2
x = [2,2]
d = [-1,-1]
print(dx)
print(grad)
输出为:
[-400*x0*(-x0**2 + x1) + 2*x0 - 2, -200*x0**2 + 200*x1] #偏导表达式,梯度
[1602,-400] #梯度值
4.3 统计维度方法的改进
- 在有些程序中,需要手动输入目标函数的维度:
n = int(input("please input the dimensions n="))
我们利用正则表达式来统计目标函数的维度:
import re
dimension_set = []
dimension_set = re.findall(r'x[0-9]\d*',str(f))
#统计x0,x1,...总数,'x[0-9]\d*'表示非负整数,\d表示单个数字字符,*表示前面的字符0次或多次
dimensions = len(set(dimension_set)) # 用set()去重
5 Python代码实现
import numpy as np
from sympy import *
import re
'''Wolfe-Powell非精确线性搜索,返回函数f在x处,方向d时的步长alpha'''
def WolfePowell(f,d,x,alpha_,rho,sigma):
maxk = 1000 #迭代上限
k = 0
phi0 = fun(x)
dimensions = dim(f_)
dx = []
grad = []
for a in range(dimensions):
dx.append(diff(f, symbols('x'+str(a),real=True))) #求偏导
item={}
for i in range(dimensions):
item[x_[i]] = x[i]
grad.append(dx[a].subs(item)) #求梯度
phi0_ = np.dot(grad,d)
print(dx)
print(grad)
a1 = 0
a2 = alpha_
alpha = (a1+a2)/2
phi1 = phi0
phi1_ = phi0_
k = 0;
for k in range(maxk): #限制迭代上限,避免时间太长
phi = fun(x + alpha * d)
if phi <= phi1 + rho * alpha * phi1_:
newx = x + alpha * d
newdx = []
newgrad = []
for a in range(dimensions):
newdx.append(diff(f, symbols('x'+str(a),real=True))) # 求偏导
newitem={}
for i in range(dimensions):
newitem[x_[i]] = newx[i]
newgrad.append(newdx[a].subs(newitem)) #求梯度
phi_ = np.dot(newgrad,d)
if phi_ >= sigma*phi0_:
break
else:
alpha_new = alpha + (alpha-a1) *phi_ / (phi1_-phi_)
a1 = alpha
alpha = alpha_new
phi1 = phi
phi1_ = phi_
else:
alpha_new = a1 + 0.5*(a1-alpha)**2*phi1_/((phi1-phi)-(a1-alpha)*phi1_)
a2 = alpha
alpha = alpha_new
k = k + 1
return alpha
'''利用正则表达式统计目标函数维度'''
def dim(f_):
dimension_set = []
dimension_set = re.findall(r'x[0-9]\d*',str(f_))
dimensions = len(set(dimension_set))
return dimensions
'''测试'''
x_ = []
for a in range(100):
x_.append(symbols('x'+str(a),real=True)) #设置符号变量
def fun(x_):
return 100*(x_[1]-x_[0]**2)**2+(1-x_[0])**2 #目标函数
f_ = fun(x_) #用于W-P
alpha_ = 1 #alpha_max
rho = 0.25 # rho∈(0,1/2)
sigma = 0.5 # sigma∈(rho,1)
x = np.random.rand(dim(f_)) #随机的初始点
d = np.array([-1,-1]) #初始方向
alpha=WolfePowell(f_,d,x,alpha_,rho,sigma)
print(alpha)
本文地址:https://blog.csdn.net/weixin_43761124/article/details/107436454
上一篇: Python画树状图