1 拟牛顿法
牛顿法:$$x_{k+1} = x_k - \nabla^2 f(x_k)^{-1} \nabla f(x_k)$$
牛顿法的缺陷在于 \(\nabla^2 f(x_k)\) 难以求取和存储,因此有拟牛顿法(Quasi-Newton methods),通过在牛顿法的迭代中加入近似求取每一步\(Hessian\)矩阵的迭代步,仅通过迭代点处的梯度信息来求取\(Hessian\)矩阵的近似值。
我们把\(Hessian\)矩阵在第\(k\)步的迭代的近似值记作\(B_k\),通过\(B_k\)代替\(\nabla^2 f(x_k)\),得到新的迭代式:$$x_{k+1} = x_k - B_k^{-1} \nabla f(x_k)$$
为了保证\(B_{k+1}\)能较好地近似第\(k+1\)步的\(Hessian\)矩阵,令\(H_k = B_k^{-1}\),可以得到以下两个式子:
\[\begin{aligned} x_{k+1} = x_k - H_k \nabla f(x_k) \\ H_{k+1} = f(H_k) \end{aligned}\]注意: 拟牛顿法是一种思想而不是具体方法,只要能满足上述迭代式并且收敛,那就是拟牛顿法
2 算法框架
令\(H_0 = I\),给出迭代初始值\(x_0\)和迭代停止阈值\(\epsilon > 0\)
while(1) {
计算梯度:\(\nabla f(x_k)\);
if(\(\| \nabla f(x_k) \|\) < \(\epsilon\))
break;
计算搜索方向:\(d_k = -H_k \nabla f(x_k)\);
更新迭代点:\(x_{k+1} = x_k - d_k\);
根据\(x_k\)处的信息更新\(H_k\);
}
为了观察\(f(x_{k+1})和f(x_k)\)的差值,使用泰勒展开:$$f(x_{k+1}) = f(x_k) + g_k^T d + o(|d|^2)$$
这里默认\(g_k^T<0\),但因为这里\(x_k\)和\(g_k\)(迭代点\(x_k\)处的梯度)固定,只能改变 d,因此问题转化为$$\max_d |g_k^Td|$$
当\(\cos<g_k, d>\)固定时,只要\(\|d\|^2\)够大,\(|g_k^Td|\)会变成\(\infty\),需要进行约束。规定\(\|d\|_{B_k}\), 定义为\(\|d^TB_kd\|\),则:$$\max_d |g_k^Td| \qquad s.t. |d|_{B_k} = 1$$
由 \(Cauchy-Schwartz\)不等式,
当且仅当\(d = -B_k^{-1}g_k\), "="
3 如何确定\(H_k\)的迭代更新?
设\(s_k = x_{k+1} - x_k, y_k = \nabla f(x_{k+1}) - \nabla f(x_k)\),则$$B_{k+1}s_k = y_k$$
\[H_{k+1}y_k = s_k \]SR1
根据\(x_k\)处的信息得到一个修正量\(\Delta H_k\)来更新:$$H_{k+1} = H_k + \Delta H_k$$
由于我们希望\(H_k\)接近\(\nabla^2 f(x_k)^{-1}\), \(H_{k+1}\)接近\(\nabla^2 f(x_{k+1})^{-1}\),有$$\Delta H_k = \nabla2f(x_{k+1}) - \nabla2f(x_k)$$
由于\(\nabla^2f(x_{k+1})^{-1} 和 \nabla^2f(x_k)^{-1}\)对称,所以\(\Delta H_k\)对称。
该算法设\(\delta H_k = \beta uu^T\),那么迭代式为$$H_{k+1} = H_k + \beta uu^T \tag{1}$$
其中\(\beta \in \mathbb{R}, u \in \mathbb{R}^n\)
在\((1)\)两边乘上\(y_k\),再根据\(H_{k+1}y_k = s_k\),有$$H_{k+1}y_k = H_ky_k + \beta uu^T y_k = s_k$$
\[\rightarrow s_k - H_ky_k = \beta (u^Ty_k)u \tag{2} \]而\(\beta (u^T y_k)\)是实数,有 $$u = \frac{1}{\beta (u^T y_k)} (s_k - H_ky_k)$$
令\(\gamma = \frac{1}{\beta (u^T y_k)}\),所以$$u = \gamma (s_k - H_ky_k) \tag{3}$$
(2)可以表示成\(s_k - H_ky_k = \beta uu^T y_k\),将(3)代入,有
注意到\(\beta \gamma^2 \in \mathbb{R}\) 和 \((s_k - H_ky_k)^Ty_k \in \mathbb{R}\),所以$$s_k - H_ky_k = [\beta \gamma^2(s_k - H_k y_k)^T y_k](s_k - H_ky_k)$$
显然当\(\beta \gamma^2(s_k - H_k y_k)^T y_k = 1\)时成立,则$$\beta \gamma^2 = \frac{1}{(s_k-H_ky_k)^Ty_k} \tag{5}$$
将(3)(5)代入(1),有
令\(H_0 = I\),给出迭代初始值\(x_0\)和迭代停止阈值\(\epsilon > 0\)
while(1) {
计算梯度:\(\nabla f(x_k)\);
if(\(\| \nabla f(x_k) \|\) < \(\epsilon\))
计算搜索方向:\(d_k = -H_k \nabla f(x_k)\);
更新迭代点:\(x_{k+1} = x_k - d_k\);
更新\(H_{k+1} = H_k + \frac{(s_k-H_ky_k)(s_k-H_ky_k)^T}{(s_k-H_ky_k)^Ty_k}\);
}
DFP
其基本思路同\(SR1\)类似,我们令\(\Delta H_k = \beta uu^T + \gamma vv^T\),有$$H_{k+1} = H_k + \beta uu^T + \gamma vv^T \tag{1} $$
\(\beta\) 和 \(\gamma\)是待定系数
在上式两边右乘\(y_k\): $$s_k = H_{k+1}y_k = H_ky_k + \beta uu^Ty_k + \gamma vv^T y_k \tag{2}$$
取\(u = s_k, v = H_ky_k\)(取的过程我不知道┭┮﹏┭┮),改写(2): $$s_k - H_ky_k = \beta uu^T y_k + \gamma vv^T y_k$$
令\(\beta uu^T y_k = s_k, \gamma vv^T y_k = - H_ky_k\),所以 $$u = s_k = \beta uu^Ty_k = (\beta u^Ty_k)u$$
为了上式恒成立,有\(\beta u^Ty_k = 1\),因此$$\beta = \frac{1}{u^Ty_k} = \frac{1}{s_k^Ty_k}$$
同理, \(\gamma = - \frac{1}{v^Ty_k} = - \frac{1}{y^T_kH_ky_k}\)
将\(\beta, \gamma, u, v\)代入(1),我们可以得到迭代式 $$H_{k+1} = H_k + \frac{s_ks_kT}{s_KTy_k} - \frac{H_ky_ky_kTH_k}{y_kTH_ky_k}$$
令\(H_0 = I\),给出迭代初始值\(x_0\)和迭代停止阈值\(\epsilon > 0\)
while(1) {
计算梯度:\(\nabla f(x_k)\);
if(\(\| \nabla f(x_k) \|\) < \(\epsilon\))
计算搜索方向:\(d_k = -H_k \nabla f(x_k)\);
更新迭代点:\(x_{k+1} = x_k - d_k\);
更新\(H_{k+1} = H_k + \frac{s_ks_k^T}{s_K^Ty_k} - \frac{H_ky_ky_k^TH_k}{y_k^TH_ky_k}\);
}
BFGS
考虑对\(B_k\) (对\(\nabla^2 f(x_k)\)的近似) 的修正,同DFP得到$$B_{k+1} = B_{k} + \frac{y_ky_kT}{y_kTs_k} - \frac{B_ks_ks_kTB_k}{s_kTB_ks_k} \tag{1} $$
引理1(矩阵求逆引理)
\[\begin{aligned} (A^{-1}+BD^{-1}C)^{-1} & \equiv A - AB(D+CAB)^{-1}CA \\ (D+CAB)^{-1} & \equiv D^{-1} - D^{-1}C(A^{-1}+BD^{-1}C)^{-1}BD^{-1} \\ AB(D+CAB)^{-1} & \equiv (A^{-1} + BD^{-1}C)^{-1}BD^{-1} \\ (D+CAB)^{-1}CA & \equiv D^{-1}C(A^{-1}+BD^{-1}C)^{-1}\end{aligned} \]证明
对一个矩阵分别进行LDU和UDL分解:
\[\left[ \begin{matrix} A^{-1} & -B \\ C & D \end{matrix} \right] = \left[ \begin{matrix} 1 & 0 \\ CA & 1 \end{matrix} \right] \left[ \begin{matrix} A^{-1} & 0 \\ 0 & D+CAB \end{matrix} \right] \left[ \begin{matrix} 1 & -AB \\ 0 & 1 \end{matrix} \right] \tag{LDU} \]\[\left[ \begin{matrix} A^{-1} & -B \\ C & D \end{matrix} \right] = \left[ \begin{matrix} 1 & -BD^{-1} \\ 0 & 1 \end{matrix} \right] \left[ \begin{matrix} A^{-1}+BD^{-1}C & 0 \\ 0 & D \end{matrix} \right] \left[ \begin{matrix} 1 & 0 \\ D^{-1}C & 1 \end{matrix} \right] \tag{ULD} \]两边同时求逆
\[\begin{aligned} (\left[ \begin{matrix} A^{-1} & -B \\ C & D \end{matrix} \right])^{-1} & = \left[ \begin{matrix} 1 & AB \\ 0 & 1 \end{matrix} \right] \left[ \begin{matrix} A & 0 \\ 0 & (D+CAB)^{-1} \end{matrix} \right] \left[ \begin{matrix} 1 & 0 \\ -CA & 1 \end{matrix} \right] \\ & = \left[ \begin{matrix} A-AB(D+CAB)^{-1}CA & AB(D+CAB)^{-1} \\ -(D+CAB)^{-1}CA & (D+CAB)^{-1} \end{matrix} \right] \end{aligned} \]\[\begin{aligned} (\left[ \begin{matrix} A^{-1} & -B \\ C & D \end{matrix} \right])^{-1} & = \left[ \begin{matrix} 1 & 0 \\ -D^{-1}C & 1 \end{matrix} \right] \left[ \begin{matrix} (A^{-1}+BD^{-1}C)^{-1} & 0 \\ 0 & D^{-1} \end{matrix} \right] \left[ \begin{matrix} 1 & BD^{-1} \\ 0 & 1 \end{matrix} \right] \\ & = \left[ \begin{matrix} (A^{-1}+BD^{-1}C)^{-1} & (A^{-1}+BD^{-1}C)^{-1}BD^{-1} \\ -D^{-1}C(A^{-1}+BD^{-1}C)^{-1} & D^{-1}-D^{-1}C(A^{-1}+BD^{-1}C)^{-1}BD^{-1} \end{matrix} \right] \end{aligned} \]参考 《机器人学中的状态估计》 @深蓝学院
递归使用上面的引理1,令\(A = B_k + \frac{y_ky_k^T}{y_k^Ts_k},u = - \frac{B_ks_k}{s_k^TB_ks_k}, v = B_ks_k\)代入SMW公式,其中\(A^{-1}\)的计算再用一次SMW,可以得到
\[\begin{aligned} H_{k+1} = B^{-1}_{k+1} & = (B_k + \frac{y_ky_k^T}{y_k^Ts_k} - \frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k})^{-1} \\ & = H_k + (1+\frac{y_k^TH_ky_k}{s_k^Ty_k}) \frac{s_ks_k^T}{s_k^Ty_k}-(\frac{s_ky_k^TH_k+H_ky_ks_k^T}{s_k^Ty_k}) \end{aligned} \]令\(H_0 = I\),给出迭代初始值\(x_0\)和迭代停止阈值\(\epsilon > 0\)
标签:begin,end,matrix,nabla,牛顿,beta,ky From: https://www.cnblogs.com/wanyy-home/p/17904350.htmlwhile(1) {
计算梯度:\(\nabla f(x_k)\);
if(\(\| \nabla f(x_k) \|\) < \(\epsilon\))
计算搜索方向:\(d_k = -H_k \nabla f(x_k)\);
更新迭代点:\(x_{k+1} = x_k - d_k\);
更新\(H_{k+1} = H_k + (1+\frac{y_k^TH_ky_k}{s_k^Ty_k}) \frac{s_ks_k^T}{s_k^Ty_k}-(\frac{s_ky_k^TH_k+H_ky_ks_k^T}{s_k^Ty_k})\);
}