最近認識了研學論壇小波版塊的一位很牛的學長--沙威，他目前在香港大學高效計算方法研究小組做博士後研究。在交流中，沙威學長介紹了近年訊號處理領域的一個新理論--壓縮感測。這個理論能有效地解決傳統訊號編碼技術在處理速度、儲存空間和抗干擾能力等方面的問題，被譽為是訊號處理領域的一個‘Big Idea’。下面這篇文章是沙威學長關於壓縮感測的引論，簡明易懂，深入淺出，有助我們快速瞭解壓縮感測的基本原理。在此轉帖全文，希望能有更多的朋友一起致力研究 n_n

注：該文包含大量的數學公式，沙威學長用了類似Latex的語法來編寫，一個更好閱讀的pdf版本放在了他的個人網頁上，並提供了相應的演算法程式碼。（程式和帖子的PDF版下載地址：

  
      到了香港做博士後，我的研究領域繼續鎖定計算電磁學。我遠離小波和計算時諧分析有很長的時間了。儘管如此，我仍然陸續收到一些年輕學者和工程人員的來信，和我探討有關的問題。對這個領域，我在理論上幾乎有很少的貢獻，唯一可以拿出手的就是幾篇網上釋出的帖子和一些簡單的入門級的程式碼。但是，從小波和其相關的理論學習中，我真正懂得了一些有趣的知識，並獲益良多。我深切地感到越來越多的學者和工程師開始使用這個工具解決實際的問題，我也發現網際網路上關於這方面的話題多了起來。但是，我們永遠不能只停留在某個階段，因為當今學術界的知識更新實在太快。就像我們學習了一代小波，就要學習二代小波；學習了二代小波，就要繼續學習方向性小波（X-let）。我也是在某個特殊的巧合下不斷地學習某方面的知識。就像最近，我的一個友人讓我幫她看看“壓縮感測”（Compressive Sensing）這個話題的時候，我的興趣又一次來了。我花了一個星期，閱讀文獻、思考問題、程式設計序、直到寫出今天的帖子。我希望這篇帖子，能對那些沒進入且迫切想進入這個領域的學者和工程師有所幫助。並且，我也希望和我一個星期前一樣，對這個訊號處理學界的“一個大想法”（A Big Idea）絲毫不瞭解的人，可以嘗試去了解它。我更希望，大家可以和我探討這個問題，因為我到現在甚至不完全確定我對壓縮感測的某些觀點是否正確，儘管我的簡單的不到50行的程式碼工作良好。在這個領域中，華裔數學家陶哲軒和斯坦福大學的統計學家David Donoho教授做出了重要的貢獻。在這個引言中，我用簡單的關鍵字，給出我們為什麼需要了解甚至是研究這個領域的原因。那是因為，我們從中可以學習到，下面的這些：矩陣分析、統計概率論、拓撲幾何、優化與運籌學、泛函分析、時諧分析（傅立葉變換、小波變換、方向性小波變換、框架小波）、訊號處理、影象處理等等。所以，我們有什麼理由，拒絕這個有意思的東西呢？讓我們開始吧。
 
傳統思路——正交變換
      對於一維的訊號$x/in R^{N/times1}$，大多數情況下，資訊是冗餘的。我們可以通過正交變換的方法來壓縮它。正變換：$y=/Psi x$，反變換$x=/Psi^{H}y$。這裡，$/Psi/Psi^{H}=/Psi^{H}/Psi=I$，$/Psi/in C^{N/timesN}$，$I$是單位矩陣。對於$y/in C^{N/times1}$，能量較$x$集中，本質上去除了$x$中的相關性。因此，我們只保留$K$個較大的分量，而把其它$N-K$個置為零。通過反變換，我們能夠近乎完美的重建原始訊號。因為，那$N-K$個變換域係數的貢獻，實在微乎其微。具有這樣性質的訊號被稱為$K$“稀疏”的。於是，我們有了如下編碼解碼的策略：

編碼：構造$/Psi$，做正變換$y=/Psi x$，保留$y$中最重要的$K$個分量，和其對應的位置。
解碼：把$K$個分量放回到對應的位置，其它位置填$0$，構造$/Psi^{H}$，反變換$/hat x=/Psi^{H}/hat y$。

而解碼能否近乎得到原始訊號呢？顯然，我們希望$||x-/hat x||_{2}=||y-/hat y||_{2}/leq/delta$，$/delta$是一個小的常數。但更有效的是用相對誤差$||y-/hat y||_{2}/||y||_{2}/leq/delta$。
 
      但這種編碼解碼方法有些缺點：1、考慮到夏農（Shannon）取樣定理，為了獲得很好的訊號解析度，取樣間隔會很小，造成了原始訊號長度會很長，因此變換過程會消耗很長的時間。2、$K$個需要保留的重要分量的位置，是隨著訊號的不同而不同的。因此，這種策略是“自適應”的，且需要分配多餘的空間儲存這些位置。3、一旦在傳輸過程中$K$個分量中的某幾個丟失了，後果可想而知。如果我們製作一個音訊裝置，1將帶來電力的消耗和使用者的不滿，2將帶來儲存空間的增加，3將帶來較差的抗干擾能力。
 
新的思路——壓縮感測
      壓縮感測（Compressive Sensing）是一個很有意思的新的方向。它也正成為訊號處理領域的“A Big Idea”。對於訊號$x/in R^{N/times 1}$，我們可以找到它的$M$個線性測量（Liner Measurement），$s=/Phi x$，$/Phi/in R^{M/times N}$。擁有了這$M$個測量和$/Phi$，我們就可以近乎完美的重構原始訊號了。聽起來“相當”傳奇，事實上，它基於如下嚴格的數學最優化（Optimization）問題：

目標函式$min ||/hat y||_{0}$，且滿足等式約束$/Phi/Psi^{H}/hat y=s$
或者，可以寫成
$min ||s-/Phi/Psi^{H}/hat y||_{2}+/lambda||/hat y||_{0}$

求解該最優化問題，得到變換域的$/hat y$，然後反變換，便可以得到時域的$/hat x$。公式中的2是我們熟悉的2-範數，而0是什麼呢？是0-範數，也就是向量$/hat y$中非零元素的個數。看起來很有道理，因為$/hat y$是待求的變換域向量，它是$K$稀疏的。使$/hat y$非零元素的個數儘量小，也就是保留了儘量少的重要的$K$個分量，顯然這幾個分量可以近乎完美重構$x$。我們回到傳統的思路，這$K$個分量是我們在變換域“自適應”找的，而該優化演算法也可以使我們找到這$K$個分量。

      這就足夠了嗎？顯然不行，我們仍然沒有探討測量矩陣$/Phi$需要滿足的性質。我們用極限分析法。如果我們把$/Phi$構造成和$/Psi^{H}$極端相似的矩陣，也就是拿出$/Psi$的前$M$行。用這個演算法求$/hat y$，我們將得到$/hat y=(s^{M/times 1};0^{(N-M)/times 1})$，這顯然是錯誤的。也就是說，你強迫的認為前$M$個變換域分量是重要的。而事實是，重要的$K$個分量的位置我們事先是不知道的，是隨著訊號的不同而不同的。當然，你可以將$/Phi$恰好構造成對應最重要分量的$K$行，得到正確的結果。而這種的做法要付出的概率代價$1/C_{N}^{K}$。也就是說，你必須窮舉$C_{N}^{K}$次，才能得到你想要的結果。但是，即使你有幸碰到了它，也並不能肯定這個結果就是對的。因此，我們選擇$/Phi$和$/Psi^{H}$極端不相似。於是，$/Phi$可以是滿足高斯分佈的白噪聲矩陣，或貝努裡分佈的 矩陣等等。除此之外，我們希望線性測量有穩定的能量性質：$1-/delta /leq||/Phi/Psi^{H}/hat y||_{2}/||/hat y||_{2}/leq 1+/delta $，也就是它要保持$K$個重要分量的長度。綜合上面的，我們有了如下編碼解碼的策略：

編碼：構造$/Phi$，生成測量$s=/Phi x$，保留$s$。
解碼：構造同樣的$/Phi$，構造任一種正交變換$/Psi^{H}$，根據$s$重構$x$。

      壓縮感測的優勢：1、非自適應的，一開始就可以傳輸長度較短的訊號，甚至突破取樣定理的極限。2、抗干擾，$s$中任何一項都是重要的，或者說不重要的。丟失了某幾項，仍然可以完美重構。它的缺點：1、實際中，$s$的長度一般是原始訊號的$1/4$，才能近乎完美重構。數學上更嚴格的，$M/approx K/log(N/K)$。2、重構（恢復）演算法是NP問題。即使將0-範數轉化為1-範數，由於其不可微性，演算法的計算複雜度仍然很高。它的應用前景廣泛：低成本數碼相機和音訊採集裝置；節電型音訊和影象採集裝置；天文（影象本身就稀疏，例如天空的星星）；網路；軍事（用很簡易的攝像機隨機記錄場景，可以完全重構軍事地圖）；超寬頻（雷達訊號處理）。這裡，值得指出的是，美國的工程學家已經設計出了實際的產品。
 
快速演算法——正交匹配追蹤
      對於0-範數的優化問題，實際上是NP問題，就是在多項式時間內難以求解，甚至無法驗證解的可靠性。於是，我們必須將0-範數換一下，變成1-範數。為什麼不是2-範數呢，那樣就會簡單多了。畢竟2-範數的優化問題可以轉化成2次型問題，而1-範數，$||/hat y||_{1}=/sum_{i}|/hat y_{i}|$，在$0$點處不可導，因此無論是梯度演算法、矩陣求導等等手段都變得相形見絀。因此，基於1-範數的優化演算法需要特殊處理，且複雜度很高。下面我們來解釋下為什麼要用1-範數，而不是2-範數。我們令恢復（Recovery）矩陣$T=/Phi/Psi^{H}$，則等式約束可重寫為：$T/hat y=s$。$/hat y$中未知數有$N$個，方程只有$M$個，且$M/ll N$。因此，方程有無窮多解。從幾何上說，$T/hat y-s=0$是一個超平面，為了簡化，在2-D問題中（$K=1$，$/hat y$只有兩個元素待求）可認為它就是一條直線。而範數約束呢？0-範數是一個十字架，因此它的最外側（範數的最小值）是4個點。所以其和直線的交點，必然在座標軸上。也就是說，能使$/hat y$產生更多的$0$，這正是我們想要的“稀疏”的結果。2範數是一個圓，因此它的最外側邊界和直線的交點（就是切線的概念），以壓倒性的概率不在座標軸上，除了直線的斜率恰好為$0$或者無窮大。其實直線的斜率恰好為$0$或者無窮大，是不可能的，因為$/Phi$和$/Psi^{H}$極端不相似。只有$/Phi$取$/Psi$的某一行時，兩者相似，才會發生斜率恰好為$0$或者無窮大的情況（因此，你的勝算只有$1/C_{2}^{1}$，但你不知道哪個是對的。）。依上所述，用2-範數優化的結果，使$/hat y$幾乎沒有$0$，這是我們不期望的。而1-範數是一個菱形，四個角都在座標軸上，因此它和直線的交點以壓倒性的概率落在座標軸上。這就是我們要用1-範數的原因。事實上，$p$範數滿足，$0/leq p<1$，它的外邊界都向中心凹。而$p>1$，外邊界向外凸。所以前者的外邊界和直線的交點無疑落在座標軸上。根據這個幾何解釋，我們可以將問題轉化成：

$min ||s-T/hat y||_{2}+/lambda||/hat y||_{1}$

      眾所周知，對於優化問題，我們一般用梯度的方法來求解。而對$||/hat y||_{1}$，在$0$點導數不存在，因為這個點正好位於兩條直線的交點上，左右導數不相等。這也正是很多數學方法的考慮。像子梯度（Subgradient）法、平滑近似法（Smooth Approximation）等等。我們暫不談這些方法，因為它需要特殊的數學背景。我們談一談工程領域最常用的正交匹配追蹤法（Orthogonal Matching Pursuit）。它的思想本質上還是來自於這個$K$“稀疏”。我們繞了一圈，還是為了找這$K$個關鍵的分量。既然是關鍵，顯然它的係數的絕對值應該比其它$N-K$個分量大得多。為了簡單起見，我們先假設$K=1$，唯一非零元素$/hat y_{q}$在$/hat y$中對應的位置在$q$。於是$T/hat y$就是恢復矩陣$T$的第$q$列$T_{q}$與$/hat y$中的非零元素$/hat y_{q}$的乘積，即$T_{q}/hat y_{q}=s$。且$||s-s_{q}||_{2}/||s||_{2}/leq/delta$。換句話說，$T$的第$q$列與$s$的相似程度最高，即$|<T_{q},s>|=|T_{q}^{H}s|/gg|T_{r}^{H}s|=|<T_{r},s>|,r/neq q$。所以，我們只要計算恢復矩陣$T$的所有列與$s$的內積，找到內積絕對值最大的那列就行了，該列對應的位置就是$q$。根據最小二乘法，$/hat y_{q}=(T_{q}^{H}T_{q})^{-1}T_{q}^{H}s$，就是使$||s-T_{q}/hat y_{q}||_{2}$最小的那個$/hat y_{q}$。呵呵，感覺到了嗎，實際上這有點或者非常像施密特（Schimidt）正交化方法。餘量$r_{n}=s-<T_{q},s>/< T_{q},T_{q}>T_{q}$，始終同$T_{q}$正交。這也是為什麼這個方法叫“正交”匹配追蹤的意思了。而匹配，顯然是你找到了最大的$|<T_{q},s>|$。對於$K>1$呢，其實是差不多的。我們找到餘量$r_{n}$同$T$中所有列向量最大的那個即可（但第一次找到的那列要排除，因為它已經保留了下來。）。於是，找到使||s-(T_{q2},T_{q1})(/hat y_{q2};/hat y_{q1})||_{2}$最小的那個$(/hat y_{q2};/hat y_{q1})$。這裡，$T_{q1}$是我們第一次找到的那一列，$T_{q2}$是我們新找到的那一列（也要記住它的列號$q_{2}$）。可以看出，$/hat y_{q}=(/hat y_{q2};/hat y_{q1})$被更新了，由原來的一個變成兩個了，也就是我們找到兩個在變換域最關鍵的元素和其在$/hat y$中對應的位置了。令$T_{q}=(T_{q2},T_{q1})$，餘量$r_{n}$又一次被寫為：$r_{n}=s-<T_{q},s>/< T_{q},T_{q}>T_{q}$。我們再一次看到了施密特正交化的影子。繼續上面的步驟，直至找到變換域所有$K$個最重要的分量。也就是說，正交匹配追蹤的迭代次數$m/geqK$。實際操作上只要滿足$||r_{n}||_{2}/||s||_{2}/leq/delta$，迭代就可以中止了。最後，假設正交變換$/Psi$是傅立葉變換，我們來分析下計算複雜度。生成恢復矩陣$T$，採用FFT快速演算法，需要計算複雜度$O(MN/log(N))$；找到$T$中匹配的列向量需要$O(mNM)$次操作，注意$m$是演算法需要的迭代次數；求解$/hat y_{q}$和$r_{n}$需要$O(m^{2}M)$次操作。考慮到$m~K$，$M~K/log(N/K)$，總計算複雜度可以控制在$O(N/log(N)_{2})$以下。到此為止，我們完成了壓縮感測所有基礎的理論和操作。

      為了大家更好的理解壓縮感測，在我的個人網站上，給出了一個沒有經過優化的簡單的操作程式碼，希望對大家有所幫助。此外，由於此帖含有大量公式，我用了類似Latex的語法敲出了它，一個更好閱讀的pdf版本也放在了我的個人網頁上。最後，如果閱讀這篇帖子的人，試圖轉貼的話，請註明這篇帖子的出處和作者，以此作為對作者工作小小的肯定和支援。（程式和帖子的PDF版下載地址：
http://www.eee.hku.hk/~wsha/Freecode/freecode.htm）