flsa.cpp

Go to the documentation of this file.
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation