1 """
2 Non-negative matrix factorization (NNMF).
3
4 Uses algorithms in
5
6 Lee, D. D. and Seung, H. S.
7 "Algorithms for Non-negative Matrix Factorization".
8 http://hebb.mit.edu/people/seung/papers/nmfconverge.pdf
9
10 Instead of simple functions, this example uses classes to simplify
11 the different algorithms.
12
13 """
14
15 import numpy as np
16
18
19 """
20 Use multiplicative
21 updates to find rank k NNMF of X, written as X=WH.
22
23 Cost function is the cost method. Updates are done using the update
24 method.
25 """
26
27 maxit = 500
28
30 n, p = self.X.shape
31 k = self.k
32 return n*p - (k * (p + k - 1) - 1)
33 df = property(_getdf)
34
35 - def __init__(self, X, k=None, initial=None):
36 """
37
38 If k is None, then k=min(X.shape).
39
40 If initial=(W,H) is supplied, it is used a starting point, and k is ignored.
41
42 If normalize, then return W,H,L where X is approximated as WLH with
43 the columns of W normalized to sum to 1, and the rows of H normalized to
44 sum to 1.
45
46 """
47 self.X = X
48 self.k = k
49
51 if initial is not None:
52 W, H = [i.copy() for i in initial]
53 self.k = min(W.shape + H.shape)
54 else:
55 if self.k is None:
56 self.k = min(self.X.shape)
57 s = np.sum(self.X)
58 W = np.random.uniform(size=(self.X.shape[0], self.k))
59 W /= np.sum(W, axis=0)
60 H = np.random.uniform(size=(self.k, self.X.shape[1]))
61 H = (H.T / np.sum(H, axis=1)).T
62 W *= np.sqrt(s / self.k)
63 H *= np.sqrt(s / self.k)
64 return W, H
65
66 - def fit(self, initial=None, tol=1.0e-06,
67 normalize=False):
68
69 l = []
70 W, H = self.initialize(initial=initial)
71
72 self.approx = np.dot(W, H)
73 err = np.nan
74 it = 0
75
76 while True:
77
78 oerr, err = err, self.cost(self.X, self.approx)
79 l.append(err)
80 if np.fabs(oerr - err) / err < tol:
81 break
82 if it >= self.maxit:
83 break
84
85 print it, np.fabs(oerr - err) / err
86 it += 1
87
88 W, H = self.update(W, H)
89
90
91 self.running_cost = l
92 self.W, self.H = W, H
93 if normalize:
94 self.normalize()
95 return self.W, self.H, self.L
96 else:
97 return self.W, self.H
98
100 L = np.ones(self.k)
101 for i in range(self.k):
102 lw = self.W[:,i].sum()
103 L[i] *= lw
104 self.W[:,i] /= lw
105
106 lh = self.H[i].sum()
107 L[i] *= lh
108 self.H[i] /= lh
109
110 self.L = L
111
112 - def cost(self, X, W, H):
113 raise NotImplementedError
114
116 raise NotImplementedError
117
119 """
120 Uses Frobenius (Euclidean) distance
121 for NNMF.
122 """
123
125 self.approx = np.dot(W, H)
126
127 A = np.dot(W.T, self.X)
128 B = np.dot(W.T, self.approx)
129 H = H * A / B
130
131 self.approx = np.dot(W, H)
132
133 C = np.dot(self.X, H.T)
134 D = np.dot(self.approx, H.T)
135 W = W * C / D
136 return W, H
137
138 - def cost(self, X, approx):
139 return np.sqrt(((X - approx)**2).sum())
140
142 """
143 Uses KL divergence
144 for NNMF.
145 """
146
148 self.approx = np.dot(W, H)
149 A = (self.X / (self.approx + np.equal(self.X, 0)))
150 B = np.dot(W.T, A).T
151 C = np.sum(W, axis=0)
152 H = H * (B / C).T
153
154 self.approx = np.dot(W, H)
155 A = (self.X / (self.approx + np.equal(self.X, 0)))
156 B = np.dot(H, A.T).T
157 C = np.sum(H, axis=1)
158 W = W * B / C
159
160 return W, H
161
162 - def cost(self, X, approx):
163 return (X * np.log(X / approx) - X + approx).sum()
164
166
167 """
168 An EM algorithm for non-negative matrix factorization (NNMF).
169
170 This is EXACTLY the same as using the KL divergence from
171 Lee & Seung.
172
173 The model is:
174
175 X[i,j] ~ Poisson(mu[i,j])
176
177 with X[i,j]'s independent with
178
179 mu[i,j] = sum(W[i,l] * H[l,j], l=1...k)
180
181 where the entries of W and H are non-negative.
182
183 The complete data for this EM is
184
185 (X,X1,...,Xk)
186
187 where
188
189 X = X1 + ... + Xk
190
191 Xl[i,j] ~ Poisson(W[i,l] * H[l,j])
192
193 """
194
196 """
197 E-step:
198
199 Conditional on X, the unobserved matrices X=X1+...+Xk have
200 entries that are independent multinomials.
201
202 Since logL is linear in the Xl[i,j]'s, we just have to compute
203 the conditional expecation of Xl[i,j]|X[i,j], l=1,...,k
204 given our current estimates of the means of the Xl[i,j]'s.
205
206 This results in l 'pseudo' X matrices.
207 """
208
209 pX = np.array([np.multiply.outer(W[:,i], H[i]) for i in range(W.shape[1])])
210 pX /= np.sum(pX, axis=0)
211 pX *= self.X
212 return pX
213
215 """
216 Fits rank one independence models to each pseudo X.
217 """
218 W = []
219 H = []
220 for i in range(pX.shape[0]):
221 n = np.sum(pX[i])
222 w = np.sum(pX[i], axis=1) / n
223 h = np.sum(pX[i], axis=0)
224 W.append(w)
225 H.append(h)
226 return np.array(W).T, np.array(H)
227
228 - def logL(self, mu):
229 """
230 log Likelihood for the Poisson model, ignoring the factorial term.
231 """
232 return (self.X * np.log(mu) - mu).sum()
233
239
240 - def fit(self, normalize=False, initial=None, iter=500):
258
259
260
262
263 l1, l2 = 12000, 10000
264 p11 = np.ones(4, np.float) / 4
265 p12 = np.arange(1,6).astype(np.float)
266 p12 /= p12.sum()
267
268 p1 = np.multiply.outer(p11, p12)
269
270 p21 = np.arange(1,5).astype(np.float)
271 p21 /= p21.sum()
272
273 p22 = np.ones(5, np.float) / 5
274 p2 = np.multiply.outer(p21, p22)
275
276 X = poissonm(p1, l=l1) + poissonm(p2, l=l2)
277 mu = p1 * l1 + p2 * l2
278
279 em = EM(X, k=2)
280 W0, H0 = em.initialize()
281 nnmf = KL(X, k=2)
282 W, H, L = nnmf.fit(normalize=True, initial=(W0,H0))
283
284 em.fit(normalize=True, initial=(W0,H0))
285 return nnmf, mu, em
286
287
288 -def X2(X, approx, k):
289 n, p = approx.shape
290 df = n*p - (k * (p + k - 1) + 1)
291 return ((X - mu)**2 / mu).sum() / df
292
293 if __name__ == "__main__":
294 f, mu, em = main()
295
296
297
298
299
300
301 import pylab
302 pylab.plot(em.ll)
303 pylab.show()
304