使用这个Matlab到python代码转换表,我能够从Matlab工具箱库重写NMF 。 我不得不分解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
usingtoarray()
方法。因此,使用numpy.sign()
功能创建了蒙版。但是,如果您具有nan
值,则可以使用numpy.isnan()
函数获得相同的结果。