统计计算第二、三节课，一些优化的方法

这是我上的统计计算课讲的主要内容，写在这可以互相交流，有些地方我不是很理解的会标出来（用加粗斜体*标出），求大佬在留言处表达自己的看法，另外如果有啥问题也可以在留言处留言，如果我看到了会回复

这次的内容相比于上次更杂，也稍微有点难度，含有一些证明，因此这次的文章结构采用三个阶段的写法，第一个阶段介绍概念，了解方法的动机、主要结果、以及优缺点，第二个阶段是最基本方法的证明，第三个阶段扩展到其他方法的证明（不过我可能写不出来）

第一阶段

背景：求解最大似然估计
主要方法：拉格朗日乘子法解带约束的凸问题，梯度下降法及其一些扩展数值求解最优点

拉格朗日乘子法

自行百度的概念：凸集，凸函数（这里指下凸函数）
凸函数的重要例子：p范数（p>=1）
凸函数的重要性质：若$f(x)$是凸函数，则有$$f(E(x))\le E(f(x))（证明放在第二阶段）$$

凸优化问题的表述

$$min{f}_0(x)$$$$s.t.f_i(x)\le 0i=1,2,...,m$$$$h_j(x)=0j=1,2,...,p$$
其中所有函数都是凸函数。

记号
最优（小）值：$p^{\ast }$
可行域：$D={\bigcap }_{i=1}^n{D}_{f_i}\cap {\bigcap }_{j=1}^p{D}_{h_j}$
最优点：$x^{\ast }$，且有$f(x^{\ast })=p^{\ast }$

凸问题中判别最优点的准则
设$f_0(x)$可微，则$x^{\ast }$是最优点，当且仅当$$▽f_0(x{)}^T(y-x)\ge 0\forall 可行域中的y$$

拉格朗日函数求解最优点

拉格朗日函数
$L(x,\mathbf{\lambda },\mathbf{\nu })=f_0(x)+{\sum }_{i=1}^n{\lambda }_i{f}_i(x)+{\sum }_{j=1}^p{\nu }_j{h}_j(x)$

对偶函数
$g(\mathbf{\lambda },\mathbf{\nu })=\underset{x\in D}{\inf }L(x,\mathbf{\lambda },\mathbf{\nu })$
我们可以观察到对偶函数有如下重要性质：
1、对偶函数始终都是凸函数（不管原问题是否是凸问题）
2、若对于$\forall$可行域中的x，若$\mathbf{\lambda }\ge 0$，则有$g(\mathbf{\lambda },\mathbf{\nu })\le L(x,\mathbf{\lambda },\mathbf{\nu })\le f_0(x)$，因此可推出$g(\mathbf{\lambda },\mathbf{\nu })\le p^{\ast }，\forall \mathbf{\lambda }\ge 0,\mathbf{\nu }\in R^p$

拉格朗日对偶问题$$ming(\mathbf{\lambda },\mathbf{\nu })$$$$s.t.\mathbf{\lambda }\ge 0$$记对偶问题的最小值为$d^{\ast }$，最小值点为$({\mathbf{\lambda }}^{\ast },{\mathbf{\nu }}^{\ast })$
若$d^{\ast }\le p^{\ast }$，则称该情况为弱对偶，称$d^{\ast }-p^{\ast }$为最优值对偶间隙
若$d^{\ast }=p^{\ast }$，则称该情况为强对偶

强对偶的一个充分条件
Slater’s Condition（还不太明白这个条件说的啥，查好后补充）*

强对偶能带来的一些结论
首先我们考察下面这个连不等式
$$\begin{array}{rl}{f}_0(x^{\ast })=& g({\mathbf{\lambda }}^{\ast },{\mathbf{\nu }}^{\ast })=\underset{x\in D}{\inf }L(x,{\mathbf{\lambda }}^{\ast },{\mathbf{\nu }}^{\ast })\\ =& \underset{x\in D}{\inf }{f}_0(x)+\sum _{i=1}^n{\lambda }_i^{\ast }{f}_i(x)+\sum _{j=1}^p{\nu }_j^{\ast }{h}_j(x)\\ \le & f_0(x^{\ast })+\sum _{i=1}^n{\lambda }_i^{\ast }{f}_i(x^{\ast })+\sum _{j=1}^p{\nu }_j^{\ast }{h}_j(x^{\ast })\\ \le & f_0(x^{\ast })\end{array}$$由于有强对偶，所以上式中的不等号全部全部取等，可得两个重要条件：
(1)：$x^{\ast }$使$L(x,{\mathbf{\lambda }}^{\ast },{\mathbf{\nu }}^{\ast })$取最小值
(2)：${\lambda }_i^{\ast }{f}_i(x^{\ast })=0$ （互补条件）

通过上面的讨论我们可以得到求解优化问题的一套流程：首先写出拉格朗日函数，算出对偶函数，由对偶函数的凸性算出对偶函数的极值点$({\mathbf{\lambda }}^{\ast },{\mathbf{\nu }}^{\ast })$，最后在该极值点下求出拉格朗日函数的极值点$x^{\ast }$。
但是我们发现这一套流程很麻烦，所以想找出一个较为容易的求解方法，这就催生了KKT条件。

KKT条件
设$(x^{\ast },{\mathbf{\lambda }}^{\ast },{\mathbf{\nu }}^{\ast })$是最优点，则KKT条件的内容为拉格朗日函数关于x的梯度在最优点处取值为0+约束条件+${\mathbf{\lambda }}^{\ast }\ge 0$+互补条件
首先由上面的讨论可以知道KKT条件是最优点的必要条件
定理：当原问题是凸问题时，KKT条件是最优点的充要条件 （不会证）

梯度下降法

上面的拉格朗日乘子法对于有解析解的问题较为好办，但如果问题没有或者很难求出解析解，数值求最优点就变得很重要了，接下来将讨论无约束问题的梯度下降法及其若干变种，每种方法将从算法流程，收敛性方面来讨论

在讨论之前，我们在这里列出三个假设，便于下文引用：
（1）f在$R^n$上凸且连续可微
（2）$▽f(x)$在$R^n$上L-李普西斯连续，即对$\forall x,y\in R^n,\Vert ▽f(x)-▽f(y)\Vert \le L\Vert x-y\Vert$
（3）f的最优值有限且最优点能取到
（4）f是m-强凸的

最原始的梯度下降法

算法（略）
收敛速度
（1,2,3）无（4）时为$O(1/k)$，有是（4）时为指数收敛

牛顿法

牛顿法产生的背景
对于梯度下降法，我们可以如下理解：假设$f(x)$有一阶导数，所以可以将$f(x)$作一阶泰勒展开并略去高阶小量，得到$f(x+v)=f(x)+f^{{}^{\prime }}(x{)}^{T}v$，我们的目标是取某个方向的单位向量$v$，使得$f(x+v)$最小，此时$v$就是最原始的梯度下降法中得到的方向

牛顿法
牛顿法就是当$f(x)$有二阶导数时，将$f(x)$作二阶泰勒展开，略去高阶小量并解上述的最优化问题得到的$v$，具体形式为$$v=-({▽}^{2}f(\mathbf{x}){)}^{-1}▽f(\mathbf{x})$$

当然我们也可以从另一方面来看牛顿法，使用一个正定矩阵新定义一个距离，在该距离下，v是下降最快的方向（与最原始的梯度下降法的想法类似）

拟牛顿法

拟牛顿法产生的背景
在牛顿法中，我们需要计算海森矩阵的逆，在矩阵较大时计算开销会很大，所以我们希望每次更新$x$时使用一个矩阵近似海森矩阵，而这个近似矩阵的逆便于求得，减少计算开销

算法
设当前所在点为$x_{n+1}$，上一步所在点为$x_n$，已经求得上一步的近似矩阵$M_n$，则算法如下：$$\begin{gathered} \Delta x_n=x_{n+1}-x_n \\ {y}_n=▽f(x_{n+1})-▽f(x_n) \\ {v}_n=y_n-M_n{x}_n \\ {M}_{n+1}=M_n+\frac{v_n{v}_n^T}{v_n^T\Delta x_n} \\ \text{计算Mn+1的逆矩阵Mn+1−1} \\ {x}_{n+2}=x_{n+1}-t{M}_{n+1}^{-1}▽f(x_{n+1}) \end{gathered}$$
其中$M_{n+1}^{-1}$的计算使用了Sherman-Morrison-Woodbury 公式$$\begin{gathered} (A+UCV{)}^{-1}=A^{-1}-A^{-1}U(C^{-1}+V{A}^{-1}U{)}^{-1}V{A}^{-1} \\ 其中A\in R^{n\times n},U\in R^{n\times d},C\in R^{d\times d},V\in R^{d\times n} \\ 其中d<n,所有矩阵满秩 \end{gathered}$$（该公式暂时不在本文中证明，可能以后会证一下）
在拟牛顿法中矩阵C为1，矩阵U为列向量$v_n$，矩阵V为行向量$v_n^T$，所以$M_{n+1}^{-1}$易求出

拟牛顿法的扩展
上面在更新$M_n$时，使用了一个秩为1的矩阵，为了更加逼近海森矩阵，我们想用一个秩为2的矩阵进行更新，这就产生了BFGS方法
$$M_{n+1}=M_n+\frac{y_n{y}_n^T}{y_n^T\Delta x_n}-\frac{M_n\Delta x_n(M_n\Delta x_n{)}^T}{(\Delta x_n{)}^T{M}_n\Delta x_n}$$

遗留的问题
在这里我没有说明为什么这样构造的$M_n$会收敛到海森矩阵，以及为什么BFGS方法比秩1方法好，我会在以后填一下这个坑
关于拟牛顿法的收敛速度，我也会通过数值试验来填一下坑

动量方法

背景
原始的梯度下降法有如下的缺点：
（1）陷入鞍点
（2）收敛速度不够快
（3）难以解决大规模数据
（4）无法求没有导数的函数的最值
这就催生了另一些改进的下降方法来求函数最值

原始的动量方法

算法
$$\begin{gathered} m_k=\mu m_{k-1}+(1-\mu )▽f(x_{n-1}) \\ {x}_n=x_{n-1}-\alpha m_k \\ 0\le \mu <1 \end{gathered}$$
直观想法
从名字就可以看出，选择下一个点的原则和物体运动类似，其方向收到上一步方向的影响，当$\mu =0$时，即为梯度下降法

优点
相比于梯度下降法，在等高线为椭圆形，出发点位于长轴附近时（偏离长轴一些），动量方法的更新方向沿长轴的分量会逐渐积累，沿短轴的分量会相互抵消，这样可以更快的向最优点收敛（如图）

Nesterov加速法（NAG）

背景
见其他大佬写的文章

https://blog.csdn.net/tsyccnh/article/details/76673073

算法
$$\begin{gathered} 令x_0=x_{-1} \\ y=x_{n-1}+\frac{n-2}{n+1}(x_{n-1}-x_{n-2}) \\ {x}_n=y-t▽f(y) \end{gathered}$$可以看出此方法和动量法差别在于需要求导数的点的位置，此方法求导的位置比当前位置靠前一点，而动量方法求导的位置就是在当前位置。

优点
在一定条件下，收敛速度在达到一阶方法的最快速度（一阶方法是什么锤子，找不到定义），下面这个定理说明一阶方法收敛的最快速度为$\frac{1}{k^2}$（本文不证明该定理）

定理：对任意的整数$\mathbf{k}<\frac{\mathbf{n}-1}{2}$(n为维数)和初始值${\mathbf{x}}_0,$,都存在满足假设（1，2，3）的函数f，使得对于任意的一阶方法，$$\mathbf{f}({\mathbf{x}}_{\mathbf{k}})-\mathbf{f}({\mathbf{x}}^{\ast })\ge \frac{3}{32}\frac{\mathbf{L}\Vert {\mathbf{x}}_0-{\mathbf{x}}^{\ast }{\Vert }^2}{(\mathbf{k}+1{)}^2}$$

收敛速度

定理：f满足假设（1,2,3)，则对于$\mathbf{t}\le \frac{1}{\mathbf{L}}$，有$$\mathbf{f}({\mathbf{x}}_{\mathbf{k}})-\mathbf{f}({\mathbf{x}}^{\ast })\le \frac{2\Vert {\mathbf{x}}_0-{\mathbf{x}}^{\ast }{\Vert }^2}{\mathbf{t}(\mathbf{k}+1{)}^2}$$
该定理说明Nesterov加速法在一定条件下，收敛速度在达到一阶方法的最快速度，我们将在第二阶段给出证明

小结
至此我们可以解决梯度下降法的缺点（2）

近端（proximal）梯度下降

背景
该方法可以求部分没有导数的函数的最值，只要这个这个函数可以分解为两个部分$f=g+h$，$g$为凸且连续可微的，$h$为凸且较为简单，这个较为简单定义得比较模糊，可以通过下面近端函数的一些例子来理解。

prox函数
定义：设h是凸函数，则关于h的近端函数为$$pro{x}_h(x)=\underset{u}{\mathrm{argmin}}(h(u)+\frac{1}{2}\Vert u-x{\Vert }_2^2)$$例如：$$\begin{gathered} （1）h(x)=0时，pro{x}_h(x)=x \\ （2）定义示性函数1_C(x)=\{\begin{array}{rl}& 0x\in C\\ & +\infty x\notin C\end{array} \\ 当h(x)=1_C(x)时，pro{x}_h(x)为点x在C上的投影 \\ (3)当h(x)=\Vert x{\Vert }_1时，pro{x}_h为软阈值算子， \\ 其中pro{x}_h(x)的第i个分量为pro{x}_h(x{)}_i=\{\begin{array}{rl}& x_i-1x_i>1\\ & 0-1\le x_i\le 1\\ & x_i+1x<-1\end{array} \end{gathered}$$通过上面的例子，我们可以找到一些在近端算子作用下的形式较为“简单”的函数。

算法
$$x_{n+1}=pro{x}_{t\cdot h}(x_n-t▽g(x_n))$$

计算上的解释
$$\begin{array}{rl}{x}_{n+1}=& \underset{u}{\mathrm{argmin}}(h(u)+\frac{1}{2t}\Vert u-x_n+t▽g(x_n){\Vert }_2^2)\\ =& \underset{u}{\mathrm{argmin}}(h(u)+g(x)+▽g(x_n{)}^T(u-x_n)+\frac{1}{2t}\Vert u-x_n{\Vert }_2^2)\end{array}$$可以看出$x_{n+1}$最小化了$h(u)$加$g(x)$的一个二次近似（不是泰勒展开，但接近）

直观上的解释
见其他大佬写的文章

https://www.cnblogs.com/raby/p/5886699.html

具体的近端梯度下降的例子
当$h(x)=0$时，即为梯度下降
当$h(x)=1_C(x)$时，即为投影梯度下降
当$h(x)=\Vert x{\Vert }_1$时，即为ISTA算法（Iterative Shrinkage-Thresholding Algorithm）$$\begin{gathered} x_{n+1}=S_t(x-t▽g(x_n)) \\ 其中S_t(u)=(\Vert u\Vert -t{)}_{+}sign(u) \end{gathered}$$

近端梯度下降法的收敛速度和梯度下降法相同，只要$\mathbf{h}(\mathbf{x})$是凸且闭的（本文不给出证明）

加速近端梯度下降法

背景
模仿NAG法

算法
$$\begin{gathered} 令x_0=x_{-1} \\ y=x_{k-1}+\frac{k-2}{k+1}(x_{k-1}-x_{k-2}) \\ {x}_k=pro{x}_{th}(y-t▽g(y)) \end{gathered}$$当$h(x)=\Vert x{\Vert }_1$时，算法又称为FISTA

收敛速度和NAG法一样（本文不给出证明）

小结
至此我们可以部分解决梯度下降法的缺点（4）

随机优化

该部分以浏览为主，理论证明部分较少，算法流程较多，会提及部分算法的想法和优缺点，可快速读过（主要是敲字敲到这太累了，不想敲了$\times$）
我参考了其他大佬的文章

https://blog.csdn.net/u014595019/article/details/52989301

随机梯度下降（SGD）

目的
逃离不稳定的极值点，摆脱鞍点，可以处理大规模数据

算法
目标优化函数仍然是$f(x)$
$$\begin{gathered} x_{k+1}=x_k-tg(x_k) \\ 其中E(g(x))=▽f(x) \end{gathered}$$

例子
太麻烦，不想敲，如果有想要的可以在评论里留言，我会填坑

收敛速度

缺点
收敛速度慢，调参（学习率）较麻烦

自适应随机梯度下降（Adaptive Stochastic Gradient Descent）

AdaGrad算法

目的
解决SGD学习率不变的缺点

想法来源
假设$f(x)$是某个随机变量的期望，那么我们可以使用Fisher信息量来调整学习率（因为Fisher信息量可以类比于函数的曲率），这启发我们使用$g({\theta }_k)g({\theta }_k{)}^T$的某种平均来调整学习率

算法
$$\begin{gathered} g({\theta }_k)仍表示随机取样本点算出的梯度 \\ r=r+g({\theta }_k)\odot g({\theta }_k)其中圈乘表示逐元素相乘 \\ {\theta }_{k+1}={\theta }_k-\frac{\eta }{ϵ+\sqrt{r}}\odot g({\theta }_k)其中所有运算均为逐元素运算 \end{gathered}$$其中超参数有$\eta ,ϵ$

缺点
$r$的累积最终导致学习率下降

RMSprop

目的
解决AdaGrad的缺点

想法
降低$r$的累积速度

算法
$$\begin{gathered} g({\theta }_k)仍表示随机取样本点算出的梯度 \\ r=\rho r+(1-\rho )g({\theta }_k)\odot g({\theta }_k)其中圈乘表示逐元素相乘 \\ {\theta }_{k+1}={\theta }_k-\frac{\eta }{ϵ+\sqrt{r}}\odot g({\theta }_k)其中所有运算均为逐元素运算 \end{gathered}$$其中超参数有$\eta ,ϵ,\rho$
可以看出此方法中有些矩方法的样子

Adam

想法
融合矩方法的优点，将$r$和梯度都使用矩方法更新

算法
不想敲了，上面的大佬链接写得很清楚

小结
至此，我们解决了梯度下降法的缺点（1）（3）.

总之，这些算法长什么形式不重要，因为一般从形式上也看不出大佬是怎么想出来的，尤其是Adam方法和NAG方法，所以阅读原始文章和其中的证明对于理解方法更重要，可以先在这里留个坑，以后用到时会挖

第二阶段

第一阶段终于写完了，感觉长得都不像第一阶段$\times$)，第二阶段就两个证明，第一个是梯度下降法的收敛速度，第二个证明是随机梯度下降法的收敛速度

第一个证明