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

通過矩陣分解補全矩陣用處很廣，如推薦系統、缺失值預處理等。本文先簡單介紹補全缺失值的一般方法，然後介紹如何用矩陣分解補全缺失值，最後給出含缺失值矩陣分解的python代碼。

目錄

• 1 常用的缺失值預處理方式
  • 1.1 不處理
  • 1.2 剔除
  • 1.3 填充
    • 1.3.1 簡單填充
    • 1.3.2 模組化填充
• 2 利用矩陣分解補全缺失值
• 3 矩陣分解補全缺失值代碼實作
• 4 通過矩陣分解補全矩陣的一些小問題
• References

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

矩陣補全可以通過矩陣分解（matrix factorization）将一個含缺失值的矩陣 X 分解為兩個（或多個）矩陣，然後這些分解後的矩陣相乘就可以得到原矩陣的近似 X'，我們用這個近似矩陣 X' 的值來填補原矩陣 X 的缺失部分。

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

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

資料缺失在實際場景中不可避免，對于一個包含 \(n\) 個 samples，每個 sample 有 \(m\) 個 features 的資料集 \(D\)，我們可以将該資料集 \(D\) 整理為一個 \(n×m\) 的矩陣 \(X\)。

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

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

對于矩陣 \(X\) 中缺失值很多的行或列，直接剔除。

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

當資料超級多時，我們甚至可以對含有缺失值的樣本直接剔除，當剔除的比例不大時，這也完全可以接受。

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

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

• 填 0
• 填 均值
• 填 衆數
• 填 中位數

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

用 EM 算法估計缺失值也可以歸為這一類。

當然，常用的缺失值處理方式還有許多，這裡就不再列舉了。可以看看部落格 SAM'S NOTE。

如果矩陣 \(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}

則梯度為：

&\frac{\partial J}{\partial u_{il}} = - 2e_{ij}v_{jl}

\\ &\frac{\partial J}{\partial v_{jl}} = - 2e_{ij}u_{il}

梯度下降更新公式：

&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}

算法到這裡其實就可以用了，但為了更加完美，可以考慮以下步驟，加入正則項和偏置。

加入正則項

加入正則項，保證矩陣 $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)

&\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}

此時梯度下降更新公式為：

&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})

加入偏置

偏置可以了解為每個樣本都有其特性，每個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\) 為：

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)

再加入 bias 後，令

\[e_{ij} = x_{ij} - \sum_{l = 1}^k u_{il}v_{jl} - b - b\_u_i - b\_v_j

\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

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)

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)
           

4.1 需不需要對 bias 進行正則化？

按照一般 deep learning 模型，是不對 bias 進行正則化的，而本文的代碼對 bias 進行了正則化，具體有沒有影響不得而知。

4.2 如果出現 "RuntimeWarning: overflow encountered in multiply" 等 Warning 造成最後的結果為 nan，怎麼辦？

可以嘗試調低 learning rate \(\alpha\)。

決策樹（decision tree）（四）——缺失值處理

【2.5】缺失值的處理 - SAM'S NOTE

Matrix Factorization: A Simple Tutorial and Implementation in Python

作者：wuliytTaotao

出處：https://www.cnblogs.com/wuliytTaotao/

本作品采用知識共享署名-非商業性使用-相同方式共享 4.0 國際許可協定進行許可，歡迎轉載，但未經作者同意必須保留此段聲明，且在文章頁面明顯位置給出原文連接配接。