作业笔记：基于二次插值的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 概述

二次插值法（抛物线法）基本思路：在极小点附近，用二次三项式$\phi (x)$逼近目标函数$f(x)$

分为三点二次插值法和二点二次插值法：

• 三点二次插值法：已知三点的函数值，求极小点
• 二点二次插值法：
  • 已知两点的函数值，和其中一点的导数值，求极小点（内插）
  • 已知两点的导数值，求极小点（外插）

1.2 二点二次插值法

（1）内插：找两点之间的极小点

令二次三项式$\phi (x)=a+b(x-x_1)+c(x-x_1{)}^2$

已知$\phi (x_1)=f(x_1)$，${\phi }^{\prime }(x_1)=f^{\prime }(x_1)$，$\phi (x_2)=f(x_2)$，其中$f^{\prime }(x_1)<0$

为保证二次插值函数$\phi (x)$有极小点，要求$f(x_2)>f(x_1)+f^{\prime }(x_1)(x_2-x_1)$或者$f^{\prime }(x_2)>0$

得到以下方程组：
$$\{\begin{array}{c}\phi (x_1)=a=f(x_1)\\ {\phi }^{\prime }(x_1)=b=f^{\prime }(x_1)\\ \phi (x_2)=a+b(x_2-x_1)+c(x_2-x_1{)}^2=f(x_2)\end{array}$$
解方程组，得$a=f(x_1)$，$b=f^{\prime }(x_1)$，$c=\frac{f(x_2)-f(x_1)-f^{\prime }(x_1)(x_2-x_1)}{(x_2-x_1{)}^2}$。

为求$\phi (x)$的极小点，则令${\phi }^{\prime }(x)=b+2c(x-x_1)=0$ ，$\Rightarrow x=x_1-\frac{b}{2c}$，将$b,c$代入得极小点
$$\bar{x}=x_1-\frac{f^{\prime }(x_1)(x_2-x_1{)}^2}{2\left[f(x_2)-f(x_1)-f^{\prime }(x_1)(x_2-x_1)\right]}$$

（2）外插：找两点之外的极小点

令二次三项式$\phi (x)=a+b(x-x_2)+c(x-x_2{)}^2$

已知${\phi }^{\prime }(x_1)=f^{\prime }(x_1)$，${\phi }^{\prime }(x_2)=f^{\prime }(x_2)$

为保证二次插值函数$\phi (x)$有极小点，要求$f^{\prime }(x_1)<f^{\prime }(x_2)<0$

得到以下方程组：
$$\{\begin{array}{c}{\phi }^{\prime }(x_1)=b+2c(x_1-x_2)=f^{\prime }(x_1)\\ {\phi }^{\prime }(x_2)=b+2c(x_2-x_2)=f^{\prime }(x_2)\end{array}$$
解方程组，得$b=f^{\prime }(x_2)$，$c=\frac{f^{\prime }(x_1)-f^{\prime }(x_2)}{2(x_2-x_1)}$。

为求$\phi (x)$的极小点，则令${\phi }^{\prime }(x)=b+2c(x-x_2)=0$ ，$\Rightarrow x=x_2-\frac{b}{2c}$，将$b,c$代入得极小点
$$\bar{x}=x_2+\frac{f^{\prime }(x_2)(x_2-x_1)}{f^{\prime }(x_1)-f^{\prime }(x_2)}$$

2 准备知识：Wolfe Conditions

$$\begin{gathered} f(x_k+{\alpha }_k{d}_k)\le f(x_k)+\rho {\alpha }_{k}g(x_k{)}^T{d}_{k}condition(1) \\ g(x_{k+1}{)}^T{d}_k\ge \sigma g(x_k{)}^T{d}_{k}condition(2) \end{gathered}$$

其中$\rho \in (0,1/2),\sigma \in (\rho ,1)$

---

3 基于二点二次插值的Wolfe-Powell法

3.1 基本思路

• 满足$condition(1)$，且满足$condition(2)$，直接输出步长$\alpha$
• 满足$condition(1)$，但不满足$condition(2)$，用外插公式增大步长$\alpha$（见Step 3）
• 不满足$condition(1)$，用内插公式缩小步长$\alpha$（见Step 2）

3.2 计算步骤

Step 1.

选取初始数据，给定初始搜索区间$[0,{\alpha }_{\max }]$，给出参数$\rho \in (0,1/2),\sigma \in (\rho ,1)$，令${\alpha }_1=0,{\alpha }_2={\alpha }_{\max }$，取$\alpha \in ({\alpha }_1,{\alpha }_2)$，计算${\varphi }_1=f(x_k)$，${\varphi }_1^{\prime }=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_：${\alpha }_{\max }$，rho：$\rho$，sigma：$\sigma$。另外，a1：${\alpha }_1$，a2：${\alpha }_2$，alpha：$\alpha$，算法中取$\alpha =({\alpha }_1+{\alpha }_2)/2$。

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

phi1：${\varphi }_1$，phi1_：${\varphi }_1^{\prime }$，求梯度的方法会在后面说明。

Step 2.

计算$\varphi =\varphi (\alpha )=f(x_k+\alpha d_k)$。

若$\varphi \le {\varphi }_1+\rho \alpha {\varphi }_1^{\prime }$（条件1）成立，转到Step 3；否则，由二点二次插值法的内插公式计算$\bar{\alpha }$：
$$\bar{\alpha }={\alpha }_1-\frac{{\varphi }_1^{\prime }(\alpha -{\alpha }_1{)}^2}{2[\varphi -{\varphi }_1-{\varphi }_1^{\prime }(\alpha -{\alpha }_1)]}$$
令${\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.

计算${\varphi }^{\prime }={\varphi }^{\prime }(\alpha )=g(x_k+\alpha d_k{)}^T{d}_k$。

若${\varphi }^{\prime }\ge \sigma {\varphi }_1^{\prime }$（条件2）成立，则令${\alpha }_k=\alpha$，输出${\alpha }_k$，停止；否则，由二点二次插值法的外插公式计算$\bar{\alpha }$：
$$\bar{\alpha }=\alpha +\frac{{\varphi }^{\prime }(\alpha -{\alpha }_1)}{{\varphi }_1^{\prime }-{\varphi }^{\prime }}$$

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)$，所以要缩小步长$\lambda$。此时步长$\lambda$太大，取$\lambda$和步长下界$a$的中值为新步长，就能使步长缩小。因为$a\ge 0$，所以步长最多缩小到原来的$1/2$。但是，这种方法得到的新步长，并不一定能使目标函数最小。采用二次插值法的内插法，可以找到$[{\alpha }_1,\alpha ]$之间的极小点，能使目标函数近似最小，提高了算法效率。

• 转到第(4)步时，满足$condition(1)$，但不满足$condition(2)$，所以要增大步长$\lambda$。此时步长$\lambda$太小，取$\lambda$和步长上界$b$的中值为新步长，就能使步长放大。因为$\lambda =\min (2\lambda ,(\lambda +b)/2)$，步长最多放大到原来的$2$倍。但是，这种方法得到的新步长，并不一定能使目标函数最小。采用二次插值法的外插法，可以找到$[{\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)