字典學習（Dictionary Learning, KSVD）詳解
注：字典學習也是一種資料降維的方法，這裡我用到SVD的知識，對SVD不太理解的地方，可以看看這篇部落格：《SVD（奇異值分解）小結 》。
1、字典學習思想
字典學習的思想應該源來實際生活中的字典的概念。字典是前輩們學習總結的精華，當我們需要學習新的知識的時候，不必與先輩們一樣去學習先輩們所有學習過的知識，我們可以參考先輩們給我們總結的字典，通過查閱這些字典，我們可以大致學會到這些知識。
為了將上述過程用準確的數學語言描述出來，我們需要將“總結字典”、“查閱字典”做出一個更為準確的描述。就從我們的常識出發：
- 我們通常會要求的我們的字典儘可能全面，也就是說總結出的字典不能漏下關鍵的知識點。
- 查字典的時候，我們想要我們查字典的過程儘可能簡潔，迅速，準確。即，查字典要快、準、狠。
- 查到的結果，要儘可能地還原出原來知識。當然，如果要完全還原出來，那麼這個字典和查字典的方法會變得非常複雜，所以我們只需要儘可能地還原出原知識點即可。
注： 以上內容，完全是自己的理解，如有不當之處，歡迎各位拍磚。
下面，我們要討論的就是如何將上述問題抽象成一個數學問題，並解決這個問題。
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\)
\[ \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{eqnarray} {\|\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{eqnarray} \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）
輸入：原始樣本，字典，稀疏矩陣輸出：字典，稀疏矩陣
初始化： 從原始樣本\(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，直到達到指定的迭代步數，或收斂到指定的誤差：
稀疏編碼： 利用字典上一步得到的字典\(\mathbf{D}^{(j)}\)，稀疏編碼，得到\(\mathbf{X}^{(j)}\in\mathbf{R}^{K \times n}\)。
- 字典更新： 逐列更新字典\(\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)