SHOGUN  5.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
epsgLasso.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 
18 #ifndef EPSGLASSO_SLEP
19 #define EPSGLASSO_SLEP
20 
21 #include <shogun/lib/config.h>
22 #ifdef USE_GPL_SHOGUN
23 
24 #include <stdlib.h>
25 #include <stdio.h>
26 #include <time.h>
27 #include <math.h>
28 #include <shogun/lib/slep/q1/epph.h> /* This is the head file that contains the implementation of the used functions*/
29 
30 
31 /*
32  Projection for sgLasso
33 
34  min 1/2 \|X - V\|_F^2 + \lambda_1 \|X\|_1 + \lambda_2 \|X\|_{p,1}
35 
36  Written by Jun Liu, January 15, 2010
37  For any problem, please contact: j.liu@asu.edu
38 
39  */
40 
41 void epsgLasso(double *X, double * normx, double * V, int k, int n, double lambda1, double lambda2, int pflag){
42  int i, j, *iter_step, nn=n*k, m;
43  double *v, *x;
44  double normValue,c0=0, c;
45 
46  v=(double *)malloc(sizeof(double)*n);
47  x=(double *)malloc(sizeof(double)*n);
48  iter_step=(int *)malloc(sizeof(int)*2);
49 
50  /*
51  initialize normx
52  */
53  normx[0]=normx[1]=0;
54 
55 
56  /*
57  X and V are k x n matrices in matlab, stored in column priority manner
58  x corresponds a row of X
59 
60  pflag=2: p=2
61  pflag=0: p=inf
62  */
63 
64  /*
65  soft thresholding
66  by lambda1
67 
68  the results are stored in X
69  */
70  for (i=0;i<nn;i++){
71  if (V[i]< -lambda1)
72  X[i]=V[i] + lambda1;
73  else
74  if (V[i]> lambda1)
75  X[i]=V[i] - lambda1;
76  else
77  X[i]=0;
78  }
79 
80  /*
81  Shrinkage or Truncating
82  by lambda2
83  */
84  if (pflag==2){
85  for(i=0; i<k; i++){
86 
87  /*
88  process the i-th row, and store it in v
89  */
90  normValue=0;
91 
92  m=n%5;
93  for(j=0;j<m;j++){
94  v[j]=X[i + j*k];
95  }
96  for(j=m;j<n;j+=5){
97  v[j ]=X[i + j*k];
98  v[j+1]=X[i + (j+1)*k ];
99  v[j+2]=X[i + (j+2)*k];
100  v[j+3]=X[i + (j+3)*k];
101  v[j+4]=X[i + (j+4)*k];
102  }
103 
104  m=n%5;
105  for(j=0;j<m;j++){
106  normValue+=v[j]*v[j];
107  }
108  for(j=m;j<n;j+=5){
109  normValue+=v[j]*v[j]+
110  v[j+1]*v[j+1]+
111  v[j+2]*v[j+2]+
112  v[j+3]*v[j+3]+
113  v[j+4]*v[j+4];
114  }
115 
116  /*
117  for(j=0; j<n; j++){
118  v[j]=X[i + j*k];
119 
120  normValue+=v[j]*v[j];
121  }
122  */
123 
124  normValue=sqrt(normValue);
125 
126  if (normValue<= lambda2){
127  for(j=0; j<n; j++)
128  X[i + j*k]=0;
129 
130  /*normx needs not to be updated*/
131  }
132  else{
133 
134  normx[1]+=normValue-lambda2;
135  /*update normx[1]*/
136 
137  normValue=(normValue-lambda2)/normValue;
138 
139  m=n%5;
140  for(j=0;j<m;j++){
141  X[i + j*k]*=normValue;
142  normx[0]+=fabs(X[i + j*k]);
143  }
144  for(j=m; j<n;j+=5){
145  X[i + j*k]*=normValue;
146  X[i + (j+1)*k]*=normValue;
147  X[i + (j+2)*k]*=normValue;
148  X[i + (j+3)*k]*=normValue;
149  X[i + (j+4)*k]*=normValue;
150 
151  normx[0]+=fabs(X[i + j*k])+
152  fabs(X[i + (j+1)*k])+
153  fabs(X[i + (j+2)*k])+
154  fabs(X[i + (j+3)*k])+
155  fabs(X[i + (j+4)*k]);
156  }
157 
158  /*
159  for(j=0; j<n; j++)
160  X[i + j*k]*=normValue;
161  */
162  }
163  }
164  }
165  else{
166  for(i=0; i<k; i++){
167 
168  /*
169  process the i-th row, and store it in v
170  */
171  normValue=0;
172  for(j=0; j<n; j++){
173  v[j]=X[i + j*k];
174 
175  normValue+=fabs(v[j]);
176  }
177 
178  if (normValue<= lambda2){
179  for(j=0; j<n; j++)
180  X[i + j*k]=0;
181  }
182  else{
183  eplb(x, &c, iter_step, v, n, lambda2, c0);
184 
185  for(j=0; j<n; j++){
186  if (X[i + j*k] > c)
187  X[i + j*k]=c;
188  else
189  if (X[i + j*k]<-c)
190  X[i + j*k]=-c;
191  }
192  }
193  }
194  }
195 
196 
197  free(v);
198  free(x);
199  free(iter_step);
200 }
201 #endif //USE_GPL_SHOGUN
202 #endif /* ----- #ifndef EPSGLASSO_SLEP ----- */
203 

SHOGUN Machine Learning Toolbox - Documentation