牛顿迭代与高斯牛顿法

一元函数下

对于一元函数 \(f(x)\) ，其泰勒二阶展开如下为

\[f(x_{k+1})=f(x_k)+f'(x_k)(x_{k+1}-x_k)+\frac{1}{2}f''(x_k)(x_{k+1}-x_k)^2+O(x^2).\]

记 \(z=x_{k+1}-x_k\) 省略余项可表示为

\[f(x_{k+1})\approx f(x_k)+f'(x_k)\cdot z+\frac{1}{2}f''(x_k)\cdot z^2=g(x).\]

  • 若要求 \(f(x)\) 的零点，即想要 \(f(x_{k+1})=0\)，省略二阶展开易得

\[f(x_k)+f'(x_k)\cdot z=0,\]

\[z=-\frac{f(x_k)}{f'(x_k)},\]

\[x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}.\]

  • 若要求 \(f(x)\) 的极值点，即想要 \(f'(x_{k+1})=0\)，也就是 \(g'(z)=0\) ，有

\[g'(z)=f'(x_k)+f''(x_k)\cdot z=0，\]

\[z=-\frac{f'(x_k)}{f''(x_k)},\]

\[x_{k+1}=x_k-\frac{f'(x_k)}{f''(x_k)}.\]

多元函数下

基本过程为：

  1. 二阶逼近
  2. 逼近式对增量求一阶导为零
  3. 求出增量计算格式
  4. 代回得到迭代式

二阶逼近为

\[f\left(\mathbf{x}_{k}+\mathbf{z}\right) \approx f\left(\mathbf{x}_{k}\right)+\left[\nabla f\left(\mathbf{x}_{k}\right)\right]^{T} \mathbf{z}+\frac{1}{2} \mathbf{z}^{T}\left[\nabla^{2} f\left(\mathbf{x}_{k}\right)\right] \mathbf{z}.\]

其中 \(\mathbf{x}=\left(x_1,x_2\cdots x_n\right)^T\)，\(\mathbf{z}=\left(z_1,z_2\cdots z_n\right)^T.\) 有

\[\nabla f\left(\mathbf{x}_{k}\right)=\left(\begin{array}{cccc} \frac{\partial f(\mathbf{x})}{\partial x_{1}} & \frac{\partial f(\mathbf{x})}{\partial x_{2}} & \cdots & \frac{\partial f(\mathbf{x})}{\partial x_{n}} \end{array}\right)^T.\]

又有

\[\nabla^{2} \mathbf{f}=\nabla \otimes \nabla \mathbf{f}=\left[\begin{array}{cccc} \frac{\partial^{2} f}{\partial^{2} x_{1}} & \frac{\partial^{2} f}{\partial x_{1} \partial x_{2}} & \cdots & \frac{\partial^{2} f}{\partial x_{1} \partial x_{n}} \\ \frac{\partial^{2} f}{\partial x_{2} \partial x_{1}} & \frac{\partial^{2} f}{\partial x_{2} \partial x_{2}} & \cdots & \frac{\partial^{2} f}{\partial x_{2} \partial x_{n}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^{2} f}{\partial x_{n} \partial x_{1}} & \cdots & \ldots & \frac{\partial^{2} f}{\partial^{2} x_{n}} \end{array}\right]\]

这就是Hessian矩阵，记为 \(\mathbf{H}_{f}(\mathbf{x})\).

令一阶导为 \(0\)，得到

\[\nabla f\left(\mathbf{x}_{k}\right)+\left[\nabla^{2} f\left(\mathbf{x}_{k}\right)\right] \mathbf{z}=0,\]

\[\mathbf{z}=-\mathbf{H}_{f}^{-1}\left(\mathbf{x}_{k}\right) \nabla f\left(\mathbf{x}_{k}\right),\]

\[\mathbf{x}_{k+1}=\mathbf{x}_{k}-\mathbf{H}_{f}^{-1}\left(\mathbf{x}_{k}\right) \nabla f\left(\mathbf{x}_{k}\right).\]

因为这里Hessian矩阵不一定可逆，最简单直接的思路是在 \(\mathbf{H}^{-1}_{f}(\mathbf{x})\) 的对角线上加上一个正数：

\[\mathbf{x}_{k+1}=\mathbf{x}_{k}-\left(\mathbf{H}_{f}^{-1}\left(\mathbf{x}_{k}\right)+\mu \mathbf{I}\right) \nabla f\left(\mathbf{x}_{k}\right).\]

求最小二乘问题

记损失函数为

\[\ell(\mathbf{w})=\frac{1}{2}\|f(\mathbf{X} , \mathbf{w})-\mathbf{Y}\|^{2}=\frac{1}{2} \mathbf{e}^{2}.\]

其中 \(\mathbf{w}\) 是全体参数组成的向量， \(\mathbb{e}\) 是误差值组成的向量。其二阶展开式如下

\[\ell\left(\mathbf{w}_{k}+\mathbf{z}\right) \approx \ell\left(\mathbf{w}_{k}\right)+\left(\nabla \ell\left(\mathbf{w}_{k}\right)\right)^{T} \mathbf{z}+\frac{1}{2} \mathbf{z}^{T}\left(\nabla^{2} \ell\left(\mathbf{w}_{k}\right)\right) \mathbf{z}.\]

因为最小二乘法的特殊性，可以简化为

\[\nabla \ell(\mathbf{w})=\frac{\partial\left(\frac{1}{2} \mathbf{e}^{2}\right)}{\partial \mathbf{w}}=\left(\frac{\partial \mathbf{e}}{\partial \mathbf{w}}\right)^{T} \mathbf{e}=\mathbf{J}_{\mathbf{e}}^{T}(\mathbf{w}) \mathbf{e}.\]

这就是Jacobian矩阵即为 \(\mathbf{J}_{\mathbf{e}}^{}(\mathbf{w})\).注意始终保证最后的结果是列向量。求二阶导为

\[\nabla^{2} \ell(\mathbf{w})=\nabla(\nabla \ell(\mathbf{w}))=\frac{\partial\left(\left(\frac{\partial \mathbf{e}}{\partial \mathbf{w}}\right)^{T} \mathbf{e}\right)}{\partial \mathbf{w}}=\left(\frac{\partial \mathbf{e}}{\partial \mathbf{w}}\right)^{T}\left(\frac{\partial \mathbf{e}}{\partial \mathbf{w}}\right)+\frac{\partial^{2} \mathbf{e}}{\partial \mathbf{w}^{2}} \mathbf{e}\]

\[=\mathbf{J}_{\mathbf{e}}^{T}(\mathbf{w}) \mathbf{J}_{\mathbf{e}}(\mathbf{w})+\frac{\partial^{2} \mathbf{e}}{\partial \mathbf{w}^{2}} \mathbf{e}\]

代回得到迭代式

\[\mathbf{w}_{k+1}=\mathbf{w}_{k}-\left.\left(\left(\mathbf{J}_{\mathbf{e}}^{T}(\mathbf{w}) \mathbf{J}_{\mathbf{e}}(\mathbf{w})+\frac{\partial^{2} \mathbf{e}}{\partial \mathbf{w}^{2}} \mathbf{e}\right)^{-1} \mathbf{J}_{\mathbf{e}}^{T}(\mathbf{w}) \mathbf{e}\right)\right|_{\mathbf{w}=\mathbf{w}_{k}}\]

高斯-牛顿法

在文章 Theoria motus corporum coelestium in sectionibus conicis solem ambientum (1809), pp. 179-180. 中，高斯提到 \(\nabla^{2} \ell(\mathbf{w})\) 的二阶项应该很小，所以直接省略，于是再对符号简化一下，这个Hessian近似就可以写成：

\[\nabla^{2} \ell(\mathbf{w})=\mathbf{J}^T\mathbf{J}.\]

这个式子就是著名的 Gauss-Newton 法的迭代公式。但 \(\mathbf{J}^T\mathbf{J}\) 通常是半正定的，但不一定可逆的于是在对角线上加上一个正数，得到 Levenberg-Marquardt 算法。

\[\mathbf{w}_{k+1}=\mathbf{w}_{k}-\left(\mathbf{J}^{T} \mathbf{J}+\mu \mathbf{I}\right)^{-1} \mathbf{J}^{T} \mathbf{e}\]

代码

一元函数

求核酸混检问题中， \[E(x)=1-(1-p)^k-\frac{1}{k}\] 的最小值点。

clc;clear;
p=0.01;
syms k; %自变量
E=@(k) 1-(1-p).^k+1./k; %每个人检测次数的期望
kDif=matlabFunction(diff(E,k,1));
kDiif=matlabFunction(diff(E,k,2));
%% 牛顿迭代法
x=[0.1];%当x过大kDiif(x)为负，函数为凸，会倾向于寻找极大值
while abs(kDif(x(end)))>1e-10
    z=-kDif(x(end))/kDiif(x(end));
    x(end+1)=abs(x(end)+z);%因为kDif(x)是偶函数
end
fprintf('ans=%d\n',x(end));
%% 画图
xx = -10:0.1:100;
yy=kDif(xx);
plot(xx,yy,'LineWidth',3);
hold on;
for i=1:size(x,2)-1;
    plot([x(i),x(i+1)],[kDif(x(i)),0]);
    plot([x(i+1),x(i+1)],[0,kDif(x(i+1))]);
end
axis([-10,100,-0.1,1]);

参考

本文几乎全部复制于知乎回答，仅作学习使用。