压缩感知还原算法

压缩感知(Compressive Sensing, CS)是一种新兴的取样方法，相比长久以来被使用的奈奎斯特取样定理(Nyquist Sampling Theorem)，能更高效的方式采样信号。压缩感知最主要利用信号的稀疏性来寻找欠定线性系统的稀疏解，因此能从较少的取样样本中还原信号。近几年有许多文献提出了许多有效的算法，大多是透过迭代的方式找到系数并还原信号，相较同时找到最大系数的方式，能更精确的还原信号。

然而，压缩感知的运作机制是建立于可同时获得取样后的样本振幅与相位资讯的前提下，方能完美重建原始稀疏讯号，因此在一些现实的应用情境中，例如X射线晶体学、全像摄影、量子光学、毫米波无线通讯…等，若其布署的感测元件或系统无法取得样本的相位资 讯时(或相位资讯不可靠)，则会严重影响压缩感知讯号重建的效果。为了克服此难题，近几年许多研究致力于将建立在线性测量模型(Linear measurement model)下的压缩感知理论与相关的讯号重建技术开发延伸至非线性测量系统模型；特别是针对只知道线性量测资料的振幅资讯这类的非线性模型进行深入探讨。由于是在不知道量测值的相位资讯条件下完成稀疏讯号还原，因此这类的研究课题被称为压缩相位回复(Compressive phase retrieval)。许多高能源效益的稀疏讯号还原方法相继被提出，例如基于凸优化理论所开发的 PhaseLift与PhaseCut还原方法，或是GESPAR以及SPARTA等贪婪式算法(Greedy algorithms)。

还原算法

对于压缩感知中的稀疏信号还原有很多不同类型的算法，这些算法大致可以被分为六类，以下将详细叙述各类型算法的特色^[1]。

凸松弛(Convex Relaxation)

这类算法最主要透过线性规划方式求解凸优化的问题来对信号进行还原。 虽然要精确地还原信号所需的测量数量很少，但是这些方法在计算上是较复杂的。主要代表的算法有基追踪(Basis Pursuit, BP)、抗噪基追踪(Basis Pursuit De-noising, BPDN)、最小角度回归(Least Angle Regression, LARS)。

基追踪优化公式：${\displaystyle \min _{x}\|x\|_{1}\quad {\mbox{subject to}}\quad y=Ax.}$

抗噪基追踪优化公式：${\displaystyle \min _{x}{\frac {1}{2}}\|y-Ax\|_{2}^{2}+\lambda \|x\|_{1},}$

贪婪迭代算法(Greedy Iterative Algorithm)

这类算法主要是透过迭代的方式逐步找到答案来还原信号。这类型的算法在每次迭代中，以贪婪的方式选择传感矩阵${\displaystyle \Theta }$中与${\displaystyle Y}$做内积最大的列，也就是选出与传感矩阵中最相关的列，接着该列对${\displaystyle Y}$的贡献将会从${\displaystyle Y}$中被减去，并且对相减之后的结果进行迭代，直到识别出正确的列集合。相比之下，最小平方误差的方法则是在每次迭代中最小化误差。这类算法的优点在于运算复杂度低，还原速度快，常见的算法有正交匹配追踪(Orthogonal Matching Pursuit, OMP)、随机梯度追踪(Stochastic Gradient Pursuit,SGP)、正则化正交匹配追踪(Regularized OMP)、压缩采样匹配追踪(Compressive Sampling Matching Pursuit, CoSaMP)和子空间追踪(Subspace Pursuit,SP)。

正交匹配追踪(Orthogonal Matching Pursuit, OMP)

• 
  正交匹配追踪(OMP)算法

一开始会有一个空集合${\displaystyle \Lambda _{0}}$，以及剩余的部分${\displaystyle r_{0}}$，每次迭代都会从${\displaystyle {\overline {\Lambda {_{l-1}}}}}$找出一个A矩阵投影到剩余${\displaystyle r_{l-1}}$有最大值的位置，把这个位置加到${\displaystyle \Lambda _{l}}$之中，并从${\displaystyle {\overline {\Lambda {_{l}}}}}$当中移除，最后再更新剩余${\displaystyle r_{l}}$。利用迭代的方式，不断地找出x向量当中最有可能非零的位置，直到达到算法停止条件。

随机梯度追踪(Stochastic Gradient Pursuit,SGP)

可视为是正交匹配追踪的改良版，在面对噪声的环境下也能更好地还原。整体算法都与正交匹配追踪相同，除了用最小二乘法(least square,LS)那一项。

最小二乘法可以用零强制线性均衡(zero forcing linear equalization,ZFLE)来实现：

• ${\displaystyle x_{ZF}=A_{\Omega }^{\dagger }y=A_{\Omega }^{\dagger }(Ax+n)=x+A_{\Omega }^{T}(A_{\Omega }A_{\Omega }^{T})^{-1}n}$,n为噪声

但是零强制线性均衡的缺点可以直接从上式看出，会造成放大噪声的影响。为了解决噪声的问题，可以利用最小均方误差(minimum mean square error,MMSE)均衡器来最小化整体的误差

• ${\displaystyle x_{MMSE}=A_{\Omega }^{T}(A_{\Omega }A_{\Omega }^{T}+{\frac {\sigma _{n}^{2}}{\sigma _{x|\Omega }^{2}}}I)^{-1}y}$

和零强制线性均衡相比，在反矩阵内部多增加了额外的项，可直接视为噪声项目。虽然最小均方误差方法的正交匹配追踪^[2]更能够对抗噪声干扰，但是${\displaystyle {\sigma _{x|\Omega }^{2}}}$和${\displaystyle {\sigma _{n}^{2}}}$并不适用压缩感知的众多应用上，所以并不能直接套用在压缩感知的还原算法上。

因为最小均方(least mean square,LMS)在不断地迭代后的平衡值会接近最小均方误差，因此随机梯度追踪利用最小均方来取代最小二乘法。最小均方的迭代将根据预估误差来调整support系数，而预估误差可表示为：

• ${\displaystyle e_{\lambda }=y_{\lambda }-{\tilde {A_{\lambda }}}{\tilde {x_{\lambda }}}}$

因此最小均方的梯度下降递归(gradient decent recursion)可以表示为：

• ${\displaystyle {\tilde {x}}_{\lambda +1}={\tilde {x}}_{\lambda }+\mu e_{\lambda }A_{\Omega }^{T}}$，${\displaystyle \mu }$为step size

总结来说随机梯度追踪算法里面包含了两个迭代循环。一个大的迭代是由正交匹配追踪所形成，而小的循环则被大的迭代所包含，是由最小均方所形成，不断地迭代来更新${\displaystyle {\tilde {x}}_{\lambda }}$的值。

压缩采样匹配追踪(Compressive Sampling Matching Pursuit, CoSaMP)

在进行该还原算法前需要额外的输入参数-稀疏性数量K(信号非零项的数量)。和一般的贪婪迭代算法一样，整个还原过程同样是由(1)匹配追踪(找出内积最大的列)、(2)还原值估计(通过最大内积列来估计还原值)、(3)剩余值(residual)更新(通过估计还原值来更新误差)不停地迭代组成。

和正交匹配追踪不同的是正交匹配追踪在每次迭代的匹配追踪上只寻找并增加一个内积最大的列；而压缩采样匹配追踪则在每次迭代的匹配追踪上寻找并增加2K个内积最大的列。

压缩采样匹配追踪过程：

参数定义：π为相关向量、r(t)为在第t次迭代的剩余值(residual)、ω为在π里面最大的值、Ω为support的集合

1. 初始化，${\displaystyle x(0)=0_{Px1},r(0)=y,\Omega (0)=\emptyset ,t=0}$
2. 计算相关向量π
   • π = |${\displaystyle A^{T}}$r(t)|
3. 从相关向量π里面寻找出2K个最高correlation值的supports
   • ω ${\displaystyle \equiv }$${\displaystyle {\underset {\omega \in [P]}{\arg max}}}$x{ | π(ω) | }
4. 把找到的support加入support 集合里面
   • Ω = Ω ᑌ supp(T(π,2K))，T(π,2K)为在π里面2K个最高的值、supp(T(π,2K))为T(π,2K)的索引集合
   • 因为初始Ω里面含有K个supports，所以在经过3.后的Ω里面将含有2K至3K个supports
5. 挑出需要运算的矩阵A
   • ${\displaystyle A_{\Omega }}$= {${\displaystyle {\vec {\theta _{j}}}|j\in \Omega }$}，${\displaystyle {\vec {\theta _{j}}}}$为第j列的矩阵A
6. 利用最小二乘法(least square,LS)计算出${\displaystyle x_{i|\Omega }(t+1)}$
   • ${\displaystyle x_{i|\Omega }(t+1)}$ = ${\displaystyle A_{\Omega }^{\dagger }y=(A_{\Omega }^{T}A_{\Omega })^{-1}A_{\Omega }^{T}y}$ ,${\displaystyle x_{i|\Omega ^{C}}(t+1)}$= 0
7. 以每一列绝对值的大小排列Ω里面的每一列，并只保留前K个最大的supports，把剩余的删去
   • Ω = supp${\displaystyle (T(|x_{i|\Omega }(t+1)|,K))}$，${\displaystyle x_{i|\Omega ^{C}}(t+1)}$= 0
8. 剩余值(residual)更新
   • r(t+1) = y - Ax(t+1)
9. 如达到终止条件，结束该还原，返还x(t)；否则回到1.继续。
   • 终止条件： ${\displaystyle (\lVert r(t)\rVert _{2}\leq THR)}$或者 (t=${\displaystyle Iter_{max}}$)

子空间追踪(Subspace Pursuit,SP)

子空间追踪还原算法和压缩采样匹配追踪还原算法相似，它们之间只有三个主要的差别。

1. 只从相关向量π里面寻找出K个最高correlation值的supports，而不是K个。
   • Ω = Ω ᑌ supp(T(π,K))
2. 在把K个supports移除过后，再进行一次最小二乘法，让预估计算的值更精准。
3. 终止条件的更改。因为压缩采样匹配追踪还原算法需要额外定义的THR来当做终止条件，故子空间追踪用前一次迭代的剩余值(residual)来取代THR。
   • ${\displaystyle \lVert r_{past}\rVert _{2}\leq \lVert r\rVert _{2}}$

利用剩余值(residual)是否有在迭代过程在降低来当成循环终止条件，因此使得子空间追踪收敛的速度比压缩采样匹配追踪来的更快。