228ee602 |
1 | #include <stdlib.h> |
2 | #include <gsl/gsl_linalg.h> |
3 | #include "utils.h" |
4 | |
5 | // Compute pseudo-inverse of a square matrix |
6 | static Real* pinv(const Real* matrix, int dim) |
7 | { |
8 | gsl_matrix* U = gsl_matrix_alloc(dim,dim); |
9 | gsl_matrix* V = gsl_matrix_alloc(dim,dim); |
10 | gsl_vector* S = gsl_vector_alloc(dim); |
11 | gsl_vector* work = gsl_vector_alloc(dim); |
12 | Real EPS = 1e-10; //threshold for singular value "== 0" |
13 | |
14 | //copy matrix into U |
15 | copyArray(matrix, U->data, dim*dim); |
16 | |
17 | //U,S,V = SVD of matrix |
18 | gsl_linalg_SV_decomp(U, V, S, work); |
19 | gsl_vector_free(work); |
20 | |
21 | // Obtain pseudo-inverse by V*S^{-1}*t(U) |
22 | Real* inverse = (Real*)malloc(dim*dim*sizeof(Real)); |
23 | for (int i=0; i<dim; i++) |
24 | { |
25 | for (int ii=0; ii<dim; ii++) |
26 | { |
27 | Real dotProduct = 0.0; |
28 | for (int j=0; j<dim; j++) |
29 | dotProduct += V->data[i*dim+j] * (S->data[j] > EPS ? 1.0/S->data[j] : 0.0) * U->data[ii*dim+j]; |
30 | inverse[i*dim+ii] = dotProduct; |
31 | } |
32 | } |
33 | |
34 | gsl_matrix_free(U); |
35 | gsl_matrix_free(V); |
36 | gsl_vector_free(S); |
37 | return inverse; |
38 | } |
39 | |
40 | // TODO: comment EMGrank purpose |
41 | void EMGrank_core( |
42 | // IN parameters |
43 | const Real* Pi, // parametre de proportion |
44 | const Real* Rho, // parametre initial de variance renormalisé |
45 | int mini, // nombre minimal d'itérations dans l'algorithme EM |
46 | int maxi, // nombre maximal d'itérations dans l'algorithme EM |
47 | const Real* X, // régresseurs |
48 | const Real* Y, // réponse |
49 | Real tau, // seuil pour accepter la convergence |
50 | const int* rank, // vecteur des rangs possibles |
51 | // OUT parameters |
52 | Real* phi, // parametre de moyenne renormalisé, calculé par l'EM |
53 | Real* LLF, // log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres |
54 | // additional size parameters |
55 | int n, // taille de l'echantillon |
56 | int p, // nombre de covariables |
57 | int m, // taille de Y (multivarié) |
58 | int k) // nombre de composantes |
59 | { |
60 | // Allocations, initializations |
61 | Real* Phi = (Real*)calloc(p*m*k,sizeof(Real)); |
62 | Real* hatBetaR = (Real*)malloc(p*m*sizeof(Real)); |
63 | int signum; |
64 | Real invN = 1.0/n; |
65 | int deltaPhiBufferSize = 20; |
66 | Real* deltaPhi = (Real*)malloc(deltaPhiBufferSize*sizeof(Real)); |
67 | int ite = 0; |
68 | Real sumDeltaPhi = 0.0; |
69 | Real* YiRhoR = (Real*)malloc(m*sizeof(Real)); |
70 | Real* XiPhiR = (Real*)malloc(m*sizeof(Real)); |
71 | Real* Xr = (Real*)malloc(n*p*sizeof(Real)); |
72 | Real* Yr = (Real*)malloc(n*m*sizeof(Real)); |
73 | Real* tXrXr = (Real*)malloc(p*p*sizeof(Real)); |
74 | Real* tXrYr = (Real*)malloc(p*m*sizeof(Real)); |
75 | gsl_matrix* matrixM = gsl_matrix_alloc(p, m); |
76 | gsl_matrix* matrixE = gsl_matrix_alloc(m, m); |
77 | gsl_permutation* permutation = gsl_permutation_alloc(m); |
78 | gsl_matrix* V = gsl_matrix_alloc(m,m); |
79 | gsl_vector* S = gsl_vector_alloc(m); |
80 | gsl_vector* work = gsl_vector_alloc(m); |
81 | |
82 | //Initialize class memberships (all elements in class 0; TODO: randomize ?) |
83 | int* Z = (int*)calloc(n, sizeof(int)); |
84 | |
85 | //Initialize phi to zero, because some M loops might exit before phi affectation |
86 | zeroArray(phi, p*m*k); |
87 | |
88 | while (ite<mini || (ite<maxi && sumDeltaPhi>tau)) |
89 | { |
90 | ///////////// |
91 | // Etape M // |
92 | ///////////// |
93 | |
94 | //M step: Mise à jour de Beta (et donc phi) |
95 | for (int r=0; r<k; r++) |
96 | { |
97 | //Compute Xr = X(Z==r,:) and Yr = Y(Z==r,:) |
98 | int cardClustR=0; |
99 | for (int i=0; i<n; i++) |
100 | { |
101 | if (Z[i] == r) |
102 | { |
103 | for (int j=0; j<p; j++) |
104 | Xr[mi(cardClustR,j,n,p)] = X[mi(i,j,n,p)]; |
105 | for (int j=0; j<m; j++) |
106 | Yr[mi(cardClustR,j,n,m)] = Y[mi(i,j,n,m)]; |
107 | cardClustR++; |
108 | } |
109 | } |
110 | if (cardClustR == 0) |
111 | continue; |
112 | |
113 | //Compute tXrXr = t(Xr) * Xr |
114 | for (int j=0; j<p; j++) |
115 | { |
116 | for (int jj=0; jj<p; jj++) |
117 | { |
118 | Real dotProduct = 0.0; |
119 | for (int u=0; u<cardClustR; u++) |
120 | dotProduct += Xr[mi(u,j,n,p)] * Xr[mi(u,jj,n,p)]; |
121 | tXrXr[mi(j,jj,p,p)] = dotProduct; |
122 | } |
123 | } |
124 | |
125 | //Get pseudo inverse = (t(Xr)*Xr)^{-1} |
126 | Real* invTXrXr = pinv(tXrXr, p); |
127 | |
128 | // Compute tXrYr = t(Xr) * Yr |
129 | for (int j=0; j<p; j++) |
130 | { |
131 | for (int jj=0; jj<m; jj++) |
132 | { |
133 | Real dotProduct = 0.0; |
134 | for (int u=0; u<cardClustR; u++) |
135 | dotProduct += Xr[mi(u,j,n,p)] * Yr[mi(u,jj,n,m)]; |
136 | tXrYr[mi(j,jj,p,m)] = dotProduct; |
137 | } |
138 | } |
139 | |
140 | //Fill matrixM with inverse * tXrYr = (t(Xr)*Xr)^{-1} * t(Xr) * Yr |
141 | for (int j=0; j<p; j++) |
142 | { |
143 | for (int jj=0; jj<m; jj++) |
144 | { |
145 | Real dotProduct = 0.0; |
146 | for (int u=0; u<p; u++) |
147 | dotProduct += invTXrXr[mi(j,u,p,p)] * tXrYr[mi(u,jj,p,m)]; |
148 | matrixM->data[j*m+jj] = dotProduct; |
149 | } |
150 | } |
151 | free(invTXrXr); |
152 | |
153 | //U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr |
154 | gsl_linalg_SV_decomp(matrixM, V, S, work); |
155 | |
156 | //Set m-rank(r) singular values to zero, and recompose |
157 | //best rank(r) approximation of the initial product |
158 | for (int j=rank[r]; j<m; j++) |
159 | S->data[j] = 0.0; |
160 | |
161 | //[intermediate step] Compute hatBetaR = U * S * t(V) |
162 | double* U = matrixM->data; //GSL require double precision |
163 | for (int j=0; j<p; j++) |
164 | { |
165 | for (int jj=0; jj<m; jj++) |
166 | { |
167 | Real dotProduct = 0.0; |
168 | for (int u=0; u<m; u++) |
169 | dotProduct += U[j*m+u] * S->data[u] * V->data[jj*m+u]; |
170 | hatBetaR[mi(j,jj,p,m)] = dotProduct; |
171 | } |
172 | } |
173 | |
174 | //Compute phi(:,:,r) = hatBetaR * Rho(:,:,r) |
175 | for (int j=0; j<p; j++) |
176 | { |
177 | for (int jj=0; jj<m; jj++) |
178 | { |
179 | Real dotProduct=0.0; |
180 | for (int u=0; u<m; u++) |
181 | dotProduct += hatBetaR[mi(j,u,p,m)] * Rho[ai(u,jj,r,m,m,k)]; |
182 | phi[ai(j,jj,r,p,m,k)] = dotProduct; |
183 | } |
184 | } |
185 | } |
186 | |
187 | ///////////// |
188 | // Etape E // |
189 | ///////////// |
190 | |
191 | Real sumLogLLF2 = 0.0; |
192 | for (int i=0; i<n; i++) |
193 | { |
194 | Real sumLLF1 = 0.0; |
195 | Real maxLogGamIR = -INFINITY; |
196 | for (int r=0; r<k; r++) |
197 | { |
198 | //Compute |
199 | //Gam(i,r) = Pi(r) * det(Rho(:,:,r)) * exp( -1/2 * (Y(i,:)*Rho(:,:,r) - X(i,:)... |
200 | //*phi(:,:,r)) * transpose( Y(i,:)*Rho(:,:,r) - X(i,:)*phi(:,:,r) ) ); |
201 | //split in several sub-steps |
202 | |
203 | //compute det(Rho(:,:,r)) [TODO: avoid re-computations] |
204 | for (int j=0; j<m; j++) |
205 | { |
206 | for (int jj=0; jj<m; jj++) |
207 | matrixE->data[j*m+jj] = Rho[ai(j,jj,r,m,m,k)]; |
208 | } |
209 | gsl_linalg_LU_decomp(matrixE, permutation, &signum); |
210 | Real detRhoR = gsl_linalg_LU_det(matrixE, signum); |
211 | |
212 | //compute Y(i,:)*Rho(:,:,r) |
213 | for (int j=0; j<m; j++) |
214 | { |
215 | YiRhoR[j] = 0.0; |
216 | for (int u=0; u<m; u++) |
217 | YiRhoR[j] += Y[mi(i,u,n,m)] * Rho[ai(u,j,r,m,m,k)]; |
218 | } |
219 | |
220 | //compute X(i,:)*phi(:,:,r) |
221 | for (int j=0; j<m; j++) |
222 | { |
223 | XiPhiR[j] = 0.0; |
224 | for (int u=0; u<p; u++) |
225 | XiPhiR[j] += X[mi(i,u,n,p)] * phi[ai(u,j,r,p,m,k)]; |
226 | } |
227 | |
228 | //compute dotProduct < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) > |
229 | Real dotProduct = 0.0; |
230 | for (int u=0; u<m; u++) |
231 | dotProduct += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]); |
232 | Real logGamIR = log(Pi[r]) + log(detRhoR) - 0.5*dotProduct; |
233 | |
234 | //Z(i) = index of max (gam(i,:)) |
235 | if (logGamIR > maxLogGamIR) |
236 | { |
237 | Z[i] = r; |
238 | maxLogGamIR = logGamIR; |
239 | } |
240 | sumLLF1 += exp(logGamIR) / pow(2*M_PI,m/2.0); |
241 | } |
242 | |
243 | sumLogLLF2 += log(sumLLF1); |
244 | } |
245 | |
246 | // Assign output variable LLF |
247 | *LLF = -invN * sumLogLLF2; |
248 | |
249 | //newDeltaPhi = max(max((abs(phi-Phi))./(1+abs(phi)))); |
250 | Real newDeltaPhi = 0.0; |
251 | for (int j=0; j<p; j++) |
252 | { |
253 | for (int jj=0; jj<m; jj++) |
254 | { |
255 | for (int r=0; r<k; r++) |
256 | { |
257 | Real tmpDist = fabs(phi[ai(j,jj,r,p,m,k)]-Phi[ai(j,jj,r,p,m,k)]) |
258 | / (1.0+fabs(phi[ai(j,jj,r,p,m,k)])); |
259 | if (tmpDist > newDeltaPhi) |
260 | newDeltaPhi = tmpDist; |
261 | } |
262 | } |
263 | } |
264 | |
265 | //update distance parameter to check algorithm convergence (delta(phi, Phi)) |
266 | //TODO: deltaPhi should be a linked list for perf. |
267 | if (ite < deltaPhiBufferSize) |
268 | deltaPhi[ite] = newDeltaPhi; |
269 | else |
270 | { |
271 | sumDeltaPhi -= deltaPhi[0]; |
272 | for (int u=0; u<deltaPhiBufferSize-1; u++) |
273 | deltaPhi[u] = deltaPhi[u+1]; |
274 | deltaPhi[deltaPhiBufferSize-1] = newDeltaPhi; |
275 | } |
276 | sumDeltaPhi += newDeltaPhi; |
277 | |
278 | // update other local variables |
279 | for (int j=0; j<m; j++) |
280 | { |
281 | for (int jj=0; jj<p; jj++) |
282 | { |
283 | for (int r=0; r<k; r++) |
284 | Phi[ai(jj,j,r,p,m,k)] = phi[ai(jj,j,r,p,m,k)]; |
285 | } |
286 | } |
287 | ite++; |
288 | } |
289 | |
290 | //free memory |
291 | free(hatBetaR); |
292 | free(deltaPhi); |
293 | free(Phi); |
294 | gsl_matrix_free(matrixE); |
295 | gsl_matrix_free(matrixM); |
296 | gsl_permutation_free(permutation); |
297 | gsl_vector_free(work); |
298 | gsl_matrix_free(V); |
299 | gsl_vector_free(S); |
300 | free(XiPhiR); |
301 | free(YiRhoR); |
302 | free(Xr); |
303 | free(Yr); |
304 | free(tXrXr); |
305 | free(tXrYr); |
306 | free(Z); |
307 | } |