SHOGUN  6.0.0
epph.cpp
Go to the documentation of this file.
1 /* This program is free software: you can redistribute it and/or modify
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
18 #ifdef USE_GPL_SHOGUN
19
20 #include <stdlib.h>
21 #include <stdio.h>
22 #include <time.h>
23 #include <math.h>
24
25 #define delta 1e-8
26
27 #define innerIter 1000
28 #define outerIter 1000
29
30 void eplb(double * x, double *root, int * steps, double * v,int n, double z, double lambda0)
31 {
32
33  int i, j, flag=0;
34  int rho_1, rho_2, rho, rho_T, rho_S;
35  int V_i_b, V_i_e, V_i;
36  double lambda_1, lambda_2, lambda_T, lambda_S, lambda;
37  double s_1, s_2, s, s_T, s_S, v_max, temp;
38  double f_lambda_1, f_lambda_2, f_lambda, f_lambda_T, f_lambda_S;
39  int iter_step=0;
40
41  /* find the maximal absolute value in v
42  * and copy the (absolute) values from v to x
43  */
44
45  if (z< 0){
46  printf("\n z should be nonnegative!");
47  return;
48  }
49
50  V_i=0;
51  if (v[0] !=0){
52  rho_1=1;
53  s_1=x[V_i]=v_max=fabs(v[0]);
54  V_i++;
55  }
56  else{
57  rho_1=0;
58  s_1=v_max=0;
59  }
60
61  for (i=1;i<n; i++){
62  if (v[i]!=0){
63  x[V_i]=fabs(v[i]); s_1+= x[V_i]; rho_1++;
64
65  if (x[V_i] > v_max)
66  v_max=x[V_i];
67  V_i++;
68  }
69  }
70
71  /* If ||v||_1 <= z, then v is the solution */
72  if (s_1 <= z){
73  flag=1; lambda=0;
74  for(i=0;i<n;i++){
75  x[i]=v[i];
76  }
77  *root=lambda;
78  *steps=iter_step;
79  return;
80  }
81
82  lambda_1=0; lambda_2=v_max;
83  f_lambda_1=s_1 -z;
84  /*f_lambda_1=s_1-rho_1* lambda_1 -z;*/
85  rho_2=0; s_2=0; f_lambda_2=-z;
86  V_i_b=0; V_i_e=V_i-1;
87
88  lambda=lambda0;
89  if ( (lambda<lambda_2) && (lambda> lambda_1) ){
90  /*-------------------------------------------------------------------
91  Initialization with the root
92  *-------------------------------------------------------------------
93  */
94
95  i=V_i_b; j=V_i_e; rho=0; s=0;
96  while (i <= j){
97  while( (i <= V_i_e) && (x[i] <= lambda) ){
98  i++;
99  }
100  while( (j>=V_i_b) && (x[j] > lambda) ){
101  s+=x[j];
102  j--;
103  }
104  if (i<j){
105  s+=x[i];
106
107  temp=x[i]; x[i]=x[j]; x[j]=temp;
108  i++; j--;
109  }
110  }
111
112  rho=V_i_e-j; rho+=rho_2; s+=s_2;
113  f_lambda=s-rho*lambda-z;
114
115  if ( fabs(f_lambda)< delta ){
116  flag=1;
117  }
118
119  if (f_lambda <0){
120  lambda_2=lambda; s_2=s; rho_2=rho; f_lambda_2=f_lambda;
121
122  V_i_e=j; V_i=V_i_e-V_i_b+1;
123  }
124  else{
125  lambda_1=lambda; rho_1=rho; s_1=s; f_lambda_1=f_lambda;
126
127  V_i_b=i; V_i=V_i_e-V_i_b+1;
128  }
129
130  if (V_i==0){
131  /*printf("\n rho=%d, rho_1=%d, rho_2=%d",rho, rho_1, rho_2);
132
133  printf("\n V_i=%d",V_i);*/
134
135  lambda=(s - z)/ rho;
136  flag=1;
137  }
138  /*-------------------------------------------------------------------
139  End of initialization
140  *--------------------------------------------------------------------
141  */
142
143  }/* end of if(!flag) */
144
145  while (!flag){
146  iter_step++;
147
148  /* compute lambda_T */
149  lambda_T=lambda_1 + f_lambda_1 /rho_1;
150  if(rho_2 !=0){
151  if (lambda_2 + f_lambda_2 /rho_2 > lambda_T)
152  lambda_T=lambda_2 + f_lambda_2 /rho_2;
153  }
154
155  /* compute lambda_S */
156  lambda_S=lambda_2 - f_lambda_2 *(lambda_2-lambda_1)/(f_lambda_2-f_lambda_1);
157
158  if (fabs(lambda_T-lambda_S) <= delta){
159  lambda=lambda_T; flag=1;
160  break;
161  }
162
163  /* set lambda as the middle point of lambda_T and lambda_S */
164  lambda=(lambda_T+lambda_S)/2;
165
166  s_T=s_S=s=0;
167  rho_T=rho_S=rho=0;
168  i=V_i_b; j=V_i_e;
169  while (i <= j){
170  while( (i <= V_i_e) && (x[i] <= lambda) ){
171  if (x[i]> lambda_T){
172  s_T+=x[i]; rho_T++;
173  }
174  i++;
175  }
176  while( (j>=V_i_b) && (x[j] > lambda) ){
177  if (x[j] > lambda_S){
178  s_S+=x[j]; rho_S++;
179  }
180  else{
181  s+=x[j]; rho++;
182  }
183  j--;
184  }
185  if (i<j){
186  if (x[i] > lambda_S){
187  s_S+=x[i]; rho_S++;
188  }
189  else{
190  s+=x[i]; rho++;
191  }
192
193  if (x[j]> lambda_T){
194  s_T+=x[j]; rho_T++;
195  }
196
197  temp=x[i]; x[i]=x[j]; x[j]=temp;
198  i++; j--;
199  }
200  }
201
202  s_S+=s_2; rho_S+=rho_2;
203  s+=s_S; rho+=rho_S;
204  s_T+=s; rho_T+=rho;
205  f_lambda_S=s_S-rho_S*lambda_S-z;
206  f_lambda=s-rho*lambda-z;
207  f_lambda_T=s_T-rho_T*lambda_T-z;
208
209  /*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);*/
210
211  if ( fabs(f_lambda)< delta ){
212  /*printf("\n lambda");*/
213  flag=1;
214  break;
215  }
216  if ( fabs(f_lambda_S)< delta ){
217  /* printf("\n lambda_S");*/
218  lambda=lambda_S; flag=1;
219  break;
220  }
221  if ( fabs(f_lambda_T)< delta ){
222  /* printf("\n lambda_T");*/
223  lambda=lambda_T; flag=1;
224  break;
225  }
226
227  /*
228  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);
229  printf("\n lambda_1=%5.6f, lambda_2=%5.6f, lambda=%5.6f",lambda_1, lambda_2, lambda);
230  printf("\n rho_1=%d, rho_2=%d, rho=%d ",rho_1, rho_2, rho);
231  */
232
233  if (f_lambda <0){
234  lambda_2=lambda; s_2=s; rho_2=rho;
235  f_lambda_2=f_lambda;
236
237  lambda_1=lambda_T; s_1=s_T; rho_1=rho_T;
238  f_lambda_1=f_lambda_T;
239
240  V_i_e=j; i=V_i_b;
241  while (i <= j){
242  while( (i <= V_i_e) && (x[i] <= lambda_T) ){
243  i++;
244  }
245  while( (j>=V_i_b) && (x[j] > lambda_T) ){
246  j--;
247  }
248  if (i<j){
249  x[j]=x[i];
250  i++; j--;
251  }
252  }
253  V_i_b=i; V_i=V_i_e-V_i_b+1;
254  }
255  else{
256  lambda_1=lambda; s_1=s; rho_1=rho;
257  f_lambda_1=f_lambda;
258
259  lambda_2=lambda_S; s_2=s_S; rho_2=rho_S;
260  f_lambda_2=f_lambda_S;
261
262  V_i_b=i; j=V_i_e;
263  while (i <= j){
264  while( (i <= V_i_e) && (x[i] <= lambda_S) ){
265  i++;
266  }
267  while( (j>=V_i_b) && (x[j] > lambda_S) ){
268  j--;
269  }
270  if (i<j){
271  x[i]=x[j];
272  i++; j--;
273  }
274  }
275  V_i_e=j; V_i=V_i_e-V_i_b+1;
276  }
277
278  if (V_i==0){
279  lambda=(s - z)/ rho; flag=1;
280  /*printf("\n V_i=0, lambda=%5.6f",lambda);*/
281  break;
282  }
283  }/* end of while */
284
285
286  for(i=0;i<n;i++){
287  if (v[i] > lambda)
288  x[i]=v[i]-lambda;
289  else
290  if (v[i]< -lambda)
291  x[i]=v[i]+lambda;
292  else
293  x[i]=0;
294  }
295  *root=lambda;
296  *steps=iter_step;
297 }
298
299 void epp1(double *x, double *v, int n, double rho)
300 {
301  int i;
302
303  /*
304  we assume rho>=0
305  */
306
307  for(i=0;i<n;i++){
308  if (fabs(v[i])<=rho)
309  x[i]=0;
310  else
311  if (v[i]< -rho)
312  x[i]=v[i]+rho;
313  else
314  x[i]=v[i]-rho;
315  }
316 }
317
318 void epp2(double *x, double *v, int n, double rho)
319 {
320  int i;
321  double v2=0, ratio;
322
323  /*
324  we assume rho>=0
325  */
326
327  for(i=0; i< n; i++){
328  v2+=v[i]*v[i];
329  }
330  v2=sqrt(v2);
331
332  if (rho >= v2)
333  for(i=0;i<n;i++)
334  x[i]=0;
335  else{
336  ratio= (v2-rho) /v2;
337  for(i=0;i<n;i++)
338  x[i]=v[i]*ratio;
339  }
340 }
341
342 void eppInf(double *x, double * c, int * iter_step, double *v, int n, double rho, double c0)
343 {
344  int i, steps;
345
346  /*
347  we assume rho>=0
348  */
349
350  eplb(x, c, &steps, v, n, rho, c0);
351
352  for(i=0; i< n; i++){
353  x[i]=v[i]-x[i];
354  }
355  iter_step[0]=steps;
356  iter_step[1]=0;
357 }
358
359 void zerofind(double *root, int * iterStep, double v, double p, double c, double x0)
360 {
361  double x, f, fprime, p1=p-1, pp;
362  int step=0;
363
364
365  if (v==0){
366  *root=0; *iterStep=0; return;
367  }
368
369  if (c==0){
370  *root=v; * iterStep=0; return;
371  }
372
373
374  if ( (x0 <v) && (x0>0) )
375  x=x0;
376  else
377  x=v;
378
379
380  pp=pow(x, p1);
381  f= x + c* pp -v;
382
383
384  /*
385  We apply the Newton's method for solving the root
386  */
387  while (1){
388  step++;
389
390  fprime=1 + c* p1 * pp / x;
391  /*
392  The derivative at the current solution x
393  */
394
395  x = x- f/fprime; /*
396  The new solution is computed by the Newton method
397  */
398
399
400
401  if (p>2){
402  if (x>v){
403  x=v;
404  }
405  }
406  else{
407  if ( (x<0) || (x>v)){
408  x=1e-30;
409
410  f= x+c* pow(x,p1)-v;
411
412  if (f>0){ /*
413  If f(x) = x + c x^{p-1} - v <0 at x=1e-30,
414  this shows that the real root is between (0, 1e-30).
415  For numerical reason, we just set x=0
416  */
417
418  *root=x;
419  * iterStep=step;
420
421  break;
422  }
423  }
424  }
425  /*
426  This makes sure that x lies in the interval [0, v]
427  */
428
429  pp=pow(x, p1);
430  f= x + c* pp -v;
431  /*
432  The function value at the new solution
433  */
434
435  if ( fabs(f) <= delta){
436  *root=x;
437  * iterStep=step;
438  break;
439  }
440
441  if (step>=innerIter){
442  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);
444  return;
445  }
446
447  }
448
449  /*
450  printf("\n x=%e, f=%e, step=%d\n",x, f, step);
451  */
452 }
453
454 double norm(double * v, double p, int n)
455 {
456  int i;
457  double t=0;
458
459
460  /*
461  we assume that v[i]>=0
462  p>1
463  */
464
465  for(i=0;i<n;i++)
466  t+=pow(v[i], p);
467
468  return( pow(t, 1/p) );
469 };
470
471 void eppO(double *x, double * cc, int * iter_step, double *v, int n, double rho, double p)
472 {
473  int i, *flag, bisStep, newtonStep=0, totoalStep=0;
474  double vq=0, epsilon, vmax=0, vmin=1e10; /* we assume that the minimal value in |v| is less than 1e10*/
475  double q=1/(1-1/p), c, c1, c2, root, f, xp;
476
477  double x_diff=0; /* this value denotes the maximal difference of the x values computed from c1 and c2*/
478  double temp;
479  int p_n=1; /* p_n indicates the previous phi(c) is positive or negative*/
480
481  flag=(int *)malloc(sizeof(int)*n);
482
483  /*
484  compute vq, the q-norm of v
485  flag denotes the sign of v:
486  flag[i]=0 denotes v[i] is non-negative
487  flag[i]=1 denotes v[i] is negative
488  vmin and vmax are the maximal and minimal value of |v| (excluding 0)
489  */
490  for(i=0; i< n; i++){
491
492  x[i]=0;
493
494  if (v[i]==0)
495  flag[i]=0;
496  else
497  {
498  if (v[i]>0)
499  flag[i]=0;
500  else
501  {
502  flag[i]=1;
503  v[i]=-v[i];/*
504  we set v[i] to its absolute value
505  */
506  }
507
508  vq+=pow(v[i], q);
509
510
511  if (v[i]>vmax)
512  vmax=v[i];
513
514  if (v[i]<vmin)
515  vmin=v[i];
516  }
517  }
518  vq=pow(vq, 1/q);
519
520  /*
521  zero solution
522  */
523  if (rho >= vq){
524  *cc=0;
525  iter_step[0]=iter_step[1]=0;
526
527
528  for(i=0;i<n;i++){
529  if (flag[i])
530  v[i]=-v[i]; /* set the value of v[i] back*/
531  }
532
533  free(flag);
534  return;
535  }
536
537  /*
538  compute epsilon
539  initialize c1 and c2, the interval where the root lies
540  */
541  epsilon=(vq -rho)/ vq;
542  if (p>2){
543
544  if ( log((1-epsilon) * vmin) - (p-1) * log( epsilon* vmin ) >= 709 )
545  {
546  /* If this contition holds, we have c2 >= 1e308, exceeding the machine precision.
547
548  In this case, the reason is that p is too large
549  and meanwhile epsilon * vmin is typically small.
550
551  For numerical stablity, we just regard p=inf, and run eppInf
552  */
553
554
555  for(i=0;i<n;i++){
556  if (flag[i])
557  v[i]=-v[i]; /* set the value of v[i] back*/
558  }
559
560  eppInf(x, cc, iter_step, v, n, rho, 0);
561
562  free(flag);
563  return;
564  }
565
566  c1= (1-epsilon) * vmax / pow(epsilon* vmax, p-1);
567  c2= (1-epsilon) * vmin / pow(epsilon* vmin, p-1);
568  }
569  else{ /*1 < p < 2*/
570
571  c2= (1-epsilon) * vmax / pow(epsilon* vmax, p-1);
572  c1= (1-epsilon) * vmin / pow(epsilon* vmin, p-1);
573  }
574
575
576  /*
577  printf("\n c1=%e, c2=%e", c1, c2);
578  */
579
580  if (fabs(c1-c2) <= delta){
581  c=c1;
582  }
583  else
584  c=(c1+c2)/2;
585
586
587  bisStep =0;
588
589  while(1){
590  bisStep++;
591
592  /*compute the root corresponding to c*/
593  x_diff=0;
594  for(i=0;i<n;i++){
595  zerofind(&root, &newtonStep, v[i], p, c, x[i]);
596
597  temp=fabs(root-x[i]);
598  if (x_diff< temp )
599  x_diff=temp; /*x_diff denotes the largest gap to the previous solution*/
600
601  x[i]=root;
602  totoalStep+=newtonStep;
603  }
604
605  xp=norm(x, p, n);
606
607  f=rho * pow(xp, 1-p) - c;
608
609  if ( fabs(f)<=delta || fabs(c1-c2)<=delta )
610  break;
611  else{
612  if (f>0){
613  if ( (x_diff <=delta) && (p_n==0) )
614  break;
615
616  c1=c; p_n=1;
617  }
618  else{
619
620  if ( (x_diff <=delta) && (p_n==1) )
621  break;
622
623  c2=c; p_n=0;
624  }
625  }
626  c=(c1+c2)/2;
627
628  if (bisStep>=outerIter){
629
630
631  if ( fabs(c1-c2) <=delta * c2 )
632  break;
633  else{
634  printf("\n The number of bisection steps exceed %d.", outerIter);
635  printf("\n c1=%e, c2=%e, x_diff=%e, f=%e",c1,c2,x_diff,f);
637  free(flag);
638
639  return;
640  }
641  }
642
643  /*
644  printf("\n c1=%e, c2=%e, f=%e, newtonStep=%d", c1, c2, f, newtonStep);
645  */
646  }
647
648  /*
649  printf("\n c1=%e, c2=%e, x_diff=%e, f=%e, bisStep=%d, totoalStep=%d",c1,c2, x_diff, f,bisStep,totoalStep);
650  */
651
652  for(i=0;i<n;i++){
653  if (flag[i]){
654  x[i]=-x[i];
655  v[i]=-v[i];
656  }
657  }
658  free(flag);
659
660  *cc=c;
661
662  iter_step[0]=bisStep;
663  iter_step[1]=totoalStep;
664 }
665
666 void epp(double *x, double * c, int * iter_step, double * v, int n, double rho, double p, double c0){
667  if (rho <0){
668  printf("\n rho should be non-negative!");
669  exit(1);
670  }
671
672  if (p==1){
673  epp1(x, v, n, rho);
674  *c=0;
675  iter_step[0]=iter_step[1]=0;
676  }
677  else
678  if (p==2){
679  epp2(x, v, n, rho);
680  *c=0;
681  iter_step[0]=iter_step[1]=0;
682  }
683  else
684  if (p>=1e6) /* when p >=1e6, we treat it as infity*/
685  eppInf(x, c, iter_step, v, n, rho, c0);
686  else
687  eppO(x, c, iter_step, v, n, rho, p);
688 }
689
690 #endif //USE_GPL_SHOGUN

SHOGUN Machine Learning Toolbox - Documentation