SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
epph.cpp
Go to the documentation of this file.
1 /* This program is free software: you can redistribute it and/or modify
2  * it under the terms of the GNU General Public License as published by
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);
441  printf("\n If you meet with this problem, please contact Jun Liu (j.liu@asu.edu). Thanks!");
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);
634  printf("\n If you meet with this problem, please contact Jun Liu (j.liu@asu.edu). Thanks!");
635 
636  return;
637  }
638  }
639 
640  /*
641  printf("\n c1=%e, c2=%e, f=%e, newtonStep=%d", c1, c2, f, newtonStep);
642  */
643  }
644 
645  /*
646  printf("\n c1=%e, c2=%e, x_diff=%e, f=%e, bisStep=%d, totoalStep=%d",c1,c2, x_diff, f,bisStep,totoalStep);
647  */
648 
649  for(i=0;i<n;i++){
650  if (flag[i]){
651  x[i]=-x[i];
652  v[i]=-v[i];
653  }
654  }
655  free(flag);
656 
657  *cc=c;
658 
659  iter_step[0]=bisStep;
660  iter_step[1]=totoalStep;
661 }
662 
663 void epp(double *x, double * c, int * iter_step, double * v, int n, double rho, double p, double c0){
664  if (rho <0){
665  printf("\n rho should be non-negative!");
666  exit(1);
667  }
668 
669  if (p==1){
670  epp1(x, v, n, rho);
671  *c=0;
672  iter_step[0]=iter_step[1]=0;
673  }
674  else
675  if (p==2){
676  epp2(x, v, n, rho);
677  *c=0;
678  iter_step[0]=iter_step[1]=0;
679  }
680  else
681  if (p>=1e6) /* when p >=1e6, we treat it as infity*/
682  eppInf(x, c, iter_step, v, n, rho, c0);
683  else
684  eppO(x, c, iter_step, v, n, rho, p);
685 }
686 

SHOGUN Machine Learning Toolbox - Documentation