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

æ°å¼ä¼åï¼Numerical Optimizationï¼(3)-çé¡¿æ³

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

假设函数 $f(x)$ 的二阶导数 $\nabla ^2 f(x)$ 连续，函数 $f(x)$ 在 $x_k$ 处的二阶泰勒展开为

$f(x) \approx f_k + \nabla f(x_k)(x-x_k) + \frac{1}{2}(x-x_k)^T \nabla ^2f(x_k) (x-x_k)$

其中 $f_k:=f(x_k)$ ，求函数的驻点那就是求导并令导数为零，即

$\nabla f(x) = \nabla f(x_k) + \nabla ^2f(x_k)(x-x_k) = 0$

如果二阶导数非奇异，可以得到下一个迭代点为（上式求出来的 $x$ 就是 $x_{k+1}$ ）

$x_{k+1}=x_k - \nabla^2f(x_k)^{-1} \nabla f(x_k)$

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

$\nabla^2 f(x_k) d = -\nabla f(x_k)$

后计算下一个迭代点 $x_{k+1} = x_k + d_k$ 。

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

初值设置：初始点以及终止准则
检验是否满足终止准则
计算二阶导数，确定搜索方向 $d_k$ : $\nabla ^2 f(x_k) d = - \nabla f(x_k)$
计算下一个迭代点 $x_{k+1}=x_k + d_k$ ，回到步骤2

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

修正牛顿法

最基础的改进是在基本牛顿法中加入线搜索方法求得步长 $\alpha_k$ 且令 $x_{k+1} = x_k + \alpha _k d_k$ ，这种方法称为阻尼牛顿法。

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

$B_k = \nabla ^2f(x_k) + E_k >0$

其中 $E_k$ 为修正阵。

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

输入：初始点，终止阈值

循环：

找到修正阵 $E_k$ 使得 $B_k$ 正定
求解线性方程组 $B_k d = -\nabla f(x_k)$ 得到下降方向 $d_k$
线搜法计算步长 $\alpha_k$
更新迭代点 $x_{k+1}=x_k + \alpha_k d_k$

可以看到修正阵 $E_k$ 的选择对算法起关键作用，针对这个有不少修正的方案，这里简要介绍一下基于 Cholesky 分解法的思想，这种算法在对二阶导数矩阵分解过程中调整对角元使得修正后的二阶导数充分正定，也就是说在矩阵的 $LDL^T$ 分解中， $D$ 的对角元不小于某一个给定常数；并且如果原矩阵正定，那么修正后的二阶导数矩阵也是就原矩阵。

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

养生的控制人：数值优化（Numerical Optimization）(2)-信赖域法38 赞同 · 7 评论文章

拟牛顿法

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

假设函数 $f$ 二次连续可微，在 $x_{k+1}$ 的二次近似为

$f(x) \approx f_{k+1} + \nabla f(x_{k+1}) (x-x_{k+1}) + \frac{1}{2}(x-x_{k+1})^T \nabla^2 f(x_{k+1}) (x-x_{k+1})$

对上式两边求导可得

$\nabla f(x) \approx \nabla f(x_{k+1}) +\nabla^2 f(x_{k+1}) (x-x_{k+1})$

如果令 $x=x_k$ 可得

$\nabla f(x_k) = \nabla f(x_{k+1}) + \nabla^2 f(x_{k+1})(x_k - x_{k+1})$

等价于

$\nabla ^2 f(x_{k+1})(x_{k+1}-x_k)=\nabla f(x_{k+1}) - \nabla f(x_k)$

假设二阶导数矩阵的逆矩阵 $\nabla ^2 f(x_{k+1})$ 近似为 $H_{k+1}$ 满足

$H_{k+1} (\nabla f(x_{k+1})- \nabla f(x_k)) = x_{k+1} -x_k$

上式也称为拟牛顿方程，可以看到 $H_{k+1}$ 和迭代点的位移 $x_{k+1}-x_k$ 和梯度差 $\nabla f_{k+1}- \nabla f_k$ 决定。

DFP

下面介绍一下第一个拟牛顿法，DFP算法，算法中假设 $H_{k+1}$ 由 $H_k$ 修正得到，且修正矩阵为秩二矩阵

$E_k = auu^T + b vv^T$

根据假设 $H_{k+1}=H_k + E_k$ 以及拟牛顿方程可以得到

$s_k = H_k y_k + auu^T y_k + bvv^T y_k$

其中 $s_k = x_{k+1}-x_k$ 为位移， $y_k = \nabla f_{k+1} - \nabla f_k$ 为梯度差。这里 $u,v$ 的选择并不是唯一的，可以取 $u=s_k,v=H_k y_k$ ，则

$s_k = H_k y_k +(au^Ty_y) s_k + (bv^T y_k) H_k y_k$

因此有

$au^T y_k = 1, bv^Ty_k =-1$

从而求出 $a,b$ 并带入 $H_{k+1}$ 的修正表达式得到

$H_{k+1}=H_k - \frac{H_k y_k y_k^T H_k}{y_k^T H_k y_k} + \frac{s_ks_k^T}{s_k^T y_k}$

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

DFP算法流程

选择初值 $x_0,H_0$ 以及收敛阈值
检验终止条件
计算搜索方向 $d_k = -H_k \nabla f_k$
确定步长 $\alpha_k$ 和下一个迭代点 $x_{k+1} = x_{k} + \alpha_k d_k$
根据 DFP 校正公式计算矩阵 $H_{k+1}$ ，令 $k=k+1$ 并返回步骤2

BFGS

BFGS算法的推导过程和 DFP 完全类似，令 $B_{k+1} = H_{k+1}^{-1}$ 可以得到拟牛顿方程的另一个表达式

$B_{k+1}s_k = y_k$

得到的更新公式为

$B_{k+1} = B_k + \frac{y_k y_k^T}{y_k^T s_k} - \frac{B_k s_k s_k^T B_k^T}{s_k^T B_k s_k}$

BFGS算法流程

设置初值 $x_0,B_0$ 和收敛阈值
求解线性方程 $B_k d_k = - \nabla f_k$
线搜索得到步长 $\alpha_k$
令 $s_k = \alpha_k d_k$ 更新 $x_{k+1} = x_k +s_k$
更新梯度差 $y_k = \nabla f_{k+1} - \nabla f_k$
计算矩阵 $B_{k+1}$

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

$B_{k+1}^{-1} = (I-\frac{s_ky_y^T}{y_k^Ts_k})B_k^{-1}(I-\frac{y_ks_k^T}{y_k^Ts_k})+\frac{s_ks_k^T}{y_k^Ts_k}$

一种更有效的计算为

$B_{k+1}^{-1}=B_{k}^{-1}+\frac{(s_k^Ty_k+y_k^TB_k^{-1}y_k)(s_ks_k^T)}{(s_k^Ty_k)^2} - \frac{B_k^{-1}y_ks_k^T+s_k y_k^TB_k^{-1}}{s_k^Ty_k}$

L-BFGS

对于 BFGS 算法，需要储存近似逆二阶导数矩阵 $B_k^{-1}$ ，对于维度较大的问题不再适用，因此有了内存受限的 BFGS 算法 （Limited-memory BFGS）。虽然 L-BFGS 不需要储存近似逆矩阵，但要保存每次迭代的中间信息，不过都是一维数组，且迭代次数不会有很多，所以对储存要求大大降低。

定义

$H_k = B_k ^{-1},\rho_k=\frac{1}{y_k^Ts_k}, V_k = I-\rho_k y_ks_k^T$

则 BFGS 的公式可以改写为

$H_{k+1} = V_k^TH_kV_k+\rho_ks_ks_k^T$

在第 $k$ 次迭代，当前点为 $x_k$ ，且存有 $\{s_i,y_i\},i=k-m,...,k-1$ 总共 $m$ 步的位移和梯度差。选择一个初始的 $H_k^0$ ，则可以推导出 $H_k$ 的表达式为

$H_k=\left( V_{k-1}^{T}...V_{k-m}^{T} \right) H_{k}^{0}\left( V_{k-m}...V_{k-1} \right) \\ +\rho _{k-m}\left( V_{k-1}^{T}...V_{k-m+1}^{T} \right) s_{k-m}s_{k-m}^{T}\left( V_{k-m+1}...V_{k-1} \right) \\ +\rho _{k-m+1}\left( V_{k-1}^{T}...V_{k-m+2}^{2} \right) s_{k-m+1}s_{k-m+1}^{T}\left( V_{k-m+2}...V_{k-1} \right) \\ +... \\ +\rho _{k-1}s_{k-1}s_{k-1}^{T}$

根据这个表达式我们可以推导出计算 $H_k \nabla f_k$ 的递归算法

令 $q=\nabla f_k$

循环1： $i=k-1,...,k-m$

$\alpha_i = \rho_i s_i^Tq$
$q=q-\alpha_i y_i$

令 $r = H_k^0 q$

循环2： $i=k-m,k-m+1,...,k-1$

$\beta = \rho_i y_i^T r$
$r =r + s_i (\alpha_i -\beta)$

输出 $H_k \nabla f_k =r$

上面的递归算法涉及了 $H_k^0$ 的选择，一种有效的方式是

$H_k^0 = \gamma_k I, \gamma_k = \frac{s_{k-1}^Ty_{k-1}}{y_{k-1}^Ty_{k-1}}$

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

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

输入：初始点 $x_0$ ，记忆步长 $m$ ，令 $k=0$

循环直到收敛：

选择 $H_k^0$
根据双循环递归算法计算 $d_k = -H_k \nabla f_k$
计算下一迭代点 $x_{k+1}=x_{k}+\alpha_k d_k$ 其中 $\alpha_k$ 的选择需要满足 Wolfe 条件
如果 $k>m$
计算并保存 $s_k = x_{k+1}-x_k$ ， $y_{k}=\nabla f_{k+1}-\nabla f_k$
令 $k=k+1$

MATLAB示例

求 Rosenbrock函数的最小点

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

计算梯度为

$\nabla f=\left[ \begin{array}{c} -2\left( 1-x \right) -400x\left( y-x^2 \right)\\ 200\left( y-x^2 \right)\\ \end{array} \right]$

计算二阶导数为

$\nabla ^2f=\left[ \begin{matrix} 2-400\left( y-x^2 \right) +800x^2& -400x\\ -400x& 200\\ \end{matrix} \right]$

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

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 函数，初始点为 $(4,14)$ ，L-BFGS的记忆步数为 30 步。

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

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

养生的控制人

修正牛顿法

拟牛顿法

DFP

BFGS

L-BFGS

MATLAB示例

为你推荐

06-darknet函数图谱

【强化学习与最优控制】笔记（七） Rollout 与 Policy Improvement

Pytorch深度学习（一）

「RL篇肆」Vanilla Policy Gradient 原理与实战

MATLAB神经网络（六）：GRNN

深入理解深度学习——语境词嵌入（Contextual Word Embedding）

关于作者

养生的控制人

14

0

0

1

高斯过程

共轭梯度法简介

卡尔曼估计两步法

相关推荐

【强化学习与最优控制】笔记（九）值函数，Q函数和策略空间的近似

Faster rcnn 训练安全帽数据集搭载服务器运行

ML1. 线性回归(Linear Regression) -2

【强化学习与最优控制】笔记（八）模型预测控制（Model Predictive Control）

深度学习 pytorch cifar10数据集训练

K近邻算法和KD树详细介绍及其原理详解

热门泡泡

30积分失眠，聊聊自己搞ROS的心得体会吧

ros学习路线

30积分 TF_REPEATED_DATA ignoring data错误

各位大佬，有什么ROS定位算法推荐吗

5积分想买能用ROS2的开发套件。或者开发板

5积分 ros中启动gazebo时报错

给作者打赏

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

养生的控制人

修正牛顿法

拟牛顿法

DFP

BFGS

L-BFGS

MATLAB示例

为你推荐

06-darknet函数图谱

【强化学习与最优控制】笔记（七） Rollout 与 Policy Improvement

Pytorch深度学习（一）

「RL篇 肆」Vanilla Policy Gradient 原理与实战

MATLAB神经网络（六）：GRNN

深入理解深度学习——语境词嵌入（Contextual Word Embedding）

评论（0）

关于作者

养生的控制人

14

0

0

1

高斯过程

共轭梯度法简介

卡尔曼估计两步法

相关推荐

【强化学习与最优控制】笔记（九）值函数，Q函数和策略空间的近似

Faster rcnn 训练 安全帽数据集搭载服务器运行

ML1. 线性回归(Linear Regression) -2

【强化学习与最优控制】笔记（八） 模型预测控制（Model Predictive Control）

深度学习 pytorch cifar10数据集训练

K近邻算法和KD树详细介绍及其原理详解

热门泡泡

30积分 失眠，聊聊自己搞ROS的心得体会吧

ros学习路线

30积分 TF_REPEATED_DATA ignoring data错误

各位大佬，有什么ROS定位算法推荐吗

5积分 想买能用ROS2的开发套件。或者开发板

5积分 ros中启动gazebo时报错

给作者打赏

忘记密码

修改头像

添加你感兴趣的标签

举报类型（必选）

举报详情（选填）

「RL篇肆」Vanilla Policy Gradient 原理与实战

Faster rcnn 训练安全帽数据集搭载服务器运行

【强化学习与最优控制】笔记（八）模型预测控制（Model Predictive Control）

30积分失眠，聊聊自己搞ROS的心得体会吧

5积分想买能用ROS2的开发套件。或者开发板