prepare structure for R package
[valse.git] / src / sources / EMGrank.c
1 #include "EMGrank.h"
2 #include <gsl/gsl_linalg.h>
3
4 // Compute pseudo-inverse of a square matrix
5 static Real* pinv(const Real* matrix, mwSize dim)
6 {
7 gsl_matrix* U = gsl_matrix_alloc(dim,dim);
8 gsl_matrix* V = gsl_matrix_alloc(dim,dim);
9 gsl_vector* S = gsl_vector_alloc(dim);
10 gsl_vector* work = gsl_vector_alloc(dim);
11 Real EPS = 1e-10; //threshold for singular value "== 0"
12
13 //copy matrix into U
14 for (mwSize i=0; i<dim*dim; i++)
15 U->data[i] = matrix[i];
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 (mwSize i=0; i<dim; i++)
24 {
25 for (mwSize ii=0; ii<dim; ii++)
26 {
27 Real dotProduct = 0.0;
28 for (mwSize 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(
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 mwSize n, // taille de l'echantillon
56 mwSize p, // nombre de covariables
57 mwSize m, // taille de Y (multivarié)
58 mwSize 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 mwSize 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 for (mwSize i=0; i<p*m*k; i++)
87 phi[i] = 0.0;
88
89 while (ite<mini || (ite<maxi && sumDeltaPhi>tau))
90 {
91 /////////////
92 // Etape M //
93 /////////////
94
95 //M step: Mise à jour de Beta (et donc phi)
96 for (mwSize r=0; r<k; r++)
97 {
98 //Compute Xr = X(Z==r,:) and Yr = Y(Z==r,:)
99 mwSize cardClustR=0;
100 for (mwSize i=0; i<n; i++)
101 {
102 if (Z[i] == r)
103 {
104 for (mwSize j=0; j<p; j++)
105 Xr[cardClustR*p+j] = X[i*p+j];
106 for (mwSize j=0; j<m; j++)
107 Yr[cardClustR*m+j] = Y[i*m+j];
108 cardClustR++;
109 }
110 }
111 if (cardClustR == 0)
112 continue;
113
114 //Compute tXrXr = t(Xr) * Xr
115 for (mwSize j=0; j<p; j++)
116 {
117 for (mwSize jj=0; jj<p; jj++)
118 {
119 Real dotProduct = 0.0;
120 for (mwSize u=0; u<cardClustR; u++)
121 dotProduct += Xr[u*p+j] * Xr[u*p+jj];
122 tXrXr[j*p+jj] = dotProduct;
123 }
124 }
125
126 //Get pseudo inverse = (t(Xr)*Xr)^{-1}
127 Real* invTXrXr = pinv(tXrXr, p);
128
129 // Compute tXrYr = t(Xr) * Yr
130 for (mwSize j=0; j<p; j++)
131 {
132 for (mwSize jj=0; jj<m; jj++)
133 {
134 Real dotProduct = 0.0;
135 for (mwSize u=0; u<cardClustR; u++)
136 dotProduct += Xr[u*p+j] * Yr[u*m+jj];
137 tXrYr[j*m+jj] = dotProduct;
138 }
139 }
140
141 //Fill matrixM with inverse * tXrYr = (t(Xr)*Xr)^{-1} * t(Xr) * Yr
142 for (mwSize j=0; j<p; j++)
143 {
144 for (mwSize jj=0; jj<m; jj++)
145 {
146 Real dotProduct = 0.0;
147 for (mwSize u=0; u<p; u++)
148 dotProduct += invTXrXr[j*p+u] * tXrYr[u*m+jj];
149 matrixM->data[j*m+jj] = dotProduct;
150 }
151 }
152 free(invTXrXr);
153
154 //U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr
155 gsl_linalg_SV_decomp(matrixM, V, S, work);
156
157 //Set m-rank(r) singular values to zero, and recompose
158 //best rank(r) approximation of the initial product
159 for (mwSize j=rank[r]; j<m; j++)
160 S->data[j] = 0.0;
161
162 //[intermediate step] Compute hatBetaR = U * S * t(V)
163 Real* U = matrixM->data;
164 for (mwSize j=0; j<p; j++)
165 {
166 for (mwSize jj=0; jj<m; jj++)
167 {
168 Real dotProduct = 0.0;
169 for (mwSize u=0; u<m; u++)
170 dotProduct += U[j*m+u] * S->data[u] * V->data[jj*m+u];
171 hatBetaR[j*m+jj] = dotProduct;
172 }
173 }
174
175 //Compute phi(:,:,r) = hatBetaR * Rho(:,:,r)
176 for (mwSize j=0; j<p; j++)
177 {
178 for (mwSize jj=0; jj<m; jj++)
179 {
180 Real dotProduct=0.0;
181 for (mwSize u=0; u<m; u++)
182 dotProduct += hatBetaR[j*m+u] * Rho[u*m*k+jj*k+r];
183 phi[j*m*k+jj*k+r] = dotProduct;
184 }
185 }
186 }
187
188 /////////////
189 // Etape E //
190 /////////////
191
192 Real sumLogLLF2 = 0.0;
193 for (mwSize i=0; i<n; i++)
194 {
195 Real sumLLF1 = 0.0;
196 Real maxLogGamIR = -INFINITY;
197 for (mwSize r=0; r<k; r++)
198 {
199 //Compute
200 //Gam(i,r) = Pi(r) * det(Rho(:,:,r)) * exp( -1/2 * (Y(i,:)*Rho(:,:,r) - X(i,:)...
201 // *phi(:,:,r)) * transpose( Y(i,:)*Rho(:,:,r) - X(i,:)*phi(:,:,r) ) );
202 //split in several sub-steps
203
204 //compute det(Rho(:,:,r)) [TODO: avoid re-computations]
205 for (mwSize j=0; j<m; j++)
206 {
207 for (mwSize jj=0; jj<m; jj++)
208 matrixE->data[j*m+jj] = Rho[j*m*k+jj*k+r];
209 }
210 gsl_linalg_LU_decomp(matrixE, permutation, &signum);
211 Real detRhoR = gsl_linalg_LU_det(matrixE, signum);
212
213 //compute Y(i,:)*Rho(:,:,r)
214 for (mwSize j=0; j<m; j++)
215 {
216 YiRhoR[j] = 0.0;
217 for (mwSize u=0; u<m; u++)
218 YiRhoR[j] += Y[i*m+u] * Rho[u*m*k+j*k+r];
219 }
220
221 //compute X(i,:)*phi(:,:,r)
222 for (mwSize j=0; j<m; j++)
223 {
224 XiPhiR[j] = 0.0;
225 for (mwSize u=0; u<p; u++)
226 XiPhiR[j] += X[i*p+u] * phi[u*m*k+j*k+r];
227 }
228
229 //compute dotProduct < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) >
230 Real dotProduct = 0.0;
231 for (mwSize u=0; u<m; u++)
232 dotProduct += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]);
233 Real logGamIR = log(Pi[r]) + log(detRhoR) - 0.5*dotProduct;
234
235 //Z(i) = index of max (gam(i,:))
236 if (logGamIR > maxLogGamIR)
237 {
238 Z[i] = r;
239 maxLogGamIR = logGamIR;
240 }
241 sumLLF1 += exp(logGamIR) / pow(2*M_PI,m/2.0);
242 }
243
244 sumLogLLF2 += log(sumLLF1);
245 }
246
247 // Assign output variable LLF
248 *LLF = -invN * sumLogLLF2;
249
250 //newDeltaPhi = max(max((abs(phi-Phi))./(1+abs(phi))));
251 Real newDeltaPhi = 0.0;
252 for (mwSize j=0; j<p; j++)
253 {
254 for (mwSize jj=0; jj<m; jj++)
255 {
256 for (mwSize r=0; r<k; r++)
257 {
258 Real tmpDist = fabs(phi[j*m*k+jj*k+r]-Phi[j*m*k+jj*k+r])
259 / (1.0+fabs(phi[j*m*k+jj*k+r]));
260 if (tmpDist > newDeltaPhi)
261 newDeltaPhi = tmpDist;
262 }
263 }
264 }
265
266 //update distance parameter to check algorithm convergence (delta(phi, Phi))
267 //TODO: deltaPhi should be a linked list for perf.
268 if (ite < deltaPhiBufferSize)
269 deltaPhi[ite] = newDeltaPhi;
270 else
271 {
272 sumDeltaPhi -= deltaPhi[0];
273 for (int u=0; u<deltaPhiBufferSize-1; u++)
274 deltaPhi[u] = deltaPhi[u+1];
275 deltaPhi[deltaPhiBufferSize-1] = newDeltaPhi;
276 }
277 sumDeltaPhi += newDeltaPhi;
278
279 // update other local variables
280 for (mwSize j=0; j<m; j++)
281 {
282 for (mwSize jj=0; jj<p; jj++)
283 {
284 for (mwSize r=0; r<k; r++)
285 Phi[j*m*k+jj*k+r] = phi[j*m*k+jj*k+r];
286 }
287 }
288 ite++;
289 }
290
291 //free memory
292 free(hatBetaR);
293 free(deltaPhi);
294 free(Phi);
295 gsl_matrix_free(matrixE);
296 gsl_matrix_free(matrixM);
297 gsl_permutation_free(permutation);
298 gsl_vector_free(work);
299 gsl_matrix_free(V);
300 gsl_vector_free(S);
301 free(XiPhiR);
302 free(YiRhoR);
303 free(Xr);
304 free(Yr);
305 free(tXrXr);
306 free(tXrYr);
307 free(Z);
308 }