Mercurial > hg > plosone_underreview
comparison scripts/nmftools.py @ 9:c4841876a8ff branch-tests
adding notebooks and trying to explain classifier coefficients
author | Maria Panteli <m.x.panteli@gmail.com> |
---|---|
date | Mon, 11 Sep 2017 19:06:40 +0100 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
8:0f3eba42b425 | 9:c4841876a8ff |
---|---|
1 # -*- coding: utf-8 -*- | |
2 """ | |
3 Created on Wed Oct 26 12:46:13 2016 | |
4 | |
5 @author: https://github.com/keik/nmftools/blob/master/nmftools/core.py | |
6 """ | |
7 | |
8 import numpy as np | |
9 | |
10 | |
11 def nmf(Y, R=3, n_iter=50, init_H=[], init_U=[], verbose=False): | |
12 """ | |
13 decompose non-negative matrix to components and activation with NMF | |
14 Y ≈ HU | |
15 Y ∈ R (m, n) | |
16 H ∈ R (m, k) | |
17 HU ∈ R (k, n) | |
18 parameters | |
19 ---- | |
20 Y: target matrix to decompose | |
21 R: number of bases to decompose | |
22 n_iter: number for executing objective function to optimize | |
23 init_H: initial value of H matrix. default value is random matrix | |
24 init_U: initial value of U matrix. default value is random matrix | |
25 return | |
26 ---- | |
27 Array of: | |
28 0: matrix of H | |
29 1: matrix of U | |
30 2: array of cost transition | |
31 """ | |
32 | |
33 eps = np.spacing(1) | |
34 | |
35 # size of input spectrogram | |
36 M = Y.shape[0] | |
37 N = Y.shape[1] | |
38 | |
39 # initialization | |
40 if len(init_U): | |
41 U = init_U | |
42 R = init_U.shape[0] | |
43 else: | |
44 U = np.random.rand(R,N); | |
45 | |
46 if len(init_H): | |
47 H = init_H; | |
48 R = init_H.shape[1] | |
49 else: | |
50 H = np.random.rand(M,R) | |
51 | |
52 # array to save the value of the euclid divergence | |
53 cost = np.zeros(n_iter) | |
54 | |
55 # computation of Lambda (estimate of Y) | |
56 Lambda = np.dot(H, U) | |
57 | |
58 # iterative computation | |
59 for i in range(n_iter): | |
60 | |
61 # compute euclid divergence | |
62 cost[i] = euclid_divergence(Y, Lambda) | |
63 | |
64 # update H | |
65 H *= np.dot(Y, U.T) / (np.dot(np.dot(H, U), U.T) + eps) | |
66 | |
67 # update U | |
68 U *= np.dot(H.T, Y) / (np.dot(np.dot(H.T, H), U) + eps) | |
69 | |
70 # recomputation of Lambda | |
71 Lambda = np.dot(H, U) | |
72 | |
73 return [H, U, cost] | |
74 | |
75 | |
76 def ssnmf(Y, R=3, n_iter=50, F=[], init_G=[], init_H=[], init_U=[], verbose=False): | |
77 """ | |
78 decompose non-negative matrix to components and activation with semi-supervised NMF | |
79 Y ≈ FG + HU | |
80 Y ∈ R (m, n) | |
81 F ∈ R (m, x) | |
82 G ∈ R (x, n) | |
83 H ∈ R (m, k) | |
84 U ∈ R (k, n) | |
85 parameters | |
86 ---- | |
87 Y: target matrix to decompose | |
88 R: number of bases to decompose | |
89 n_iter: number for executing objective function to optimize | |
90 F: matrix as supervised base components | |
91 init_W: initial value of W matrix. default value is random matrix | |
92 init_H: initial value of W matrix. default value is random matrix | |
93 return | |
94 ---- | |
95 Array of: | |
96 0: matrix of F | |
97 1: matrix of G | |
98 2: matrix of H | |
99 3: matrix of U | |
100 4: array of cost transition | |
101 """ | |
102 | |
103 eps = np.spacing(1) | |
104 | |
105 # size of input spectrogram | |
106 M = Y.shape[0]; | |
107 N = Y.shape[1]; | |
108 X = F.shape[1] | |
109 | |
110 # initialization | |
111 if len(init_G): | |
112 G = init_G | |
113 X = init_G.shape[1] | |
114 else: | |
115 G = np.random.rand(X, N) | |
116 | |
117 if len(init_U): | |
118 U = init_U | |
119 R = init_U.shape[0] | |
120 else: | |
121 U = np.random.rand(R, N) | |
122 | |
123 if len(init_H): | |
124 H = init_H | |
125 R = init_H.shape[1] | |
126 else: | |
127 H = np.random.rand(M, R) | |
128 | |
129 # array to save the value of the euclid divergence | |
130 cost = np.zeros(n_iter) | |
131 | |
132 # computation of Lambda (estimate of Y) | |
133 Lambda = np.dot(F, G) + np.dot(H, U) | |
134 | |
135 # iterative computation | |
136 for it in range(n_iter): | |
137 | |
138 # compute euclid divergence | |
139 cost[it] = euclid_divergence(Y, Lambda + eps) | |
140 | |
141 # update of H | |
142 H *= (np.dot(Y, U.T) + eps) / (np.dot(np.dot(H, U) + np.dot(F, G), U.T) + eps) | |
143 | |
144 # update of U | |
145 U *= (np.dot(H.T, Y) + eps) / (np.dot(H.T, np.dot(H, U) + np.dot(F, G)) + eps) | |
146 | |
147 # update of G | |
148 G *= (np.dot(F.T, Y) + eps)[np.arange(G.shape[0])] / (np.dot(F.T, np.dot(H, U) + np.dot(F, G)) + eps) | |
149 | |
150 # recomputation of Lambda (estimate of V) | |
151 Lambda = np.dot(H, U) + np.dot(F, G) | |
152 | |
153 return [F, G, H, U, cost] | |
154 | |
155 | |
156 def euclid_divergence(V, Vh): | |
157 d = 1 / 2 * (V ** 2 + Vh ** 2 - 2 * V * Vh).sum() | |
158 return d |