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

稀疏表示 分爲 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

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

    end

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

題記：從今年九月份開始看壓縮感知方面的文獻，本身的感受是要求的數學知識太多了，上研時只學了一個矩陣理論，學的感受還算紮實，但學完到如今基本都忘掉了，什麼優化之類的根據就沒學過，因此如今看文獻真心很費勁，不過不會能夠慢慢學嘛，之後我會根據本身的學習狀況陸續發一些壓縮感知經常使用到的數學知識。

本次說三個問題：

一、稀疏

二、範數

三、符號arg min

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

 

=========================如下爲正文=========================

一、稀疏：什麼是K稀疏呢？

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

爲何要談稀疏這個問題呢？由於若是信號是稀疏的，則它是可壓縮的，也就是說裏面那麼多零，我只記錄那些非零值及它的位置就行了。

固然，現實中的信號自己通常並非稀疏的，但通過一個變換後，在一組基上面是稀疏的，這就是信號的稀疏表示。

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

 

二、範數||x||p

常見的有l0範數、l1範數、l2範數，常常要將l0範數等價爲l1範數去求解，由於l1範數求解是一個凸優化問題，而l0範數求解是一個NP難問題，這些後面慢慢再說。

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

l1範數指的是x中全部元素模值的和

l2範數指的是x中全部元素模值平方的和 再平方，這個帶公式就能夠了，它表明着距離的概念

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

 

三、符號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.

     

 

b)石光明，劉丹華，高大化，劉哲，林傑，王良君.壓縮感知理論及其進展[J].電子學報，2009，37(5)：1070-1081.

c)楊海蓉，張成，丁大爲，韋穗.壓縮傳感理論與重構算法[J].電子學報，2011，39(1)：142-148.

在壓縮感知理論方面，無論是用min仍是argmin（文獻ab與文獻c區別），無論min下面有沒有變量（文獻a與文獻b區別），其實表達的意思都是同樣的：

若是用0範數，則是求得知足後面約束條件的最稀疏的x(或θ)；

若是用1範數，則是求得知足後面約束條件的元素模值和最小的x(或θ)；

固然兩種求法在知足必定條件下(RIP)是等價的，RIP又是另外一回事了，慢慢之後再說吧。