SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
epsp.h
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 #ifndef EPSP_SLEP
18 #define EPSP_SLEP
19 
20 #include <shogun/lib/config.h>
21 
22 #include <stdlib.h>
23 #include <stdio.h>
24 #include <time.h>
25 #include <math.h>
26 
27 #define delta 1e-12
28 
29 /*
30  Euclidean Projection onto the simplex (epsp)
31 
32  min 1/2 ||x- y||_2^2
33  s.t. ||x||_1 = z, x >=0
34 
35 which is converted to the following zero finding problem
36 
37  f(lambda)= sum( max( x-lambda,0) )-z=0
38 
39  Usage:
40  [x, lambda, iter_step]=epsp(y, n, z, lambda0);
41 
42  */
43 
44 void epsp(double * x, double *root, int * stepsp, double * v,
45 int n, double z, double lambda0)
46 {
47 
48  int i, j, flag=0;
49  int rho_1, rho_2, rho, rho_T, rho_S;
50  int V_i_b, V_i_e, V_i;
51  double lambda_1, lambda_2, lambda_T, lambda_S, lambda;
52  double s_1, s_2, s, s_T, s_S, v_max, temp;
53  double f_lambda_1, f_lambda_2, f_lambda, f_lambda_T, f_lambda_S;
54  int iter_step=0;
55 
56 
57 
58  if (z< 0){
59  printf("\n z should be nonnegative!");
60  exit(-1);
61  }
62 
63 
64  /*
65  * find the maximal value in v
66  */
67  v_max=v[0];
68  for (i=1;i<n; i++){
69  if (v[i] > v_max)
70  v_max=v[i];
71  }
72 
73 
74  lambda_1=v_max - z; lambda_2=v_max;
75  /*
76  * copy v to x
77  * compute f_lambda_1, rho_1, s_1
78  */
79  V_i=0;s_1=0; rho_1=0;
80  for (i=0;i<n;i++){
81  if (v[i] > lambda_1){
82  x[V_i]=v[i];
83 
84  s_1+=x[V_i]; rho_1++;
85 
86  V_i++;
87  }
88  }
89  f_lambda_1=s_1-rho_1* lambda_1 -z;
90 
91  rho_2=0; s_2=0; f_lambda_2=-z;
92  V_i_b=0; V_i_e=V_i-1;
93 
94  lambda=lambda0;
95  if ( (lambda<lambda_2) && (lambda> lambda_1) ){
96  /*-------------------------------------------------------------------
97  Initialization with the root
98  *-------------------------------------------------------------------
99  */
100 
101  i=V_i_b; j=V_i_e; rho=0; s=0;
102  while (i <= j){
103  while( (i <= V_i_e) && (x[i] <= lambda) ){
104  i++;
105  }
106  while( (j>=V_i_b) && (x[j] > lambda) ){
107  s+=x[j];
108  j--;
109  }
110  if (i<j){
111  s+=x[i];
112 
113  temp=x[i]; x[i]=x[j]; x[j]=temp;
114  i++; j--;
115  }
116  }
117 
118  rho=V_i_e-j; rho+=rho_2; s+=s_2;
119  f_lambda=s-rho*lambda-z;
120 
121  if ( fabs(f_lambda)< delta ){
122  flag=1;
123  }
124 
125  if (f_lambda <0){
126  lambda_2=lambda; s_2=s; rho_2=rho; f_lambda_2=f_lambda;
127 
128  V_i_e=j; V_i=V_i_e-V_i_b+1;
129  }
130  else{
131  lambda_1=lambda; rho_1=rho; s_1=s; f_lambda_1=f_lambda;
132 
133  V_i_b=i; V_i=V_i_e-V_i_b+1;
134  }
135 
136  if (V_i==0){
137  /*printf("\n rho=%d, rho_1=%d, rho_2=%d",rho, rho_1, rho_2);*/
138 
139  /*printf("\n V_i=%d",V_i);*/
140 
141  lambda=(s - z)/ rho;
142  flag=1;
143  }
144  /*-------------------------------------------------------------------
145  End of initialization
146  *--------------------------------------------------------------------
147  */
148 
149  }/* end of if(!flag) */
150 
151  while (!flag){
152  iter_step++;
153 
154  /* compute lambda_T */
155  lambda_T=lambda_1 + f_lambda_1 /rho_1;
156  if(rho_2 !=0){
157  if (lambda_2 + f_lambda_2 /rho_2 > lambda_T)
158  lambda_T=lambda_2 + f_lambda_2 /rho_2;
159  }
160 
161  /* compute lambda_S */
162  lambda_S=lambda_2 - f_lambda_2 *(lambda_2-lambda_1)/(f_lambda_2-f_lambda_1);
163 
164  if (fabs(lambda_T-lambda_S) <= delta){
165  lambda=lambda_T; flag=1;
166  break;
167  }
168 
169  /* set lambda as the middle point of lambda_T and lambda_S */
170  lambda=(lambda_T+lambda_S)/2;
171 
172  s_T=s_S=s=0;
173  rho_T=rho_S=rho=0;
174  i=V_i_b; j=V_i_e;
175  while (i <= j){
176  while( (i <= V_i_e) && (x[i] <= lambda) ){
177  if (x[i]> lambda_T){
178  s_T+=x[i]; rho_T++;
179  }
180  i++;
181  }
182  while( (j>=V_i_b) && (x[j] > lambda) ){
183  if (x[j] > lambda_S){
184  s_S+=x[j]; rho_S++;
185  }
186  else{
187  s+=x[j]; rho++;
188  }
189  j--;
190  }
191  if (i<j){
192  if (x[i] > lambda_S){
193  s_S+=x[i]; rho_S++;
194  }
195  else{
196  s+=x[i]; rho++;
197  }
198 
199  if (x[j]> lambda_T){
200  s_T+=x[j]; rho_T++;
201  }
202 
203  temp=x[i]; x[i]=x[j]; x[j]=temp;
204  i++; j--;
205  }
206  }
207 
208  s_S+=s_2; rho_S+=rho_2;
209  s+=s_S; rho+=rho_S;
210  s_T+=s; rho_T+=rho;
211  f_lambda_S=s_S-rho_S*lambda_S-z;
212  f_lambda=s-rho*lambda-z;
213  f_lambda_T=s_T-rho_T*lambda_T-z;
214 
215  /*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);*/
216 
217  if ( fabs(f_lambda)< delta ){
218  /*printf("\n lambda");*/
219  flag=1;
220  break;
221  }
222  if ( fabs(f_lambda_S)< delta ){
223  /* printf("\n lambda_S");*/
224  lambda=lambda_S; flag=1;
225  break;
226  }
227  if ( fabs(f_lambda_T)< delta ){
228  /* printf("\n lambda_T");*/
229  lambda=lambda_T; flag=1;
230  break;
231  }
232 
233  /*
234  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);
235  printf("\n lambda_1=%5.6f, lambda_2=%5.6f, lambda=%5.6f",lambda_1, lambda_2, lambda);
236  printf("\n rho_1=%d, rho_2=%d, rho=%d ",rho_1, rho_2, rho);
237  */
238 
239  if (f_lambda <0){
240  lambda_2=lambda; s_2=s; rho_2=rho;
241  f_lambda_2=f_lambda;
242 
243  lambda_1=lambda_T; s_1=s_T; rho_1=rho_T;
244  f_lambda_1=f_lambda_T;
245 
246  V_i_e=j; i=V_i_b;
247  while (i <= j){
248  while( (i <= V_i_e) && (x[i] <= lambda_T) ){
249  i++;
250  }
251  while( (j>=V_i_b) && (x[j] > lambda_T) ){
252  j--;
253  }
254  if (i<j){
255  x[j]=x[i];
256  i++; j--;
257  }
258  }
259  V_i_b=i; V_i=V_i_e-V_i_b+1;
260  }
261  else{
262  lambda_1=lambda; s_1=s; rho_1=rho;
263  f_lambda_1=f_lambda;
264 
265  lambda_2=lambda_S; s_2=s_S; rho_2=rho_S;
266  f_lambda_2=f_lambda_S;
267 
268  V_i_b=i; j=V_i_e;
269  while (i <= j){
270  while( (i <= V_i_e) && (x[i] <= lambda_S) ){
271  i++;
272  }
273  while( (j>=V_i_b) && (x[j] > lambda_S) ){
274  j--;
275  }
276  if (i<j){
277  x[i]=x[j];
278  i++; j--;
279  }
280  }
281  V_i_e=j; V_i=V_i_e-V_i_b+1;
282  }
283 
284  if (V_i==0){
285  lambda=(s - z)/ rho; flag=1;
286  /*printf("\n V_i=0, lambda=%5.6f",lambda);*/
287  break;
288  }
289  }/* end of while */
290 
291 
292  for(i=0;i<n;i++){
293  if (v[i] > lambda)
294  x[i]=v[i]-lambda;
295  else
296  x[i]=0;
297  }
298  *root=lambda;
299  *stepsp=iter_step;
300 }
301 #endif /* ----- #ifndef EPSP_SLEP ----- */
302 

SHOGUN Machine Learning Toolbox - Documentation