KSVD去噪

Posted 2023-04-22

参考技术A       在解释KSVD去噪原理之前先解释几个名词，首先：

       原子：信号的基本构成成分，比如一个长为N的列向量；

       字典：许多原子的排序集合，一个N*T的矩阵，如果T>N（列数大于行数），则为过完备或冗余字典。信号稀疏之前在压缩感知学习的时候有介绍过，就是信号的非零值很少，这个概念主要应用与信号处理领域，自然信号中主要是低频信息，高频信息大多就是噪声了，在图像中应用滤波器滤除高频成分也就是去噪了。具体怎么完成，就要依靠下面介绍的KSVD。

       KSVD的目的就是找到最稀疏的稀疏矩阵X，使得原始矩阵（Y）有最稀疏的表示。假设现在有了一个N*T的过完备字典 D，一个要表示的对象y（要还原的图像），求一套系数x，使得y=Dx，这里y是一个已知的长为N的列向量，x是一个未知的长为T的列向量，解方程。这是一个T个未知数，N个方程的方程组，T>N，所以是有无穷多解，但是针对问题目标我们会给这个方程添加约束条件，比如在图像去噪中，我们想要x最稀疏，就是非零值最少，这个在已知D和y求解x的过程就是稀疏编码。

       稀疏编码求解的模型就是：x = argminxnorm(y − Dx, 2)2,  s.t.norm(x, 1)≤ε。这里又被分为了两种可能，D已知情况下可以用OMP算法（大意是先找到D和y最接近的一个原子D(m)，求出合适的系数x(m)，新的y'=D(m) * x(m)，再找下一个最接近的原子，直到找完合适的x）；当D未知时就变成了矩阵分解问题，比如MOD算法的路子：Sparse Coding和Dictionary Update，两步走，第一步固定D,更新x：x = argminxnorm(y − Ax, 2)2, s.t.norm(x, 1)≤k；第二步更新D: D = argminxnorm(y − Ax, 2)2很像EM算法。

       KSVD和MOD最大的不同就是K每次只更新D当中的一个原子，就是D矩阵当中的某一列，因为矩阵相乘可以看做是前一个矩阵的列向量分别乘后一个矩阵的行向量。Loss函数在这里可以进行转化：

      而SVD就用在对E的分解，和上面的MOD类似，也是迭代进行就是每次更改D的原子。就比如说KSVD用于图像去噪的话，假如有一个零均值高斯白噪声，即 n ∼ N(0, σ) ，σ是噪声的标准差，有噪声的图像为 z = y + n ，目的是从信号 z 中恢复出原始无噪信号 y，通过最大后验概率，求得目标函数的解，即可恢复出y：x = argminx∥z−Dx∥22,   s.t.∥x∥0 ≤ T 。其中T依赖于 ε 和 σ 。为方便优化计算，实际操作中往往转化成： x = argminx∥z−Dx∥22+ μ∥x∥0选取恰当的μ可以让上面两式等价。

字典学习（Dictionary Learning, KSVD）详解

注：字典学习也是一种数据降维的方法，这里我用到SVD的知识，对SVD不太理解的地方，可以看看这篇博客：《SVD（奇异值分解）小结 》。

1、字典学习思想

字典学习的思想应该源来实际生活中的字典的概念。字典是前辈们学习总结的精华，当我们需要学习新的知识的时候，不必与先辈们一样去学习先辈们所有学习过的知识，我们可以参考先辈们给我们总结的字典，通过查阅这些字典，我们可以大致学会到这些知识。

为了将上述过程用准确的数学语言描述出来，我们需要将“总结字典”、“查阅字典”做出一个更为准确的描述。就从我们的常识出发：

1. 我们通常会要求的我们的字典尽可能全面，也就是说总结出的字典不能漏下关键的知识点。
2. 查字典的时候，我们想要我们查字典的过程尽可能简洁，迅速，准确。即，查字典要快、准、狠。
3. 查到的结果，要尽可能地还原出原来知识。当然，如果要完全还原出来，那么这个字典和查字典的方法会变得非常复杂，所以我们只需要尽可能地还原出原知识点即可。

注： 以上内容，完全是自己的理解，如有不当之处，欢迎各位拍砖。

下面，我们要讨论的就是如何将上述问题抽象成一个数学问题，并解决这个问题。

2、字典学习数学模型

2.1 数学描述

我们将上面的所提到的关键点用几个数学符号表示一下：

• “以前的知识”，更专业一点，我们称之为原始样本，用矩阵(mathbf{Y})表示；
• “字典”，我们称之为字典矩阵，用(mathbf{D})表示，“字典”中的词条，我们称之为原子（atom），用列向量(mathbf{d}_k)表示；
• “查字典的方法”，我们称为稀疏矩阵，用(mathbf{X})；
• “查字典的过程”，我们可以用矩阵的乘法来表示，即(mathbf{DX})。

用数学语言描述，字典学习的主要思想是，利用包含(K)个原子(mathbf{d}_k)的字典矩阵(mathbf{D}in mathbf{R}^{m imes K})，稀疏线性表示原始样本(mathbf{Y} in mathbf{R}^{m imes n})（其中(m)表示样本数，(n)表示样本的属性），即有(mathbf{Y=DX})（这只是我们理想的情况），其中(mathbf{X} in mathbf{R}^{K imes n})为稀疏矩阵，可以将上述问题用数学语言描述为如下优化问题

[ min_{mathbf{D, X}}{|mathbf{Y}-mathbf{DX}|^2_F},quad ext{s.t.} forall i, |mathbf{x}_i|_0 le T_0 ag{2-1} ]

或者

[ min_{mathbf{D, X}}sum_i|mathbf{x}_i|_0, quad ext{s.t.} min_{mathbf{D, X}}{|mathbf{Y}-mathbf{DX}|^2_F} le epsilon, ag{2-2} ]

上式中(mathbf{X})为稀疏编码的矩阵，(mathbf{x}_i, (i=1,2,cdots,K))为该矩阵中的行向量，代表字典矩阵的系数。

注： (|mathbf{x}_i|_0)表示零阶范数，它表示向量中不为0的数的个数。

2.2 求解问题

式（2-1）的目标函数表示，我们要最小化查完的字典与原始样本的误差，即要尽可能还原出原始样本；它的限的制条件(|mathbf{x}_i|_0 le T_0)，表示查字典的方式要尽可能简单，即(mathbf{X})要尽可能稀疏。式（2-2）同理。

式（2-1）或式（2-2）是一个带有约束的优化问题，可以利用拉格朗日乘子法将其转化为无约束优化问题

[ min_{mathbf{D, X}}{|mathbf{Y}-mathbf{DX}|^2_F}+lambda|mathbf{x}_i|_1 ag{2-3} ]

注： 我们将(|mathbf{x}_i|_0)用(|mathbf{x}_i|_1)代替，主要是(|mathbf{x}_i|_1)更加便于求解。

这里有两个优化变量(mathbf{D, X})，为解决这个优化问题，一般是固定其中一个优化变量，优化另一个变量，如此交替进行。式（2-3）中的稀疏矩阵(mathbf{X})可以利用已有经典算法求解，如Lasso（Least Absolute Shrinkage and Selection Operator）、OMP（Orthogonal Matching Pursuit），这里我重点讲述如何更新字典(mathbf{D})，对更新(mathbf{X})不多做讨论。

假设(mathbf{X})是已知的，我们逐列更新字典。下面我们仅更新字典的第(k)列，记(mathbf{d}_k)为字典(mathbf{D})的第(k)列向量，记(mathbf{x}^k_T)为稀疏矩阵(mathbf{X})的第(k)行向量，那么对式（2-1），我们有

[ egin{eqnarray} {|mathbf{Y}-mathbf{DX}|^2_F} &=&left|mathbf{Y}-sum^K_{j=1}mathbf{d}_jmathbf{x}^j_T ight|^2_F &=&left|left(mathbf{Y}-sum_{j e k}mathbf{d}_jmathbf{x}^j_T ight)-mathbf{d}_kmathbf{x}^k_T ight|^2_F &=&left|mathbf{E}_k - mathbf{d}_kmathbf{x}_T^k ight|^2_F end{eqnarray} ag{2-4} ]

上式中残差(mathbf{E}_k=mathbf{Y}-sum_{j e k}mathbf{d}_jmathbf{x}^j_T)，

此时优化问题可描述为

[ min_{mathbf{d}_k, mathbf{x}^k_T}left|mathbf{E}_k - mathbf{d}_kmathbf{x}_T^k ight|^2_F ]

因此我们需要求出最优的(mathbf{d}_k, mathbf{x}_T^k)，这是一个最小二乘问题，可以利用最小二乘的方法求解，或者可以利用SVD进行求解，这里利用SVD的方式求解出两个优化变量。

但是，在这里我人需要注意的是，不能直接利用(mathbf{E}_k)进行求解，否则求得的新的(mathbf{x}_k^T)不稀疏。因此我们需要将(mathbf{E}_k)中对应的(mathbf{x}_T^k)不为0的位置提取出来，得到新的(mathbf{E}_k^{'})，这个过程如图2-1所示，这样描述更加清晰。

图2-1 提取部分残差

如上图，假设我们要更新第0列原子，我们将(mathbf{x}_T^k)中为零的位置找出来，然后把(mathbf{E}_k)对应的位置删除，得到(mathbf{E}_k^{'})，此时优化问题可描述为

[ min_{mathbf{d}_k, mathbf{x}^k_T}left|mathbf{E}_k^{'} - mathbf{d}_kmathbf{x{'}}_T^{k} ight|^2_F ag{2-5} ]

因此我们需要求出最优的(mathbf{d}_k, mathbf{x^{'}}_T^k)

[ mathbf{E}_k^{'}=USigma V^T ag{2-6} ]

取左奇异矩阵(U)的第1个列向量(mathbf{u}_1=U(cdot,1))作为(mathbf{d}_k)，即(mathbf{d}_k=mathbf{u}_1)，取右奇异矩阵的第1个行向量与第1个奇异值的乘积作为(mathbf{x{'}}_T^k)，即(mathbf{x{'}}^k_T=Sigma(1,1)V^T(1,cdot))。得到(mathbf{x{'}}^k_T)后，将其对应地更新到原(mathbf{x}_T^k)。

注： 式（2-6）所求得的奇异值矩阵(Sigma)中的奇异值应从大到小排列；同样也有(mathbf{x{'}}^k_T=Sigma(1,1)V(cdot,1)^T)，这与上面(mathbf{x{'}}^k_T)的求法是相等的。

2.3 字典学习算法实现

据2.2小节，利用稀疏算法求解得到稀疏矩阵(mathbf{X})后，逐列更新字典，有如下算法1.1。

算法1.1：字典学习（K-SVD）

输入：原始样本，字典，稀疏矩阵
输出：字典，稀疏矩阵

1. 初始化： 从原始样本(Y in mathbf{R}^{m imes n})随机取(K)个列向量或者取它的左奇异矩阵的前(K)个列向量({mathbf{d}_1,mathbf{d}_2,cdots,mathbf{d}_K})作为初始字典的原子，得到字典(mathbf{D}^{(0)} in mathbf{R}^{m imes K})。令(j=0)，重复下面步骤2-3，直到达到指定的迭代步数，或收敛到指定的误差：

2. 稀疏编码： 利用字典上一步得到的字典(mathbf{D}^{(j)})，稀疏编码，得到(mathbf{X}^{(j)}inmathbf{R}^{K imes n})。

3. 字典更新： 逐列更新字典(mathbf{D}^{(j)})，字典的列(mathbf{d}_k in {mathbf{d}_1,mathbf{d}_2,cdots,mathbf{d}_K})

   • 当更新(mathbf{d}_k)时，计算误差矩阵(mathbf{E}_k)

     [ mathbf{E}_k=mathbf{Y}-sum_{j e k}mathbf{d}_jmathbf{x}^j_T. ]

   • 取出稀疏矩阵第(k)个行向量(mathbf{x}^k_T)不为0的索引的集合(omega_k = {i|1le ile n, mathbf{x}_T^k(i) e 0})，(mathbf{x'}_T^{k} = {mathbf{x}_T^k(i)|1le ile n, mathbf{x}_T^k(i) e 0})

   • 从(mathbf{E}_k)取出对应(omega_k)不为0的列，得到(mathbf{E}_k^{'}).

   • 对(mathbf{E}_k^{'})作奇异值分解(mathbf{E}_k=USigma V^T)，取(U)的第1列更新字典的第(k)列，即(mathbf{d}_k=U(cdot,1))；令(mathbf{x'}^k_T=Sigma(1,1)V(cdot,1)^T)，得到(mathbf{x{'}}^k_T)后，将其对应地更新到原(mathbf{x}_T^k)。

   • (j = j + 1)

3、字典学习Python实现

以下实验的运行环境为python3.6+jupyter5.4。

载入数据

import numpy as np
import pandas as pd
from scipy.io import loadmat

train_data_mat = loadmat("../data/train_data2.mat")

train_data = train_data_mat["Data"]
train_label = train_data_mat["Label"]

print(train_data.shape, train_label.shape)

注： 上面的数据集，可以随便使用一个，也可以随便找一个张图片。

初始化字典

u, s, v = np.linalg.svd(train_data)
n_comp = 50
dict_data = u[:, :n_comp]

字典更新

def dict_update(y, d, x, n_components):
    """
    使用KSVD更新字典的过程
    """
    for i in range(n_components):
        index = np.nonzero(x[i, :])[0]
        if len(index) == 0:
            continue
        # 更新第i列
        d[:, i] = 0
        # 计算误差矩阵
        r = (y - np.dot(d, x))[:, index]
        # 利用svd的方法，来求解更新字典和稀疏系数矩阵
        u, s, v = np.linalg.svd(r, full_matrices=False)
        # 使用左奇异矩阵的第0列更新字典
        d[:, i] = u[:, 0]
        # 使用第0个奇异值和右奇异矩阵的第0行的乘积更新稀疏系数矩阵
        for j,k in enumerate(index):
            x[i, k] = s[0] * v[0, j]
    return d, x

注： 上面代码的16~17需要注意python的numpy中的普通索引和花式索引的区别，花式索引会产生一个原数组的副本，所以对花式索引的操作并不会改变原数据，因此不能像第10行一样，需利用直接索引更新x。

迭代更新求解

可以指定迭代更新的次数，或者指定收敛的误差。

from sklearn import linear_model

max_iter = 10
dictionary = dict_data

y = train_data
tolerance = 1e-6

for i in range(max_iter):
    # 稀疏编码
    x = linear_model.orthogonal_mp(dictionary, y)
    e = np.linalg.norm(y - np.dot(dictionary, x))
    if e < tolerance:
        break
    dict_update(y, dictionary, x, n_comp)

sparsecode = linear_model.orthogonal_mp(dictionary, y)

train_restruct = dictionary.dot(sparsecode)