线搜索技术

线搜索技术

求解一维函数的极值点的数值迭代法称为一维搜索方法，又称为线搜索技术。实际工程优化问题的目标函数是一元函数的并不多，但是，在求解多元函数最优解的过程中，常常伴随着一系列的线搜索技术。线搜索方法对整个算法的收敛速度、精度都有较大的影响。可以说线搜索技术是优化方法的重要支柱。
在前面的优化技术理论基础中介绍了下降迭代算法的基本迭代公式
$$x_{k+1}=x_k+{\alpha }_k{d}_k,(k=0,1,2,\cdots )$$

当第k步迭代的初始点$x_k$及搜索方向$d_k$确定后，迭代所得到的的新点$x_{k+1}$取决于步长因子${\alpha }_k$，不同的${\alpha }_k$会得到不同的点$x_{k+1}$和不同的函数值$f(x_{k+1})$。因此，在多维优化的过程中一维线搜索的目的是在既定的$x_k$和$d_k$下，寻求最优步长因子${\alpha }_k$，使第k步迭代产生的新点$x_{k+1}$的函数值$f(x_{k+1})$为最小，即
$$minf(x_{k+1})=minf(x_k+\alpha d_k)=minf(\alpha )$$
显然，求解上式就是一维优化问题，如果一元函数是次数较低的简单函数，可以用一阶导数等于零来求解。若是高次复杂函数，其一阶导数的次数仍然比较高，难以求解，则需采用数值迭代的方法求近似解，这就是本文将要讨论的线搜索技术。
先搜索技术包括精确线搜索技术和非精确线搜索技术，线搜索技术的主要步骤分为两步：
（1）确定搜索区间（单峰区间），这是各种线搜索技术的共性问题；
（2）在搜索区间内求步长因子${\alpha }_k$，即选用具体的优化方法.

一、搜索区间的确定
（1）单峰区间
在进行一维搜索时，首先要在给定的方向上确定包含函数值最小值点的搜索区间，而且在该区间内，函数$f(\alpha )$有唯一最小值点${\alpha }^{\ast }$。这种只有一个函数峰值的区间称为单峰区间，如下图所示，

在极小点${\alpha }^{\ast }$的左边函数值是严格下降的，而在极小值点${\alpha }^{\ast }$的右边函数值是严格上升的，其函数图形呈“高-低-高”的形状。
（2）单峰区间的确定
在进行一维搜索时，必须首先找到单峰区间。若能实现估计出极小值点的大致位置，可以直接得到单峰区间；否则，需要用试探的方法寻找单峰区间，其常用的方法是进退法；
进退法的基本思想是：按照一定的规律给出三个试算点，依次比较各试算点函数值的大小，指导找出函数值按“大-小-大”变化的区间为止，其具体算法步骤如下

算法1：进退法
步1 给定初始点$a_1$和初始步长$h$;
步2 令$a_2=a_1+h$，得到两个试算点$a_1,a_2$，计算它们的函数值$f_1=f(a_1),f_2=f(a_2)$;
步3 比较$f_1$与$f_2$的大小，存在如下两种情况
（1）若$f_1>f_2$，则作前进运算，为了加速计算过程，步长增大一倍$h=2h$，取第三个试算点$a_3=a_2+h$,计算函数值$f_3=f(a_3)$，并比较$f_2$与$f_3$的大小；
若$f_3>f_2$,则找到了连续三个试算点$a_1,a_2,a_3$的函数值按照“大-小-大”变化，故搜索区间为$[a,b]=[a_1,a_3]$;
若$f_3<f_2$,则抛弃$a_1$，令$a_1=a_2,f_1=f_2,a_2=a_3,f_2=f_3$，返回到步3中的（1），重复该过程，直至找到连续三个试算点$a_1,a_2,a_3$，其函数值按“大-小-大”变化要求为止；
（2）若$f_1<f_2$，则作后退运算，令$h=-h$，将$a_1,f_1$与$a_2,f_2$对调，并取第三个试算点$a_3=a_2+h$，计算$f_3=f(a_3)$，比较对调后的函数值$f_2$与$f_3$:
若$f_3>f_2$，则搜索区间$[a,b]=[a_3,a_1]$;
若$f_3<f_2$，则步长加倍，抛弃点$a_1$，继续作后退运算，直至找到单峰区间为止。

二、区间消去原理
各种一维搜索方法的共同点是，在找到单峰区间后，再将单峰区间逐渐缩小，从而找到极小值点的近似解；
区间消去法的基本思想是：在单峰区间$[a,b]$内取两个点，通过比较该两点的函数值的大小，抛掉函数值较大的段，保留极小值点所在的段，反复比较舍弃，指导极小值点所在的单峰区间缩小到满足精度要求的长度，即$\overline{ab}\le ϵ$，这里$ϵ$为给定的精度，从而找到近似的极小值点，如下图所示

在单峰区间[a,b]内取两个点$x_1,x_2$，并计算对应的函数值$f_1,f_2$。若$f_1>f_2$，如上图a)所示，显然极小值点不在$a{x}_1$段内，取$a=x_1$，即消去了$a{x}_1$段，得到一个较小的单峰区间；同理如图b)所示若$f_1<f_2$，取$b=x_2$，即消去了$x_{2}b$段；如图c)所示，$f_1=f_2$，可取$a=x_1,b=x_2$，同时消去了$a{x}_1$段和$b=x_2$段；但往往为了简化程序，图c)所示的情况归并到前两种情况之一中去了。点$x_1,x_2$选取的方法不同，就得到不同的一维线搜索方法。

三、黄金分割法精确线搜索
0.618法又称黄金分割法，它是一种等比例缩短区间的搜索方法，其基本思想是：在单峰区间内取两队对称点，逐步把区间缩小，直到极小值点所在的区间缩小到给定的精度范围内，从而得到近似的最优解。
方法原理：设函数为$f(x)$，初始单峰区间为$[a,b]$，其长$L=b-a$，在$[a,b]$内取两对称点$x_1,x_2,a{x}_2=x_{1}b=\lambda L$，计算对应的函数值$f_1=f(x_1),f_2=f(x_2)$，比较它们的大小：
若$f_1>f_2$，显然，极小值点必在区间$[x_1,b]$内，因而抛掉区间$[a,x_1]$，即取$a=x_1$，则新取件为$[x_1,b]$，为了程序循环迭代，更新区间$a=x_1,x_1=x_2,f_1=f_2$，然后计算新点$x_2=a+\lambda (b-a)$;
令区间缩短率为
$$\frac{新区间长}{原区间长}=\frac{\lambda L}{L}=\lambda$$
在得到的新区间内继续比较$f_1,f_2$的大小，若$f_1>f_2$，抛掉区间$[a,x_1]$，即取$a=x_1$，则新区间为$[x_1,b]$，作相应的更新：$a=x_1,x_1=x_2,f_1=f_2$，然后计算新点$x_2=a+\lambda (b-a)$和函数值$f_2=f(x_2)$，则新区建缩短率为
$$\frac{新区间长}{原区间长}=\frac{(1-\lambda )L}{\lambda L}=\frac{1-\lambda }{\lambda }$$
黄金分割法的特点是取每次的缩短率相等，于是有
$$\lambda =\frac{1-\lambda }{\lambda }$$
即
$${\lambda }^2+\lambda -1=0$$
由上式求解可知$\lambda =\frac{\sqrt{5}-1}{2}\approx 0.618$，即每次缩小区间后，所得到的的新区间是原区间的0.618倍，舍弃的区间是原区间的0.382倍。在反复迭代计算过程中，除初始区间要找两个计算点外，在以后每次缩短的新区建内只需要再计算一个新点的函数值即可；对于上述过程，若$f_1<f_2$，应消去区间$[x_2,b]$新区间为$[a,x_2]$,作相应的更新$b=x_2,x_2=x_1,f_2=f_1$，然后计算新点$x_1=b-\lambda (b-a)$和函数值$f_1=f(x_1)$，比较$f_1,f_2$的大小，重复迭代

算法2 黄金分割法
步1 给定单峰区间$[a,b],\lambda =0.618$及收敛精度$ϵ$;
步2 在该区间内取两对称点并计算对应的函数值：
$$\begin{gathered} x_1=b-\lambda (b-a),f_1=f(x_1) \\ {x}_2=a+\lambda (b-a),f_2=f(x_2) \end{gathered}$$
步3 比较函数$f_1$与$f_2$的大小：
若$f_1>f_2$，则取区间[x_1,b]为新区间，而$x_2$则作为新区间内的第一个试算点，即令
$$a=x_1,x_1=x_2,f_1=f_2$$
而另一个试算点要重新计算
$$x_2=a+\lambda (b-a),f_2=f(x_2)$$
若$f_1\le f_2$，则取$[a,x_2]$为新区间，而$x_1$则作为新区间的第二试算点，即令
$$b=x_2,x_2=x_1,f_2=f_1$$
而另外一个试算点要重新计算
$$x_1=b-\lambda (b-a),f_1=f(x_1)$$

步4 若满足迭代精度$b-a\le ϵ$，则转到下一步，否则返回步骤3，进行下一次迭代计算，进一步缩短区间
步5 输出最优解$x^{\ast }=0.5(b+a),f^{\ast }=f(x^{\ast })$

黄金分割法的特点：
（1）每次区间的缩短率为0.618
（2）每一次只需要计算一个新点$x_1$或$x_2$的函数值，另外利用上次旧迭代点的信息，节省计算时间
（3）方法简单，对函数没有特殊要求，适应性强
（4）没有利用函数自身的数学特性，收敛速度慢

四、Armijo准则
线搜索技术是求解许多优化问题下降算法的基本组成部分，但精确线搜索往往需要计算很多的函数值甚至梯度值，对于复杂函数而言，这种函数值以及梯度值的计算耗费较多的计算资源。特别是当迭代点远离最优点时，精确线搜索通常不是十分有效与合理的。对于许多优化算法，其收敛速度并不依赖于精确线搜索过程。因此，既能保证目标函数具有可接受的下降量又能使最终形成的迭代序列收敛的非精确线搜索变得越来越流行，本小节介绍的Armijo准则以及Wolfe准则就是非精确线搜索常用的准则。

Armijo准则：设$d_k$是点$x_k$处的下降方向，若
$$f(x_k+\alpha d_k)\le f(x_k)+c_1\alpha \nabla f(x_k{)}^T{d}_k$$
则称步长$\alpha$满足Armijo准则，其中$c_1\in (0,1)$是一个常数.
Armijo准则有非常直观的几何含义，它指的是点$(\alpha ,\varphi (\alpha ))$必须在直线$l(\alpha )=\varphi (0)+c_1\alpha \nabla f(x_k{)}^T{d}_k$的下方

如上图所示，区间$[0,{\alpha }_1]$中的点均满足Armijo准则。在实际应用中，参数$c_1$通常选为一个很小的正数，这使得Armijo准则非常容易得到满足。但是仅仅使用Armijo准则并不能保证迭代的收敛性，这是因为$\alpha =0$也满足Armijo准则，但是这意味着迭代序列中的点固定不变，研究这样的步长没有意义。为此Armijo准则需要配合其他准则共同使用。在优化算法的视线中，寻找一个满足Armijo准则的步长是比较容易的，一个最常用的算法是回退法。给定初值$\hat{\alpha }$，回退法通过不断以指数方式缩小试探步长，找到第一个满足Armijo准则的点。具体来说，回退法选取${\alpha }_k={\gamma }^{j_0}\hat{\alpha }$，其中
j 0 = m i n { j = 0 , 1 , ⋯ ∣ f ( x k + γ j α ^ d k ) ≤ f ( x k ) + c 1 γ j α ^ ∇ f ( x k ) T d k } j_0=min \{j=0,1,\cdots|f(x_k+\gamma^j\hat{\alpha}d_k)\le f(x_k)+c_1\gamma^j\hat{\alpha}\nabla f(x_k)^Td_k\} j0​=min{j=0,1,⋯∣f(xk​+γjα^dk​)≤f(xk​)+c1​γjα^∇f(xk​)Tdk​}
参数$\gamma \in (0,1)$为一给定的实数。回退法的基本过程如下

armijo回退法：
步1 选择初始步长$\hat{\alpha }$，参数$\gamma \in (0,1)$，初始化$\alpha =\hat{\alpha }$
步2 while $f(x_k+\alpha d_k)>f(x_k)+c\alpha \nabla f(x_k{)}^T{d}_k$ do
令$\alpha =\gamma \alpha ,k=k+1$
end while
步3 输出${\alpha }_k=\alpha$

该算法被称为回退法是因为$\alpha$的试验值是由大到小的，它可以确保输出的${\alpha }_k$能尽量地大。此外回退法不会无限进行下去，因为$d_k$是一个下降方向，当$\alpha$充分小时，Armijo准则总是成立的，在实际应用中我们通常会给$\alpha$设置一个下界，防止步长过小。

五、Wolfe准则
Wolfe准则：设$d_k$是点$x_k$处的下降方向，若
$$\begin{array}{rl}& f(x_k+\alpha d_k)\le f(x_k)+c_1\alpha \nabla f(x_k{)}^T{d}_k,\\ & \nabla f(x_k+\alpha d_k{)}^T{d}_k\ge c_2\nabla f(x_k{)}^T{d}_k,\end{array}$$
则称步长$\alpha$满足Wolfe准则，其中$c_1,c_2\in (0,1)$为给定的常数且$c_1<c_2$.
在该准则中，第一个不等式记为Armijo准则，而第二个不等式则是Wolfe准则的本质要求，Wolfe准则要求$\varphi (\alpha )$在点$\alpha$处的斜率不能小于其在零点的斜率的$c_2$倍。注意到在$\varphi (\alpha )$的极小值点$\alpha \ast$处有${\varphi }^{{}^{\prime }}({\alpha }^{\ast })=\nabla f(x_k+{\alpha }^{\ast }{d}_k{)}^T{d}_k=0$，因此$\alpha \ast$永远满足Wolfe准则的第二式，而选择较小的$c_1$可使${\alpha }^{\ast }$同时满足第一式，即Wolfe准则在绝大多数情况下会包含线搜索子问题的精确解。在实际应用中，参数$c_2$通常选取为0.9

六、线搜索收敛性
线搜索算法的框架：
步0 初始化，选取有关参数的初始迭代点$x_0$，设定容许误差$ϵ$，迭代步数$k=0$；
步1 检验终止判别主责，计算$g_k=\nabla f(x_k)$，若$|||g_k|\le ϵ$，输出$x^{\ast }\approx x_k$，停止计算；
步2 确定下降方向$d_k$，使满足$g_k^T{d}_k<0$;
步3 确定步长因子${\alpha }_k$，可以在系列“精确”与“非精确”两种线搜索技术中选用其一：
（1）用前面介绍的黄金分割法或者其他精确线搜索技术求
$${\alpha }_k=arg\underset{\alpha >0}{\min }f(x_k+\alpha d_k)$$
（2）用前面介绍的Armijo准则或者Wolfe准则等非精确线搜索技术求${\alpha }_k$;
步4 迭代更新点，令$x_{k+1}=x_k+{\alpha }_k{d}_k,k=k+1$，转步1

值得说明的是，为了保证这类算法的收敛性，除了在步3要选用适当的线搜索技术之外，搜索方向$d_k$也需要满足一定的条件，即对于所有的$k,d_k$与$-g_k$的夹角满足
$$0\le {\theta }_k\le \frac{\pi }{2}-\mu ,\mu \in (0,\frac{\pi }{2})$$即要求$d_k$与负梯度方向夹角为锐角

七、线搜索代码实现及举例
（1）进退法
进退法代码实现如下

function [a, b] = forward_backward_method(func, x0, h)

a1 = x0;
a2 = a1 + h;
f1 = feval(func, a1);
f2 = feval(func, a2);

maxIter = 1000;
iIter = 0;

while iIter < maxIter
    if f1 < f2
        h = -h;
        a3 = a1; f3 = f1;
        a1 = a2; f1 = f2;
        a2 = a3; f2 = f3;
    end
    h = 2 * h; a3 = a2 + h; f3 = feval(func, a3);
    
    if f3 > f2
        if h > 0
            a = a1; b = a3;
            break;
        else
            a = a3; b = a1;
            break
        end
    else
        a1 = a2; f1 = f2;
        a2 = a3; f2 = f3;
    end
    
    iIter = iIter + 1;
end

if iIter == 1000
    fprintf('reach maximum iteration, and not found suitable interval!\n');
    a = 0;  b = 0;
else
    fprintf('find interval: [a, b] = [%d, %d]\n', a, b);
end

end

进退法实验1
$f(x)=x^2-\sin (x)$，初始点及步长为$x0=-1,h=0.2$

>> x0 = -1;
>> h = 0.2;
>> func = inline('x.^2 - sin(x)')
>> [a, b] = forward_backward_method(func, x0, h);
find interval: [a, b] = [-4.000000e-01, 2]

进退法实验2
$f(x)=\frac{1}{x}$

> x0 = 1;
>> h = 0.01;
>> [a, b] = forward_backward_method(func, x0, h);
reach maximum iteration, and not found suitable interval!

（2）黄金分割法
黄金分割法线搜索代码实现如下：

function [xmin, fmin] = golden_line_search(func, a, b, epsilon)

lambda = ( sqrt(5) - 1 ) / 2;
L = b - a;
x1 = b - lambda * L; f1 = feval(func, x1);
x2 = a + lambda * L; f2 = feval(func, x2);

iIter = 0;
iterMax = 1000;

while iIter < iterMax
    iIter = iIter + 1;
    if f1 > f2
        a = x1; x1 = x2; f1 = f2;
        x2 = a + lambda * (b - a); f2 = feval(func, x2);
    else
        b = x2; x2 = x1; f2 = f1;
        x1 = b - lambda * (b - a); f1 = feval(func, x1);
    end
    
    if abs(b - a) <= epsilon
        break;
    end
    
end

if iIter == iterMax
    fprintf('reach maimum iteration, and not found xmin!\n');
    xmin = 0; fmin = 0;
else
    xmin = 0.5 * (a + b); fmin = feval(func, xmin);
end

end

黄金分割法实验1：
求$f(x)=x+\frac{20}{x}$在$[0.2,1]$上的最小值

function f = mostCheapVelocity(x)

f = x + 20 / x;

end

% golden line search test
x1 = 0.2;
x2 = 1;
epsilon = 0.0001;
[fmin, xmin] = golden_line_search('mostCheapVelocity', x1, x2, epsilon);
fprintf('fmin = %f, xmin = %f\n', fmin, xmin);


[x, fval] = fminbnd('mostCheapVelocity', x1, x2);
fprintf('matlab buildin 1d search: fmin = %f, xmin = %f\n', fval, x);

[fmin, xmin] = squareInterpSearch('mostCheapVelocity', x1, x2, epsilon);
fprintf('square interp search: fmin = %f, xmin = %f\n', fmin, xmin);

计算结果：

fmin = 0.996748, xmin = 21.062005
matlab buildin 1d search: fmin = 21.000992, xmin = 0.999948
square interp search: fmin = 23.193544, xmin = 0.897000

从实验结果来看，实现代码计算的结果与matlab内置的线搜索的结果一致

（3）Armijo准则
Armijo准则代码实现：

function [fnew, xnew, alpha] = armijo_rule(func, gfunc, x0, alpha0, dk)

c1 = 1e-3;
alpha = alpha0;
gamma = 0.8;

iIter = 0;
iterMax = 200;
alphaMin = 1e-5;

while iIter < iterMax 
    
    xnew = x0 + alpha * dk;
    fnew = feval(func, xnew);
    f0 = feval(func, x0);
    if fnew <= f0 + c1 * feval(gfunc, x0)' * dk
        break;
    end
    
    if alpha < alphaMin
        break;
    end
    
    alpha = gamma * alpha;
    iIter = iIter + 1;   
end

if iIter == iterMax
    alpha = alphaMin;
    fprintf('reach maximum iteration, and not found suitable alpha.\n');
end

xnew = x0 + alpha * dk;
fnew = feval(func, xnew);

end

（4）Wolfe准则
Wolfe准则代码实现：

function [fnew, xnew, alpha] = armijo_rule(func, gfunc, x0, alpha0, dk)

c1 = 1e-3;
alpha = alpha0;
gamma = 0.8;

iIter = 0;
iterMax = 200;
alphaMin = 1e-5;

while iIter < iterMax 
    
    xnew = x0 + alpha * dk;
    fnew = feval(func, xnew);
    f0 = feval(func, x0);
    if fnew <= f0 + c1 * feval(gfunc, x0)' * dk
        break;
    end
    
    if alpha < alphaMin
        break;
    end
    
    alpha = gamma * alpha;
    iIter = iIter + 1;   
end

if iIter == iterMax
    alpha = alphaMin;
    fprintf('reach maximum iteration, and not found suitable alpha.\n');
end

xnew = x0 + alpha * dk;
fnew = feval(func, xnew);

end

考虑无约束优化问题
$$\underset{x\in R^2}{\min }f(x)=100(x_1^2-x_2{)}^2+(x_1-1{)}^2$$设当前迭代点$x_k=(-1,1{)}^T$，下降方向$d_k=(1,-2{)}^T$，利用Armijo和Wolfe准则求解步长因子${\alpha }_k$
函数计算

function f = func(x)
f = 100 * (x(1)^2 - x(2))^2 + (x(1) - 1)^2;
end

梯度计算

function f = func(x)
f = 100 * (x(1)^2 - x(2))^2 + (x(1) - 1)^2;
end

调用armijo准则以及wolfe准则计算下一步迭代步长：

% armijo rule test
xk = [-1, 1]';
dk = [1, -2]';
alpha0 = 10;


[fNew, xNew, alpha] = armijo_rule('func', 'gfunc', xk, alpha0, dk);
fOld = feval('func', xk);
fprintf('armijo rule: fOld = %f, fNew = %f, xNew = [%f, %f], alpha = %f\n', fOld, fNew, xNew(1), xNew(2), alpha);

[fNew, xNew, alpha] = wolfe_rule('func', 'gfunc', xk, alpha0, dk);
fOld = feval('func', xk);
fprintf('wolfe rule: fOld = %f, fNew = %f, xNew = [%f, %f], alpha = %f\n', fOld, fNew, xNew(1), xNew(2), alpha);

结果为：

armijo rule: fOld = 4.000000, fNew = 3.581038, xNew = [-0.718525, 0.437050], alpha = 0.281475
wolfe rule: fOld = 4.000000, fNew = 3.581038, xNew = [-0.718525, 0.437050], alpha = 0.281475