最优化学习笔记1——关于拟牛顿法推导

从2020年3月份到现在，一年的时间里断断续续自学最优化，结合豆瓣读书、知乎等网站上的推荐，翻阅了以下书和课程：

• Stephen Boyd和Lieven Vandenberghe的《Convex Optimization》，这本书被公认为学习凸优化的必读经典著作，其中大量的篇幅讲应用，讲理论和算法的篇幅稍少一些。在B网站上有Boyd的视频课程，油管上也有，不过只有英文字幕，如果是零基础接触，直接看Boyd的视频稍显费劲，B站上也有凌青老师以《Convex Optimization》为教材的讲解，推荐作为入门观看；
• Jorge Nocedal和Stephen J. Wright的《Numerical Optimization》，这本书更多地介绍最优化算法，在算法层面和《Convex Optimization》正好形成互补；
• B站上崔雪婷老师的最优化课程（貌似是因为疫情原因录制的），综合了《Convex Optimization》、《Numerical Optimization》等教材的内容，推导非常详细，推荐观看，崔老师自己也写了一本书《最优化基础理论与方法》，书很薄，适合快速入门（貌似B站上的视频已经被删掉了）。

关于最优化的书和课程非常多，其他书和课程包括但不限于：

• Dimitri P. Bertsekas的《Nonlinear Programming》（作者比较苦口婆心，内容详尽）
• Yurii Nesterov的《Lectures On Convex Optimization》，被誉为思路很巧妙，不过目前没有细读，所以也没有体会到
• Aharon Ben-Taly和Arkadi Nemirovski的《Lectures On Modern Convex Optimization》，主要作为课程的配套讲义
• B站上方述诚老师的非线性规划课（只有讲义，没有对应的教材，方述诚老师另有关于线性锥优化方面的专著）
• David G. Luenberger和Yinyu Ye的《Linear and Nonlinear Programming》
• Mokhtar S. Bazaraa，Hanif D. Sherali和C. M. Shetty的《Nonlinear Programming Theory and Algorithms》
• 袁亚湘和孙文瑜的《最优化理论与方法》，毕竟是院士的著作，国产为数不多的关于最优化的好书

侧重于应用方面的有：

• Andreas Antoniou和Wu-Sheng Lu的《Practical Optimization: Algorithms and Engineering Applications》，侧重介绍具有工程实用性的算法
• Lorenz T. Biegler的《Nonlinear programming: Concepts, Algorithms, and Applications to Chemical Processes》，虽然书中是以化工过程的优化举例，但是也同样适用于其他动态优化（或者最优控制）问题

值得一提的是，最优化和控制学科里面的最优控制联系紧密，是最优控制中直接法的数学理论基石，除了上述Lorenz T. Biegler的著作涉及最优控制外，侧重于最优控制求解的还有：

• John T. Betts的《Practical methods for optimal control and estimation using nonlinear programming》
• S´ebastien Gros和Moritz Diehl的《Numerical Optimal Control》（与课程对应的讲义，科学上网可以看到视频，可惜没有配套字幕）。

个人觉得，上述这些著作全部看完不必要也不现实，应该先选择1至2本入门，然后根据自己的研究需要选读。

扯完前面的废话，谈谈拟牛顿法，主要参考《Numerical Optimization》中的内容，根据第3章的描述，拟牛顿法最开始主要用于克服牛顿法计算无约束优化问题中的搜索方向的不足（有无约束和拟牛顿的思想没有太大关系）。以线搜索为例（对于信赖域法类似），我们想要
$$\underset{x}{\min }f(x)\text{\hspace{1em}(1)}$$

的方式是通过下述迭代过程：
$$x_{k+1}=x_k+{\alpha }_k{p}_k\text{\hspace{1em}(2)}$$

式中：${\alpha }_k$为步长，$p_k$为搜索方向。

这里不考虑步长${\alpha }_k$的选取，下面主要考虑搜索方向$p_k$的选取。首先看牛顿法，需要求$f(x)$对应的Hessian矩阵${\nabla }^{2}f$，给出的搜索方向为
$$p_k=-{\nabla }^2{f}_k^{-1}\nabla f_k\text{\hspace{1em}(3)}$$

式中：${\nabla }^2{f}_k^{-1}$为${\nabla }^2{f}_k$的逆矩阵，考虑到${\nabla }^2{f}_k^{-1}$并不容易计算，拟牛顿法提出用下式代替牛顿法的搜索方向：
$$p_k=-B_k^{-1}\nabla f_k\text{\hspace{1em}(4)}$$

式中：$B_k$一般选为对称正定矩阵，若能够直接获得$B_k$的计算方式，那么就免去了求Hessian矩阵的繁琐工作。那么如何得到$B_k$呢？

1. 拟牛顿法中的割线方程推导

考虑到牛顿法其实可以从对目标函数的二阶近似获得，即考虑在当前解$x_k$附近对目标函数的二阶近似：
$$m_k(p)=f_k+\nabla f_k^{\mathrm{T}}p+\frac{1}{2}{p}^{\mathrm{T}}{\nabla }^2{f}_{k}p\text{\hspace{1em}(5)}$$
式中：$p=x-x_k$。选择$p$使得式(5)中的目标函数取极值，就得到了牛顿法的公式，那么一个自然的想法是，式(4)可看作是使得下述目标函数取极值：
$$m_k(p)=f_k+\nabla f_k^{\mathrm{T}}p+\frac{1}{2}{p}^{\mathrm{T}}{B}_{k}p\text{\hspace{1em}(6)}$$

即用对称正定矩阵$B_k$代替Hessian矩阵${\nabla }^2{f}_k$。由于直接获得$B_k$的表达式不太容易，转而考虑迭代计算$B_k$，即需要推导$B_{k+1}$和$B_k$的关系，根据式(6)，到了$x_{k+1}$附近时，有
$$m_{k+1}(p)=f_{k+1}+\nabla f_{k+1}^{\mathrm{T}}p+\frac{1}{2}{p}^{\mathrm{T}}{B}_{k+1}p\text{\hspace{1em}(7)}$$

式中：$p=x-x_{k+1}$。根据式(7)，若$p=0$，即停留在$x_{k+1}$时，有$\nabla m_{k+1}(0)=\nabla f_{k+1}$，于是想到后退到$x_k$(对应$p=-{\alpha }_k{p}_k$)时，$m_{k+1}$的梯度应与$\nabla f_k$相等(即$m_{k+1}$应能保证从$x_k$到$x_{k+1}$的过程中对$\nabla f$足够近似)，从而有
$$\nabla m_{k+1}(-{\alpha }_k{p}_k)=\nabla f_{k+1}-{\alpha }_k{B}_{k+1}{p}_k=\nabla f_k\text{\hspace{1em}(8)}$$
也即
$$B_{k+1}{\alpha }_k{p}_k=\nabla f_{k+1}-\nabla f_k\text{\hspace{1em}(9)}$$

定义$s_k=x_{k+1}-x_k={\alpha }_k{p}_k$, $y_k=\nabla f_{k+1}-\nabla f_k$，则式(9)可写为
$$B_{k+1}{s}_k=y_k\text{\hspace{1em}(10)}$$

式(10)又称为割线方程(secant equation)。令$H_k=B_k^{-1}$，则有互补形式的割线方程：
$$H_{k+1}{y}_k=s_k\text{\hspace{1em}(11)}$$

另一方面，式(8)也可这么理解（参考《最优化理论与方法》），即直接考虑以$x$为自变量（而不是引入额外的变量$p$）:
$$f(x)\approx f(x_{k+1})+\nabla f_{k+1}(x-x_{k+1})+\frac{1}{2}(x-x_{k+1}{)}^{\mathrm{T}}{B}_{k+1}(x-x_{k+1})\text{\hspace{1em}(12)}$$

对上式两边求导，有
$$\nabla f(x)\approx \nabla f_{k+1}+B_{k+1}(x-x_{k+1})\text{\hspace{1em}(13)}$$

上式中令$x=x_k$即可得到式(8)。

2.拟牛顿法中的矩阵递推

拟牛顿法中的矩阵递推关系具有比较强的构造痕迹，典型的方法包括秩一和秩二校正。在介绍具体内容之前，首先引入矩阵求逆的Sherman–Morrison公式：对于非奇异方阵$A$，设$\bar{A}=A+a{b}^{\mathrm{T}}$，其中$a,b\in {\mathbb{R}}^n$，若$\bar{A}$非奇异，则有
$${\bar{A}}^{-1}=A^{-1}-\frac{A^{-1}a{b}^{\mathrm{T}}{A}^{-1}}{1+b^{\mathrm{T}}{A}^{-1}a}\text{\hspace{1em}(14)}$$

2.1 秩一校正

为了保证对称性，首先考虑构造如下关系（秩一校正）：
$$B_{k+1}=B_k+\sigma v{v}^{\mathrm{T}}\text{\hspace{1em}(15)}$$

式中：$\sigma$取值为$1$或者$-1$，且选择$\sigma$和$v$使得$B_{k+1}$满足割线方程(10)，于是有
$$y_k=B_k{s}_k+\left[\sigma v^{\mathrm{T}}{s}_k\right]v\Rightarrow \left[\sigma v^{\mathrm{T}}{s}_k\right]v=y_k-B_k{s}_k\text{\hspace{1em}(16)}$$

由于$\sigma v^{\mathrm{T}}{s}_k$为标量，因此$v$必定为$y_k-B_k{s}_k$的倍数，即存在标量$\delta$使得$v=\delta \left(y_k-B_k{s}_k\right)$，则有
$$y_k-B_k{s}_k=\sigma {\delta }^2\left[s_k^{\mathrm{T}}\left(y_k-B_k{s}_k\right)\right]\left(y_k-B_k{s}_k\right)\text{\hspace{1em}(17)}$$

为确保上式成立，可选取
$$\sigma =\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left[s_k^{\mathrm{T}}\left(y_k-B_k{s}_k\right)\right],\delta =\pm {\left|s_k^{\mathrm{T}}\left(y_k-B_k{s}_k\right)\right|}^{-\frac{1}{2}}\text{\hspace{1em}(18)}$$

则$v=\pm {\left|s_k^{\mathrm{T}}\left(y_k-B_k{s}_k\right)\right|}^{-\frac{1}{2}}\left(y_k-B_k{s}_k\right)$，结合式(15)和(18)可得
$$\begin{array}{rl}{B}_{k+1}& =B_k+\sigma v{v}^{\mathrm{T}}\\ & =B_k+\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left[s_k^{\mathrm{T}}\left(y_k-B_k{s}_k\right)\right]\times {\left|s_k^{\mathrm{T}}\left(y_k-B_k{s}_k\right)\right|}^{-1}\left(y_k-B_k{s}_k\right){\left(y_k-B_k{s}_k\right)}^{\mathrm{T}}\\ & =B_k+\frac{\left(y_k-B_k{s}_k\right){\left(y_k-B_k{s}_k\right)}^{\mathrm{T}}}{{\left(y_k-B_k{s}_k\right)}^{\mathrm{T}}{s}_k}\end{array}\text{\hspace{1em}(19)}$$

根据式(19)，取$a=\frac{y_k-B_k{s}_k}{{\left(y_k-B_k{s}_k\right)}^{\mathrm{T}}{s}_k}$，$b=y_k-B_k{s}_k$，且考虑到$H_k$和$B_k$为互逆关系，且均为对称矩阵，则由Sherman–Morrison公式(14)有
$$\begin{array}{rl}{H}_{k+1}& =H_k-\frac{H_k\frac{y_k-B_k{s}_k}{{\left(y_k-B_k{s}_k\right)}^{\mathrm{T}}{s}_k}{\left(y_k-B_k{s}_k\right)}^{\mathrm{T}}{H}_k}{1+{\left(y_k-B_k{s}_k\right)}^{\mathrm{T}}{H}_k\frac{y_k-B_k{s}_k}{{\left(y_k-B_k{s}_k\right)}^{\mathrm{T}}{s}_k}}\\ & =H_k-\frac{\left(H_k{y}_k-s_k\right)(H}{{\left(y_k-B_k{s}_k\right)}^{\mathrm{T}}{s}_k+{\left(y_k-B_k{s}_k\right)}^{\mathrm{T}}\left(H_k{y}_k-s_k\right)}\\ & \end{array}$$