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

SHOGUN Machine Learning Toolbox - Documentation