SHOGUN  v2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
invCov.cpp
Go to the documentation of this file.
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation