SHOGUN
v2.0.0
|
00001 /* This program is free software: you can redistribute it and/or modify 00002 * it under the terms of the GNU General Public License as published by 00003 * the Free Software Foundation, either version 3 of the License, or 00004 * (at your option) any later version. 00005 * 00006 * This program is distributed in the hope that it will be useful, 00007 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00008 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00009 * GNU General Public License for more details. 00010 * 00011 * You should have received a copy of the GNU General Public License 00012 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00013 * 00014 * Copyright (C) 2009 - 2012 Jun Liu and Jieping Ye 00015 */ 00016 00017 #include <shogun/lib/slep/overlapping/overlapping.h> 00018 00019 void identifySomeZeroEntries(double * u, int * zeroGroupFlag, int *entrySignFlag, 00020 int *pp, int *gg, 00021 double *v, double lambda1, double lambda2, 00022 int p, int g, double * w, double *G){ 00023 00024 int i, j, newZeroNum, iterStep=0; 00025 double twoNorm, temp; 00026 00027 /* 00028 * process the L1 norm 00029 * 00030 * generate the u>=0, and assign values to entrySignFlag 00031 * 00032 */ 00033 for(i=0;i<p;i++){ 00034 if (v[i]> lambda1){ 00035 u[i]=v[i]-lambda1; 00036 00037 entrySignFlag[i]=1; 00038 } 00039 else{ 00040 if (v[i] < -lambda1){ 00041 u[i]= -v[i] -lambda1; 00042 00043 entrySignFlag[i]=-1; 00044 } 00045 else{ 00046 u[i]=0; 00047 00048 entrySignFlag[i]=0; 00049 } 00050 } 00051 } 00052 00053 /* 00054 * Applying Algorithm 1 for identifying some sparse groups 00055 * 00056 */ 00057 00058 /* zeroGroupFlag denotes whether the corresponding group is zero */ 00059 for(i=0;i<g;i++) 00060 zeroGroupFlag[i]=1; 00061 00062 while(1){ 00063 00064 iterStep++; 00065 00066 if (iterStep>g+1){ 00067 00068 printf("\n Identify Zero Group: iterStep= %d. The code might have a bug! Check it!", iterStep); 00069 return; 00070 } 00071 00072 /*record the number of newly detected sparse groups*/ 00073 newZeroNum=0; 00074 00075 for (i=0;i<g;i++){ 00076 00077 if(zeroGroupFlag[i]){ 00078 00079 /*compute the two norm of the */ 00080 00081 twoNorm=0; 00082 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00083 temp=u[ (int) G[j]]; 00084 twoNorm+=temp*temp; 00085 } 00086 twoNorm=sqrt(twoNorm); 00087 00088 /* 00089 printf("\n twoNorm=%2.5f, %2.5f",twoNorm,lambda2 * w[3*i+2]); 00090 */ 00091 00092 /* 00093 * Test whether this group should be sparse 00094 */ 00095 if (twoNorm<= lambda2 * w[3*i+2] ){ 00096 zeroGroupFlag[i]=0; 00097 00098 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++) 00099 u[ (int) G[j]]=0; 00100 00101 newZeroNum++; 00102 00103 /* 00104 printf("\n zero group=%d", i); 00105 */ 00106 } 00107 } /*end of if(!zeroGroupFlag[i]) */ 00108 00109 } /*end of for*/ 00110 00111 if (newZeroNum==0) 00112 break; 00113 } 00114 00115 *pp=0; 00116 /* zeroGroupFlag denotes whether the corresponding entry is zero */ 00117 for(i=0;i<p;i++){ 00118 if (u[i]==0){ 00119 entrySignFlag[i]=0; 00120 *pp=*pp+1; 00121 } 00122 } 00123 00124 *gg=0; 00125 for(i=0;i<g;i++){ 00126 if (zeroGroupFlag[i]==0) 00127 *gg=*gg+1; 00128 } 00129 } 00130 00131 void xFromY(double *x, double *y, 00132 double *u, double *Y, 00133 int p, int g, int *zeroGroupFlag, 00134 double *G, double *w){ 00135 00136 int i,j; 00137 00138 00139 for(i=0;i<p;i++) 00140 x[i]=u[i]; 00141 00142 for(i=0;i<g;i++){ 00143 if(zeroGroupFlag[i]){ /*this group is non-zero*/ 00144 00145 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00146 x[ (int) G[j] ] -= Y[j]; 00147 } 00148 } 00149 }/*end of for(i=0;i<g;i++) */ 00150 00151 for(i=0;i<p;i++){ 00152 if (x[i]>=0){ 00153 y[i]=0; 00154 } 00155 else{ 00156 y[i]=x[i]; 00157 x[i]=0; 00158 } 00159 } 00160 } 00161 00162 void YFromx(double *Y, 00163 double *xnew, double *Ynew, 00164 double lambda2, int g, int *zeroGroupFlag, 00165 double *G, double *w){ 00166 00167 int i, j; 00168 double twoNorm, temp; 00169 00170 for(i=0;i<g;i++){ 00171 if(zeroGroupFlag[i]){ /*this group is non-zero*/ 00172 00173 twoNorm=0; 00174 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00175 temp=xnew[ (int) G[j] ]; 00176 00177 Y[j]=temp; 00178 00179 twoNorm+=temp*temp; 00180 } 00181 twoNorm=sqrt(twoNorm); /* two norm for x_{G_i}*/ 00182 00183 if (twoNorm > 0 ){ /*if x_{G_i} is non-zero*/ 00184 temp=lambda2 * w[3*i+2] / twoNorm; 00185 00186 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++) 00187 Y[j] *= temp; 00188 } 00189 else /*if x_{G_j} =0, we let Y^i=Ynew^i*/ 00190 { 00191 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++) 00192 Y[j]=Ynew[j]; 00193 } 00194 } 00195 }/*end of for(i=0;i<g;i++) */ 00196 } 00197 00198 void dualityGap(double *gap, double *penalty2, 00199 double *x, double *Y, int g, int *zeroGroupFlag, 00200 double *G, double *w, double lambda2){ 00201 00202 int i,j; 00203 double temp, twoNorm, innerProduct; 00204 00205 *gap=0; *penalty2=0; 00206 00207 for(i=0;i<g;i++){ 00208 if(zeroGroupFlag[i]){ /*this group is non-zero*/ 00209 00210 twoNorm=0;innerProduct=0; 00211 00212 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00213 temp=x[ (int) G[j] ]; 00214 00215 twoNorm+=temp*temp; 00216 00217 innerProduct+=temp * Y[j]; 00218 } 00219 00220 twoNorm=sqrt(twoNorm)* w[3*i +2]; 00221 00222 *penalty2+=twoNorm; 00223 00224 twoNorm=lambda2 * twoNorm; 00225 if (twoNorm > innerProduct) 00226 *gap+=twoNorm-innerProduct; 00227 } 00228 }/*end of for(i=0;i<g;i++) */ 00229 } 00230 00231 void overlapping_gd(double *x, double *gap, double *penalty2, 00232 double *v, int p, int g, double lambda1, double lambda2, 00233 double *w, double *G, double *Y, int maxIter, int flag, double tol){ 00234 00235 int YSize=(int) w[3*(g-1) +1]+1; 00236 double *u=(double *)malloc(sizeof(double)*p); 00237 double *y=(double *)malloc(sizeof(double)*p); 00238 00239 double *xnew=(double *)malloc(sizeof(double)*p); 00240 double *Ynew=(double *)malloc(sizeof(double)* YSize ); 00241 00242 int *zeroGroupFlag=(int *)malloc(sizeof(int)*g); 00243 int *entrySignFlag=(int *)malloc(sizeof(int)*p); 00244 int pp, gg; 00245 int i, j, iterStep; 00246 double twoNorm,temp, L=1, leftValue, rightValue, gapR, penalty2R; 00247 int nextRestartStep=0; 00248 00249 /* 00250 * call the function to identify some zero entries 00251 * 00252 * entrySignFlag[i]=0 denotes that the corresponding entry is definitely zero 00253 * 00254 * zeroGroupFlag[i]=0 denotes that the corresponding group is definitely zero 00255 * 00256 */ 00257 00258 identifySomeZeroEntries(u, zeroGroupFlag, entrySignFlag, 00259 &pp, &gg, 00260 v, lambda1, lambda2, 00261 p, g, w, G); 00262 00263 penalty2[1]=pp; 00264 penalty2[2]=gg; 00265 /*store pp and gg to penalty2[1] and penalty2[2]*/ 00266 00267 00268 /* 00269 *------------------- 00270 * Process Y 00271 *------------------- 00272 * We make sure that Y is feasible 00273 * and if x_i=0, then set Y_{ij}=0 00274 */ 00275 for(i=0;i<g;i++){ 00276 00277 if(zeroGroupFlag[i]){ /*this group is non-zero*/ 00278 00279 /*compute the two norm of the group*/ 00280 twoNorm=0; 00281 00282 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00283 00284 if (! u[ (int) G[j] ] ) 00285 Y[j]=0; 00286 00287 twoNorm+=Y[j]*Y[j]; 00288 } 00289 twoNorm=sqrt(twoNorm); 00290 00291 if (twoNorm > lambda2 * w[3*i+2] ){ 00292 temp=lambda2 * w[3*i+2] / twoNorm; 00293 00294 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++) 00295 Y[j]*=temp; 00296 } 00297 } 00298 else{ /*this group is zero*/ 00299 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++) 00300 Y[j]=0; 00301 } 00302 } 00303 00304 /* 00305 * set Ynew to zero 00306 * 00307 * in the following processing, we only operator Y and Ynew in the 00308 * possibly non-zero groups by "if(zeroGroupFlag[i])" 00309 * 00310 */ 00311 for(i=0;i<YSize;i++) 00312 Ynew[i]=0; 00313 00314 /* 00315 * ------------------------------------ 00316 * Gradient Descent begins here 00317 * ------------------------------------ 00318 */ 00319 00320 /* 00321 * compute x=max(u-Y * e, 0); 00322 * 00323 */ 00324 xFromY(x, y, u, Y, p, g, zeroGroupFlag, G, w); 00325 00326 00327 /*the main loop */ 00328 00329 for(iterStep=0;iterStep<maxIter;iterStep++){ 00330 00331 00332 /* 00333 * the gradient at Y is 00334 * 00335 * omega'(Y)=-x e^T 00336 * 00337 * where x=max(u-Y * e, 0); 00338 * 00339 */ 00340 00341 00342 /* 00343 * line search to find Ynew with appropriate L 00344 */ 00345 00346 while (1){ 00347 /* 00348 * compute 00349 * Ynew = proj ( Y + x e^T / L ) 00350 */ 00351 for(i=0;i<g;i++){ 00352 if(zeroGroupFlag[i]){ /*this group is non-zero*/ 00353 00354 twoNorm=0; 00355 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00356 Ynew[j]= Y[j] + x[ (int) G[j] ] / L; 00357 00358 twoNorm+=Ynew[j]*Ynew[j]; 00359 } 00360 twoNorm=sqrt(twoNorm); 00361 00362 if (twoNorm > lambda2 * w[3*i+2] ){ 00363 temp=lambda2 * w[3*i+2] / twoNorm; 00364 00365 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++) 00366 Ynew[j]*=temp; 00367 } 00368 } 00369 }/*end of for(i=0;i<g;i++) */ 00370 00371 /* 00372 * compute xnew=max(u-Ynew * e, 0); 00373 * 00374 *void xFromY(double *x, double *y, 00375 * double *u, double *Y, 00376 * int p, int g, int *zeroGroupFlag, 00377 * double *G, double *w) 00378 */ 00379 xFromY(xnew, y, u, Ynew, p, g, zeroGroupFlag, G, w); 00380 00381 /* test whether L is appropriate*/ 00382 leftValue=0; 00383 for(i=0;i<p;i++){ 00384 if (entrySignFlag[i]){ 00385 temp=xnew[i]-x[i]; 00386 00387 leftValue+= temp * ( 0.5 * temp + y[i]); 00388 } 00389 } 00390 00391 rightValue=0; 00392 for(i=0;i<g;i++){ 00393 if(zeroGroupFlag[i]){ /*this group is non-zero*/ 00394 00395 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00396 temp=Ynew[j]-Y[j]; 00397 00398 rightValue+=temp * temp; 00399 } 00400 } 00401 }/*end of for(i=0;i<g;i++) */ 00402 rightValue=rightValue/2; 00403 00404 if ( leftValue <= L * rightValue){ 00405 00406 temp= L * rightValue / leftValue; 00407 00408 if (temp >5) 00409 L=L*0.8; 00410 00411 break; 00412 } 00413 else{ 00414 temp=leftValue / rightValue; 00415 00416 if (L*2 <= temp) 00417 L=temp; 00418 else 00419 L=2*L; 00420 00421 00422 if ( L / g - 2* g ){ 00423 00424 if (rightValue < 1e-16){ 00425 break; 00426 } 00427 else{ 00428 00429 printf("\n GD: leftValue=%e, rightValue=%e, ratio=%e", leftValue, rightValue, temp); 00430 00431 printf("\n L=%e > 2 * %d * %d. There might be a bug here. Otherwise, it is due to numerical issue.", L, g, g); 00432 00433 break; 00434 } 00435 } 00436 } 00437 } 00438 00439 /* compute the duality gap at (xnew, Ynew) 00440 * 00441 * void dualityGap(double *gap, double *penalty2, 00442 * double *x, double *Y, int g, int *zeroGroupFlag, 00443 * double *G, double *w, double lambda2) 00444 * 00445 */ 00446 dualityGap(gap, penalty2, xnew, Ynew, g, zeroGroupFlag, G, w, lambda2); 00447 00448 /* 00449 * flag =1 means restart 00450 * 00451 * flag =0 means with restart 00452 * 00453 * nextRestartStep denotes the next "step number" for 00454 * initializing the restart process. 00455 * 00456 * This is based on the fact that, the result is only beneficial when 00457 * xnew is good. In other words, 00458 * if xnew is not good, then the 00459 * restart might not be helpful. 00460 */ 00461 00462 if ( (flag==0) || (flag==1 && iterStep < nextRestartStep )){ 00463 00464 /* copy Ynew to Y, and xnew to x */ 00465 memcpy(x, xnew, sizeof(double) * p); 00466 memcpy(Y, Ynew, sizeof(double) * YSize); 00467 00468 /* 00469 printf("\n iterStep=%d, L=%2.5f, gap=%e", iterStep, L, *gap); 00470 */ 00471 00472 } 00473 else{ 00474 /* 00475 * flag=1 00476 * 00477 * We allow the restart of the program. 00478 * 00479 * Here, Y is constructed as a subgradient of xnew, based on the 00480 * assumption that Y might be a better choice than Ynew, provided 00481 * that xnew is good enough. 00482 * 00483 */ 00484 00485 /* 00486 * compute the restarting point Y with xnew and Ynew 00487 * 00488 *void YFromx(double *Y, 00489 * double *xnew, double *Ynew, 00490 * double lambda2, int g, int *zeroGroupFlag, 00491 * double *G, double *w) 00492 */ 00493 YFromx(Y, xnew, Ynew, lambda2, g, zeroGroupFlag, G, w); 00494 00495 /*compute the solution with the starting point Y 00496 * 00497 *void xFromY(double *x, double *y, 00498 * double *u, double *Y, 00499 * int p, int g, int *zeroGroupFlag, 00500 * double *G, double *w) 00501 * 00502 */ 00503 xFromY(x, y, u, Y, p, g, zeroGroupFlag, G, w); 00504 00505 /*compute the duality at (x, Y) 00506 * 00507 * void dualityGap(double *gap, double *penalty2, 00508 * double *x, double *Y, int g, int *zeroGroupFlag, 00509 * double *G, double *w, double lambda2) 00510 * 00511 */ 00512 dualityGap(&gapR, &penalty2R, x, Y, g, zeroGroupFlag, G, w, lambda2); 00513 00514 if (*gap< gapR){ 00515 /*(xnew, Ynew) is better in terms of duality gap*/ 00516 /* copy Ynew to Y, and xnew to x */ 00517 memcpy(x, xnew, sizeof(double) * p); 00518 memcpy(Y, Ynew, sizeof(double) * YSize); 00519 00520 /*In this case, we do not apply restart, as (x,Y) is not better 00521 * 00522 * We postpone the "restart" by giving a 00523 * "nextRestartStep" 00524 */ 00525 00526 /* 00527 * we test *gap here, in case *gap=0 00528 */ 00529 if (*gap <=tol) 00530 break; 00531 else{ 00532 nextRestartStep=iterStep+ (int) sqrt(gapR / *gap); 00533 } 00534 } 00535 else{ /*we use (x, Y), as it is better in terms of duality gap*/ 00536 *gap=gapR; 00537 *penalty2=penalty2R; 00538 } 00539 00540 /* 00541 printf("\n iterStep=%d, L=%2.5f, gap=%e, gapR=%e", iterStep, L, *gap, gapR); 00542 */ 00543 00544 } 00545 00546 /* 00547 * if the duality gap is within pre-specified parameter tol 00548 * 00549 * we terminate the algorithm 00550 */ 00551 if (*gap <=tol) 00552 break; 00553 } 00554 00555 penalty2[3]=iterStep; 00556 00557 penalty2[4]=0; 00558 for(i=0;i<g;i++){ 00559 if (zeroGroupFlag[i]==0) 00560 penalty2[4]=penalty2[4]+1; 00561 else{ 00562 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00563 if (x[ (int) G[j] ] !=0) 00564 break; 00565 } 00566 00567 if (j>(int) w[3*i +1]) 00568 penalty2[4]=penalty2[4]+1; 00569 } 00570 } 00571 00572 /* 00573 * assign sign to the solution x 00574 */ 00575 for(i=0;i<p;i++){ 00576 if (entrySignFlag[i]==-1){ 00577 x[i]=-x[i]; 00578 } 00579 } 00580 00581 free (u); 00582 free (y); 00583 free (xnew); 00584 free (Ynew); 00585 free (zeroGroupFlag); 00586 free (entrySignFlag); 00587 } 00588 00589 void gradientDescentStep(double *xnew, double *Ynew, 00590 double *LL, double *u, double *y, int *entrySignFlag, double lambda2, 00591 double *x, double *Y, int p, int g, int * zeroGroupFlag, 00592 double *G, double *w){ 00593 00594 double twoNorm, temp, L=*LL, leftValue, rightValue; 00595 int i,j; 00596 00597 00598 00599 while (1){ 00600 00601 /* 00602 * compute 00603 * Ynew = proj ( Y + x e^T / L ) 00604 */ 00605 for(i=0;i<g;i++){ 00606 if(zeroGroupFlag[i]){ /*this group is non-zero*/ 00607 00608 twoNorm=0; 00609 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00610 Ynew[j]= Y[j] + x[ (int) G[j] ] / L; 00611 00612 twoNorm+=Ynew[j]*Ynew[j]; 00613 } 00614 twoNorm=sqrt(twoNorm); 00615 00616 if (twoNorm > lambda2 * w[3*i+2] ){ 00617 temp=lambda2 * w[3*i+2] / twoNorm; 00618 00619 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++) 00620 Ynew[j]*=temp; 00621 } 00622 } 00623 }/*end of for(i=0;i<g;i++) */ 00624 00625 /* 00626 * compute xnew=max(u-Ynew * e, 0); 00627 * 00628 *void xFromY(double *x, double *y, 00629 * double *u, double *Y, 00630 * int p, int g, int *zeroGroupFlag, 00631 * double *G, double *w) 00632 */ 00633 xFromY(xnew, y, u, Ynew, p, g, zeroGroupFlag, G, w); 00634 00635 /* test whether L is appropriate*/ 00636 leftValue=0; 00637 for(i=0;i<p;i++){ 00638 if (entrySignFlag[i]){ 00639 temp=xnew[i]-x[i]; 00640 00641 leftValue+= temp * ( 0.5 * temp + y[i]); 00642 } 00643 } 00644 00645 rightValue=0; 00646 for(i=0;i<g;i++){ 00647 if(zeroGroupFlag[i]){ /*this group is non-zero*/ 00648 00649 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00650 temp=Ynew[j]-Y[j]; 00651 00652 rightValue+=temp * temp; 00653 } 00654 } 00655 }/*end of for(i=0;i<g;i++) */ 00656 rightValue=rightValue/2; 00657 00658 /* 00659 printf("\n leftValue =%e, rightValue=%e, L=%e", leftValue, rightValue, L); 00660 */ 00661 00662 if ( leftValue <= L * rightValue){ 00663 00664 temp= L * rightValue / leftValue; 00665 00666 if (temp >5) 00667 L=L*0.8; 00668 00669 break; 00670 } 00671 else{ 00672 temp=leftValue / rightValue; 00673 00674 if (L*2 <= temp) 00675 L=temp; 00676 else 00677 L=2*L; 00678 00679 if ( L / g - 2* g >0 ){ 00680 00681 if (rightValue < 1e-16){ 00682 break; 00683 } 00684 else{ 00685 00686 printf("\n One Gradient Step: leftValue=%e, rightValue=%e, ratio=%e", leftValue, rightValue, temp); 00687 00688 printf("\n L=%e > 2 * %d * %d. There might be a bug here. Otherwise, it is due to numerical issue.", L, g, g); 00689 00690 break; 00691 } 00692 } 00693 } 00694 } 00695 00696 *LL=L; 00697 } 00698 00699 void overlapping_agd(double *x, double *gap, double *penalty2, 00700 double *v, int p, int g, double lambda1, double lambda2, 00701 double *w, double *G, double *Y, int maxIter, int flag, double tol){ 00702 00703 int YSize=(int) w[3*(g-1) +1]+1; 00704 double *u=(double *)malloc(sizeof(double)*p); 00705 double *y=(double *)malloc(sizeof(double)*p); 00706 00707 double *xnew=(double *)malloc(sizeof(double)*p); 00708 double *Ynew=(double *)malloc(sizeof(double)* YSize ); 00709 00710 double *xS=(double *)malloc(sizeof(double)*p); 00711 double *YS=(double *)malloc(sizeof(double)* YSize ); 00712 00713 /*double *xp=(double *)malloc(sizeof(double)*p);*/ 00714 double *Yp=(double *)malloc(sizeof(double)* YSize ); 00715 00716 int *zeroGroupFlag=(int *)malloc(sizeof(int)*g); 00717 int *entrySignFlag=(int *)malloc(sizeof(int)*p); 00718 00719 int pp, gg; 00720 int i, j, iterStep; 00721 double twoNorm,temp, L=1, leftValue, rightValue, gapR, penalty2R; 00722 int nextRestartStep=0; 00723 00724 double alpha, alphap=0.5, beta, gamma; 00725 00726 /* 00727 * call the function to identify some zero entries 00728 * 00729 * entrySignFlag[i]=0 denotes that the corresponding entry is definitely zero 00730 * 00731 * zeroGroupFlag[i]=0 denotes that the corresponding group is definitely zero 00732 * 00733 */ 00734 00735 identifySomeZeroEntries(u, zeroGroupFlag, entrySignFlag, 00736 &pp, &gg, 00737 v, lambda1, lambda2, 00738 p, g, w, G); 00739 00740 penalty2[1]=pp; 00741 penalty2[2]=gg; 00742 /*store pp and gg to penalty2[1] and penalty2[2]*/ 00743 00744 /* 00745 *------------------- 00746 * Process Y 00747 *------------------- 00748 * We make sure that Y is feasible 00749 * and if x_i=0, then set Y_{ij}=0 00750 */ 00751 for(i=0;i<g;i++){ 00752 00753 if(zeroGroupFlag[i]){ /*this group is non-zero*/ 00754 00755 /*compute the two norm of the group*/ 00756 twoNorm=0; 00757 00758 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00759 00760 if (! u[ (int) G[j] ] ) 00761 Y[j]=0; 00762 00763 twoNorm+=Y[j]*Y[j]; 00764 } 00765 twoNorm=sqrt(twoNorm); 00766 00767 if (twoNorm > lambda2 * w[3*i+2] ){ 00768 temp=lambda2 * w[3*i+2] / twoNorm; 00769 00770 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++) 00771 Y[j]*=temp; 00772 } 00773 } 00774 else{ /*this group is zero*/ 00775 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++) 00776 Y[j]=0; 00777 } 00778 } 00779 00780 /* 00781 * set Ynew and Yp to zero 00782 * 00783 * in the following processing, we only operate, Yp, Y and Ynew in the 00784 * possibly non-zero groups by "if(zeroGroupFlag[i])" 00785 * 00786 */ 00787 for(i=0;i<YSize;i++) 00788 YS[i]=Yp[i]=Ynew[i]=0; 00789 00790 00791 /* 00792 * --------------- 00793 * 00794 * we first do a gradient descent step for determing the value of an approporate L 00795 * 00796 * Also, we initialize gamma 00797 * 00798 * with Y, we compute a new Ynew 00799 * 00800 */ 00801 00802 00803 /* 00804 * compute x=max(u-Y * e, 0); 00805 */ 00806 xFromY(x, y, u, Y, p, g, zeroGroupFlag, G, w); 00807 00808 /* 00809 * compute (xnew, Ynew) from (x, Y) 00810 * 00811 * 00812 * gradientDescentStep(double *xnew, double *Ynew, 00813 double *LL, double *u, double *y, int *entrySignFlag, double lambda2, 00814 double *x, double *Y, int p, int g, int * zeroGroupFlag, 00815 double *G, double *w) 00816 */ 00817 00818 gradientDescentStep(xnew, Ynew, 00819 &L, u, y,entrySignFlag,lambda2, 00820 x, Y, p, g, zeroGroupFlag, 00821 G, w); 00822 00823 /* 00824 * we have finished one gradient descent to get 00825 * 00826 * (x, Y) and (xnew, Ynew), where (xnew, Ynew) is 00827 * 00828 * a gradient descent step based on (x, Y) 00829 * 00830 * we set (xp, Yp)=(x, Y) 00831 * 00832 * (x, Y)= (xnew, Ynew) 00833 */ 00834 00835 /*memcpy(xp, x, sizeof(double) * p);*/ 00836 memcpy(Yp, Y, sizeof(double) * YSize); 00837 00838 /*memcpy(x, xnew, sizeof(double) * p);*/ 00839 memcpy(Y, Ynew, sizeof(double) * YSize); 00840 00841 gamma=L; 00842 00843 /* 00844 * ------------------------------------ 00845 * Accelerated Gradient Descent begins here 00846 * ------------------------------------ 00847 */ 00848 00849 00850 for(iterStep=0;iterStep<maxIter;iterStep++){ 00851 00852 00853 while (1){ 00854 00855 00856 /* 00857 * compute alpha as the positive root of 00858 * 00859 * L * alpha^2 = (1-alpha) * gamma 00860 * 00861 */ 00862 00863 alpha= ( - gamma + sqrt( gamma * gamma + 4 * L * gamma ) ) / 2 / L; 00864 00865 beta= gamma * (1-alphap)/ alphap / (gamma + L * alpha); 00866 00867 /* 00868 * compute YS= Y + beta * (Y - Yp) 00869 * 00870 */ 00871 for(i=0;i<g;i++){ 00872 if(zeroGroupFlag[i]){ /*this group is non-zero*/ 00873 00874 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00875 00876 YS[j]=Y[j] + beta * (Y[j]-Yp[j]); 00877 00878 } 00879 } 00880 }/*end of for(i=0;i<g;i++) */ 00881 00882 00883 /* 00884 * compute xS 00885 */ 00886 xFromY(xS, y, u, YS, p, g, zeroGroupFlag, G, w); 00887 00888 00889 /* 00890 * 00891 * Ynew = proj ( YS + xS e^T / L ) 00892 * 00893 */ 00894 for(i=0;i<g;i++){ 00895 if(zeroGroupFlag[i]){ /*this group is non-zero*/ 00896 00897 twoNorm=0; 00898 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00899 00900 Ynew[j]= YS[j] + xS[ (int) G[j] ] / L; 00901 00902 twoNorm+=Ynew[j]*Ynew[j]; 00903 } 00904 twoNorm=sqrt(twoNorm); 00905 00906 if (twoNorm > lambda2 * w[3*i+2] ){ 00907 temp=lambda2 * w[3*i+2] / twoNorm; 00908 00909 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++) 00910 Ynew[j]*=temp; 00911 } 00912 } 00913 }/*end of for(i=0;i<g;i++) */ 00914 00915 /* 00916 * compute xnew=max(u-Ynew * e, 0); 00917 * 00918 *void xFromY(double *x, double *y, 00919 * double *u, double *Y, 00920 * int p, int g, int *zeroGroupFlag, 00921 * double *G, double *w) 00922 */ 00923 00924 xFromY(xnew, y, u, Ynew, p, g, zeroGroupFlag, G, w); 00925 00926 /* test whether L is appropriate*/ 00927 leftValue=0; 00928 for(i=0;i<p;i++){ 00929 if (entrySignFlag[i]){ 00930 temp=xnew[i]-xS[i]; 00931 00932 leftValue+= temp * ( 0.5 * temp + y[i]); 00933 } 00934 } 00935 00936 rightValue=0; 00937 for(i=0;i<g;i++){ 00938 if(zeroGroupFlag[i]){ /*this group is non-zero*/ 00939 00940 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 00941 temp=Ynew[j]-YS[j]; 00942 00943 rightValue+=temp * temp; 00944 } 00945 } 00946 }/*end of for(i=0;i<g;i++) */ 00947 rightValue=rightValue/2; 00948 00949 if ( leftValue <= L * rightValue){ 00950 00951 temp= L * rightValue / leftValue; 00952 00953 if (temp >5) 00954 L=L*0.8; 00955 00956 break; 00957 } 00958 else{ 00959 temp=leftValue / rightValue; 00960 00961 if (L*2 <= temp) 00962 L=temp; 00963 else 00964 L=2*L; 00965 00966 00967 00968 if ( L / g - 2* g >0 ){ 00969 00970 if (rightValue < 1e-16){ 00971 break; 00972 } 00973 else{ 00974 00975 printf("\n AGD: leftValue=%e, rightValue=%e, ratio=%e", leftValue, rightValue, temp); 00976 00977 printf("\n L=%e > 2 * %d * %d. There might be a bug here. Otherwise, it is due to numerical issue.", L, g, g); 00978 00979 break; 00980 } 00981 } 00982 } 00983 } 00984 00985 /* compute the duality gap at (xnew, Ynew) 00986 * 00987 * void dualityGap(double *gap, double *penalty2, 00988 * double *x, double *Y, int g, int *zeroGroupFlag, 00989 * double *G, double *w, double lambda2) 00990 * 00991 */ 00992 dualityGap(gap, penalty2, 00993 xnew, Ynew, g, zeroGroupFlag, 00994 G, w, lambda2); 00995 00996 00997 /* 00998 * if the duality gap is within pre-specified parameter tol 00999 * 01000 * we terminate the algorithm 01001 */ 01002 if (*gap <=tol){ 01003 01004 memcpy(x, xnew, sizeof(double) * p); 01005 memcpy(Y, Ynew, sizeof(double) * YSize); 01006 01007 break; 01008 } 01009 01010 01011 01012 /* 01013 * flag =1 means restart 01014 * 01015 * flag =0 means with restart 01016 * 01017 * nextRestartStep denotes the next "step number" for 01018 * initializing the restart process. 01019 * 01020 * This is based on the fact that, the result is only beneficial when 01021 * xnew is good. In other words, 01022 * if xnew is not good, then the 01023 * restart might not be helpful. 01024 */ 01025 01026 if ( (flag==0) || (flag==1 && iterStep < nextRestartStep )){ 01027 01028 01029 /*memcpy(xp, x, sizeof(double) * p);*/ 01030 memcpy(Yp, Y, sizeof(double) * YSize); 01031 01032 /*memcpy(x, xnew, sizeof(double) * p);*/ 01033 memcpy(Y, Ynew, sizeof(double) * YSize); 01034 01035 gamma=gamma * (1-alpha); 01036 01037 alphap=alpha; 01038 01039 /* 01040 printf("\n iterStep=%d, L=%2.5f, gap=%e", iterStep, L, *gap); 01041 */ 01042 01043 } 01044 else{ 01045 /* 01046 * flag=1 01047 * 01048 * We allow the restart of the program. 01049 * 01050 * Here, Y is constructed as a subgradient of xnew, based on the 01051 * assumption that Y might be a better choice than Ynew, provided 01052 * that xnew is good enough. 01053 * 01054 */ 01055 01056 /* 01057 * compute the restarting point YS with xnew and Ynew 01058 * 01059 *void YFromx(double *Y, 01060 * double *xnew, double *Ynew, 01061 * double lambda2, int g, int *zeroGroupFlag, 01062 * double *G, double *w) 01063 */ 01064 YFromx(YS, xnew, Ynew, lambda2, g, zeroGroupFlag, G, w); 01065 01066 /*compute the solution with the starting point YS 01067 * 01068 *void xFromY(double *x, double *y, 01069 * double *u, double *Y, 01070 * int p, int g, int *zeroGroupFlag, 01071 * double *G, double *w) 01072 * 01073 */ 01074 xFromY(xS, y, u, YS, p, g, zeroGroupFlag, G, w); 01075 01076 /*compute the duality at (xS, YS) 01077 * 01078 * void dualityGap(double *gap, double *penalty2, 01079 * double *x, double *Y, int g, int *zeroGroupFlag, 01080 * double *G, double *w, double lambda2) 01081 * 01082 */ 01083 dualityGap(&gapR, &penalty2R, xS, YS, g, zeroGroupFlag, G, w, lambda2); 01084 01085 if (*gap< gapR){ 01086 /*(xnew, Ynew) is better in terms of duality gap*/ 01087 01088 /*In this case, we do not apply restart, as (xS,YS) is not better 01089 * 01090 * We postpone the "restart" by giving a 01091 * "nextRestartStep" 01092 */ 01093 01094 /*memcpy(xp, x, sizeof(double) * p);*/ 01095 memcpy(Yp, Y, sizeof(double) * YSize); 01096 01097 /*memcpy(x, xnew, sizeof(double) * p);*/ 01098 memcpy(Y, Ynew, sizeof(double) * YSize); 01099 01100 gamma=gamma * (1-alpha); 01101 01102 alphap=alpha; 01103 01104 nextRestartStep=iterStep+ (int) sqrt(gapR / *gap); 01105 } 01106 else{ 01107 /*we use (xS, YS), as it is better in terms of duality gap*/ 01108 01109 *gap=gapR; 01110 *penalty2=penalty2R; 01111 01112 if (*gap <=tol){ 01113 01114 memcpy(x, xS, sizeof(double) * p); 01115 memcpy(Y, YS, sizeof(double) * YSize); 01116 01117 break; 01118 }else{ 01119 /* 01120 * we do a gradient descent based on (xS, YS) 01121 * 01122 */ 01123 01124 /* 01125 * compute (x, Y) from (xS, YS) 01126 * 01127 * 01128 * gradientDescentStep(double *xnew, double *Ynew, 01129 * double *LL, double *u, double *y, int *entrySignFlag, double lambda2, 01130 * double *x, double *Y, int p, int g, int * zeroGroupFlag, 01131 * double *G, double *w) 01132 */ 01133 gradientDescentStep(x, Y, 01134 &L, u, y, entrySignFlag,lambda2, 01135 xS, YS, p, g, zeroGroupFlag, 01136 G, w); 01137 01138 /*memcpy(xp, xS, sizeof(double) * p);*/ 01139 memcpy(Yp, YS, sizeof(double) * YSize); 01140 01141 gamma=L; 01142 01143 alphap=0.5; 01144 01145 } 01146 01147 01148 } 01149 01150 /* 01151 * printf("\n iterStep=%d, L=%2.5f, gap=%e, gapR=%e", iterStep, L, *gap, gapR); 01152 */ 01153 01154 }/* flag =1*/ 01155 01156 } /* main loop */ 01157 01158 01159 01160 penalty2[3]=iterStep+1; 01161 01162 /* 01163 * get the number of nonzero groups 01164 */ 01165 01166 penalty2[4]=0; 01167 for(i=0;i<g;i++){ 01168 if (zeroGroupFlag[i]==0) 01169 penalty2[4]=penalty2[4]+1; 01170 else{ 01171 for(j=(int) w[3*i] ; j<= (int) w[3*i +1]; j++){ 01172 if (x[ (int) G[j] ] !=0) 01173 break; 01174 } 01175 01176 if (j>(int) w[3*i +1]) 01177 penalty2[4]=penalty2[4]+1; 01178 } 01179 } 01180 01181 01182 /* 01183 * assign sign to the solution x 01184 */ 01185 for(i=0;i<p;i++){ 01186 if (entrySignFlag[i]==-1){ 01187 x[i]=-x[i]; 01188 } 01189 } 01190 01191 free (u); 01192 free (y); 01193 01194 free (xnew); 01195 free (Ynew); 01196 01197 free (xS); 01198 free (YS); 01199 01200 /*free (xp);*/ 01201 free (Yp); 01202 01203 free (zeroGroupFlag); 01204 free (entrySignFlag); 01205 } 01206 01207 void overlapping(double *x, double *gap, double *penalty2, 01208 double *v, int p, int g, double lambda1, double lambda2, 01209 double *w, double *G, double *Y, int maxIter, int flag, double tol){ 01210 01211 switch(flag){ 01212 case 0: 01213 case 1: 01214 overlapping_gd(x, gap, penalty2, 01215 v, p, g, lambda1, lambda2, 01216 w, G, Y, maxIter, flag,tol); 01217 break; 01218 case 2: 01219 case 3: 01220 01221 overlapping_agd(x, gap, penalty2, 01222 v, p, g, lambda1, lambda2, 01223 w, G, Y, maxIter, flag-2,tol); 01224 01225 break; 01226 default: 01227 /* printf("\n Wrong flag! The value of flag should be 0,1,2,3. The program uses flag=2.");*/ 01228 01229 overlapping_agd(x, gap, penalty2, 01230 v, p, g, lambda1, lambda2, 01231 w, G, Y, maxIter, 0,tol); 01232 break; 01233 } 01234 01235 01236 }