【python】最优化方法之一维搜索（黄金分割法+斐波那契法）

1.概念

\qquad 一维搜索是最优化方法最简单的一种，即求一个在(a,b)内，连续下单峰函数$f(x)$的极小值。所谓下单峰函数就是只有一个极小值的函数。

2.遍历搜索

\qquad 这个问题看似简单，我们只需要指定步长，从区间左端点搜索到右端点，每次根据步长新产生一个值，比对它和前两个值的大小。记录三个值从左到由依次为$x_1$,$x_2$,$x_3$，若出现$x_1>x_2<x_3$的情况，则说明极小值点必定落在$（x_1,x_3）$之间。因此如果指定了最终区间不得超过eps的话，那么搜索步长不得大于0.5eps。
我们很快就会发现，这种遍历法虽然算法简单，但是迭代次数多，浪费资源严重，因此我们需要更快地得到最终区间的方法。在最优化方法这个数学分支中，可以证明斐波拉契法是压缩比（在同样迭代次数下开始区间和最终区间的比值）最高的方法，黄金分割比相比它压缩比小一些，但是算法会比较简单，占用计算资源也会略少。

3.优化算法

3.1.一维搜索原则

\qquad 首先在初始区间[a,b]中取两个点$x_1,x_2(x_1<x_2)$，计算$f(a),f(x_1),f(x_2),f(b)$的值。由于$f(x)$在(a,b)为下单峰函数，即有唯一的极小值，因此a或b必定不是极小值点。

1.若$f(x_1)<f(x_2)$则说明$(x_1,x_2)$我下降段，极小值$x^{\ast }$点应该在$x_1$的右侧，因此修改查找区间为$(x_1,b)$
2.若$f(x_1)>f(x_2)$则说明$(x_1,x_2)$我上升段，极小值$x^{\ast }$点应该在$x_2$的左侧，因此修改查找区间为$(a,x_2)$
一维搜索算法的核心就是解决取$x_1和x_2$的问题

3.2.黄金分割法

\qquad 设黄金分割比为k，取$x_1=b-(b-a)\ast k,x_2=a+(b-a)\ast k$。
在一维搜索中，如果满足条件1，那么需要在$(x_1,b)$取值。设$\rho =1-k$则
$\overline{CD}=1-\rho ,\overline{DE}=\rho ,\overline{EF}=1-\rho$由于在黄金分割比的作用下，$$\frac{1-\rho }{1}=\frac{1-2\rho }{1-\rho }$$，因此$$\frac{DE}{GH}=\frac{DE}{DF}=\frac{DF}{CF}=k\approx 0.618$$
因此若将E点平移至GH线段，将G的横坐标当做a，H的横坐标当做b，它就相当于是$x_1$的位置，因此只需要计算$x_2=a+(b-a)\ast k=b-(b-a)\ast \rho$即可
同理，如果满足条件2，$$\frac{DE}{IJ}=\frac{DE}{CE}=\frac{CE}{CF}=k\approx 0.618$$D点相当于IJ线段的$x_2$的位置，因此只需要计算$x_1=b-(b-a)\ast k=a+(b-a)\ast \rho$

Code Block

编写用户自定义函数的代码讲解参见以下链接：

[动态创建数学函数(CSDN])
以下为黄金分割法进行一维搜索的代码：

from math import *
import matplotlib.pyplot as plt
from pylab import *
# 通用函数f(x)靠用户录入
def function(x):
    fx = str_fx.replace("x", "%(x)f")  # 所有的"x"换为"%(x)function"
    return eval(fx % {"x": x})  # 字典类型的格式化字符串，将所有的"x"替换为变量x


# 绘图函数：给定闭区间（绘图间隔），绘图间隔默认为0.05，若区间较小，请自行修改
def drawf(a,b,interp=0.05):
    x = [a+ele*interp for ele in range(0, int((b-a)/interp))]
    y = [function(ele) for ele in x]
    plt.figure(1)
    plt.plot(x, y)
    xlim(a, b)
    title(init_str, color="b")
    plt.show()

# 黄金分割法进行一维搜索的函数
def gold_div_search(a,b,esp):
    data=list()
    x1=a+rou*(b-a)
    x2=b-rou*(b-a)
    data.append([a,x1,x2,b])
    while(b-a>esp):
        if function(x1)>function(x2):  #如果f(x1)>function(x2)，则在区间(x1,b)内搜索
            a=x1
            x1=x2
            x2=b-rou*(b-a)
            plt.plot(x2,function(x2),'r*')
        elif function(x1)<function(x2):  #如果f(x1)<function(x2),则在区间(a,x2)内搜索
            b=x2
            x2=x1
            x1=a+rou*(b-a)
            plt.plot(x1,function(x1),'r*')
        else:  #如果f(x1)=function(x2)，则在区间(x1,x2)内搜索
            a=x1
            b=x2
            x1=a+rou*(b-a)
            x2=b-rou*(b-a)
            plt.plot(x1,function(x1),'r*',x2,function(x2),'r*')
        data.append([a,x1,x2,b])
    with open("一维搜索（黄金分割法）.txt",mode="w",encoding="utf-8")as a_file: 
    # 保存的txt文件在程序的同目录下
        for i in range(0,len(data)):
            a_file.write("%d：\t"%(i+1))
            for j in range(0,4):
                a_file.write("function(%.3f)=%.3f\t"%(data[i][j],function(data[i][j])))
            a_file.write("\n")
    print("写入文件成功！")
    return [a,b]

rou = 1-(sqrt(5)-1)/2  # 1-rou为黄金分割比
init_str = input("请输入一个函数，默认变量为x：\n")  # 输入的最初字符串
para=input("请依次输入一维搜索的区间a,b和最终区间的精确值（用空格分隔）").split() # 导入区间
para=[float(ele) for ele in para]
a,b,esp=para
str_fx = init_str.replace("^", "**")  # 将所有的“^"替换为python的幂形式"**"
gold_div_search(a,b,esp)  # 调用黄金分割法并保存文件
drawf(a,b,(b-a)/2000)  # 绘制函数图形

以下为交互式界面的效果：（请在Shell或Pycharm中的Python Console里面运行）

    请输入一个函数，默认变量为x：
    atan(x)-log(2*x+1)+4*x**2-3*x
    请依次输入一维搜索的区间a,b和最终区间的精确值（用空格分隔）0 1 0.001
    写入文件成功！

文本文档数据：（如用此程序运行，保存的文本文档与程序同目录）

1： f(0.000)=0.0000000 f(0.382)=-0.7649875 f(0.618)=-0.5773825 f(1.000)=0.6867859
2： f(0.000)=0.0000000 f(0.236)=-0.6401822 f(0.382)=-0.7649875 f(0.618)=-0.5773825
3： f(0.236)=-0.6401822 f(0.382)=-0.7649875 f(0.472)=-0.7485370 f(0.618)=-0.5773825
4： f(0.236)=-0.6401822 f(0.326)=-0.7399126 f(0.382)=-0.7649875 f(0.472)=-0.7485370
5： f(0.326)=-0.7399126 f(0.382)=-0.7649875 f(0.416)=-0.7669244 f(0.472)=-0.7485370
6： f(0.382)=-0.7649875 f(0.416)=-0.7669244 f(0.438)=-0.7630202 f(0.472)=-0.7485370
7： f(0.382)=-0.7649875 f(0.403)=-0.7673941 f(0.416)=-0.7669244 f(0.438)=-0.7630202
8： f(0.382)=-0.7649875 f(0.395)=-0.7669382 f(0.403)=-0.7673941 f(0.416)=-0.7669244
9： f(0.395)=-0.7669382 f(0.403)=-0.7673941 f(0.408)=-0.7673906 f(0.416)=-0.7669244
10： f(0.395)=-0.7669382 f(0.400)=-0.7672874 f(0.403)=-0.7673941 f(0.408)=-0.7673906
11： f(0.400)=-0.7672874 f(0.403)=-0.7673941 f(0.405)=-0.7674185 f(0.408)=-0.7673906
12： f(0.403)=-0.7673941 f(0.405)=-0.7674185 f(0.406)=-0.7674176 f(0.408)=-0.7673906
13： f(0.403)=-0.7673941 f(0.404)=-0.7674129 f(0.405)=-0.7674185 f(0.406)=-0.7674176
14： f(0.404)=-0.7674129 f(0.405)=-0.7674185 f(0.406)=-0.7674196 f(0.406)=-0.7674176
15： f(0.405)=-0.7674185 f(0.406)=-0.7674196 f(0.406)=-0.7674194 f(0.406)=-0.7674176
16： f(0.405)=-0.7674185 f(0.405)=-0.7674193 f(0.406)=-0.7674196 f(0.406)=-0.7674194

可以看出我们的最终区间为(0.405,0.406)，精度为0.001.
黄金分割法查找的点会被记录在matplotlib绘制的函数图像上，一目了然，可以看出，在最终区间精度为0.001的情况下，我们的算法收敛速度很快。


\qquad 绘图可以看出函数在极小点附近非常平坦，极小值点差不多在0.405的位置，与文本文档的数据一致。

3.3.斐波拉契法

\qquad 斐波拉契法是一维搜索中压缩比最高的搜索算法。斐波拉契法基于斐波拉契数列产生比例值，斐波拉契数列 { F n } \lbrace F_n \rbrace {Fn​}的定义如下：
$$F_1=F_2=1,F_{n+2}=F_{n+1}+F_n$$
取分点时，$x_1=a+\frac{F_{n-2}}{F_n}(b-a),x_2=a+\frac{F_{n-1}}{F_n}(b-a)$

$$\{\begin{array}{l}{x}_1-a=\frac{F_{n-2}}{F_n}(b-a)\\ x_2-a=\frac{F_{n-1}}{F_n}(b-a)\\ x_2-x_1=\frac{F_{n-1}-F_{n-2}}{F_n}(b-a)=\frac{F_{n-3}}{F_n}(b-a)\\ b-x_2=\frac{F_n-F_{n-1}}{F_n}(b-a)=\frac{F_{n-2}}{F_n}(b-a)\end{array}$$

1.若$x_1＜x_2$，则取新区间为$[a,x_2]$，则$$\begin{gathered} \overline{DE}=\frac{F_{n-3}}{F_n}(b-a)，\overline{IJ}=\frac{F_{n-1}}{F_n}(b-a),\overline{CD}=\frac{F_{n-2}}{F_n}(b-a), \\ \frac{CD}{IJ}=\frac{F_{n-2}}{F_{n-1}}，因此D相当于IJ中的x_2的位置 \\ 只需要再取x_1=a+\frac{F_{n-3}}{F_{n-1}}(b-a)即可 \end{gathered}$$

2.若$x_1>x_2$，则取新区间为$[x_1,b]$，则$$\begin{gathered} \overline{DE}=\frac{F_{n-3}}{F_n}(b-a)，\overline{GH}=\frac{F_{n-1}}{F_n}(b-a),\overline{CD}=\frac{F_{n-2}}{F_n}(b-a), \\ \frac{DE}{GH}=\frac{F_{n-3}}{F_{n-1}}，因此D相当于GH中的x_1的位置 \\ 只需要再取x_2=a+\frac{F_{n-2}}{F_{n-1}}(b-a)即可 \end{gathered}$$

对于要求精度esp，取压缩比 C = b − a e s p ， F n = m i n { x ∣ x ≥ C , x ∈ N + } C=\frac{b-a}{esp}，F_n=min\lbrace x|x≥C,x∈N^+\rbrace C=espb−a​，Fn​=min{x∣x≥C,x∈N+}，按照上述方法迭代至分母为$F_2$为止，最后按照以下方法进行：

Code Block

首先我们需要一个生产斐波那契数列（Fabonaci）的函数：

# 在本案例中，使用一次性生成斐波拉契列表的方法算法复杂度更低
# 生成的斐波拉契数列F(0)=0,F(1)=F(2)=1,剩余项和普通斐波拉契数列类似
def Fab_list(Fmax):
    a,b=0,1
    Fablist = [a,b]  # 返回一个列表
    while Fablist[-1]< Fmax:
        a,b = b,a+b
        Fablist.append(b)
    return Fablist  

再按照斐波那契法迭代编写算法即可：
通用函数的编写仍然参考以下链接：
[动态创建数学函数(CSDN)]
注：需要安装matplotlib模块、pylab模块和Pillow模块，版本为python3

步骤：
Windows键+R，在【运行】窗口输入cmd并回车，在命令提示行下输入：
pip install matplotlib
pip install pylab
pip install Pillow
即可

from math import *
import matplotlib.pyplot as plt  # 绘图模块
from pylab import *  # 绘图辅助模块



# 通用函数f(x)靠用户录入
def function(x):
    fx = str_fx.replace("x", "%(x)f")  # 所有的"x"换为"%(x)function"
    return eval(fx % {"x": x})  # 字典类型的格式化字符串，将所有的"x"替换为变量x


# 绘图函数：给定闭区间（绘图间隔），绘图间隔默认为0.05，若区间较小，请自行修改
def drawf(a, b, interp=0.05):
    x = [a + ele * interp for ele in range(0, int((b - a) / interp))]
    y = [function(ele) for ele in x]
    plt.figure(1)
    plt.plot(x, y)
    xlim(a, b)
    title(init_str, color="b")  # 标注函数名
    plt.show()
    return "绘图完成！"


# 在本案例中，使用一次性生成斐波拉契列表的方法算法复杂度更低
# 生成的斐波拉契数列F(0)=0,F(1)=F(2)=1,剩余项和普通斐波拉契数列类似
def Fab_list(Fmax):
    a,b=0,1
    Fablist = [a,b]  # 返回一个列表
    while Fablist[-1]< Fmax:
        a,b = b,a+b
        Fablist.append(b)
    Fablist.pop()  # 舍掉最后一个元素
    return Fablist    


def Fabonaci_search(a, b, Fab):
    n = len(Fab)-1  # 获取斐波那契数列的长度
    if n < 3:
        return "精度过低，无法进行斐波那契一维搜索"
    data = list()
    data.append([a, b])
    x1 = a + Fab[n - 2]/Fab[n] * (b - a)
    x2 = a + Fab[n - 1]/Fab[n] * (b - a)
    t = n
    while (t > 3):
        if function(x1) > function(x2):  # 如果f(x1)>f(x2)，则在区间(x1,b)内搜索
            a = x1
            x1 = x2
            x2 = a + Fab[t - 1]/Fab[t] * (b - a)
            plt.plot(x2, function(x2), 'r*')
            data.append([x1, b])
        elif function(x1) < function(x2):  # 如果f(x1)<(x2),则在区间(a,x2)内搜索
            b = x2
            x2 = x1
            x1 = a + Fab[t - 2]/Fab[t]* (b - a)
            plt.plot(x1, function(x1), 'r*')
            data.append([a, x2])
        else:  # 如果f(x1)=function(x2)，则在区间(x1,x2)内搜索
            a = x1
            b = x2
            x1 = a + Fab[t - 2]/Fab[t] * (b - a)
            x2 = a + Fab[t - 1]/Fab[t] * (b - a)
            plt.plot(x1, function(x1), 'r*', x2, function(x2), 'r*')
        data.append([x1, x2])
        t -= 1
    x1 = a + 0.5 * (b - a)  # 斐波那契数列第3项和第2项的比
    x2 = x1 + 0.1 * (b - a)  # 偏离一定距离，人工构造的点
    if function(x1) > function(x2):  # 如果f(x1)>function(x2)，则在区间(x1,b)内搜索
        plt.plot(x2, function(x2), 'r*')
        data.append([x1, b])
        a = x1
    elif function(x1) < function(x2):  # 如果f(x1)<function(x2),则在区间(a,x2)内搜索
        plt.plot(x1, function(x1), 'r*')
        data.append([a, x2])
        b = x2
    else:  # 如果f(x1)=function(x2)，则在区间(x1,x2)内搜索
        plt.plot(x1, function(x1), 'r*', x2, function(x2), 'r*')
        data.append([x1, x2])
    with open(r"C:\Users\ouni\桌面\一维搜索（黄金分割法）.txt", mode="w", encoding="utf-8")as a_file:
        for i in range(0, len(data)):
            a_file.write("%d：\t" % (i + 1))
            for j in range(0, 2):
                a_file.write("f(%.3f)=%.7f\t" % (data[i][j], function(data[i][j])))
            a_file.write("\n")
    print("写入文件成功！")
    return [a, b]


init_str = input("请输入一个函数，默认变量为x：\n")  # 输入的最初字符串
para = input("请依次输入一维搜索的区间a,b和最终区间的精确值（用空格分隔）").split()  # 导入区间
para = [float(ele) for ele in para]  # 将输入的字符串转换为浮点数
low, high, esp = para  # 输入参数列表（最小值、最大值和最终精度）
str_fx = init_str.replace("^", "**")  # 将所有的“^"替换为python的幂形式"**"
print(Fabonaci_search(low, high, Fab_list(int(ceil((high-low)/esp)))))  # 传入区间和斐波拉契列表
drawf(low, high, (high - low) / 2000)  # 默认精度是2000个点

以下代码请在交互式界面运行：

请输入一个函数，默认变量为x：
e^(-x)+x^2
请依次输入一维搜索的区间a,b和最终区间的精确值（用空格分隔）0 1 0.01
写入文件成功！
[0.3402663805075117, 0.35279320792829183]
关闭绘图框以继续交互式命令

文本文档如下：

1： f(0.000)=1.0000000 f(1.000)=1.3678794
2： f(0.000)=1.0000000 f(0.382)=0.8284208
3： f(0.236)=0.8454509 f(0.382)=0.8284208
4： f(0.382)=0.8284208 f(0.618)=0.9209301
5： f(0.382)=0.8284208 f(0.472)=0.8465897
6： f(0.236)=0.8454509 f(0.382)=0.8284208
7： f(0.326)=0.8280571 f(0.382)=0.8284208
8： f(0.236)=0.8454509 f(0.326)=0.8280571
9： f(0.292)=0.8320851 f(0.326)=0.8280571
10：f(0.326)=0.8280571 f(0.382)=0.8284208
11： f(0.326)=0.8280571 f(0.347)=0.8272109
12： f(0.347)=0.8272109 f(0.382)=0.8284208
13： f(0.347)=0.8272109 f(0.361)=0.8273036
14： f(0.326)=0.8280571 f(0.347)=0.8272109
15： f(0.340)=0.8273620 f(0.347)=0.8272109
16： f(0.347)=0.8272109 f(0.361)=0.8273036
17： f(0.347)=0.8272109 f(0.354)=0.8271921
18： f(0.340)=0.8273620 f(0.353)=0.8271855
}

图形如下：


可以看出函数的极小值确实在0.35附近，与文本文档标定的一致。
注：以上程序请在Shell或Cmd中运行，如需当做模块运行，请在主要程序的前面加上if __name__=="__main__"。