牛顿迭代法

Goodness2020

已于 2022-03-23 08:49:59 修改

阅读量4k

点赞数 3

分类专栏：优化文章标签：机器学习

于 2022-03-19 09:32:59 首次发布

本文链接：https://blog.csdn.net/Goodness2020/article/details/123589164

版权

优化专栏收录该内容

8 篇文章 0 订阅

订阅专栏

一、应用方向

1、求方程的根

并不是所有的方程 f(x)=0 都有求根公式，或者求根公式很复杂，这时可以利用牛顿法，迭代求解。

原理是利用泰勒公式在 $x_{0}$ 处展开，且展开到一阶，即 $f(x)=f(x_{0})+(x-x_{0})f^{'}(x_{0})$ ，求解 f(x)=0 ，即 $f(x_{0})+(x-x_{0})f^{'}(x_{0})=0$ ，化简得： $x=x_{0}-\tfrac{f(x_{0})}{f^{'}(x_{0})}$ ，因此获得迭代公式如下：

$x_{k+1}=x_{k}-\frac{f(x_{k})}{f^{'}(x_{k})}$

2、最优化

在最优化问题中，线性最优化可以用行法求解，但对于非线性最优化问题，牛顿法提供了一种求解的办法。假设任务是优化一个目标函数 f(x) ，求函数的极大极小问题，可以转化为求解函数的导数 $f^{'}(x)=0$ 的问题，这样就可以把优化问题转化为方程求解问题。

若要求方程 $f^{'}(x)=0$ 的根，则迭代公式如下：

$x_{k+1}=x_{k}-\frac{f^{'}(x_{k})}{f^{''}(x_{k})}$

二、牛顿法

考虑无约束最优化问题： min f(x) ，其中， $x\in \mathbf{R}^{n}$ ， $x^{*}$ 为目标函数的极小点。

假设 f(x) 有二阶连续偏导数，若第k次迭代值为 $x^{(k)}$ ，则可将在附近进行二阶泰勒展开：

$f(x)=f(x^{(k)})+g(x^{(k)})^{T}(x-x^{(k)})+\frac{1}{2}(x-x^{(k)})^{T}H(x^{(k)})(x-x^{(k)})$

这里， $g(x^{(k)})=\triangledown f(x^{(k)})$ 是 f(x) 的梯度向量在点 $x^{(k)}$ 的值； $H(x^{(k)})$ 是 f(x) 的黑塞矩阵

$H(x) = \left [ \frac{\partial ^{2} f}{\partial x_{i} \partial x_{j}} \right ]_{n\times n}$

在点 $x^{(k)}$ 的值。函数有极值点的必要条件是在极值点处一阶导数为0，即梯度向量为0， $\triangledown f(x)=0$ 。特别是当 H(x) 在极值点处是正定矩阵时，函数 f(x) 的极值为极小值。

由二阶泰勒展开式，两边对求导，得 $\triangledown f(x)=g(x^{(k)})+H(x^{(k)})(x-x^{(k)})$ ；再假设 $x^{(k+1)}$ 满足 $\triangledown f(x^{(k+1)})=0$ ，则 $g(x^{(k)})+H(x^{(k)})(x^{(k+1)}-x^{(k)})=0$ ，因此，迭代公式如下：

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

其中， $-g(x^{(k)})=-\triangledown f(x^{(k)})$ 为迭代方向， $H(x^{(k)})$ 为迭代步长。

完整算法如下：

三、代码实现

matlab实现求正定二次函数的极小值点

% 牛顿法---收敛速度很快；但是每一步都需要计算Hessian矩阵的逆，计算复杂
% 牛顿法---基于二阶导数，基于当前位置的切线来确定下一次的位置

function [x,m] = newton(A,b,x,ep,n)
% A:正定矩阵  b:向量  x:初始向量  ep:停止误差  n:最大迭代次数
% 求二次函数 f(x)=(1/2)x'*A*x+b'x+c 的极小值，A为正定矩阵
% A为f(x)的Hessian矩阵  f'(x)=Ax+b为梯度向量
m = 0; % 统计迭代次数
for k = 0:n
    g = A * x + b; % 梯度向量
    if (g'*g <= ep) % 终止条件
        break;
    end
    H = A; % hessian矩阵
    x = x - H^(-1) * g; % 迭代公式
    m = m + 1;
end

% A = [2,1;1,2];
% b = [1;2];
% x = [0;0];
% ep = 0.00001;
% n = 1000;
% [y,m] = newton(A,b,x,ep,n);  % [0;-1]