00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <shogun/lib/slep/SpInvCoVa/invCov.h>
00018
00019 void m_Ax(double *Ax, double *A, double *x, int n, int ith)
00020 {
00021 int i, j;
00022 double t;
00023 for(i=0;i<n;i++){
00024 if (i==ith)
00025 continue;
00026 t=0;
00027 for(j=0;j<n;j++){
00028 if (j==ith)
00029 continue;
00030 t+=A[i*n+j]* x[j];
00031 Ax[i]=t;
00032 }
00033 }
00034 }
00035
00036
00037 int lassoCD(double *Theta, double *W, double *S, double lambda, int n, int ith, int flag, int maxIter, double fGap, double xGap)
00038 {
00039 int iter_step, i,j;
00040 double * Ax, * x;
00041 double u, v, s_v, t=0, x_new;
00042 double fun_new,fun_old=-100;
00043 double x_change;
00044
00045 Ax= (double *)malloc(sizeof(double)*n);
00046 if (Ax==NULL)
00047 {
00048 printf("\n Memory allocation failure!");
00049 return (-1);
00050 }
00051
00052 x= (double *)malloc(sizeof(double)*n);
00053 if (x==NULL)
00054 {
00055 printf("\n Memory allocation failure!");
00056 free(Ax);
00057 return (-1);
00058 }
00059
00060
00061 for(i=0;i<n;i++){
00062 if (i==ith)
00063 continue;
00064 x[i]=Theta[i*n+ith];
00065 }
00066
00067
00068 m_Ax(Ax, W, x, n, ith);
00069
00070 for (iter_step=0;iter_step<maxIter; iter_step++){
00071
00072
00073
00074 x_change=0;
00075
00076 for (i=0;i<n;i++){
00077 if(i==ith)
00078 continue;
00079
00080 u=W[i*n + i];
00081
00082 v=Ax[i]-x[i]*u;
00083
00084 s_v=S[i*n+ ith]-v;
00085
00086 if(s_v > lambda)
00087 x_new= (s_v-lambda) / u;
00088 else{
00089 if(s_v < -lambda)
00090 x_new= (s_v + lambda) / u;
00091 else
00092 x_new=0;
00093 }
00094 if (x[i]!=x_new){
00095 for(j=0;j<n;j++){
00096 if (j==ith)
00097 continue;
00098 Ax[j]+= W[j*n+ i]*(x_new - x[i]);
00099 }
00100 x_change+=fabs(x[i]-x_new);
00101
00102 x[i]=x_new;
00103 }
00104 }
00105
00106 fun_new=0;
00107 t=0;
00108 for(i=0;i<n;i++){
00109 if (i==ith)
00110 continue;
00111 t+= Ax[i]*x[i] ;
00112 fun_new+=- S[i*n+ith]*x[i]+ lambda* fabs(x[i]);
00113 }
00114 fun_new+=0.5*t;
00115
00116
00117
00118
00119
00120
00121
00122
00123 if ( (fabs(fun_new-fun_old) <=fGap) || x_change <=xGap){
00124
00125
00126
00127 break;
00128 }
00129 else{
00130
00131
00132
00133
00134 fun_old=fun_new;
00135 }
00136 }
00137
00138
00139
00140 if (flag){
00141 t=1/(W[ith*n+ith]-t);
00142 Theta[ith*n + ith]=t;
00143
00144 for(i=0;i<n;i++){
00145 if (i==ith)
00146 continue;
00147 W[i*n+ ith]=W[ith*n +i]=Ax[i];
00148 Theta[i*n+ ith]=Theta[ith*n +i]=-x[i]*t;
00149 }
00150 }
00151 else{
00152 for(i=0;i<n;i++){
00153 if (i==ith)
00154 continue;
00155 W[i*n+ ith]=W[ith*n +i]=Ax[i];
00156 Theta[i*n+ ith]=Theta[ith*n +i]=x[i];
00157 }
00158 }
00159
00160
00161 free(Ax); free(x);
00162
00163 return(iter_step);
00164 }
00165
00166
00167 void invCov(double *Theta, double *W, double *S, double lambda, double sum_S, int n,
00168 int LassoMaxIter, double fGap, double xGap,
00169 int maxIter, double xtol)
00170 {
00171 int iter_step, i,j, ith;
00172 double * W_old;
00173 double gap;
00174 int flag=0;
00175
00176 W_old= (double *)malloc(sizeof(double)*n*n);
00177
00178
00179 if ( W_old==NULL ){
00180 printf("\n Memory allocation failure!");
00181 exit (-1);
00182 }
00183
00184 for(i=0;i<n;i++)
00185 for(j=0;j<n;j++){
00186 if (i==j)
00187 W_old[i*n+j]=W[i*n+j]=S[i*n+j]+lambda;
00188 else
00189 W_old[i*n+j]=W[i*n+j]=S[i*n+j];
00190
00191 Theta[i*n+j]=0;
00192 }
00193
00194 for (iter_step=0;iter_step<=maxIter; iter_step++){
00195 for(ith=0;ith<n;ith++)
00196 lassoCD(Theta, W, S, lambda, n, ith, flag, LassoMaxIter,fGap, xGap);
00197
00198 if (flag)
00199 break;
00200
00201 gap=0;
00202 for(i=0;i<n;i++)
00203 for(j=0;j<n;j++){
00204 gap+=fabs(W[i*n+j]-W_old[i*n+j]);
00205 W_old[i*n+j]=W[i*n+j];
00206 }
00207
00208
00209
00210
00211 if ( (gap <= xtol) || (iter_step==maxIter-1) ){
00212 flag=1;
00213 }
00214
00215
00216
00217
00218 }
00219
00220 free(W_old);
00221
00222
00223 }
00224