欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页  >  IT编程

作业笔记:基于二次插值的Wolfe-Powell非精确线搜索算法及Python代码实现

程序员文章站 2022-10-04 08:39:49
1 准备知识:二次插值法1.1 概述**二次插值法(抛物线法)基本思路:**在极小点附近,用二次三项式φ(x)\varphi(x)φ(x)逼近目标函数f(x)f(x)f(x)分为三点二次插值法和二点二次插值法:三点二次插值法:已知三点的函数值,求极小点二点二次插值法:已知两点的函数值,和其中一点的导数值,求极小点(内插)已知两点的导数值,求极小点(外插)1.2 二点二次插值法(1)内插令二次三项式φ(x)=a+b(x−x1)+c(x−x1)2\varphi(x)=a+b(x-x...

写在前面的说明:

  • 查阅了很多资料,发现资料里对于非精确线搜索的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(xx1)+c(xx1)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)(x2x1)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(x2x1)+c(x2x1)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)(x2x1)(x2x1)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(xx1)=0\varphi'(x)=b+2c(x-x_1)=0x=x1b2c\Rightarrow x=x_1-\frac{b}{2c},将b,cb,c代入得极小点
xˉ=x1f(x1)(x2x1)22[f(x2)f(x1)f(x1)(x2x1)] \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(xx2)+c(xx2)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(x1x2)=f(x1)φ(x2)=b+2c(x2x2)=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(x2x1)c=\frac{f'(x_1)-f'(x_2)}{2(x_2-x_1)}

为求φ(x)\varphi(x)的极小点,则令φ(x)=b+2c(xx2)=0\varphi'(x)=b+2c(x-x_2)=0x=x2b2c\Rightarrow x=x_2-\frac{b}{2c},将b,cb,c代入得极小点
xˉ=x2+f(x2)(x2x1)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ρ\rhosigmaσ\sigma。另外,a1α1\alpha_1a2α2\alpha_2alphaα\alpha,算法中取α=(α1+α2)/2\alpha=(\alpha_1+\alpha_2)/2

phi1 = fun(x) 
phi1_ = np.dot(grad,d)

phi1ϕ1\phi_1phi1_ϕ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程序中,采用如下的方法计算步长α
    作业笔记:基于二次插值的Wolfe-Powell非精确线搜索算法及Python代码实现

本文的算法中,二点二次插值法的内插和外插,分别对应上面第(3)步和第(4)步。但不同之处在于对步长放缩的方式:

  • 转到第(3)步时,不满足condition(1)condition(1),所以要缩小步长λ\lambda。此时步长λ\lambda太大,取λ\lambda和步长下界aa的中值为新步长,就能使步长缩小。因为a0a\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

相关标签: 算法 python