牛顿迭代与高斯牛顿法

请注意,文章中部分理解可能已被笔者废弃,本文最后更新于:6 个月前

一元函数下

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

f(xk+1)=f(xk)+f(xk)(xk+1xk)+12f(xk)(xk+1xk)2+O(x2).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=xk+1xkz=x_{k+1}-x_k 省略余项可表示为

f(xk+1)f(xk)+f(xk)z+12f(xk)z2=g(x).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)零点,即想要 f(xk+1)=0f(x_{k+1})=0,省略二阶展开易得

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

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

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

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

g(z)=f(xk)+f(xk)z=0g'(z)=f'(x_k)+f''(x_k)\cdot z=0,

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

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

多元函数下

基本过程为:

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

二阶逼近为

f(xk+z)f(xk)+[f(xk)]Tz+12zT[2f(xk)]z.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}.

其中 x=(x1,x2xn)T\mathbf{x}=\left(x_1,x_2\cdots x_n\right)^Tz=(z1,z2zn)T.\mathbf{z}=\left(z_1,z_2\cdots z_n\right)^T.

f(xk)=(f(x)x1f(x)x2f(x)xn)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.

又有

2f=f=[2f2x12fx1x22fx1xn2fx2x12fx2x22fx2xn2fxnx12f2xn]\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矩阵,记为 Hf(x)\mathbf{H}_{f}(\mathbf{x}).

令一阶导为 00,得到

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

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

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

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

xk+1=xk(Hf1(xk)+μI)f(xk).\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).

求最小二乘问题

记损失函数为

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

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

(wk+z)(wk)+((wk))Tz+12zT(2(wk))z.\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}.

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

(w)=(12e2)w=(ew)Te=JeT(w)e.\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矩阵即为 Je(w)\mathbf{J}_{\mathbf{e}}^{}(\mathbf{w}).注意始终保证最后的结果是列向量。求二阶导为

2(w)=((w))=((ew)Te)w=(ew)T(ew)+2ew2e\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}

=JeT(w)Je(w)+2ew2e=\mathbf{J}_{\mathbf{e}}^{T}(\mathbf{w}) \mathbf{J}_{\mathbf{e}}(\mathbf{w})+\frac{\partial^{2} \mathbf{e}}{\partial \mathbf{w}^{2}} \mathbf{e}

代回得到迭代式

wk+1=wk((JeT(w)Je(w)+2ew2e)1JeT(w)e)w=wk\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. 中,高斯提到 2(w)\nabla^{2} \ell(\mathbf{w}) 的二阶项应该很小,所以直接省略,于是再对符号简化一下,这个Hessian近似就可以写成:

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

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

wk+1=wk(JTJ+μI)1JTe\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(1p)k1kE(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]);

参考

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