Python非负矩阵分解，如何处理零值和缺失数据？
我在找一个可以用Python操作的NMF（非负矩阵分解）实现，并且能够处理缺失数据和零值。
我不想在开始分解之前先填补缺失值，我希望这些缺失值在最小化的过程中被忽略。
看起来scikit-learn、nimfa、graphlab和mahout都没有提供这样的选项。
谢谢！
3 个回答
0
下面是一个基于梯度下降法和拟牛顿法的相对快速且简单的解决方案：
import numpy.random as rd
from autograd_minimize import minimize
import numpy as np
import tensorflow as tf
def low_rank_matrix_factorization(X: np.array, rank: int = 5, l2_reg: float = 0.,
positive: bool = False, opt_kwargs:
dict = {'method': 'L-BFGS-B', 'tol': 1e-12}) -> np.array:
"""Factorises a matrix X into the product of two matrices.
The matrix can have nans.
:param X: Input matrix
:type X: np.array
:param rank: rank of matrix decomposition, defaults to 5
:type rank: int, optional
:param l2_reg: l2 regularisation for the submatrices, defaults to 0.
:type l2_reg: float, optional
:param positive: if true, the matrices must have positive coefficients, defaults to False
:type positive: bool, optional
:param opt_kwargs: parameters for the optimizer, defaults to {'method': 'L-BFGS-B', 'tol': 1e-12}
:type opt_kwargs: dict, optional
:return: completed matrix
:rtype: np.array
"""
mask = tf.constant(~np.isnan(X), dtype=tf.float32)
X_ = np.nan_to_num(X.copy(), 0)
X_t = tf.constant(X_, dtype=tf.float32)
Npos = tf.reduce_sum(mask)
def model(U=None, V=None):
return tf.reduce_sum(((U @ V) - X_t)**2 * mask)/Npos + (tf.reduce_mean(U**2) + tf.reduce_mean(V**2)) * l2_reg
x0 = {'U': rd.random((X_t.shape[0], rank)),
'V': rd.random((rank, X_t.shape[1]))}
if positive:
opt_kwargs['bounds'] = {'U': (0, None), 'V': (0, None)}
res = minimize(model, x0, **opt_kwargs, backend='tf')
return np.where(np.isnan(X), res.x['U'] @ res.x['V'], X)
11
Scipy有一个方法可以解决非负最小二乘问题（NNLS）。在这个回答中，我将重现我在博客上关于使用scipy的NNLS进行非负矩阵分解的内容。你可能也会对我其他的博客文章感兴趣，这些文章使用了autograd、Tensorflow和CVXPY来进行NNMF。
目标：我们的目标是给定一个矩阵A，将其分解为两个非负因子，如下所示：
A (M×N) ≈ W (M×K) × H (K×N), such that W (M×K) ≥ 0 and H (K×N) ≥ 0
概述
我们的解决方案分为两个步骤。首先，我们固定W，学习H，给定A。接下来，我们固定H，学习W，给定A。我们会反复进行这个过程。固定一个变量并学习另一个变量（在这种情况下）被称为交替最小二乘法，因为这个问题被简化为一个最小二乘问题。然而，需要注意的是，由于我们希望W和H都是非负的，所以我们使用NNLS而不是普通的最小二乘法。
步骤1：在给定A和W的情况下学习H
通过上面的图示，我们可以使用A的相应列和矩阵W来学习H的每一列。
H[:,j]=NNLS(W,A[:,j])
处理A中的缺失项
在协同过滤的问题中，A通常是用户-物品矩阵，并且有很多缺失项。这些缺失项对应的是没有对物品进行评分的用户。我们可以修改我们的公式来考虑这些缺失项。假设A中有M' ≤ M个条目是有观察数据的，我们现在将上述方程修改为：
H[:,j]=NNLS(W[mask],A[:,j][mask])
其中，掩码是通过仅考虑M′条目来找到的。
步骤2：在给定A和H的情况下学习W
代码示例
导入库
import numpy as np
import pandas as pd
创建要分解的矩阵
M, N = 20, 10
np.random.seed(0)
A_orig = np.abs(np.random.uniform(low=0.0, high=1.0, size=(M,N)))
print pd.DataFrame(A_orig).head()
0 1 2 3 4 5 6 7 8 9
0 0.548814 0.715189 0.602763 0.544883 0.423655 0.645894 0.437587 0.891773 0.963663 0.383442
1 0.791725 0.528895 0.568045 0.925597 0.071036 0.087129 0.020218 0.832620 0.778157 0.870012
2 0.978618 0.799159 0.461479 0.780529 0.118274 0.639921 0.143353 0.944669 0.521848 0.414662
3 0.264556 0.774234 0.456150 0.568434 0.018790 0.617635 0.612096 0.616934 0.943748 0.681820
4 0.359508 0.437032 0.697631 0.060225 0.666767 0.670638 0.210383 0.128926 0.315428 0.363711
掩盖一些条目
A = A_orig.copy()
A[0, 0] = np.NAN
A[3, 1] = np.NAN
A[6, 3] = np.NAN
A_df = pd.DataFrame(A)
print A_df.head()
0 1 2 3 4 5 6 7 8 9
0 NaN 0.715189 0.602763 0.544883 0.423655 0.645894 0.437587 0.891773 0.963663 0.383442
1 0.791725 0.528895 0.568045 0.925597 0.071036 0.087129 0.020218 0.832620 0.778157 0.870012
2 0.978618 0.799159 0.461479 0.780529 0.118274 0.639921 0.143353 0.944669 0.521848 0.414662
3 0.264556 NaN 0.456150 0.568434 0.018790 0.617635 0.612096 0.616934 0.943748 0.681820
4 0.359508 0.437032 0.697631 0.060225 0.666767 0.670638 0.210383 0.128926 0.315428 0.363711
定义矩阵W和H
K = 4
W = np.abs(np.random.uniform(low=0, high=1, size=(M, K)))
H = np.abs(np.random.uniform(low=0, high=1, size=(K, N)))
W = np.divide(W, K*W.max())
H = np.divide(H, K*H.max())
pd.DataFrame(W).head()
0 1 2 3
0 0.078709 0.175784 0.095359 0.045339
1 0.006230 0.016976 0.171505 0.114531
2 0.135453 0.226355 0.250000 0.054753
3 0.167387 0.066473 0.005213 0.191444
4 0.080785 0.096801 0.148514 0.209789
pd.DataFrame(H).head()
0 1 2 3 4 5 6 7 8 9
0 0.074611 0.216164 0.157328 0.003370 0.088415 0.037721 0.250000 0.121806 0.126649 0.162827
1 0.093851 0.034858 0.209333 0.048340 0.130195 0.057117 0.024914 0.219537 0.247731 0.244654
2 0.230833 0.197093 0.084828 0.020651 0.103694 0.059133 0.033735 0.013604 0.184756 0.002910
3 0.196210 0.037417 0.020248 0.022815 0.171121 0.062477 0.107081 0.141921 0.219119 0.185125
定义我们想要最小化的成本
def cost(A, W, H):
from numpy import linalg
WH = np.dot(W, H)
A_WH = A-WH
return linalg.norm(A_WH, 'fro')
然而，由于A中有缺失项，我们必须根据A中存在的条目来定义成本。
def cost(A, W, H):
from numpy import linalg
mask = pd.DataFrame(A).notnull().values
WH = np.dot(W, H)
WH_mask = WH[mask]
A_mask = A[mask]
A_WH_mask = A_mask-WH_mask
# Since now A_WH_mask is a vector, we use L2 instead of Frobenius norm for matrix
return linalg.norm(A_WH_mask, 2)
让我们先看看我们随机分配的W和H的初始值的成本。
cost(A, W, H)
7.3719938519859509
交替NNLS过程
num_iter = 1000
num_display_cost = max(int(num_iter/10), 1)
from scipy.optimize import nnls
for i in range(num_iter):
if i%2 ==0:
# Learn H, given A and W
for j in range(N):
mask_rows = pd.Series(A[:,j]).notnull()
H[:,j] = nnls(W[mask_rows], A[:,j][mask_rows])[0]
else:
for j in range(M):
mask_rows = pd.Series(A[j,:]).notnull()
W[j,:] = nnls(H.transpose()[mask_rows], A[j,:][mask_rows])[0]
WH = np.dot(W, H)
c = cost(A, W, H)
if i%num_display_cost==0:
print i, c
0 4.03939072472
100 2.38059096458
200 2.35814781954
300 2.35717011529
400 2.35711130357
500 2.3571079918
600 2.35710729854
700 2.35710713129
800 2.35710709085
900 2.35710708109
A_pred = pd.DataFrame(np.dot(W, H))
A_pred.head()
0 1 2 3 4 5 6 7 8 9
0 0.564235 0.677712 0.558999 0.631337 0.536069 0.621925 0.629133 0.656010 0.839802 0.545072
1 0.788734 0.539729 0.517534 1.041272 0.119894 0.448402 0.172808 0.658696 0.493093 0.825311
2 0.749886 0.575154 0.558981 0.931156 0.270149 0.502035 0.287008 0.656178 0.588916 0.741519
3 0.377419 0.743081 0.370408 0.637094 0.071684 0.529433 0.767696 0.628507 0.832910 0.605742
4 0.458661 0.327143 0.610012 0.233134 0.685559 0.377750 0.281483 0.269960 0.468756 0.114950
让我们查看掩盖条目的值。
A_pred.values[~pd.DataFrame(A).notnull().values]
array([ 0.56423481, 0.74308143, 0.10283106])
Original values were:
A_orig[~pd.DataFrame(A).notnull().values]
array([ 0.5488135 , 0.77423369, 0.13818295])
11
我使用了这个从Matlab到Python的代码转换表，成功把NMF从Matlab工具箱库中重写了。
我需要对一个大小为40k X 1k的矩阵进行分解，这个矩阵的稀疏度是0.7%。使用500个潜在特征，我的机器花了20分钟完成100次迭代。
这里是方法：
import numpy as np
from scipy import linalg
from numpy import dot
def nmf(X, latent_features, max_iter=100, error_limit=1e-6, fit_error_limit=1e-6):
"""
Decompose X to A*Y
"""
eps = 1e-5
print 'Starting NMF decomposition with {} latent features and {} iterations.'.format(latent_features, max_iter)
X = X.toarray() # I am passing in a scipy sparse matrix
# mask
mask = np.sign(X)
# initial matrices. A is random [0,1] and Y is A\X.
rows, columns = X.shape
A = np.random.rand(rows, latent_features)
A = np.maximum(A, eps)
Y = linalg.lstsq(A, X)[0]
Y = np.maximum(Y, eps)
masked_X = mask * X
X_est_prev = dot(A, Y)
for i in range(1, max_iter + 1):
# ===== updates =====
# Matlab: A=A.*(((W.*X)*Y')./((W.*(A*Y))*Y'));
top = dot(masked_X, Y.T)
bottom = (dot((mask * dot(A, Y)), Y.T)) + eps
A *= top / bottom
A = np.maximum(A, eps)
# print 'A', np.round(A, 2)
# Matlab: Y=Y.*((A'*(W.*X))./(A'*(W.*(A*Y))));
top = dot(A.T, masked_X)
bottom = dot(A.T, mask * dot(A, Y)) + eps
Y *= top / bottom
Y = np.maximum(Y, eps)
# print 'Y', np.round(Y, 2)
# ==== evaluation ====
if i % 5 == 0 or i == 1 or i == max_iter:
print 'Iteration {}:'.format(i),
X_est = dot(A, Y)
err = mask * (X_est_prev - X_est)
fit_residual = np.sqrt(np.sum(err ** 2))
X_est_prev = X_est
curRes = linalg.norm(mask * (X - X_est), ord='fro')
print 'fit residual', np.round(fit_residual, 4),
print 'total residual', np.round(curRes, 4)
if curRes < error_limit or fit_residual < fit_error_limit:
break
return A, Y
在这里，我使用了Scipy的稀疏矩阵作为输入，并且把缺失的值用0替换，方法是使用toarray()。因此，我使用numpy.sign()函数创建了一个掩码。不过，如果你有nan值的话，也可以用numpy.isnan()函数得到相同的结果。