最优化理论-拟牛顿法[Matlab]

努力上进的math boy

已于 2024-07-18 17:12:37 修改

阅读量375

点赞数 6

分类专栏：运筹学文章标签： matlab 开发语言算法

于 2024-07-17 02:05:32 首次发布

本文链接：https://blog.csdn.net/qq_52176723/article/details/140472600

版权

运筹学专栏收录该内容

4 篇文章 0 订阅

订阅专栏

1.牛顿法的局限性

牛顿法是一种具有较高实用性的优化问题求解方法。如果牛顿法收敛，其收敛阶数至少是2。但牛顿法有其局限性，这里列举两个有代表性的缺陷，其中一个缺陷是引入拟牛顿法的基本思路。

牛顿法的基本思路是在每次迭代中，利用二次型函数来局部近似目标函数 $f$ ，并求解近似函数的极小点作为下一个迭代点，迭代公式为：

$x^{(k+1)}=x^{(k)}-F(x^{(k)})^{-1}g^{(k)}$

该公式在初始点 $x^{(0)}$ 不足够接近极小点时，牛顿法可能不具有下降特性[即对于某些k，有 $f(x^{(k+1)})\nless f(x^{(k)})$ ]

这便是牛顿法的第一个缺陷，对于该问题，可以引入一个适当的步长 $\alpha_{k}$ ，使得牛顿法每一步迭代具有下降特性，即

$x^{(k+1)}=x^{(k)}-\alpha_{k}F(x^{(k)})^{-1}g^{(k)}$

其中， $\alpha_{k}=argmin_{\alpha\geqslant 0}f(x^{(k)}-\alpha_{}F(x^{(k)})^{-1}g^{(k)})$ ，可以使用精确线搜索或非精确线搜索（Armijo条件，Wolfe条件等），合理地确定步长，使得

$f(x^{(k+1)})< f(x^{(k)})$

但牛顿法还有另外一个比较严重的缺陷，就是必须要计算第k步函数的二阶黑塞矩阵 $F(x^{(k)})^{-1}$ ，这十分消耗计算资源，为了避免这种求逆运算，可以设计 $F(x^{(k)})^{-1}$ 的近似矩阵来替代 $F(x^{(k)})^{-1}$ ，这便是拟牛顿法的基本思路。 $F(x^{(k)})^{-1}$ 的近似矩阵随着迭代不断更新，使其至少拥有一部分 $F(x^{(k)})^{-1}$ 的部分性质。即令 $H_{k}=F(x^{(k)})^{-1}$ ，得到

$x^{(k+1)}=x^{(k)}-\alpha_{k}H_{k}g^{(k)}$

其中的 $H_k$ 只需要利用目标函数值和梯度，因此，只要确定了合适的近似矩阵 $H_k$ 的构造方法，在迭代过程中不需要任何涉及到黑塞矩阵以及求逆的计算工作。

2.黑塞矩阵逆矩阵的近似

令

$x^{(k+1)}=x^{(k)}-\alpha_{}H_{k}g^{(k)}$

在 $x^{(k)}$ 处对 $f$ 进行一阶泰勒展开，并将 $x^{(k+1)}$ 代入，可得

$f(x^{(k+1)})=f(x^{(k)})+g^{(k)T}(x^{(k+1)}-x^{(k)})+o(||x^{(k+1)}-x^{(k)}||)$

代入 $x^{(k+1)}=x^{(k)}-\alpha_{}H_{k}g^{(k)}$ 得

$f(x^{(k+1)})=f(x^{(k)})-\alpha_{}g^{(k)T}H_{k}g^{(k)}+o(||H_{k}g^{(k)}||\alpha)$

要使得 $f(x^{(k+1)})< f(x^{(k)})$ 成立，当 $\alpha$ 非常小的时候，需要有 $g^{(k)T}H_kg^{(k)}>0$ 成立，最简单的方法就是使得 $H_k$ 是正定的。于是我们得到了黑塞矩阵近似的第一个条件： $H_k$ 正定。

假定目标函数 $f$ 是二次型函数，即黑塞矩阵Q是常正定对称矩阵（保证二次型函数一定有最小值），此时

$g^{(k+1)}-g^{(k)}=Q(x^{(k+1)}-x^{(k)})$

规定

$\Delta g^{(k)}=g^{(k+1)}-g^{(k)}$

$\Delta_{}x^{(k)}=x^{(k+1)}-x^{(k)}$

可得

$\Delta_{}g^{(k)}=Q\Delta_{}x^{(k)}$

则 $Q^{-1}$ 应该满足

$Q^{-1}\Delta_{}g^{(i)}=\Delta_{}x^{(i)}\ \ \ \ \ \ 0\leq i\leq k$

因此，近似矩阵 $H_{k+1}$ 应该满足

$H_{k+1}\Delta_{}g^{(i)}=\Delta_{}x^{(i)}\ \ \ \ \ \ 0\leq i\leq k$

如果迭代n次，总共产生n个迭代方向 $\Delta_{}x^{(0)}$ ， $\Delta_{}x^{(1)}$ ，...， $\Delta_{}x^{(n-1)}$ 。由此可得:

$H_{n}\Delta_{}g^{(0)}=\Delta_{}x^{(0)}$

$H_{n}\Delta_{}g^{(1)}=\Delta_{}x^{(1)}$

...

$H_{n}\Delta_{}g^{(n-1)}=\Delta_{}x^{(n-1)}$

如果 $\Delta_{}g^{(0)}$ ， $\Delta_{}g^{(1)}$ ，...， $\Delta_{}g^{(n-1)}$ 线性无关，则可唯一确定 $H_{n}$ ，此时又 $Q^{-1}$ 满足上式，由解的唯一性可知 $H_{n}=Q^{-1}$ ，即此时拟牛顿法成为了真正的牛顿法。由此可得近似矩阵 $H_{k+1}$ 的第二个条件，满足拟牛顿方程

$H_{k+1}\Delta_{}g^{(i)}=\Delta_{}x^{(i)}\ \ \ \ \ \ 0\leq i\leq k$

可以用数学归纳法证明，拟牛顿法也是一种共轭方向法，即证 $d^{(0)},d^{(1)},...,d^{(k+1)}$ 是关于Q共轭的即可，这里不作证明。由此可知拟牛顿法在解决二次型函数的最小值问题时具有二次终止性，即至多n步迭代可得最优值。

第三个要求是一个跟实际计算相关的要求，即近似矩阵 $H_{k}$ 到 $H_{k+1}$ 的修改幅度不能过大，在 $H_{k}$ 的基础上加一个低秩的矩阵进行修正即可，这样做的目的是为了减少计算复杂度。

从矩阵 $H_{k}$ 所满足的方程来看，矩阵 $H_{k}$ 不能唯一确定。因此，有许多算法可以合理设计出 $H_{k}$ ，本文给出三种修正公式。

3.秩一修正算法

在秩一修正算法中，修正项为 $\alpha_{k}z^{(k)}z^{(k)T}$ ， $\alpha_{k}\in\mathbb{R}$ ， $z^{(k)}\in\mathbb{R}^{n}$ ，是一个对称矩阵，近似矩阵的更新方程为

$H_{k+1}=H_{k}+\alpha_{k}z^{(k)}z^{(k)T}$

其中，由代数知识容易证明， $rank(z^{(k)}z^{(k)T})=1$ ，故称为秩1修正算法，显然，如果 $H_k$ 对称，则 $H_{(k+1)}$ 亦对称。下面给出秩1修正算法的推导过程

令

$H_{k+1}\Delta_{}g^{(k)}=(H_k+\alpha_kz^{(k)}z^{(k)T})\Delta_{}g^{(k)}=\Delta_{}x^{(k)}$

整理得

$z^{(k)}=\frac{\Delta_{}x^{(k)}-H_kg^{(k)}}{\alpha_k(z^{(k)T}\Delta_{}g^{(k)}}$

据此可得

$\alpha_kz^{(k)}z^{(k)T}=\frac{(\Delta_{}x^{(k)}-H_kg^{(k)})(\Delta_{}x-H_kg^{(k)})^T}{\alpha_k(z^{(k)T}\Delta_{}g^{(k)})^2}$

接下来把 $\alpha_k(z^{(k)T}\Delta_{}g^{(k)})^2$ 用 $H_k,\Delta_{}g^{(k)},\Delta_{}x^{(k)}$ 表示即可，利用 $H_{k+1}\Delta_{}g^{(k)}=(H_k+\alpha_kz^{(k)}z^{(k)T})\Delta_{}g^{(k)}=\Delta_{}x^{(k)}$ ,可得

$\Delta_{}x^{(k)}-H_k\Delta_{}g^{(k)}=(\alpha_kz^{(k)T}\Delta_{}g^{(k)})z^{(k)}$

两边同左乘 $\Delta_{}g^{(k)T}$ ，可得

$\Delta_{}g^{(k)T}\Delta_{}x^{(k)}-\Delta_{}g^{(k)T}H_k\Delta_{}g^{(k)}=\alpha_k(z^{(k)T}\Delta_{}g^{(k)})^2$

替换 $\alpha_kz^{(k)}z^{(k)T}=\frac{(\Delta_{}x^{(k)}-H_kg^{(k)})(\Delta_{}x-H_kg^{(k)})^T}{\alpha_k(z^{(k)T}\Delta_{}g^{(k)})^2}$ 中的分母可得最终的秩一修正算法

$H_{k+1}=H_{k}+\frac{(\Delta_{}x^{(k)}-H_k\Delta_{}g^{(k)})(\Delta_{}x^{(k)}-H_k\Delta_{}g^{(k)})^T}{\Delta_{}g^{(k)T}(\Delta_{}x^{(k)}-H_k\Delta_{}g^{(k)})}$

最终我们得出了秩1算法中近似矩阵 $H_k$ 的更新公式，最初近似矩阵 $H_0$ 使用单位矩阵 $I_n$ ，这样第一步就是最速下降法，往后才是秩1修正算法。

秩1算法步骤如下：

1.令 $k=0$ ;选择初始点 $x^{(0)}$ ，任选一个对称实矩阵 $H_0$ （一般选择位矩阵 $I_n$ ）

2.如果 $g^{(k)}=0$ ，停止迭代；否则，令 $d^{(k)}=-H_kg^{(k)}$ 。

3.计算

$\alpha_k=argmin_{\alpha\geq 0}f(x^{(k)}+\alpha_{}d^{(k)})$

$x^{(k+1)}=x^{(k)}+\alpha_kd^{(k)}$

4.计算

$\Delta_{}x^{(k)}=\alpha_kd^{(k)}$

$\Delta_{}g^{(k)}=g^{(k+1)}-g^{(k)}$

$H_{k+1}=H_{k}+\frac{(\Delta_{}x^{(k)}-H_k\Delta_{}g^{(k)})(\Delta_{}x^{(k)}-H_k\Delta_{}g^{(k)})^T}{\Delta_{}g^{(k)T}(\Delta_{}x^{(k)}-H_k\Delta_{}g^{(k)})}$

5.令 $k=k+1$ ，转第二步

这里有个问题，秩1算法是由 $H_{k+1}\Delta_{}g^{(k)}=\Delta_{}x^{(k)}$ 推导出来的，但希望满足 $H_{k+1}\Delta_{}g^{(i)}=\Delta_{}x^{(i)}\ \ \ \ \ \ 0\leq i\leq k$ ,幸运的是，秩1算法具有继承性，其近似矩阵 $H_k$ 能自动地满足上述方程。由数学归纳法可以轻松证明之，这里不过多赘述。

秩1算法的matlab代码如下：

function [xmin,fval] = rand1_QNM(f,g,x0,tolerance)
%此函数使用秩一修正公式求解二次型函数的极小值，
%f是待求的二次型函数，g为二次型函数的梯度，x0是初值，tolerance是退出循环的最小范围，Q是黑塞矩阵
k=1;
n=length(x0);
x=x0;   %储存初始点x0
d=zeros(n,1);   %储存方向向量d_k
H=eye(n,n);    %初始拟牛顿矩阵选择单位矩阵，即用最速下降法选取方式
y=[f(x0)];
while norm(g(x),2)>=tolerance
    d=-H*g(x);
    alpha_1=10;   %设置alpha_1和alpha_2，以用割线法求出一个极小值
    alpha_2=9;
    while abs(g(x+alpha_2*d)'*d)>abs(g(x)'*d)*0.5  %使用强Wolfe条件选取步长alpha
        g_1=g(x+alpha_1*d)'*d;
        g_2=g(x+alpha_2*d)'*d;
        delta_a=alpha_2-alpha_1;
        if g_2-g_1~=0  %若g_2-g_1不等于零时，直接使用割线法公式
            alpha_1=alpha_2;
            alpha_2=alpha_2-g_2*delta_a/(g_2-g_1);  
        else   
            alpha_1=alpha_2;
            alpha_2=alpha_2-g_2*delta_a/(g_2-g_1+10^(-5));   %若g_2-g_1等于零，则在分母上加上一个非常小的常数，继续迭代
        end
    end
    alpha=alpha_2;
    x_2=x+alpha*d;
    if mod(k,6)==0 %再启动rand1方法，效率会降低，但能保证稳定性
        H=eye(n,n);
    else
        delta_x=x_2-x;
        delta_g=g(x_2)-g(x);
       if delta_g'*(delta_x-H*delta_g)>0
            H=H+(delta_x-H*delta_g)*(delta_x-H*delta_g)'/(delta_g'*(delta_x-H*delta_g)); %秩一修正公式
        elseif delta_g'*(delta_x-H*delta_g)==0
            H=H+(delta_x-H*delta_g)*(delta_x-H*delta_g)'/(delta_g'*(delta_x-H*delta_g)+10^-5);  %若分母为零，则添加一个非常小的常数epsilon，保证分母有意义
        else
            H=eye(n,n);   %若delta_g'*(delta_x-H*delta_g)<0,求得的H_(k+1)有可能非正定，此时使用最速下降法，即H_(k+1)为单位矩阵。
       end
    end
    x=x_2;
    y(k+1)=f(x);
    k=k+1;
end
xmin=x;
fval=f(xmin);
disp(['使用步数k=',num2str(k-1)])
semilogy(0:1:k-1,y,'-*r')
title('秩一修正优化效果图')
xlabel ('迭代次数')
ylabel ('优化函数值')
end

测试函数:

$f(x)=\frac{1}{2}x^T\begin{bmatrix} 2 &0 \\ 0& 1 \end{bmatrix}x+3$

初始值 $x^{(0)}=\begin{bmatrix} 1 & 2 \end{bmatrix}^T$

测试结果：

需要注意的是，秩1修正算法不能完全令人满意。首先，该算法产生的矩阵 $H_{k+1}$ ，并不能保证一定是正定的，这将导致 $d^{(k+1)}$ 不再是下降方向。即使对于二次型问题，这种情况仍有可能发生，其次如果 $\Delta_{}g^{(k)T}(\Delta_{}x^{(k)}-H_k\Delta_{}g^{(k)})$ 接近于0，也会带来计算的困难。

例如：

$f(x)=(x_2-x_1)^4+12x_1x_2-x_1+x_2-3$

在初始点

$x_0=\begin{pmatrix} -0.5262 &0.6014 \end{pmatrix}^T$

初始近似矩阵

$H_0=\begin{bmatrix} 0.1186 &-0.0376 \\ -0.0376 & 0.1191 \end{bmatrix}$

情况下，得到的 $H_1$ 就是非正定的。

下列介绍的DFP算法和BFGS算法都可以在一维搜索精确的情况下，保证每个近似矩阵 $H_k$ 都正定。

4.DFP算法

DFP算法步骤如下：

1.令 $k=0$ ;选择初始点 $x^{(0)}$ ，任选一个对称实矩阵 $H_0$ （一般选择位矩阵 $I_n$ ）

2.如果 $g^{(k)}=0$ ，停止迭代；否则，令 $d^{(k)}=-H_kg^{(k)}$ 。

3.计算

$\alpha_k=argmin_{\alpha\geq 0}f(x^{(k)}+\alpha_{}d^{(k)})$

$x^{(k+1)}=x^{(k)}+\alpha_kd^{(k)}$

4.计算

$\Delta_{}x^{(k)}=\alpha_kd^{(k)}$

$\Delta_{}g^{(k)}=g^{(k+1)}-g^{(k)}$

$H_{k+1}=H_k+\frac{\Delta_{}x^{(k)}\Delta_{}x^{(k)T}}{\Delta_{}x^{(k)T}\Delta_{}g^{(k)}}-\frac{[H_k\Delta_{}g^{(k)}][H_k\Delta_{}g^{(k)}]^T}{\Delta_{}g^{(k)T}H_k\Delta_{}g^{(k)}}$

5.令 $k=k+1$ ，转第二步

同秩1修正法，利用数学归纳法可以轻松证明DFP算法满足拟牛顿方程

$H_{k+1}\Delta_{}g^{(i)}=\Delta_{}x^{(i)}\ \ \ \ \ \ 0\leq i\leq k$

并且可以证明，只要 $H_k$ 正定， $H_{k+1}$ 一定正定只需要利用

$H_{k+1}=H_k+\frac{\Delta_{}x^{(k)}\Delta_{}x^{(k)T}}{\Delta_{}x^{(k)T}\Delta_{}g^{(k)}}-\frac{[H_k\Delta_{}g^{(k)}][H_k\Delta_{}g^{(k)}]^T}{\Delta_{}g^{(k)T}H_k\Delta_{}g^{(k)}}$

再对任意非零向量 $x$ ，有 $x^TH_{k+1}x>0$ 即可，在此不作证明。

DFP算法的maltab代码如下：

function [x_min,fval] = DFP_QNM(f,g,x0,tolerance)
%此函数使用DFP算法求解函数极小值，在一维搜索的时候使用精确线搜索，以割线法求得线搜索方向的极小值
n=length(x0);   %获取向量x0长度，以构建对应规模的矩阵
x=x0;
y=[f(x0)];
H=eye(n,n);  %初始H0使用单位矩阵，即使用最速下降法作为初始拟牛顿矩阵
k=1;
while norm(g(x),2)>=tolerance   %停机准则为x_k处的梯度二范数小于给定容忍值
    d=-H*g(x);  %求出第k次迭代的方向d_k
    alpha_1=10;   %设置alpha_1和alpha_2，以用割线法求出一个极小值
    alpha_2=9;
    while abs(g(x+alpha_2*d)'*d)>abs(g(x)'*d)*0.5  %使用强Wolfe条件选取步长alpha
        g_1=g(x+alpha_1*d)'*d;
        g_2=g(x+alpha_2*d)'*d;
        delta_a=alpha_2-alpha_1;
        if g_2-g_1~=0  %若g_2-g_1不等于零时，直接使用割线法公式
            alpha_1=alpha_2;
            alpha_2=alpha_2-g_2*delta_a/(g_2-g_1);  
        else   
            alpha_1=alpha_2;
            alpha_2=alpha_2-g_2*delta_a/(g_2-g_1+10^(-5));   %若g_2-g_1等于零，则在分母上加上一个非常小的常数，继续迭代
        end
    end
    alpha=alpha_2;
    x_2=x+alpha*d;   %确定x_k+1的值
    %if mod(k,6)==0  再启动DFP方法，但效率会降低
    %    H=eye(n,n);
    %else
        delta_x=x_2-x;
        delta_g=g(x_2)-g(x);
        H=H+delta_x*delta_x'/(delta_x'*delta_g)-(H*delta_g)*(H*delta_g)'/(delta_g'*H*delta_g);   %以DFP算法确定H_k+1
    %end
    x=x_2;
    y(k+1)=f(x);
    k=k+1;
end
x_min=x;
fval=f(x);
disp(['迭代次数k=',num2str(k-1)]);
semilogy(0:k-1,y,'-*r')
title('DFP算法求解函数极小值')
xlabel('迭代次数')
ylabel('优化函数值')
end

测试函数：

$f(x)=\frac{1}{2}x^T\begin{bmatrix} 4 & 2\\ 2& 2 \end{bmatrix}x-x^T\begin{bmatrix} -1\\1 \end{bmatrix}$

初始点为 $x^{(0)}=\begin{bmatrix} 0 &0 \end{bmatrix}^T$

DFP算法能够使得矩阵 $H_k$ 保持正定，因此优于秩1算法。但是，当处理一些较大的非二次型问题时，DFP有时会被“卡住”，使得迭代无法展开。产生这一现象的原因是矩阵 $H_k$ 接近成为奇异矩阵了。但BFGS算法可以有效避免这个问题。

5.BFGS算法

在DFP算法中，根据 $H_{k+1}\Delta_{}g^{(i)}=\Delta_{}x^{(i)}\ \ \ \ \ \ 0\leq i\leq k$ 和秩2的条件，得到

$H_{k+1}=H_k+\frac{\Delta_{}x^{(k)}\Delta_{}x^{(k)T}}{\Delta_{}x^{(k)T}\Delta_{}g^{(k)}}-\frac{[H_k\Delta_{}g^{(k)}][H_k\Delta_{}g^{(k)}]^T}{\Delta_{}g^{(k)T}H_k\Delta_{}g^{(k)}}$

这里的 $H_k$ 是为了近似黑塞矩阵 $Q^{-1}$ 而得到的，但我们还可以构造矩阵Q的近似矩阵。即令 $B_k$ 为第k次迭代中关于矩阵Q的近似，则 $B_{k+1}$ 应该满足

$\Delta_{}g^{(i)}=B_{k+1}\Delta_{}x^{(i)}\ \ \ \ \ 0\leq i\leq k$

可以看到，这与 $H_{k+1}\Delta_{}g^{(i)}=\Delta_{}x^{(i)}\ \ \ \ \ \ 0\leq i\leq k$ 非常相似，唯一区别在于 $\Delta_{}x^{(i)}$ 与 $\Delta_{}g^{(i)}$ 交换了位置。因此，给定关于 $H_k$ 的更新公式，交换 $\Delta_{}x^{(i)}$ 与 $\Delta_{}g^{(i)}$ 的位置，并将 $H_k$ 替换为 $B_k$ 即得到 $B_k$ 的更新公式。

$H^{DFP}_{k+1}=H_k+\frac{\Delta_{}x^{(k)}\Delta_{}x^{(k)T}}{\Delta_{}x^{(k)T}\Delta_{}g^{(k)}}-\frac{[H_k\Delta_{}g^{(k)}][H_k\Delta_{}g^{(k)}]^T}{\Delta_{}g^{(k)T}H_k\Delta_{}g^{(k)}}$

$B_{k+1}=B_k+\frac{\Delta_{}g^{(k)}\Delta_{}g^{(k)T}}{\Delta_{}g^{(k)T}\Delta_{}x^{(k)}}-\frac{[B_k\Delta_{}x^{(k)}][B_k\Delta_{}x^{(k)}]^T}{\Delta_{}x^{(k)T}B_k\Delta_{}x^{(k)}}$

基于上述公式，可知在BFGS算法中，为获得黑塞矩阵逆矩阵的近似矩阵更新公式，只需对矩阵 $B_{k+1}$ 求逆即可：

$H^{BFGS}_{k+1}=B_{k+1}^{-1}=(B_k+\frac{\Delta_{}g^{(k)}\Delta_{}g^{(k)T}}{\Delta_{}g^{(k)T}\Delta_{}x^{(k)}}-\frac{[B_k\Delta_{}x^{(k)}][B_k\Delta_{}x^{(k)}]^T}{\Delta_{}x^{(k)T}B_k\Delta_{}x^{(k)}})^{-1}$

为化简上式，引入公式：

引理（Sherman-Morrison公式）：如果矩阵A非奇异，u和v为列向量，满足 $1+v^TA^{-1}u\neq 0$ ，那么 $A+uv^T$ 非奇异，其逆矩阵可以用 $A^{-1}$ 表示：

$(A+uv^T)^{-1}=A^{-1}-\frac{(A^{-1}u)(v^TA^{-1})}{1+v^TA^{-1}u}$

利用该引理，可以将上述式子化简，最终形式如下：

$H^{BFGS}_{k+1}=H_{k}+(1+\frac{\Delta_{}g^{(k)T}H_k\Delta_{}g^{(k)}}{\Delta_{}g^{(k)T}\Delta_{}x^{(k)}})\frac{\Delta_{}x^{(k)}\Delta_{}x^{(k)T}}{\Delta_{}x^{(k)T}\Delta_{}g^{(k)}}-\frac{H_k\Delta_{}g^{(k)}\Delta_{}x^{(k)T}+(H_k\Delta_{}g^{(k)}\Delta_{}x^{(k)T})^T}{\Delta_{}g^{(k)T}\Delta_{}x^{(k)}}$

BFGS算法计算步骤如下：

1.令 $k=0$ ;选择初始点 $x^{(0)}$ ，任选一个对称实矩阵 $H_0$ （一般选择位矩阵 $I_n$ ）

2.如果 $g^{(k)}=0$ ，停止迭代；否则，令 $d^{(k)}=-H_kg^{(k)}$ 。

3.计算

$\alpha_k=argmin_{\alpha\geq 0}f(x^{(k)}+\alpha_{}d^{(k)})$

$x^{(k+1)}=x^{(k)}+\alpha_kd^{(k)}$

4.计算

$\Delta_{}x^{(k)}=\alpha_kd^{(k)}$

$\Delta_{}g^{(k)}=g^{(k+1)}-g^{(k)}$

5.令 $k=k+1$ ，转第二步

BFGS算法maltab代码如下：

function [x_min,fval] = BFGS_QNM(f,g,x0,tolerance)
%此函数使用BFGS公式求解函数极小值，在一维搜索使用精确线搜索，使用割线法求得最优步长alpha
n=length(x0);   %获取向量x0长度，以构建对应规模的矩阵
x=x0;
y=[f(x0)];
H=eye(n,n);  %初始H0使用单位矩阵，即使用最速下降法作为初始拟牛顿矩阵
k=1;
while norm(g(x),2)>=tolerance   %停机准则为x_k处的梯度二范数小于给定容忍值
    d=-H*g(x);  %求出第k次迭代的方向d_k
    alpha_1=10;   %设置alpha_1和alpha_2，以用割线法求出一个极小值
    alpha_2=9;
    while abs(g(x+alpha_2*d)'*d)>abs(g(x)'*d)*0.5  %使用强Wolfe条件选取步长alpha
        g_1=g(x+alpha_1*d)'*d;
        g_2=g(x+alpha_2*d)'*d;
        delta_a=alpha_2-alpha_1;
        if g_2-g_1~=0  %若g_2-g_1不等于零时，直接使用割线法公式
            alpha_1=alpha_2;
            alpha_2=alpha_2-g_2*delta_a/(g_2-g_1);  
        else   
            alpha_1=alpha_2;
            alpha_2=alpha_2-g_2*delta_a/(g_2-g_1+10^(-5));   %若g_2-g_1等于零，则在分母上加上一个非常小的常数，继续迭代
        end
    end
    alpha=alpha_2;
    x_2=x+alpha*d;   %确定x_k+1的值
    %if mod(k,6)==0  再启动BFGS方法，但效率会降低
    %    H=eye(n,n);
    %else
        delta_x=x_2-x;
        delta_g=g(x_2)-g(x);
        H=H+(1+(delta_g'*H*delta_g)/(delta_g'*delta_x))*(delta_x*delta_x')/(delta_x'*delta_g)-((H*delta_g*delta_x')+(H*delta_g*delta_x')')/(delta_g'*delta_x);    %BFGS的拟牛顿矩阵更新公式
    %end
        x=x_2;
        y(k+1)=f(x);
        k=k+1; 
end
x_min=x;
fval=f(x);
disp(['迭代次数k=',num2str(k-1)]);
semilogy(0:1:k-1,y,'-*r')
title('BFGS算法求解函数极小值')
xlabel('迭代次数')
ylabel('优化函数值')
end

测试函数：

$f(x)=\frac{1}{2}x^T\begin{bmatrix} 4 & 2\\ 2& 2 \end{bmatrix}x-x^T\begin{bmatrix} -1\\1 \end{bmatrix}$

初始点为 $x^{(0)}=\begin{bmatrix} 0 &0 \end{bmatrix}^T$

由于BFGS算法是由DFP算法对偶而来，因此保持了DFP算法正定的性质，以及拟牛顿方程具有继承性，共轭方向等特性。并且，当迭代过程中一维搜索的精度不高时，BFGS算法仍然比较稳健。

上述几种算法的实际表现符合预期，对于二次型问题可以做到至多n步找到最优值。但需要指出的是，在实际设计算法时，要考虑一维搜索停机准则，建议使用非精确的线搜索以减少计算量，秩1算法还需要考虑 $H_k$ 非正定的情况，这种情况可以使用最速下降法的方式保证可以继续迭代，这种情况下并不能保证拟牛顿法的二次终止性，但能极大节省计算资源和避免错误。总之，在实际运用时，还需要考虑计算复杂度和分母为零等各种情况，需要根据情况进行改进设计。

6.补充说明

更新了秩1修正、DFP和BFGS算法，进一步地减少了其空间复杂度和时间复杂度。对于秩1修正算法，由于其不能保证拟牛顿矩阵 $H_k$ 的正定性，当出现这种情况是，使用最速下降法替代，使得每一步具有下降性质，并每隔6次使用一次最速下降法以保证其算法的稳定性。这里使用了Rosenbrock函数进行测试，结果如下：

测试函数:

$f(x,y)=100(y-x^2)^2+(1-x)^2$