import numpy as np

X = np.array([[1, 1], [2, 1], [3, 1.2], [4, 1], [5, 0.8], [6, 1]])

n_components, n_samples, n_features,  = (2,) + X.shape
W = np.random.uniform(size = (n_samples, n_components))
H = np.random.uniform(size = (n_components, n_features))

eps = 1e-4

# NMF
for i in range(100):
    # update B
    A = X.T.dot(W)
    B = W.T.dot(W)
    for j in range(n_components):
        tmp = H[j, :] + A[:, j] - H.T.dot(B[:, j])
        H[j, :] = np.maximum(tmp, eps)

    # update A
    C = X.dot(H.T)
    D = H.dot(H.T)
    for j in range(n_components):
        tmp = W[:, j] * D[j, j] + C[:, j] - W.dot(D[:, j])
        W[:, j] = np.maximum(tmp, eps)
        norm = np.linalg.norm(W[:, j])
        if norm > 0:
            W[:, j] /= norm

print(W)
print(H)
print(W.dot(H))