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 #include <shogun/lib/slep/flsa/flsa.h> 00018 #include <shogun/lib/slep/flsa/sfa.h> 00019 #include <stdlib.h> 00020 #include <stdio.h> 00021 #include <time.h> 00022 #include <math.h> 00023 00024 void flsa(double *x, double *z, double *infor, 00025 double * v, double *z0, 00026 double lambda1, double lambda2, int n, 00027 int maxStep, double tol, int tau, int flag) 00028 { 00029 00030 int i, nn=n-1, m; 00031 double zMax, temp; 00032 double *Av, *g, *s; 00033 int iterStep, numS; 00034 double gap; 00035 double *zz = NULL; /*to replace z0, so that z0 shall not revised after */ 00036 00037 00038 Av=(double *) malloc(sizeof(double)*nn); 00039 00040 /* 00041 Compute Av= A*v (n=4, nn=3) 00042 A= [ -1 1 0 0; 00043 0 -1 1 0; 00044 0 0 -1 1] 00045 */ 00046 00047 for (i=0;i<nn; i++) 00048 Av[i]=v[i+1]-v[i]; 00049 00050 /* 00051 Sovlve the linear system via Thomas's algorithm or Rose's algorithm 00052 B * z0 = Av 00053 */ 00054 00055 Thomas(&zMax, z, Av, nn); 00056 00057 /* 00058 We consider two cases: 00059 1) lambda2 >= zMax, which leads to a solution with same entry values 00060 2) lambda2 < zMax, which needs to first run sfa, and then perform soft thresholding 00061 */ 00062 00063 /* 00064 First case: lambda2 >= zMax 00065 */ 00066 if (lambda2 >= zMax){ 00067 00068 temp=0; 00069 m=n%5; 00070 if (m!=0){ 00071 for (i=0;i<m;i++) 00072 temp+=v[i]; 00073 } 00074 for (i=m;i<n;i+=5){ 00075 temp += v[i] + v[i+1] + v[i+2] + v[i+3] + v[i+4]; 00076 } 00077 temp/=n; 00078 /* temp is the mean value of v*/ 00079 00080 00081 /* 00082 soft thresholding by lambda1 00083 */ 00084 if (temp> lambda1) 00085 temp= temp-lambda1; 00086 else 00087 if (temp < -lambda1) 00088 temp= temp+lambda1; 00089 else 00090 temp=0; 00091 00092 m=n%7; 00093 if (m!=0){ 00094 for (i=0;i<m;i++) 00095 x[i]=temp; 00096 } 00097 for (i=m;i<n;i+=7){ 00098 x[i] =temp; 00099 x[i+1] =temp; 00100 x[i+2] =temp; 00101 x[i+3] =temp; 00102 x[i+4] =temp; 00103 x[i+5] =temp; 00104 x[i+6] =temp; 00105 } 00106 00107 gap=0; 00108 00109 free(Av); 00110 00111 if (infor) 00112 { 00113 infor[0]= gap; 00114 infor[1]= 0; 00115 infor[2]=zMax; 00116 infor[3]=0; 00117 } 00118 00119 return; 00120 } 00121 00122 00123 /* 00124 Second case: lambda2 < zMax 00125 00126 We need to call sfa for computing x, and then do soft thresholding 00127 00128 Before calling sfa, we need to allocate memory for g and s, 00129 and initialize z and z0. 00130 */ 00131 00132 00133 /* 00134 Allocate memory for g and s 00135 */ 00136 00137 g =(double *) malloc(sizeof(double)*nn), 00138 s =(double *) malloc(sizeof(double)*nn); 00139 00140 00141 00142 m=flag /10; 00143 /* 00144 00145 If m=0, then this shows that, z0 is a "good" starting point. (m=1-6) 00146 00147 Otherwise (m=11-16), we shall set z as either the solution to the linear system. 00148 or the zero point 00149 00150 */ 00151 if (m==0){ 00152 for (i=0;i<nn;i++){ 00153 if (z0[i] > lambda2) 00154 z[i]=lambda2; 00155 else 00156 if (z0[i]<-lambda2) 00157 z[i]=-lambda2; 00158 else 00159 z[i]=z0[i]; 00160 } 00161 } 00162 else{ 00163 if (lambda2 >= 0.5 * zMax){ 00164 for (i=0;i<nn;i++){ 00165 if (z[i] > lambda2) 00166 z[i]=lambda2; 00167 else 00168 if (z[i]<-lambda2) 00169 z[i]=-lambda2; 00170 } 00171 } 00172 else{ 00173 for (i=0;i<nn;i++) 00174 z[i]=0; 00175 00176 } 00177 } 00178 00179 flag=flag %10; /* 00180 flag is now in [1:6] 00181 00182 for sfa, i.e., flag in [1:4], we need initialize z0 with zero 00183 */ 00184 00185 if (flag>=1 && flag<=4){ 00186 zz =(double *) malloc(sizeof(double)*nn); 00187 00188 for (i=0;i<nn;i++) 00189 zz[i]=0; 00190 } 00191 00192 /* 00193 call sfa, sfa_one, or sfa_special to compute z, for finding the subgradient 00194 and x 00195 */ 00196 00197 if (flag==6) 00198 iterStep=sfa_one(x, &gap, &numS, 00199 z, v, Av, 00200 lambda2, nn, maxStep, 00201 s, g, 00202 tol, tau); 00203 else 00204 if (flag==5) 00205 iterStep=sfa_special(x, &gap, &numS, 00206 z, v, Av, 00207 lambda2, nn, maxStep, 00208 s, g, 00209 tol, tau); 00210 else{ 00211 iterStep=sfa(x, &gap, &numS, 00212 z, zz, v, Av, 00213 lambda2, nn, maxStep, 00214 s, g, 00215 tol,tau, flag); 00216 00217 free (zz); 00218 /*free the variable zz*/ 00219 } 00220 00221 00222 /* 00223 soft thresholding by lambda1 00224 */ 00225 00226 for(i=0;i<n;i++) 00227 if (x[i] > lambda1) 00228 x[i]-=lambda1; 00229 else 00230 if (x[i]<-lambda1) 00231 x[i]+=lambda1; 00232 else 00233 x[i]=0; 00234 00235 00236 free(Av); 00237 free(g); 00238 free(s); 00239 00240 if (infor) 00241 { 00242 infor[0]=gap; 00243 infor[1]=iterStep; 00244 infor[2]=zMax; 00245 infor[3]=numS; 00246 } 00247 } 00248