[綜] Sparse Representation 稀疏表示壓縮感知

本文轉載自查看原文 2015-06-28 00:40 2038 Matlab 應用

稀疏表示分為 2個過程：1. 獲得字典（訓練優化字典；直接給出字典），其中字典學習又分為2個步驟：Sparse Coding和Dictionary Update；2. 用得到超完備字典后，對測試數據進行稀疏編碼Sparse Coding，求出稀疏矩陣。

1. 訓練字典的方法：MOD，K-SVD，Online ...

MOD (Method of Optimal Direction):

Sparse Coding其采用的方法是OMP貪婪算法;

Dictionary Update采用的是最小二乘法，即D=argmin norm(y-D*x,2)^2 解的形式是D=Y*x'*inv(x*x’).

2. OMP(Orthogonal Matching Pursuit)

計算機視覺小菜鳥的專欄

稀疏表示

http://blog.csdn.net/carson2005/article/details/7989675

幫助快速了解稀疏表示。

稀疏表示是最近幾年信號處理領域的熱點之一，簡單來說，它其實是一種對原始信號的分解過程，該分解過程借助一個事先得到的字典（也有人稱之為過完備基，overcomplete basis，后面會介紹到），將輸入信號表示為字典的線性近似的過程。即：

秋楚驛的博客

[轉載]稀疏表示step by step(轉)

http://blog.sina.com.cn/s/blog_60542f6601018udw.html

系統講解了稀疏表示。

稀疏表示step by step(1)

聲明：本人屬於絕對的新手，剛剛接觸“稀疏表示”這個領域。之所以寫下以下的若干個連載，是鼓勵自己不要急功近利，而要步步為贏！所以下文肯定有所紕漏，敬請指出，我們共同進步！

踏入“稀疏表達”（Sparse Representation）這個領域，純屬偶然中的必然。之前一直在研究壓縮感知（Compressed Sensing）中的重構問題。照常理來講，首先會找一維的稀疏信號（如下圖）來驗證CS理論中的一些原理，性質和算法，如測量矩陣為高斯隨機矩陣，貝努利矩陣，亞高斯矩陣時使用BP,MP,OMP等重構算法的異同和效果。然后會找來二維稀疏信號來驗證一些問題。當然，就像你所想的，這些都太簡單。是的，接下來你肯定會考慮對於二維的稠密信號呢，如一幅lena圖像？我們知道CS理論之所以能突破乃奎斯特采樣定律，使用更少的采樣信號來精確的還原原始信號，其中一個重要的先驗知識就是該信號的稀疏性，不管是本身稀疏，還是在變換域稀疏的。因此我們需要對二維的稠密信號稀疏化之后才能使用CS的理論完成重構。問題來了，對於lena圖像這樣一個二維的信號，其怎樣稀疏表示，在哪個變換域上是稀疏的，稀疏后又是什么？於是竭盡全力的google…后來發現了馬毅的“Image Super-Resolution via Sparse Representation”（IEEE Transactions on Image Processing，Nov.2010）這篇文章，於是與稀疏表達的緣分開始啦！ [break] 談到稀疏表示就不能不提下面兩位的團隊，Yi Ma AND Elad Michael，國內很多高校（像TSinghua，USTC）的學生直奔兩位而去。（下圖是Elad M的團隊，后來知道了CS界大牛Donoho是Elad M的老師，怪不得…）其實對於馬毅，之前稍有了解，因為韋穗老師，我們實驗室的主任從前兩年開始着手人臉識別這一領域並且取得了不錯的成績，人臉識別這個領域馬毅算是大牛了…因此每次開會遇到相關的問題，韋老師總會提到馬毅，於是通過各種渠道也了解了一些有關他科研和個人的信息。至於Elad.M，恕我直言，我在踏入這個領域其實真的完全不知道，只是最近文章看的比較多，發現看的文章中大部分的作者都有Elad，於是乎，好奇心驅使我了解了這位大牛以及他的團隊成員…也深深的了解到了一個團隊對一個領域的貢獻，從Elad.M那兒畢業的學生現在都成了這個領域中的佼佼者…不禁感嘆到：一個好的導師是多么的重要！！下面舉個簡單的例子，說說二維信號的稀疏性，也為后面將稀疏表示做個鋪墊。我們以一幅大小為256×256的Lena圖像為例，過完備字典（Dictionary，具體含義見后文，先理解為基吧，其實不完全等同）選擇離散余弦變換DCT，字典大小選擇64×256，對圖像進行分塊處理，由於僅僅為了說明稀疏性的概念，所以不進行重疊處理，每塊大小8×8（pixel）…簡單的稀疏表示后的稀疏如下圖。可以看出，絕大多數的稀疏集中在0附近。當然，這里僅僅是簡單的說明一下。后面我們有更好的選擇，比如說，字典的選擇，圖像塊的選擇等等…

本文固定鏈接: http://www.win7soft.com/justplus/srstepbystep1 | JustPlus

-----------------------------------------------------

稀疏表示step by step(2)

壓縮感知（CS），或許你最近聽說的比較多，不錯，CS最近比較火，什么問題不管三七二十一就往上粘連，先試試能不能解決遇到的問題，能的話就把文章發出來忽悠大家，這就是中國學術浮躁的表現…我們沒有時間去思考的更多，因為你一思考，別人可能就把“你的東西”搶先發表了…不扯了，反正也干預不了…稀疏表示的現狀有點像CS，能做很多事，也不能做很多事…但是它確實是解決一些棘手問題的方法，至少能提供一種思路…目前用稀疏表示解決的問題主要集中在圖像去噪（Denoise），代表性paper：Image Denoise Via Sparse and Redundant Representations Over Learned Dictionaries(Elad M. and Aharon M. IEEE Trans. on Image Processing,Dec,2006);Image Sequence Denoising Via Sparse and Redundant Representations(Protter M. and Elad M.IEEE Trans. on Image Processing,Jan,2009)，還有超分辨率（Super-Resolution OR Scale-Up），代表性paper：Image Super-Resolution via Sparse Representation（Jianchao Yang, John Wright, Thomas Huang, and Yi Ma，IEEE Transactions on Image Processing, Nov,2010）,A Shrinkage Learning Approach for Single Image Super-Resolution with Overcomplete Representations( A. Adler, Y. Hel-Or, and M. Elad,ECCV,Sep,2010)…. 另外還有inpait，deblur，Face Recognition，compression等等..更多應用參考Elad M的書，google能找到電子檔，這里不提
供下載地址

當然Elad M.和Yi Ma的團隊不僅僅在應用上大做文章，在理論上也是不斷革新…就拿字典學習為例，一開始將固定字典（如過完備 DCT，Contourlet，Wavelet字典）用在去噪，超分辨率上，效果不明顯，可能某些情況下還不如空域或者頻域的去噪效果好（我拿Lena圖像和DCT實驗了，以PSNR為標准，比時域的去噪方法要好，但是比小波去噪相比，稍稍遜色），但是速度快是它的優點。於是他們開始研究自適應的字典學習算法，開始使用很多圖像進行學習，后來采用單幅圖像進行學習來提高運算速度（使用很多圖像進行學習屬於半自適應的學習，對於自然圖像的處理需要學習自然圖像，對遙感圖像的處理需要學習遙感圖像，但是對自然圖像或遙感圖像的去噪，超分辨率處理，都可以使用已經訓練好的相應的字典）；同時學習的方法也不盡相同，開始使用MOD，后來就是一直比較流行的K-SVD，最近又出來了Online，總體而言Online比較快。下面是我提到的幾種字典的例子，所有的字典都是64×256大小的，依次為DCT，globally（訓練圖像是：標准圖像 lena，boat，house，barbara，perppers），K-SVD（單幅含躁lena，噪聲標准差為25），online（單幅含躁 lena，噪聲標准差為25），其中globally的訓練方法是將訓練圖像分成8×8的overlap patch，平均取，共取10000塊，K-SVD和online也是分成相同的重疊塊，取所有可能的塊。

總之，Elad M.和Yi Ma為稀疏表示這個領域作出了很大的貢獻…向大牛們致敬！！最后稍微說一下國內的研究現狀，國內的很多研究還沒浮出水面，不知道是不是我想的的這樣（我比較疑惑的是，為什么國外研究了十幾年，國內還沒大動靜？），至少從google學術以及IEEE的文章搜索上來看是這樣的…不過還是有幾位教授在這方面作出了很大的貢獻的…

-------------------------------------------------

稀疏表示step by step(3)

我們來考慮信號的稀疏表示問題，假如我們有了過完備字典D，如何求出信號x在這個過完備字典上的稀疏表示？先來回顧一下在壓縮感知中常常會遇到的問題，信號x在經過測量矩陣A后得到測量值y，即y=A*x，其中測量矩陣A_mxn(m遠小於n)，那么怎樣從y中精確的恢復出x呢？

由於m遠小於n，用m個方程求解n個未知數，因此y=A*x是個欠定方程，有無窮多個解。就像我們解優化問題一樣，如果我們加上適當的限定條件，或者叫正則項，問題的解會變得明朗一些！這里我們加上的正則項是norm(x,0)，即使重構出的信號 x盡可能的稀疏（零范數：值為0的元素個數），后來Donoho和Elad這對師徒證明了如果A滿足某些條件，那么argmin norm(x,0) s.t.y=A*x 這個優化問題即有唯一解！唯一性確定了，仍然不能求解出該問題，后來就嘗試使用l1和l2范數來替代l0范數，華裔科學家陶哲軒和candes合作證明了在A滿足UUP原則這樣一個條件下，l0范數可以使用l1范數替代，所以優化問題變成argmin norm(x,1) s.t.y=A*x這樣一個凸優化問題，可以通過線性優化的問題來解決！（參考文獻：Stable signal recovery from incomplete and inaccurate measurements（E. J. Candès, J. Romberg and T. Tao））至此，稀疏表示的理論已經初步成型。至於之后的優化問題，都是一些變形，像lasso模型，TV模型等….這里推薦一本stanford的凸優化教材convex optimization，我准備抽個時間好好看一看，搞稀疏表達這一塊的，優化問題少不了…最近一直在感嘆：數學不好的人哪傷不起啊！！有木有！！[break]

我們再回到我們一開始提到的問題上來，一開始說字典D已知，求出y在過完備字典D上的稀疏表示x，這個在稀疏表示里面被稱作為Sparse Coding…問題的模型是 x=argmin norm(y-D*x,2)^2 s.t.norm(x,1)<=k; 我們下面來介紹一下解決這個問題的常用方法OMP(Orthogonal Matching Pursuit) .我們主要目標是找出x中最主要的K個分量（即x滿足K稀疏），不妨從第1個系數找起，假設x中僅有一個非零元x(m)，那么 y0=D(:,m)*x(m)即是在只有一個主元的情況下最接近y的情況，norm(y-y0,2)/norm(y,2)<=sigma,換句話說在只有一個非零元的情況下，D的第m列與y最“匹配”，要確定m的值，只要從D的所有列與y的內積中找到最大值所對應的D的列數即可，然后通過最小二乘法即可確定此時的稀疏系數。考慮非零元大於1的情況，其實是類似的，只要將余量r=y-y0與D的所有列做內積，找到最大值所對應D的列即可。…下面是代碼和示例結果[其中圖1是lena圖像某個8x8的圖像塊在其自身訓練得到的字典上的稀疏表示，稀疏值k=30，相對誤差norm(y- y0,2)/norm(y,2)=1.8724，圖2是相同的塊在相同的字典下的稀疏表述，稀疏值k=50，相對誤差為0.0051]

function A=OMP(D,X,L)

% 輸入參數: % D - 過完備字典，注意：必須字典的各列必須經過了規范化

% X - 信號

%L - 系數中非零元個數的最大值（可選，默認為D的列數，速度可能慢）

% 輸出參數: % A - 稀疏系數

if nargin==2

L=size(D,2);

end

P=size(X,2);

K=size(D,2);

for k=1:1:P,

a=[];

x=X(:,k);

residual=x;

indx=zeros(L,1);

for j=1:1:L,

proj=D'*residual;

[maxVal,pos]=max(abs(proj));

pos=pos(1);

indx(j)=pos;

a=pinv(D(:,indx(1:j)))*x;

residual=x-D(:,indx(1:j))*a;

if sum(residual.^2) < 1e-6

break;

end

end;

temp=zeros(K,1);

temp(indx(1:j))=a;

A(:,k)=sparse(temp);

end;

return;

end

----------------------------------------------------

稀疏表示step by step(4)

回顧一下前面所說的OMP算法，前提條件是字典D已知，求一個信號在這個字典上的稀疏表示…那假如我們不想使用過完備 DCT，wavelet呢？因為它們沒有自適應的能力，不能隨着信號的變化作出相應的變化，一經選定，對所有的信號一視同仁，當然不是我們想要的。我們想通過訓練學習的方法獲取字典D來根據輸入信號的不同作出自適應的變化。那現在我們的未知量有兩個，過完備字典D，稀疏系數x，已知量是輸入信號y，當然先驗知識是輸入信號在字典D上可以稀疏表示…我們再次列出sparse-land模型： [D,x]=argmin norm(y-D*x,2)^2 s.t.norm(x,1)<=k。如何同時獲取字典D和稀疏系數x呢？方法是將該模型分解：第一步將D固定，求出x的值，這就是你常聽到的稀疏分解（Sparse Coding），也就是上一節提到的字典D固定，求信號y在D上稀疏表示的問題；第二步是使用上一步得到的x來更新字典D，即字典更新（Dictionary Update）。如此反復迭代幾次即可得到優化的D和x。

Sparse Coding：x=argmin norm(y-D*x,2) s.t.norm(x,1)<=k
Dictinary Update:D=argmin norm(y-D*x,2)^2

[break]我們主要通過實例介紹三種方法：MOD，K-SVD,Online… 首先是MOD(Method of Optimal Direction)。Sparse Coding其采用的方法是OMP貪婪算法，Dictionary Update采用的是最小二乘法，即D=argmin norm(y-D*x,2)^2 解的形式是D=Y*x'*inv(x*x’)。因此MOD算法的流程如下：

初始化：字典D_mxn可以初始化為隨機分布的mxn的矩陣，也可以從輸入信號中隨機的選取n個列向量，下面的實驗我們選取后者。注意OMP要求字典的各列必須規范化，因此這一步我們要將字典規范化。根據輸入信號確定原子atoms的個數，即字典的列數。還有迭代次數。
主循環： Sparse Coding使用OMP算法； Dictionary Update采用最小二乘法。注意這一步得到的字典D可能會有列向量的二范數接近於0，此時為了下一次迭代應該忽略該列原子，重新選取一個服從隨機分布的原子。

function [A,x]= MOD(y,codebook_size,errGoal)

%==============================

%input parameter

% y - input signal

% codebook_size - count of atoms

%output parameter

% A - dictionary

% x - coefficent

%==============================

if(size(y,2)<codebook_size)

disp('codebook_size is too large or training samples is too small');

return;

end

% initialization

[rows,cols]=size(y);

r=randperm(cols);

A=y(:,r(1:codebook_size));

A=A./repmat(sqrt(sum(A.^2,1)),rows,1);

mod_iter=10;

% main loop

for k=1:mod_iter

% sparse coding

if nargin==2

x=OMP(A,y,5.0/6*rows);

elseif nargin==3

x=OMPerr(A,y,errGoal);

end

% update dictionary

A=y*x'/(x*x');

sumdictcol=sum(A,1);

zeroindex=find(abs(sumdictcol)<eps);

A(zeroindex)=randn(rows,length(zeroindex));

A=A./repmat(sqrt(sum(A.^2,1)),rows,1);

if(sum((y-A*x).^2,1)<=1e-6)

break;

end

----------------------------------------------------

稀疏表示step by step(5)

回憶一下前面提到字典學習的方法之一MOD，其分為兩個步驟：Sparse Coding和Dictionary Update。現在來看一種比較流行的方法K-SVD，至少到目前來說比較流行，雖然速度有點慢（迭代次數+收斂速度的影響）。之所以叫K-SVD，估計Aharon M.和Elad M.他們是從K-Means中得到了靈感，K-Means中的K是指要迭代K次，每次都要求一次均值，所以叫K均值，K-SVD也是類似，要迭代K次，每次都要計算一次SVD分解。其實在K-SVD出來之前，字典學習的方法已經有很多種，Aharon M.和Elad M.的文章中有提到，其中包括最大似然ML，MOD，最大后驗概率MAP等等，並對它們進行了算法的靈活性，復雜度的比較。 K-SVD同MOD一樣也分為Sparse Coding和Dictionary Update兩個步驟，Sparse Coding沒有什么特殊的，也是固定過完備字典D，使用各種迭代算法求信號在字典上的稀疏系數。同MOD相比，K-SVD最大的不同在字典更新這一步，K-SVD每次更新一個原子（即字典的一列）和其對應的稀疏系數，直到所有的原子更新完畢，重復迭代幾次即可得到優化的字典和稀疏系數。下面我們來具體的解釋這句話。

如上圖（左上），現在我們要更新第k個原子，即d_k..那我們需要知道在上一步迭代之后哪些信號使用了該原子，即稀疏系數不為0的部分是哪些？從左上圖中很容易看出，x'的第k列T(x_k)，也就是x的第k行中不為0的那部分所對應的T(y_k)即是我們要找的信號，結果見左下圖藍色部分。我們用d_k和稀疏系數x(k)'來重構這部分使用了d(k)的信號，它和T(y_k)的差值即E，右上圖中綠色部分，接下來我們要使用右下圖這個約束來更新x_k和d_k這兩個值…如此反復，直到過完備字典D的所有原子更新完畢為止…求解這個x_k和d_k，直接對E進行SVD分解即可。

下面是K-SVD的matlab代碼，由於是深化對原理的理解，所以沒有經過任何的優化和改進，優化的代碼可以參考K-SVD toolboxwriten by Elad M.

function [A,x]= KSVD(y,codebook_size,errGoal)

%==============================

%input parameter

% y - input signal

% codebook_size - count of atoms

%output parameter

% A - dictionary

% x - coefficent

%reference:K-SVD:An Algorithm for Designing of Overcomplete Dictionaries % for Sparse Representation,Aharon M.,Elad M.etc

%==============================

if(size(y,2)<codebook_size)

disp('codebook_size is too large or training samples is too small');

return;

end

% initialization

[rows,cols]=size(y);

r=randperm(cols);

A=y(:,r(1:codebook_size));

A=A./repmat(sqrt(sum(A.^2,1)),rows,1);

ksvd_iter=10;

for k=1:ksvd_iter % sparse coding

if nargin==2

x=OMP(A,y,5.0/6*rows);

elseif nargin==3

x=OMPerr(A,y,errGoal);

end

% update dictionary

for m=1:codebook_size

mindex=find(x(m,:));

if ~isempty(mindex)

mx=x(:,mindex);

mx(m,:)=0;

my=A*mx;

resy=y(:,mindex);

mE=resy-my;

[u,s,v]=svds(mE,1);

A(:,m)=u;

x(m,mindex)=s*v';

end

三個公開的matlab工具箱：

SPArse Modeling Software

SPAMS：

http://spams-devel.gforge.inria.fr/index.html

SparseLab：

http://sparselab.stanford.edu/

Michael Elad:

http://www.cs.technion.ac.il/~elad/

The SMALLbox, a framework for sparse representation and dictionary learning:

http://www.small-project.eu/software-data/smallbox

SPAMS稀疏建模工具箱

https://chunqiu.blog.ustc.edu.cn/?p=570

彬彬有禮的專欄

壓縮感知中的數學知識：稀疏、范數、符號arg min

http://blog.csdn.net/jbb0523/article/details/40262629

題記：從今年九月份開始看壓縮感知方面的文獻，自己的感覺是要求的數學知識太多了，上研時只學了一個矩陣理論，學的感覺還算扎實，但學完到現在基本都忘掉了，什么優化之類的根據就沒學過，所以現在看文獻真心很費勁，不過不會可以慢慢學嘛，以后我會根據自己的學習情況陸續發一些壓縮感知常用到的數學知識。

本次說三個問題：

1、稀疏

2、范數

3、符號arg min

前面兩個問題從矩陣理論的書籍中應該可以找到，最后一個問題從最優化類的書籍中應該可以找到。

=========================以下為正文=========================

1、稀疏：什么是K稀疏呢？

在壓縮感知里經常提到 “K稀疏” 的概念，這個是很容易理解的：也就是對於長度為N的向量（實際上是指一個N維離散離值信號）來說，它的N個元素值只有K個是非零的，其中K<<N，這時我們稱這個向量是K稀疏的或者說是嚴格K稀疏的；實際中要做到嚴格K稀疏不容易，一般來說，只要除了這K個值其它的值很小很小，我們就認為向量是稀疏的，這時區別於嚴格K稀疏且就叫它K稀疏吧。

為什么要談稀疏這個問題呢？因為如果信號是稀疏的，則它是可壓縮的，也就是說里面那么多零，我只記錄那些非零值及它的位置就好了。

當然，現實中的信號本身一般並不是稀疏的，但經過一個變換后，在一組基上面是稀疏的，這就是信號的稀疏表示。

稀疏性是壓縮感知的前提。

2、范數||x||p

常見的有l0范數、l1范數、l2范數，經常要將l0范數等價為l1范數去求解，因為l1范數求解是一個凸優化問題，而l0范數求解是一個NP難問題，這些后面慢慢再說。

l0范數指的是x中非零元素的個數，即x的稀疏度，如果x是K稀疏的，則l0范數等於K；

l1范數指的是x中所有元素模值的和

l2范數指的是x中所有元素模值平方的和再平方，這個帶公式就可以了，它代表着距離的概念

還有無窮范數，指的是x中元素模的最大值

3、符號arg min

看壓縮感知的參考文獻最讓我難受的是很多數學符號都不認識，更難受的是還不知道這些符號從什么書里可以找到。

壓縮感知中常見如下表示：

s.t. 表示受約束於，是“subject to”的縮寫。

為了說明argmin的含義，可以參見Wikipedia中對argmax的解釋：

argmax : In mathematics, arg max stands for the argument of the maximum, that is to say, the set of points of the given argument for which the given function attains its maximum value.

舉三個例子自己體會一下就可以了：

argmin與其類似，琢磨一下就是了。

下面轉一段話：（max 和 argmax的區別）

y = f(t) 是一般常見的函數式，如果給定一個t值，f（t）函數式會賦一個值給y。
y = max f(t) 代表：y 是f(t)函式所有的值中最大的output。
y = arg max f(t) 代表：y 是f(t)函式中，會產生最大output的那個參數t。
看起來很模糊，舉個例子應該比較好理解：
假設有一個函式 f(t)，t 的可能范圍是 {0,1,2}，f(t=0) = 10 ; f(t=1) = 20 ; f(t=2) = 7，那分別對應的y如下：
y = max f(t) = 20
y= arg max f(t) = 1

這一塊要好好說一說，因為這是壓縮感知最基本的表示，是最常見的，但在不同的論文里面表示是不統一的：

a)焦李成，楊淑媛，劉芳，侯彪.壓縮感知回顧與展望[J].電子學報，2011，39(7)：1651-1662.