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

目录

• 1 常用的缺失值预处理方式
  • 1.1 不处理
  • 1.2 剔除
  • 1.3 填充
• 2 利用矩阵分解补全缺失值
• 3 矩阵分解补全缺失值代码实现
• 4 通过矩阵分解补全矩阵的一些小问题
• References

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

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

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

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

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

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

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

1.1 不处理

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

1.2 剔除

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

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

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

1.3 填充

1.3.1 简单填充

在矩阵 XX 的每个缺失位置上填上一个数，来代替缺失值。填一个数也不能乱来，如果 feature 代表年龄，那么肯定要填正数；如果 feature 代表性别，那么填 0 或 1 更合适（0 代表男，1 代表女）。

一般有以下几种简单的填充值：（均值和众数都是在一个 feature 下计算，即在矩阵 XX 的每一列中计算均值和众数）

• 填 0
• 填 均值
• 填 众数

1.3.2 建模填充

这种方式通过观察缺失的 feature 和其它已有的 features 之间的联系，建立一个统计模型或者回归模型，然后然后预测缺失 feature 的值应该是什么。

用 EM 算法估计缺失值也可以归为这一类。

当然，常用的缺失值处理方式还有许多，这里就不再列举了。可以看看博客 SAM'S NOTE。

2 利用矩阵分解补全缺失值

如果矩阵 XX 不含缺失值，那么矩阵分解可以将矩阵 XX 分解成两个矩阵 UU (大小 m×km×k )、VV (大小 m×km×k )，其中 k<min{m,n}k<min{m,n} ，则：

 

X=UV⊤X=UV⊤

因为 k<min{m,n}k<min{m,n} ，所以 rank(U)≤krank(U)≤k 、rank(V)≤krank(V)≤k ，该矩阵分解又叫做低秩矩阵分解（low-rank matrix factorization）。

那么为什么 k<min{m,n}k<min{m,n} ？

• 在 samples 和 features 之间存在 k 个关系，每个关系的具体含义不得而知，但如果 k≥min{m,n}k≥min{m,n} ，那么意味着每个 sample 和 feature 之间可以构建一个的关系，而其它的 samples 或者 features 可以和该关系基本无关，体现在矩阵 UU （或 VV ）中就是某一列仅有一个元素不为0，这是不合理的。（参考矩阵分解用在推荐系统方面的解释）
• 当 k 越大，计算量也会越大。

如果矩阵 XX 是完整的，那么矩阵分解 X=UV⊤X=UV⊤ 完全没问题，但现在 XX 中含了缺失值，故没有办法用线性代数的知识直接进行矩阵分解，我们需要一种近似的解法——梯度下降法。

这个时候我们令 X≈X^=UV⊤X≈X^=UV⊤ ，∥X−X^∥2F‖X−X^‖F2 表示含缺失值的原矩阵 XX 和 还原后的近似矩阵 X^X^ 之间误差的平方（Square error），或者称之为 reconstruction error，当然 ∥X−X^∥2F‖X−X^‖F2 的计算只能在不含缺失值的项上。（∥⋅∥F‖⋅‖F 表示 Frobenius norm。）

文献中一般会将 reconstruction error ∥X−X^∥2F‖X−X^‖F2 记为 ∥∥RΩ(X−X^)∥∥2F‖RΩ(X−X^)‖F2 ，其中 [RΩ(X−X^)]ij={xij−x^ij0 if (i,j)∈Ω otherwise [RΩ(X−X^)]ij={xij−x^ij if (i,j)∈Ω0 otherwise  ，ΩΩ 表示非缺失值矩阵元素下标的集合。这里为了简便，直接使用 ∥X−X^∥2F‖X−X^‖F2 ，知道只在不含缺失值的项上计算平方和即可。

我们的目标的是找到矩阵 XX 的近似矩阵 X^X^ ，通过 X^X^ 中对应的值来填充 XX 中缺失的部分。而想要找到 X^X^ ，就是要找到矩阵 UU 和 VV 。当然 X^X^ 要尽可能像 XX ，体现在函数上就是 min∥X−X^∥2Fmin‖X−X^‖F2 。

NOTE：以下矩阵的范数都默认为 Frobenius norm。

Loss function JJ 为:

 

J=∥X−X^∥2=∥X−UV⊤∥2=∑i,j,xij≠nan(xij−∑l=1kuilvjl)2J=‖X−X^‖2=‖X−UV⊤‖2=∑i,j,xij≠nan(xij−∑l=1kuilvjl)2

其中，i,ji,j 分别表示矩阵 XX 的行和列，要求 xij≠nanxij≠nan ，否则没有办法求最小值了。上式中，未知的就是 uil,vjluil,vjl ，也是我们想要求的。

随机初始化矩阵 U,VU,V ，loss function JJ 就可以得到一个误差，基于该误差计算梯度，而想要更新 U,VU,V ，只需要按照梯度下降的公式来即可。
令:

 

eij=xij−∑l=1kuilvjleij=xij−∑l=1kuilvjl

则梯度为：

 

∂J∂uil=−2eijvjl∂J∂vjl=−2eijuil∂J∂uil=−2eijvjl∂J∂vjl=−2eijuil

梯度下降更新公式：

 

uil=uil−α∂J∂uil=uil+2αeijvjlvjl=vjl−α∂J∂vjl=uil+2αeijuiluil=uil−α∂J∂uil=uil+2αeijvjlvjl=vjl−α∂J∂vjl=uil+2αeijuil

算法到这里其实就可以用了，但为了更加完美，可以考虑以下步骤，加入正则项和偏置。

加入正则项

加入正则项，保证矩阵 U,VU,V 中元素不要太大，此时 loss function JJ 如下所示：

 

J=∥X−X^∥2+β2(∥U∥2+∥V∥2)=∑i,j,xij≠nan(xij−∑l=1kuilvjl)2+β2(∑i,lu2il+∑j,lv2jl)J=‖X−X^‖2+β2(‖U‖2+‖V‖2)=∑i,j,xij≠nan(xij−∑l=1kuilvjl)2+β2(∑i,luil2+∑j,lvjl2)

则梯度为：

 

∂J∂uil=−2eijvjl+βuil∂J∂vjl=−2eijuil+βvjl∂J∂uil=−2eijvjl+βuil∂J∂vjl=−2eijuil+βvjl

此时梯度下降更新公式为：

 

uil=uil−α∂J∂uil=uil+α(2eijvjl−βuil)vjl=vjl−α∂J∂vjl=uil+α(2eijuil−βvjl)uil=uil−α∂J∂uil=uil+α(2eijvjl−βuil)vjl=vjl−α∂J∂vjl=uil+α(2eijuil−βvjl)

加入偏置

偏置可以理解为每个样本都有其特性，每个feature也有其特点，故可以加入 bias 来控制。bias 分为三种，第一种是矩阵 XX 整体的的 bias，记为 bb ，那么 b=mean(X)b=mean(X) ，即可以用矩阵 XX 中存在元素的均值来赋值；第二种是 sample 的 bias，记为 b_uib_ui ；第三种是 feature 的 bias，记为 b_vjb_vj 。
则：

 

x^ij=∑l=1kuilvjl+(b+b_ui+b_vj)x^ij=∑l=1kuilvjl+(b+b_ui+b_vj)

其中，b=∑i,j,xij≠nanxijNb=∑i,j,xij≠nanxijN ，NN 表示分子求和元素的个数。
则 loss function JJ 为：

 

J=∥X−X^∥2+β2(∥U∥2+∥V∥2+b_u2+b_v2)=∑i,j,xij≠nan(xij−∑l=1kuilvjl−b−b_ui−b_vj)2+β2(∑i,lu2il+∑j,lv2jl+∑ib_u2i+∑jb_v2j)J=‖X−X^‖2+β2(‖U‖2+‖V‖2+b_u2+b_v2)=∑i,j,xij≠nan(xij−∑l=1kuilvjl−b−b_ui−b_vj)2+β2(∑i,luil2+∑j,lvjl2+∑ib_ui2+∑jb_vj2)

再加入 bias 后，令

 

eij=xij−∑l=1kuilvjl−b−b_ui−b_vjeij=xij−∑l=1kuilvjl−b−b_ui−b_vj

则梯度为：

 

∂J∂uil∂J∂vjl∂J∂b_ui∂J∂b_vj=−2eijvjl+βuil=−2eijuil+βvjl=−2eij+βb_ui=−2eij+βb_vj∂J∂uil=−2eijvjl+βuil∂J∂vjl=−2eijuil+βvjl∂J∂b_ui=−2eij+βb_ui∂J∂b_vj=−2eij+βb_vj

此时梯度下降更新公式为：

 

uilvjlb_uib_vj=uil+α(2eijvjl−βuil)=uil+α(2eijuil−βvjl)=b_ui+α(2eij−βb_ui)=b_vj+α(2eij−βb_vj)uil=uil+α(2eijvjl−βuil)vjl=uil+α(2eijuil−βvjl)b_ui=b_ui+α(2eij−βb_ui)b_vj=b_vj+α(2eij−βb_vj)

3 矩阵分解补全缺失值代码实现

1. import numpy as np

2.  

3. class MF():

4.  

5. def __init__(self, X, k, alpha, beta, iterations):

6. """

7. Perform matrix factorization to predict np.nan entries in a matrix.

8. Arguments

9. - X (ndarray) : sample-feature matrix

10. - k (int) : number of latent dimensions

11. - alpha (float) : learning rate

12. - beta (float) : regularization parameter

13. """

14.  

15. self.X = X

16. self.num_samples, self.num_features = X.shape

17. self.k = k

18. self.alpha = alpha

19. self.beta = beta

20. self.iterations = iterations

21. # True if not nan

22. self.not_nan_index = (np.isnan(self.X) == False)

23.  

24. def train(self):

25. # Initialize factorization matrix U and V

26. self.U = np.random.normal(scale=1./self.k, size=(self.num_samples, self.k))

27. self.V = np.random.normal(scale=1./self.k, size=(self.num_features, self.k))

28.  

29. # Initialize the biases

30. self.b_u = np.zeros(self.num_samples)

31. self.b_v = np.zeros(self.num_features)

32. self.b = np.mean(self.X[np.where(self.not_nan_index)])

33. # Create a list of training samples

34. self.samples = [

35. (i, j, self.X[i, j])

36. for i in range(self.num_samples)

37. for j in range(self.num_features)

38. if not np.isnan(self.X[i, j])

39. ]

40.  

41. # Perform stochastic gradient descent for number of iterations

42. training_process = []

43. for i in range(self.iterations):

44. np.random.shuffle(self.samples)

45. self.sgd()

46. # total square error

47. se = self.square_error()

48. training_process.append((i, se))

49. if (i+1) % 10 == 0:

50. print("Iteration: %d ; error = %.4f" % (i+1, se))

51.  

52. return training_process

53.  

54. def square_error(self):

55. """

56. A function to compute the total square error

57. """

58. predicted = self.full_matrix()

59. error = 0

60. for i in range(self.num_samples):

61. for j in range(self.num_features):

62. if self.not_nan_index[i, j]:

63. error += pow(self.X[i, j] - predicted[i, j], 2)

64. return error

65.  

66. def sgd(self):

67. """

68. Perform stochastic graident descent

69. """

70. for i, j, x in self.samples:

71. # Computer prediction and error

72. prediction = self.get_x(i, j)

73. e = (x - prediction)

74.  

75. # Update biases

76. self.b_u[i] += self.alpha * (2 * e - self.beta * self.b_u[i])

77. self.b_v[j] += self.alpha * (2 * e - self.beta * self.b_v[j])

78.  

79. # Update factorization matrix U and V

80. """

81. If RuntimeWarning: overflow encountered in multiply,

82. then turn down the learning rate alpha.

83. """

84. self.U[i, :] += self.alpha * (2 * e * self.V[j, :] - self.beta * self.U[i,:])

85. self.V[j, :] += self.alpha * (2 * e * self.U[i, :] - self.beta * self.V[j,:])

86.  

87. def get_x(self, i, j):

88. """

89. Get the predicted x of sample i and feature j

90. """

91. prediction = self.b + self.b_u[i] + self.b_v[j] + self.U[i, :].dot(self.V[j, :].T)

92. return prediction

93.  

94. def full_matrix(self):

95. """

96. Computer the full matrix using the resultant biases, U and V

97. """

98. return self.b + self.b_u[:, np.newaxis] + self.b_v[np.newaxis, :] + self.U.dot(self.V.T)

99.  

100. def replace_nan(self, X_hat):

101. """

102. Replace np.nan of X with the corresponding value of X_hat

103. """

104. X = np.copy(self.X)

105. for i in range(self.num_samples):

106. for j in range(self.num_features):

107. if np.isnan(X[i, j]):

108. X[i, j] = X_hat[i, j]

109. return X

110.  

111. if __name__ == '__main__':

112. X = np.array([

113. [5, 3, 0, 1],

114. [4, 0, 0, 1],

115. [1, 1, 0, 5],

116. [1, 0, 0, 4],

117. [0, 1, 5, 4],

118. ], dtype=np.float)

119. # replace 0 with np.nan

120. X[X == 0] = np.nan

121. print(X)

122. # np.random.seed(1)

123. mf = MF(X, k=2, alpha=0.1, beta=0.1, iterations=100)

124. mf.train()

125. X_hat = mf.full_matrix()

126. X_comp = mf.replace_nan(X_hat)

127.  

128. print(X_hat)

129. print(X_comp)

130. print(X)

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

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

4.2 如果出现 "RuntimeWarning: overflow encountered in multiply" 等 Warning 造成最后的结果为 nan，怎么办？
可以尝试调低 learning rate αα 。

References

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

转载于:https://www.cnblogs.com/wuliytTaotao/p/10814770.html