从稀疏表示到K-SVD，再到图像去噪
初学入门，将来有了更深的理解回来纠正错误。
MathJax需要加载，公式如未显示请稍候。
稀疏表示
稀疏性的理解
最初稀疏性产生于信号处理领域，自然界信号低频居多，高频主要是噪声，图像处理中的频率域滤波是个典型例子。
假设有一个干净没有噪声的图像，经过传输，收到的是一个受到干扰有了噪声的图像，而噪声主要是高频分量，对图片做二维傅里叶变换，对低频的波形保持，高频的一刀切，还原回来的图像就平滑了许多，大部分高频噪声就去除了。这个假设的场景就是个“信号恢复”的过程。如果把所有的频率的波都看作一个个相互正交的向量，恢复数据就是给这些向量找到一组系数，它们一乘、一合并，得到原始信号，频率从大到小是有无穷多个的，而由于自然界信号的“稀疏性”，对于图像而言，就是指有用的频率主要是低频，那么对应高频的系数基本都是0了，低频部分也不见得全是非0，这一系列系数0很多，非0很少，就很“稀疏”。
稀疏表示的概念
稀疏表示的目的就是在给定的超完备字典中用尽可能少的原子来表示信号。
意义在于降维，可以是压缩，可以用于机器学习特征提取，还有很多我也不知道的事情。。。
原子：信号的基本构成成分，比如一个长为N的列向量。
字典：许多原子的排序集合，一个N*T的矩阵，如果T>N，则为过完备或冗余字典。
咦，线性代数又出现了，假设列向量两两线性无关，
N*N的矩阵的秩就是N了，再增加向量也不会增加额外的信息。
稀疏表示的一个场景
假设现在有了一个N*T的过完备字典 D（比如前面所述图像傅里叶变换的所有频率的波），一个要表示的对象y（要还原的图像），求一套系数x，使得y=Dx，这里y是一个已知的长为N的列向量，x是一个未知的长为T的列向量，解方程。
这是一个T个未知数，N个方程的方程组，T>N，所以是有无穷多解的，线性代数中这样的方程很熟悉了。 \[
\begin{bmatrix}
1\\
2\\
3\\
4\\
5
\end{bmatrix} = \begin{bmatrix}
1 & 2 & 3 & 4 & 5 & 6 & 7 & 8\\
0 & 1 & 0 & 0 & 0 & 1 & 0 & 0\\
0 & 0 & 1 & 0 & 0 & 0 & 0 & 1\\
0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0
\end{bmatrix}\times
\begin{bmatrix}
x1\\
x2\\
x3\\
x4\\
x5\\
x6\\
x7\\
x8
\end{bmatrix}
\] 上面我就随便举了个N=5, T=8的例子，用来随便感受下。
这里可以引出一个名词，ill-posed problem（不适定问题），即有多个满足条件的解，无法判断哪个解更加合适，这是更“落地”的应用场景，inverse problem（逆问题），比如图像去噪，从噪声图中提取干净图。于是需要做一个约束。
增加限制条件，要求x尽可能稀疏，怎么“稀疏”呢？就是x的0尽可能多，即norm(x, 0)（零范数：非0元素个数）尽可能小。这样就有唯一解了吗？也还不是，如何能“约束”出各位合适的解，如何解，正是稀疏表示所研究的重点问题。比如后来有证明D满足一定条件情况下x满足norm(x,1)即可还原原始数据等，这有不少大神开启这个领域的故事这里就不讲了。
奥卡姆剃刀原理的思想：如果两个模型的解释力相同，选择较简洁的那个。稀疏表达就符合这一点。
针对这个例子我有个疑问，
x都比y还大了，这哪里压缩了。这个问题应该容易解答，例子里x是比y大，但是如果每个原子不是长度为N的列向量，而是个矩阵，或者更复杂的东西呢，x却依然只是一列系数。
求解
字典D已知，求y在过完备字典D上的稀疏表示x，被称作稀疏编码，模型是： \[
x = \mathop{\arg\min}_{x} norm(y - Dx, 2)^{2}, s.t. norm(x, 1) \le \varepsilon
\] 如何求解x？D又是怎么来的？先说D，D可以是前面说的傅里叶变换的一系列波啊，也可以是DCT的，也可以是小波的。但是科学家为了特定问题能有更具适应性的字典，让D也变成一个设计出来的量了，手工设计是不行，那么D也成了一个需要求的未知量。
已知D求x有OMP算法，大意是先找到D和y最接近的一个原子D(m)，求出合适的系数x(m)，新的y'=D(m) * x(m)，再找下一个最接近的原子，直到找完合适的x，如何确定最接近，如何计算x(m)，这里不再细说。
当任务是同时求出一个好的字典D，并得到一个满足稀疏约束的x，两步求解算法：先固定D，求个x出来，再固定x更新D，交替进行。
两步求解好熟悉的老套路。也许可以这么理解，求
x不再是个稀罕问题，而训练一个好的字典渐渐成为解决应用问题的关键，也是研究重点。
做这个求解的方法就有MOD、K-SVD等等一系列了。MOD分为两个步骤：Sparse Coding和Dictionary Update。
Sparse Coding: \[ x = \mathop{\arg\min}_{x} norm(y - Ax, 2)^{2}, s.t. norm(x, 1) \le k \]
Dictinary Update: \[ D = \mathop{\arg\min}_{x} norm(y - Ax, 2)^{2} \]
x的更新类似OMP，字典D的更新使用最小二乘法。
K-SVD
迭代K次，每次计算一下SVD分解的算法。SVD即奇异值，在了解SVD之前先复习一下矩阵的特征值。
特征值分解
特征值分解和奇异值分解是机器学习领域常见的方法。线性代数中我们熟悉的特征值 \(\lambda\) ，设 v是矩阵A的特征向量，则 \[
Av = \lambda v
\] v是 \(\lambda\) 对应的特征值。 矩阵的一组特征向量是相互正交的，特征值分解 将矩阵分解成如下形式： \[
A = Q\Sigma Q^{-1}
\]
其中Q是矩阵A的特征向量组成的矩阵， \(\Sigma\) 是一个对角阵，对角线上每个元素就是一个特征值，由大到小排列。
这像什么？\(\Sigma\) 里的一串 \(\lambda\) 就像前面
y=Dx里的x，而特征值矩阵Q就像字典D啊。特征向量的大小描述了每个特征值的权重，它们一起组合成了矩阵A，也就是那个y。
然而，特征值分解有个严重的局限——A必须是个方阵。
SVD（奇异值）分解
类比特征值分解，奇异值分解定义成这样： \[
A = U\Sigma V'
\] 假设A是一个N*M的矩阵，则U是N*N的方阵，\(\Sigma\) 是N*M的矩阵，V是M*M的方阵，于是奇异值分解就是求 \(\Sigma\)、U、V，V'是V的转置。于是套公式就可以求出奇异值、U、V，公式就不堆这里了。
说来也巧，我之前在一个量子计算的书里看过SVD分解，当时还做了一下习题手算了半天SVD，这部分理解起来舒服多了
SVD怎么和稀疏性搭边呢？因为奇异值矩阵 \(\Sigma\) （虽然不是方阵，但也是按45°角放着一串奇异值的类似对角阵的东西）的“对角线”上大部分数值是0或接近0的，类比特征值分解，这些奇异值就是“权重”嘛，如果把接近0的这些丢掉，是不是清（稀）爽（疏）很多？部分奇异值分解，A还是N*M，但假设取比较大的r个奇异值， \(\Sigma\) 变成r*r，部分SVD分解即：
\[ A_{N\times M} \approx U_{N\times r} \Sigma_{r\times r} V^{-1}_{r\times M} \]
如果r很小，而等式左右又能很接近，那数据就被压缩的相当不错，保存U、\(\Sigma\) 、V 要比保存A本身节省空间多了。
K-SVD字典学习
K-SVD和MOD最大的不同在于，每次只更新字典的一个原子（即D的一列），而不是每次用一个x更新整个D。
回忆下前面的y=Dx，但是学一个字典，当然不能只用一个数据，现在来升级版： \[
Y=DX
\] 哈？小写变大写？意思是一组y和其对应的一组x，那么Y和X指的矩阵。
\[ \begin{bmatrix} y_{11} & y_{21} & y_{31} & y_{41} \\ y_{12} & y_{22} & y_{32} & y_{42} \\ y_{13} & y_{23} & y_{33} & y_{43} \\ y_{14} & y_{24} & y_{34} & y_{44} \\ y_{15} & y_{25} & y_{35} & y_{45} \end{bmatrix} = \begin{bmatrix} d_{11} & d_{21} & d_{31} & d_{41} & d_{51} & d_{61} & d_{71} & d_{81}\\ d_{12} & d_{22} & d_{32} & d_{42} & d_{52} & d_{62} & d_{72} & d_{82}\\ d_{13} & d_{23} & d_{33} & d_{43} & d_{53} & d_{63} & d_{73} & d_{83}\\ d_{14} & d_{24} & d_{34} & d_{44} & d_{54} & d_{64} & d_{74} & d_{84}\\ d_{15} & d_{25} & d_{35} & d_{45} & d_{55} & d_{65} & d_{75} & d_{85} \end{bmatrix}\times \begin{bmatrix} x_{11} & x_{21} & x_{31} & x_{41}\\ x_{12} & x_{22} & x_{32} & x_{42}\\ x_{13} & x_{23} & x_{33} & x_{43}\\ x_{14} & x_{24} & x_{34} & x_{44}\\ x_{15} & x_{25} & x_{35} & x_{45}\\ x_{16} & x_{26} & x_{36} & x_{46}\\ x_{17} & x_{27} & x_{37} & x_{47}\\ x_{18} & x_{28} & x_{38} & x_{48} \end{bmatrix} \]
现在要更新字典D的第k个原子，也就是第k列，它能影响到的是Y的第k行，同样对应D的第k列的系数，也是X的第k行。
目标函数的转化：
\[
\begin{align*}
& \quad\ \left \|Y - DX\right \|^{2}_{F} \\
& =\left \|Y-\Sigma_{j=1}^{K}d_{j}x^{T}_{j}\right \|^{2}_{F} \\
& =\left \|(Y-\Sigma_{j\neq k}d_{j}x^{T}_{j})-d_{k}x^{T}_{k}\right \|^{2}_{F} \\
& =\left \|E_{k}-d_{k}x^{T}_{k}\right \|^{2}_{F}
\end{align*}
\] \(E_{k}\) 是去掉原子 \(d_{k}\) 的 \(D\) 中的误差，于是目标函数转化为 \(D\) 的其他列固定，要更新的 \(d_{k}\) 使全局误差（ \(\left \|Y-DX\right \|\) ）最小化。 即可得到字典的第k个原子。求解这里的 \(D_{k}, x^{T}_{k}\) ，就用到对 \(E\) 的SVD分解了。
但是直接分解 \(E\) 得到的 \(x^{T}_{k}\) 并不稀疏。
更新字典和稀疏系数是迭代进行的，在“本次”迭代中，找到“上次”迭代中哪些Y用到了字典D的原子k，也就是X的第k行哪些元素不为0，\(x_{1k},x_{2k},x_{3k},x_{4k}\) 里，假设 \(x_{1k},x_{3k}\) 不为0，那么对应的Y的1,3列就是用到了D的原子k的信号（Y的每列是一个信号）。现在把它们拆出来： \[
\begin{align*}
& Y^{temp}_{k} = \begin{bmatrix}
y_{11} & y_{31} \\
y_{12} & y_{32} \\
y_{13} & y_{33} \\
y_{14} & y_{34} \\
y_{15} & y_{35}
\end{bmatrix} \\
& X^{temp}_{k} = \begin{bmatrix}
x_{1k} & x_{3k}
\end{bmatrix} \\
& D^{temp}_{k} = \begin{bmatrix}
d_{k1}\\
d_{k2}\\
d_{k3}\\
d_{k4}\\
d_{k5}
\end{bmatrix} \\
\end{align*}
\]
这样得到只保留非零位置的X、D计算目标函数后得到的只保留对应位置的 \(E^{temp}_{k}\) ，对这个 \(E^{temp}_{k}\) 再做SVD分解，\(E^{temp}_{k} = U \Sigma V^{T}\)， \(U\) 的第一列即为新的 \(\widetilde{d}_{k}\)， \(V\) 的第一列与 \(\Sigma(1, 1)\) 的乘积为新的 \(\widetilde{x}^{T}_{k}\) 。
逐列更新得到新字典 \(\widetilde{D}\) 。
K-SVD图像去噪
上面提到过稀疏表示基本的目标函数是： \[ x = \mathop{\arg\min}_{x} \left \|y - Dx\right \|^{2}_{2}, s.t. \left \|x \right \|_{0} \le \varepsilon \]
这里暂时用0范数来说明。
假设现在有一个零均值高斯白噪声，即 \(n\sim N(0,\sigma)\) ， \(\sigma\) 是噪声的标准差，有噪声的图像为 \(z=y+n\) ，目的是从信号 z 中恢复出原始无噪信号 y，通过最大后验概率，求得目标函数的解，即可恢复出y：
\[
x = \mathop{\arg\min}_{x} \left \|z - Dx\right \|^{2}_{2}, s.t. \left \|x \right \|_{0} \le T
\] > 这里最大后验概率能恢复y是为什么暂时没管，肯定有证明了。。
其中T依赖于 \(\varepsilon\) 和 \(\sigma\) 。为方便优化计算，实际操作中往往转化成： \[
x = \mathop{\arg\min}_{x} \left \|z - Dx\right \|^{2}_{2}+\mu \left \|x \right \|_{0}
\] 选取恰当的\(\mu\)可以让上面两式等价。
优化属于NP难问题，有一系列研究，可通过OMP、BP、FOCUSS等算法来获得近似解。如果 x 足够稀疏，近似解就足够接近精确解。