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

SHOGUN Machine Learning Toolbox - Documentation