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

時間 2020-08-02

標籤字典學習 dictionary learning ksvd 詳解简体版

原文原文鏈接

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

一、字典學習思想

字典學習的思想應該源來實際生活中的<span style="color:#007acc;font-weight:normal">字典</span>的概念。字典是前輩們學習總結的精華，當咱們須要學習新的知識的時候，沒必要與先輩們同樣去學習先輩們全部學習過的知識，咱們能夠參考先輩們給咱們總結的字典，經過查閱這些字典，咱們能夠大體學會到這些知識。python

爲了將上述過程用準確的數學語言描述出來，咱們須要將<span style="color:#007acc;font-weight:normal">「總結字典」、「查閱字典」</span>作出一個更爲準確的描述。就從咱們的常識出發：算法

咱們一般會要求的咱們的字典儘量全面，也就是說總結出的字典不能漏下關鍵的知識點。
查字典的時候，咱們想要咱們查字典的過程儘量簡潔，迅速，準確。即，查字典要快、準、狠。
查到的結果，要儘量地還原出原來知識。固然，若是要徹底還原出來，那麼這個字典和查字典的方法會變得很是複雜，因此咱們只須要儘量地還原出原知識點便可。

注：以上內容，徹底是本身的理解，若有不當之處，歡迎各位拍磚。數組

下面，咱們要討論的就是如何將上述問題抽象成一個數學問題，並解決這個問題。less

二、字典學習數學模型

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} $$優化

或者ui

$$ \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。

初始化： 從原始樣本$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$

三、字典學習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

相關標籤/搜索

每日一句

每一个你不满意的现在，都有一个你没有努力的曾经。