【最优化】精确一维线搜索—黄金分割法在Matlab上的实现

1 前言

  在线搜索迭代算法中，当迭代方向$d^k$确定后，需要看沿$d^k$ 走多远对算法有利，即需要寻找迭代步长。迭代步长通常由精确一维线搜索、非精确一维线搜索、非单调一维线搜索等方式获得。实现精确一维线搜索的方法有很多，其中包括$\text{Fibonacci}$ 法、黄金分割法、进退法、平分法、抛物线法等。下面只介绍黄金分割法。黄金分割法又称为$0.618$法。

2 算法内容

  算法（黄金分割法） 给定目标函数$f(x),$初始搜索区间$[a_0,b_0]\subseteq [0,+\infty )$和精度要求$\epsilon >0.$ 若$b_0-a_0<\epsilon ,$那么算法终止$,$输出${\lambda }_{\ast }:=\frac{a_0+b_0}{2}.$

  步1 置$t_0=a_0+0.382(b_0-a_0),s_0=a_0+0.618(b_0-a_0),$计算$f(t_0),f(s_0).$置$k:=0$.

  步2 若$f(t_k)>f(s_k),$则转步3$;$否则转步4.

  步3 若$b_k-t_k<\epsilon ,$则算法终止$,$输出${\lambda }_{\ast }:=s_k;$否则置$a_{k+1}:=t_k,b_{k+1}:=b_k,t_{k+1}$$:=s_k,f(t_{k+1}):=f(s_k),$计算$s_{k+1}:=a_{k+1}+0.618(b_{k+1}-a_{k+1})$和$f(s_{k+1}),$转步5.

  步4 若$s_k-a_k<\epsilon ,$则算法终止$,$输出${\lambda }_{\ast }:=t_k;$否则置$a_{k+1}:=a_k,b_{k+1}:=s_k,s_{k+1}$$:=t_k,f(s_{k+1}):=f(t_k),$计算$t_{k+1}:=a_{k+1}+0.382(b_{k+1}-a_{k+1})$和$f(t_{k+1}),$转步5.

  步5 置$k:=k+1,$转步2.

3 在Matlab上的实现

3.1 代码

3.1.1 Goldsection.m（函数.m文件）

function [target_value,target_x,k,lambda] = Goldsection(f,epsilon,x,d,upper,lower)
% 黄金分割法
k=0;
kk=k+1;
if (upper-lower)<epsilon
    lambda=(lower(kk)+upper(kk))/2;
else
    t=lower+0.382*(upper-lower);   s=lower+0.618*(upper-lower);
    theta_t=f(x+t*d);      theta_s=f(x+s*d);
    while 1==1
        if theta_t(kk)>theta_s(kk)
            if (upper(kk)-t(kk))<epsilon
                lambda=s(kk);
                break
            else
                lower(kk+1)=t(kk);   upper(kk+1)=upper(kk);
                t(kk+1)=s(kk);
                s(kk+1)=lower(kk+1)+0.618*(upper(kk+1)-lower(kk+1));
                theta_t(kk+1)=f(x+t(kk+1)*d);
                theta_s(kk+1)=f(x+s(kk+1)*d);
                k=k+1;
                kk=k+1;
            end
        else
            if (s(kk)-lower(kk))<epsilon
                lambda=t(kk);
                break
            else
                lower(kk+1)=lower(kk);  upper(kk+1)=s(kk);
                s(kk+1)=t(kk);
                t(kk+1)=lower(kk+1)+0.382*(upper(kk+1)-lower(kk+1));
                theta_t(kk+1)=f(x+t(kk+1)*d);
                theta_s(kk+1)=f(x+s(kk+1)*d);
                k=k+1;
                kk=k+1;
            end
        end
    end
end
target_x=x+lambda*d;
target_value=f(x+lambda*d);
fprintf('算法迭代次数为%d，输出的λ=%f，目标函数值为：%f\n',k,lambda,target_value);  %输出结果
target_x
end

main.m

clear all
clc

%% 数据输入
% 创建目标函数的句柄函数
f=@(x)(1-x(1))^2+100*(x(2)-x(1)^2)^2;
%输入起始点
x=[0,0]';
%输入迭代方向
d=[1,0]';
%输入精度
epsilon=0.001;
%初始区间下界
lower=0;
%初始区间上界
upper=1;    

%% 通过黄金分割法计算
[target_value,target_x,k,lambda] = Goldsection(f,epsilon,x,d,upper,lower);

%target_value 储存当前精度下目标函数极小值
%target_x 储存当前精度下目标函数为极小值时x的位置，通常是向量
%  k    储存迭代次数
%lambda 储存迭代步长

3.2 使用说明

以下题为例：
  求Rosenbrock函数
$$\min f(x)=(1-x_1{)}^2+100(x_2-x_1^2{)}^2$$
在点$x_k=(0,0{)}^T$处沿方向$d^k=(1,0{)}^T$的近似步长。其中黄金分割法的初始区间$[0,1]$，精度$ϵ=0.001$.

1. 保证main.m和Goldsection.m文件在同一目录下.
2. 打开main.m，
   输入目标函数句柄函数：f=@(x)(1-x(1))$^2+100*(x(2)-x(1)^2)^2;
   修改为起始点x=[0,0]';
   修改迭代方向d=[1,0]';
   修改精度epsilon=0.001;
   修改初始区间下界lower=0;
   修改初始区间上界upper=1;
3. 运行，即可得到
   
   其中target_x为当前精度下目标函数为极小值时x的位置。