SHOGUN  5.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
invCov.cpp
Go to the documentation of this file.
1 /* This program is free software: you can redistribute it and/or modify
2  * it under the terms of the GNU General Public License as published by
3  * the Free Software Foundation, either version 3 of the License, or
4  * (at your option) any later version.
5  *
6  * This program is distributed in the hope that it will be useful,
7  * but WITHOUT ANY WARRANTY; without even the implied warranty of
8  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
9  * GNU General Public License for more details.
10  *
11  * You should have received a copy of the GNU General Public License
12  * along with this program. If not, see <http://www.gnu.org/licenses/>.
13  *
14  * Copyright (C) 2009 - 2012 Jun Liu and Jieping Ye
15  */
16 
17 
19 #ifdef USE_GPL_SHOGUN
21 
22 #include <stdlib.h>
23 #include <time.h>
24 
25 void m_Ax(double *Ax, double *A, double *x, int n, int ith)
26 {
27  int i, j;
28  double t;
29  for(i=0;i<n;i++){
30  if (i==ith)
31  continue;
32  t=0;
33  for(j=0;j<n;j++){
34  if (j==ith)
35  continue;
36  t+=A[i*n+j]* x[j];
37  Ax[i]=t;
38  }
39  }
40 }
41 
42 
43 int lassoCD(double *Theta, double *W, double *S, double lambda, int n, int ith, int flag, int maxIter, double fGap, double xGap)
44 {
45  int iter_step, i,j;
46  double * Ax, * x;
47  double u, v, s_v, t=0, x_new;
48  double fun_new,fun_old=-100;
49  double x_change;
50 
51  Ax= (double *)malloc(sizeof(double)*n);
52  if (Ax==NULL)
53  {
54  printf("\n Memory allocation failure!");
55  return (-1);
56  }
57 
58  x= (double *)malloc(sizeof(double)*n);
59  if (x==NULL)
60  {
61  printf("\n Memory allocation failure!");
62  free(Ax);
63  return (-1);
64  }
65 
66  /* give x an intialized value, from previously Theta*/
67  for(i=0;i<n;i++){
68  if (i==ith)
69  continue;
70  x[i]=Theta[i*n+ith];
71  }
72 
73  /* Ax contains the derivative*/
74  m_Ax(Ax, W, x, n, ith);
75 
76  for (iter_step=0;iter_step<maxIter; iter_step++){
77 
78  /*printf("\n Iter: %d",iter_step);*/
79 
80  x_change=0;
81 
82  for (i=0;i<n;i++){
83  if(i==ith)
84  continue;
85 
86  u=W[i*n + i];
87 
88  v=Ax[i]-x[i]*u;
89 
90  s_v=S[i*n+ ith]-v;
91 
92  if(s_v > lambda)
93  x_new= (s_v-lambda) / u;
94  else{
95  if(s_v < -lambda)
96  x_new= (s_v + lambda) / u;
97  else
98  x_new=0;
99  }
100  if (x[i]!=x_new){
101  for(j=0;j<n;j++){
102  if (j==ith)
103  continue;
104  Ax[j]+= W[j*n+ i]*(x_new - x[i]);
105  }
106  x_change+=fabs(x[i]-x_new);
107 
108  x[i]=x_new;
109  }
110  }
111 
112  fun_new=0;
113  t=0;
114  for(i=0;i<n;i++){
115  if (i==ith)
116  continue;
117  t+= Ax[i]*x[i] ;
118  fun_new+=- S[i*n+ith]*x[i]+ lambda* fabs(x[i]);
119  }
120  fun_new+=0.5*t;
121 
122 
123  /*
124  the Lasso terminates either
125  the change of the function value is less than fGap
126  or the change of the solution in terms of L1 norm is less than xGap
127  or the maximal iteration maxIter has achieved
128  */
129  if ( (fabs(fun_new-fun_old) <=fGap) || x_change <=xGap){
130  /*printf("\n %d, Fun value: %2.5f",iter_step, fun_new);
131  printf("\n The objective gap between adjacent solutions is less than %e",1e-6);
132  */
133  break;
134  }
135  else{
136  /*
137  if(iter_step%10 ==0)
138  printf("\n %d, Fun value: %2.5f",iter_step, fun_new);
139  */
140  fun_old=fun_new;
141  }
142  }
143 
144  /*printf("\n %d, Fun value: %2.5f",iter_step, fun_new);*/
145 
146  if (flag){
147  t=1/(W[ith*n+ith]-t);
148  Theta[ith*n + ith]=t;
149 
150  for(i=0;i<n;i++){
151  if (i==ith)
152  continue;
153  W[i*n+ ith]=W[ith*n +i]=Ax[i];
154  Theta[i*n+ ith]=Theta[ith*n +i]=-x[i]*t;
155  }
156  }
157  else{
158  for(i=0;i<n;i++){
159  if (i==ith)
160  continue;
161  W[i*n+ ith]=W[ith*n +i]=Ax[i];
162  Theta[i*n+ ith]=Theta[ith*n +i]=x[i];
163  }
164  }
165 
166 
167  free(Ax); free(x);
168 
169  return(iter_step);
170 }
171 
172 
173 void invCov(double *Theta, double *W, double *S, double lambda, double sum_S, int n,
174  int LassoMaxIter, double fGap, double xGap, /*for the Lasso (inner iteration)*/
175  int maxIter, double xtol) /*for the outer iteration*/
176 {
177  int iter_step, i,j, ith;
178  double * W_old;
179  double gap;
180  int flag=0;
181 
182  W_old= (double *)malloc(sizeof(double)*n*n);
183 
184 
185  if ( W_old==NULL ){
186  printf("\n Memory allocation failure!");
187  exit (-1);
188  }
189 
190  for(i=0;i<n;i++)
191  for(j=0;j<n;j++){
192  if (i==j)
193  W_old[i*n+j]=W[i*n+j]=S[i*n+j]+lambda;
194  else
195  W_old[i*n+j]=W[i*n+j]=S[i*n+j];
196 
197  Theta[i*n+j]=0;
198  }
199 
200  for (iter_step=0;iter_step<=maxIter; iter_step++){
201  for(ith=0;ith<n;ith++)
202  lassoCD(Theta, W, S, lambda, n, ith, flag, LassoMaxIter,fGap, xGap);
203 
204  if (flag)
205  break;
206 
207  gap=0;
208  for(i=0;i<n;i++)
209  for(j=0;j<n;j++){
210  gap+=fabs(W[i*n+j]-W_old[i*n+j]);
211  W_old[i*n+j]=W[i*n+j];
212  }
213 
214  /* printf("\n Outer Loop: %d, gap %e\n",iter_step,gap); */
215 
216 
217  if ( (gap <= xtol) || (iter_step==maxIter-1) ){
218  flag=1;
219  }
220  /*
221  The outer loop terminates either the difference between ajacent solution in terms of L1 norm is less than xtol,
222  or the maximal iterations has achieved
223  */
224  }
225 
226  free(W_old);
227 
228  /*return (iter_step);*/
229 }
230 
231 #endif //USE_GPL_SHOGUN

SHOGUN Machine Learning Toolbox - Documentation