K-SVD算法总结

这几天看了稀疏表示的一些文章，对字典学习方法K-SVD[1]查阅了相关资料，特此总结如下，如有理解上不正确的地方，还望指正，本人还处于初学者的状态。一、概述K-SVD是一种迭代算法，是K-means算法的扩展，一般是用来在稀疏表示问题中的字典训练方面。这里的“字典”是一个过完备的矩阵，由其，使得一个信号向量可以表示成字典中原子（字典的列向量）的稀疏线性组合。 K-SVD和K-means方法本质上都属于一种压缩的思想，都主要包含以下两个步骤： 1)稀疏编码2)字典更新在K-means方法中，K-means 先随机选择K个初始点作为字典，K个初始点就代表K类。然后在每次迭代过程中，稀疏编码阶段：计算数据集中每个数据点与这K个点的距离，距离最短则代表改点属于该类；字典更新：每一类中所有点的均值作为新的字典。 而在K-SVD中，稀疏编码可以采用任何基方法（MP、OMP、BP）；字典更新采用SVD奇异值分解。文章原文引用如下：The K-SVD algorithm is an iterative method that alternates between sparse coding of theexamples based on the current dictionary and an update process for the dictionary atoms so as tobetter fit the data, generalizing the K-means algorithm. The update of the dictionary columns is donejointly with an update the sparse representation coefficients related to it, resulting in acceleratedconvergence.二、K-SVD方法6.5（双系统）it+CUDA6.5配置说明这里给出文章中对K-SVD方法的描述 以下简要描述该算法： 用数学表达式表示稀疏表示问题，K-SVD即是求解满足以下式子的字典D 其中是信号 ,  是过完备字典矩阵,  为系数矩阵 。 1、初始化操作：初始化字典矩阵D，D中原子(列向量)必须是单位向量，2、重复以下步骤直到满足收敛条件：1）稀疏编码：这里使用OMP(正交匹配基)来得到每一个信号样本对应的稀疏系数向量的近似逼近解。 - OMP算法 本质思想：以贪婪迭代的方法选择字典原子，使得每次迭代的过程中原子与信号最大相关（内积最大），从信号向量减去相关部分得到残差，残差按照上述方法反复迭代，直到迭代次数达到稀疏度K 或者残差小于误差阈值，停止迭代。 OMP算法步骤： 输入：字典矩阵D, 采样信号 , 稀疏度K 输出：稀疏系数  的 K-稀疏的逼近() 初始化：残差 , 索引集 , 循环执行步骤1-5： 步骤1：找出残差  和字典内积最大的原子索引，即 步骤2：更新索引集，记录字典的重建原子集合步骤3：由最小二乘得到。OMP要求在每次迭代过程中残差与之前所选择的原子都正交，正交是通过最小二乘实现（见下图，b=y, A=D, y、Dx、残差组成正交三角形，避免了重复选择原子的现象，每一步都使残差最小，同时加快收敛速度) 6.5（双系统）it+CUDA6.5配置说明{}subjectto∀i,≤minD,X∥Y−DX∥2F∥∥xi0T0Y(n×N)D(n×K)X(K×N)t=1yixiyxx^=yr0=∅Λ0t=1rλ=arg⟨,⟩λimaxj=1,...,N∣∣ri−1dj∣∣=⋃{}ΛtΛt−1λt=[,]DtDt−1dλt=argminx^t∥y−∥Dtx^2步骤4：更新残差步骤5：判断是否满足, 若满足，则停止迭代；若不满足，则执行步骤1.这里给出OMP的Matlab代码，以助于理解rt=y−Dtx^t，t=t+1 tK2）字典更新：采用逐次更新字典原子向量，通过K次迭代完成字典的一次更新。具体如下： 首先 将 DX表示成D中原子（列）与X中每行的乘积, 然后剥离第k个原子，上述表达式会产生一个“空洞”，字典更新的目的就是寻找新的  和  来填补这个“空洞”来更加趋于收敛情况，所用方法就是SVD. 奇异值分解SVD是为了提取出一个矩阵最重要的特征, 适用于任何普通矩阵. function[A]=OMP(D,X,L);%=============================================%Sparsecodingofagroupofsignalsbasedonagiven%dictionaryandspecifiednumberofatomstouse.%inputarguments:%D-thedictionary(itscolumnsMUSTbenormalized).%X-thesignalstorepresent%L-themax.numberofcoefficientsforeachsignal.%outputarguments:%A-sparsecoefficientmatrix.%=============================================[n,P]=size(X);[n,K]=size(D);fork=1:1:P,a=[];x=X(:,k);residual=x;indx=zeros(L,1);forj=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;ifsum(residual.^2)1e-6break;endend;temp=zeros(K,1);temp(indx(1:j))=a;A(:,k)=sparse(temp);end;return;1234567891011121314151617181920212223242526272829303132333412345678910111213141516171819202122232425262728293031323334DX=∑Ki=1dixTidixi的列矢量均是正交基，Λ 是对角矩阵。若Λ的对角元素从大到小排列，则表示Ek的能量分量主轴在相应几个正交方向上由大到小分配，如此我们取U的第一个列向量来表示, 取V的第一个列向量与,这样就完成了字典一个原子的更新。 然而上述操作会引发一个问题，X 本是一个稀疏矩阵，而上述操作得到的xk 不一定稀疏。解决方法是将 中非， as a matrix of size N×|ωk| , with ones on theth entries and zeros elsewhere.) 然后进行下式变换： 对 进行SVD分解即可。式子中，是为了仅保留字典中有用原子对应的误差矩阵中有用分量，从EkΩk 而避免直接SVD分解不稀疏的情况。上式中的 EK 是误差矩阵, 对其做SVD分解，Ek=UΛVT , 其中U和VdkΛ 的第一个元素的乘积表示xkxk零元素构建一个新的矩阵Ω(N×|wk|) |ωk| 为xk 中非零元素的个数。(Define Ωk (wk(i),i)ERkER=k