La factorisation matricielle non négative ne converge pas

Question

J'essaie de mettre en œuvre une factorisation matricielle non négative en utilisant la divergence de Kullback-Liebler comme mesure de similarité. L'algorithme est décrit dans: http://hebb.mit.edu/people/seung/papers/nmfconverge.pdf. Voici ma mise en œuvre python/numpy, avec un exemple de matrice pour l'exécuter. En un mot, l'algorithme est supposé apprendre les matrices W (n par r) et H (r par m) de sorte que V (n par m) soit approximativement WH. Vous commencez avec des valeurs aléatoires dans W et H, et en suivant les règles de mise à jour décrites dans le document Seung and Lee, vous êtes supposé vous approcher de plus en plus des bonnes approximations pour W et H.

réduire de façon monotone la mesure de divergence, mais ce n'est pas ce qui se passe dans ma mise en œuvre. Au lieu de cela, il s'installe dans une alternance entre deux valeurs de divergence. Si vous regardez W et H, vous pouvez voir que la factorisation résultante n'est pas particulièrement bonne. Je me suis demandé si utiliser l'ancien H mis à jour lors du calcul de la mise à jour pour W. J'ai essayé dans les deux sens, et cela ne change pas le comportement de l'implémentation.

J'ai vérifié mon implémentation contre le papier plusieurs fois, et je ne vois pas ce que je fais de mal. Quelqu'un peut-il faire la lumière sur la question?

import numpy as np 

def update(V, W, H, r, n, m): 
    n,m = V.shape 
    WH = W.dot(H) 

    # equation (5) 
    H_coeff = np.zeros(H.shape) 
    for a in range(r): 
     for mu in range(m): 
      for i in range(n): 
       H_coeff[a, mu] += W[i, a] * V[i, mu]/WH[i, mu] 
      H_coeff[a, mu] /= sum(W)[a] 
    H = H * H_coeff 

    W_coeff = np.zeros(W.shape) 
    for i in range(n): 
     for a in range(r): 
      for mu in range(m): 
       W_coeff[i, a] += H[a, mu] * V[i, mu]/WH[i, mu] 
      W_coeff[i, a] /= sum(H.T)[a] 
    W = W * W_coeff 

    return W, H 


def factor(V, r, iterations=100): 
    n, m = V.shape 
    avg_V = sum(sum(V))/n/m 
    W = np.random.random(n*r).reshape(n,r)*avg_V 
    H = np.random.random(r*m).reshape(r,m)*avg_V 

    for i in range(iterations): 
     WH = W.dot(H) 
     divergence = sum(sum(V * np.log(V/WH) - V + WH)) # equation (3) 
     print "At iteration " + str(i) + ", the Kullback-Liebler divergence is", divergence 
     W,H = update(V, W, H, r, n, m) 

    return W, H 


V = np.arange(0.01,1.01,0.01).reshape(10,10) 

W, H = factor(V, 6) 

Answer 1

Comment éliminer l'effet de l'alternance:

La dernière ligne de la démonstration du théorème 2 lit,

En inversant les rôles de H et W, la règle de mise à jour pour W, il peut être montré que n'est pas croissante. Nous pouvons donc supposer que la mise à jour H peut être effectuée indépendamment de la mise à jour W.

Cela signifie que la mise à jour après H:

H = H * H_coeff 

nous devons également mettre à jour la valeur intermédiaire WH avant la mise à jour W:

WH = W.dot(H) 
W = W * W_coeff 

Les deux mises à jour diminuent la divergence. Essayez-le: Coller WH = W.dot(H) avant le calcul pour W_coeff, et l'effet d'alternance disparaît.

---

Simplifier le code:

Lorsque vous traitez avec des tableaux numpy, utiliser leurs méthodes mean et sum, et évitez d'utiliser la fonction Python sum:

avg_V = sum(sum(V))/n/m 

peut être écrit comme

avg_V = V.mean() 

et

divergence = sum(sum(V * np.log(V/WH) - V + WH)) # equation (3) 

peut être écrit comme

divergence = ((V * np.log(V_over_WH)) - V + WH).sum() 

Évitez la fonction builtin Python sum parce

• il est plus lent que la méthode NumPy sum et
• il n'est pas aussi polyvalent que la méthode NumPy sum. (Il ne vous permet pas de spécifier l'axe sur lequel la somme. Nous avons réussi à éliminer deux appels à Python sum avec un appel à la méthode sum ou mean de NumPy.)

---

éliminons le triple -loop:

Mais une plus grande amélioration à la fois la vitesse et la lisibilité peut être obtenue en remplaçant

H_coeff = np.zeros(H.shape) 
for a in range(r): 
    for mu in range(m): 
     for i in range(n): 
      H_coeff[a, mu] += W[i, a] * V[i, mu]/WH[i, mu] 
     H_coeff[a, mu] /= sum(W)[a] 
H = H * H_coeff 

avec

V_over_WH = V/WH 
H *= (np.dot(V_over_WH.T, W)/W.sum(axis=0)).T 

---

Explication:

Si vous regardez l'équation 5 règle de mise à jour pour H, premier avis que les indices pour V et (W H) sont identiques. Ainsi, vous pouvez remplacer V/(W H) avec

V_over_WH = V/WH 

Ensuite, notons que dans le numérateur, nous sommant sur l'indice i, qui est le premier indice dans les deux W et V_over_WH. Nous pouvons exprimer que la multiplication matrice:

np.dot(V_over_WH.T, W).T 

Et le dénominateur est simplement:

W.sum(axis=0).T 

Si l'on divise le numérateur et le dénominateur

(np.dot(V_over_WH.T, W)/W.sum(axis=0)).T 

nous obtenons une matrice indexée par les deux indices restants, alpha et mu, dans cet ordre.C'est la même chose que les indices pour H. Nous voulons donc multiplier H par ce ratio par élément. Parfait. NumPy multiplie les tableaux par élément par défaut.

Ainsi, nous pouvons exprimer la règle entière de mise à jour pour H comme

H *= (np.dot(V_over_WH.T, W)/W.sum(axis=0)).T 

---

Donc, mettre tous ensemble:

import numpy as np 
np.random.seed(1) 


def update(V, W, H, WH, V_over_WH): 
    # equation (5) 
    H *= (np.dot(V_over_WH.T, W)/W.sum(axis=0)).T 

    WH = W.dot(H) 
    V_over_WH = V/WH 
    W *= np.dot(V_over_WH, H.T)/H.sum(axis=1) 

    WH = W.dot(H) 
    V_over_WH = V/WH 
    return W, H, WH, V_over_WH 


def factor(V, r, iterations=100): 
    n, m = V.shape 
    avg_V = V.mean() 
    W = np.random.random(n * r).reshape(n, r) * avg_V 
    H = np.random.random(r * m).reshape(r, m) * avg_V 
    WH = W.dot(H) 
    V_over_WH = V/WH 

    for i in range(iterations): 
     W, H, WH, V_over_WH = update(V, W, H, WH, V_over_WH) 
     # equation (3) 
     divergence = ((V * np.log(V_over_WH)) - V + WH).sum() 
     print("At iteration {i}, the Kullback-Liebler divergence is {d}".format(
      i=i, d=divergence)) 
    return W, H 

V = np.arange(0.01, 1.01, 0.01).reshape(10, 10) 
# V = np.arange(1,101).reshape(10,10).astype('float') 
W, H = factor(V, 6)