矩阵补全（Matrix Completion）和缺失值预处理

目录

矩阵补全（Matrix Completion），就是补上一个含缺失值矩阵的缺失部分。

矩阵补全可以通过矩阵分解（matrix factorization）将一个含缺失值的矩阵 X 分解为两个（或多个）矩阵，然后这些分解后的矩阵相乘就可以得到原矩阵的近似 X'，我们用这个近似矩阵 X' 的值来填补原矩阵 X 的缺失部分。

矩阵补全有很多方面的应用，如推荐系统、缺失值预处理。

除了 EM 算法、树模型，机器学习中的大多数算法都需要输入的数据是不含缺失值的。在 deep learning 模型中，通过梯度的计算公式就可以发现，如果 feature 中含有缺失值，那么梯度也会含缺失值，梯度也就未知了。对缺失值的处理是在模型训练开始前就应该完成的，故也称为预处理。

数据缺失在实际场景中不可避免，对于一个包含 \(n\) 个 samples，每个 sample 有 \(m\) 个 features 的数据集 \(D\)，我们可以将该数据集 \(D\) 整理为一个 \(n×m\) 的矩阵 \(X\)。

通过矩阵分解补全矩阵是一种处理缺失值的方式，但在介绍之前，先介绍一些简单常用的缺失值预处理方式。

1 常用的缺失值预处理方式

1.1 不处理

不进行缺失值预处理，缺了就缺了，找一个对缺失值不敏感的算法（如“树模型”），直接训练。

1.2 剔除

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

缺失值较多的行，即一个 sample 的很多 features 都缺失了；缺失值较多的列，即大部分 samples 都没有该 feature。剔除这些 samples 或 features，而不是填充它们，避免引入过多的噪声。

当数据超级多时，我们甚至可以对含有缺失值的样本直接剔除，当剔除的比例不大时，这也完全可以接受。

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

posted @ 2019-05-06 17:06  wuliytTaotao  阅读(18088)  评论(27)  编辑  收藏  举报