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

SHOGUN Machine Learning Toolbox - Documentation