---

typora-root-url: md_pic

---

字典学习（Dictionary Learning）和稀疏表示（Sparse Representation）在学术界的正式称谓应该是稀疏字典学习（Sparse Dictionary Learning）。该算法理论包含两个阶段：字典构建阶段（Dictionary Generate）和利用字典（稀疏的）表示样本阶段（Sparse coding with a precomputed dictionary）。这两个阶段（如下图）的每个阶段都有许多不同算法可供选择，每种算法的诞生时间都不一样，以至于稀疏字典学习的理论提出者已变得不可考。笔者尝试找了Wikipedia和Google Scolar都无法找到这一系列理论的最早发起人。

​ 使用字典学习的目的：（1）是对于庞大数据集进行降维表示；（2）字典学习总是尝试学习蕴藏在样本背后最质朴的特征（假如样本最质朴的特征就是样本最好的特征），这两条原因同时也是这两年深度学习之风日盛的情况下字典学习也开始随之升温的原因。稀疏表示的本质：用尽可能少的资源表示尽可能多的知识，这种表示还能带来一个附加的好处，即计算速度快。

稀疏字典学习的Python实现

首先是各种工具包的导入和测试样例的导入

from time import time
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
from sklearn.decomposition import MiniBatchDictionaryLearning
from sklearn.feature_extraction.image import extract_patches_2d
from sklearn.feature_extraction.image import reconstruct_from_patches_2d
from sklearn.utils.testing import SkipTest
from sklearn.utils.fixes import sp_version
if sp_version < (0, 12):
    raise SkipTest("Skipping because SciPy version earlier than 0.12.0 and "
                   "thus does not include the scipy.misc.face() image.")
try:
    from scipy import misc
    face = misc.face(gray=True)
except AttributeError:
    # Old versions of scipy have face in the top level package
    face = sp.face(gray=True)

通过测试样例计算字典V

# Convert from uint8 representation with values between 0 and 255 to
# a floating point representation with values between 0 and 1.
face = face / 255.0
# downsample for higher speed
face = face[::2, ::2] + face[1::2, ::2] + face[::2, 1::2] + face[1::2, 1::2]
face = face / 4.0
height, width = face.shape
# Distort the right half of the image
print('Distorting image...')
distorted = face.copy()
distorted[:, width // 2:] += 0.075 * np.random.randn(height, width // 2)
# Extract all reference patches from the left half of the image
print('Extracting reference patches...')
t0 = time()
patch_size = (7, 7)
data = extract_patches_2d(distorted[:, :width // 2], patch_size)
data = data.reshape(data.shape[0], -1)
data -= np.mean(data, axis=0)
data /= np.std(data, axis=0)
print('done in %.2fs.' % (time() - t0))
print('Learning the dictionary...')
t0 = time()
dico = MiniBatchDictionaryLearning(n_components=100, alpha=1, n_iter=500)
V = dico.fit(data).components_
dt = time() - t0
print('done in %.2fs.' % dt)
plt.figure(figsize=(4.2, 4))
for i, comp in enumerate(V[:100]):
    plt.subplot(10, 10, i + 1)
    plt.imshow(comp.reshape(patch_size), cmap=plt.cm.gray_r,
               interpolation='nearest')
    plt.xticks(())
    plt.yticks(())
plt.suptitle('Dictionary learned from face patches\n' +
             'Train time %.1fs on %d patches' % (dt, len(data)),
             fontsize=16)
plt.subplots_adjust(0.08, 0.02, 0.92, 0.85, 0.08, 0.23)#left, right, bottom, top, wspace, hspace

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

字典学习也是一种数据降维的方法。

字典学习思想

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

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

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

字典学习数据模型

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

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

​ 字典学习的主要思想是，利用包含$K$K个原子$d_k$dk​的字典矩阵$\mathbf{D}\in\mathbf{R}^{m\times K}$D∈Rm×K ，稀疏线性表示原始样本$\mathbf{Y} \in \mathbf{R}^{m \times n}$Y∈Rm×n（其中$m$m表示样本数，$n$n表示样本的属性），即有$\mathbf{Y}=\mathbf{D} \mathbf{S}$Y=DS（这只是我们理想的情况），其中$\mathbf{S} \in \mathbf{R}^{K \times n}$S∈RK×n为稀疏矩阵，可以将上述问题用数学语言描述为如下优化问题
$\min _{\mathbf{D}, \mathbf{X}}\|\mathbf{Y}-\mathbf{D} \mathbf{S}\|_{F}^{2}, \quad \text { s.t. } \forall i,\left\|\mathbf{s}_{i}\right\|_{0} \leq T_{0}$D,Xmin​∥Y−DS∥F2​, s.t. ∀i,∥si​∥0​≤T0​
上式中，$\mathbf{S}$S为稀疏编码的矩阵，$\mathbf{s}_{i}(i=1,2, \cdots, K)$si​(i=1,2,⋯,K)为该矩阵中的行向量，代表字典矩阵的系数。$\left\|\mathbf{s}_{i}\right\|_{0}$∥si​∥0​表示零阶范数，它表示向量中不为0的数的个数。

求解问题

式（2-1）的目标函数表示，我们要最小化查完的字典与原始样本的误差，即要尽可能还原出原始样本；它的限的制条件$\left\|\mathbf{s}_{i}\right\|_{0} \leq T_{0}$∥si​∥0​≤T0​，表示查字典的方式要尽可能简单，即$\mathbf{X}$X要尽可能稀疏。

式（2-1）或式（2-2）是一个带有约束的优化问题，可以利用拉格朗日乘子法将其转化为无约束优化问题
$\min _{\mathbf{D}, \mathbf{S}}\|\mathbf{Y}-\mathbf{D} \mathbf{S}\|_{F}^{2}+\lambda\left\|\mathbf{s}_{i}\right\|_{1}$D,Smin​∥Y−DS∥F2​+λ∥si​∥1​
注： 我们将$\left\|\mathbf{s}_{i}\right\|_{0}$∥si​∥0​用$\left\|\mathbf{s}_{i}\right\|_{1}$∥si​∥1​代替，主要是更加便于$\left\|\mathbf{s}_{i}\right\|_{1}$∥si​∥1​求解。

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

Hello World

假设$\mathbf{X}$X是已知的，我们逐列更新字典。下面我们仅更新字典的第$k$k列，记$d_k$dk​为字典$\mathbf{D}$D的第$k$k列向量，记$\mathbf{S}_{T}^{k}$STk​为稀疏矩阵$\mathbf{S}$S的第$k$k行向量，那么对式（2-1），我们有
$\begin{aligned}\|\mathbf{Y}-\mathbf{D} \mathbf{S}\|_{F}^{2} &=\left\|\mathbf{Y}-\sum_{j=1}^{K} \mathbf{d}_{j} \mathbf{s}_{T}^{j}\right\|_{F}^{2} \\ &=\left\|\left(\mathbf{Y}-\sum_{j \neq k} \mathbf{d}_{j} \mathbf{s}_{T}^{j}\right)-\mathbf{d}_{k} \mathbf{s}_{T}^{k}\right\|_{F}^{2} \\ &=\left\|\mathbf{E}_{k}-\mathbf{d}_{k} \mathbf{s}_{T}^{k}\right\|_{F}^{2} \end{aligned}$∥Y−DS∥F2​​=∥∥∥∥∥​Y−j=1∑K​dj​sTj​∥∥∥∥∥​F2​=∥∥∥∥∥∥​⎝⎛​Y−j̸​=k∑​dj​sTj​⎠⎞​−dk​sTk​∥∥∥∥∥∥​F2​=∥∥​Ek​−dk​sTk​∥∥​F2​​
上式中，残差为$\mathbf{E}_{k}=\mathbf{Y}-\sum_{j \neq k} \mathbf{d}_{j} \mathbf{s}_{T}^{j}$Ek​=Y−∑j̸​=k​dj​sTj​。此时优化问题可描述为
$\min _{\mathbf{d}_{k}, \mathbf{s}_{T}^{k}}\left\|\mathbf{E}_{k}-\mathbf{d}_{k} \mathbf{s}_{T}^{k}\right\|_{F}^{2}$dk​,sTk​min​∥∥​Ek​−dk​sTk​∥∥​F2​
​ 因此，我们需要求出最优的$\mathbf{d}_{k}$dk​, $\mathbf{s}_{T}^{k}$sTk​，这是一个最小二乘问题，可以利用最小二乘的方法求解，或者可以利用SVD进行求解，这里利用SVD求解出两个优化变量。

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

[外链图片转存失败(img-S8uc5Drn-1566026922960)(/sparseDictionaryLearning(extractionEx)].png)

如上图，假设我们要更新第0列原子，我们将$\mathbf{s}_{T}^{k}$sTk​中为零的位置找出来，然后把对$\mathbf{E}_{k}$Ek​应的位置删除，得到$\mathbf{E}_{k}^{\prime}$Ek′​，此时优化问题可描述为
$\min _{\mathrm{d}_{k}, \mathrm{s}_{F}^{k}}\left\|\mathbf{E}_{k}^{\prime}-\mathbf{d}_{k} \mathbf{s}_{T}^{\prime k}\right\|_{F}^{2}$dk​,sFk​min​∥∥​Ek′​−dk​sT′k​∥∥​F2​
因此我们需要求出最优的$\mathbf{d}_{k}$dk​, $\mathbf{s}_{T}^{\prime k}$sT′k​
$\mathbf{E}_{k}^{\prime}=U \Sigma V^{T}$Ek′​=UΣVT
取左奇异矩阵$U$U的第1个列向量$\mathbf{u}_{1}=U(\cdot, 1)$u1​=U(⋅,1)作为$\mathbf{d}_{k}$dk​，即$\mathbf{d}_{k}=\mathbf{u}_{1}$dk​=u1​，取右奇异矩阵的第1个行向量与第1个奇异值的乘积作为$\mathbf{x}_{T}^{\prime k}$xT′k​，即$\mathbf{s}_{T}^{\prime k}=\Sigma(1,1) V^{T}(1, \cdot)$sT′k​=Σ(1,1)VT(1,⋅)。得到$\mathbf{S}_{T}^{\prime k}$ST′k​后，将其对应地更新到原$\mathbf{S}_{T}^{k}$STk​。

字典学习算法实现

好的呀

利用稀疏算法求解得到稀疏矩阵$\mathbf{X}$X后，逐列更新字典，有如下算法1.1。

算法1.1 字典学习（K-SVD）
输入：原始样本$Y$Y，字典，稀疏矩阵
输出：字典，稀疏矩阵
1. 初始化：从原始样本$Y \in \mathbf{R}^{m \times n}$Y∈Rm×n随机取$K$K个列向量或者取它的左奇异矩阵的前$K$K个列向量$\left\{\mathbf{d}_{1}, \mathbf{d}_{2}, \cdots, \mathbf{d}_{K}\right\}${d1​,d2​,⋯,dK​}作为初始字典的原子，得到字典$\mathbf{D}^{(0)} \in \mathbf{R}^{m \times K}$D(0)∈Rm×K。令$j=0$j=0，重复下面步骤2-3，直到达到指定的迭代步数，或收敛到指定的误差。
2. 稀疏编码：利用字典上一步得到的字典$\mathbf{D}^{(j)}$D(j)，稀疏编码，得到$\mathbf{S}^{(j)} \in \mathbf{R}^{K \times n}$S(j)∈RK×n。
3. 字典更新：逐列更新字典$\mathbf{D}^{(j)}$D(j)，以及字典的列$\mathbf{d}_{k} \in\left\{\mathbf{d}_{1}, \mathbf{d}_{2}, \cdots, \mathbf{d}_{K}\right\}$dk​∈{d1​,d2​,⋯,dK​}：
(1) 当更新$\mathrm{d}_{k}$dk​时，计算误差矩阵$\mathbf{E}_{k}$Ek​
$\mathbf{E}_{k}=\mathbf{Y}-\sum_{j \neq k} \mathbf{d}_{j} \mathbf{s}_{T}^{j}$Ek​=Y−j̸​=k∑​dj​sTj​
(2) 取出稀疏矩阵第$k$k个行向量$\mathbf{s}_{T}^{k}$sTk​不为0的索引的集合$\omega_{k}=\left{i
(3) 从$\mathbf{E}_k$Ek​取出对应$\omega_k$ωk​不为0的列，得到$\mathbf{E}_k^{\prime}$Ek′​。
(4) 对$\mathbf{E}_k^{\prime}$Ek′​作奇异值分解$\mathbf{E}_{k}=U \Sigma V^{T}$Ek​=UΣVT，取的第1列更新字典的第$k$k列，即$\mathbf{d}_{k}=U(\cdot, 1)$dk​=U(⋅,1)；令$\mathbf{s}_{T}^{\prime k}=\Sigma(1,1) V(\cdot, 1)^{T}$sT′k​=Σ(1,1)V(⋅,1)T，得到$\mathbf{s}_{T}^{\prime k}$sT′k​后，将其对应地更新到原$\mathbf{s}_{T}^{k}$sTk​。
$j=j+1$j=j+1

字典学习Python实现

import numpy as np
import matplotlib.pyplot as plt
import scipy.io
import sklearn.linear_model

train_data_mat = scipy.io.loadmat("Indian_pines_corrected.mat")
train_data = train_data_mat['indian_pines_corrected']
train_label_dat = scipy.io.loadmat("Indian_pines_gt.mat")
train_label = train_label_dat['indian_pines_gt']
train_data = np.swapaxes(train_data, 0, 2)
train_data = np.swapaxes(train_data, 1, 2)
train_data = train_data[0, :, :]
print(train_data.shape, train_label.shape)

# Initialize the dictionary
n_comp = 50
u_tmp, s_tmp, v_tmp = np.linalg.svd(train_data)
dict_data = u_tmp[:, :n_comp]


def dict_update(y, d, x, n_components):
    """
    Update the dictionary using K-SVD
    """
    for jj in np.arange(n_components):
        index = np.nonzero(x[jj, :])[0]
        if len(index) == 0:
            continue
        # update the jj-th column
        d[:, jj] = 0
        # calculate the error matrix
        r = (y - np.dot(d, x))[:, index]
        # Update the dictionary and the sparse coefficiency based on svd
        u_r, s_r, v_r = np.linalg.svd(r, full_matrices=False)
        # 使用左奇异矩阵的第0列更新字典
        d[:, jj] = u_r[:, 0]
        # 使用第0个奇异值和右奇异矩阵的第0行的乘积更新稀疏系数矩阵
        for j, k in enumerate(index):
            x[jj, k] = s_r[0] * v_r[0, j]
    return d, x


max_iter = 10
dictionary = dict_data

e_tolerance = 1e-6

for i in range(max_iter):
    # sparse coding
    sc = sklearn.linear_model.orthogonal_mp(dictionary, train_data)
    error = np.linalg.norm(train_data - np.dot(dictionary, sc))
    if error < e_tolerance:
        break
    dictionary, sc = dict_update(train_data, dictionary, sc, n_comp)

sparsecode = sklearn.linear_model.orthogonal_mp(dictionary, train_data)
train_restruct = dictionary.dot(sparsecode)

plt.figure(figsize=(5, 3.3))
plt.subplot(1, 2, 1)
plt.title('Original Image')
plt.imshow(train_data.astype(np.uint16), cmap=plt.cm.gray,
           interpolation='nearest')
plt.subplot(1, 2, 2)
plt.title('Image')
plt.imshow(train_restruct.astype(np.uint16), cmap=plt.cm.gray,
           interpolation='nearest')
plt.show()