矩陣補全（Matrix Completion）和缺失值預處理

時間 2019-12-05

原文原文鏈接

目錄python

矩陣補全（Matrix Completion），就是補上一個含缺失值矩陣的缺失部分。算法

矩陣補全能夠經過矩陣分解（matrix factorization）將一個含缺失值的矩陣 X 分解爲兩個（或多個）矩陣，而後這些分解後的矩陣相乘就能夠獲得原矩陣的近似 X'，咱們用這個近似矩陣 X' 的值來填補原矩陣 X 的缺失部分。app

矩陣補全有不少方面的應用，如推薦系統、缺失值預處理。dom

除了 EM 算法、樹模型，機器學習中的大多數算法都須要輸入的數據是不含缺失值的。在 deep learning 模型中，經過梯度的計算公式就能夠發現，若是 feature 中含有缺失值，那麼梯度也會含缺失值，梯度也就未知了。對缺失值的處理是在模型訓練開始前就應該完成的，故也稱爲預處理。機器學習

數據缺失在實際場景中不可避免，對於一個包含 \(n\) 個 samples，每一個 sample 有 \(m\) 個 features 的數據集 \(D\)，咱們能夠將該數據集 \(D\) 整理爲一個 \(n×m\) 的矩陣 \(X\)。ide

經過矩陣分解補全矩陣是一種處理缺失值的方式，但在介紹以前，先介紹一些簡單經常使用的缺失值預處理方式。函數

1 經常使用的缺失值預處理方式

1.1 不處理

不進行缺失值預處理，缺了就缺了，找一個對缺失值不敏感的算法（如「樹模型」），直接訓練。post

1.2 剔除

對於矩陣 \(X\) 中缺失值不少的行或列，直接剔除。學習

缺失值較多的行，即一個 sample 的不少 features 都缺失了；缺失值較多的列，即大部分 samples 都沒有該 feature。剔除這些 samples 或 features，而不是填充它們，避免引入過多的噪聲。spa

當數據超級多時，咱們甚至能夠對含有缺失值的樣本直接剔除，當剔除的比例不大時，這也徹底能夠接受。

1.3 填充

1.3.1 簡單填充

在矩陣 \(X\) 的每一個缺失位置上填上一個數，來代替缺失值。填一個數也不能亂來，若是 feature 表明年齡，那麼確定要填正數；若是 feature 表明性別，那麼填 0 或 1 更合適（0 表明男，1 表明女）。

通常有如下幾種簡單的填充值：（均值和衆數都是在一個 feature 下計算，即在矩陣 \(X\) 的每一列中計算均值和衆數）

填 0
填均值
填衆數

1.3.2 建模填充

這種方式經過觀察缺失的 feature 和其它已有的 features 之間的聯繫，創建一個統計模型或者回歸模型，而後而後預測缺失 feature 的值應該是什麼。

用 EM 算法估計缺失值也能夠歸爲這一類。

固然，經常使用的缺失值處理方式還有許多，這裏就再也不列舉了。能夠看看博客 SAM'S NOTE。

2 利用矩陣分解補全缺失值

若是矩陣 \(X\) 不含缺失值，那麼矩陣分解能夠將矩陣 \(X\) 分解成兩個矩陣 \(U\) (大小 \(m×k\))、\(V\) (大小 \(m×k\))，其中 \(k < \min\{m, n\}\)，則：
\[ X = UV^{\top} \]

由於 \(k < \min\{m, n\}\)，因此 \(rank(U) \le k\)、\(rank(V) \le k\)，該矩陣分解又叫作低秩矩陣分解（low-rank matrix factorization）。

那麼爲何 \(k < \min\{m, n\}\)？

在 samples 和 features 之間存在 k 個關係，每一個關係的具體含義不得而知，但若是 \(k \ge \min\{m, n\}\)，那麼意味着每一個 sample 和 feature 之間能夠構建一個的關係，而其它的 samples 或者 features 能夠和該關係基本無關，體如今矩陣 \(U\)（或 \(V\)）中就是某一列僅有一個元素不爲0，這是不合理的。（參考矩陣分解用在推薦系統方面的解釋）
當 k 越大，計算量也會越大。

若是矩陣 \(X\) 是完整的，那麼矩陣分解 \(X = UV^{\top}\) 徹底沒問題，但如今 \(X\) 中含了缺失值，故沒有辦法用線性代數的知識直接進行矩陣分解，咱們須要一種近似的解法——梯度降低法。

這個時候咱們令 \(X \approx \hat X = UV^{\top}\)，\(\|X - \hat X\|_{\mathrm{F}}^2\) 表示含缺失值的原矩陣 \(X\) 和還原後的近似矩陣 \(\hat X\) 之間偏差的平方（Square error），或者稱之爲 reconstruction error，固然 \(\|X - \hat X\|_{\mathrm{F}}^2\) 的計算只能在不含缺失值的項上。（\(\|\cdot\|_{\mathrm{F}}\) 表示 Frobenius norm。）

文獻中通常會將 reconstruction error \(\|X - \hat X\|_{\mathrm{F}}^2\) 記爲 \(\left\|\mathcal{R}_{\Omega}(X - \hat{X})\right\|_{\mathrm{F}}^{2}\)，其中 \(\left[\mathcal{R}_{\Omega}(X - \hat X)\right]_{ij}=\left\{\begin{array}{ll}{x_{ij} - \hat{x}_{ij}} & {\text { if }(i, j) \in \Omega} \\ {0} & {\text { otherwise }}\end{array}\right.\)，\(\Omega\) 表示非缺失值矩陣元素下標的集合。這裏爲了簡便，直接使用 \(\|X - \hat X\|_{\mathrm{F}}^2\)，知道只在不含缺失值的項上計算平方和便可。

咱們的目標的是找到矩陣 \(X\) 的近似矩陣 \(\hat X\)，經過 \(\hat X\) 中對應的值來填充 \(X\) 中缺失的部分。而想要找到 \(\hat X\)，就是要找到矩陣 \(U\) 和 \(V\)。固然 \(\hat X\) 要儘量像 \(X\)，體如今函數上就是 \(\min \|X - \hat X\|_{\mathrm{F}}^2\)。

NOTE：如下矩陣的範數都默認爲 Frobenius norm。

Loss function \(J\) 爲:
\[ \begin{split} J &= \|X - \hat X\|^2 \\ &= \|X - UV^{\top}\|^2 \\ &= \sum_{i, j, x_{ij} \not = nan} (x_{ij} - \sum_{l = 1}^k u_{il}v_{jl})^2 \end{split} \]

其中，\(i,j\) 分別表示矩陣 \(X\) 的行和列，要求 \(x_{ij} \not = nan\)，不然沒有辦法求最小值了。上式中，未知的就是 \(u_{il}, v_{jl}\)，也是咱們想要求的。

隨機初始化矩陣 \(U, V\)，loss function \(J\) 就能夠獲得一個偏差，基於該偏差計算梯度，而想要更新 \(U, V\)，只須要按照梯度降低的公式來便可。
令:
\[ e_{ij} = x_{ij} - \sum_{l = 1}^k u_{il}v_{jl} \]

則梯度爲：
\[ \begin{split} &\frac{\partial J}{\partial u_{il}} = - 2e_{ij}v_{jl} \\ &\frac{\partial J}{\partial v_{jl}} = - 2e_{ij}u_{il} \end{split} \]

梯度降低更新公式：
\[ \begin{split} &u_{il} = u_{il} - \alpha\frac{\partial J}{\partial u_{il}} = u_{il} + 2\alpha e_{ij}v_{jl} \\ &v_{jl} = v_{jl} - \alpha\frac{\partial J}{\partial v_{jl}} = u_{il} + 2\alpha e_{ij}u_{il} \end{split} \]

算法到這裏其實就能夠用了，但爲了更加完美，能夠考慮如下步驟，加入正則項和偏置。

加入正則項

加入正則項，保證矩陣 \(U,V\) 中元素不要太大，此時 loss function \(J\) 以下所示：
\[ \begin{split} J &= \|X - \hat X\|^2 + \frac{\beta}{2}(\|U\|^2 + \|V\|^2) \\ &=\sum_{i, j, x_{ij} \not = nan} (x_{ij} - \sum_{l = 1}^k u_{il}v_{jl})^2 + \frac{\beta}{2}(\sum_{i,l} u_{il}^2 + \sum_{j, l} v_{jl}^2) \end{split} \]

則梯度爲：
\[ \begin{split} &\frac{\partial J}{\partial u_{il}} = - 2e_{ij}v_{jl} + \beta u_{il} \\ &\frac{\partial J}{\partial v_{jl}} = - 2e_{ij}u_{il} + \beta v_{jl} \end{split} \]

此時梯度降低更新公式爲：
\[ \begin{split} &u_{il} = u_{il} - \alpha\frac{\partial J}{\partial u_{il}} = u_{il} + \alpha(2e_{ij}v_{jl} - \beta u_{il}) \\ &v_{jl} = v_{jl} - \alpha\frac{\partial J}{\partial v_{jl}} = u_{il} + \alpha(2e_{ij}u_{il} - \beta v_{jl}) \end{split} \]

加入偏置

偏置能夠理解爲每一個樣本都有其特性，每一個feature也有其特色，故能夠加入 bias 來控制。bias 分爲三種，第一種是矩陣 \(X\) 總體的的 bias，記爲 \(b\)，那麼 \(b = mean(X)\)，便可以用矩陣 \(X\) 中存在元素的均值來賦值；第二種是 sample 的 bias，記爲 \(b\_u_{i}\)；第三種是 feature 的 bias，記爲 \(b\_v_j\)。
則：
\[ \hat x_{ij} = \sum_{l = 1}^{k} u_{il}v_{jl} + (b + b\_u_i + b\_v_j) \]

其中，\(b = \frac{\sum_{i, j, x_{ij} \not = nan} x_{ij}}{N}\)，\(N\) 表示分子求和元素的個數。
則 loss function \(J\) 爲：
\[ \begin{split} J &= \|X - \hat X\|^2 + \frac{\beta}{2}(\|U\|^2 + \|V\|^2 + b\_u^2 + b\_v^2) \\ &=\sum_{i, j, x_{ij} \not = nan} (x_{ij} - \sum_{l = 1}^k u_{il}v_{jl} - b - b\_u_i - b\_v_j)^2 \\ &+ \frac{\beta}{2}(\sum_{i,l} u_{il}^2 + \sum_{j, l} v_{jl}^2 + \sum_{i} b\_u_i^2 +\sum_{j} b\_v_j^2) \end{split} \]

再加入 bias 後，令
\[ e_{ij} = x_{ij} - \sum_{l = 1}^k u_{il}v_{jl} - b - b\_u_i - b\_v_j \]

則梯度爲：
\[ \begin{split} \frac{\partial J}{\partial u_{il}} &= - 2e_{ij}v_{jl} + \beta u_{il} \\ \frac{\partial J}{\partial v_{jl}} &= - 2e_{ij}u_{il} + \beta v_{jl} \\ \frac{\partial J}{\partial b\_u_i} &= -2e_{ij} + \beta b\_u_i \\ \frac{\partial J}{\partial b\_v_j} &= -2e_{ij} + \beta b\_v_j \end{split} \]

此時梯度降低更新公式爲：
\[ \begin{split} u_{il} &= u_{il} + \alpha(2e_{ij}v_{jl} - \beta u_{il}) \\ v_{jl} &= u_{il} + \alpha(2e_{ij}u_{il} - \beta v_{jl}) \\ b\_u_i &= b\_u_i + \alpha(2e_{ij} - \beta b\_u_i) \\ b\_v_j &= b\_v_j + \alpha(2e_{ij} - \beta b\_v_j) \end{split} \]

3 矩陣分解補全缺失值代碼實現

import numpy as np

class MF():

    def __init__(self, X, k, alpha, beta, iterations):
        """
        Perform matrix factorization to predict np.nan entries in a matrix.
        Arguments
        - X (ndarray)   : sample-feature matrix
        - k (int)       : number of latent dimensions
        - alpha (float) : learning rate
        - beta (float)  : regularization parameter
        """

        self.X = X
        self.num_samples, self.num_features = X.shape
        self.k = k
        self.alpha = alpha
        self.beta = beta
        self.iterations = iterations
        # True if not nan
        self.not_nan_index = (np.isnan(self.X) == False)

    def train(self):
        # Initialize factorization matrix U and V
        self.U = np.random.normal(scale=1./self.k, size=(self.num_samples, self.k))
        self.V = np.random.normal(scale=1./self.k, size=(self.num_features, self.k))

        # Initialize the biases
        self.b_u = np.zeros(self.num_samples)
        self.b_v = np.zeros(self.num_features)
        self.b = np.mean(self.X[np.where(self.not_nan_index)])
        # Create a list of training samples
        self.samples = [
            (i, j, self.X[i, j])
            for i in range(self.num_samples)
            for j in range(self.num_features)
            if not np.isnan(self.X[i, j])
        ]

        # Perform stochastic gradient descent for number of iterations
        training_process = []
        for i in range(self.iterations):
            np.random.shuffle(self.samples)
            self.sgd()
            # total square error
            se = self.square_error()
            training_process.append((i, se))
            if (i+1) % 10 == 0:
                print("Iteration: %d ; error = %.4f" % (i+1, se))

        return training_process

    def square_error(self):
        """
        A function to compute the total square error
        """
        predicted = self.full_matrix()
        error = 0
        for i in range(self.num_samples):
            for j in range(self.num_features):
                if self.not_nan_index[i, j]:
                    error += pow(self.X[i, j] - predicted[i, j], 2)
        return error

    def sgd(self):
        """
        Perform stochastic graident descent
        """
        for i, j, x in self.samples:
            # Computer prediction and error
            prediction = self.get_x(i, j)
            e = (x - prediction)

            # Update biases
            self.b_u[i] += self.alpha * (2 * e - self.beta * self.b_u[i])
            self.b_v[j] += self.alpha * (2 * e - self.beta * self.b_v[j])

            # Update factorization matrix U and V
            """
            If RuntimeWarning: overflow encountered in multiply,
            then turn down the learning rate alpha.
            """
            self.U[i, :] += self.alpha * (2 * e * self.V[j, :] - self.beta * self.U[i,:])
            self.V[j, :] += self.alpha * (2 * e * self.U[i, :] - self.beta * self.V[j,:])

    def get_x(self, i, j):
        """
        Get the predicted x of sample i and feature j
        """
        prediction = self.b + self.b_u[i] + self.b_v[j] + self.U[i, :].dot(self.V[j, :].T)
        return prediction

    def full_matrix(self):
        """
        Computer the full matrix using the resultant biases, U and V
        """
        return self.b + self.b_u[:, np.newaxis] + self.b_v[np.newaxis, :] + self.U.dot(self.V.T)

    def replace_nan(self, X_hat):
        """
        Replace np.nan of X with the corresponding value of X_hat
        """
        X = np.copy(self.X)
        for i in range(self.num_samples):
            for j in range(self.num_features):
                if np.isnan(X[i, j]):
                    X[i, j] = X_hat[i, j]
        return X

if __name__ == '__main__':
    X = np.array([
        [5, 3, 0, 1],
        [4, 0, 0, 1],
        [1, 1, 0, 5],
        [1, 0, 0, 4],
        [0, 1, 5, 4],
    ], dtype=np.float)
    # replace 0 with np.nan
    X[X == 0] = np.nan
    print(X)
    # np.random.seed(1)
    mf = MF(X, k=2, alpha=0.1, beta=0.1, iterations=100)
    mf.train()
    X_hat = mf.full_matrix()
    X_comp = mf.replace_nan(X_hat)

    print(X_hat)
    print(X_comp)
    print(X)