pkg/src/EMGrank.c

   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 renormalise
  45   int mini, // nombre minimal d'iterations dans l'algorithme EM
  46   int maxi, // nombre maximal d'iterations dans l'algorithme EM
  47   const Real* X, // regresseurs
  48   const Real* Y, // reponse
  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 renormalise, calcule par l'EM
  53   Real* LLF, // log vraisemblance associe a cet echantillon, pour les valeurs estimees des parametres
  54   // additional size parameters
  55   int n, // taille de l'echantillon
  56   int p, // nombre de covariables
  57   int m, // taille de Y (multivarie)
  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 a 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       if (p >= m)
 155         gsl_linalg_SV_decomp(matrixM, V, S, work);
 156       else
 157       {
 158         gsl_matrix* matrixM_ = gsl_matrix_alloc(m, p);
 159         for (int j=0; j<m; j++)
 160         {
 161           for (int jj=0; jj<p; jj++)
 162             matrixM_->data[jj*m+j] = matrixM->data[j*p +jj];
 163         }
 164         gsl_matrix* V_ = gsl_matrix_alloc(p,p);
 165         gsl_vector* S_ = gsl_vector_alloc(p);
 166         gsl_vector* work_ = gsl_vector_alloc(p);
 167         gsl_linalg_SV_decomp(matrixM_, V_, S_, work_);
 168         gsl_vector_free(work_);
 169         for (int j=0; j<m; j++)
 170         {
 171           if (j < p)
 172             S->data[j] = S_->data[j];
 173           else
 174             S->data[j] = 0.0;
 175           for (int jj=0; jj<m; jj++)
 176           {
 177             if (j < p && jj < p)
 178               V->data[jj * m + j] = V_->data[jj * m + j];
 179             else
 180               V->data[jj * m + j] = 0.0;
 181           }
 182         }
 183         for (int j=0; j<m; j++)
 184         {
 185           for (int jj=0; jj<p; jj++)
 186             matrixM->data[j*p+jj] = matrixM_->data[jj*m +j];
 187         }
 188         gsl_matrix_free(matrixM_);
 189         gsl_matrix_free(V_);
 190         gsl_vector_free(S_);
 191       }
 192
 193       //Set m-rank(r) singular values to zero, and recompose
 194       //best rank(r) approximation of the initial product
 195       for (int j=rank[r]; j<m; j++)
 196         S->data[j] = 0.0;
 197
 198       //[intermediate step] Compute hatBetaR = U * S * t(V)
 199       double* U = matrixM->data; //GSL require double precision
 200       for (int j=0; j<p; j++)
 201       {
 202         for (int jj=0; jj<m; jj++)
 203         {
 204           Real dotProduct = 0.0;
 205           for (int u=0; u<m; u++)
 206             dotProduct += U[j*m+u] * S->data[u] * V->data[jj*m+u];
 207           hatBetaR[mi(j,jj,p,m)] = dotProduct;
 208         }
 209       }
 210
 211       //Compute phi(:,:,r) = hatBetaR * Rho(:,:,r)
 212       for (int j=0; j<p; j++)
 213       {
 214         for (int jj=0; jj<m; jj++)
 215         {
 216           Real dotProduct=0.0;
 217           for (int u=0; u<m; u++)
 218             dotProduct += hatBetaR[mi(j,u,p,m)] * Rho[ai(u,jj,r,m,m,k)];
 219           phi[ai(j,jj,r,p,m,k)] = dotProduct;
 220         }
 221       }
 222     }
 223
 224     /////////////
 225     // Etape E //
 226     /////////////
 227
 228     Real sumLogLLF2 = 0.0;
 229     for (int i=0; i<n; i++)
 230     {
 231       Real sumLLF1 = 0.0;
 232       Real maxLogGamIR = -INFINITY;
 233       for (int r=0; r<k; r++)
 234       {
 235         //Compute
 236         //Gam(i,r) = Pi(r) * det(Rho(:,:,r)) * exp( -1/2 * (Y(i,:)*Rho(:,:,r) - X(i,:)...
 237         //*phi(:,:,r)) * transpose( Y(i,:)*Rho(:,:,r) - X(i,:)*phi(:,:,r) ) );
 238         //split in several sub-steps
 239
 240         //compute det(Rho(:,:,r)) [TODO: avoid re-computations]
 241         for (int j=0; j<m; j++)
 242         {
 243           for (int jj=0; jj<m; jj++)
 244             matrixE->data[j*m+jj] = Rho[ai(j,jj,r,m,m,k)];
 245         }
 246         gsl_linalg_LU_decomp(matrixE, permutation, &signum);
 247         Real detRhoR = gsl_linalg_LU_det(matrixE, signum);
 248
 249         //compute Y(i,:)*Rho(:,:,r)
 250         for (int j=0; j<m; j++)
 251         {
 252           YiRhoR[j] = 0.0;
 253           for (int u=0; u<m; u++)
 254             YiRhoR[j] += Y[mi(i,u,n,m)] * Rho[ai(u,j,r,m,m,k)];
 255         }
 256
 257         //compute X(i,:)*phi(:,:,r)
 258         for (int j=0; j<m; j++)
 259         {
 260           XiPhiR[j] = 0.0;
 261           for (int u=0; u<p; u++)
 262             XiPhiR[j] += X[mi(i,u,n,p)] * phi[ai(u,j,r,p,m,k)];
 263         }
 264
 265         //compute dotProduct < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) >
 266         Real dotProduct = 0.0;
 267         for (int u=0; u<m; u++)
 268           dotProduct += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]);
 269         Real logGamIR = log(Pi[r]) + log(detRhoR) - 0.5*dotProduct;
 270
 271         //Z(i) = index of max (gam(i,:))
 272         if (logGamIR > maxLogGamIR)
 273         {
 274           Z[i] = r;
 275           maxLogGamIR = logGamIR;
 276         }
 277         sumLLF1 += exp(logGamIR) / pow(2*M_PI,m/2.0);
 278       }
 279
 280       sumLogLLF2 += log(sumLLF1);
 281     }
 282
 283     // Assign output variable LLF
 284     *LLF = -invN * sumLogLLF2;
 285
 286     //newDeltaPhi = max(max((abs(phi-Phi))./(1+abs(phi))));
 287     Real newDeltaPhi = 0.0;
 288     for (int j=0; j<p; j++)
 289     {
 290       for (int jj=0; jj<m; jj++)
 291       {
 292         for (int r=0; r<k; r++)
 293         {
 294           Real tmpDist = fabs(phi[ai(j,jj,r,p,m,k)]-Phi[ai(j,jj,r,p,m,k)])
 295             / (1.0+fabs(phi[ai(j,jj,r,p,m,k)]));
 296           if (tmpDist > newDeltaPhi)
 297             newDeltaPhi = tmpDist;
 298         }
 299       }
 300     }
 301
 302     //update distance parameter to check algorithm convergence (delta(phi, Phi))
 303     //TODO: deltaPhi should be a linked list for perf.
 304     if (ite < deltaPhiBufferSize)
 305       deltaPhi[ite] = newDeltaPhi;
 306     else
 307     {
 308       sumDeltaPhi -= deltaPhi[0];
 309       for (int u=0; u<deltaPhiBufferSize-1; u++)
 310         deltaPhi[u] = deltaPhi[u+1];
 311       deltaPhi[deltaPhiBufferSize-1] = newDeltaPhi;
 312     }
 313     sumDeltaPhi += newDeltaPhi;
 314
 315     // update other local variables
 316     for (int j=0; j<m; j++)
 317     {
 318       for (int jj=0; jj<p; jj++)
 319       {
 320         for (int r=0; r<k; r++)
 321           Phi[ai(jj,j,r,p,m,k)] = phi[ai(jj,j,r,p,m,k)];
 322       }
 323     }
 324     ite++;
 325   }
 326
 327   //free memory
 328   free(hatBetaR);
 329   free(deltaPhi);
 330   free(Phi);
 331   gsl_matrix_free(matrixE);
 332   gsl_matrix_free(matrixM);
 333   gsl_permutation_free(permutation);
 334   gsl_vector_free(work);
 335   gsl_matrix_free(V);
 336   gsl_vector_free(S);
 337   free(XiPhiR);
 338   free(YiRhoR);
 339   free(Xr);
 340   free(Yr);
 341   free(tXrXr);
 342   free(tXrYr);
 343   free(Z);
 344 }