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

SHOGUN Machine Learning Toolbox - Documentation