压缩感知与稀疏模型——Proximal Gradient
实际上这节课讲的内容也比较多，对应的是convex optimization，这里主要介绍一下lasso回归以及对应使用的Proximal Gradient Decent（PGD）。
Lasso Regression
之前提到，对于$l^0$范数问题，用来求解最稀疏的解，我们可以将最小化$l^0$范数放松到最小化$l^1$范数。假如现在有下面的问题：
\text{minimize }\Vert x \Vert_1\\ \text{subject to } Ax = y.
这个问题是著名的BP(Basic Pursuit)问题，它的另外一个版本是考虑了高斯噪声：$y = Ax_0 + z$，将这个问题的限制用一个惩罚函数$f(x) = \frac {1}{2} \Vert y-Ax \Vert_2^2$代替，可以得到下面的解决策略：
\min_x \frac{1}{2} \Vert y - Ax \Vert ^2_2 + \lambda \Vert x\Vert_1
这就是著名的Lasso回归。
之前学习的优化方法大多数是对于$l^2$范数的优化，因为它处处可导。而$l^1$范数虽然也是凸函数，但是却不是平滑函数，在局部会出现不可导的情况。
除此之外，有着类似形式的问题还有低秩矩阵的恢复问题。我们希望恢复低秩矩阵$L_0$，能观察到的量为$Y = L_0 + S_0$。一个常用的方法是解决PCP（principal component pursuit）问题：
\begin{array}{cl}{ {\operatorname{minimize} } }& {\|\boldsymbol{L}\|_ {*}+\lambda\|\boldsymbol{S}\|_{1} } \\ {\text { subject to } } & {\boldsymbol{L}+\boldsymbol{S}=\boldsymbol{Y} } \end{array}
如果数据中包含噪声，那么一个更稳定的版本为:
\operatorname{minimize} \frac{1}{2}\|\boldsymbol{Y}-\boldsymbol{L}-\boldsymbol{S}\|_{F}^{2}+\lambda_{1}\|\boldsymbol{L}\|_{*}+\lambda_{2}\|\boldsymbol{S}\|_{1}
Proximal Gradient Methods
实际上，我们会遇到很多类似于lasso回归一样的优化问题，它们都有着下面的形式：
F(x) = f(x) + g(x)
其中$f(x)$是一个平滑的凸函数，而$g(x)$是一个凸函数但是却不是平滑的。在上面提到的Lasso问题中：
f(x) = \frac{1}{2} \Vert y - Ax \Vert ^2_2, g(x) = \lambda \Vert x\Vert_1.
我们希望能找到一种快速的方法来解决这样的优化问题。
首先，由于局部不可导，因此我们没法使用梯度下降，最容易想到的是使用subgradient（次梯度）。
使用次梯度下降：
x_{k+1} = x_k - \gamma_k g_k, g_k \in \partial F(x_k)
实际上是可以解决Lasso问题的，但是主要缺点是收敛得太慢了。一般来说，次梯度的收敛速率对于非平滑的目标函数是：
F(x_k) - F(x_*) = O(\frac{1}{\sqrt k})
convergence rates of first order methods
这部分内容介绍一下一般梯度下降（使用一阶导数）的收敛速度是$O(\frac{1}{k})$。对于convex函数来说，一阶导数估计实际上提供的是一个下界。
我们知道梯度下降每次迭代规律为：
x_{k+1} = x_k - \gamma_k \nabla f(x_k)
而迭代后的代入$f$真实值，实际上是大于根据梯度估计的值的：
f(x’) \ge f(x) + \langle\nabla f(x),x’-x\rangle
梯度越大，步子就可以迈得更大一点，因为对于平滑函数来说，这意味着距离最低点还更远一点。引入新的概念：L-Lipschitz gradient，如果一个函数满足：
\Vert \nabla f(x_1) - \nabla f(x_2)\Vert \leq L\Vert x_1 - x_2\Vert, \forall x_1,x_2 \in \mathbb R^n
$L>0$,那么这个$L$被称为Lipschitz常量。有了这个条件，我们可以给线性下界补充一个二次上界。
引理：
如果$f$是可导函数，而且$\nabla f$满足L-Lipschitz，那么对于任意$x,x’$有
f(x’) \leq f(x) + \langle\nabla f(x),x’-x\rangle + \frac{L}{2}\Vert x’ - x \Vert _2^2
证明上面的结论：
\begin{aligned} f\left(\boldsymbol{x}^{\prime}\right) &=f\left.\left(\boldsymbol{x}+t\left(\boldsymbol{x}^{\prime}-\boldsymbol{x}\right)\right)\right|_{t=1} \\ &=f(\boldsymbol{x})+\int_{t=0}^{1} \frac{d}{d t} f\left(\boldsymbol{x}+t\left(\boldsymbol{x}^{\prime}-\boldsymbol{x}\right)\right) d t \\ &=f(\boldsymbol{x})+\int_{t=0}^{1}\left\langle\nabla f\left(\boldsymbol{x}+t\left(\boldsymbol{x}^{\prime}-\boldsymbol{x}\right)\right), \boldsymbol{x}^{\prime}-\boldsymbol{x}\right\rangle d t \\ &=f(\boldsymbol{x})+\left\langle\nabla f(\boldsymbol{x}), \boldsymbol{x}^{\prime}-\boldsymbol{x}\right\rangle \\ & \quad+\int_{t=0}^{1}\left\langle\nabla f\left(\boldsymbol{x}+t\left(\boldsymbol{x}^{\prime}-\boldsymbol{x}\right)\right)-\nabla f(\boldsymbol{x}), \boldsymbol{x}^{\prime}-\boldsymbol{x}\right\rangle d t \\ & \leq f(\boldsymbol{x})+\left\langle\nabla f(\boldsymbol{x}), \boldsymbol{x}^{\prime}-\boldsymbol{x}\right\rangle+\frac{L}{2}\left\|\boldsymbol{x}^{\prime}-\boldsymbol{x}\right\|_{2}^{2} \end{aligned}
所以我们可以得到：
f\left(\boldsymbol{x}^{\prime}\right) \leq \hat{f}\left(\boldsymbol{x}^{\prime}, \boldsymbol{x}\right) \doteq f(\boldsymbol{x})+\left\langle\nabla f(\boldsymbol{x}), \boldsymbol{x}^{\prime}-\boldsymbol{x}\right\rangle+\frac{L}{2}\left\|\boldsymbol{x}^{\prime}-\boldsymbol{x}\right\|_{2}^{2}
如果我们想要最小化这个上界，可以得到：
\arg\min _{x’} \hat f(x’,x) = x - \frac{1}{L} \nabla f(x).
这就是一个很简单的梯度下降迭代。而且可以保证，这样的迭代不会增加目标函数的值，因为$\hat f(x,x) = f(x)$：
f\left(\boldsymbol{x}_{\star}^{\prime}\right) \leq \hat{f}\left(\boldsymbol{x}_{\star}^{\prime}, \boldsymbol{x}\right) \leq \hat{f}(\boldsymbol{x}, \boldsymbol{x})=f(\boldsymbol{x})
可以证明的是，使用这种迭代方法到最佳值的收敛速度是：
f\left(\boldsymbol{x}_{k}\right)-f\left(\boldsymbol{x}_{\star}\right) \leq \frac{L\left\|\boldsymbol{x}_{0}-\boldsymbol{x}_{\star}\right\|_{2}^{2} }{2 k}=O(1 / k)
这个收敛速度会比次梯度要快很多。
From gradient to proximal gradient
实际上，我们可以从上面的梯度下降算法来获得一点启发，来设计一个算法优化lasso回归这类的问题：
F(x) = f(x) + g(x).
因为原函数是不可导的，因此梯度下降算法不能直接使用。尽管如此，如果目标函数中平滑的那一项$f(x)$的导数$\nabla f$是Lipschitz，我们依然可以构造出一个$F$的简单上界：对平滑项$f$放大到上界，对于不可导项就什么也不做：
\hat{F}\left(\boldsymbol{x}^{\prime}, \boldsymbol{x}\right)=f(\boldsymbol{x})+\left\langle\nabla f(\boldsymbol{x}), \boldsymbol{x}^{\prime}-\boldsymbol{x}\right\rangle+\frac{L}{2}\left\|\boldsymbol{x}^{\prime}-\boldsymbol{x}\right\|_{2}^{2}+g\left(\boldsymbol{x}^{\prime}\right)
不断地最小化上式，可以得到类似于之前的梯度下降算法，这样的算法比次梯度下降有更好的收敛率。因此，我们要做的就是最小化$\hat F$:
\boldsymbol{x}_{k+1}=\arg \min _{\boldsymbol{x} } \hat{F}\left(\boldsymbol{x}, \boldsymbol{x}_{k}\right)
将$x_k$带入后补充全式，我们可以得到：
\hat{F}\left(\boldsymbol{x}, \boldsymbol{x}_{k}\right)=\frac{L}{2}\left\|\boldsymbol{x}-\left(\boldsymbol{x}_{k}-\frac{1}{L} \nabla f\left(\boldsymbol{x}_{k}\right)\right)\right\|_{2}^{2}+g(\boldsymbol{x})+h\left(\boldsymbol{x}_{k}\right)
这里的$h(x_k)$只与$x_k$有关，因此每次迭代就变成了：
\begin{aligned} \boldsymbol{x}_{k+1} &=\arg \min _{\boldsymbol{x} } \frac{L}{2}\left\|\boldsymbol{x}-\left(\boldsymbol{x}_{k}-\frac{1}{L} \nabla f\left(\boldsymbol{x}_{k}\right)\right)\right\|_{2}^{2}+g(\boldsymbol{x}) \\ &=\arg \min _{\boldsymbol{x} } g(\boldsymbol{x})+\frac{L}{2}\left\|\boldsymbol{x}-\boldsymbol{w}_{k}\right\|_{2}^{2} \end{aligned}
这里的$\boldsymbol{w}_{k}=\boldsymbol{x}_{k}-(1 / L) \nabla f\left(\boldsymbol{x}_{k}\right)$。也就是，我们希望让$g$尽量小，同时$x$又不要离$w_k$太远。由于二范数是凸函数，因此它总会有一个唯一最优解。
分析到这里，这个式子还是一个$F(x) = f(x) + g(x)$的形式，但是这里的$f(x)$是二次方（quadratic）的。这样形式在凸优化中经常会遇到，它有个特殊的名字——近端（Proximal Operator）：
\operatorname{prox}_{g}(\boldsymbol{w}) \doteq \arg \min _{\boldsymbol{x} }\left\{g(\boldsymbol{x})+\frac{1}{2}\|\boldsymbol{x}-\boldsymbol{w}\|_{2}^{2}\right\}
因此迭代步骤可以写为：
\boldsymbol{x}_{k+1}=\operatorname{prox}_{g / L}\left(\boldsymbol{w}_{k}\right).
比较幸运的是，一般我们遇到的凸函数，他们的近端都是有着固定的形式,或者可以很高效地被计算出来，下面将会举几个例子：
$g(x) = I_D$，指示函数，也就是对于一个凸集$D$,如果$x\in D$，$I_D(x) = 0$，否则$I_D(x) = \infty$。那么: \operatorname{prox}_{g}(\boldsymbol{w})=\arg \min _{\boldsymbol{x} \in \mathcal{D} }\|\boldsymbol{x}-\boldsymbol{w}\|_{2}^{2}=\mathcal{P}_{\mathcal{D} }[\boldsymbol{w}]
$g(\boldsymbol{x})=\lambda|\boldsymbol{x}|_{1}$，$l^1$范数：
\operatorname{prox}_{g}\left(u_{i}\right)=\operatorname{soft}\left(u_{i}, \lambda\right) \doteq \operatorname{sign}\left(u_{i}\right) \max \left(\left|u_{i}\right|-\lambda, 0\right)
$g(\boldsymbol{x})=\lambda|\boldsymbol{X}|_{*}$，矩阵核范数，那么：
\operatorname{prox}_{g}(\boldsymbol{W})=\mathcal{S}(\boldsymbol{W}, \lambda) \doteq \boldsymbol{U} \operatorname{soft}(\boldsymbol{\Sigma}, \lambda) \boldsymbol{V}^{*}
$(\boldsymbol{U}, \boldsymbol{\Sigma}, \boldsymbol{V})$是$\mathbf W$的SVD分解。
由以上这些，我们可以得到一个收敛速率为$O(\frac{1}{K})$的Proximal Gradient算法。更正式的定理如下：
$F(x) = f(x) + g(x)$，其中$f$是一个可导凸函数，并且导数符合Lipschitz，$g$是一个凸函数，考虑下面的迭代方式：
\boldsymbol{w}_{k} \leftarrow \boldsymbol{x}_{k}-\frac{1}{L} \nabla f\left(\boldsymbol{x}_{k}\right), \quad \boldsymbol{x}_{k+1} \leftarrow \operatorname{prox}_{g / L}\left(\boldsymbol{w}_{k}\right).
假设$F(x)$在$x^*$有最小值，那么对于任何$k>1$，都有：
F\left(\boldsymbol{x}_{k}\right)-F\left(\boldsymbol{x}_{\star}\right) \leq \frac{L\left\|\boldsymbol{x}_{0}-\boldsymbol{x}_{\star}\right\|_{2}^{2} }{2 k}
因此，我们得到近端梯度下降：
对于不同的目标函数，因为有着不同的近端形式，它们的近端梯度算法在具体形式上也略微有所不同，因此下面介绍几种针对不同问题的具体的近端梯度下降算法。
- Proximal Gradient for the Lasso
Lasso的$g$是$l^1$范数形式，而它的$L$可以是矩阵$A^TA$的最大特征值，s 可以提前计算得到的。有时候对于Lasso的PGD算法又被称为ISTA(iterative soft-thresholding algorithm)，具体算法流程如下图：
- Proximal Gradient for Stable PCP
观察之前不同形式的近端操作，可以看到对核范数$\Vert X\Vert_*$也是有y一个简单的近端形式的。因此我们可以用PGD来解决恢复低秩矩阵的问题，比如稳定主成分追踪（stable principal
component pursuit,PCP）: \min _{L, S}\|\boldsymbol{L}\|_{*}+\lambda\|\boldsymbol{S}\|_{1}+\frac{\mu}{2}\|\boldsymbol{Y}-\boldsymbol{L}-\boldsymbol{S}\|_{F}^{2}上式中有两个不可导项： \|\boldsymbol{L}\|_{*}+\lambda\|\boldsymbol{S}\|_{1}各自都有比较简单的近端形式。我们还需要知道下面的内容，假设$\boldsymbol{x}=\left[\boldsymbol{x}_{1} ; \boldsymbol{x}_{2}\right] \text {, } g(\boldsymbol{x})=g_{1}\left(\boldsymbol{x}_{1}\right)+g_{2}\left(\boldsymbol{x}_{2}\right)$，那么: \operatorname{prox}_{g}(\boldsymbol{w})=\left(\operatorname{prox}_{g_{1} }\left(\boldsymbol{w}_{1}\right), \operatorname{prox}_{g_{2} }\left(\boldsymbol{w}_{2}\right)\right)对于该问题的具体优化算法如下图：
Accelerate Proximal Gradient(APG)
之前的PG算法可以让收敛速度达到$O(\frac{1}{k})$，但是如果我们需要解决的问题有更特殊的结构，可以生成更高效快速的梯度算法。在1970～1980年这段日子里，有很多俄国的优化理论学者，再思考一个问题：对于优化一个平滑函数$f$，梯度方法是不是一阶方法里最佳的优化方法（这里局限于一阶方法，二阶方法有更快的收敛速度，但是我们也知道，二阶方法需要求海森矩阵，而高维情况下对海森矩阵的求解需要很大的代价）？
为了回答这个问题，首先需要一个计算模型。他们提出一个黑盒模型：模型中，算法会产生一个迭代序列$x_0,….,x_k$。对应于每个迭代，又会有对应的值$f(x_i)$以及梯度$\nabla f(x_i)$。通过这些来产生下一个迭代值，也就是：
\boldsymbol{x}_{k+1}=\varphi_{k}\left(\boldsymbol{x}_{0}, \ldots, \boldsymbol{x}_{k}, f\left(\boldsymbol{x}_{0}\right), \ldots, f\left(\boldsymbol{x}_{k}\right), \nabla f\left(\boldsymbol{x}_{0}\right), \ldots, \nabla f\left(\boldsymbol{x}_{k}\right)\right)
从上述模型中，我们可以从最坏的角度来研究。首先，固定一类函数$\mathcal F$，然后观察这个算法在最差的函数下能做到什么样子（这里更像是最好的情况下，sup上界，inf为下界）：
\sup _{f \in \mathcal{F}, \boldsymbol{x}_{0} }\left\{f\left(\boldsymbol{x}_{k}\right)-\inf _{\boldsymbol{x} } f(\boldsymbol{x})\right\}
对于我们感兴趣的一类函数，也就是可导凸函数并且满足Lipschitz条件:
\mathcal{F}_{L, r} \doteq\left\{f : \mathfrak{B}(0, r) \rightarrow \mathbb{R} |\left\|\nabla f(\boldsymbol{x})-\nabla f\left(\boldsymbol{x}^{\prime}\right)\right\| \leq L\left\|\boldsymbol{x}-\boldsymbol{x}^{\prime}\right\| \forall \boldsymbol{x}, \boldsymbol{x}^{\prime} \in \mathfrak{B}(0, r)\right\}
仅仅使用梯度下降也就是$O(\frac 1 k)$的收敛速度，我们有：
\sup _{f \in \mathcal{F}_{L, r}, \boldsymbol{x}_{0} }\left\{f\left(\boldsymbol{x}_{k}\right)-\inf _{\boldsymbol{x} } f(\boldsymbol{x})\right\} \leq \frac{C L r^{2} }{k}
而对于这类函数，可以证明的最好的下界是：
\sup _{f \in \mathcal{F}_{L, r}, x_{0} }\left\{f\left(\boldsymbol{x}_{k}\right)-\inf _{\boldsymbol{x} } f(\boldsymbol{x})\right\} \geq \frac{c L r^{2} }{k^{2} }
可以看到这之间还是有一个gap的。因此更快的收敛速度还是有可能做到的。
首先考虑比较简单的迭代方式，依然使用梯度下降：
x_{k+1} = x_k - \alpha \nabla f(x_k),
这里的$\alpha$是学习率，对于梯度符合Lipschitz条件的函数，取$\alpha = \frac{1}{L}$是一个很好的选择，而梯度$\nabla f(x_k)$是下降速度最快的方向。但是朝着当前最快的方向走并不一定就是最优的，因为梯度是局部的，因此梯度下降实际上是贪心算法，但是贪心算法不一定是最优的，另外一种保守的做法是参考之前迭代的方向，保留之前迭代的动量（momentum），被称为重量球方法：
\boldsymbol{x}_{k+1}=\boldsymbol{x}_{k}-\alpha \nabla f\left(\boldsymbol{x}_{k}\right)+\beta\left(\boldsymbol{x}_{k}-\boldsymbol{x}_{k-1}\right)
重量球方法往往会带来更平滑的迭代策略以及更快的收敛。
1983年，Yuri Nesterov提出了达到收敛速度为$O(\frac{1}{k^2})$的算法。在算法中，他也使用到了动量这个步骤。它引入了一个新的辅助点（auxiliary point）$p_{k+1}$:
\boldsymbol{p}_{k+1} \doteq \boldsymbol{x}_{k}+\beta_{k}\left(\boldsymbol{x}_{k}-\boldsymbol{x}_{k-1}\right)
而每一次迭代步骤就是：
\boldsymbol{x}_{k+1}=\boldsymbol{p}_{k+1}-\alpha \nabla f\left(\boldsymbol{p}_{k+1}\right)
这里$\alpha =\frac{1}{L}$，而$\beta _k$的值要仔细选择，来达到最佳的收敛率：
t_{k}=\frac{1+\sqrt{1+4 t_{k-1}^{2} } }{2}, \quad \beta_{k}=\frac{t_{k-1}-1}{t_{k} }
使用这种迭代策略，可以达到$O(\frac{1}{k^2})$的收敛率。
到目前位置，我们说的都是建立在平滑函数上的，而对于有非平滑函数项的目标函数，还需要做别的改进，也就是APG算法。
之前的PG算法，我们可以得到目标函数的一个上界：
\begin{aligned} \hat{F}\left(\boldsymbol{x}, \boldsymbol{x}_{k}\right) & \doteq f\left(\boldsymbol{x}_{k}\right)+\left\langle\nabla f\left(\boldsymbol{x}_{k}\right), \boldsymbol{x}-\boldsymbol{x}_{k}\right\rangle+\frac{L}{2}\left\|\boldsymbol{x}-\boldsymbol{x}_{k}\right\|_{2}^{2}+g(\boldsymbol{x}) \\ & \doteq \frac{L}{2}\left\|\boldsymbol{x}-\boldsymbol{w}_{k}\right\|_{2}^{2}+g(\boldsymbol{x})+\text { terms that do not depend on } \boldsymbol{x}, \end{aligned}
每一次迭代，实际上是最小化这个上界。
而对于Nesterov momentum来说，也可以按照类似的方法来得到它对应的迭代策略：
\boldsymbol{x}_{k+1}=\operatorname{prox}_{g}\left(\boldsymbol{p}_{k+1}-\frac{1}{L} \nabla f\left(\boldsymbol{p}_{k+1}\right)\right).
这样我们就得到了APG算法的迭代策略：
APG算法有$O(\frac{1}{k^2})$的收敛速率，比普通的梯度方法会好很多，更正式的定理为：
序列${x_k}$是通过APG策略对凸函数$F(x) = f(x) + g(x)$进行优化的迭代步骤，而且$f$的梯度是Lipschitz的，如果$F(x)$有最小值对应的$x_*$，那么对于任意的$k>1$:
F\left(\boldsymbol{x}_{k}\right)-F\left(\boldsymbol{x}_{\star}\right) \leq \frac{2 L\left\|\boldsymbol{x}_{0}-\boldsymbol{x}_{\star}\right\|_{2}^{2} }{(k+1)^{2} }
APG for Basis Pursuit Denoising
APG for Stable Principal Component Pursuit
对于收敛性的证明，这里都不多赘述了。