注:字典學習也是一種數據降維的方法,這里我用到SVD的知識,對SVD不太理解的地方,可以看看這篇博客:《SVD(奇異值分解)小結 》;數據集:https://pan.baidu.com/s/1ZmpUSIscy4VltcimwwIWew
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\)表示樣本數,\(n\)表示樣本的屬性),即有\(\mathbf{Y=DX}\)(這只是我們理想的情況),其中\(\mathbf{X} \in \mathbf{R}^{K \times n}\)為稀疏矩陣,可以將上述問題用數學語言描述為如下優化問題
或者
上式中\(\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)是一個帶有約束的優化問題,可以利用拉格朗日乘子法將其轉化為無約束優化問題
注: 我們將\(\|\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),我們有
上式中殘差\(\mathbf{E}_k=\mathbf{Y}-\sum_{j\ne k}\mathbf{d}_j\mathbf{x}^j_T\),
此時優化問題可描述為
因此我們需要求出最優的\(\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所示,這樣描述更加清晰。

如上圖,假設我們要更新第0列原子,我們將\(\mathbf{x}_T^k\)中為零的位置找出來,然后把\(\mathbf{E}_k\)對應的位置刪除,得到\(\mathbf{E}_k^{'}\),此時優化問題可描述為
因此我們需要求出最優的\(\mathbf{d}_k,\ \mathbf{x^{'}}_T^k\)
取左奇異矩陣\(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)