2 #include <gsl/gsl_linalg.h>
5 // Compute pseudo-inverse of a square matrix
6 static float* pinv(const float* matrix
, int dim
)
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 float EPS
= 1e-10; //threshold for singular value "== 0"
15 copyArray(matrix
, U
->data
, dim
*dim
);
17 //U,S,V = SVD of matrix
18 gsl_linalg_SV_decomp(U
, V
, S
, work
);
19 gsl_vector_free(work
);
21 // Obtain pseudo-inverse by V*S^{-1}*t(U)
22 float* inverse
= (float*)malloc(dim
*dim
*sizeof(float));
23 for (int i
=0; i
<dim
; i
++)
25 for (int ii
=0; ii
<dim
; ii
++)
27 float 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
;
40 // TODO: comment EMGrank purpose
43 const float* Pi
, // parametre de proportion
44 const float* 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 float* X
, // régresseurs
48 const float* Y
, // réponse
49 float tau
, // seuil pour accepter la convergence
50 const int* rank
, // vecteur des rangs possibles
52 float* phi
, // parametre de moyenne renormalisé, calculé par l'EM
53 float* 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
60 // Allocations, initializations
61 float* Phi
= (float*)calloc(p
*m
*k
,sizeof(float));
62 float* hatBetaR
= (float*)malloc(p
*m
*sizeof(float));
65 int deltaPhiBufferSize
= 20;
66 float* deltaPhi
= (float*)malloc(deltaPhiBufferSize
*sizeof(float));
68 float sumDeltaPhi
= 0.0;
69 float* YiRhoR
= (float*)malloc(m
*sizeof(float));
70 float* XiPhiR
= (float*)malloc(m
*sizeof(float));
71 float* Xr
= (float*)malloc(n
*p
*sizeof(float));
72 float* Yr
= (float*)malloc(n
*m
*sizeof(float));
73 float* tXrXr
= (float*)malloc(p
*p
*sizeof(float));
74 float* tXrYr
= (float*)malloc(p
*m
*sizeof(float));
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
);
82 //Initialize class memberships (all elements in class 0; TODO: randomize ?)
83 int* Z
= (int*)calloc(n
, sizeof(int));
85 //Initialize phi to zero, because some M loops might exit before phi affectation
86 zeroArray(phi
, p
*m
*k
);
88 while (ite
<mini
|| (ite
<maxi
&& sumDeltaPhi
>tau
))
94 //M step: Mise à jour de Beta (et donc phi)
95 for (int r
=0; r
<k
; r
++)
97 //Compute Xr = X(Z==r,:) and Yr = Y(Z==r,:)
99 for (int i
=0; i
<n
; i
++)
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
)];
113 //Compute tXrXr = t(Xr) * Xr
114 for (int j
=0; j
<p
; j
++)
116 for (int jj
=0; jj
<p
; jj
++)
118 float 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
;
125 //Get pseudo inverse = (t(Xr)*Xr)^{-1}
126 float* invTXrXr
= pinv(tXrXr
, p
);
128 // Compute tXrYr = t(Xr) * Yr
129 for (int j
=0; j
<p
; j
++)
131 for (int jj
=0; jj
<m
; jj
++)
133 float dotProduct
= 0.0;
134 for (int u
=0; u
<cardClustR
; u
++)
135 dotProduct
+= Xr
[mi(u
,j
,n
,p
)] * Yr
[mi(u
,j
,n
,m
)];
136 tXrYr
[mi(j
,jj
,p
,m
)] = dotProduct
;
140 //Fill matrixM with inverse * tXrYr = (t(Xr)*Xr)^{-1} * t(Xr) * Yr
141 for (int j
=0; j
<p
; j
++)
143 for (int jj
=0; jj
<m
; jj
++)
145 float 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
;
153 //U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr
154 gsl_linalg_SV_decomp(matrixM
, V
, S
, work
);
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
++)
161 //[intermediate step] Compute hatBetaR = U * S * t(V)
162 double* U
= matrixM
->data
;
163 for (int j
=0; j
<p
; j
++)
165 for (int jj
=0; jj
<m
; jj
++)
167 float 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
;
174 //Compute phi(:,:,r) = hatBetaR * Rho(:,:,r)
175 for (int j
=0; j
<p
; j
++)
177 for (int jj
=0; jj
<m
; jj
++)
179 float 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
;
191 float sumLogLLF2
= 0.0;
192 for (int i
=0; i
<n
; i
++)
195 float maxLogGamIR
= -INFINITY
;
196 for (int r
=0; r
<k
; r
++)
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
203 //compute det(Rho(:,:,r)) [TODO: avoid re-computations]
204 for (int j
=0; j
<m
; j
++)
206 for (int jj
=0; jj
<m
; jj
++)
207 matrixE
->data
[j
*m
+jj
] = Rho
[ai(j
,jj
,r
,m
,m
,k
)];
209 gsl_linalg_LU_decomp(matrixE
, permutation
, &signum
);
210 float detRhoR
= gsl_linalg_LU_det(matrixE
, signum
);
212 //compute Y(i,:)*Rho(:,:,r)
213 for (int j
=0; j
<m
; j
++)
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
)];
220 //compute X(i,:)*phi(:,:,r)
221 for (int j
=0; j
<m
; j
++)
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
)];
228 //compute dotProduct < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) >
229 float dotProduct
= 0.0;
230 for (int u
=0; u
<m
; u
++)
231 dotProduct
+= (YiRhoR
[u
]-XiPhiR
[u
]) * (YiRhoR
[u
]-XiPhiR
[u
]);
232 float logGamIR
= log(Pi
[r
]) + log(detRhoR
) - 0.5*dotProduct
;
234 //Z(i) = index of max (gam(i,:))
235 if (logGamIR
> maxLogGamIR
)
238 maxLogGamIR
= logGamIR
;
240 sumLLF1
+= exp(logGamIR
) / pow(2*M_PI
,m
/2.0);
243 sumLogLLF2
+= log(sumLLF1
);
246 // Assign output variable LLF
247 *LLF
= -invN
* sumLogLLF2
;
249 //newDeltaPhi = max(max((abs(phi-Phi))./(1+abs(phi))));
250 float newDeltaPhi
= 0.0;
251 for (int j
=0; j
<p
; j
++)
253 for (int jj
=0; jj
<m
; jj
++)
255 for (int r
=0; r
<k
; r
++)
257 float 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
;
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
;
271 sumDeltaPhi
-= deltaPhi
[0];
272 for (int u
=0; u
<deltaPhiBufferSize
-1; u
++)
273 deltaPhi
[u
] = deltaPhi
[u
+1];
274 deltaPhi
[deltaPhiBufferSize
-1] = newDeltaPhi
;
276 sumDeltaPhi
+= newDeltaPhi
;
278 // update other local variables
279 for (int j
=0; j
<m
; j
++)
281 for (int jj
=0; jj
<p
; jj
++)
283 for (int r
=0; r
<k
; r
++)
284 Phi
[ai(j
,jj
,r
,p
,m
,k
)] = phi
[ai(j
,jj
,r
,p
,m
,k
)];
294 gsl_matrix_free(matrixE
);
295 gsl_matrix_free(matrixM
);
296 gsl_permutation_free(permutation
);
297 gsl_vector_free(work
);