SHOGUN
v3.2.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 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 /* give x an intialized value, from previously Theta*/ 00061 for(i=0;i<n;i++){ 00062 if (i==ith) 00063 continue; 00064 x[i]=Theta[i*n+ith]; 00065 } 00066 00067 /* Ax contains the derivative*/ 00068 m_Ax(Ax, W, x, n, ith); 00069 00070 for (iter_step=0;iter_step<maxIter; iter_step++){ 00071 00072 /*printf("\n Iter: %d",iter_step);*/ 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 the Lasso terminates either 00119 the change of the function value is less than fGap 00120 or the change of the solution in terms of L1 norm is less than xGap 00121 or the maximal iteration maxIter has achieved 00122 */ 00123 if ( (fabs(fun_new-fun_old) <=fGap) || x_change <=xGap){ 00124 /*printf("\n %d, Fun value: %2.5f",iter_step, fun_new); 00125 printf("\n The objective gap between adjacent solutions is less than %e",1e-6); 00126 */ 00127 break; 00128 } 00129 else{ 00130 /* 00131 if(iter_step%10 ==0) 00132 printf("\n %d, Fun value: %2.5f",iter_step, fun_new); 00133 */ 00134 fun_old=fun_new; 00135 } 00136 } 00137 00138 /*printf("\n %d, Fun value: %2.5f",iter_step, fun_new);*/ 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, /*for the Lasso (inner iteration)*/ 00169 int maxIter, double xtol) /*for the outer iteration*/ 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 /* printf("\n Outer Loop: %d, gap %e\n",iter_step,gap); */ 00209 00210 00211 if ( (gap <= xtol) || (iter_step==maxIter-1) ){ 00212 flag=1; 00213 } 00214 /* 00215 The outer loop terminates either the difference between ajacent solution in terms of L1 norm is less than xtol, 00216 or the maximal iterations has achieved 00217 */ 00218 } 00219 00220 free(W_old); 00221 00222 /*return (iter_step);*/ 00223 } 00224