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 概述

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

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

• 三点二次插值法：已知三点的函数值，求极小点
• 二点二次插值法：

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

1.2 二点二次插值法

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

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

已知$\varphi(x_1)=f(x_1)$，$\varphi$，$\varphi(x_2)=f(x_2)$，其中$f$

为保证二次插值函数$\varphi(x)$有极小点，要求$f(x_2)>f(x_1)+f$或者$f$

得到以下方程组：
$\left\{\begin{matrix} \varphi(x_1)=a=f(x_1) \\ \varphi$
解方程组，得$a=f(x_1)$，$b=f$，$c=\frac{f(x_2)-f(x_1)-f$。

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

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

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

已知$\varphi$，$\varphi$

为保证二次插值函数$\varphi(x)$有极小点，要求$f$

得到以下方程组：
$\left\{\begin{matrix} \varphi$
解方程组，得$b=f$，$c=\frac{f$。

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

2 准备知识：Wolfe Conditions

$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)$

其中$\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)$，计算$\phi_1=f(x_k)$，$\phi_1$。

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：$\phi_1$，phi1_：$\phi_1$，求梯度的方法会在后面说明。

Step 2.

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

若$\phi\leq\phi_1+\rho\alpha\phi_1$（条件1）成立，转到Step 3；否则，由二点二次插值法的内插公式计算$\bar \alpha$：
$\bar \alpha=\alpha_1-\frac{\phi_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.

计算$\phi$。

若$\phi$（条件2）成立，则令$\alpha_k=\alpha$，输出$\alpha_k$，停止；否则，由二点二次插值法的外插公式计算$\bar \alpha$：
$\bar \alpha=\alpha+\frac{\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)$，所以要缩小步长$\lambda$。此时步长$\lambda$太大，取$\lambda$和步长下界$a$的中值为新步长，就能使步长缩小。因为$a\geq0$，所以步长最多缩小到原来的$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)