字典學習（Dictionary Learning, KSVD）詳解

注：字典學習也是一種數據降維的方法，這里我用到SVD的知識，對SVD不太理解的地方，可以看看這篇博客：《SVD（奇異值分解）小結 》；數據集：https://pan.baidu.com/s/1ZmpUSIscy4VltcimwwIWew

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 \times K}\)，稀疏線性表示原始樣本\(\mathbf{Y} \in \mathbf{R}^{m \times n}\)（其中\(m\)表示樣本數，\(n\)表示樣本的屬性），即有\(\mathbf{Y=DX}\)（這只是我們理想的情況），其中\(\mathbf{X} \in \mathbf{R}^{K \times n}\)為稀疏矩陣，可以將上述問題用數學語言描述為如下優化問題

\[\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}\)可以利用已有經典算法求解，如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），我們有

\[\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所示，這樣描述更加清晰。

圖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}_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\)中的奇異值應從大到小排列；同樣也有\(\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 \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\)

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)