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