SHOGUN  4.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
19 
20 #include <stdlib.h>
21 #include <time.h>
22 
23 void m_Ax(double *Ax, double *A, double *x, int n, int ith)
24 {
25  int i, j;
26  double t;
27  for(i=0;i<n;i++){
28  if (i==ith)
29  continue;
30  t=0;
31  for(j=0;j<n;j++){
32  if (j==ith)
33  continue;
34  t+=A[i*n+j]* x[j];
35  Ax[i]=t;
36  }
37  }
38 }
39 
40 
41 int lassoCD(double *Theta, double *W, double *S, double lambda, int n, int ith, int flag, int maxIter, double fGap, double xGap)
42 {
43  int iter_step, i,j;
44  double * Ax, * x;
45  double u, v, s_v, t=0, x_new;
46  double fun_new,fun_old=-100;
47  double x_change;
48 
49  Ax= (double *)malloc(sizeof(double)*n);
50  if (Ax==NULL)
51  {
52  printf("\n Memory allocation failure!");
53  return (-1);
54  }
55 
56  x= (double *)malloc(sizeof(double)*n);
57  if (x==NULL)
58  {
59  printf("\n Memory allocation failure!");
60  free(Ax);
61  return (-1);
62  }
63 
64  /* give x an intialized value, from previously Theta*/
65  for(i=0;i<n;i++){
66  if (i==ith)
67  continue;
68  x[i]=Theta[i*n+ith];
69  }
70 
71  /* Ax contains the derivative*/
72  m_Ax(Ax, W, x, n, ith);
73 
74  for (iter_step=0;iter_step<maxIter; iter_step++){
75 
76  /*printf("\n Iter: %d",iter_step);*/
77 
78  x_change=0;
79 
80  for (i=0;i<n;i++){
81  if(i==ith)
82  continue;
83 
84  u=W[i*n + i];
85 
86  v=Ax[i]-x[i]*u;
87 
88  s_v=S[i*n+ ith]-v;
89 
90  if(s_v > lambda)
91  x_new= (s_v-lambda) / u;
92  else{
93  if(s_v < -lambda)
94  x_new= (s_v + lambda) / u;
95  else
96  x_new=0;
97  }
98  if (x[i]!=x_new){
99  for(j=0;j<n;j++){
100  if (j==ith)
101  continue;
102  Ax[j]+= W[j*n+ i]*(x_new - x[i]);
103  }
104  x_change+=fabs(x[i]-x_new);
105 
106  x[i]=x_new;
107  }
108  }
109 
110  fun_new=0;
111  t=0;
112  for(i=0;i<n;i++){
113  if (i==ith)
114  continue;
115  t+= Ax[i]*x[i] ;
116  fun_new+=- S[i*n+ith]*x[i]+ lambda* fabs(x[i]);
117  }
118  fun_new+=0.5*t;
119 
120 
121  /*
122  the Lasso terminates either
123  the change of the function value is less than fGap
124  or the change of the solution in terms of L1 norm is less than xGap
125  or the maximal iteration maxIter has achieved
126  */
127  if ( (fabs(fun_new-fun_old) <=fGap) || x_change <=xGap){
128  /*printf("\n %d, Fun value: %2.5f",iter_step, fun_new);
129  printf("\n The objective gap between adjacent solutions is less than %e",1e-6);
130  */
131  break;
132  }
133  else{
134  /*
135  if(iter_step%10 ==0)
136  printf("\n %d, Fun value: %2.5f",iter_step, fun_new);
137  */
138  fun_old=fun_new;
139  }
140  }
141 
142  /*printf("\n %d, Fun value: %2.5f",iter_step, fun_new);*/
143 
144  if (flag){
145  t=1/(W[ith*n+ith]-t);
146  Theta[ith*n + ith]=t;
147 
148  for(i=0;i<n;i++){
149  if (i==ith)
150  continue;
151  W[i*n+ ith]=W[ith*n +i]=Ax[i];
152  Theta[i*n+ ith]=Theta[ith*n +i]=-x[i]*t;
153  }
154  }
155  else{
156  for(i=0;i<n;i++){
157  if (i==ith)
158  continue;
159  W[i*n+ ith]=W[ith*n +i]=Ax[i];
160  Theta[i*n+ ith]=Theta[ith*n +i]=x[i];
161  }
162  }
163 
164 
165  free(Ax); free(x);
166 
167  return(iter_step);
168 }
169 
170 
171 void invCov(double *Theta, double *W, double *S, double lambda, double sum_S, int n,
172  int LassoMaxIter, double fGap, double xGap, /*for the Lasso (inner iteration)*/
173  int maxIter, double xtol) /*for the outer iteration*/
174 {
175  int iter_step, i,j, ith;
176  double * W_old;
177  double gap;
178  int flag=0;
179 
180  W_old= (double *)malloc(sizeof(double)*n*n);
181 
182 
183  if ( W_old==NULL ){
184  printf("\n Memory allocation failure!");
185  exit (-1);
186  }
187 
188  for(i=0;i<n;i++)
189  for(j=0;j<n;j++){
190  if (i==j)
191  W_old[i*n+j]=W[i*n+j]=S[i*n+j]+lambda;
192  else
193  W_old[i*n+j]=W[i*n+j]=S[i*n+j];
194 
195  Theta[i*n+j]=0;
196  }
197 
198  for (iter_step=0;iter_step<=maxIter; iter_step++){
199  for(ith=0;ith<n;ith++)
200  lassoCD(Theta, W, S, lambda, n, ith, flag, LassoMaxIter,fGap, xGap);
201 
202  if (flag)
203  break;
204 
205  gap=0;
206  for(i=0;i<n;i++)
207  for(j=0;j<n;j++){
208  gap+=fabs(W[i*n+j]-W_old[i*n+j]);
209  W_old[i*n+j]=W[i*n+j];
210  }
211 
212  /* printf("\n Outer Loop: %d, gap %e\n",iter_step,gap); */
213 
214 
215  if ( (gap <= xtol) || (iter_step==maxIter-1) ){
216  flag=1;
217  }
218  /*
219  The outer loop terminates either the difference between ajacent solution in terms of L1 norm is less than xtol,
220  or the maximal iterations has achieved
221  */
222  }
223 
224  free(W_old);
225 
226  /*return (iter_step);*/
227 }
228 

SHOGUN Machine Learning Toolbox - Documentation