三种非导数区间消去法理论推导和优劣对比-Matlab

        在处理连续函数最小值问题中，有很多种方法可以求得一维函数最小值点，并且收敛速度很快。但大多都需要函数的导数信息，以作为迭代信息使用，而在考虑连续不可导函数的问题时，这些方法就不是特别好用，需要做一些处理（如近似导数、割线法等）。而有一类区间消去法，仅仅使用函数连续的性质，即可求出函数最小值点的近似解。本文将介绍三种不用导数信息，而仅仅使用函数本身函数值来压缩最小值区间求得近似解，并给出理论上证明和基于Matlab实验上的测试。并且是在假设函数是一个单谷函数的情况下来讨论三种方法。

一：方法介绍

        二分法：最容易实现的区间分割法，理论上可以证明收敛于最小值点，且代码复杂度不高，是最常用的一个区间消去法。二分法本质基于贪婪算法，每一步取点都是当前步骤的最优选择，这使得二分法在不限制使用函数值次数的情况下，收敛速度比之其他两种方法最快。但贪婪算法的缺点很明显，在有限次数的情况下，如果只考虑眼前最优而非全局最优，最终结果会丧失部分精度。其原因在于二分法在第k步使用的函数值不能作为第k+1步的函数值使用，会浪费掉上一步求得的函数值，在实际科研中，往往取得一个实验数据的难度很大，所以在有限的情况下，二分法是较劣的方法。

        黄金分割法(0.618法):在二分法的基础上加以改进，其改进思路是第k步的函数值可以作为第k+1步的函数值来使用，这样保证了函数值能够被充分利用。在有限步数的情况下，黄金分割法表现得比二分法更好。例如，如果仅计算10次函数值，若考虑二者运气最差的情况(即每次消去的区间都是较少的那部分，二分法每次都是消去1/2)，则二分法只能迭代5次，最后得到的区间缩减为原来的(1/2)^5,而黄金分割法除了第一步需要使用2个函数值，其余步骤都只需要一个函数值，故而区间缩减为原来的(2/3)^9，其缩减区间长度更多。故在有限步数(函数值)情况下，黄金分割法优于二分法。

        斐波那契分割法：黄金分割法是从每一步都为下一步着想来设计的算法，但还可以有一些改进空间。因为对于最后一个函数值，放在黄金分割点的位置为下一步考虑，但实际上我们已经没有下一步了，所以这个处理方法显然不太可行。而斐波那契法可以考虑到全局，即既可以每一步都为下一步做准备，也能够在最后一步做出最优的选择，即放在中点附近消去1/2的长度。从理论上来推导，仅用10个函数值，斐波那契分割法可以将区间缩减到原来的1/89，这个缩减的区间长度是要比黄金分割法和二分法都要多的。可以证明，在有限步数的情况下，斐波那契分割法是最优的区间分割法，没有方法比它更好。

二、原理解释

以下用手写推导简要写出了二分法，黄金分割法和斐波那契分割法的过程。

1.二分法

2.黄金分割法

斐波那契分割法

三、代码实现

1.二分法:

   

function [f1]=Dichotomy(f,a,b,expect)
        %y0,y1，画图时作为初始区间。因为后续a，b会改变，所以提前存入变量值
        %delta是a，b之间的长度
        %二分法不使用导数，使用中点附近的函数值做差商，近似于导数来判断消去的区间，故而每一次迭代都需要两个函数值
        y0=a;
        y1=b;
        delta=b-a;
        h=0.1*expect*delta;   %h是选取中点附近点的距离，采用从中点往左右两边走h步长的方法近似导数来判断要消去的区间
        tic
        count=0;    %count用于计数，判断二分法要达到所给精度需要用的函数值个数
        while abs(b-a)>expect*delta    %当b-a长度大于用户期待的区间长度时，进入循环
            c=(a+b)/2;
            if f(c+h)-f(c-h)>0      %类似于导数，如果大于零，则消去[c,b]部分，留下[a,c]部分
                b=c;
            elseif f(c+h)-f(c-h)<0  %同上述
                a=c;
            else
                break;
            end
            count=count+2;  %每次迭代都需要用两个函数值，故次数加二
        end
        toc
        c=(a+b)/2;%取最终得到的区间的中点值作为近似得到的零点
        disp(['用了' num2str(count) '次函数值'])
        disp(['函数极小值点是:' num2str(c)]);
        disp(['区间缩减精度为原来的:' num2str(100*abs(b-a)/delta) '%']);  %打印出相比于最初的区间，区间缩短了多少
        disp(['与函数真正极小值点距离:' num2str(abs(fminbnd(f,a,b)-c))]);  %用内置函数求出最小值点，看二分法是否收敛于最小值点
        f1=f(c);
        xx=y0:0.01:y1;  %用图像直观显示出所求点的位置
        y1=f(xx);
        plot(xx,y1,'b-');
        hold on
        plot(c,f(c),'r*')
        legend 函数曲线 近似极小值点
        hold off
    end

2.黄金分割法：

function [f2]=Golden(f,a,b,expect)  %此为黄金分割法函数
        y0=a;   %此处作用跟二分法一致
        y1=b;
        delta=b-a;  %区间初始长度
        x1=a+0.382*abs(b-a);    %x1为区间左侧的黄金分割点
        x2=a+0.618*abs(b-a);    %x2为区间右侧的黄金分割点
        tic    %计时查看函数效率
        count=1;    %黄金分割法除最开始使用两个函数点，其余每次都只用一个点，经测试count＝1时与实际符合
        while abs(b-a)>delta*expect    %当区间长度大于预期长度，进入循环
            if f(x2)-f(x1)>0    %此时消去[x2,b]，保留[a,x1]，此时原先的x1变为x2，x1另外计算，故每次只用计算一个函数值
                b=x2;
                x2=x1;
                x1=a+0.382*(b-a);
            elseif f(x2)-f(x1)<0
                a=x1;
                x1=x2;
                x2=a+0.618*(b-a);
            else
                break;
            end
            count=count+1;  %每次迭代只多用了一个函数值
        end
        c=(a+b)/2;  %取a，b中点作为零点近似值
        f2=f(c);
        toc
        disp(['用了' num2str(count) '次函数值'])
        disp(['函数极小值点是:' num2str(c)]);
        disp(['区间缩减为原来的:' num2str(100*abs(b-a)/delta) '%']);
        disp(['与函数真正极小值点距离:' num2str(abs(fminbnd(f,0,2)-c))]);
        g=y0:0.01:y1;
        y2=f(g);
        plot(g,y2);
        hold on
        plot(c,f(c),'r*')
        legend 函数曲线 近似极小值点
        hold off
    end

3.斐波那契分割法：

    function[f3]=Fibonacci(f,a,b,expect)    %使用斐波那契法求函数最小值
        epsilon=0.05;   %设置迭代最后一步的步长epsilon，避免x0,x1都取中点
        y0=a;
        y1=b;
        Fib_num=2;  %构造斐波那契数列
        F=[1,1];
        delta=b-a;
        while (1+2*epsilon)/F(Fib_num)>expect    %当斐波那契数列中最后一个数的倒数大于用户输入精度，进入循环
            F(Fib_num+1)=F(Fib_num)+F(Fib_num-1);
            Fib_num=Fib_num+1;  
        end
        F(1)=0.9;   %设置F(1),F(2)不为1，目的是偏离中点，取中点附近的值来比较，避免两个点重合，所以斐波那契法会丢失一点点分割精度
        F(2)=1.1;
        x1=a+F(Fib_num-2)/F(Fib_num)*(b-a);    %确定斐波那契左侧点x1
        x2=a+F(Fib_num-1)/F(Fib_num)*(b-a);    %确定斐波那契右侧点x2
        count=1;
        tic
        while Fib_num>=3
            if f(x2)>f(x1)  %左侧函数值小于右侧函数值，此时消去[x2,b]区间，原来x1的位置变为x2,x1需要另算，这样每次只计算一次函数值
                Fib_num=Fib_num-1;  %每次迭代，斐波那契数列往前推一个数，直到下标3结束
                b=x2;
                x2=x1;
                x1=a+F(Fib_num-2)/F(Fib_num)*(b-a);  %x1的位置由此刻斐波那契数的前面的第二个数决定，x2的位置由斐波那契数前面的第一个数决定
            else
                Fib_num=Fib_num-1;
                a=x1;
                x1=x2;
                x2=a+F(Fib_num-1)/F(Fib_num)*(b-a);
            end
            %disp(['第' num2str(count) '次迭代区间缩减为原来的' num2str((b-a)/4*100) '%'])
            count=count+1;
        end
        toc
        c=(a+b)/2;
        disp(['用了' num2str(count) '次函数值'])
        disp(['函数极小值点是:' num2str(c)]);
        disp(['区间缩减为原来的:' num2str(100*abs(b-a)/delta) '%']);
        disp(['与函数真正极小值点距离:' num2str(abs(fminbnd(f,0,2)-c))]);
        f3=f(c);
        cci=y0:0.01:y1;
        y3=f(cci);
        plot(cci,y3,'b-');
        hold on
        plot(c,f(c),'r*')
        legend 函数曲线 近似极小值点
        hold off
    end

各代码运行结果图示：

1.二分法

2.黄金分割法

3.斐波那契分割法(理论上使用10次函数值，未加修正因子epsilon的结果，修正因子解释见斐波那契代码注释)

4.斐波那契分割法(实际上使用11次函数值，因为最后一步避免两个检测点x0，x1取中点，所以加了修正因子epsilon)

注意：笔者所用函数是整合了三种方法，并且可以处理各种异常情况（如区间不合法，单峰函数和多极值函数等），在下一章，笔者将分享整合的函数代码和分析讨论代码实现时的细节问题。