24 void identifySomeZeroEntries(
double * u,
int * zeroGroupFlag,
int *entrySignFlag,
26 double *v,
double lambda1,
double lambda2,
27 int p,
int g,
double * w,
double *G){
29 int i, j, newZeroNum, iterStep=0;
73 printf(
"\n Identify Zero Group: iterStep= %d. The code might have a bug! Check it!", iterStep);
87 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
91 twoNorm=sqrt(twoNorm);
100 if (twoNorm<= lambda2 * w[3*i+2] ){
103 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++)
131 if (zeroGroupFlag[i]==0)
136 void xFromY(
double *x,
double *y,
137 double *u,
double *Y,
138 int p,
int g,
int *zeroGroupFlag,
139 double *G,
double *w){
148 if(zeroGroupFlag[i]){
150 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
151 x[ (int) G[j] ] -= Y[j];
167 void YFromx(
double *Y,
168 double *xnew,
double *Ynew,
169 double lambda2,
int g,
int *zeroGroupFlag,
170 double *G,
double *w){
173 double twoNorm, temp;
176 if(zeroGroupFlag[i]){
179 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
180 temp=xnew[ (int) G[j] ];
186 twoNorm=sqrt(twoNorm);
189 temp=lambda2 * w[3*i+2] / twoNorm;
191 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++)
196 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++)
203 void dualityGap(
double *gap,
double *penalty2,
204 double *x,
double *Y,
int g,
int *zeroGroupFlag,
205 double *G,
double *w,
double lambda2){
208 double temp, twoNorm, innerProduct;
213 if(zeroGroupFlag[i]){
215 twoNorm=0;innerProduct=0;
217 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
218 temp=x[ (int) G[j] ];
222 innerProduct+=temp * Y[j];
225 twoNorm=sqrt(twoNorm)* w[3*i +2];
229 twoNorm=lambda2 * twoNorm;
230 if (twoNorm > innerProduct)
231 *gap+=twoNorm-innerProduct;
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){
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);
244 double *xnew=(
double *)malloc(
sizeof(
double)*p);
245 double *Ynew=(
double *)malloc(
sizeof(
double)* YSize );
247 int *zeroGroupFlag=(
int *)malloc(
sizeof(
int)*g);
248 int *entrySignFlag=(
int *)malloc(
sizeof(
int)*p);
251 double twoNorm,temp, L=1, leftValue, rightValue, gapR, penalty2R;
252 int nextRestartStep=0;
263 identifySomeZeroEntries(u, zeroGroupFlag, entrySignFlag,
282 if(zeroGroupFlag[i]){
287 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
289 if (! u[ (
int) G[j] ] )
294 twoNorm=sqrt(twoNorm);
296 if (twoNorm > lambda2 * w[3*i+2] ){
297 temp=lambda2 * w[3*i+2] / twoNorm;
299 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++)
304 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++)
329 xFromY(x, y, u, Y, p, g, zeroGroupFlag, G, w);
334 for(iterStep=0;iterStep<maxIter;iterStep++){
357 if(zeroGroupFlag[i]){
360 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
361 Ynew[j]= Y[j] + x[ (int) G[j] ] / L;
363 twoNorm+=Ynew[j]*Ynew[j];
365 twoNorm=sqrt(twoNorm);
367 if (twoNorm > lambda2 * w[3*i+2] ){
368 temp=lambda2 * w[3*i+2] / twoNorm;
370 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++)
384 xFromY(xnew, y, u, Ynew, p, g, zeroGroupFlag, G, w);
389 if (entrySignFlag[i]){
392 leftValue+= temp * ( 0.5 * temp + y[i]);
398 if(zeroGroupFlag[i]){
400 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
403 rightValue+=temp * temp;
407 rightValue=rightValue/2;
409 if ( leftValue <= L * rightValue){
411 temp= L * rightValue / leftValue;
419 temp=leftValue / rightValue;
429 if (rightValue < 1e-16){
434 printf(
"\n GD: leftValue=%e, rightValue=%e, ratio=%e", leftValue, rightValue, temp);
436 printf(
"\n L=%e > 2 * %d * %d. There might be a bug here. Otherwise, it is due to numerical issue.", L, g, g);
451 dualityGap(gap, penalty2, xnew, Ynew, g, zeroGroupFlag, G, w, lambda2);
467 if ( (flag==0) || (flag==1 && iterStep < nextRestartStep )){
470 memcpy(x, xnew,
sizeof(
double) * p);
471 memcpy(Y, Ynew,
sizeof(
double) * YSize);
498 YFromx(Y, xnew, Ynew, lambda2, g, zeroGroupFlag, G, w);
508 xFromY(x, y, u, Y, p, g, zeroGroupFlag, G, w);
517 dualityGap(&gapR, &penalty2R, x, Y, g, zeroGroupFlag, G, w, lambda2);
522 memcpy(x, xnew,
sizeof(
double) * p);
523 memcpy(Y, Ynew,
sizeof(
double) * YSize);
537 nextRestartStep=iterStep+ (int) sqrt(gapR / *gap);
560 penalty2[3]=iterStep;
564 if (zeroGroupFlag[i]==0)
565 penalty2[4]=penalty2[4]+1;
567 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
568 if (x[ (
int) G[j] ] !=0)
572 if (j>(
int) w[3*i +1])
573 penalty2[4]=penalty2[4]+1;
581 if (entrySignFlag[i]==-1){
590 free (zeroGroupFlag);
591 free (entrySignFlag);
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){
599 double twoNorm, temp, L=*LL, leftValue, rightValue;
611 if(zeroGroupFlag[i]){
614 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
615 Ynew[j]= Y[j] + x[ (int) G[j] ] / L;
617 twoNorm+=Ynew[j]*Ynew[j];
619 twoNorm=sqrt(twoNorm);
621 if (twoNorm > lambda2 * w[3*i+2] ){
622 temp=lambda2 * w[3*i+2] / twoNorm;
624 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++)
638 xFromY(xnew, y, u, Ynew, p, g, zeroGroupFlag, G, w);
643 if (entrySignFlag[i]){
646 leftValue+= temp * ( 0.5 * temp + y[i]);
652 if(zeroGroupFlag[i]){
654 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
657 rightValue+=temp * temp;
661 rightValue=rightValue/2;
667 if ( leftValue <= L * rightValue){
669 temp= L * rightValue / leftValue;
677 temp=leftValue / rightValue;
684 if ( L / g - 2* g >0 ){
686 if (rightValue < 1e-16){
691 printf(
"\n One Gradient Step: leftValue=%e, rightValue=%e, ratio=%e", leftValue, rightValue, temp);
693 printf(
"\n L=%e > 2 * %d * %d. There might be a bug here. Otherwise, it is due to numerical issue.", L, g, g);
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){
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);
712 double *xnew=(
double *)malloc(
sizeof(
double)*p);
713 double *Ynew=(
double *)malloc(
sizeof(
double)* YSize );
715 double *xS=(
double *)malloc(
sizeof(
double)*p);
716 double *YS=(
double *)malloc(
sizeof(
double)* YSize );
719 double *Yp=(
double *)malloc(
sizeof(
double)* YSize );
721 int *zeroGroupFlag=(
int *)malloc(
sizeof(
int)*g);
722 int *entrySignFlag=(
int *)malloc(
sizeof(
int)*p);
726 double twoNorm,temp, L=1, leftValue, rightValue, gapR, penalty2R;
727 int nextRestartStep=0;
729 double alpha, alphap=0.5, beta, gamma;
740 identifySomeZeroEntries(u, zeroGroupFlag, entrySignFlag,
758 if(zeroGroupFlag[i]){
763 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
765 if (! u[ (
int) G[j] ] )
770 twoNorm=sqrt(twoNorm);
772 if (twoNorm > lambda2 * w[3*i+2] ){
773 temp=lambda2 * w[3*i+2] / twoNorm;
775 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++)
780 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++)
793 YS[i]=Yp[i]=Ynew[i]=0;
811 xFromY(x, y, u, Y, p, g, zeroGroupFlag, G, w);
823 gradientDescentStep(xnew, Ynew,
824 &L, u, y,entrySignFlag,lambda2,
825 x, Y, p, g, zeroGroupFlag,
841 memcpy(Yp, Y,
sizeof(
double) * YSize);
844 memcpy(Y, Ynew,
sizeof(
double) * YSize);
855 for(iterStep=0;iterStep<maxIter;iterStep++){
868 alpha= ( - gamma + sqrt( gamma * gamma + 4 * L * gamma ) ) / 2 / L;
870 beta= gamma * (1-alphap)/ alphap / (gamma + L * alpha);
877 if(zeroGroupFlag[i]){
879 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
881 YS[j]=Y[j] + beta * (Y[j]-Yp[j]);
891 xFromY(xS, y, u, YS, p, g, zeroGroupFlag, G, w);
900 if(zeroGroupFlag[i]){
903 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
905 Ynew[j]= YS[j] + xS[ (int) G[j] ] / L;
907 twoNorm+=Ynew[j]*Ynew[j];
909 twoNorm=sqrt(twoNorm);
911 if (twoNorm > lambda2 * w[3*i+2] ){
912 temp=lambda2 * w[3*i+2] / twoNorm;
914 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++)
929 xFromY(xnew, y, u, Ynew, p, g, zeroGroupFlag, G, w);
934 if (entrySignFlag[i]){
937 leftValue+= temp * ( 0.5 * temp + y[i]);
943 if(zeroGroupFlag[i]){
945 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
948 rightValue+=temp * temp;
952 rightValue=rightValue/2;
954 if ( leftValue <= L * rightValue){
956 temp= L * rightValue / leftValue;
964 temp=leftValue / rightValue;
973 if ( L / g - 2* g >0 ){
975 if (rightValue < 1e-16){
980 printf(
"\n AGD: leftValue=%e, rightValue=%e, ratio=%e", leftValue, rightValue, temp);
982 printf(
"\n L=%e > 2 * %d * %d. There might be a bug here. Otherwise, it is due to numerical issue.", L, g, g);
997 dualityGap(gap, penalty2,
998 xnew, Ynew, g, zeroGroupFlag,
1009 memcpy(x, xnew,
sizeof(
double) * p);
1010 memcpy(Y, Ynew,
sizeof(
double) * YSize);
1031 if ( (flag==0) || (flag==1 && iterStep < nextRestartStep )){
1035 memcpy(Yp, Y,
sizeof(
double) * YSize);
1038 memcpy(Y, Ynew,
sizeof(
double) * YSize);
1040 gamma=gamma * (1-alpha);
1069 YFromx(YS, xnew, Ynew, lambda2, g, zeroGroupFlag, G, w);
1079 xFromY(xS, y, u, YS, p, g, zeroGroupFlag, G, w);
1088 dualityGap(&gapR, &penalty2R, xS, YS, g, zeroGroupFlag, G, w, lambda2);
1100 memcpy(Yp, Y,
sizeof(
double) * YSize);
1103 memcpy(Y, Ynew,
sizeof(
double) * YSize);
1105 gamma=gamma * (1-alpha);
1109 nextRestartStep=iterStep+ (int) sqrt(gapR / *gap);
1115 *penalty2=penalty2R;
1119 memcpy(x, xS,
sizeof(
double) * p);
1120 memcpy(Y, YS,
sizeof(
double) * YSize);
1138 gradientDescentStep(x, Y,
1139 &L, u, y, entrySignFlag,lambda2,
1140 xS, YS, p, g, zeroGroupFlag,
1144 memcpy(Yp, YS,
sizeof(
double) * YSize);
1165 penalty2[3]=iterStep+1;
1173 if (zeroGroupFlag[i]==0)
1174 penalty2[4]=penalty2[4]+1;
1176 for(j=(
int) w[3*i] ; j<= (int) w[3*i +1]; j++){
1177 if (x[ (
int) G[j] ] !=0)
1181 if (j>(
int) w[3*i +1])
1182 penalty2[4]=penalty2[4]+1;
1191 if (entrySignFlag[i]==-1){
1208 free (zeroGroupFlag);
1209 free (entrySignFlag);
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){
1219 overlapping_gd(x, gap, penalty2,
1220 v, p, g, lambda1, lambda2,
1221 w, G, Y, maxIter, flag,tol);
1226 overlapping_agd(x, gap, penalty2,
1227 v, p, g, lambda1, lambda2,
1228 w, G, Y, maxIter, flag-2,tol);
1234 overlapping_agd(x, gap, penalty2,
1235 v, p, g, lambda1, lambda2,
1236 w, G, Y, maxIter, 0,tol);
1242 #endif //USE_GPL_SHOGUN