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

SHOGUN Machine Learning Toolbox - Documentation