SHOGUN
v3.2.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 <stdlib.h> 00018 #include <stdio.h> 00019 #include <time.h> 00020 #include <math.h> 00021 #include <shogun/lib/slep/q1/epph.h> 00022 00023 #define delta 1e-8 00024 00025 #define innerIter 1000 00026 #define outerIter 1000 00027 00028 void eplb(double * x, double *root, int * steps, double * v,int n, double z, double lambda0) 00029 { 00030 00031 int i, j, flag=0; 00032 int rho_1, rho_2, rho, rho_T, rho_S; 00033 int V_i_b, V_i_e, V_i; 00034 double lambda_1, lambda_2, lambda_T, lambda_S, lambda; 00035 double s_1, s_2, s, s_T, s_S, v_max, temp; 00036 double f_lambda_1, f_lambda_2, f_lambda, f_lambda_T, f_lambda_S; 00037 int iter_step=0; 00038 00039 /* find the maximal absolute value in v 00040 * and copy the (absolute) values from v to x 00041 */ 00042 00043 if (z< 0){ 00044 printf("\n z should be nonnegative!"); 00045 return; 00046 } 00047 00048 V_i=0; 00049 if (v[0] !=0){ 00050 rho_1=1; 00051 s_1=x[V_i]=v_max=fabs(v[0]); 00052 V_i++; 00053 } 00054 else{ 00055 rho_1=0; 00056 s_1=v_max=0; 00057 } 00058 00059 for (i=1;i<n; i++){ 00060 if (v[i]!=0){ 00061 x[V_i]=fabs(v[i]); s_1+= x[V_i]; rho_1++; 00062 00063 if (x[V_i] > v_max) 00064 v_max=x[V_i]; 00065 V_i++; 00066 } 00067 } 00068 00069 /* If ||v||_1 <= z, then v is the solution */ 00070 if (s_1 <= z){ 00071 flag=1; lambda=0; 00072 for(i=0;i<n;i++){ 00073 x[i]=v[i]; 00074 } 00075 *root=lambda; 00076 *steps=iter_step; 00077 return; 00078 } 00079 00080 lambda_1=0; lambda_2=v_max; 00081 f_lambda_1=s_1 -z; 00082 /*f_lambda_1=s_1-rho_1* lambda_1 -z;*/ 00083 rho_2=0; s_2=0; f_lambda_2=-z; 00084 V_i_b=0; V_i_e=V_i-1; 00085 00086 lambda=lambda0; 00087 if ( (lambda<lambda_2) && (lambda> lambda_1) ){ 00088 /*------------------------------------------------------------------- 00089 Initialization with the root 00090 *------------------------------------------------------------------- 00091 */ 00092 00093 i=V_i_b; j=V_i_e; rho=0; s=0; 00094 while (i <= j){ 00095 while( (i <= V_i_e) && (x[i] <= lambda) ){ 00096 i++; 00097 } 00098 while( (j>=V_i_b) && (x[j] > lambda) ){ 00099 s+=x[j]; 00100 j--; 00101 } 00102 if (i<j){ 00103 s+=x[i]; 00104 00105 temp=x[i]; x[i]=x[j]; x[j]=temp; 00106 i++; j--; 00107 } 00108 } 00109 00110 rho=V_i_e-j; rho+=rho_2; s+=s_2; 00111 f_lambda=s-rho*lambda-z; 00112 00113 if ( fabs(f_lambda)< delta ){ 00114 flag=1; 00115 } 00116 00117 if (f_lambda <0){ 00118 lambda_2=lambda; s_2=s; rho_2=rho; f_lambda_2=f_lambda; 00119 00120 V_i_e=j; V_i=V_i_e-V_i_b+1; 00121 } 00122 else{ 00123 lambda_1=lambda; rho_1=rho; s_1=s; f_lambda_1=f_lambda; 00124 00125 V_i_b=i; V_i=V_i_e-V_i_b+1; 00126 } 00127 00128 if (V_i==0){ 00129 /*printf("\n rho=%d, rho_1=%d, rho_2=%d",rho, rho_1, rho_2); 00130 00131 printf("\n V_i=%d",V_i);*/ 00132 00133 lambda=(s - z)/ rho; 00134 flag=1; 00135 } 00136 /*------------------------------------------------------------------- 00137 End of initialization 00138 *-------------------------------------------------------------------- 00139 */ 00140 00141 }/* end of if(!flag) */ 00142 00143 while (!flag){ 00144 iter_step++; 00145 00146 /* compute lambda_T */ 00147 lambda_T=lambda_1 + f_lambda_1 /rho_1; 00148 if(rho_2 !=0){ 00149 if (lambda_2 + f_lambda_2 /rho_2 > lambda_T) 00150 lambda_T=lambda_2 + f_lambda_2 /rho_2; 00151 } 00152 00153 /* compute lambda_S */ 00154 lambda_S=lambda_2 - f_lambda_2 *(lambda_2-lambda_1)/(f_lambda_2-f_lambda_1); 00155 00156 if (fabs(lambda_T-lambda_S) <= delta){ 00157 lambda=lambda_T; flag=1; 00158 break; 00159 } 00160 00161 /* set lambda as the middle point of lambda_T and lambda_S */ 00162 lambda=(lambda_T+lambda_S)/2; 00163 00164 s_T=s_S=s=0; 00165 rho_T=rho_S=rho=0; 00166 i=V_i_b; j=V_i_e; 00167 while (i <= j){ 00168 while( (i <= V_i_e) && (x[i] <= lambda) ){ 00169 if (x[i]> lambda_T){ 00170 s_T+=x[i]; rho_T++; 00171 } 00172 i++; 00173 } 00174 while( (j>=V_i_b) && (x[j] > lambda) ){ 00175 if (x[j] > lambda_S){ 00176 s_S+=x[j]; rho_S++; 00177 } 00178 else{ 00179 s+=x[j]; rho++; 00180 } 00181 j--; 00182 } 00183 if (i<j){ 00184 if (x[i] > lambda_S){ 00185 s_S+=x[i]; rho_S++; 00186 } 00187 else{ 00188 s+=x[i]; rho++; 00189 } 00190 00191 if (x[j]> lambda_T){ 00192 s_T+=x[j]; rho_T++; 00193 } 00194 00195 temp=x[i]; x[i]=x[j]; x[j]=temp; 00196 i++; j--; 00197 } 00198 } 00199 00200 s_S+=s_2; rho_S+=rho_2; 00201 s+=s_S; rho+=rho_S; 00202 s_T+=s; rho_T+=rho; 00203 f_lambda_S=s_S-rho_S*lambda_S-z; 00204 f_lambda=s-rho*lambda-z; 00205 f_lambda_T=s_T-rho_T*lambda_T-z; 00206 00207 /*printf("\n %d & %d & %5.6f & %5.6f & %5.6f & %5.6f & %5.6f \\\\ \n \\hline ", iter_step, V_i, lambda_1, lambda_T, lambda, lambda_S, lambda_2);*/ 00208 00209 if ( fabs(f_lambda)< delta ){ 00210 /*printf("\n lambda");*/ 00211 flag=1; 00212 break; 00213 } 00214 if ( fabs(f_lambda_S)< delta ){ 00215 /* printf("\n lambda_S");*/ 00216 lambda=lambda_S; flag=1; 00217 break; 00218 } 00219 if ( fabs(f_lambda_T)< delta ){ 00220 /* printf("\n lambda_T");*/ 00221 lambda=lambda_T; flag=1; 00222 break; 00223 } 00224 00225 /* 00226 printf("\n\n f_lambda_1=%5.6f, f_lambda_2=%5.6f, f_lambda=%5.6f",f_lambda_1,f_lambda_2, f_lambda); 00227 printf("\n lambda_1=%5.6f, lambda_2=%5.6f, lambda=%5.6f",lambda_1, lambda_2, lambda); 00228 printf("\n rho_1=%d, rho_2=%d, rho=%d ",rho_1, rho_2, rho); 00229 */ 00230 00231 if (f_lambda <0){ 00232 lambda_2=lambda; s_2=s; rho_2=rho; 00233 f_lambda_2=f_lambda; 00234 00235 lambda_1=lambda_T; s_1=s_T; rho_1=rho_T; 00236 f_lambda_1=f_lambda_T; 00237 00238 V_i_e=j; i=V_i_b; 00239 while (i <= j){ 00240 while( (i <= V_i_e) && (x[i] <= lambda_T) ){ 00241 i++; 00242 } 00243 while( (j>=V_i_b) && (x[j] > lambda_T) ){ 00244 j--; 00245 } 00246 if (i<j){ 00247 x[j]=x[i]; 00248 i++; j--; 00249 } 00250 } 00251 V_i_b=i; V_i=V_i_e-V_i_b+1; 00252 } 00253 else{ 00254 lambda_1=lambda; s_1=s; rho_1=rho; 00255 f_lambda_1=f_lambda; 00256 00257 lambda_2=lambda_S; s_2=s_S; rho_2=rho_S; 00258 f_lambda_2=f_lambda_S; 00259 00260 V_i_b=i; j=V_i_e; 00261 while (i <= j){ 00262 while( (i <= V_i_e) && (x[i] <= lambda_S) ){ 00263 i++; 00264 } 00265 while( (j>=V_i_b) && (x[j] > lambda_S) ){ 00266 j--; 00267 } 00268 if (i<j){ 00269 x[i]=x[j]; 00270 i++; j--; 00271 } 00272 } 00273 V_i_e=j; V_i=V_i_e-V_i_b+1; 00274 } 00275 00276 if (V_i==0){ 00277 lambda=(s - z)/ rho; flag=1; 00278 /*printf("\n V_i=0, lambda=%5.6f",lambda);*/ 00279 break; 00280 } 00281 }/* end of while */ 00282 00283 00284 for(i=0;i<n;i++){ 00285 if (v[i] > lambda) 00286 x[i]=v[i]-lambda; 00287 else 00288 if (v[i]< -lambda) 00289 x[i]=v[i]+lambda; 00290 else 00291 x[i]=0; 00292 } 00293 *root=lambda; 00294 *steps=iter_step; 00295 } 00296 00297 void epp1(double *x, double *v, int n, double rho) 00298 { 00299 int i; 00300 00301 /* 00302 we assume rho>=0 00303 */ 00304 00305 for(i=0;i<n;i++){ 00306 if (fabs(v[i])<=rho) 00307 x[i]=0; 00308 else 00309 if (v[i]< -rho) 00310 x[i]=v[i]+rho; 00311 else 00312 x[i]=v[i]-rho; 00313 } 00314 } 00315 00316 void epp2(double *x, double *v, int n, double rho) 00317 { 00318 int i; 00319 double v2=0, ratio; 00320 00321 /* 00322 we assume rho>=0 00323 */ 00324 00325 for(i=0; i< n; i++){ 00326 v2+=v[i]*v[i]; 00327 } 00328 v2=sqrt(v2); 00329 00330 if (rho >= v2) 00331 for(i=0;i<n;i++) 00332 x[i]=0; 00333 else{ 00334 ratio= (v2-rho) /v2; 00335 for(i=0;i<n;i++) 00336 x[i]=v[i]*ratio; 00337 } 00338 } 00339 00340 void eppInf(double *x, double * c, int * iter_step, double *v, int n, double rho, double c0) 00341 { 00342 int i, steps; 00343 00344 /* 00345 we assume rho>=0 00346 */ 00347 00348 eplb(x, c, &steps, v, n, rho, c0); 00349 00350 for(i=0; i< n; i++){ 00351 x[i]=v[i]-x[i]; 00352 } 00353 iter_step[0]=steps; 00354 iter_step[1]=0; 00355 } 00356 00357 void zerofind(double *root, int * iterStep, double v, double p, double c, double x0) 00358 { 00359 double x, f, fprime, p1=p-1, pp; 00360 int step=0; 00361 00362 00363 if (v==0){ 00364 *root=0; *iterStep=0; return; 00365 } 00366 00367 if (c==0){ 00368 *root=v; * iterStep=0; return; 00369 } 00370 00371 00372 if ( (x0 <v) && (x0>0) ) 00373 x=x0; 00374 else 00375 x=v; 00376 00377 00378 pp=pow(x, p1); 00379 f= x + c* pp -v; 00380 00381 00382 /* 00383 We apply the Newton's method for solving the root 00384 */ 00385 while (1){ 00386 step++; 00387 00388 fprime=1 + c* p1 * pp / x; 00389 /* 00390 The derivative at the current solution x 00391 */ 00392 00393 x = x- f/fprime; /* 00394 The new solution is computed by the Newton method 00395 */ 00396 00397 00398 00399 if (p>2){ 00400 if (x>v){ 00401 x=v; 00402 } 00403 } 00404 else{ 00405 if ( (x<0) || (x>v)){ 00406 x=1e-30; 00407 00408 f= x+c* pow(x,p1)-v; 00409 00410 if (f>0){ /* 00411 If f(x) = x + c x^{p-1} - v <0 at x=1e-30, 00412 this shows that the real root is between (0, 1e-30). 00413 For numerical reason, we just set x=0 00414 */ 00415 00416 *root=x; 00417 * iterStep=step; 00418 00419 break; 00420 } 00421 } 00422 } 00423 /* 00424 This makes sure that x lies in the interval [0, v] 00425 */ 00426 00427 pp=pow(x, p1); 00428 f= x + c* pp -v; 00429 /* 00430 The function value at the new solution 00431 */ 00432 00433 if ( fabs(f) <= delta){ 00434 *root=x; 00435 * iterStep=step; 00436 break; 00437 } 00438 00439 if (step>=innerIter){ 00440 printf("\n The number of steps exceed %d, in finding the root for f(x)= x + c x^{p-1} - v, 0< x< v.", innerIter); 00441 printf("\n If you meet with this problem, please contact Jun Liu (j.liu@asu.edu). Thanks!"); 00442 return; 00443 } 00444 00445 } 00446 00447 /* 00448 printf("\n x=%e, f=%e, step=%d\n",x, f, step); 00449 */ 00450 } 00451 00452 double norm(double * v, double p, int n) 00453 { 00454 int i; 00455 double t=0; 00456 00457 00458 /* 00459 we assume that v[i]>=0 00460 p>1 00461 */ 00462 00463 for(i=0;i<n;i++) 00464 t+=pow(v[i], p); 00465 00466 return( pow(t, 1/p) ); 00467 }; 00468 00469 void eppO(double *x, double * cc, int * iter_step, double *v, int n, double rho, double p) 00470 { 00471 int i, *flag, bisStep, newtonStep=0, totoalStep=0; 00472 double vq=0, epsilon, vmax=0, vmin=1e10; /* we assume that the minimal value in |v| is less than 1e10*/ 00473 double q=1/(1-1/p), c, c1, c2, root, f, xp; 00474 00475 double x_diff=0; /* this value denotes the maximal difference of the x values computed from c1 and c2*/ 00476 double temp; 00477 int p_n=1; /* p_n indicates the previous phi(c) is positive or negative*/ 00478 00479 flag=(int *)malloc(sizeof(int)*n); 00480 00481 /* 00482 compute vq, the q-norm of v 00483 flag denotes the sign of v: 00484 flag[i]=0 denotes v[i] is non-negative 00485 flag[i]=1 denotes v[i] is negative 00486 vmin and vmax are the maximal and minimal value of |v| (excluding 0) 00487 */ 00488 for(i=0; i< n; i++){ 00489 00490 x[i]=0; 00491 00492 if (v[i]==0) 00493 flag[i]=0; 00494 else 00495 { 00496 if (v[i]>0) 00497 flag[i]=0; 00498 else 00499 { 00500 flag[i]=1; 00501 v[i]=-v[i];/* 00502 we set v[i] to its absolute value 00503 */ 00504 } 00505 00506 vq+=pow(v[i], q); 00507 00508 00509 if (v[i]>vmax) 00510 vmax=v[i]; 00511 00512 if (v[i]<vmin) 00513 vmin=v[i]; 00514 } 00515 } 00516 vq=pow(vq, 1/q); 00517 00518 /* 00519 zero solution 00520 */ 00521 if (rho >= vq){ 00522 *cc=0; 00523 iter_step[0]=iter_step[1]=0; 00524 00525 00526 for(i=0;i<n;i++){ 00527 if (flag[i]) 00528 v[i]=-v[i]; /* set the value of v[i] back*/ 00529 } 00530 00531 free(flag); 00532 return; 00533 } 00534 00535 /* 00536 compute epsilon 00537 initialize c1 and c2, the interval where the root lies 00538 */ 00539 epsilon=(vq -rho)/ vq; 00540 if (p>2){ 00541 00542 if ( log((1-epsilon) * vmin) - (p-1) * log( epsilon* vmin ) >= 709 ) 00543 { 00544 /* If this contition holds, we have c2 >= 1e308, exceeding the machine precision. 00545 00546 In this case, the reason is that p is too large 00547 and meanwhile epsilon * vmin is typically small. 00548 00549 For numerical stablity, we just regard p=inf, and run eppInf 00550 */ 00551 00552 00553 for(i=0;i<n;i++){ 00554 if (flag[i]) 00555 v[i]=-v[i]; /* set the value of v[i] back*/ 00556 } 00557 00558 eppInf(x, cc, iter_step, v, n, rho, 0); 00559 00560 free(flag); 00561 return; 00562 } 00563 00564 c1= (1-epsilon) * vmax / pow(epsilon* vmax, p-1); 00565 c2= (1-epsilon) * vmin / pow(epsilon* vmin, p-1); 00566 } 00567 else{ /*1 < p < 2*/ 00568 00569 c2= (1-epsilon) * vmax / pow(epsilon* vmax, p-1); 00570 c1= (1-epsilon) * vmin / pow(epsilon* vmin, p-1); 00571 } 00572 00573 00574 /* 00575 printf("\n c1=%e, c2=%e", c1, c2); 00576 */ 00577 00578 if (fabs(c1-c2) <= delta){ 00579 c=c1; 00580 } 00581 else 00582 c=(c1+c2)/2; 00583 00584 00585 bisStep =0; 00586 00587 while(1){ 00588 bisStep++; 00589 00590 /*compute the root corresponding to c*/ 00591 x_diff=0; 00592 for(i=0;i<n;i++){ 00593 zerofind(&root, &newtonStep, v[i], p, c, x[i]); 00594 00595 temp=fabs(root-x[i]); 00596 if (x_diff< temp ) 00597 x_diff=temp; /*x_diff denotes the largest gap to the previous solution*/ 00598 00599 x[i]=root; 00600 totoalStep+=newtonStep; 00601 } 00602 00603 xp=norm(x, p, n); 00604 00605 f=rho * pow(xp, 1-p) - c; 00606 00607 if ( fabs(f)<=delta || fabs(c1-c2)<=delta ) 00608 break; 00609 else{ 00610 if (f>0){ 00611 if ( (x_diff <=delta) && (p_n==0) ) 00612 break; 00613 00614 c1=c; p_n=1; 00615 } 00616 else{ 00617 00618 if ( (x_diff <=delta) && (p_n==1) ) 00619 break; 00620 00621 c2=c; p_n=0; 00622 } 00623 } 00624 c=(c1+c2)/2; 00625 00626 if (bisStep>=outerIter){ 00627 00628 00629 if ( fabs(c1-c2) <=delta * c2 ) 00630 break; 00631 else{ 00632 printf("\n The number of bisection steps exceed %d.", outerIter); 00633 printf("\n c1=%e, c2=%e, x_diff=%e, f=%e",c1,c2,x_diff,f); 00634 printf("\n If you meet with this problem, please contact Jun Liu (j.liu@asu.edu). Thanks!"); 00635 00636 return; 00637 } 00638 } 00639 00640 /* 00641 printf("\n c1=%e, c2=%e, f=%e, newtonStep=%d", c1, c2, f, newtonStep); 00642 */ 00643 } 00644 00645 /* 00646 printf("\n c1=%e, c2=%e, x_diff=%e, f=%e, bisStep=%d, totoalStep=%d",c1,c2, x_diff, f,bisStep,totoalStep); 00647 */ 00648 00649 for(i=0;i<n;i++){ 00650 if (flag[i]){ 00651 x[i]=-x[i]; 00652 v[i]=-v[i]; 00653 } 00654 } 00655 free(flag); 00656 00657 *cc=c; 00658 00659 iter_step[0]=bisStep; 00660 iter_step[1]=totoalStep; 00661 } 00662 00663 void epp(double *x, double * c, int * iter_step, double * v, int n, double rho, double p, double c0){ 00664 if (rho <0){ 00665 printf("\n rho should be non-negative!"); 00666 exit(1); 00667 } 00668 00669 if (p==1){ 00670 epp1(x, v, n, rho); 00671 *c=0; 00672 iter_step[0]=iter_step[1]=0; 00673 } 00674 else 00675 if (p==2){ 00676 epp2(x, v, n, rho); 00677 *c=0; 00678 iter_step[0]=iter_step[1]=0; 00679 } 00680 else 00681 if (p>=1e6) /* when p >=1e6, we treat it as infity*/ 00682 eppInf(x, c, iter_step, v, n, rho, c0); 00683 else 00684 eppO(x, c, iter_step, v, n, rho, p); 00685 } 00686