无约束优化方法之powell方法

无约束优化方法是优化技术中心极为重要也是最基本的内容之一。无约束优化问题的数学模型为：
$$\underset{x\in R^n}{\min }f(x)$$
无约束优化方法不仅可解决设计中无约束优化问题，更重要的是通过对无约束优化方法的研究，给求解约束优化方法提供良好的概念和理论基础，而且约束优化问题通过数学处理可化为无约束优化问题，利用无约束优化方法来求解。所以无约束优化方法在优化设计中有十分重要的作用；
根据确定搜索方向所使用的信息和方法，无约束优化方法可以分为两大类：
（1）直接法：这列方法仅利用计算目标函数值的信息来构造搜索方向。如坐标轮换法、powell法和单纯形法等。由于只需要计算和比较目标函数值，对于那些无法求导数或求导数很困难的函数，这类方法就有适用性广和可靠性高的显著优越性，但是这类方法一般收敛速度比较慢。
（2）间接法：这类方法利用函数的一阶导数甚至二阶导数的信息来构造搜索方向。如梯度法、牛顿法、变尺度法等。由于需要计算函数的一阶导数甚至二阶导数，因此这类方法计算量较大，但收敛速度较快。

本篇重点讨论共轭方向及powell法。

一、共轭方向的定义
在优化方法中，共轭方向的概念有着重要的意义，一些有效的无约束方法大多都是以共轭方向为搜索方向而形成的。首先考察二维的二次函数
$$f(x)=\frac{1}{2}{x}^{T}Ax+b^{T}x+c$$

如上图，设由任一点$x_0$出发，沿下降方向$s_0$进行一维搜索得到$x_1$，且
$$x_1=x_0+{\alpha }_0{s}_0$$
再由点$x_1$出发，沿直指极值点$x^{\ast }$的方向$s_1$进行一维搜索得极值点$x^{\ast }$.
显然
$$x^{\ast }=x_1+{\alpha }_1{s}_1$$
下面分析直指极值点$x^{\ast }$的方向$s_1$与方向$s_0$的关系。
对二次型求导可知
$$\nabla f(x)=Ax+b$$
由极值的必要条件可知$\nabla f(x^{\ast })=A{x}^{\ast }+b=0$
于是有
$$\nabla f(x^{\ast })=A(x_1+{\alpha }_1{s}_1)+b=(A{x}_1+b)+{\alpha }_{1}A{s}_1=0$$
改写上式有
$$\nabla f(x^{\ast })=\nabla f(x_1)+{\alpha }_{1}A{s}_1=0$$
两边左乘$s_0$，有
$$s_0^T\nabla f(x_1)+\alpha s_0^{T}A{s}_1=0$$
由梯度的性质可知：$\nabla f(x_1)$与$s_0$正交，又由于最优步长因此不等于零，故有
$$s_0^{T}A{s}_1=0$$
这就是两个向量共轭的数学表达式；
定义：设A是nxn阶实对称正定矩阵，如果有两个n为非零向量$s_1$与$s_2$，若满足：
$$s_1^{T}A{s}_2=0$$
则称$s_1$与$s_2$关于对称正定矩阵A是共轭的，或简称向量$s_1$与$s_2$关于A共轭；
若非零向量组$s_1,s_2,\cdots ,s_k$满足
$$s_i^{T}A{s}_j=0,(i\ne j;i,j=1,2,\cdots ,k)$$
则称向量组$s_1,s_2,\cdots ,s_k$关于矩阵A是共轭的。A的一个特例是目标函数的海塞矩阵；若A是单位矩阵I，则向量组互相正交；若A不是单位矩阵，则可理解为向量$s_2$经过矩阵A线性变换后的向量$(s_2{)}^{{}^{\prime }}=A{s}_2$与向量$s_1$正交。这说明政教向量是共轭的一种特例，或者说向量共轭是正交的推广。

二、共轭方向的性质
这里不加证明的给出共轭方向的性质；
可以证明，对于一般函数，共轭方向具有以下重要性质，这里设A为nxn阶实对称正定矩阵：
（1）若$s_1,s_2,\cdots ,s_n$为关于矩阵A共轭的n哥非零向量，则这一组向量线性无关。
由此可得推论：在n维空间中相互共轭的非零向量的个数不超过n个
（2）设向量$s_i^1(i=1,2,\cdots ,n)$是一组线性无关的非零向量，则可构造出n个非零向量$s_i^2(i=1,2,\cdots ,n)$，使满足
$$(s_i^2{)}^{T}A{S}_j^1=0,(i\ne j;i,j=1,2,\cdots ,n)$$
这说明共轭向量由线性无关的向量构造。由于坐标轴方向是线性无关的，故可得推论：可有坐标轴方向来构造共轭向量。
（3）若$s_1,s_2,\cdots ,s_n$为关于矩阵A共轭的n哥非零向量，则对于正定的n维二次函数
$$f(x)=\frac{1}{2}{x}^{T}Ax+b^{T}x+c$$
从任一点$x_0$出发，依次沿$s_i,(i=1,2,\cdots ,n)$方向进行一维搜索，之多n次变可收敛到极小值点$x^{\ast }$。这说明沿共轭方向搜索具有二次收敛速度。而对于非二次函数上述的极小值点还不是该函数的极小值点，需要进一步搜索，得到新的共轭方向组，反复迭代，直至最后逼近最小值点。

n维函数共轭方向的形成：沿n个线性无关的方向进行一轮一维搜索，连接该轮迭代的起始点$x_0^k$和终止点$x_n^k$的向量，即为新产生的共轭方向。

三、共轭方向法
（1）基本思想
根据上面的共轭方向的概念和性质，共轭方向法的基本思想：首先采用坐标轴方向作为第一轮搜索方向进行一轮一维搜索迭代。然后将第一轮迭代的初始点与最末一个极小值点连接构成新的方向，以该新的方向作为最末一个方向，淘汰第一个方向，得下一轮迭代的n个搜索方向。依此下去，直至逼近极值点。
（2）共轭方向法的步骤

---

步1 给定初始点$x_0$，迭代轮数$k=1$，令$x_0^k=x_0$，收敛精度$ϵ$
步2 输入坐标轴方向作为第一轮的搜索方向，即
$$s_i^k=e_i,e_i=[0,\cdots ,1,0,\cdots ,0{]}^T,(i=1,2,\cdots ,n)$$
步3 从$x_0^k$开始，以此沿方向$s_i^k$进行n次一维搜索得到n个一维极小值点$x_i^k,(i=1,2,\cdots ,n)$，连接点$x_0^k$与$x_n^k$，构造第一个共轭方向：
$$s_{n+1}^k=x_n^k-x_0^k$$
步4 从$x_n^k$出发沿$s_{n+1}^k$进行一维搜索，得极小值点$x_{n+1}^k$(完成了一轮迭代)
步5 收敛精度判断，若$||x_{n+1}^k-x_0^k||<ϵ$，则输出$x^{\ast }=x_{n+1}^k$，结束；否则，淘汰$s_1^k$，增加$s_{n+1}^k$，构成第k+1轮的搜索方向：$s_i^{k+1}=s_{i+1}^k,(i=1,2,\cdots ,n)$
步6 将$x_{k+1}^k\to x_0^{k+1}$，令k=k+1，返回步3

---

（3）共轭方向法存在的问题
在上述共轭方向基本算法中，方向组的替换采用了固定的格式，运算比较简便。但由此形成的方向组中，有可能出现几个方向线性相关或近似线性相关的现象。也就是说，在新的方向组中有可能存在两个方向平行的情况。显然，这将导致新的方向组共轭方向个数减小，从而使迭代运算退化到一个较低维空间中进行，即降维搜索，因此无法得到真正的极小值点。
以三维空间为例：从$x_0^1$出发分别沿三个方向作一维搜索，得到三个方向上的极小值点为：
$$\begin{array}{rl}{x}_1^1& =x_0^1+{\alpha }_1^1{s}_1^1\\ x_2^1& =x_1^1+{\alpha }_2^1{s}_2^1=x_0^1+{\alpha }_1^1{s}_1^1+{\alpha }_2^1{s}_2^1\\ x_3^1& =x_2^1+{\alpha }_3^1{s}_3^1=x_0^1+{\alpha }_1^1{s}_1^1+{\alpha }_2^1{s}_2^1+{\alpha }_3^1{s}_3^1\end{array}$$
由首末两点构造的新方向为：
$$s_4^1=x_3^1-x_0^1={\alpha }_1^1{s}_1^1+{\alpha }_2^1{s}_2^1+{\alpha }_3^1{s}_3^1$$
若其中的${\alpha }_1^1=0$或者${\alpha }_1^1\approx 0$，也就是说，在方向$s_1^1$上函数无下降或者下降得很小，则上式为
$$s_4^1={\alpha }_2^1{s}_2^1+{\alpha }_3^1{s}_3^1$$
若以$s_4^1$替换$s_1^1$，则由上式看出，由此形成的三个方向线性相关，因此，以后的迭代将始终局限在该平面内，使三维的问题退化到二维问题。显然在二维平面内无法搜索到三维空间中的真正极小值点。鉴于以上问题，M.J.D.Powell与1964年对共轭方向基本算法进行了改进，形成了Powell法

四、共轭方向阀的改进—Powell法
（1）基本思想
powell法与上述为改进的共轭方向法迭代过程基本相同，其不同的点是：
a）在（k+1）轮迭代中，是否选用新产生的方向$s_{n+1}^k$要进行判断；
b）若选用，也不一定淘汰第k轮的第一个方向$s_1^k$，淘汰哪一个方向也要进行判断，其目的是避免新一轮的搜索方向线性相关，出现退化现象，导致搜索降维。
（2）共轭方向$s_{k+1}^n$的选取与淘汰方向的判断
在得到新的方向$s_{n+1}^k=x_n^k-x_0^k$之后，沿方向$s_{n+1}^k$找到点$x_0^k$关于点$x_n^k$的反射点$x_{n+2}^k$，如下图所示

$$x_{n+1}^k=x_n^k+(x_n^k-x_0^k)=2x_n^k-x_0^k$$
分别计算三个点的函数值
$$f_1=f(x_0^k),f_2=f(x_n^k),f_3=f(x_{n+2}^k)$$
然后找出前一轮一维搜索过程中函数值下降最大的方向$s_m^k$及其最大的下降量${\Delta }_m$，即${\Delta }_m^k={\max }_{i=1,2,\cdots ,n}|f(x_{i-1}^k)-f(x_i^k)|$，Powell的研究结果，若以下关系
$$\begin{array}{rl}& f_3<f_1\\ & (f_1+f_3-2f_2)(f_1-f_3-{\Delta }_m^k{)}^2<\frac{1}{2}{\Delta }_m^k(f_1-f_3{)}^2\end{array}$$
同时满足，则表明方向$s_{n+1}^k$与园方向线性无关，因此可以用来替换。替换的方向就是${\Delta }_m$所对应的方向$s_m^k$，即淘汰原方向$s_m^k$，而把新的方向$s_{n+1}^k$加到方向组的末尾。替换的通式为
$$\begin{array}{r}{s}_i^{k+1}=s_i^k,i<m\\ s_i^{k+1}=s_{i+1}^k,i>m\end{array}$$
若Powell判别式不满足，则表明$s_{n+1}^k$与原方向组中的某些方向线性相关，因此不能选用新方向，扔用第k轮中的迭代方向作为第k+1轮的搜索方向，即令
$$s_i^{k+1}=s_i^k,(i=1,2,\cdots ,n)$$

Powell法的具体步骤

---

步1 给定初始点$x_0$，迭代论数k=1，令$x_0^k=x_0$，维数n，收敛精度${ϵ}_1,{ϵ}_2$
步2 输入坐标轴方向作为第一轮的搜索方向，即
$$s_i^k=e_i,e_i=[0,\cdots ,1,0,\cdots ,0{]}^T,(i=1,2,\cdots ,n)$$
步3 从$x_0^k$开始，依次沿方向$s_i^k$进行一轮（n次）一维搜索得到n个一维极小值点：
$$x_i^k=x_{i-1}^k+{\alpha }_i^k{s}_i^k,(i=1,2,\cdots ,n)$$
计算：${\Delta }_m^k={\max }_{i=1,2,\cdots ,n}|f(x_{i-1}^k)-f(x_i^k)|$，得相应的方向$s_m^k$;
步4 收敛精度判断：若满足$||x_n^k-x_0^k||\le ϵ1$，并且$|\frac{f(x_0^k)-f(x_n^k)}{f(x_0^k)}|\le {ϵ}_2$，则令$x^{\ast }=x_n^k,f(x^{\ast })=f(x_n^k)$，结束迭代；否则转下一步；
步5 构造共轭方向，连接点$x_0^k$与$x_n^k$构造共轭方向$s_{n+1}^k=x_n^k-x_0^k$
步6 求反射点 $x_{n+2}^k=2x_n^k-x_0^k$，计算函数值$f_1=f(x_0^k),f_2=f(x_n^k),f_3=f(x_{n+2}^k)$;
步7 Powell条件判别：若满足判别式
$$\begin{array}{rl}& f_3<f_1\\ & (f_1+f_3-2f_2)(f_1-f_3-{\Delta }_m^k{)}^2<\frac{1}{2}{\Delta }_m^k(f_1-f_3{)}^2\end{array}$$
则取共轭方向$s_{n+1}^k$，转下一步；步满足，则转步10
步8 从点$x_n^k$出发，沿方向$s_{n+1}^k$作一维搜索取得极小值点$x_{n+1}^k=x_n^k+{\alpha }_n^k{s}_{n+1}^k$
步9 进行方向替换，构成新的方向组$s_i^{k+1},(i=1,2,\cdots ,n)$，令
$$\begin{array}{r}{s}_i^{k+1}=s_i^k,i<m\\ s_i^{k+1}=s_{i+1}^k,i>m\end{array}$$
$x_0^{k+1}=x_{n+1}^k,k=k+1$，转步3；
步10 不进行方向替换，扔用第k轮的搜索方向组，令$s_i^{k+1}=s_i^k,(i=1,2,\cdots ,n)$，若$f_2<f_2$时，取$x_0^{k+1}=x_n^k$；否则，取$x_0^{k+1}=x_{n+2}^k,k=k+1$，转步3；

---

五、powell法的代码实现

function [fmin, xmin] = powell_method(func, func1D, x0, eps1, eps2, N)

%k = 1;
%x0_old = x0;
s = diag(ones(1, N));
x = zeros(N, N+1);
x(:, 1) = x0;

iIter = 0;
iterMax = 5000;
epsilon = 1e-5;
deltaMax = zeros(N, 1);

while iIter < iterMax
    for i = 2:N+1
        h = 0.5;
        alpha0 = 0;
        func1DTemp = @(y)(func1D(y, s(:, i-1)));
        [alphaMin, alphaMax] = forward_backward_method(func1DTemp, alpha0, h);
        [alpha, ~] = golden_line_search(func, alphaMin, alphaMax, epsilon);
        x(:, i) = x(:, i-1) + alpha * s(:, i-1);
    end
    
    % find max
    for i = 2:N+1
        deltaMax(i-1, 1) = abs( feval(func, x(:, i-1)) - feval(func, x(:, i)) );
    end
    
    [deltaMaxK, kMax] = max(dletaMax);
    
    if norm(x(:,N+1) - x(:, 1), 2) <= eps1
        fxN = feval(func, x(:, N+1));
        fx0 = feval(func, x(:, 1));
        if abs((fxN - fx0) / fx0) <= eps2
            xmin = x(:, N+1);
            fmin = feval(func, xmin);
            break;
        end
    end
    
    sNP1 =  x(:, N+1) - x(:, 1);
    xNP2 =  2 * x(:, N+1) - x(:, 1);
    f1 = feval(func, x(:, 1));
    f2 = feval(func, x(:, N+1));
    f3 = feval(func, xNP2);
    y1 =  (f1 + f3 - 2 * f2) * (f1 - f2 - deltaMaxK)^2;
    y2 = 0.5 * deltaMaxK * (f1 - f3)^3;
    
    if (f3 < f1) && (y1 < y2)
        h = 0.5;
        alpha0 = 0;
        func1DTemp = @(y)(func1D(y, sNP1));
        [alphaMin, alphaMax] = forward_backward_method(func1DTemp, alpha0, h);
        [alpha, ~] = golden_line_search(func, alphaMin, alphaMax, epsilon);
        xNP1 = x(:, N+1) + alpha * sNP1;
        
        for i = kMax:N-1
            s(:, i) = s(:, i+1);
        end
        
        s(:, N) = sNP1;
        x(:, 1) = xNP1; 
    else
        if f2 < f3
            x(:, 1) = x(:, N+1);
        else
            x(:, 1) = xNP2;
        end
    end
    
    iIter = iIter + 1;
    %k = k + 1;
end

if iIter == iterMax
    fprintf('reach maximum iteration, and not found suitable xmin.\n');
    xmin = zeros(N, 1);
    fmin = feval(func, xmin);
end  
    
end