稀疏表示与字典学习

参考：
1.稀疏表示：https://www.cnblogs.com/yifdu25/p/8128028.html
2.K-SVD : https://blog.csdn.net/chlele0105/article/details/16886795
3.OMP : https://www.cnblogs.com/yifdu25/p/8385204.html
4.https://www.cnblogs.com/yifdu25/p/8128028.html
5.MP的实例 https://korediantousman.staff.telkomuniversity.ac.id/files/2017/08/main.pdf
6.K-SVD https://www.cnblogs.com/salan668/p/3555871.html
7.OMP与MP算法： https://blog.csdn.net/scucj/article/details/7467955

1. 稀疏表示

稀疏表示(sparse representation)假设自然信号可以被一个字典的原子(基)稀疏的线性组合近似表示。即y≈Dα，这里$y\in R^n$表示信号，$\alpha \in R^K$表示稀疏系数，$D\in R^{n\times K}$表示是过完备字典。向量α只具有很少的非0值，这便是稀疏的含义。之所以需要字典，是因为一般的自然信号x本身并不是稀疏的，往往需要在某种稀疏基上才可以进行稀疏表示。所以稀疏表示的关键是在某种基(字典)下的系数稀疏，而这些字典与图像本身密切相关，可以将字典看做是组成图像的基本信息。也就是说在稀疏表示中我们把字典当成一种变换域，在该变换中信号的表示是稀疏的。该问题用公式表述为：

其中$Y\in R^{n\times M}$表示样本矩阵，每一列代表一个patch，M个样本， D = { d 1 , d 2 , . . . d K } ∈ R n × K D=\{d_{1},d_{2},...d_{K}\}∈R^{n×K} D={d1​,d2​,...dK​}∈Rn×K表示字典，每一列代表一个原子，K个原子，$X\in R^{K\times M}$表示系数矩阵，每一列$x_i$对应一个样本在该字典上的系数，ε表示重构误差。对于表示一个信号使用字典稀疏表示的情况如下。

在这里，稀疏性表现在$x_i$的系数大部分都是0，而只有很少的非0值。
由于求解0范数是非凸的问题，而在一定条件下，0范数问题可以转化为1范数求解，所以(1)式也可以转化为：

1.1 字典的生成方式

通常有2种生成字典的方式，一种基于分析模型(pre-constructed)，比如利用傅立叶基，小波基(WT)，离散余弦变换(DCT)基，Gabor基，Contourlet基，等等。依靠这种方法生成的字典结构良好，数值计算快速。但是这种字典具有局限性，只能表示某种类型的信号，过度依赖于图像的几何特征，并且只有少数图像块可以用稀疏的原子组合来表示。另一种基于学习模型，将图像分成一个个小patch，通过训练样本得到字典。由这种方法得到的字典具有很好的适应性，能够捕捉图像的几何特征。相较于第一种往往表现更好。

1.2 字典的种类

而从字典的形式上看，通俗的说就是从字典的长和宽的对比，字典又可以分成3种：过完备字典(n<K)，即原子个数要远远大于信号的维度，这种情况在稀疏表示中最常见，完备字典(n=K)例如傅里叶变换和DCT变换都是这种情况和欠完备字典(n>K)[4]。通常我们使用的都是过完备字典。那么为什么通常都要求过完备字典呢？参考别人的说法，如果n<K，这个方程要么有无穷多个解（当增广矩阵[D,y]的秩=矩阵[D]的秩时），要么无解（当增广矩阵[D,y]的秩>矩阵[D]的秩时）。可能在方程有无穷多解的时候更容易找到更适合的稀疏系数。这大概就是过完备性字典的优点。

1.3 稀疏表示的相关算法

从上面的讲解来看稀疏表示其实就是在求解式(1)，但由于其中只有Y是已知的，而D和X都是未知的，这里就涉及到2个问题，1)如何对信号进行稀疏表示，即如何求稀疏系数$x_i$？2)如何根据样本生成字典，即如何求解D？概括地说，我们可以使用任何追踪算法来求解稀疏系数$x_i$，例如MP，OMP(正交匹配追踪算法)，BasisPursuit (BP)，FocalUnderdetermined System Solver (FOCUSS)。而更新字典最常用的就是K-SVD算法。

2. 如何稀疏表示

2.1 MP算法

这里借用[3]的博客的内容来讲解MP算法。假设有这样一个问题：要使用OB,OC,OD来表示OA。这里OB,OC,OD就相当于字典，OA相当于要稀疏表示的信号。那么应该如何求解才能达到使用最少的向量近似表示OA呢？答案是每次将OA投影到各个向量上，选择OA投影最大的那个向量表示原信号会在当前情况下保存信号最多的信息。
1)如图1，将OA向OB,OC,OD分别进行投影，发现OA在OB上的投影最大，这样OA就被分成投影向量OM和残差向量MA。假设$\frac{OM}{OB}=a$，那么OA=a OB+MA；2)然后继续对残差MA进行分解，和第一步一样将MA向各个向量上投影，如图2，发现MA在OC向量上的投影最大，则选择OC来表示MA，并且假设$\frac{ON}{OC}=b$，此时$MA=b\ast OC+N{A}^{{}^{\prime }}$，3)假设此时会有2种情况， [1]发现残差仍然大于阈值ε，那么就按照第一步继续分解，直到残差小于阈值。[2]发现残差已经小于阈值了，那么此时$OA=a\ast OB+b\ast OC+N{A}^{{}^{\prime }}$。
以上就是追踪算法MP的思想。稀疏表示记为OA ≈ aOB +b*OC 。a,b就相当于系数向量里面的元素。可以看出如果残差值可以忽略，则信号OA就是这些原子的线性组合。
那如何使用代码来实现投影最大呢？也就是如何用代码寻找最接近原子？又该如何计算残差呢？

1. 选择最接近残差的原子：MP里用向量内积定义原子与残差的距离。也就是残差与某原子内积最大，即表示残差在该原子上的投影最大。用R表示残差，$d_i$表示原子，则每次都在寻找$Max[Dist(R,d_j)]$，假设此时最大的原子下标为$d_i$。
2. 残差更新：$R=R-<R,d_i>d_i$;继续选择下一个原子，直至收敛。
   使用二维空间上的向量表示该过程。
   
   红色向量r代表当前残差，绿色向量则表示r投影到$d_3$之后的残差更新。从上图中也可以看出当前残差r是直角三角形的斜边，而新的残差绿色向量则是一条直角边。即$|r_{k+1}|<|r_k|$，所以该算法是收敛的。
   值得注意的有2点，1) MP算法使用的各个原子是归一化的，即$||d_i|{|}^2=1$。因为选取$max<R,d_i>$时，如果$d_i$长度不统一，不能得出最好的投影。2) MP算法不是最优的，得到的解是次优解。举个例子，如果字典中只有两个向量d1,d2,那么MP算法会在这两个向量间交叉迭代投影，也就是f=a1d1+a2d2+a3d1+a4d2+……;可以看到之前投影过的原子方向 $d_1$，之后还有可能投影。换句话说，MP的方向选择不是最优的，是次优的。
   理论上说，假设MP算法第k次迭代的结果为$y={\sum }_{i=1}^k{a}_i{d}_i+R_{k}f=f_k+R_{k}f$。由于MP算法仅能保证$R_{k}f\perp d_k$，所以$f_k$一般情况下是次优的。这是什么意思呢？$f_k$是字典中k个原子对y的线性表示，这个组合的值作为近似值，只有在第k个残差和$f_k$里所有原子都正交，才是最优的。如果第k个残差与$f_k$正交，意味着第k个残差与$f_k$的任意一项都线性无关，那么其在后面的分解中，不可能出现$f_k$中已经出现的原子，这才是最优的。而一般情况下，MP不能满足这个条件，它一般只能满足第k个残差和$d_k$正交。也正是因为MP的这个特点，所以MP算法需要更多次迭代才能收敛。

MP算法流程

算法流程总结如下，需要说明的是这里使用x表示信号，α表示系数向量。与上面的符号表示略有差异，主要是因为直接摘抄别人的东西。

2.2 OMP算法

OMP算法即正交的MP算法。MP算法的次最优性来源其残差只与当前投影方向垂直，这样在接下来的投影中，因为残差与已选到的原子不是线性无关，所以残差很有可能会再次投影到原来的方向。于是，在投影时，如果使得残差$R_{k+1}$与$x_1到x_{k+1}$(已选到的原子)的所有向量垂直，则可以克服这个问题。这正是OMP算法改进的地方，OMP算法在分解的每一步对所选择的全部原子进行正交化处理，这使得在精度要求相同的情况下，OMP算法的收敛速度更快。更多信息请参考https://blog.csdn.net/scucj/article/details/7467955

3. 如何生成字典–K-SVD

在进行稀疏表示的时候，我们往往会先初始化一个字典，但是它并不是最优的。常用K-SVD算法对字典进行更新。但是K-SVD算法并不是独立的，它需要和上面讲述的追踪算法结合在一起使用才能对字典进行更新。需要注意的是，这里

3.1 稀疏表示

首先要初始化一个字典D，使用该字典对数据进行稀疏表示，得到系数矩阵X。然后把DX看做是D中每列$d_i$与X中每行$x_i$(在这里表示X的一行，前面该符号表示X的一列)的乘积，即将DX的结果分片。如前所述 D = { d 1 , d 2 , . . . d K } ∈ R n × K D=\{d_{1},d_{2},...d_{K}\}∈R^{n×K} D={d1​,d2​,...dK​}∈Rn×K表示字典，每一列$d_i\in R^{n\times 1}$代表一个原子，有K个原子，$X\in R^{K\times M}$表示系数矩阵，每一行$x_i\in R^{1\times M}$。$d_i$与$x_i$相乘的结果也是一个n×M的矩阵。表示该原子在整个重构矩阵中的贡献。DX分片的结果如下，

3.2 字典更新

由于初始的字典往往并不是最优的，所以需要对初始字典进行更新。而更新的动力就来源于残差矩阵。K-SVD采用逐列更新字典的方法，通过K次迭代完成一次字典更新。之所以需要K次主要是因为字典D有K个原子。而小写k表示每次更新字典所正在操作的原子。具体来说每次都提取出一个原子$d_k$对整个重构矩阵DX的影响，因为缺少了$d_k$，上述表达式会产生一个"空洞"。K-SVD算法的目的就是如何选择一个新的原子填补空洞，使得结果更加逼近Y，减小整体误差。提取出一个原子$d_k$的影响之后，当前误差矩阵如下。

误差值为$$E_n=||E|{|}_F^2$$
正如3.1所述，在开始更新字典之前就已经得到了字典D和其对应的系数矩阵X。言下之意，在更新字典的时候，D,X都是固定的，但因为重构的结果不好所以要更新。假设此时将要更新的是系数矩阵的$x_k$行和字典的$d_k$列，则有，

上式中的$E_k$是误差矩阵，对$E_k$做SVD分解，得到$E_k=U\Lambda V^T$，其中U和V的列矢量均是正交基，Λ是对角矩阵。若Λ的对角元素从大到小排列，则表示$E_k$的能量分量主轴在相应几个正交方向上由大到小分配，如此取U的第一个列向量来表示$d_i$，取V的第一个列向量与Λ的第一个元素的乘积表示$x_i$，这样就完成了字典一个条目的更新。
值得注意的是，如果直接用上面的方法来更新$x_k$和$d_k$，则会导致$x_k$不稀疏，出现"发散"。换句话说，更新之后的$x_{knew}$与更新$d_k$前的$x_k$的非零元素所处位置和value不一样。
那应该如何更新$x_k$和$d_k$？
处理方法是在进行SVD分解之前，我们使用如下公式对$E_k$和$x_k^T$进行变换。即只保留系数中的非零值，再进行SVD分解就不会出现这种现象了。

其中$\Omega \in R^{M\times W}$，W是$x_k^T$中的非零元素个数，M是样本个数。其实质就是保留$x_k^T$中的非零元素，$E_k$中只保留$d_k$和$x_k^T$中非零元素乘积后的那些项的贡献，形成$E_k^R$。换句话说，乘以新的矩阵Ω之后的结果就是$E_k$,$x_k^T$去掉0输入之后的收缩结果。
那如何生成矩阵Ω？定义集合 w k = { i ∣ 1 ≤ i ≤ M , x k T ( i ) ≠ 0 } w_{k}=\{i|1≤i≤M,x_{k}^{T}(i)≠0\} wk​={i∣1≤i≤M,xkT​(i)̸​=0}表示$x_k^T(i)̸=0$的点的索引。定义Ω为$M\times |w_k|$矩阵，它在$(w_k(i),i)$处值为1，其余都是0。
此时对$E_k^R$做SVD分解得到$E_k^R=U\Delta V^T$，按照上面的思路使用U的第一列更新$d_k$，使用V的第一列和Δ(1,1)的乘积更新$x_k^T$，这样更新字典就正确了。

K-SVD算法流程

总得来说，整个字典学习的流程如下。它分为2个阶段，1)对数据使用OMP算法(也可以使用其他算法)进行稀疏表示，得到对应字典的系数矩阵X。2) 依据得到系数矩阵X，根据K-SVD算法对字典逐列进行更新，同时也对系数矩阵进行更新，经过K次迭代完成一次字典更新。直到算法收敛。下面是K-SVD算法的流程。