SHOGUN  4.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
overlapping.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 #include <stdlib.h>
20 #include <string.h>
21 
22 void identifySomeZeroEntries(double * u, int * zeroGroupFlag, int *entrySignFlag,
23  int *pp, int *gg,
24  double *v, double lambda1, double lambda2,
25  int p, int g, double * w, double *G){
26 
27  int i, j, newZeroNum, iterStep=0;
28  double twoNorm, temp;
29 
30  /*
31  * process the L1 norm
32  *
33  * generate the u>=0, and assign values to entrySignFlag
34  *
35  */
36  for(i=0;i<p;i++){
37  if (v[i]> lambda1){
38  u[i]=v[i]-lambda1;
39 
40  entrySignFlag[i]=1;
41  }
42  else{
43  if (v[i] < -lambda1){
44  u[i]= -v[i] -lambda1;
45 
46  entrySignFlag[i]=-1;
47  }
48  else{
49  u[i]=0;
50 
51  entrySignFlag[i]=0;
52  }
53  }
54  }
55 
56  /*
57  * Applying Algorithm 1 for identifying some sparse groups
58  *
59  */
60 
61  /* zeroGroupFlag denotes whether the corresponding group is zero */
62  for(i=0;i<g;i++)
63  zeroGroupFlag[i]=1;
64 
65  while(1){
66 
67  iterStep++;
68 
69  if (iterStep>g+1){
70 
71  printf("\n Identify Zero Group: iterStep= %d. The code might have a bug! Check it!", iterStep);
72  return;
73  }
74 
75  /*record the number of newly detected sparse groups*/
76  newZeroNum=0;
77 
78  for (i=0;i<g;i++){
79 
80  if(zeroGroupFlag[i]){
81 
82  /*compute the two norm of the */
83 
84  twoNorm=0;
85  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
86  temp=u[ (int) G[j]];
87  twoNorm+=temp*temp;
88  }
89  twoNorm=sqrt(twoNorm);
90 
91  /*
92  printf("\n twoNorm=%2.5f, %2.5f",twoNorm,lambda2 * w[3*i+2]);
93  */
94 
95  /*
96  * Test whether this group should be sparse
97  */
98  if (twoNorm<= lambda2 * w[3*i+2] ){
99  zeroGroupFlag[i]=0;
100 
101  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++)
102  u[ (int) G[j]]=0;
103 
104  newZeroNum++;
105 
106  /*
107  printf("\n zero group=%d", i);
108  */
109  }
110  } /*end of if(!zeroGroupFlag[i]) */
111 
112  } /*end of for*/
113 
114  if (newZeroNum==0)
115  break;
116  }
117 
118  *pp=0;
119  /* zeroGroupFlag denotes whether the corresponding entry is zero */
120  for(i=0;i<p;i++){
121  if (u[i]==0){
122  entrySignFlag[i]=0;
123  *pp=*pp+1;
124  }
125  }
126 
127  *gg=0;
128  for(i=0;i<g;i++){
129  if (zeroGroupFlag[i]==0)
130  *gg=*gg+1;
131  }
132 }
133 
134 void xFromY(double *x, double *y,
135  double *u, double *Y,
136  int p, int g, int *zeroGroupFlag,
137  double *G, double *w){
138 
139  int i,j;
140 
141 
142  for(i=0;i<p;i++)
143  x[i]=u[i];
144 
145  for(i=0;i<g;i++){
146  if(zeroGroupFlag[i]){ /*this group is non-zero*/
147 
148  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
149  x[ (int) G[j] ] -= Y[j];
150  }
151  }
152  }/*end of for(i=0;i<g;i++) */
153 
154  for(i=0;i<p;i++){
155  if (x[i]>=0){
156  y[i]=0;
157  }
158  else{
159  y[i]=x[i];
160  x[i]=0;
161  }
162  }
163 }
164 
165 void YFromx(double *Y,
166  double *xnew, double *Ynew,
167  double lambda2, int g, int *zeroGroupFlag,
168  double *G, double *w){
169 
170  int i, j;
171  double twoNorm, temp;
172 
173  for(i=0;i<g;i++){
174  if(zeroGroupFlag[i]){ /*this group is non-zero*/
175 
176  twoNorm=0;
177  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
178  temp=xnew[ (int) G[j] ];
179 
180  Y[j]=temp;
181 
182  twoNorm+=temp*temp;
183  }
184  twoNorm=sqrt(twoNorm); /* two norm for x_{G_i}*/
185 
186  if (twoNorm > 0 ){ /*if x_{G_i} is non-zero*/
187  temp=lambda2 * w[3*i+2] / twoNorm;
188 
189  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++)
190  Y[j] *= temp;
191  }
192  else /*if x_{G_j} =0, we let Y^i=Ynew^i*/
193  {
194  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++)
195  Y[j]=Ynew[j];
196  }
197  }
198  }/*end of for(i=0;i<g;i++) */
199 }
200 
201 void dualityGap(double *gap, double *penalty2,
202  double *x, double *Y, int g, int *zeroGroupFlag,
203  double *G, double *w, double lambda2){
204 
205  int i,j;
206  double temp, twoNorm, innerProduct;
207 
208  *gap=0; *penalty2=0;
209 
210  for(i=0;i<g;i++){
211  if(zeroGroupFlag[i]){ /*this group is non-zero*/
212 
213  twoNorm=0;innerProduct=0;
214 
215  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
216  temp=x[ (int) G[j] ];
217 
218  twoNorm+=temp*temp;
219 
220  innerProduct+=temp * Y[j];
221  }
222 
223  twoNorm=sqrt(twoNorm)* w[3*i +2];
224 
225  *penalty2+=twoNorm;
226 
227  twoNorm=lambda2 * twoNorm;
228  if (twoNorm > innerProduct)
229  *gap+=twoNorm-innerProduct;
230  }
231  }/*end of for(i=0;i<g;i++) */
232 }
233 
234 void overlapping_gd(double *x, double *gap, double *penalty2,
235  double *v, int p, int g, double lambda1, double lambda2,
236  double *w, double *G, double *Y, int maxIter, int flag, double tol){
237 
238  int YSize=(int) w[3*(g-1) +1]+1;
239  double *u=(double *)malloc(sizeof(double)*p);
240  double *y=(double *)malloc(sizeof(double)*p);
241 
242  double *xnew=(double *)malloc(sizeof(double)*p);
243  double *Ynew=(double *)malloc(sizeof(double)* YSize );
244 
245  int *zeroGroupFlag=(int *)malloc(sizeof(int)*g);
246  int *entrySignFlag=(int *)malloc(sizeof(int)*p);
247  int pp, gg;
248  int i, j, iterStep;
249  double twoNorm,temp, L=1, leftValue, rightValue, gapR, penalty2R;
250  int nextRestartStep=0;
251 
252  /*
253  * call the function to identify some zero entries
254  *
255  * entrySignFlag[i]=0 denotes that the corresponding entry is definitely zero
256  *
257  * zeroGroupFlag[i]=0 denotes that the corresponding group is definitely zero
258  *
259  */
260 
261  identifySomeZeroEntries(u, zeroGroupFlag, entrySignFlag,
262  &pp, &gg,
263  v, lambda1, lambda2,
264  p, g, w, G);
265 
266  penalty2[1]=pp;
267  penalty2[2]=gg;
268  /*store pp and gg to penalty2[1] and penalty2[2]*/
269 
270 
271  /*
272  *-------------------
273  * Process Y
274  *-------------------
275  * We make sure that Y is feasible
276  * and if x_i=0, then set Y_{ij}=0
277  */
278  for(i=0;i<g;i++){
279 
280  if(zeroGroupFlag[i]){ /*this group is non-zero*/
281 
282  /*compute the two norm of the group*/
283  twoNorm=0;
284 
285  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
286 
287  if (! u[ (int) G[j] ] )
288  Y[j]=0;
289 
290  twoNorm+=Y[j]*Y[j];
291  }
292  twoNorm=sqrt(twoNorm);
293 
294  if (twoNorm > lambda2 * w[3*i+2] ){
295  temp=lambda2 * w[3*i+2] / twoNorm;
296 
297  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++)
298  Y[j]*=temp;
299  }
300  }
301  else{ /*this group is zero*/
302  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++)
303  Y[j]=0;
304  }
305  }
306 
307  /*
308  * set Ynew to zero
309  *
310  * in the following processing, we only operator Y and Ynew in the
311  * possibly non-zero groups by "if(zeroGroupFlag[i])"
312  *
313  */
314  for(i=0;i<YSize;i++)
315  Ynew[i]=0;
316 
317  /*
318  * ------------------------------------
319  * Gradient Descent begins here
320  * ------------------------------------
321  */
322 
323  /*
324  * compute x=max(u-Y * e, 0);
325  *
326  */
327  xFromY(x, y, u, Y, p, g, zeroGroupFlag, G, w);
328 
329 
330  /*the main loop */
331 
332  for(iterStep=0;iterStep<maxIter;iterStep++){
333 
334 
335  /*
336  * the gradient at Y is
337  *
338  * omega'(Y)=-x e^T
339  *
340  * where x=max(u-Y * e, 0);
341  *
342  */
343 
344 
345  /*
346  * line search to find Ynew with appropriate L
347  */
348 
349  while (1){
350  /*
351  * compute
352  * Ynew = proj ( Y + x e^T / L )
353  */
354  for(i=0;i<g;i++){
355  if(zeroGroupFlag[i]){ /*this group is non-zero*/
356 
357  twoNorm=0;
358  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
359  Ynew[j]= Y[j] + x[ (int) G[j] ] / L;
360 
361  twoNorm+=Ynew[j]*Ynew[j];
362  }
363  twoNorm=sqrt(twoNorm);
364 
365  if (twoNorm > lambda2 * w[3*i+2] ){
366  temp=lambda2 * w[3*i+2] / twoNorm;
367 
368  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++)
369  Ynew[j]*=temp;
370  }
371  }
372  }/*end of for(i=0;i<g;i++) */
373 
374  /*
375  * compute xnew=max(u-Ynew * e, 0);
376  *
377  *void xFromY(double *x, double *y,
378  * double *u, double *Y,
379  * int p, int g, int *zeroGroupFlag,
380  * double *G, double *w)
381  */
382  xFromY(xnew, y, u, Ynew, p, g, zeroGroupFlag, G, w);
383 
384  /* test whether L is appropriate*/
385  leftValue=0;
386  for(i=0;i<p;i++){
387  if (entrySignFlag[i]){
388  temp=xnew[i]-x[i];
389 
390  leftValue+= temp * ( 0.5 * temp + y[i]);
391  }
392  }
393 
394  rightValue=0;
395  for(i=0;i<g;i++){
396  if(zeroGroupFlag[i]){ /*this group is non-zero*/
397 
398  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
399  temp=Ynew[j]-Y[j];
400 
401  rightValue+=temp * temp;
402  }
403  }
404  }/*end of for(i=0;i<g;i++) */
405  rightValue=rightValue/2;
406 
407  if ( leftValue <= L * rightValue){
408 
409  temp= L * rightValue / leftValue;
410 
411  if (temp >5)
412  L=L*0.8;
413 
414  break;
415  }
416  else{
417  temp=leftValue / rightValue;
418 
419  if (L*2 <= temp)
420  L=temp;
421  else
422  L=2*L;
423 
424 
425  if ( L / g - 2* g ){
426 
427  if (rightValue < 1e-16){
428  break;
429  }
430  else{
431 
432  printf("\n GD: leftValue=%e, rightValue=%e, ratio=%e", leftValue, rightValue, temp);
433 
434  printf("\n L=%e > 2 * %d * %d. There might be a bug here. Otherwise, it is due to numerical issue.", L, g, g);
435 
436  break;
437  }
438  }
439  }
440  }
441 
442  /* compute the duality gap at (xnew, Ynew)
443  *
444  * void dualityGap(double *gap, double *penalty2,
445  * double *x, double *Y, int g, int *zeroGroupFlag,
446  * double *G, double *w, double lambda2)
447  *
448  */
449  dualityGap(gap, penalty2, xnew, Ynew, g, zeroGroupFlag, G, w, lambda2);
450 
451  /*
452  * flag =1 means restart
453  *
454  * flag =0 means with restart
455  *
456  * nextRestartStep denotes the next "step number" for
457  * initializing the restart process.
458  *
459  * This is based on the fact that, the result is only beneficial when
460  * xnew is good. In other words,
461  * if xnew is not good, then the
462  * restart might not be helpful.
463  */
464 
465  if ( (flag==0) || (flag==1 && iterStep < nextRestartStep )){
466 
467  /* copy Ynew to Y, and xnew to x */
468  memcpy(x, xnew, sizeof(double) * p);
469  memcpy(Y, Ynew, sizeof(double) * YSize);
470 
471  /*
472  printf("\n iterStep=%d, L=%2.5f, gap=%e", iterStep, L, *gap);
473  */
474 
475  }
476  else{
477  /*
478  * flag=1
479  *
480  * We allow the restart of the program.
481  *
482  * Here, Y is constructed as a subgradient of xnew, based on the
483  * assumption that Y might be a better choice than Ynew, provided
484  * that xnew is good enough.
485  *
486  */
487 
488  /*
489  * compute the restarting point Y with xnew and Ynew
490  *
491  *void YFromx(double *Y,
492  * double *xnew, double *Ynew,
493  * double lambda2, int g, int *zeroGroupFlag,
494  * double *G, double *w)
495  */
496  YFromx(Y, xnew, Ynew, lambda2, g, zeroGroupFlag, G, w);
497 
498  /*compute the solution with the starting point Y
499  *
500  *void xFromY(double *x, double *y,
501  * double *u, double *Y,
502  * int p, int g, int *zeroGroupFlag,
503  * double *G, double *w)
504  *
505  */
506  xFromY(x, y, u, Y, p, g, zeroGroupFlag, G, w);
507 
508  /*compute the duality at (x, Y)
509  *
510  * void dualityGap(double *gap, double *penalty2,
511  * double *x, double *Y, int g, int *zeroGroupFlag,
512  * double *G, double *w, double lambda2)
513  *
514  */
515  dualityGap(&gapR, &penalty2R, x, Y, g, zeroGroupFlag, G, w, lambda2);
516 
517  if (*gap< gapR){
518  /*(xnew, Ynew) is better in terms of duality gap*/
519  /* copy Ynew to Y, and xnew to x */
520  memcpy(x, xnew, sizeof(double) * p);
521  memcpy(Y, Ynew, sizeof(double) * YSize);
522 
523  /*In this case, we do not apply restart, as (x,Y) is not better
524  *
525  * We postpone the "restart" by giving a
526  * "nextRestartStep"
527  */
528 
529  /*
530  * we test *gap here, in case *gap=0
531  */
532  if (*gap <=tol)
533  break;
534  else{
535  nextRestartStep=iterStep+ (int) sqrt(gapR / *gap);
536  }
537  }
538  else{ /*we use (x, Y), as it is better in terms of duality gap*/
539  *gap=gapR;
540  *penalty2=penalty2R;
541  }
542 
543  /*
544  printf("\n iterStep=%d, L=%2.5f, gap=%e, gapR=%e", iterStep, L, *gap, gapR);
545  */
546 
547  }
548 
549  /*
550  * if the duality gap is within pre-specified parameter tol
551  *
552  * we terminate the algorithm
553  */
554  if (*gap <=tol)
555  break;
556  }
557 
558  penalty2[3]=iterStep;
559 
560  penalty2[4]=0;
561  for(i=0;i<g;i++){
562  if (zeroGroupFlag[i]==0)
563  penalty2[4]=penalty2[4]+1;
564  else{
565  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
566  if (x[ (int) G[j] ] !=0)
567  break;
568  }
569 
570  if (j>(int) w[3*i +1])
571  penalty2[4]=penalty2[4]+1;
572  }
573  }
574 
575  /*
576  * assign sign to the solution x
577  */
578  for(i=0;i<p;i++){
579  if (entrySignFlag[i]==-1){
580  x[i]=-x[i];
581  }
582  }
583 
584  free (u);
585  free (y);
586  free (xnew);
587  free (Ynew);
588  free (zeroGroupFlag);
589  free (entrySignFlag);
590 }
591 
592 void gradientDescentStep(double *xnew, double *Ynew,
593  double *LL, double *u, double *y, int *entrySignFlag, double lambda2,
594  double *x, double *Y, int p, int g, int * zeroGroupFlag,
595  double *G, double *w){
596 
597  double twoNorm, temp, L=*LL, leftValue, rightValue;
598  int i,j;
599 
600 
601 
602  while (1){
603 
604  /*
605  * compute
606  * Ynew = proj ( Y + x e^T / L )
607  */
608  for(i=0;i<g;i++){
609  if(zeroGroupFlag[i]){ /*this group is non-zero*/
610 
611  twoNorm=0;
612  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
613  Ynew[j]= Y[j] + x[ (int) G[j] ] / L;
614 
615  twoNorm+=Ynew[j]*Ynew[j];
616  }
617  twoNorm=sqrt(twoNorm);
618 
619  if (twoNorm > lambda2 * w[3*i+2] ){
620  temp=lambda2 * w[3*i+2] / twoNorm;
621 
622  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++)
623  Ynew[j]*=temp;
624  }
625  }
626  }/*end of for(i=0;i<g;i++) */
627 
628  /*
629  * compute xnew=max(u-Ynew * e, 0);
630  *
631  *void xFromY(double *x, double *y,
632  * double *u, double *Y,
633  * int p, int g, int *zeroGroupFlag,
634  * double *G, double *w)
635  */
636  xFromY(xnew, y, u, Ynew, p, g, zeroGroupFlag, G, w);
637 
638  /* test whether L is appropriate*/
639  leftValue=0;
640  for(i=0;i<p;i++){
641  if (entrySignFlag[i]){
642  temp=xnew[i]-x[i];
643 
644  leftValue+= temp * ( 0.5 * temp + y[i]);
645  }
646  }
647 
648  rightValue=0;
649  for(i=0;i<g;i++){
650  if(zeroGroupFlag[i]){ /*this group is non-zero*/
651 
652  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
653  temp=Ynew[j]-Y[j];
654 
655  rightValue+=temp * temp;
656  }
657  }
658  }/*end of for(i=0;i<g;i++) */
659  rightValue=rightValue/2;
660 
661  /*
662  printf("\n leftValue =%e, rightValue=%e, L=%e", leftValue, rightValue, L);
663  */
664 
665  if ( leftValue <= L * rightValue){
666 
667  temp= L * rightValue / leftValue;
668 
669  if (temp >5)
670  L=L*0.8;
671 
672  break;
673  }
674  else{
675  temp=leftValue / rightValue;
676 
677  if (L*2 <= temp)
678  L=temp;
679  else
680  L=2*L;
681 
682  if ( L / g - 2* g >0 ){
683 
684  if (rightValue < 1e-16){
685  break;
686  }
687  else{
688 
689  printf("\n One Gradient Step: leftValue=%e, rightValue=%e, ratio=%e", leftValue, rightValue, temp);
690 
691  printf("\n L=%e > 2 * %d * %d. There might be a bug here. Otherwise, it is due to numerical issue.", L, g, g);
692 
693  break;
694  }
695  }
696  }
697  }
698 
699  *LL=L;
700 }
701 
702 void overlapping_agd(double *x, double *gap, double *penalty2,
703  double *v, int p, int g, double lambda1, double lambda2,
704  double *w, double *G, double *Y, int maxIter, int flag, double tol){
705 
706  int YSize=(int) w[3*(g-1) +1]+1;
707  double *u=(double *)malloc(sizeof(double)*p);
708  double *y=(double *)malloc(sizeof(double)*p);
709 
710  double *xnew=(double *)malloc(sizeof(double)*p);
711  double *Ynew=(double *)malloc(sizeof(double)* YSize );
712 
713  double *xS=(double *)malloc(sizeof(double)*p);
714  double *YS=(double *)malloc(sizeof(double)* YSize );
715 
716  /*double *xp=(double *)malloc(sizeof(double)*p);*/
717  double *Yp=(double *)malloc(sizeof(double)* YSize );
718 
719  int *zeroGroupFlag=(int *)malloc(sizeof(int)*g);
720  int *entrySignFlag=(int *)malloc(sizeof(int)*p);
721 
722  int pp, gg;
723  int i, j, iterStep;
724  double twoNorm,temp, L=1, leftValue, rightValue, gapR, penalty2R;
725  int nextRestartStep=0;
726 
727  double alpha, alphap=0.5, beta, gamma;
728 
729  /*
730  * call the function to identify some zero entries
731  *
732  * entrySignFlag[i]=0 denotes that the corresponding entry is definitely zero
733  *
734  * zeroGroupFlag[i]=0 denotes that the corresponding group is definitely zero
735  *
736  */
737 
738  identifySomeZeroEntries(u, zeroGroupFlag, entrySignFlag,
739  &pp, &gg,
740  v, lambda1, lambda2,
741  p, g, w, G);
742 
743  penalty2[1]=pp;
744  penalty2[2]=gg;
745  /*store pp and gg to penalty2[1] and penalty2[2]*/
746 
747  /*
748  *-------------------
749  * Process Y
750  *-------------------
751  * We make sure that Y is feasible
752  * and if x_i=0, then set Y_{ij}=0
753  */
754  for(i=0;i<g;i++){
755 
756  if(zeroGroupFlag[i]){ /*this group is non-zero*/
757 
758  /*compute the two norm of the group*/
759  twoNorm=0;
760 
761  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
762 
763  if (! u[ (int) G[j] ] )
764  Y[j]=0;
765 
766  twoNorm+=Y[j]*Y[j];
767  }
768  twoNorm=sqrt(twoNorm);
769 
770  if (twoNorm > lambda2 * w[3*i+2] ){
771  temp=lambda2 * w[3*i+2] / twoNorm;
772 
773  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++)
774  Y[j]*=temp;
775  }
776  }
777  else{ /*this group is zero*/
778  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++)
779  Y[j]=0;
780  }
781  }
782 
783  /*
784  * set Ynew and Yp to zero
785  *
786  * in the following processing, we only operate, Yp, Y and Ynew in the
787  * possibly non-zero groups by "if(zeroGroupFlag[i])"
788  *
789  */
790  for(i=0;i<YSize;i++)
791  YS[i]=Yp[i]=Ynew[i]=0;
792 
793 
794  /*
795  * ---------------
796  *
797  * we first do a gradient descent step for determing the value of an approporate L
798  *
799  * Also, we initialize gamma
800  *
801  * with Y, we compute a new Ynew
802  *
803  */
804 
805 
806  /*
807  * compute x=max(u-Y * e, 0);
808  */
809  xFromY(x, y, u, Y, p, g, zeroGroupFlag, G, w);
810 
811  /*
812  * compute (xnew, Ynew) from (x, Y)
813  *
814  *
815  * gradientDescentStep(double *xnew, double *Ynew,
816  double *LL, double *u, double *y, int *entrySignFlag, double lambda2,
817  double *x, double *Y, int p, int g, int * zeroGroupFlag,
818  double *G, double *w)
819  */
820 
821  gradientDescentStep(xnew, Ynew,
822  &L, u, y,entrySignFlag,lambda2,
823  x, Y, p, g, zeroGroupFlag,
824  G, w);
825 
826  /*
827  * we have finished one gradient descent to get
828  *
829  * (x, Y) and (xnew, Ynew), where (xnew, Ynew) is
830  *
831  * a gradient descent step based on (x, Y)
832  *
833  * we set (xp, Yp)=(x, Y)
834  *
835  * (x, Y)= (xnew, Ynew)
836  */
837 
838  /*memcpy(xp, x, sizeof(double) * p);*/
839  memcpy(Yp, Y, sizeof(double) * YSize);
840 
841  /*memcpy(x, xnew, sizeof(double) * p);*/
842  memcpy(Y, Ynew, sizeof(double) * YSize);
843 
844  gamma=L;
845 
846  /*
847  * ------------------------------------
848  * Accelerated Gradient Descent begins here
849  * ------------------------------------
850  */
851 
852 
853  for(iterStep=0;iterStep<maxIter;iterStep++){
854 
855 
856  while (1){
857 
858 
859  /*
860  * compute alpha as the positive root of
861  *
862  * L * alpha^2 = (1-alpha) * gamma
863  *
864  */
865 
866  alpha= ( - gamma + sqrt( gamma * gamma + 4 * L * gamma ) ) / 2 / L;
867 
868  beta= gamma * (1-alphap)/ alphap / (gamma + L * alpha);
869 
870  /*
871  * compute YS= Y + beta * (Y - Yp)
872  *
873  */
874  for(i=0;i<g;i++){
875  if(zeroGroupFlag[i]){ /*this group is non-zero*/
876 
877  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
878 
879  YS[j]=Y[j] + beta * (Y[j]-Yp[j]);
880 
881  }
882  }
883  }/*end of for(i=0;i<g;i++) */
884 
885 
886  /*
887  * compute xS
888  */
889  xFromY(xS, y, u, YS, p, g, zeroGroupFlag, G, w);
890 
891 
892  /*
893  *
894  * Ynew = proj ( YS + xS e^T / L )
895  *
896  */
897  for(i=0;i<g;i++){
898  if(zeroGroupFlag[i]){ /*this group is non-zero*/
899 
900  twoNorm=0;
901  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
902 
903  Ynew[j]= YS[j] + xS[ (int) G[j] ] / L;
904 
905  twoNorm+=Ynew[j]*Ynew[j];
906  }
907  twoNorm=sqrt(twoNorm);
908 
909  if (twoNorm > lambda2 * w[3*i+2] ){
910  temp=lambda2 * w[3*i+2] / twoNorm;
911 
912  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++)
913  Ynew[j]*=temp;
914  }
915  }
916  }/*end of for(i=0;i<g;i++) */
917 
918  /*
919  * compute xnew=max(u-Ynew * e, 0);
920  *
921  *void xFromY(double *x, double *y,
922  * double *u, double *Y,
923  * int p, int g, int *zeroGroupFlag,
924  * double *G, double *w)
925  */
926 
927  xFromY(xnew, y, u, Ynew, p, g, zeroGroupFlag, G, w);
928 
929  /* test whether L is appropriate*/
930  leftValue=0;
931  for(i=0;i<p;i++){
932  if (entrySignFlag[i]){
933  temp=xnew[i]-xS[i];
934 
935  leftValue+= temp * ( 0.5 * temp + y[i]);
936  }
937  }
938 
939  rightValue=0;
940  for(i=0;i<g;i++){
941  if(zeroGroupFlag[i]){ /*this group is non-zero*/
942 
943  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
944  temp=Ynew[j]-YS[j];
945 
946  rightValue+=temp * temp;
947  }
948  }
949  }/*end of for(i=0;i<g;i++) */
950  rightValue=rightValue/2;
951 
952  if ( leftValue <= L * rightValue){
953 
954  temp= L * rightValue / leftValue;
955 
956  if (temp >5)
957  L=L*0.8;
958 
959  break;
960  }
961  else{
962  temp=leftValue / rightValue;
963 
964  if (L*2 <= temp)
965  L=temp;
966  else
967  L=2*L;
968 
969 
970 
971  if ( L / g - 2* g >0 ){
972 
973  if (rightValue < 1e-16){
974  break;
975  }
976  else{
977 
978  printf("\n AGD: leftValue=%e, rightValue=%e, ratio=%e", leftValue, rightValue, temp);
979 
980  printf("\n L=%e > 2 * %d * %d. There might be a bug here. Otherwise, it is due to numerical issue.", L, g, g);
981 
982  break;
983  }
984  }
985  }
986  }
987 
988  /* compute the duality gap at (xnew, Ynew)
989  *
990  * void dualityGap(double *gap, double *penalty2,
991  * double *x, double *Y, int g, int *zeroGroupFlag,
992  * double *G, double *w, double lambda2)
993  *
994  */
995  dualityGap(gap, penalty2,
996  xnew, Ynew, g, zeroGroupFlag,
997  G, w, lambda2);
998 
999 
1000  /*
1001  * if the duality gap is within pre-specified parameter tol
1002  *
1003  * we terminate the algorithm
1004  */
1005  if (*gap <=tol){
1006 
1007  memcpy(x, xnew, sizeof(double) * p);
1008  memcpy(Y, Ynew, sizeof(double) * YSize);
1009 
1010  break;
1011  }
1012 
1013 
1014 
1015  /*
1016  * flag =1 means restart
1017  *
1018  * flag =0 means with restart
1019  *
1020  * nextRestartStep denotes the next "step number" for
1021  * initializing the restart process.
1022  *
1023  * This is based on the fact that, the result is only beneficial when
1024  * xnew is good. In other words,
1025  * if xnew is not good, then the
1026  * restart might not be helpful.
1027  */
1028 
1029  if ( (flag==0) || (flag==1 && iterStep < nextRestartStep )){
1030 
1031 
1032  /*memcpy(xp, x, sizeof(double) * p);*/
1033  memcpy(Yp, Y, sizeof(double) * YSize);
1034 
1035  /*memcpy(x, xnew, sizeof(double) * p);*/
1036  memcpy(Y, Ynew, sizeof(double) * YSize);
1037 
1038  gamma=gamma * (1-alpha);
1039 
1040  alphap=alpha;
1041 
1042  /*
1043  printf("\n iterStep=%d, L=%2.5f, gap=%e", iterStep, L, *gap);
1044  */
1045 
1046  }
1047  else{
1048  /*
1049  * flag=1
1050  *
1051  * We allow the restart of the program.
1052  *
1053  * Here, Y is constructed as a subgradient of xnew, based on the
1054  * assumption that Y might be a better choice than Ynew, provided
1055  * that xnew is good enough.
1056  *
1057  */
1058 
1059  /*
1060  * compute the restarting point YS with xnew and Ynew
1061  *
1062  *void YFromx(double *Y,
1063  * double *xnew, double *Ynew,
1064  * double lambda2, int g, int *zeroGroupFlag,
1065  * double *G, double *w)
1066  */
1067  YFromx(YS, xnew, Ynew, lambda2, g, zeroGroupFlag, G, w);
1068 
1069  /*compute the solution with the starting point YS
1070  *
1071  *void xFromY(double *x, double *y,
1072  * double *u, double *Y,
1073  * int p, int g, int *zeroGroupFlag,
1074  * double *G, double *w)
1075  *
1076  */
1077  xFromY(xS, y, u, YS, p, g, zeroGroupFlag, G, w);
1078 
1079  /*compute the duality at (xS, YS)
1080  *
1081  * void dualityGap(double *gap, double *penalty2,
1082  * double *x, double *Y, int g, int *zeroGroupFlag,
1083  * double *G, double *w, double lambda2)
1084  *
1085  */
1086  dualityGap(&gapR, &penalty2R, xS, YS, g, zeroGroupFlag, G, w, lambda2);
1087 
1088  if (*gap< gapR){
1089  /*(xnew, Ynew) is better in terms of duality gap*/
1090 
1091  /*In this case, we do not apply restart, as (xS,YS) is not better
1092  *
1093  * We postpone the "restart" by giving a
1094  * "nextRestartStep"
1095  */
1096 
1097  /*memcpy(xp, x, sizeof(double) * p);*/
1098  memcpy(Yp, Y, sizeof(double) * YSize);
1099 
1100  /*memcpy(x, xnew, sizeof(double) * p);*/
1101  memcpy(Y, Ynew, sizeof(double) * YSize);
1102 
1103  gamma=gamma * (1-alpha);
1104 
1105  alphap=alpha;
1106 
1107  nextRestartStep=iterStep+ (int) sqrt(gapR / *gap);
1108  }
1109  else{
1110  /*we use (xS, YS), as it is better in terms of duality gap*/
1111 
1112  *gap=gapR;
1113  *penalty2=penalty2R;
1114 
1115  if (*gap <=tol){
1116 
1117  memcpy(x, xS, sizeof(double) * p);
1118  memcpy(Y, YS, sizeof(double) * YSize);
1119 
1120  break;
1121  }else{
1122  /*
1123  * we do a gradient descent based on (xS, YS)
1124  *
1125  */
1126 
1127  /*
1128  * compute (x, Y) from (xS, YS)
1129  *
1130  *
1131  * gradientDescentStep(double *xnew, double *Ynew,
1132  * double *LL, double *u, double *y, int *entrySignFlag, double lambda2,
1133  * double *x, double *Y, int p, int g, int * zeroGroupFlag,
1134  * double *G, double *w)
1135  */
1136  gradientDescentStep(x, Y,
1137  &L, u, y, entrySignFlag,lambda2,
1138  xS, YS, p, g, zeroGroupFlag,
1139  G, w);
1140 
1141  /*memcpy(xp, xS, sizeof(double) * p);*/
1142  memcpy(Yp, YS, sizeof(double) * YSize);
1143 
1144  gamma=L;
1145 
1146  alphap=0.5;
1147 
1148  }
1149 
1150 
1151  }
1152 
1153  /*
1154  * printf("\n iterStep=%d, L=%2.5f, gap=%e, gapR=%e", iterStep, L, *gap, gapR);
1155  */
1156 
1157  }/* flag =1*/
1158 
1159  } /* main loop */
1160 
1161 
1162 
1163  penalty2[3]=iterStep+1;
1164 
1165  /*
1166  * get the number of nonzero groups
1167  */
1168 
1169  penalty2[4]=0;
1170  for(i=0;i<g;i++){
1171  if (zeroGroupFlag[i]==0)
1172  penalty2[4]=penalty2[4]+1;
1173  else{
1174  for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){
1175  if (x[ (int) G[j] ] !=0)
1176  break;
1177  }
1178 
1179  if (j>(int) w[3*i +1])
1180  penalty2[4]=penalty2[4]+1;
1181  }
1182  }
1183 
1184 
1185  /*
1186  * assign sign to the solution x
1187  */
1188  for(i=0;i<p;i++){
1189  if (entrySignFlag[i]==-1){
1190  x[i]=-x[i];
1191  }
1192  }
1193 
1194  free (u);
1195  free (y);
1196 
1197  free (xnew);
1198  free (Ynew);
1199 
1200  free (xS);
1201  free (YS);
1202 
1203  /*free (xp);*/
1204  free (Yp);
1205 
1206  free (zeroGroupFlag);
1207  free (entrySignFlag);
1208 }
1209 
1210 void overlapping(double *x, double *gap, double *penalty2,
1211  double *v, int p, int g, double lambda1, double lambda2,
1212  double *w, double *G, double *Y, int maxIter, int flag, double tol){
1213 
1214  switch(flag){
1215  case 0:
1216  case 1:
1217  overlapping_gd(x, gap, penalty2,
1218  v, p, g, lambda1, lambda2,
1219  w, G, Y, maxIter, flag,tol);
1220  break;
1221  case 2:
1222  case 3:
1223 
1224  overlapping_agd(x, gap, penalty2,
1225  v, p, g, lambda1, lambda2,
1226  w, G, Y, maxIter, flag-2,tol);
1227 
1228  break;
1229  default:
1230  /* printf("\n Wrong flag! The value of flag should be 0,1,2,3. The program uses flag=2.");*/
1231 
1232  overlapping_agd(x, gap, penalty2,
1233  v, p, g, lambda1, lambda2,
1234  w, G, Y, maxIter, 0,tol);
1235  break;
1236  }
1237 
1238 
1239 }

SHOGUN Machine Learning Toolbox - Documentation