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

**注：**字典学习也是一种数据降维的方法，这里我用到SVD的知识，对SVD不太理解的地方，能够看看这篇博客：《SVD（奇异值分解）小结 》；数据集：https://pan.baidu.com/s/1ZmpUSIscy4VltcimwwIWew

一、字典学习思想

字典学习的思想应该源来实际生活中的<span style="color:#007acc;font-weight:normal">字典</span>的概念。字典是前辈们学习总结的精华，当咱们须要学习新的知识的时候，没必要与先辈们同样去学习先辈们全部学习过的知识，咱们能够参考先辈们给咱们总结的字典，经过查阅这些字典，咱们能够大体学会到这些知识。

为了将上述过程用准确的数学语言描述出来，咱们须要将<span style="color:#007acc;font-weight:normal">“总结字典”、“查阅字典”</span>作出一个更为准确的描述。就从咱们的常识出发：

1. 咱们一般会要求的咱们的字典尽量全面，也就是说总结出的字典不能漏下关键的知识点。
2. 查字典的时候，咱们想要咱们查字典的过程尽量简洁，迅速，准确。即，查字典要快、准、狠。
3. 查到的结果，要尽量地还原出原来知识。固然，若是要彻底还原出来，那么这个字典和查字典的方法会变得很是复杂，因此咱们只须要尽量地还原出原知识点便可。

注： 以上内容，彻底是本身的理解，若有不当之处，欢迎各位拍砖。

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

二、字典学习数学模型

2.1 数学描述

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

• “之前的知识”，更专业一点，咱们称之为<span style="color:#319595;font-weight:bold">原始样本</span>，用矩阵$\mathbf{Y}$表示；
• “字典”，咱们称之为<span style="color:#319595;font-weight:bold">字典矩阵</span>，用$\mathbf{D}$表示，“字典”中的词条，咱们称之为<span style="color:#319595;font-weight:bold">原子（atom）</span>，用列向量$\mathbf{d}_k$表示；
• “查字典的方法”，咱们称为<span style="color:#319595;font-weight:bold">稀疏矩阵</span>，用$\mathbf{X}$；
• “查字典的过程”，咱们能够用矩阵的乘法来表示，即$\mathbf{DX}$。

用数学语言描述，字典学习的主要思想是，利用包含$K$个原子$\mathbf{d}_k$的字典矩阵$\mathbf{D}\in \mathbf{R}^{m \times K}$，稀疏线性表示原始样本$\mathbf{Y} \in \mathbf{R}^{m \times n}$（其中$m$表示样本数，$n$表示样本的属性），即有$\mathbf{Y=DX}$（这只是咱们理想的状况），其中$\mathbf{X} \in \mathbf{R}^{K \times n}$为<span style="color:red;font-weight:bold">稀疏矩阵</span>，能够将上述问题用数学语言描述为以下优化问题

$$ \min_{\mathbf{D,\ X}}{|\mathbf{Y}-\mathbf{DX}|^2_F},\quad \text{s.t.}\ \forall i,\ |\mathbf{x}_i|_0 \le T_0 \tag{2-1} $$

或者

$$ \min_{\mathbf{D,\ X}}\sum_i|\mathbf{x}_i|0, \quad \text{s.t.}\ \min{\mathbf{D,\ X}}{|\mathbf{Y}-\mathbf{DX}|^2_F} \le \epsilon, \tag{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 \tag{2-3} $$

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

这里有两个优化变量$\mathbf{D,\ X}$，为解决这个优化问题，通常是固定其中一个优化变量，优化另外一个变量，如此交替进行。式（2-3）中的稀疏矩阵$\mathbf{X}$能够利用已有经典算法求解，如<span style="color:red;font-weight:normal">Lasso（Least Absolute Shrinkage and Selection Operator）、OMP（Orthogonal Matching Pursuit）</span>，这里我重点讲述如何更新字典$\mathbf{D}$，对更新$\mathbf{X}$很少作讨论。

假设$\mathbf{X}$是已知的，咱们<span style="color:red;font-weight:normal">逐列</span>更新字典。下面咱们仅更新字典的第$k$列，记$\mathbf{d}_k$为字典$\mathbf{D}$的第$k$列向量，记$\mathbf{x}^k_T$为稀疏矩阵$\mathbf{X}$的第$k$行向量，那么对式（2-1），咱们有

$$ \begin{aligned} {|\mathbf{Y}-\mathbf{DX}|^2_F} =&\left|\mathbf{Y}-\sum^K_{j=1}\mathbf{d}j\mathbf{x}^j_T\right|^2_F \ =&\left|\left(\mathbf{Y}-\sum{j\ne k}\mathbf{d}_j\mathbf{x}^j_T\right)-\mathbf{d}_k\mathbf{x}^k_T\right|^2_F\ =&\left|\mathbf{E}_k - \mathbf{d}_k\mathbf{x}_T^k \right|^2_F \end{aligned} \tag{2-4} $$

上式中残差$\mathbf{E}k=\mathbf{Y}-\sum{j\ne k}\mathbf{d}_j\mathbf{x}^j_T$，

此时优化问题可描述为

$$ \min_{\mathbf{d}_k,\ \mathbf{x}^k_T}\left|\mathbf{E}_k - \mathbf{d}_k\mathbf{x}_T^k \right|^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所示，这样描述更加清晰。

<div style="text-align:center;line-height:120%;padding:-1px 0;font-size:90%"> <img alt="fig1-1" src="http://images.cnblogs.com/cnblogs_com/endlesscoding/1356090/o_dict_note_1.png" style="width:150%;"/><br /> <div style="margin:0px 0px 0px 0">图2-1 提取部分残差</div> </div>

如上图，假设咱们要更新第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}_k\mathbf{x{'}}_T^{k} \right|^2_F \tag{2-5} $$

所以咱们须要求出最优的$\mathbf{d}_k,\ \mathbf{x^{'}}_T^k$

$$ \mathbf{E}_k^{'}=U\Sigma V^T \tag{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$中的奇异值应<span style="color:red;font-weight:normal">从大到小排列</span>；一样也有$\mathbf{x{'}}^k_T=\Sigma(1,1)V(\cdot,1)^T$，这与上面$\mathbf{x{'}}^k_T$的求法是相等的。

2.3 字典学习算法实现

据<span style="color:#007acc;font-weight:normal">2.2小节</span>，利用稀疏算法求解获得稀疏矩阵$\mathbf{X}$后，逐列更新字典，有以下算法1.1。

<div style="border-top:2px solid;"> <p style="color:#007acc;text-indent:0.3em;font-weight:bold;border-bottom:1px solid;margin:5px 0;font-size:1.1em">算法1.1：字典学习（K-SVD）</p> <span style="font-weight:bold;padding:0 5px">输入：</span>原始样本，字典，稀疏矩阵</br> <span style="font-weight:bold;padding:0 5px">输出：</span>字典，稀疏矩阵</br> </div>

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

2. 稀疏编码： 利用字典上一步获得的字典$\mathbf{D}^{(j)}$，稀疏编码，获得$\mathbf{X}^{(j)}\in\mathbf{R}^{K \times 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\ne k}\mathbf{d}_j\mathbf{x}^j_T. $$

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

   • 从$\mathbf{E}_k$取出对应$\omega_k$不为0的列，获得$\mathbf{E}_k^{'}$.

   • 对$\mathbf{E}_k^{'}$做奇异值分解$\mathbf{E}_k=U\Sigma 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$

<div style="border-top:2px solid;margin:-5px 0"></div>

三、字典学习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)

原文出处：https://www.cnblogs.com/endlesscoding/p/10090866.html