[综] 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又是另外一回事了，慢慢之后再说吧。