牛顿迭代法➨牛顿法➨拟牛顿法
2️⃣

牛顿迭代法➨牛顿法➨拟牛顿法

牛顿迭代法
牛顿迭代法是求解方程f(x)=0f(x)=0的近似根的一种方法。
牛顿法
牛顿法是求解无约束最优化问题min⁡f(x)\min f(x)的最优解的一种方法。
拟牛顿法
通过正定矩阵近似Hessian矩阵的逆矩阵或Hessian矩阵，简化了计算牛顿法的计算。
 

牛顿迭代法

牛顿迭代法是求解方程f(x)=0f(x)=0的近似根的一种方法。
选用泰勒展开f(x)=f(x0)+f′(x0)(x−x0)+o(x−x0)2f(x)=f(x_0)+f'(x_0)(x-x_0)+o(x-x_0)^2的前两项来近似方程，即
f(x)≈f(x0)+f′(x0)(x−x0)f(x)\approx f(x_0)+f'(x_0)(x-x_0)
非线性方程线性化求解f(x)=0f(x)=0的近似根 ⟹ f(x0)+f′(x0)(x−x0)≈0f(x_0)+f'(x_0)(x-x_0)\approx 0 ⟹ x≈x0−f(x0)f′(x0)x\approx x_0-\frac{f(x_0)}{f'(x_0)}
牛顿迭代法求
牛顿迭代法求f(x)=0f(x)=0x(k+1)=x(k)−f(x(k))f′(x(k))x^{(k+1)}=x^{(k)}-\frac{f(x^{(k)})}{f'(x^{(k)})}
👁️‍🗨️已经证明：如果f(x)f(x)是连续的，且待求零点是孤立的，那么在零点附近存在一个区域，只要初始点在这个区域内，那么牛顿法必定收敛，并且有平方收敛的能力。
➣➣牛顿迭代法流程➣➣
  • 设置初始近似值x0x_0
  • x(1)=x(0)−f(x(0))f′(x(0))x^{(1)}=x^{(0)}-\frac{f(x^{(0)})}{f'(x^{(0)})} ⇤ 解的 1次近似
  • x(2)=x(1)−f(x(1))f′(x(1))x^{(2)}=x^{(1)}-\frac{f(x^{(1)})}{f'(x^{(1)})} ⇤ 解的2次近似
  • ⋯\cdots
  • x(k+1)=x(k)−f(x(k))f′(x(k))x^{(k+1)}=x^{(k)}-\frac{f(x^{(k)})}{f'(x^{(k)})} ⇤ 解的k+1次近似
牛顿迭代法可看作对原函数的线性逼近
牛顿迭代法可看作对原函数的线性逼近
 
 

牛顿迭代法 ➠ 牛顿法

牛顿法是求解无约束最优化问题min⁡f(x)\min f(x)的最优解的一种方法
一维的公式推导（x∈Rx\in \mathbb{R}）
最优化问题min⁡f(x)\min f(x) ⟹ g(x)=f′(x)=0g(x)=f'(x)=0 ⟹ 牛顿迭代法求解g(x)=0g(x)=0
⟹ 迭代公式x(k+1)=x(k)−g(x(k))g′(x(k))=x(k)−f′(x(k))f′′(x(k))x^{(k+1)}=x^{(k)}-\frac{g(x^{(k)})}{g'(x^{(k)})}=x^{(k)}-\frac{f'(x^{(k)})}{f''(x^{(k)})}
n维的公式推导（x∈Rn,f:Rn→R\mathbf{x}\in \mathbb{R^n},f:\mathbb{R}^n\rightarrow\mathbb{R}）
泰勒展开 f(x(k+1))=f(x(k)+d)=f(x(k)))+∇f(x(k))Td+12dTH(x(k))d+o(∣∣d∣∣2),x(k+1))=x(k)+df(\mathbf{x}^{(k+1)})=f(\mathbf{x}^{(k)}+d)=f(\mathbf{x}^{(k)}))+\nabla f(\mathbf{x}^{(k)})^{T}d+\frac{1}{2}d^TH(\mathbf{x}^{(k)})d+o(||d||^2),\mathbf{x}^{(k+1)})=\mathbf{x}^{(k)}+d
舍去高阶项，f(x(k+1))≈f(x(k))+∇f(x(k))Td+12dTH(x(k))d:=h(d)f(\mathbf{x}^{(k+1)})\approx f(\mathbf{x}^{(k)})+\nabla f(\mathbf{x}^{(k)})^{T}d+\frac{1}{2}d^TH(\mathbf{x}^{(k)})d:=h(d)
min⁡x(k+1)f(x(k+1))≈min⁡dh(d)\min_{\mathbf{x}^{(k+1)}} f(\mathbf{x}^{(k+1)}) \approx \min_d h(d) ⟹ h′(d)=∇f(x(k))T+H(x(k))d=0h'(d)=\nabla f(\mathbf{x}^{(k)})^{T}+H(\mathbf{x}^{(k)})d=0
⟹ d=−H(x(k))−1∇f(x(k))d=-H(\mathbf{x}^{(k)})^{-1}\nabla f(\mathbf{x}^{(k)}) ⟹ x(k+1)=x(k)+d=x(k)−H(x(k))−1∇f(x(k))\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}+d=\mathbf{x}^{(k)}-H(\mathbf{x}^{(k)})^{-1}\nabla f(\mathbf{x}^{(k)})
其中∇f(x)\nabla f(\mathbf{x})为梯度向量，H(x)H(\mathbf{x})为Hessian矩阵：
∇f(x)=(∂f(x)∂x1∂f(x)∂x2⋮∂f(x)∂xn),H(x)=[∂2f(x)∂x12∂2f(x)∂x1∂x2⋯∂2f(x)∂x1∂xn∂2f(x)∂x2∂x1∂2f(x)∂x22⋯∂2f(x)∂x2∂xn⋮⋮⋱⋮∂2f(x)∂xn∂x1∂2f(x)∂xn∂x2⋯∂2f(x)∂xn2],x=(x1x2⋮xn)\nabla f(\mathbf{x})=\left(\begin{array}{cccc} \frac{\partial f(\mathbf{x})}{\partial x_1}\\\frac{\partial f(\mathbf{x})}{\partial x_2}\\\vdots\\\frac{\partial f(\mathbf{x})}{\partial x_n}\end{array}\right),H\left(\mathbf{x}\right)=\left[\begin{array}{cccc} \frac{\partial^2 f\left(\mathbf{x}\right)}{\partial x_1^2} & \frac{\partial^2 f\left(\mathbf{x}\right)}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2 f\left(\mathbf{x}\right)}{\partial x_1 \partial x_n} \\ \frac{\partial^2 f\left(\mathbf{x}\right)}{\partial x_2 \partial x_1} & \frac{\partial^2 f\left(\mathbf{x}\right)}{\partial x_2^2} & \cdots & \frac{\partial^2 f\left(\mathbf{x}\right)}{\partial x_2 \partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2 f\left(\mathbf{x}\right)}{\partial x_n \partial x_1} & \frac{\partial^2 f\left(\mathbf{x}\right)}{\partial x_n \partial x_2} & \cdots & \frac{\partial^2 f\left(\mathbf{x}\right)}{\partial x_n^2} \end{array}\right],\mathbf{x}=\left(\begin{array}{cccc} x_1\\x_2\\\vdots\\x_n\end{array}\right)
牛顿法
牛顿法
假设f(x)f(x)具有连续二阶偏导数，其求无约束问题min⁡f(x)\min f(x)的牛顿法迭代公式为
一维：x(k+1)=x(k)−f′(x(k))])f′′(x(k)), x∈R一维：x^{(k+1)}=x^{(k)}-\frac{f'(x^{(k)})])}{f''(x^{(k)})},\ x\in\mathbb{R} 多维：x(k+1)=x(k+1)H(x(k))−1∇f(x(k)), x∈Rn多维： \mathbf{x}^{(k+1)}=\mathbf{x}^{(k+1)}H(\mathbf{x}^{(k)})^{-1}\nabla f(\mathbf{x}^{(k)}),\ \mathbf{x}\in \mathbb{R}^n
牛顿法可看作对原函数的二次逼近，或对其导函数的线性逼近
牛顿法可看作对原函数的二次逼近，或对其导函数的线性逼近
 
如何理解牛顿法的公式，以及为何牛顿法收敛速度快？
由于目标点是导数为0的点，迭代公式每次迭代改变的步长为d=f′(x(k))f′′(x(k))d=\frac{f'(x^{(k)})}{f''(x^{(k)})}。可以看出，当分子一阶导越大，步长则越长，符合实际意义，因为我们的目标是快速找到导数为0的点，而此时一阶导数还很大，所以改变量可以大一些；另一方面，若分母二阶导越小，步长则越大，也符合实际意义，二阶导越小，意味着一阶导的变化速度慢，故一阶导为0的位置距离当前位置还很远，故该变量可以大一些。基于以上两点，牛顿法的收敛速度很快。
➣➣牛顿法算法流程➣➣
输入：目标函数f(x)f(\mathbf{x})，梯度g(x)=∇f(x)g(\mathbf{x})=\nabla f(\mathbf{x})，Hessian矩阵H(x)H(\mathbf{x})，精度要求ϵ\epsilon;
输出：f(x)f(\mathbf{x})的极小值点x∗\mathbf{x}^*.
（1）取初值点x(0)\mathbf{x}^{(0)}，置k=0k=0;
（2）计算gk=g(x(k))g_k=g(\mathbf{x}^{(k)});
（3）若∣∣gk∣∣≤ϵ||g_k||\leq\epsilon，返回x∗=x(k)\mathbf{x}^*=\mathbf{x}^{(k)};
（4）计算Hk=H(x(k))H_k=H(\mathbf{x}^{(k)})，并计算pk:Hkpk=−gkp_k:H_kp_k=-g_k，即pk=−Hk−1gkp_k=-H_k^{-1}g_k;
（5）置x(k+1)=x(k)+pk,k←k+1\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}+p_k,k\leftarrow k+1，转(2)
 

牛顿法VS梯度下降法

  • 算法思路
      • 梯度下降法是一种基于一阶导数的优化算法，核心思想是通过不断沿着负梯度方向迭代，使得目标函数的值不断逼近最小值点；
      牛顿法是一种基于二阶导数的优化算法，核心思想是通过不断更新当前点的位置，使得目标函数的值不断逼近最小值点。
  • 迭代效率
      • 梯度下降法由于只需计算目标函数的一阶导数，每次迭代的计算量相对较小，适用于处理大规模数据集。
      牛顿法每次迭代都需要同时计算目标函数的一阶和二阶导数，因此每次迭代的计算量较大，这可能导致其在处理大型数据集时的迭代速度受限。
  • 算法收敛性
      • 梯度下降法在多数情况下都能得到收敛的结果，但收敛速度较慢。
      牛顿法具有较快的收敛速度，尤其在处理较少次的迭代时，可以得到很好的结果。
      为什么牛顿法收敛更快，或迭代次数更少？
      红色曲线是利用牛顿法迭代求解，绿色曲线是利用梯度下降法求解。
      红色曲线是利用牛顿法迭代求解，绿色曲线是利用梯度下降法求解。
      通俗来说梯度下降法每次只从你当前所处位置选一个坡度最大的方向走一步，牛顿法在选择方向时，不仅会考虑坡度是否够大，还会考虑你走了一步之后，坡度是否会变得更大。所以，可以说牛顿法比梯度下降法看得更远一点，能更快地走到最底部。
  • 初始点的影响
      • 梯度下降法对初始点的选择不太敏感，算法会沿着负梯度方向不断迭代，最终找到局部最优解。
      牛顿法对初始点的选择较为敏感，不同的初始点可能会导致算法找到不同的极值点。
总结起来，牛顿法收敛更快，但每次的迭代的计算量较大。故牛顿法在处理非凸优化问题和需要快速收敛的情况下更为适用，而梯度下降法则在处理大规模数据和大规模问题的场景下表现出更高的可扩展性和稳定性。
 

牛顿法 ➠ 拟牛顿法

  • 相同：牛顿法和拟牛顿法的收敛速度都比较快；
  • 不同：牛顿法每一步需要求解目标函数Hessian矩阵的逆矩阵，计算比较复杂；
  • 不同：拟牛顿法通过正定矩阵近似Hessian矩阵的逆矩阵或Hessian矩阵，简化了计算。
基本思想：考虑用一个n阶矩阵Gk=G(x(k))G_k=G(\mathbf{x}^{(k)})来近似替代Hk−1=H−1(x(k))H_k^{-1}=H^{-1}(\mathbf{x}^{(k)})
牛顿法中Hk−1H_k^{-1}满足的性质
将f(x)f(\mathbf{x})在x(k+1)\mathbf{x}^{(k+1)}附近作泰勒展开f(x)≈f(x(k+1))+gk+1(x−x(k+1))+12(x−x(k+1))THk+1(x−x(k+1))f(\mathbf{x})\approx f(\mathbf{x}^{(k+1)})+g_{k+1}(\mathbf{x}-\mathbf{x}^{(k+1)})+\frac{1}{2}(\mathbf{x}-\mathbf{x}^{(k+1)})^TH_{k+1}(\mathbf{x}-\mathbf{x}^{(k+1)})，
则∇f(x)=gk+1+Hk+1(x−x(k+1))\nabla f(\mathbf{x})=g_{k+1}+H_{k+1}(\mathbf{x}-\mathbf{x}^{(k+1)})，
记yk=gk+1−gk,δk=x(k+1)−x(k)y_k=g_{k+1}-g_k,\delta_k=\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}，则有yk=Hk+1δk或Hk+1−1yk=δk⏟拟牛顿条件\underbrace{y_k=H_{k+1}\delta_k \text{或}H_{k+1}^{-1}y_k=\delta_k}_{拟牛顿条件}
定理：如果HkH_k是正定的（那么Hk−1H^{-1}_k也是正定的），则可以保持牛顿法的搜索方向pk=−Hk−1gkp_k=-H_k^{-1}g_k一定是下降方向。
证明：若x=x(k)+λpk=x(k)−λHk−1gk\mathbf{x}=\mathbf{x}^{(k)}+\lambda p_k=\mathbf{x}^{(k)}-\lambda H_k^{-1}g_k，则f(x)f(\mathbf{x})在x(k)\mathbf{x}^{(k)}附近（λ\lambda充分小）的泰勒展开可以近似为f(x)≈f(x(k))−λgkTHk−1gkf(\mathbf{x})\approx f(\mathbf{x}^{(k)})-\lambda g_{k}^TH_k^{-1}g_k，由于HkH_k是正定的，那么Hk−1H^{-1}_k也是正定的，则gkTHk−1gk>0g_{k}^TH_k^{-1}g_k>0.当λ\lambda充分小时，总有f(x)<f(x(k))f(\mathbf{x})<f(\mathbf{x}^{(k)})，也就是说pk=−Hk−1gkp_k=-H_k^{-1}g_k是下降方向。
牛顿法中满足的性质
牛顿法中Hk−1H_k^{-1}满足的性质δk=Hk+1−1yk, 其中yk=gk+1−gk,δk=x(k+1)−x(k) \delta_k=H_{k+1}^{-1}y_k,\ 其中y_k=g_{k+1}-g_k,\delta_k=\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}
如果HkH_k是正定的，则可以保持牛顿法的搜索方向pk=−Hk−1gkp_k=-H_k^{-1}g_k一定是下降方向。
拟牛顿法将GkG_k作为Hk−1H_k^{-1}的近似，要求GkG_k满足同样的条件
拟牛顿条件
拟牛顿条件
  1. GkG_k是正定的；
  1. GkG_k满足Gk+1yk=δk⏟拟牛顿条件\underbrace{G_{k+1}y_k=\delta_k}_{拟牛顿条件}
 
GkG_k的选择具有灵活性，有许多种具体方法。
更新:在每次迭代中基于拟牛顿条件更新
更新Gk+1G_{k+1}:在每次迭代中基于拟牛顿条件更新Gk+1=Gk+ΔGkG_{k+1}=G_k+\Delta G_k
 
DFP算法
假设Gk+1=Gk+Pk+QkG_{k+1}=G_k+P_k+Q_k，此时Gk+1yk=Gkyk+Pkyk+QkykG_{k+1}y_k=G_ky_k+P_ky_k+Q_ky_k
为了满足拟牛顿条件，则可使Pkyk=δk, Qkyk=−GkykP_ky_k=\delta_k,\ Q_ky_k=-G_ky_k，例如取Pk=δkδkTδkTyk, Qk=−GkykykTGkykTGkykP_k=\frac{\delta_k\delta_k^T}{\delta_k^T y_k},\ Q_k=-\frac{G_ky_ky_k^TG_k}{y_k^T G_k y_k}，则
Gk+1=Gk+δkδkTδkTyk−GkykykTGkykTGkykG_{k+1}=G_k+\frac{\delta_k\delta_k^T}{\delta_k^T y_k}-\frac{G_ky_ky_k^TG_k}{y_k^T G_k y_k}
可以证明，若初始G0G_0正定，则对任意k都有GkG_k正定。
➣➣DFP算法➣➣
输入：目标函数f(x)f(\mathbf{x})，梯度g(x)=∇f(x)g(\mathbf{x})=\nabla f(\mathbf{x})，精度要求ϵ\epsilon;
输出：f(x)f(\mathbf{x})的极小值点x∗\mathbf{x}^*.
（1）取初值点x(0)\mathbf{x}^{(0)}，取G0G_0为正定对称阵，置k=0k=0;
（2）计算gk=g(x(k))g_k=g(\mathbf{x}^{(k)})，若∣∣gk∣∣≤ϵ||g_k||\leq\epsilon，返回x∗=x(k)\mathbf{x}^*=\mathbf{x}^{(k)}；否则，转（3）
（3）置pk=−Gk−1gkp_k=-G_k^{-1}g_k;
（4）一维搜索：求λk\lambda_k使得f(x(k)+λkpk)=min⁡λf(x(k)+λpk)f(\mathbf{x}^{(k)}+\lambda_kp_k)=\min_\lambda f(\mathbf{x}^{(k)}+\lambda p_k);
（5）置x(k+1)=x(k)+λkpk\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}+\lambda_k p_k;
（6）计算gk+1=g(x(k+1))g_{k+1}=g(\mathbf{x}^{(k+1)})，若∣∣gk+1∣∣≤ϵ||g_{k+1}||\leq\epsilon，返回x∗=x(k+1)\mathbf{x}^*=\mathbf{x}^{(k+1)}；否则，更新Gk+1=Gk+δkδkTδkTyk−GkykykTGkykTGkykG_{k+1}=G_k+\frac{\delta_k\delta_k^T}{\delta_k^T y_k}-\frac{G_ky_ky_k^TG_k}{y_k^T G_k y_k}；
（7）置k←k+1k\leftarrow k+1，转（3）。
 
BFGS算法
可以将GkG_k作为Hk−1H_k^{-1}的近似，也可以将BkB_k作为HkH_k的近似。此时，拟牛顿条件为yk=Bk+1δky_k=B_{k+1}\delta_k。
假设Bk+1=Bk+Pk+QkB_{k+1}=B_k+P_k+Q_k，此时Bk+1δk=Bkδk+Pkδk+QkδkB_{k+1}\delta_k=B_k\delta_k+P_k\delta_k+Q_k\delta_k
为了满足拟牛顿条件，则可使Pkδk=yk, Qkδk=−BkδkP_k\delta_k=y_k,\ Q_k\delta_k=-B_k\delta_k，例如取Pk=ykykTykTδk, Qk=−BkδkδkTBkδkTBkδkP_k=\frac{y_ky_k^T}{y_k^T \delta_k},\ Q_k=-\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^T B_k \delta_k}，则
Bk+1=Bk+ykykTykTδk−BkδkδkTBkδkTBkδkB_{k+1}=B_k+\frac{y_ky_k^T}{y_k^T \delta_k}-\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^T B_k \delta_k}
可以证明，若初始B0B_0正定，则对任意k都有BkB_k正定。
➣➣BFGS算法➣➣
输入：目标函数f(x)f(\mathbf{x})，梯度g(x)=∇f(x)g(\mathbf{x})=\nabla f(\mathbf{x})，精度要求ϵ\epsilon;
输出：f(x)f(\mathbf{x})的极小值点x∗\mathbf{x}^*.
（1）取初值点x0\mathbf{x}_0，取B0B_0为正定对称阵，置k=0k=0;
（2）计算gk=g(x(k))g_k=g(\mathbf{x}^{(k)})，若∣∣gk∣∣≤ϵ||g_k||\leq\epsilon，返回x∗=x(k)\mathbf{x}^*=\mathbf{x}^{(k)}；否则，转（3）
（3）求pk:Bkpk=−gkp_k:B_kp_k=-g_k;
（4）一维搜索：求λk\lambda_k使得f(x(k)+λkpk)=min⁡λf(x(k)+λpk)f(\mathbf{x}^{(k)}+\lambda_kp_k)=\min_\lambda f(\mathbf{x}^{(k)}+\lambda p_k);
（5）置x(k+1)=x(k)+λkpk\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}+\lambda_k p_k;
（6）计算gk+1=g(xk+1)g_{k+1}=g(\mathbf{x}_{k+1})，若∣∣gk+1∣∣≤ϵ||g_{k+1}||\leq\epsilon，返回x∗=x(k+1)\mathbf{x}^*=\mathbf{x}^{(k+1)}；否则，更新Bk+1=Bk+ykykTykTδk−BkδkδkTBkδkTBkδkB_{k+1}=B_k+\frac{y_ky_k^T}{y_k^T \delta_k}-\frac{B_k\delta_k\delta_k^TB_k}{\delta_k^T B_k \delta_k}；
（7）置k←k+1k\leftarrow k+1，转（3）。
 
Broyden算法
可以从BFGS算法中BkB_k的迭代公式中得到BFGS算法关于GkG_k的迭代公式。
记Gk=Bk−1,Gk+1=Bk+1−1G_k=B_k^{-1},G_{k+1}=B_{k+1}^{-1}，可以利用两次Sherman-Morrison公式，得到
BFGS算法关于的迭代公式：
BFGS算法关于GkG_k的迭代公式：Gk+1=(I−δkykTδkTyk)Gk(I−δkykTδkTyk)T+δkykTδkTykG_{k+1}=(I-\frac{\delta_ky_k^T}{\delta_k^Ty_k})G_k(I-\frac{\delta_ky_k^T}{\delta_k^Ty_k})^T+\frac{\delta_ky_k^T}{\delta_k^Ty_k}
Sherman-Morrison公式：假设A为n阶可逆矩阵，u, v为n维向量，且A+uvTA+uv^T也可逆，则
(A+uvT)−1=A−1−A−1uvTA−11+vTA−1u(A+uv^T)^{-1}=A^{-1}-\frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}
将DFP算法中的Gk+1G_{k+1}记为Gk+1DFPG^{DFP}_{k+1}，将BFGS算法中的Gk+1G_{k+1}记为Gk+1BFGSG^{BFGS}_{k+1}，它们都满足拟牛顿条件，则它们的线性组合也满足拟牛顿条件，并且正定。这样就得到了一类拟牛顿法，称为Broden类。
Broden类
Broden类Gk+1=αGk+1DFP+(1−α)Gk+1BFGSG_{k+1}=\alpha G^{DFP}_{k+1}+(1-\alpha)G^{BFGS}_{k+1}