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