一.前言

上上上篇文章介绍了牛顿法和修正牛顿法。想看的话可以往后翻。牛顿法有二阶的收敛速度，但Hess阵必须要正定，因为只有正定才能保证它的下降方向是正确的。虽然修正牛顿法克服了这个缺点，但是它的修正参数uk的选取很难，选得太大太小都会影响到算法的收敛速度。最重要的是，这两种方法都要用到Hess阵，计算量很大，所以都不实用。本文介绍的拟牛顿法很好的克服了这些缺点。

二.拟牛顿法的基本思想

拟牛顿法的基本思想是把牛顿法中用到的Hess阵用一个Hk矩阵来代替。那么Hk是什么呢？Hk的三个特点如下：

（1）Hk近似等于牛顿法中的Hess阵，这样可以保证拟牛顿发所产生的方向与牛顿反向近似，从而保证了拟牛顿法的收敛速度。

（2）Hk是正定的

（3）Hk的更新规则有两种，即用秩1或秩2的矩阵矫正。这篇文章介绍的是秩1算法，秩2算法下篇文章再写吧。（秩1秩2法其实就是两个更新Hk的不同公式，应该是记住就行了吧，反正推导我也看不懂）

三.秩1矫正Hk公式

在做matlab仿真时，H0通常用一个单位矩阵来代替。经过一次迭代后，Hk的秩1矫正公式如下：

$H^{_{k+1}}=H_{k}+\frac{(S_{k }-H_{k}y_{k})(S_{k }-H_{k}y_{k})'}{(S_{k}-H_{k}y_{k})'y_{k}}$

其中Sk=Xk+1-Xk(就是本次的迭代点减去上一个迭代点)，yk=gk+1-gk(就是本次的梯度值减去上一次的梯度值)

四.算法步骤

步0：确定终止误差e=（0~1），设初始点x0， $\delta$ =(0~1)， $\sigma$ =（0,0.5），初始对称正定阵H0=I（单位阵），令k=0

步1：计算gk= $\bigtriangledown$ f(xk).若||gk||<=e,停算，输出xk作为最优解。否则，转步2

步2：计算搜索方向：dk=-Hk*gk

步3：用Armjio搜索技术求步长 $\alpha$ k= $\delta$ ^mk，m的值从0开始

若f(xk+ $\delta$ ^m*dk)<=f(xk)+ $\alpha$ * $\delta$ ^m*gk'dk

则 mk=m，步长 $\alpha$ k= $\delta$ ^mk，若不满足上式，则m=m+1，直到满足上述不等式为止

步4：令Xk+1=xk+ $\alpha$ k*dk，由矫正公式来确定Hk+1

$H^{_{k+1}}=H_{k}+\frac{(S_{k }-H_{k}y_{k})(S_{k }-H_{k}y_{k})'}{(S_{k}-H_{k}y_{k})'y_{k}}$

步5：k=k+1，转步1

五.代码实现

1.秩1矫正算法

function [x,val,k]=sr1(fun,gfun,x0)
%功能：用秩1拟牛顿法求无约束问题 mini f(x)
%输入：fun,gfun分别是目标函数和梯度，x0是初始点
%输出：x,val分别是近似最优点和最优值，k表示迭代次数
k=0;
maxk=500;
rho=0.55;
sigma=0.4;
e=1e-5;%精度
n=length(x0);
Hk=eye(n);
while(k<maxk)gk=feval(gfun,x0);dk=-Hk*gk;if(norm(gk)<=e),break;end
m=0;
mk=0;
while(m<20)if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk'*dk)mk=m;break;endm=m+1;end
x=x0+dk*rho^mk;
sk=x-x0;
yk=feval(gfun,x)-gk;
Hk=Hk+(sk-Hk*yk)*(sk-Hk*yk)'/((sk-Hk*yk)'*yk);
k=k+1;
x0=x;
end
x0=x;
val=feval(fun,x0);
end

2.目标函数

function f= fun(x)
%目标函数
f=100*(x(1)^2-x(2))^2+(x(1)-1)^2;
end

3.目标函数梯度

function  g=gfun(x)
%目标函数的梯度
g=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1),-200*(x(1)^2-x(2))]';
end

4.主函数

%这个问题的精确值是x=(1,1)',f(x)=0;
clear all
clc
x0=[-1.2,1]';
[x,val,k]=sr1('fun','gfun',x0);
disp('迭代次数:k=')
disp(k)
disp(['最优解：x = '])
disp(x)
disp(['此时: f(x) = ',num2str(val)])