矩阵补全（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)

4 经过矩阵分解补全矩阵的一些小问题

4.1 需不须要对 bias 进行正则化？
按照通常 deep learning 模型，是不对 bias 进行正则化的，而本文的代码对 bias 进行了正则化，具体有没有影响不得而知。

4.2 若是出现 "RuntimeWarning: overflow encountered in multiply" 等 Warning 形成最后的结果为 nan，怎么办？
能够尝试调低 learning rate \(\alpha\)。

References

决策树（decision tree）（四）——缺失值处理
 【2.5】缺失值的处理 - SAM'S NOTE
Matrix Factorization: A Simple Tutorial and Implementation in Python