SHOGUN
v2.0.0
|
00001 /* This program is free software: you can redistribute it and/or modify 00002 * it under the terms of the GNU General Public License as published by 00003 * the Free Software Foundation, either version 3 of the License, or 00004 * (at your option) any later version. 00005 * 00006 * This program is distributed in the hope that it will be useful, 00007 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00008 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00009 * GNU General Public License for more details. 00010 * 00011 * You should have received a copy of the GNU General Public License 00012 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00013 * 00014 * Copyright (C) 2009 - 2012 Jun Liu and Jieping Ye 00015 */ 00016 00017 #ifndef EPSGLASSO_SLEP 00018 #define EPSGLASSO_SLEP 00019 00020 #include <stdlib.h> 00021 #include <stdio.h> 00022 #include <time.h> 00023 #include <math.h> 00024 #include <shogun/lib/slep/q1/epph.h> /* This is the head file that contains the implementation of the used functions*/ 00025 00026 00027 /* 00028 Projection for sgLasso 00029 00030 min 1/2 \|X - V\|_F^2 + \lambda_1 \|X\|_1 + \lambda_2 \|X\|_{p,1} 00031 00032 Written by Jun Liu, January 15, 2010 00033 For any problem, please contact: j.liu@asu.edu 00034 00035 */ 00036 00037 void epsgLasso(double *X, double * normx, double * V, int k, int n, double lambda1, double lambda2, int pflag){ 00038 int i, j, *iter_step, nn=n*k, m; 00039 double *v, *x; 00040 double normValue,c0=0, c; 00041 00042 v=(double *)malloc(sizeof(double)*n); 00043 x=(double *)malloc(sizeof(double)*n); 00044 iter_step=(int *)malloc(sizeof(int)*2); 00045 00046 /* 00047 initialize normx 00048 */ 00049 normx[0]=normx[1]=0; 00050 00051 00052 /* 00053 X and V are k x n matrices in matlab, stored in column priority manner 00054 x corresponds a row of X 00055 00056 pflag=2: p=2 00057 pflag=0: p=inf 00058 */ 00059 00060 /* 00061 soft thresholding 00062 by lambda1 00063 00064 the results are stored in X 00065 */ 00066 for (i=0;i<nn;i++){ 00067 if (V[i]< -lambda1) 00068 X[i]=V[i] + lambda1; 00069 else 00070 if (V[i]> lambda1) 00071 X[i]=V[i] - lambda1; 00072 else 00073 X[i]=0; 00074 } 00075 00076 /* 00077 Shrinkage or Truncating 00078 by lambda2 00079 */ 00080 if (pflag==2){ 00081 for(i=0; i<k; i++){ 00082 00083 /* 00084 process the i-th row, and store it in v 00085 */ 00086 normValue=0; 00087 00088 m=n%5; 00089 for(j=0;j<m;j++){ 00090 v[j]=X[i + j*k]; 00091 } 00092 for(j=m;j<n;j+=5){ 00093 v[j ]=X[i + j*k]; 00094 v[j+1]=X[i + (j+1)*k ]; 00095 v[j+2]=X[i + (j+2)*k]; 00096 v[j+3]=X[i + (j+3)*k]; 00097 v[j+4]=X[i + (j+4)*k]; 00098 } 00099 00100 m=n%5; 00101 for(j=0;j<m;j++){ 00102 normValue+=v[j]*v[j]; 00103 } 00104 for(j=m;j<n;j+=5){ 00105 normValue+=v[j]*v[j]+ 00106 v[j+1]*v[j+1]+ 00107 v[j+2]*v[j+2]+ 00108 v[j+3]*v[j+3]+ 00109 v[j+4]*v[j+4]; 00110 } 00111 00112 /* 00113 for(j=0; j<n; j++){ 00114 v[j]=X[i + j*k]; 00115 00116 normValue+=v[j]*v[j]; 00117 } 00118 */ 00119 00120 normValue=sqrt(normValue); 00121 00122 if (normValue<= lambda2){ 00123 for(j=0; j<n; j++) 00124 X[i + j*k]=0; 00125 00126 /*normx needs not to be updated*/ 00127 } 00128 else{ 00129 00130 normx[1]+=normValue-lambda2; 00131 /*update normx[1]*/ 00132 00133 normValue=(normValue-lambda2)/normValue; 00134 00135 m=n%5; 00136 for(j=0;j<m;j++){ 00137 X[i + j*k]*=normValue; 00138 normx[0]+=fabs(X[i + j*k]); 00139 } 00140 for(j=m; j<n;j+=5){ 00141 X[i + j*k]*=normValue; 00142 X[i + (j+1)*k]*=normValue; 00143 X[i + (j+2)*k]*=normValue; 00144 X[i + (j+3)*k]*=normValue; 00145 X[i + (j+4)*k]*=normValue; 00146 00147 normx[0]+=fabs(X[i + j*k])+ 00148 fabs(X[i + (j+1)*k])+ 00149 fabs(X[i + (j+2)*k])+ 00150 fabs(X[i + (j+3)*k])+ 00151 fabs(X[i + (j+4)*k]); 00152 } 00153 00154 /* 00155 for(j=0; j<n; j++) 00156 X[i + j*k]*=normValue; 00157 */ 00158 } 00159 } 00160 } 00161 else{ 00162 for(i=0; i<k; i++){ 00163 00164 /* 00165 process the i-th row, and store it in v 00166 */ 00167 normValue=0; 00168 for(j=0; j<n; j++){ 00169 v[j]=X[i + j*k]; 00170 00171 normValue+=fabs(v[j]); 00172 } 00173 00174 if (normValue<= lambda2){ 00175 for(j=0; j<n; j++) 00176 X[i + j*k]=0; 00177 } 00178 else{ 00179 eplb(x, &c, iter_step, v, n, lambda2, c0); 00180 00181 for(j=0; j<n; j++){ 00182 if (X[i + j*k] > c) 00183 X[i + j*k]=c; 00184 else 00185 if (X[i + j*k]<-c) 00186 X[i + j*k]=-c; 00187 } 00188 } 00189 } 00190 } 00191 00192 00193 free(v); 00194 free(x); 00195 free(iter_step); 00196 } 00197 #endif /* ----- #ifndef EPSGLASSO_SLEP ----- */ 00198