SHOGUN
v2.0.0
|
00001 /* This program is free software: you can redistribute it and/or modify 00002 * it under the terms of the GNU General Public License as published by 00003 * the Free Software Foundation, either version 3 of the License, or 00004 * (at your option) any later version. 00005 * 00006 * This program is distributed in the hope that it will be useful, 00007 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00008 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00009 * GNU General Public License for more details. 00010 * 00011 * You should have received a copy of the GNU General Public License 00012 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00013 * 00014 * Copyright (C) 2009 - 2012 Jun Liu and Jieping Ye 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 x= (double *)malloc(sizeof(double)*n); 00047 00048 00049 if ( (Ax==NULL) || (x==NULL) ){ 00050 printf("\n Memory allocation failure!"); 00051 return (-1); 00052 } 00053 00054 /* give x an intialized value, from previously Theta*/ 00055 for(i=0;i<n;i++){ 00056 if (i==ith) 00057 continue; 00058 x[i]=Theta[i*n+ith]; 00059 } 00060 00061 /* Ax contains the derivative*/ 00062 m_Ax(Ax, W, x, n, ith); 00063 00064 for (iter_step=0;iter_step<maxIter; iter_step++){ 00065 00066 /*printf("\n Iter: %d",iter_step);*/ 00067 00068 x_change=0; 00069 00070 for (i=0;i<n;i++){ 00071 if(i==ith) 00072 continue; 00073 00074 u=W[i*n + i]; 00075 00076 v=Ax[i]-x[i]*u; 00077 00078 s_v=S[i*n+ ith]-v; 00079 00080 if(s_v > lambda) 00081 x_new= (s_v-lambda) / u; 00082 else{ 00083 if(s_v < -lambda) 00084 x_new= (s_v + lambda) / u; 00085 else 00086 x_new=0; 00087 } 00088 if (x[i]!=x_new){ 00089 for(j=0;j<n;j++){ 00090 if (j==ith) 00091 continue; 00092 Ax[j]+= W[j*n+ i]*(x_new - x[i]); 00093 } 00094 x_change+=fabs(x[i]-x_new); 00095 00096 x[i]=x_new; 00097 } 00098 } 00099 00100 fun_new=0; 00101 t=0; 00102 for(i=0;i<n;i++){ 00103 if (i==ith) 00104 continue; 00105 t+= Ax[i]*x[i] ; 00106 fun_new+=- S[i*n+ith]*x[i]+ lambda* fabs(x[i]); 00107 } 00108 fun_new+=0.5*t; 00109 00110 00111 /* 00112 the Lasso terminates either 00113 the change of the function value is less than fGap 00114 or the change of the solution in terms of L1 norm is less than xGap 00115 or the maximal iteration maxIter has achieved 00116 */ 00117 if ( (fabs(fun_new-fun_old) <=fGap) || x_change <=xGap){ 00118 /*printf("\n %d, Fun value: %2.5f",iter_step, fun_new); 00119 printf("\n The objective gap between adjacent solutions is less than %e",1e-6); 00120 */ 00121 break; 00122 } 00123 else{ 00124 /* 00125 if(iter_step%10 ==0) 00126 printf("\n %d, Fun value: %2.5f",iter_step, fun_new); 00127 */ 00128 fun_old=fun_new; 00129 } 00130 } 00131 00132 /*printf("\n %d, Fun value: %2.5f",iter_step, fun_new);*/ 00133 00134 if (flag){ 00135 t=1/(W[ith*n+ith]-t); 00136 Theta[ith*n + ith]=t; 00137 00138 for(i=0;i<n;i++){ 00139 if (i==ith) 00140 continue; 00141 W[i*n+ ith]=W[ith*n +i]=Ax[i]; 00142 Theta[i*n+ ith]=Theta[ith*n +i]=-x[i]*t; 00143 } 00144 } 00145 else{ 00146 for(i=0;i<n;i++){ 00147 if (i==ith) 00148 continue; 00149 W[i*n+ ith]=W[ith*n +i]=Ax[i]; 00150 Theta[i*n+ ith]=Theta[ith*n +i]=x[i]; 00151 } 00152 } 00153 00154 00155 free(Ax); free(x); 00156 00157 return(iter_step); 00158 } 00159 00160 00161 void invCov(double *Theta, double *W, double *S, double lambda, double sum_S, int n, 00162 int LassoMaxIter, double fGap, double xGap, /*for the Lasso (inner iteration)*/ 00163 int maxIter, double xtol) /*for the outer iteration*/ 00164 { 00165 int iter_step, i,j, ith; 00166 double * W_old; 00167 double gap; 00168 int flag=0; 00169 00170 W_old= (double *)malloc(sizeof(double)*n*n); 00171 00172 00173 if ( W_old==NULL ){ 00174 printf("\n Memory allocation failure!"); 00175 exit (-1); 00176 } 00177 00178 for(i=0;i<n;i++) 00179 for(j=0;j<n;j++){ 00180 if (i==j) 00181 W_old[i*n+j]=W[i*n+j]=S[i*n+j]+lambda; 00182 else 00183 W_old[i*n+j]=W[i*n+j]=S[i*n+j]; 00184 00185 Theta[i*n+j]=0; 00186 } 00187 00188 for (iter_step=0;iter_step<=maxIter; iter_step++){ 00189 for(ith=0;ith<n;ith++) 00190 lassoCD(Theta, W, S, lambda, n, ith, flag, LassoMaxIter,fGap, xGap); 00191 00192 if (flag) 00193 break; 00194 00195 gap=0; 00196 for(i=0;i<n;i++) 00197 for(j=0;j<n;j++){ 00198 gap+=fabs(W[i*n+j]-W_old[i*n+j]); 00199 W_old[i*n+j]=W[i*n+j]; 00200 } 00201 00202 /* printf("\n Outer Loop: %d, gap %e\n",iter_step,gap); */ 00203 00204 00205 if ( (gap <= xtol) || (iter_step==maxIter-1) ){ 00206 flag=1; 00207 } 00208 /* 00209 The outer loop terminates either the difference between ajacent solution in terms of L1 norm is less than xtol, 00210 or the maximal iterations has achieved 00211 */ 00212 } 00213 00214 free(W_old); 00215 00216 /*return (iter_step);*/ 00217 } 00218