matlab牛顿法_数值优化（Numerical Optimization）(3)-牛顿法

3基本牛顿法

牛顿法的基本思想是用迭代点的梯度信息和二阶导数对目标函数进行二次函数逼近，然后把二次函数的极小值作为新的迭代点，并不断重复这一过程，直到求出极小点。

假设函数

的二阶导数

连续，函数

在

处的二阶泰勒展开为

其中

，求函数的

驻点那就是求导并令 导数为零，即

如果二阶导数非奇异，可以得到下一个迭代点为（上式求出来的

就是

）

如果二阶导数奇异，那么可以求解下面线性方程确定搜索方向

后计算下一个迭代点

。

基本牛顿法可以归结为以下四步

1. 初值设置：初始点以及终止准则
2. 检验是否满足终止准则
3. 计算二阶导数，确定搜索方向
   :
4. 计算下一个迭代点
   ，回到步骤2

注意：牛顿法的好处在于收敛速度快，缺点在于计算二阶导数的计算量大以及求解线性方程组确定搜索方向可能是病态的。

修正牛顿法

最基础的改进是在基本牛顿法中加入线搜索方法求得步长

且令

，这种方法称为

阻尼牛顿法。

牛顿法面临的一个主要困难是二阶导数不正定，在这一情况下，下降方向就很难获得。Goldfeld 修正法在二阶导数不正定时对其进行修正

其中

为修正阵。

带有线搜索的修正牛顿法可以表述为

输入：初始点，终止阈值

循环：

• 找到修正阵
  使得
  正定
• 求解线性方程组
  得到下降方向
• 线搜法计算步长
• 更新迭代点

可以看到修正阵

的选择对算法起关键作用，针对这个有不少修正的方案，这里简要介绍一下

基于 Cholesky 分解法的思想，这种算法在对二阶导数矩阵分解过程中调整对角元使得修正后的二阶导数充分正定，也就是说在矩阵的

分解中，

的对角元不小于某一个给定常数；并且如果原矩阵正定，那么修正后的二阶导数矩阵也是就原矩阵。

Ps: 信赖域牛顿法在上一篇博文仿真案例用的就是了，感兴趣的可以去看一下。

养生的控制人：数值优化（Numerical Optimization）(2)-信赖域法​zhuanlan.zhihu.com

拟牛顿法

拟牛顿法的思想是模拟牛顿方向的生成路径，利用相邻两个点的位移和一阶导数信息构造与二阶导数阵相似的正定矩阵。所需的计算量比牛顿法少，收敛速度达到超线性。

假设函数

二次连续可微，在

的二次近似为

对上式两边求导可得

如果令

可得

等价于

假设二阶导数矩阵的逆矩阵

近似为

满足

上式也称为拟牛顿方程，可以看到

和迭代点的

位移

和

梯度差

决定。

DFP

下面介绍一下第一个拟牛顿法，DFP算法，算法中假设

由

修正得到，且修正矩阵为秩二矩阵

根据假设

以及拟牛顿方程可以得到

其中

为位移，

为梯度差。这里

的选择并不是唯一的，可以取

，则

因此有

从而求出

并带入

的修正表达式得到

这个公式也称为 DFP校正公式。

DFP算法流程

1. 选择初值
   以及收敛阈值
2. 检验终止条件
3. 计算搜索方向
4. 确定步长
   和下一个迭代点
5. 根据 DFP 校正公式计算矩阵
   ，令
   并返回步骤2

BFGS

BFGS算法的推导过程和 DFP 完全类似，令

可以得到拟牛顿方程的另一个表达式

得到的更新公式为

BFGS算法流程

1. 设置初值
   和收敛阈值
2. 求解线性方程
3. 线搜索得到步长
4. 令
   更新
5. 更新梯度差
6. 计算矩阵

在算法中

可以初始化为单位阵，第一步的矩阵求逆可以根据 Sherman-Morrison 公式进行转化后得到

一种更有效的计算为

L-BFGS

对于 BFGS 算法，需要储存近似逆二阶导数矩阵

，对于维度较大的问题不再适用，因此有了

内存受限的 BFGS 算法 （Limited-memory BFGS）。虽然 L-BFGS 不需要储存近似逆矩阵，但要保存每次迭代的中间信息，不过都是一维数组，且迭代次数不会有很多，所以对储存要求大大降低。

定义

则 BFGS 的公式可以改写为

在第

次迭代，当前点为

，且存有

总共

步的位移和梯度差。选择一个初始的

，则可以推导出

的表达式为

根据这个表达式我们可以推导出计算

的递归算法

令

循环1：

令

循环2：

输出

上面的递归算法涉及了

的选择，一种有效的方式是

这里的

叫做尺度因子用来估计沿最近搜索方向的真实二阶导数矩阵的大小。

综上，完整的 L-BFGS 算法可以描述为

输入：初始点

，记忆步长

，令

循环直到收敛：

• 选择
• 根据双循环递归算法计算
• 计算下一迭代点
  其中
  的选择需要满足 Wolfe 条件
• 如果
• 计算并保存
  ，
• 令

MATLAB示例

求 Rosenbrock函数的最小点

计算梯度为

计算二阶导数为

首先把这个函数写成输入输出的函数形式

function [f,g,H] = RosenFunc(x0)

x = x0(1);
y = x0(2);
f = (1-x)^2 + 100*(y-x^2)^2;
if (nargout > 1)
    g = [-2 + 2*x - 400*x*y + 400*x^3;
        200*(y-x^2)];
end
if (nargout > 2)
   H = [2 - 400*(y-x^2) + 800*x^2, -400*x;
       -400*x, 200];
end

另外也附上 Branin 函数用于测试

function [f,g,H] = BraninFunc(x0)
x1 = x0(1);
x2 = x0(2);
f = (x2-0.129*x1^2+1.6*x1-6)^2+6.07*cos(x1)+10;
if (nargout > 1)
    g = [2*(x2-0.129*x1^2+1.6*x1-6)*(-0.258*x1+1.6)-6.07*sin(x1);...
        2*(x2-0.129*x1^2+1.6*x1-6)];
end
if (nargout > 2)
    H = [2*(-0.258*x1+1.6)^2-0.516*(x2-0.128*x1^2+1.6*x1-6)-6.07*cos(x1),...
    -0.516*x1+3.2;-0.516*x1+3.2,2];
end

阻尼牛顿法（线搜索用回溯法）

function [x_opt,x_eval,f_eval] = Newton_basic(fun,x0,epsilon,iter)
x_eval = [];
f_eval = [];
xk = x0;
[fk,gk,Hk] = fun(xk);
for k = 1:iter  
    if norm(gk,2) <= epsilon
        break
    end
    x_eval = [x_eval,xk];
    f_eval = [f_eval,fk];
    dk = -inv(Hk)*gk;
    alpha = backstracking_linesearch(fun,xk,dk);
    xk = xk + alpha*dk;
    [fk,gk,Hk] = fun(xk);
end
x_opt = xk;

BFGS法（线搜索用回溯法）

function [x_opt,x_eval,f_eval] = Newton_BFGS(fun,x0,epsilon,iter)
x_eval = [];
f_eval = [];
xk = x0;
n = length(x0);
Hk = eye(n);
[fk,gk] = fun(xk);
for k = 1:iter
    if norm(gk,2) <= epsilon
        break
    end
    x_eval = [x_eval,xk];
    f_eval = [f_eval,fk];
    dk = -Hk*gk;
    % line search
    alpha = backstracking_linesearch(fun,xk,dk);
    x_next = xk + alpha*dk;
    [f_next,g_next] = fun(x_next);
    % BFGS
    sk = alpha*dk;
    yk = g_next - gk;
    rho = 1/(yk'*sk);
    Hk = (eye(n) - rho*sk*yk')*Hk*(eye(n)-rho*yk*sk') + rho*sk*sk';
    % 更新
    gk = g_next;
    xk = x_next;
    fk = f_next;
end
x_opt = xk;

L-BFGS法（线搜索用回溯法）

function [x_opt,x_eval,f_eval] = Newton_LBFGS(fun,x0,m,epsilon,iter)
x_eval = [];
f_eval = [];
n = length(x0);
S = zeros(n,m);
Y = zeros(n,m);
xk = x0;
[fk,gk] = fun(x0);
for k = 1:iter
    if norm(gk,2) <= epsilon
        break
    end
    x_eval = [x_eval,xk];
    f_eval = [f_eval,fk];
    if k > 1
        gamma = (S(:,1)'*Y(:,1))/(Y(:,1)'*Y(:,1));
        H0 = gamma*eye(n);
        % two-loop
        q = gk;
        rho = zeros(m,1);
        ALPHA = zeros(m,1);
        for i = 1:min(k-1,m)
            rho(i) = 1/(Y(:,i)'*S(:,i));
            ALPHA(i) = rho(i)*S(:,i)'*q;
            q = q - ALPHA(i)*Y(:,i);
        end
        dk = H0*q;
        for j = 1:min(k-1,m)-1
            i = min(k,m)-j+1;
            beta = rho(i)*Y(:,i)'*dk;
            dk = dk +S(:,i)*(ALPHA(i)-beta);
        end
        dk = -dk;
    else
        H0 = eye(n);
        dk = -H0*gk;
    end
    % line search
    alpha = backstracking_linesearch(fun,xk,dk);
    x_next = xk + alpha*dk;
    [f_next,g_next] = fun(x_next);
    % 更新位移和梯度差
    S(:,2:m) = S(:,1:m-1);
    Y(:,2:m) = Y(:,1:m-1);
    sk = x_next - xk;
    yk = g_next - gk;
    S(:,1) = sk;
    Y(:,1) = yk;
    % 更新
    xk = x_next;
    fk = f_next;
    gk = g_next;
end
x_opt = xk;

选择一个比较有意思的对比展示一下，采用 Branin 函数，初始点为

，L-BFGS的记忆步数为 30 步。

三种算法跑到了三个地方去了，两个局部最小点，基本牛顿法停留在一个在非常平的地段。对应的函数值也不一样。