查看“︁壓縮感知還原演算法”︁的源代码

{{cleanup-jargon|time=2017-12-19T06:06:21+00:00}}
{{copyedit|time=2017-12-19T06:06:21+00:00}}
'''[[壓縮感知]]'''(Compressive Sensing, CS)是一種新興的取樣方法，相比長久以來被使用的[[奈奎斯特采样定理|奈奎斯特取樣定理]](Nyquist Sampling Theorem)，能更高效的方式採樣信號。壓縮感知最主要利用信號的稀疏性來尋找欠定線性系統的稀疏解，因此能從較少的取樣樣本中還原信號。近幾年有許多文獻提出了許多有效的算法，大多是透過迭代的方式找到係數並還原信號，相較同時找到最大係數的方式，能更精確的還原信號。

然而，壓縮感知的運作機制是建立於可同時獲得取樣後的樣本振幅與相位資訊的前提下，方能完美重建原始稀疏訊號，因此在一些現實的應用情境中，例如X射線晶體學、全像攝影、量子光學、毫米波無線通訊…等，若其佈署的感測元件或系統無法取得樣本的相位資
訊時(或相位資訊不可靠)，則會嚴重影響壓縮感知訊號重建的效果。為了克服此難題，近幾年許多研究致力於將建立在線性測量模型(Linear measurement model)下的壓縮感知理論與相關的訊號重建技術開發延伸至非線性測量系統模型；特別是針對只知道線性量測資料的振幅資訊這類的非線性模型進行深入探討。由於是在不知道量測值的相位資訊條件下完成稀疏訊號還原，因此這類的研究課題被稱為壓縮相位回復(Compressive phase retrieval)。許多高能源效益的稀疏訊號還原方法相繼被提出，例如基於凸優化理論所開發的 PhaseLift與PhaseCut還原方法，或是GESPAR以及SPARTA等貪婪式演算法(Greedy algorithms)。 

==還原演算法 ==
對於壓縮感知中的稀疏信號還原有很多不同類型的演算法，這些算法大致可以被分為六類，以下將詳細敘述各類型演算法的特色<ref name=":0">{{Cite journal|title=Compressive sensing: From theory to applications, a survey|url=http://ieeexplore.ieee.org/document/6674179/|last=Qaisar|first=S.|last2=Bilal|first2=R. M.|date=October 2013|journal=Journal of Communications and Networks|issue=5|doi=10.1109/JCN.2013.000083|volume=15|pages=443–456|issn=1229-2370|last3=Iqbal|first3=W.|last4=Naureen|first4=M.|last5=Lee|first5=S.|access-date=2017-12-19|archive-date=2018-11-28|archive-url=https://web.archive.org/web/20181128150501/https://ieeexplore.ieee.org/document/6674179/|dead-url=no}}</ref>。

===凸鬆弛(Convex Relaxation)===
這類算法最主要透過線性規劃方式求解凸優化的問題來對信號進行還原。 雖然要精確地還原信號所需的測量數量很少，但是這些方法在計算上是較複雜的。主要代表的演算法有基追蹤(Basis Pursuit, BP)、抗噪基追蹤(Basis Pursuit De-noising, BPDN)、最小角度回歸(Least Angle Regression, LARS)。

基追蹤優化公式：<math>\min_x \|x\|_1 \quad \mbox{subject to} \quad y = Ax.</math>

抗噪基追蹤優化公式：<math>\min_x \frac{1}{2}\|y-Ax\|^2_2+\lambda\|x\|_1,</math>

===貪婪迭代演算法(Greedy Iterative Algorithm)===
這類算法主要是透過迭代的方式逐步找到答案來還原信號。這類型的演算法在每次迭代中，以貪婪的方式選擇傳感矩陣<math>\Theta</math>中與<math>Y</math>做內積最大的列，也就是選出與傳感矩陣中最相關的列，接著該列對<math>Y</math>的貢獻將會從<math>Y</math>中被減去，並且對相減之後的結果進行迭代，直到識別出正確的列集合。相比之下，最小平方誤差的方法則是在每次迭代中最小化誤差。這類演算法的優點在於運算複雜度低，還原速度快，常見的演算法有正交匹配追蹤(Orthogonal Matching Pursuit, OMP)、隨機梯度追蹤(Stochastic Gradient Pursuit,SGP)、正則化正交匹配追蹤(Regularized OMP)、壓縮採樣匹配追蹤(Compressive Sampling Matching Pursuit, CoSaMP)和子空間追蹤(Subspace Pursuit,SP)。

==== 正交匹配追蹤(Orthogonal Matching Pursuit, OMP) ====
<gallery widths="430" heights="250">
File:Orthogonal Matching Pursuit.png|正交匹配追蹤(OMP)演算法
</gallery>一開始會有一個空集合<math>\Lambda_0</math>，以及剩餘的部分<math>r_0</math>，每次疊代都會從<math>\overline{\Lambda{_{l-1}}}</math>找出一個A矩陣投影到剩餘<math>r_{l-1}</math>有最大值的位置，把這個位置加到<math>\Lambda_l</math>之中，並從<math>\overline{\Lambda{_{l}}}</math>當中移除，最後再更新剩餘<math>r_{l}</math>。利用疊代的方式，不斷地找出x向量當中最有可能非零的位置，直到達到演算法停止條件。

==== 隨機梯度追蹤(Stochastic Gradient Pursuit,SGP) ====
可視為是正交匹配追蹤的改良版，在面對雜訊的環境下也能更好地還原。整體演算法都與正交匹配追蹤相同，除了用最小二乘法(least square,LS)那一項。

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

* <math>x_{ZF}=A_\Omega^\dagger y = A_\Omega^\dagger(Ax+n)=x+A_\Omega^T(A_\Omega A_\Omega^T)^{-1}n</math>,'''n'''為雜訊

但是零强制线性均衡的缺點可以直接從上式看出，會造成放大雜訊的影響。為了解決雜訊的問題，可以利用最小均方误差(minimum mean square error,MMSE)均衡器來最小化整體的誤差

* <math>x_{MMSE}=A_\Omega^T(A_\Omega A_\Omega^T+\frac{\sigma_n^2}{\sigma_{x|\Omega}^2}I)^{-1}y</math>

和零强制线性均衡相比，在反矩陣內部多增加了額外的項，可直接視為雜訊項目。雖然最小均方误差方法的正交匹配追蹤<ref>{{Cite journal|title=MMSE-based version of OMP for recovery of discrete-valued sparse signals|url=http://dx.doi.org/10.1049/el.2015.0924|last=Sparrer|first=S.|last2=Fischer|first2=R.F.H.|date=2016-01-08|journal=Electronics Letters|issue=1|doi=10.1049/el.2015.0924|volume=52|pages=75–77|issn=0013-5194}}</ref>更能夠對抗雜訊干擾，但是<math>{\sigma_{x|\Omega}^2}</math>和<math>{\sigma_n^2}</math>並不適用壓縮感知的眾多應用上，所以並不能直接套用在壓縮感知的還原演算法上。

因為最小均方(least mean square,LMS)在不斷地迭代後的平衡值會接近最小均方误差，因此隨機梯度追蹤利用最小均方來取代最小二乘法。最小均方的迭代將根據預估誤差來調整support係數，而預估誤差可表示為：

* <math>e_\lambda=y_\lambda - \tilde{A_\lambda}\tilde{x_\lambda}</math>

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

* <math>\tilde{x}_{\lambda+1}=\tilde{x}_\lambda+\mu e_\lambda A_\Omega^T</math>，<math>\mu</math>為step size

總結來說隨機梯度追蹤演算法裡面包含了兩個迭代循環。一個大的迭代是由正交匹配追蹤所形成，而小的循環則被大的迭代所包含，是由最小均方所形成，不斷地迭代來更新<math>\tilde{x}_{\lambda}</math>的值。

==== 壓縮採樣匹配追蹤(Compressive Sampling Matching Pursuit, CoSaMP) ====
在進行該還原演算法前需要額外的輸入參數-稀疏性數量K(信號非零項的數量)。和一般的貪婪迭代演算法一樣，整個還原過程同樣是由(1)匹配追蹤(找出內積最大的列)、(2)還原值估計(通過最大內積列來估計還原值)、(3)剩餘值(residual)更新(通過估計還原值來更新誤差)不停地迭代組成。

和正交匹配追蹤不同的是正交匹配追蹤在每次迭代的匹配追蹤上只尋找並增加一個內積最大的列；而壓縮採樣匹配追蹤則在每次迭代的匹配追蹤上尋找並增加2K個內積最大的列。

壓縮採樣匹配追蹤過程：

參數定義：'''''π'''''為相關向量、'''r'''(t)為在第t次迭代的剩餘值(residual)、ω為在'''''π'''''裡面最大的值、'''Ω'''為support的集合

# 初始化，<math>x(0) = 0_{Px1},r(0)=y,\Omega(0)=\emptyset,t=0</math>
# 計算相關向量'''''π'''''
#* '''''π''''' = |<math>A^T</math>'''r'''(t)|
# 從相關向量'''''π'''''裡面尋找出2K個最高correlation值的supports
#* ω <math>\equiv</math><math>\underset{\omega\in[P]}{\arg max}</math>x{ | '''''π'''''(ω) | }
# 把找到的support加入support 集合裡面
#* '''Ω''' = '''Ω''' ᑌ supp(T('''''π''''',2K))，T('''''π''''',2K)為在'''''π'''''裡面2K個最高的值、supp(T('''''π''''',2K))為T('''''π''''',2K)的索引集合
#* 因為初始'''Ω'''裡面含有K個supports，所以在經過3.後的'''Ω'''裡面將含有2K至3K個supports
# 挑出需要運算的矩陣'''A'''
#* <math>A_\Omega</math>= {<math>\vec{\theta_j} | j\in\Omega</math>}，<math>\vec{\theta_j}</math>為第j列的矩陣'''A'''
# 利用最小二乘法(least square,LS)計算出<math>x_{i|\Omega}(t+1)</math>
#* <math>x_{i|\Omega}(t+1)</math> = <math>A_\Omega^\dagger y = (A_\Omega^TA_\Omega)^{-1}A_\Omega^Ty</math> ,<math>x_{i|\Omega^C}(t+1)</math>= '''0'''
# 以每一列絕對值的大小排列'''Ω'''裡面的每一列，並只保留前K個最大的supports，把剩餘的刪去
#* '''Ω''' = supp<math>(T(|x_{i|\Omega}(t+1)|,K))</math>，<math>x_{i|\Omega^C}(t+1)</math>= '''0'''
# 剩餘值(residual)更新
#* '''r'''(t+1) = '''y''' - '''Ax'''(t+1)
# 如達到終止條件，結束該還原，返還'''x'''(t)；否則回到1.繼續。
#* 終止條件： <math>(\lVert r(t) \rVert_2 \leq THR)</math>或者 (t=<math>Iter_{max}</math>)

==== 子空間追蹤(Subspace Pursuit,SP) ====
子空間追蹤還原演算法和壓縮採樣匹配追蹤還原演算法相似，它們之間只有三個主要的差別。

# 只從相關向量'''''π'''''裡面尋找出K個最高correlation值的supports，而不是K個。
#* '''Ω''' = '''Ω''' ᑌ supp(T('''''π''''',K))
# 在把K個supports移除過後，再進行一次最小二乘法，讓預估計算的值更精準。
# 終止條件的更改。因為壓縮採樣匹配追蹤還原演算法需要額外定義的THR來當做終止條件，故子空間追蹤用前一次迭代的剩餘值(residual)來取代THR。
#* <math>\lVert r_{past} \rVert_2 \leq \lVert r \rVert_2</math>

利用剩餘值(residual)是否有在迭代過程在降低來當成循環終止條件，因此使得子空間追蹤收斂的速度比壓縮採樣匹配追蹤來的更快。

===迭代閾值演算法(Iterative Thresholding Algorithm) ===
相較於凸鬆弛演算法，迭代方法在壓縮感知還原信號的速度更快。對於這類算法，主要透過過軟閾值或硬閾值來正確的重建信號，而閾值的大小則會根據迭代次數進行調整。 最近提出了擴展匹配追踪(Expander Matching Pursuits)、稀疏匹配追踪(Sparse Matching Pursuit)和序列稀疏匹配追蹤(Sequential Sparse Matching Pursuits)，實現了近線性還原時間。

迭代閾值算法是從Relaxation Method啟發而來，Relaxation Method是用於解決具有稀疏性的線性系統。

在迭代閾值算法中，問題一樣是被定義為：<div><math>\min \|x\|_1\ \ subject\ to\ y=Ax </math></blockquote></div>
而在此算法中，每一次的迭代主要分成兩部分：1、尋找新的解<math>x^{t+1} </math>。2、更新殘值<math>r^t </math>。

<div>Step 1: 主要是利用投影的方式找到更新的向量方向，再透過取閾值的方式來進行優化。</div>
<math>x^{t+1}=\eta_{thr}(x^t+\mathcal{K}\times(A^Tr^t)) </math></blockquote></div>
而<math>\mathcal{K} </math>是relaxation parameter，且<math>\mathcal{K}\in(0,1) </math>。

而<math>\eta_{thr}(.) </math>就是閾值函數(thresholding function)，主要就是為了使迭代的過程中，逐漸逼近具有稀疏性性質的解，而目前主要有兩大類取閾值的方式：硬閾值函數(Hard Thresholding Function)與軟閾值函數(Soft Thresholding Function)。

硬閾值函數就是將小於閾值(thr)的解，一律通通過濾成0。
<div><math>\eta_{thr}^H(x)=\begin{cases} x, \text{if }|x|>thr\\0,\text{otherwise} \end{cases}    </math></div>
<div>軟閾值函數則是使用較為平滑的方式，將值逼近於稀疏性的解</div>
<math>\eta_{thr}^S(x)=\sgn(x)(|x|-thr)_+\text{, where }(z)_+=z\times\mathbb{I}(z\geq0)    </math>，而<math>\mathbb{I}(.)    </math>是[[指示函數]]
<div>Step 2: 利用新算出來的<math>x^{t+1} </math>來進行殘差更新</div>
<math>r^{t+1}=y-Ax^{t+1}    </math>
<div>之後一直等到規定的迴圈數抵達即完成還原。

而此算法是Landweber iteration的變形，若沒有加入閾值函數的話，就會收斂到<math>\min \|y-Ax\|_2^2 </math>。

===組合/次線性演算法(Combinatorial/Sublinear Algorithms) ===
這類算法通過分組測試來還原稀疏信號。與凸鬆弛演算法或貪婪演算法相比，它們有高速和高效的好處，然而對於測量矩陣<math>\Phi</math>有特定模式稀疏的要求。代表性的演算法有傅里葉採樣算法(Fourier Sampling Algorithm)和鍊式追蹤(Chaining Pursuits)。

===非凸最小化演算法(Non Convex Minimization Algorithms) ===
非凸局部最小化技術通過用<math>l_p</math>規範替換<math>l_1</math>規範來從更少的測量中恢復CS信號，其中<math>p\leq1</math>。非凸優化主要應用於醫學影像斷層掃描、網絡狀態推斷、資料串流壓縮。文獻中提出了許多使用這種技術的演算法，例如聚焦欠定系統解法(Focal Underdetermined System Solution, FOCUSS)、迭代重新加權最小平方法(Iterative Re-weighted Least Squares)、稀疏貝葉斯學習算法(Sparse Bayesian Learning Algorithms)、基於蒙特卡羅(Monte-Carlo based Algorithms)的演算法。

===布里格曼迭代演算法(Bregman Iterative Algorithm) ===
這些算法提供了一種解決最小化問題的簡單而有效的方法。文獻<ref>{{Cite journal|title=Bregman Iterative Algorithms for $\ell_1$-Minimization with Applications to Compressed Sensing|url=http://epubs.siam.org/doi/abs/10.1137/070703983|last=Yin|first=W.|last2=Osher|first2=S.|date=2008-01-01|journal=SIAM Journal on Imaging Sciences|issue=1|doi=10.1137/070703983|volume=1|pages=143–168|last3=Goldfarb|first3=D.|last4=Darbon|first4=J.}}</ref>給出了一個新的想法，它透過布里格曼迭代正則化方法產生一系列無約束子問題給出約束問題的精確解。當應用於壓縮感知問題時，使用布里格曼的迭代方法能在四到六次的迭代中實現還原。與其他現有演算法相比，這些算法的計算速度特別吸引人。

== 複雜度比較 ==
{| class="wikitable"
|+'''各類型壓縮感知還原演算法複雜度比較<ref name=":0" />'''
!種類
!演算法
!複雜度
!最小測量次數
|-
|凸鬆弛
|基追踪
|<math>O(n^3)</math>
|<math>O(s\log n)</math>
|-
| rowspan="3" |貪婪迭代演算法
|正交匹配追踪
|<math>O(smn)</math>
|<math>O(s\log n)</math>
|-
|正則化正交匹配追踪
|<math>O(smn)</math>
|<math>O(s\log^2 n)</math>
|-
|壓縮採樣匹配追踪
|<math>O(mn)</math>
|<math>O(s\log n)</math>
|-
|迭代閾值演算法
|擴展匹配追踪
|<math>O(n\log (n/s))</math>
|<math>O(n\log (n/s))</math>
|-
|組合/次線性演算法
|鍊式追踪
|<math>O(s\log^2 n\log^2s)</math>
|<math>O(s\log^2 n)</math>
|}

== 參考資料 ==
{{refbegin}}
* E. J. Candes, T. Strohmer, and V. Voroninski, “PhaseLift:  Exact and Stable Signal Recovery from Magnitude Measurements via Convex Programming,” Commun. Pure Appl. Math., vol. 66, no. 8, pp. 1241–1274, 2013

{{refbegin}}
* I. Waldspurger, A. d’Aspremont, and S. Mallat, “Phase Recovery, MaxCut, and Complex Semidefinite Programming,”Math. Program. A, vol. 149, no. 1, pp.47–81, 2015.

{{refbegin}}
* Y. Shechtman, A. Beck, and Y. C. Eldar, “GESPAR: Efficient Phase Retrieval of Sparse Signals,” IEEE Trans. Signal Processing, vol. 62, no. 4, pp. 928–938, Feb. 2014.

{{refbegin}}
* G. Wang, L. Zhang, G. B. Giannakis, M. Akc¸akaya and J. Chen, “Sparse Phase Retrieval via Truncated Amplitude Flow,” IEEE Trans. Signal Processing, vol. 66, no. 2, pp. 479–491, Jan. 2018