非负矩阵分解

from numpy import *
from random import *
def difcost(a,b):
    dif=0
    for i in range(shape(a)[0]):
        for j in range(shape(a)[1]):
            dif+=pow(a[i,j]-b[i,j],2)
    return dif

def factorize(v,pc=10,iter=50):
    ic=shape(v)[0]
    fc=shape(v)[1]
    
    w=matrix([[random() for j in range(pc)] for i in range(ic)])
    h=matrix([[random() for i in range(fc)] for j in range(pc)])
    
    for i in range(iter):
        wh=w*h
        cost=difcost(v,wh)
        print cost
        if cost==0:break
        hn=(transpose(w)*v)
        hd=(transpose(w)*w*h)
        h=matrix(array(h)*array(hn)/array(hd))
        wn=(v*transpose(h))
        wd=(w*h*transpose(h))
        w=matrix(array(w)*array(wn)/array(wd))
    return w,h

if __name__=='__main__':
    m1=matrix([[1,2,3],[4,5,6]])
    m2=matrix([[1,2],[3,4],[5,6]])
    print m1*m2
    w,h=factorize(m1*m2,pc=2,iter=200)
    print w*h
    print w
    print h