# source:trunk/sinco/sinco.cpp@10

Last change on this file since 10 was 10, checked in by katyas, 5 years ago

Rename the sparse inverse covariance to SINCO and made a first real version available

File size: 13.0 KB
Line
1#include <stdlib.h>
2#include <vector>
3#include <cstdio>
4#include <stdio.h>
5#include <sstream>
6#include <iostream>
7#include <fstream>
8#include <cmath>
9
10#include "sinco.hpp"
11
12/* This is the C++ implementation to solve the sparse inverse covariance
13   problem of the form
14
15max _C[ K*log(det(C)-tr(AC)-\lambda ||S.*C||_1]
16
17K and lambda are positive scalars, A is a given symmetric matrix, S is a given nonnegative matrix
18(not necessarily symmetric), S.*C is a elemwhise product and ||*||_1 is the sum
19of all absolute values.
20
21   */
22
23double funvalue_update(double alpha,  double K, double Wij, double  Wii, double Wjj,
24                     double Aij, double Sij,double  Sji, double update) {
25  double detratio, fchange;
26  /* This routine computes the change in the objective function value when a step of length\ alpha is made in the direction e_ie_j^T+e_je_i^T  */
27if (update ==1) {
28  detratio=(1+alpha*Wij)*(1+alpha*Wij-alpha*alpha*Wii*Wjj/(1+alpha*Wij));
29  fchange=K*(log(detratio))-2*alpha*Aij-alpha*(Sij+Sji);
30 }
31 else {
32   detratio=(1-alpha*Wij)*(1-alpha*Wij-alpha*alpha*Wii*Wjj/(1-alpha*Wij));
33  fchange=K*(log(detratio))+2*alpha*Aij-alpha*(Sij+Sji);
34 }
35 return fchange;
36}
37
38
39void invupdate(Matrix& Gradp,  Matrix& Gradpp,  vector<int>& UpInd, Matrix& W, const double theta, const double K, const int p, int i, int j, int update)
40{
41 /* This routine computes the change in the dual matrix W (inverse of C) when a step of length\ alpha is made in the direction e_ie_j^T+e_je_i^T  */
42double a,b,c,d, wij, wjj, wii;
43 int ii, jj, upindind;
44Matrix UpW(p,p);
45wii=W(i,i);
46wjj=W(j,j);
47wij=W(i,j);
48if (update==1){
49 d=theta*theta*(wii*wjj-wij*wij)-1-2*theta*wij;
50 a=-(1+theta*wij)/d;
51 b=theta*wjj/d;
52 c=theta*wii/d;
53 for (int ii=p-1; ii>=0; ii--) {
54   for (int jj=ii; jj>=0; jj--){
55     upindind=ii*p+jj;
56     UpW(ii,jj)=-theta*(a*W(i,ii)*W(j,jj)+ c*W(j,ii)*W(j,jj)+b*W(i,ii)*W(i,jj)+a*W(j,ii)*W(i,jj));
57     if (fabs(UpW(ii,jj))>zerotol){UpInd[upindind]=1;}
58     else {UpInd[ii*p+jj]=0;}
59     UpW(jj,ii)=UpW(ii,jj);
60     UpInd[jj*p+ii]=UpInd[upindind];
61   }
62 }
63}
64 else{
65 d=theta*theta*(-wii*wjj+wij*wij)+1-2*theta*wij;
66 a=(1-theta*wij)/d;
67 b=theta*wjj/d;
68 c=theta*wii/d;
69
70 for (int ii=p-1; ii>=0; ii--) {
71   for (int jj=ii; jj>=0; jj--) {
72     upindind=ii*p+jj;
73     UpW(ii,jj)=theta*(a*W(i,ii)*W(j,jj)+ c*W(j,ii)*W(j,jj)+b*W(i,ii)*W(i,jj)+a*W(j,ii)*W(i,jj));
74     if (fabs(UpW(ii,jj))>zerotol){UpInd[upindind]=1;}
75     else {UpInd[ii*p+jj]=0;}
76     UpW(jj,ii)=UpW(ii,jj);
77     UpInd[jj*p+ii]=UpInd[upindind];
78   }
79 }
80}
81W.sum(UpW);
84
85}
86
87
88
89double findposstep(double K, double Wij, double Wii, double Wjj,
90                   double Aij, double Sij, double Sji, int update)  {
91 /* This routine computes the optimal length of  a step of length in the direction e_ie_j^T+e_je_i^T for the positive component of C */
92  double aux1, aux2, aux3, aux4, aux5;
93  double a, b, c, D, alpha, alpha1, alpha2;
94
95  if ( update ==1 ) {
96    aux1=2*(K*Wij-Aij)-Sji-Sij;
97    aux2=(Wii*Wjj-Wij*Wij);
98    aux3=2*Wij;
99    aux4=-2*K*(Wii*Wjj+Wij*Wij);
100    aux5=2*K*Wij*aux2;
101  }
102  else {
103    aux1=2*(-K*Wij+Aij)-Sij-Sji;
104    aux2=(-Wii*Wjj+Wij*Wij);
105    aux3=2*Wij;
106    aux4=2*K*(Wii*Wjj+Wij*Wij);
107    aux5=-2*K*Wij*aux2;
108  }
109  a=aux2*aux1-aux5;
110  b=-aux3*aux1-aux4;
111  c=-update*aux1;
112  if (fabs(a)>zerotol) {
113    D=b*b-4*a*c;
114    if (D<0) {
115      printf ( "negative discriminant, \n");
116      abort;
117    }
118    alpha1=fmin((-b-sqrt(D))/(2*a),(-b+sqrt(D))/(2*a));
119    alpha2=fmax((-b-sqrt(D))/(2*a),(-b+sqrt(D))/(2*a));
120    if (alpha1>=0){ alpha=alpha1;}
121    else  {
122      if (alpha2>=0){ alpha=alpha2;}
123       else {
124        printf ( "unbounded direction!!!, \n");
125        abort;
126      }
127    }
128  }
129  else if (-c/b>0)
130    { alpha=-c/b;}
131  else {
132    printf ( "unbounded direction!!!,\n");
133    abort;
134  }
135
136  /* Check correctness
137double theta=alpha;
138 double fprime;
139 if ( update ==1 ) {
140    fprime=(K*Wij+K*Wij-Aij-Aij-Sij-Sji)-
141     K*theta/(theta*theta*(Wii*Wjj-Wij*Wij)-1-2*theta*Wij)*(-2*(Wii*Wjj+Wij*Wij)
142   +2*theta*Wij*(Wii*Wjj-Wij*Wij));
143 }
144 else {
145   fprime=(-K*Wij-K*Wij+Aij+Aij-Sij-Sji)-K*theta/(theta*theta*(-Wii*Wjj+Wij*Wij)
146           +1-2*theta*Wij)*(2*(Wii*Wjj+Wij*Wij)
147           +2*theta*Wij*(Wii*Wjj-Wij*Wij));
148 }
149 if (fabs(fprime) > 10e-6){
151   fflush(stdout);
152    abort;
153    } */
154  return alpha;
155}
156
157
158double findnegstep(double K, bool diag, double Wij, double Wii, double Wjj,
159                   double Aij, double Sij, double Sji, double Cpij, double Cppij,  int update)  {
160/* This routine computes the optimal length of  a step of length in the direction e_ie_j^T+e_je_i^T for the negative component of C */
161  double aux1, aux2, aux3, aux4, aux5;
162  double a, b, c, D, maxstep, alpha, alpha1, alpha2;
163
164  if (update ==1) {
165   aux1=2*(K*Wij-Aij)-Sji-Sij;
166   aux2=(Wii*Wjj-Wij*Wij);
167   aux3=2*Wij;
168   aux4=-2*K*(Wii*Wjj+Wij*Wij);
169   aux5=2*K*Wij*aux2;
170   maxstep=Cpij;}
171  else {
172   aux1=2*(-K*Wij+Aij)-Sij-Sji;
173   aux2=(-Wii*Wjj+Wij*Wij);
174   aux3=2*Wij;
175   aux4=2*K*(Wii*Wjj+Wij*Wij);
176   aux5=-2*K*Wij*aux2;
177   maxstep=Cppij;
178}
179if (diag) {
180 maxstep=maxstep/2;
181}
182a=aux2*aux1-aux5;
183b=-aux3*aux1-aux4;
184c=-update*aux1;
185 if (fabs(a)>zerotol) {
186  D=b*b-4*a*c;
187  if (D<0) {
188      printf ( "negative discriminant, \n");
189      abort;
190    }
191  double sqrootD=sqrt(D);
192  alpha1=fmin(((-b-sqrootD)/(2*a)),((-b+sqrootD)/(2*a)));
193  alpha2=fmax(((-b-sqrootD)/(2*a)),((-b+sqrootD)/(2*a)));
194  if ( alpha2<0 ){
195    alpha=fmax(alpha2, -maxstep);}
196  else if ( alpha1<0 ){
197    alpha=fmax(alpha1, -maxstep);}
198  else {
199    alpha=-maxstep;
200  }
201 }
202 else if (-c/b<0){
203   alpha=fmax(-c/b, -maxstep);}
204 else {
205   alpha=-maxstep;
206 }
207 /* Check correctness
208
209double theta=alpha;
210 double fprime;
211 if ( update ==1 ) {
212    fprime=(K*Wij+K*Wij-Aij-Aij-Sij-Sji)-
213     K*theta/(theta*theta*(Wii*Wjj-Wij*Wij)-1-2*theta*Wij)*(-2*(Wii*Wjj+Wij*Wij)
214   +2*theta*Wij*(Wii*Wjj-Wij*Wij));
215 }
216 else {
217   fprime=(-K*Wij-K*Wij+Aij+Aij-Sij-Sji)-K*theta/(theta*theta*(-Wii*Wjj+Wij*Wij)
218           +1-2*theta*Wij)*(2*(Wii*Wjj+Wij*Wij)
219           +2*theta*Wij*(Wii*Wjj-Wij*Wij));
220 }
221 if ((abs(theta)<maxstep) && (fabs(fprime) > 10e-6)){
223   fflush(stdout);
224    abort;
225    } */
226  return alpha;
227}
228
229
230void ICS::sinco(SOL& ICS_sol, const SOL& ICS_start, double tol ){
231  /*This is the main routine for computing the sparse inverse covariance.
232   ICS class contains the problem data, that is A, K, p, lambda and S.
233   The input for this routine is a starting point ICS_start which contains
234   intial C, initial W=C^{-1} and initial objective function value. The output
235of the routine is final C, W and function value. Tolerance is set in tol and
236if set too high may significantly slow down the solver, since convergence
237 is slow in the end.  */
238  vector<int> UpInd;
239  double alphamax, funmax, K, fnew;
240  int imax, jmax, updatemax,  iter=0;
241  UpInd.resize(p*p, 1);
242  bool stop=false;
243  double f=ICS_start.f;
244  Steps AllSteps(p,p);
245  K=N/2;
246  ICS_sol.W=ICS_start.W;
247  ICS_sol.C=ICS_start.C;
248  Matrix Cp(p,p,0);
249  Matrix Cpp(p,p,0);
252  for (int i=p-1; i>=0; i--) {
253    for (int j=i; j>=0; j--) {
254            if (ICS_start.C(i,j)>0) {
255               Cp(i,j)=ICS_start.C(i,j);
256               Cp(j,i)=ICS_start.C(j,i); }
257            else {
258               Cpp(i,j)=-ICS_start.C(i,j);
259               Cpp(j,i)=-ICS_start.C(j,i);
260            }
261     }
262  }
269
270
271  while (!stop && (iter<100000)) {
272    // %find the largest positive derivative
273      alphamax=0;
274      funmax=f;
275      iter+=1;
276      //      printf("\n iteration [%4d] [%6f]", iter, f);
277      //    fflush(stdout);
278
279      for (int i=p-1; i>=0; i--) {
280        for (int j=i; j>=0; j--) {
281          if (UpInd[i*p+j] || UpInd[i*p+i] || UpInd[j*p+j]){
282            AllSteps(i,j).alpha=0;
284              AllSteps(i,j).update=1;
285              AllSteps(i,j).alpha=findposstep(K, ICS_sol.W(i,j),ICS_sol.W(i,i),
286                                              ICS_sol.W(j,j), A(i,j), lambda*S(i,j),
287                                                lambda*S(j,i), AllSteps(i,j).update);
288            }
290              AllSteps(i,j).update=-1;
291              AllSteps(i,j).alpha=findposstep(K, ICS_sol.W(i,j),ICS_sol.W(i,i),
292                                              ICS_sol.W(j,j), A(i,j), lambda*S(i,j),
293                                                lambda*S(j,i), AllSteps(i,j).update);
294            }
296              AllSteps(i,j).update=1;
297              AllSteps(i,j).alpha=findnegstep(K, (i==j), ICS_sol.W(i,j),ICS_sol.W(i,i),
298                                              ICS_sol.W(j,j), A(i,j), lambda*S(i,j),
299                                              lambda*S(j,i), Cp(i,j), Cpp(i,j),
300                                              AllSteps(i,j).update);
301            }
303              AllSteps(i,j).update=-1;
304              AllSteps(i,j).alpha=findnegstep(K, (i==j), ICS_sol.W(i,j),ICS_sol.W(i,i),
305                                      ICS_sol.W(j,j), A(i,j), lambda*S(i,j),
306                                       lambda*S(j,i), Cp(i,j), Cpp(i,j), AllSteps(i,j).update);
307            }
308            if ( fabs(AllSteps(i,j).alpha)> tol) {
309               AllSteps(i,j).fchange=funvalue_update(AllSteps(i,j).alpha,
310                                              K,ICS_sol.W(i,j),ICS_sol.W(i,i),
311                                              ICS_sol.W(j,j), A(i,j), lambda*S(i,j),
312                                              lambda*S(j,i),  AllSteps(i,j).update);
313            }
314            else {
315               AllSteps(i,j).fchange=0;
316            }
317          }
318          fnew=f+AllSteps(i,j).fchange;
319          if (fnew > funmax ) {
320            funmax=fnew;
321            imax=i;
322            jmax=j;
323            alphamax=AllSteps(i,j).alpha;
324            updatemax=AllSteps(i,j).update;
325          }
326        }
327      }
328      /*      for (int i=0; i<p; i++) {
329      for (int j=0; j<p; j++) {
330        printf("\n [%6f] [%6f] [%6d] [%6d]", AllSteps(i,j).fchange, AllSteps(i,j).alpha, i, j);
331       }
332       } */
333    if ( (alphamax==0) || ((funmax-f)/(fabs(f)+1))<tol ) { stop=true;}
334    else {
335      if (updatemax==1) {
336        Cp(imax,jmax)=Cp(imax,jmax)+alphamax;
337        Cp(jmax,imax)=Cp(jmax,imax)+alphamax;
338      }
339      else {
340        Cpp(imax,jmax)=Cpp(imax,jmax)+alphamax;
341        Cpp(jmax,imax)=Cpp(jmax,imax)+alphamax;
342      }
343      if ((Cp(imax,jmax)<0) || (Cpp(imax, jmax)<0)) {
344        printf("\n illegal step");
345          }
346      f=funmax;
348      /*    for (int i=0; i<p; i++) {
349       for (int j=0; j<p; j++) {
351        }
352        } */
353    }
354  }
355
356  ICS_sol.C=Cp;
357  ICS_sol.C.sumwithscal(-1, Cpp);
358  ICS_sol.f=f;
359}  // end of the sinco routine.
360
361
362
363int main()
364{
365  /* the main routine for testing in C++ environment */
366  double v1, v2;
367  ICS test;
368  test.p=200;
369  int p=test.p;
370  Matrix A(p,p,1);
371  A.dim1=p;
372  A.dim2=p;
373  test.S=A;
374  Matrix Zeros(p,p,0);
375  Zeros.dim1=p;
376  Zeros.dim2=p;
377  SOL Init_Sol;
378  SOL  Opt_Sol;
379  Init_Sol.C=Zeros;
380  Init_Sol.W=Zeros;
381  Opt_Sol.C=Zeros;
382  Opt_Sol.W=Zeros;
383  test.A=Zeros;
384  test.lambda=200;
385  test.N=50;
386  for (int i=0; i<p; i++) {
387    for (int j=0; j<p; j++) {
388      A(i,j)=(1.0*(rand() % 100))/100.0;
389    }
390  }
391  ifstream myfile ("ECmat.txt");
392  /*  myfile.open(); */
393  for (int i=0; i<p; i++) {
394    for (int j=0; j<p; j++) {
395      myfile>>test.A(i,j);
396    }
397  }
398  /* myfile.close(); */
399  /*  for (int i=0; i<p; i++) {
400    for (int j=0; j<p; j++) {
401      for (int l=0; l<p; l++){
402        test.A(i,j)+=A(i,l)*A(j,l);
403      }
404      test.A(i,j)*=test.N/2;
405      printf("\n [%16f]", test.A(i,j));
406    }
407    } */
408
411  v1=Init_Sol.C.dotmultiply(test.A);
412  v2=Init_Sol.C.normS(test.S);
413  /*    for (int i=0; i<p; i++) {
414    for (int j=0; j<p; j++) {
415      printf("\n [%6f]  [%6f]  [%6f]  [%6d] [%6d]", Init_Sol.C(i,j),
416             test.A(i,j), test.S(i,j), i, j);
417    }
418    }
419    fflush(stdout); */
420  Init_Sol.f=-v1-test.lambda*v2;
421  double tol=1e-6;
422  for (int iter=1; iter<=20; iter++){
423    test.sinco(Opt_Sol, Init_Sol, tol );
424    printf("\n finished the run!\n ");
425    test.lambda-=10;
426    Init_Sol=Opt_Sol;
427    v2=Init_Sol.C.normS(test.S);
428    Init_Sol.f-=10*v2;
429  }
430   for (int i=0; i<p; i++) {
431    for (int j=0; j<p; j++) {
432      if ( Opt_Sol.C(i,j) !=0) {
433        printf("\n [%6f] [%4d] [%4d]", Opt_Sol.C(i,j), i, j);
434      }
435    }
436   }
437}
Note: See TracBrowser for help on using the repository browser.