写在前面的说明:
- 查阅了很多资料,发现资料里对于非精确线搜索的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 概述
二次插值法(抛物线法)基本思路:在极小点附近,用二次三项式φ(x)\varphi(x)逼近目标函数f(x)f(x)
分为三点二次插值法和二点二次插值法:
- 三点二次插值法:已知三点的函数值,求极小点
- 二点二次插值法:
- 已知两点的函数值,和其中一点的导数值,求极小点(内插)
- 已知两点的导数值,求极小点(外插)
1.2 二点二次插值法
(1)内插:找两点之间的极小点
令二次三项式φ(x)=a+b(x−x1)+c(x−x1)2\varphi(x)=a+b(x-x_1)+c(x-x_1)^2
已知φ(x1)=f(x1)\varphi(x_1)=f(x_1),φ′(x1)=f′(x1)\varphi'(x_1)=f'(x_1),φ(x2)=f(x2)\varphi(x_2)=f(x_2),其中f′(x1)<0f'(x_1)<0
为保证二次插值函数φ(x)\varphi(x)有极小点,要求f(x2)>f(x1)+f′(x1)(x2−x1)f(x_2)>f(x_1)+f'(x_1)(x_2-x_1)或者f′(x2)>0f'(x_2)>0
得到以下方程组:
{φ(x1)=a=f(x1)φ′(x1)=b=f′(x1)φ(x2)=a+b(x2−x1)+c(x2−x1)2=f(x2)
\left\{\begin{matrix} \varphi(x_1)=a=f(x_1)
\\ \varphi'(x_1)=b=f'(x_1)
\\ \varphi(x_2)=a+b(x_2-x_1)+c(x_2-x_1)^2=f(x_2)
\end{matrix}\right.
解方程组,得a=f(x1)a=f(x_1),b=f′(x1)b=f'(x_1),c=f(x2)−f(x1)−f′(x1)(x2−x1)(x2−x1)2c=\frac{f(x_2)-f(x_1)-f'(x_1)(x_2-x_1)}{(x_2-x_1)^2}。
为求φ(x)\varphi(x)的极小点,则令φ′(x)=b+2c(x−x1)=0\varphi'(x)=b+2c(x-x_1)=0 ,⇒x=x1−b2c\Rightarrow x=x_1-\frac{b}{2c},将b,cb,c代入得极小点
xˉ=x1−f′(x1)(x2−x1)22[f(x2)−f(x1)−f′(x1)(x2−x1)]
\bar x=x_1-\frac{f'(x_1)(x_2-x_1)^2}{2\left [f(x_2)-f(x_1)-f'(x_1)(x_2-x_1)\right ]}
(2)外插:找两点之外的极小点
令二次三项式φ(x)=a+b(x−x2)+c(x−x2)2\varphi(x)=a+b(x-x_2)+c(x-x_2)^2
已知φ′(x1)=f′(x1)\varphi'(x_1)=f'(x_1),φ′(x2)=f′(x2)\varphi'(x_2)=f'(x_2)
为保证二次插值函数φ(x)\varphi(x)有极小点,要求f′(x1)<f′(x2)<0f'(x_1)<f'(x_2)<0
得到以下方程组:
{φ′(x1)=b+2c(x1−x2)=f′(x1)φ′(x2)=b+2c(x2−x2)=f′(x2)
\left\{\begin{matrix} \varphi'(x_1)=b+2c(x_1-x_2)=f'(x_1)
\\ \varphi'(x_2)=b+2c(x_2-x_2)=f'(x_2)
\end{matrix}\right.
解方程组,得b=f′(x2)b=f'(x_2),c=f′(x1)−f′(x2)2(x2−x1)c=\frac{f'(x_1)-f'(x_2)}{2(x_2-x_1)}。
为求φ(x)\varphi(x)的极小点,则令φ′(x)=b+2c(x−x2)=0\varphi'(x)=b+2c(x-x_2)=0 ,⇒x=x2−b2c\Rightarrow x=x_2-\frac{b}{2c},将b,cb,c代入得极小点
xˉ=x2+f′(x2)(x2−x1)f′(x1)−f′(x2)
\bar x=x_2+\frac{f'(x_2)(x_2-x_1)}{f'(x_1)-f'(x_2)}
2 准备知识:Wolfe Conditions
f(xk+αkdk)≤f(xk)+ραkg(xk)Tdk condition(1)g(xk+1)Tdk≥σg(xk)Tdk condition(2)
f(x_k+\alpha_k d_k)\leq f(x_k)+\rho \alpha_k g(x_k)^T d_k \ \ \ \ \ condition(1) \\
g(x_{k+1})^T d_k \geq \sigma g(x_k)^T d_k \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ condition(2)
其中ρ∈(0,1/2),σ∈(ρ,1)\rho \in(0,1/2),\sigma \in (\rho,1)
3 基于二点二次插值的Wolfe-Powell法
3.1 基本思路
- 满足condition(1)condition(1),且满足condition(2)condition(2),直接输出步长α\alpha
- 满足condition(1)condition(1),但不满足condition(2)condition(2),用外插公式增大步长α\alpha(见Step 3)
- 不满足condition(1)condition(1),用内插公式缩小步长α\alpha(见Step 2)
3.2 计算步骤
Step 1.
选取初始数据,给定初始搜索区间[0,αmax][0,\alpha_{\max}],给出参数ρ∈(0,1/2),σ∈(ρ,1)\rho \in(0,1/2),\sigma \in (\rho,1),令α1=0,α2=αmax\alpha_1=0,\alpha_2=\alpha_{\max},取α∈(α1,α2)\alpha \in(\alpha_1,\alpha_2),计算ϕ1=f(xk)\phi_1=f(x_k),ϕ1′=g(xk)Tdk\phi_1'=g(x_k)^T d_k。
def WolfePowell(f,d,x,alpha_,rho,sigma):
a1 = 0
a2 = alpha_
alpha = (a1+a2)/2
输入的初始数据有,f
:目标函数,d
:初始方向,x
:初始点,alpha_
:αmax\alpha_{\max},rho
:ρ\rho,sigma
:σ\sigma。另外,a1
:α1\alpha_1,a2
:α2\alpha_2,alpha
:α\alpha,算法中取α=(α1+α2)/2\alpha=(\alpha_1+\alpha_2)/2。
phi1 = fun(x)
phi1_ = np.dot(grad,d)
phi1
:ϕ1\phi_1,phi1_
:ϕ1′\phi_1',求梯度的方法会在后面说明。
Step 2.
计算ϕ=ϕ(α)=f(xk+αdk)\phi=\phi(\alpha)=f(x_k+\alpha d_k)。
若ϕ≤ϕ1+ραϕ1′\phi\leq\phi_1+\rho\alpha\phi_1'(条件1)成立,转到Step 3;否则,由二点二次插值法的内插公式计算αˉ\bar \alpha:
αˉ=α1−ϕ1′(α−α1)22[ϕ−ϕ1−ϕ1′(α−α1)]
\bar \alpha=\alpha_1-\frac{\phi_1'(\alpha-\alpha_1)^2}{2[\phi-\phi_1-\phi_1'(\alpha-\alpha_1)]}
令α2=α\alpha_2=\alpha,α=αˉ\alpha=\bar \alpha,转到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.
计算ϕ′=ϕ′(α)=g(xk+αdk)Tdk\phi'=\phi'(\alpha)=g(x_k+\alpha d_k)^T d_k。
若ϕ′≥σϕ1′\phi'\geq\sigma\phi_1'(条件2)成立,则令αk=α\alpha_k=\alpha,输出αk\alpha_k,停止;否则,由二点二次插值法的外插公式计算αˉ\bar \alpha:
αˉ=α+ϕ′(α−α1)ϕ1′−ϕ′
\bar \alpha=\alpha+\frac{\phi'(\alpha-\alpha_1)}{\phi_1'-\phi'}
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)步时,不满足condition(1)condition(1),所以要缩小步长λ\lambda。此时步长λ\lambda太大,取λ\lambda和步长下界aa的中值为新步长,就能使步长缩小。因为a≥0a\geq0,所以步长最多缩小到原来的1/21/2。但是,这种方法得到的新步长,并不一定能使目标函数最小。采用二次插值法的内插法,可以找到[α1,α][\alpha_1,\alpha]之间的极小点,能使目标函数近似最小,提高了算法效率。
-
转到第(4)步时,满足condition(1)condition(1),但不满足condition(2)condition(2),所以要增大步长λ\lambda。此时步长λ\lambda太小,取λ\lambda和步长上界bb的中值为新步长,就能使步长放大。因为λ=min(2λ,(λ+b)/2)\lambda=\min{(2\lambda,(\lambda+b)/2)},步长最多放大到原来的22倍。但是,这种方法得到的新步长,并不一定能使目标函数最小。采用二次插值法的外插法,可以找到[α1,α][\alpha_1,\alpha]之外(大于α\alpha一侧)的极小点,能使目标函数近似最小,提高了算法效率。
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