SHOGUN  v2.0.0
sfa.cpp
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
18 #include <stdlib.h>
19 #include <stdio.h>
20 #include <time.h>
21 #include <math.h>
22
23 #define delta 1e-10
24
25 /*
26  Revision History
27
28  First Version available on October 10, 2009
29
30  A runnable version on October 15, 2009
31
32  Major revision on October 29, 2009
33  (Some functions appearing in a previous version have deleted, please refer to the previous version for the old functions.
34  Some new functions have been added as well)
35
36 */
37
38 /*
39
40  Files contained in this header file sfa.h:
41
42  1. Algorithms for solving the linear system A A^T z0 = Av (see the description of A from the following context)
43
44  void Thomas(double *zMax, double *z0,
45  double * Av, int nn)
46
47  void Rose(double *zMax, double *z0,
48  double * Av, int nn)
49
50  int supportSet(double *x, double *v, double *z,
51  double *g, int * S, double lambda, int nn)
52
53  void dualityGap(double *gap, double *z,
54  double *g, double *s, double *Av,
55  double lambda, int nn)
56
57  void dualityGap2(double *gap, double *z,
58  double *g, double *s, double *Av,
59  double lambda, int nn)
60
61
62  2. The Subgraident Finding Algorithm (SFA) for solving problem (4) (refer to the description of the problem for detail)
63
64  int sfa(double *x, double *gap,
65  double *z, double *z0, double * v, double * Av,
66  double lambda, int nn, int maxStep,
67  double *s, double *g,
68  double tol, int tau, int flag)
69
70  int sfa_special(double *x, double *gap,
71  double *z, double * v, double * Av,
72  double lambda, int nn, int maxStep,
73  double *s, double *g,
74  double tol, int tau)
75
76  int sfa_one(double *x, double *gap,
77  double *z, double * v, double * Av,
78  double lambda, int nn, int maxStep,
79  double *s, double *g,
80  double tol, int tau)
81
82
83 */
84
85
86 /*
87
88  Some mathematical background.
89
90  In this file, we discuss how to solve the following subproblem,
91
92  min_x 1/2 \|x-v\|^2 + lambda \|A x\|_1, (1)
93
94  which is a key problem used in the Fused Lasso Signal Approximator (FLSA).
95
96  Also, note that, FLSA is a building block for solving the optimation problmes with fused Lasso penalty.
97
98  In (1), x and v are n-dimensional vectors,
99  and A is a matrix with size (n-1) x n, and is defined as follows (e.g., n=4):
100  A= [ -1 1 0 0;
101  0 -1 1 0;
102  0 0 -1 1]
103
104  The above problem can be reformulated as the following equivalent min-max optimization problem
105
106  min_x max_z 1/2 \|x-v\|^2 + <A x, z>
107  subject to \|z\|_{infty} \leq lambda (2)
108
109
110  It is easy to get that, at the optimal point
111
112  x = v - AT z, (3)
113
114  where z is the optimal solution to the following optimization problem
115
116  min_z 1/2 z^T A AT z - < z, A v>,
117  subject to \|z\|_{infty} \leq lambda (4)
118
119
120
121  Let B=A A^T. It is easy to get that B is a (n-1) x (n-1) tridiagonal matrix.
122  When n=5, B is defined as:
123  B= [ 2 -1 0 0;
124  -1 2 -1 0;
125  0 -1 2 -1;
126  0 0 -1 2]
127
128  Let z0 be the solution to the linear system:
129
130  A A^T * z0 = A * v (5)
131
132  The problem (5) can be solve by the Thomas Algorithm, in about 5n multiplications and 4n additions.
133
134  It can also be solved by the Rose's Algorithm, in about 2n multiplications and 2n additions.
135
136  Moreover, considering the special structure of the matrix A (and B),
137  it can be solved in about n multiplications and 3n additions
138
139  If lambda \geq \|z0\|_{infty}, x_i= mean(v), for all i,
140  the problem (1) admits near analytical solution
141
142
143  We have also added the restart technique, please refer to our paper for detail!
144
145 */
146
147
148 /*
150 */
151
152 void Thomas(double *zMax, double *z0, double * Av, int nn){
153
154  /*
155
156  We apply the Tomas algorithm for solving the following linear system
157  B * z0 = Av
158  Thomas algorithm is also called the tridiagonal matrix algorithm
159
160  B=[ 2 -1 0 0;
161  -1 2 -1 0;
162  0 -1 2 -1;
163  0 0 -1 2]
164
165  z0 is the result, Av is unchanged after the computation
166
167
168  c is a precomputed nn dimensional vector
169  c=[-1/2, -2/3, -3/4, -4/5, ..., -nn/(nn+1)]
170
171  c[i]=- (i+1) / (i+2)
172  c[i-1]=- i / (i+1)
173
174  z0 is an nn dimensional vector
175
176 */
177
178  int i;
179  double tt, z_max;
180
181  /*
182  Modify the coefficients in Av (copy to z0)
183  */
184  z0[0]=Av[0]/2;
185  for (i=1;i < nn; i++){
186  tt=Av[i] + z0[i-1];
187  z0[i]=tt - tt / (i+2);
188  }
189
190  /*z0[i]=(Av[i] + z0[i-1]) * (i+1) / (i+2);*/
191
192  /*z0[i]=(Av[i] + z0[i-1])/ ( 2 - i / (i+1));*/
193
194
195  /*
196  Back substitute (obtain the result in z0)
197  */
198  z_max= fabs(z0[nn-1]);
199
200  for (i=nn-2; i>=0; i--){
201
202  z0[i]+= z0[i+1] - z0[i+1]/ (i+2);
203
204  /*z0[i]+= z0[i+1] * (i+1) / (i+2);*/
205
206  tt=fabs(z0[i]);
207
208  if (tt > z_max)
209  z_max=tt;
210
211  }
212  *zMax=z_max;
213
214 }
215
216
217
218
219 /*
221 */
222
223 void Rose(double *zMax, double *z0, double * Av, int nn){
224
225  /*
226  We use the Rose algorithm for solving the following linear system
227  B * z0 = Av
228
229
230  B=[ 2 -1 0 0;
231  -1 2 -1 0;
232  0 -1 2 -1;
233  0 0 -1 2]
234
235  z0 is the result, Av is unchanged after the computation
236
237  z0 is an nn dimensional vector
238
239 */
240
241  int i, m;
242  double s=0, z_max;
243
244
245  /*
246  We follow the style in CLAPACK
247  */
248  m= nn % 5;
249  if (m!=0){
250  for (i=0;i<m; i++)
251  s+=Av[i] * (i+1);
252  }
253  for(i=m;i<nn;i+=5)
254  s+= Av[i] * (i+1)
255  + Av[i+1] * (i+2)
256  + Av[i+2] * (i+3)
257  + Av[i+3] * (i+4)
258  + Av[i+4] * (i+5);
259  s/=(nn+1);
260
261
262  /*
263  from nn-1 to 0
264  */
265  z0[nn-1]=Av[nn-1]- s;
266  for (i=nn-2;i >=0; i--){
267  z0[i]=Av[i] + z0[i+1];
268  }
269
270  /*
271  from 0 to nn-1
272  */
273  z_max= fabs(z0[0]);
274  for (i=0; i<nn; i++){
275
276  z0[i]+= z0[i-1];
277
278  s=fabs(z0[i]);
279
280  if (s > z_max)
281  z_max=s;
282
283  }
284  *zMax=z_max;
285
286 }
287
288
289
290 /*
292
293 x=omega(z)
294
295 v: the vector to be projected
296 z: the approximate solution
297 g: the gradient at z (g should be computed before calling this function
298
299 nn: the length of z, g, and S (maximal length for S)
300
301 n: the length of x and v
302
303 S: records the indices of the elements in the support set
304 */
305
306 int supportSet(double *x, double *v, double *z, double *g, int * S, double lambda, int nn){
307
308  int i, j, n=nn+1, numS=0;
309  double temp;
310
311
312  /*
313  we first scan z and g to obtain the support set S
314  */
315
316  /*numS: number of the elements in the support set S*/
317  for(i=0;i<nn; i++){
318  if ( ( (z[i]==lambda) && (g[i] < delta) ) || ( (z[i]==-lambda) && (g[i] >delta) )){
319  S[numS]=i;
320  numS++;
321  }
322  }
323
324  /*
325  printf("\n %d",numS);
326  */
327
328  if (numS==0){ /*this shows that S is empty*/
329  temp=0;
330  for (i=0;i<n;i++)
331  temp+=v[i];
332
333  temp=temp/n;
334  for(i=0;i<n;i++)
335  x[i]=temp;
336
337  return numS;
338  }
339
340
341  /*
342  Next, we deal with numS >=1
343  */
344
345  /*process the first block
346
347  j=0
348  */
349  temp=0;
350  for (i=0;i<=S[0]; i++)
351  temp+=v[i];
352  /*temp =sum (v [0: s[0] ]*/
353  temp=( temp + z[ S[0] ] ) / (S[0] +1);
354  for (i=0;i<=S[0]; i++)
355  x[i]=temp;
356
357
358  /*process the middle blocks
359
360  If numS=1, it belongs the last block
361  */
362  for (j=1; j < numS; j++){
363  temp=0;
364  for (i= S[j-1] +1; i<= S[j]; i++){
365  temp+=v[i];
366  }
367
368  /*temp =sum (v [ S[j-1] +1: s[j] ]*/
369
370  temp=(temp - z[ S[j-1] ] + z[ S[j] ])/ (S[j]- S[j-1]);
371
372  for (i= S[j-1] +1; i<= S[j]; i++){
373  x[i]=temp;
374  }
375  }
376
377  /*process the last block
378  j=numS-1;
379  */
380  temp=0;
381  for (i=S[numS-1] +1 ;i< n; i++)
382  temp+=v[i];
383  /*temp =sum (v [ (S[numS-1] +1): (n-1) ]*/
384
385  temp=( temp - z[ S[numS-1] ] ) / (nn - S[numS-1]); /*S[numS-1] <= nn-1*/
386
387  for (i=S[numS-1] +1 ;i< n; i++)
388  x[i]=temp;
389
390  return numS;
391
392 }
393
394
395
396 /*
397
399
400 we compute the duality corresponding the solution z
401
402 z: the approximate solution
403 g: the gradient at z (we recompute the gradient)
404 s: an auxiliary variable
405 Av: A*v
406
407 nn: the lenght for z, g, s, and Av
408
409 The variables g and s shall be revised.
410
411 The variables z and Av remain unchanged.
412 */
413
414 void dualityGap(double *gap, double *z, double *g, double *s, double *Av, double lambda, int nn){
415
416  int i, m;
417  double temp;
418
419
420  g[0]=z[0] + z[0] - z[1] - Av[0];
421  for (i=1;i<nn-1;i++){
422  g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
423  }
424  g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
425
426
427  for (i=0;i<nn;i++)
428  if (g[i]>0)
429  s[i]=lambda + z[i];
430  else
431  s[i]=-lambda + z[i];
432
433
434  temp=0;
435  m=nn%5;
436
437  if (m!=0){
438  for(i=0;i<m;i++)
439  temp+=s[i]*g[i];
440  }
441
442  for(i=m;i<nn;i+=5)
443  temp=temp + s[i] *g[i]
444  + s[i+1]*g[i+1]
445  + s[i+2]*g[i+2]
446  + s[i+3]*g[i+3]
447  + s[i+4]*g[i+4];
448  *gap=temp;
449 }
450
451
452 /*
453  Similar to dualityGap,
454
455  The difference is that, we assume that g has been computed.
456  */
457
458 void dualityGap2(double *gap, double *z, double *g, double *s, double *Av, double lambda, int nn){
459
460  int i, m;
461  double temp;
462
463
464  /*
465  g[0]=z[0] + z[0] - z[1] - Av[0];
466  for (i=1;i<nn-1;i++){
467  g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
468  }
469  g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
470
471 */
472
473  for (i=0;i<nn;i++)
474  if (g[i]>0)
475  s[i]=lambda + z[i];
476  else
477  s[i]=-lambda + z[i];
478
479
480  temp=0;
481  m=nn%5;
482
483  if (m!=0){
484  for(i=0;i<m;i++)
485  temp+=s[i]*g[i];
486  }
487
488  for(i=m;i<nn;i+=5)
489  temp=temp + s[i] *g[i]
490  + s[i+1]*g[i+1]
491  + s[i+2]*g[i+2]
492  + s[i+3]*g[i+3]
493  + s[i+4]*g[i+4];
494  *gap=temp;
495 }
496
497
498 /*
499 generateSolution:
500
501 generate the solution x based on the information of z and g
502 (!!!!we assume that g has been computed as the gradient of z!!!!)
503
504 */
505
506 int generateSolution(double *x, double *z, double *gap,
507  double *v, double *Av,
508  double *g, double *s, int *S,
509  double lambda, int nn){
510
511  int i, m, numS, n=nn+1;
512  double temp, funVal1, funVal2;
513
514  /*
515  z is the appropriate solution,
516  and g contains its gradient
517  */
518
519
520  /*
521  We assume that n>=3, and thus nn>=2
522
523  We have two ways for recovering x.
524  The first way is x = v - A^T z
525  The second way is x =omega(z)
526  */
527
528  temp=0;
529  m=nn%5;
530  if (m!=0){
531  for (i=0;i<m;i++)
532  temp+=z[i]*(g[i] + Av[i]);
533  }
534  for (i=m;i<nn;i+=5)
535  temp=temp + z[i] *(g[i] + Av[i])
536  + z[i+1]*(g[i+1] + Av[i+1])
537  + z[i+2]*(g[i+2] + Av[i+2])
538  + z[i+3]*(g[i+3] + Av[i+3])
539  + z[i+4]*(g[i+4] + Av[i+4]);
540  funVal1=temp /2;
541
542  temp=0;
543  m=nn%5;
544  if (m!=0){
545  for (i=0;i<m;i++)
546  temp+=fabs(g[i]);
547  }
548  for (i=m;i<nn;i+=5)
549  temp=temp + fabs(g[i])
550  + fabs(g[i+1])
551  + fabs(g[i+2])
552  + fabs(g[i+3])
553  + fabs(g[i+4]);
554  funVal1=funVal1+ temp*lambda;
555
556
557  /*
558  we compute the solution by the second way
559  */
560
561  numS= supportSet(x, v, z, g, S, lambda, nn);
562
563  /*
564  we compute the objective function of x computed in the second way
565  */
566
567  temp=0;
568  m=n%5;
569  if (m!=0){
570  for (i=0;i<m;i++)
571  temp+=(x[i]-v[i]) * (x[i]-v[i]);
572  }
573  for (i=m;i<n;i+=5)
574  temp=temp + (x[i]- v[i]) * ( x[i]- v[i])
575  + (x[i+1]-v[i+1]) * (x[i+1]-v[i+1])
576  + (x[i+2]-v[i+2]) * (x[i+2]-v[i+2])
577  + (x[i+3]-v[i+3]) * (x[i+3]-v[i+3])
578  + (x[i+4]-v[i+4]) * (x[i+4]-v[i+4]);
579  funVal2=temp/2;
580
581  temp=0;
582  m=nn%5;
583  if (m!=0){
584  for (i=0;i<m;i++)
585  temp+=fabs( x[i+1]-x[i] );
586  }
587  for (i=m;i<nn;i+=5)
588  temp=temp + fabs( x[i+1]-x[i] )
589  + fabs( x[i+2]-x[i+1] )
590  + fabs( x[i+3]-x[i+2] )
591  + fabs( x[i+4]-x[i+3] )
592  + fabs( x[i+5]-x[i+4] );
593  funVal2=funVal2 + lambda * temp;
594
595
596  /*
597  printf("\n funVal1=%e, funVal2=%e, diff=%e\n", funVal1, funVal2, funVal1-funVal2);
598  */
599
600
601
602
603  if (funVal2 > funVal1){ /*
604  we compute the solution by the first way
605  */
606  x[0]=v[0] + z[0];
607  for(i=1;i<n-1;i++)
608  x[i]= v[i] - z[i-1] + z[i];
609  x[n-1]=v[n-1] - z[n-2];
610  }
611  else{
612
613  /*
614  the solution x is computed in the second way
615  the gap can be further reduced
616  (note that, there might be numerical error)
617  */
618
619  *gap=*gap - (funVal1- funVal2);
620  if (*gap <0)
621  *gap=0;
622  }
623
624  return (numS);
625 }
626
627
628 void restartMapping(double *g, double *z, double * v,
629  double lambda, int nn)
630 {
631
632  int i, n=nn+1;
633  double temp;
634  int* S=(int *) malloc(sizeof(int)*nn);
635  double *x=(double *)malloc(sizeof(double)*n);
636  double *s=(double *)malloc(sizeof(double)*nn);
637  double *Av=(double *)malloc(sizeof(double)*nn);
638  //int numS=-1;
639
640  /*
641  for a given input z,
642  we compute the z0 after restarting
643
644  The elements in z lie in [-lambda, lambda]
645
646  The returned value is g
647  */
648
649
650  for (i=0;i<nn; i++)
651  Av[i]=v[i+1]-v[i];
652
653
654
655  g[0]=z[0] + z[0] - z[1] - Av[0];
656  for (i=1;i<nn-1;i++){
657  g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
658  }
659  g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
660
661
662  //numS = supportSet(x, v, z, g, S, lambda, nn);
663
664
665  /*With x, we compute z via
666  AA^T z = Av - Ax
667  */
668
669  /*
670  compute s= Av -Ax
671  */
672
673  for (i=0;i<nn; i++)
674  s[i]=Av[i] - x[i+1] + x[i];
675
676
677  /*
678  Apply Rose Algorithm for solving z
679  */
680
681  Thomas(&temp, g, s, nn);
682
683  /*
684  Rose(&temp, g, s, nn);
685  */
686
687  /*
688  project g to [-lambda, lambda]
689  */
690
691  for(i=0;i<nn;i++){
692  if (g[i]>lambda)
693  g[i]=lambda;
694  else
695  if (g[i]<-lambda)
696  g[i]=-lambda;
697  }
698
699
700  free (S);
701  free (x);
702  free (s);
703  free (Av);
704
705 }
706
707
708
709 /*
710
712
713 Our objective is to solve the fused Lasso signal approximator (flsa) problem:
714
715 min_x g(x) 1/2 \|x-v\|^2 + lambda \|A x\|_1, (1)
716
717 Let x* be the solution (which is unique), it satisfies
718
719 0 in x* - v + A^T * lambda *SGN(Ax*) (2)
720
721 To solve x*, it suffices to find
722
723 y* in A^T * lambda *SGN(Ax*) (3)
724 that satisfies
725
726 x* - v + y* =0 (4)
727 which leads to
728 x*= v - y* (5)
729
730 Due to the uniqueness of x*, we conclude that y* is unique.
731
732 As y* is a subgradient of lambda \|A x*\|_1,
733 we name our method as Subgradient Finding Algorithm (sfa).
734
735 y* in (3) can be further written as
736
737 y*= A^T * z* (6)
738 where
739
740 z* in lambda* SGN (Ax*) (7)
741
742 From (6), we have
743 z* = (A A^T)^{-1} A * y* (8)
744
745 Therefore, from the uqniueness of y*, we conclude that z* is also unique.
746 Next, we discuss how to solve this unique z*.
747
748 The problem (1) can reformulated as the following equivalent problem:
749
750 min_x max_z f(x, z)= 1/2 \|x-v\|^2 + <A x, z>
751 subject to \|z\|_{infty} \leq lambda (9)
752
753 At the saddle point, we have
754
755 x = v - AT z, (10)
756
757 which somehow concides with (5) and (6)
758
759 Plugging (10) into (9), we obtain the problem
760
761 min_z 1/2 z^T A AT z - < z, A v>,
762 subject to \|z\|_{infty} \leq lambda, (11)
763
764 In this program, we apply the Nesterov's method for solving (11).
765
766
767 Duality gap:
768
769 At a given point z0, we compute x0= v - A^T z0.
770 It is easy to show that
771 min_x f(x, z0) = f(x0, z0) <= max_z f(x0, z) (12)
772
773 Moreover, we have
774 max_z f(x0, z) - min_x f(x, z0)
775 <= lambda * \|A x0\|_1 - < z0, Av - A A^T z0> (13)
776
777 It is also to get that
778
779 f(x0, z0) <= f(x*, z*) <= max_z f(x0, z) (14)
780
781 g(x*)=f(x*, z*) (15)
782
783 g(x0)=max_z f(x0, z) (17)
784
785  Therefore, we have
786
787 g(x0)-g(x*) <= lambda * \|A x0\|_1 - < z0, Av - A A^T z0> (18)
788
789
790  We have applied a restarting technique, which is quite involved; and thus, we do not explain here.
791
793  */
794
795
796  /*
798
799  For sfa, the stepsize of the Nesterov's method is fixed to 1/4, so that no line search is needed.
800
801
802
803  Explanation of the parameters:
804
805  Output parameters
806  x: the solution to the primal problem
807  gap: the duality gap (pointer)
808
809  Input parameters
810  z: the solution to the dual problem (before calling this function, z contains a starting point)
811  !!!!we assume that the starting point has been successfully initialized in z !!!!
812  z0: a variable used for multiple purposes:
813  1) the previous solution z0
814  2) the difference between z and z0, i.e., z0=z- z0
815
816  lambda: the regularization parameter (and the radius of the infity ball, see (11)).
817  nn: the length of z, z0, Av, g, and s
818  maxStep: the maximal number of iterations
819
820  v: the point to be projected (not changed after the program)
821  Av: A*v (not changed after the program)
822
823  s: the search point (used for multiple purposes)
824  g: the gradient at g (and it is also used for multiple purposes)
825
826  tol: the tolerance of the gap
827  tau: the duality gap or the restarting technique is done every tau steps
828  flag: if flag=1, we apply the resart technique
829  flag=2, just run the SFA algorithm, terminate it when the absolution change is less than tol
830  flag=3, just run the SFA algorithm, terminate it when the duality gap is less than tol
831  flag=4, just run the SFA algorithm, terminate it when the relative duality gap is less than tol
832
833
834  We would like to emphasis that the following assumptions
835  have been checked in the functions that call this function:
836  1) 0< lambda < z_max
837  2) nn >=2
838  3) z has been initialized with a starting point
839  4) z0 has been initialized with all zeros
840
841  The termination condition is checked every tau iterations.
842
843  For the duality gap, please refer to (12-18)
844  */
845
846  int sfa(double *x, double *gap, int * activeS,
847  double *z, double *z0, double * v, double * Av,
848  double lambda, int nn, int maxStep,
849  double *s, double *g,
850  double tol, int tau, int flag){
851
852  int i, iterStep, m, tFlag=0, n=nn+1;
853  double alphap=0, alpha=1, beta=0, temp;
854  int* S=(int *) malloc(sizeof(int)*nn);
855  double gapp=-1, gappp=-1; /*gapp denotes the previous gap*/
856  int numS=-1, numSp=-2, numSpp=-3;;
857  /*
858  numS denotes the number of elements in the Support Set S
859  numSp denotes the number of elements in the previous Support Set S
860  */
861
862  *gap=-1; /*initial a value -1*/
863
864  /*
865  The main algorithm by Nesterov's method
866
867  B is an nn x nn tridiagonal matrix.
868
869  The nn eigenvalues of B are 2- 2 cos (i * PI/ n), i=1, 2, ..., nn
870  */
871
872  for (iterStep=1; iterStep<=maxStep; iterStep++){
873
874
875  /*------------- Step 1 ---------------------*/
876
877  beta=(alphap -1 ) / alpha;
878  /*
879  compute search point
880
881  s= z + beta * z0
882
883  We follow the style of CLAPACK
884  */
885  m=nn % 5;
886  if (m!=0){
887  for (i=0;i<m; i++)
888  s[i]=z[i]+ beta* z0[i];
889  }
890  for (i=m;i<nn;i+=5){
891  s[i] =z[i] + beta* z0[i];
892  s[i+1] =z[i+1] + beta* z0[i+1];
893  s[i+2] =z[i+2] + beta* z0[i+2];
894  s[i+3] =z[i+3] + beta* z0[i+3];
895  s[i+4] =z[i+4] + beta* z0[i+4];
896  }
897
898  /*
899  s and g are of size nn x 1
900
901  compute the gradient at s
902
903  g= B * s - Av,
904
905  where B is an nn x nn tridiagonal matrix. and is defined as
906
907  B= [ 2 -1 0 0;
908  -1 2 -1 0;
909  0 -1 2 -1;
910  0 0 -1 2]
911
912  We assume n>=3, which leads to nn>=2
913  */
914  g[0]=s[0] + s[0] - s[1] - Av[0];
915  for (i=1;i<nn-1;i++){
916  g[i]= - s[i-1] + s[i] + s[i] - s[i+1] - Av[i];
917  }
918  g[nn-1]= -s[nn-2] + s[nn-1] + s[nn-1] - Av[nn-1];
919
920
921  /*
922  z0 stores the previous -z
923  */
924  m=nn%7;
925  if (m!=0){
926  for (i=0;i<m;i++)
927  z0[i]=-z[i];
928  }
929  for (i=m; i <nn; i+=7){
930  z0[i] = - z[i];
931  z0[i+1] = - z[i+1];
932  z0[i+2] = - z[i+2];
933  z0[i+3] = - z[i+3];
934  z0[i+4] = - z[i+4];
935  z0[i+5] = - z[i+5];
936  z0[i+6] = - z[i+6];
937  }
938
939
940  /*
941  do a gradient step based on s to get z
942  */
943  m=nn%5;
944  if (m!=0){
945  for(i=0;i<m; i++)
946  z[i]=s[i] - g[i]/4;
947  }
948  for (i=m;i<nn; i+=5){
949  z[i] = s[i] - g[i] /4;
950  z[i+1] = s[i+1] - g[i+1]/4;
951  z[i+2] = s[i+2] - g[i+2]/4;
952  z[i+3] = s[i+3] - g[i+3]/4;
953  z[i+4] = s[i+4] - g[i+4]/4;
954  }
955
956  /*
957  project z onto the L_{infty} ball with radius lambda
958
959  z is the new approximate solution
960  */
961  for (i=0;i<nn; i++){
962  if (z[i]>lambda)
963  z[i]=lambda;
964  else
965  if (z[i]<-lambda)
966  z[i]=-lambda;
967  }
968
969  /*
970  compute the difference between the new solution
971  and the previous solution (stored in z0=-z_p)
972
973  the difference is written to z0
974  */
975
976  m=nn%5;
977  if (m!=0){
978  for (i=0;i<m;i++)
979  z0[i]+=z[i];
980  }
981  for(i=m;i<nn; i+=5){
982  z0[i] +=z[i];
983  z0[i+1]+=z[i+1];
984  z0[i+2]+=z[i+2];
985  z0[i+3]+=z[i+3];
986  z0[i+4]+=z[i+4];
987  }
988
989
990  alphap=alpha;
991  alpha=(1+sqrt(4*alpha*alpha+1))/2;
992
993  /*
994  check the termination condition
995  */
996  if (iterStep%tau==0){
997
998
999  /*
1000  The variables g and s can be modified
1001
1002  The variables x, z0 and z can be revised for case 0, but not for the rest
1003  */
1004  switch (flag){
1005  case 1:
1006
1007  /*
1008
1009  terminate the program once the "duality gap" is smaller than tol
1010
1011  compute the duality gap:
1012
1013  x= v - A^T z
1014  Ax = Av - A A^T z = -g,
1015  where
1016  g = A A^T z - A v
1017
1018
1019  the duality gap= lambda * \|Ax\|-1 - <z, Ax>
1020  = lambda * \|g\|_1 + <z, g>
1021
1022  In fact, gap=0 indicates that,
1023  if g_i >0, then z_i=-lambda
1024  if g_i <0, then z_i=lambda
1025  */
1026
1027  gappp=gapp;
1028  gapp=*gap; /*record the previous gap*/
1029  numSpp=numSp;
1030  numSp=numS; /*record the previous numS*/
1031
1032  dualityGap(gap, z, g, s, Av, lambda, nn);
1033  /*g is computed as the gradient of z in this function*/
1034
1035
1036  /*
1037  printf("\n Iteration: %d, gap=%e, numS=%d", iterStep, *gap, numS);
1038  */
1039
1040  /*
1041  If *gap <=tol, we terminate the iteration
1042  Otherwise, we restart the algorithm
1043  */
1044
1045  if (*gap <=tol){
1046  tFlag=1;
1047  break;
1048
1049  } /* end of *gap <=tol */
1050  else{
1051
1052  /* we apply the restarting technique*/
1053
1054  /*
1055  we compute the solution by the second way
1056  */
1057  numS = supportSet(x, v, z, g, S, lambda, nn);
1058  /*g, the gradient of z should be computed before calling this function*/
1059
1060  /*With x, we compute z via
1061  AA^T z = Av - Ax
1062  */
1063
1064  /*
1065  printf("\n iterStep=%d, numS=%d, gap=%e",iterStep, numS, *gap);
1066  */
1067
1068
1069  m=1;
1070  if (nn > 1000000)
1071  m=10;
1072  else
1073  if (nn > 100000)
1074  m=5;
1075
1076  if ( abs(numS-numSp) < m){
1077
1078  numS=generateSolution(x, z, gap, v, Av,
1079  g, s, S, lambda, nn);
1080  /*g, the gradient of z should be computed before calling this function*/
1081
1082
1083  if (*gap <tol){
1084  tFlag=2; /*tFlag =2 shows that the result is already optimal
1085  There is no need to call generateSolution for recomputing the best solution
1086  */
1087  break;
1088  }
1089
1090  if ( (*gap ==gappp) && (numS==numSpp) ){
1091
1092  tFlag=2;
1093  break;
1094
1095  }
1096
1097  /*we terminate the program is *gap <1
1098  numS==numSP
1099  and gapp==*gap
1100  */
1101  }
1102
1103  /*
1104  compute s= Av -Ax
1105  */
1106  for (i=0;i<nn; i++)
1107  s[i]=Av[i] - x[i+1] + x[i];
1108
1109  /*
1110  apply Rose Algorithm for solving z
1111  */
1112
1113  Thomas(&temp, z, s, nn);
1114
1115  /*
1116  Rose(&temp, z, s, nn);
1117  */
1118
1119  /*
1120  printf("\n Iteration: %d, %e", iterStep, temp);
1121  */
1122
1123  /*
1124  project z to [-lambda2, lambda2]
1125  */
1126  for(i=0;i<nn;i++){
1127  if (z[i]>lambda)
1128  z[i]=lambda;
1129  else
1130  if (z[i]<-lambda)
1131  z[i]=-lambda;
1132  }
1133
1134
1135
1136  m=nn%7;
1137  if (m!=0){
1138  for (i=0;i<m;i++)
1139  z0[i]=0;
1140  }
1141  for (i=m; i<nn; i+=7){
1142  z0[i] = z0[i+1]
1143  = z0[i+2]
1144  = z0[i+3]
1145  = z0[i+4]
1146  = z0[i+5]
1147  = z0[i+6]
1148  =0;
1149  }
1150
1151
1152  alphap=0; alpha=1;
1153
1154  /*
1155  we restart the algorithm
1156  */
1157
1158  }
1159
1160  break; /*break case 1*/
1161
1162  case 2:
1163
1164  /*
1165  The program is terminated either the summation of the absolution change (denoted by z0)
1166  of z (from the previous zp) is less than tol * nn,
1167  or the maximal number of iteration (maxStep) is achieved
1168 Note: tol indeed measures the averaged per element change.
1169 */
1170  temp=0;
1171  m=nn%5;
1172  if (m!=0){
1173  for(i=0;i<m;i++)
1174  temp+=fabs(z0[i]);
1175  }
1176  for(i=m;i<nn;i+=5)
1177  temp=temp + fabs(z0[i])
1178  + fabs(z0[i+1])
1179  + fabs(z0[i+2])
1180  + fabs(z0[i+3])
1181  + fabs(z0[i+4]);
1182  *gap=temp / nn;
1183
1184  if (*gap <=tol){
1185
1186  tFlag=1;
1187  }
1188
1189  break;
1190
1191  case 3:
1192
1193  /*
1194
1195  terminate the program once the "duality gap" is smaller than tol
1196
1197  compute the duality gap:
1198
1199  x= v - A^T z
1200  Ax = Av - A A^T z = -g,
1201  where
1202  g = A A^T z - A v
1203
1204
1205  the duality gap= lambda * \|Ax\|-1 - <z, Ax>
1206  = lambda * \|g\|_1 + <z, g>
1207
1208  In fact, gap=0 indicates that,
1209  if g_i >0, then z_i=-lambda
1210  if g_i <0, then z_i=lambda
1211  */
1212
1213
1214  g[0]=z[0] + z[0] - z[1] - Av[0];
1215  for (i=1;i<nn-1;i++){
1216  g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1217  }
1218
1219  g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1220
1221  for (i=0;i<nn;i++)
1222  if (g[i]>0)
1223  s[i]=lambda + z[i];
1224  else
1225  s[i]=-lambda + z[i];
1226
1227  temp=0;
1228  m=nn%5;
1229  if (m!=0){
1230  for(i=0;i<m;i++)
1231  temp+=s[i]*g[i];
1232  }
1233  for(i=m;i<nn;i+=5)
1234  temp=temp + s[i] *g[i]
1235  + s[i+1]*g[i+1]
1236  + s[i+2]*g[i+2]
1237  + s[i+3]*g[i+3]
1238  + s[i+4]*g[i+4];
1239  *gap=temp;
1240
1241  /*
1242  printf("\n %e", *gap);
1243  */
1244
1245
1246  if (*gap <=tol)
1247  tFlag=1;
1248
1249  break;
1250
1251  case 4:
1252
1253  /*
1254
1255  terminate the program once the "relative duality gap" is smaller than tol
1256
1257
1258  compute the duality gap:
1259
1260  x= v - A^T z
1261  Ax = Av - A A^T z = -g,
1262  where
1263  g = A A^T z - A v
1264
1265
1266  the duality gap= lambda * \|Ax\|-1 - <z, Ax>
1267  = lambda * \|g\|_1 + <z, g>
1268
1269  In fact, gap=0 indicates that,
1270  if g_i >0, then z_i=-lambda
1271  if g_i <0, then z_i=lambda
1272
1273
1274  Here, the "relative duality gap" is defined as:
1275  duality gap / - 1/2 \|A^T z\|^2 + < z, Av>
1276
1277  We efficiently compute - 1/2 \|A^T z\|^2 + < z, Av> using the following relationship
1278
1279  - 1/2 \|A^T z\|^2 + < z, Av>
1280  = -1/2 <z, A A^T z - Av -Av>
1281  = -1/2 <z, g - Av>
1282  */
1283
1284
1285  g[0]=z[0] + z[0] - z[1] - Av[0];
1286  for (i=1;i<nn-1;i++){
1287  g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1288  }
1289
1290  g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1291
1292  for (i=0;i<nn;i++)
1293  if (g[i]>0)
1294  s[i]=lambda + z[i];
1295  else
1296  s[i]=-lambda + z[i];
1297
1298  temp=0;
1299  m=nn%5;
1300  if (m!=0){
1301  for(i=0;i<m;i++)
1302  temp+=s[i]*g[i];
1303  }
1304  for(i=m;i<nn;i+=5)
1305  temp=temp + s[i] *g[i]
1306  + s[i+1]*g[i+1]
1307  + s[i+2]*g[i+2]
1308  + s[i+3]*g[i+3]
1309  + s[i+4]*g[i+4];
1310  *gap=temp;
1311  /*
1312  Now, *gap contains the duality gap
1313  Next, we compute
1314  - 1/2 \|A^T z\|^2 + < z, Av>
1315  =-1/2 <z, g - Av>
1316  */
1317
1318  temp=0;
1319  m=nn%5;
1320  if (m!=0){
1321  for(i=0;i<m;i++)
1322  temp+=z[i] * (g[i] - Av[i]);
1323  }
1324  for(i=m;i<nn;i+=5)
1325  temp=temp + z[i] * (g[i] - Av[i])
1326  + z[i+1]* (g[i+1]- Av[i+1])
1327  + z[i+2]* (g[i+2]- Av[i+2])
1328  + z[i+3]* (g[i+3]- Av[i+3])
1329  + z[i+4]* (g[i+4]- Av[i+4]);
1330  temp=fabs(temp) /2;
1331
1332  if (temp <1)
1333  temp=1;
1334
1335  *gap/=temp;
1336  /*
1337  *gap now contains the relative gap
1338  */
1339
1340
1341  if (*gap <=tol){
1342  tFlag=1;
1343  }
1344
1345  break;
1346
1347  default:
1348
1349  /*
1350  The program is terminated either the summation of the absolution change (denoted by z0)
1351  of z (from the previous zp) is less than tol * nn,
1352  or the maximal number of iteration (maxStep) is achieved
1353 Note: tol indeed measures the averaged per element change.
1354 */
1355  temp=0;
1356  m=nn%5;
1357  if (m!=0){
1358  for(i=0;i<m;i++)
1359  temp+=fabs(z0[i]);
1360  }
1361  for(i=m;i<nn;i+=5)
1362  temp=temp + fabs(z0[i])
1363  + fabs(z0[i+1])
1364  + fabs(z0[i+2])
1365  + fabs(z0[i+3])
1366  + fabs(z0[i+4]);
1367  *gap=temp / nn;
1368
1369  if (*gap <=tol){
1370
1371  tFlag=1;
1372  }
1373
1374  break;
1375
1376  }/*end of switch*/
1377
1378  if (tFlag)
1379  break;
1380
1381  }/* end of the if for checking the termination condition */
1382
1383  /*-------------- Step 3 --------------------*/
1384
1385  }
1386
1387  /*
1388  for the other cases, except flag=1, compute the solution x according the first way (the primal-dual way)
1389  */
1390
1391  if ( (flag !=1) || (tFlag==0) ){
1392  x[0]=v[0] + z[0];
1393  for(i=1;i<n-1;i++)
1394  x[i]= v[i] - z[i-1] + z[i];
1395  x[n-1]=v[n-1] - z[n-2];
1396  }
1397
1398  if ( (flag==1) && (tFlag==1)){
1399
1400  /*
1401  We assume that n>=3, and thus nn>=2
1402
1403  We have two ways for recovering x.
1404  The first way is x = v - A^T z
1405  The second way is x =omega(z)
1406  */
1407
1408  /*
1409  We first compute the objective function value of the first choice in terms f(x), see our paper
1410  */
1411
1412  /*
1413  for numerical reason, we do a gradient descent step
1414  */
1415
1416  /*
1417  ---------------------------------------------------
1418  A gradient step begins
1419  */
1420  g[0]=z[0] + z[0] - z[1] - Av[0];
1421  for (i=1;i<nn-1;i++){
1422  g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1423  }
1424  g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1425
1426
1427  /*
1428  do a gradient step based on z to get the new z
1429  */
1430  m=nn%5;
1431  if (m!=0){
1432  for(i=0;i<m; i++)
1433  z[i]=z[i] - g[i]/4;
1434  }
1435  for (i=m;i<nn; i+=5){
1436  z[i] = z[i] - g[i] /4;
1437  z[i+1] = z[i+1] - g[i+1]/4;
1438  z[i+2] = z[i+2] - g[i+2]/4;
1439  z[i+3] = z[i+3] - g[i+3]/4;
1440  z[i+4] = z[i+4] - g[i+4]/4;
1441  }
1442
1443  /*
1444  project z onto the L_{infty} ball with radius lambda
1445
1446  z is the new approximate solution
1447  */
1448  for (i=0;i<nn; i++){
1449  if (z[i]>lambda)
1450  z[i]=lambda;
1451  else
1452  if (z[i]<-lambda)
1453  z[i]=-lambda;
1454  }
1455
1456  /*
1457  ---------------------------------------------------
1458  A gradient descent step ends
1459  */
1460
1461  /*compute the gradient at z*/
1462
1463  g[0]=z[0] + z[0] - z[1] - Av[0];
1464  for (i=1;i<nn-1;i++){
1465  g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1466  }
1467  g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1468
1469
1470  numS=generateSolution(x, z, gap, v, Av,
1471  g, s, S, lambda, nn);
1472  /*g, the gradient of z should be computed before calling this function*/
1473
1474  }
1475
1476  free (S);
1477  /*
1478  free the variables S
1479  */
1480
1481  *activeS=numS;
1482  return (iterStep);
1483
1484  }
1485
1486
1487 /*
1488
1489  Refer to sfa for the defintions of the variables
1490
1491  In this file, we restart the program every step, and neglect the gradient step.
1492
1493  It seems that, this program does not converge.
1494
1495  This function shows that the gradient step is necessary.
1496  */
1497
1498 int sfa_special(double *x, double *gap, int * activeS,
1499  double *z, double * v, double * Av,
1500  double lambda, int nn, int maxStep,
1501  double *s, double *g,
1502  double tol, int tau){
1503
1504  int i, iterStep;
1505  //int tFlag=0;
1506  //int n=nn+1;
1507  double temp;
1508  int* S=(int *) malloc(sizeof(int)*nn);
1509  double gapp=-1; /*gapp denotes the previous gap*/
1510  int numS=-1, numSp=-1;
1511  /*
1512  numS denotes the number of elements in the Support Set S
1513  numSp denotes the number of elements in the previous Support Set S
1514  */
1515
1516  *gap=-1; /*initialize *gap a value*/
1517
1518  for (iterStep=1; iterStep<=maxStep; iterStep++){
1519
1520
1521  g[0]=z[0] + z[0] - z[1] - Av[0];
1522  for (i=1;i<nn-1;i++){
1523  g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1524  }
1525  g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1526
1527  numSp=numS; /*record the previous numS*/
1528  numS = supportSet(x, v, z, g, S, lambda, nn);
1529
1530
1531  /*With x, we compute z via
1532  AA^T z = Av - Ax
1533  */
1534
1535  /*
1536  compute s= Av -Ax
1537  */
1538
1539  for (i=0;i<nn; i++)
1540  s[i]=Av[i] - x[i+1] + x[i];
1541
1542
1543  /*
1544  Apply Rose Algorithm for solving z
1545  */
1546
1547  Thomas(&temp, z, s, nn);
1548
1549  /*
1550  Rose(&temp, z, s, nn);
1551  */
1552
1553  /*
1554  project z to [-lambda, lambda]
1555  */
1556
1557  for(i=0;i<nn;i++){
1558  if (z[i]>lambda)
1559  z[i]=lambda;
1560  else
1561  if (z[i]<-lambda)
1562  z[i]=-lambda;
1563  }
1564
1565
1566  if (iterStep%tau==0){
1567  gapp=*gap; /*record the previous gap*/
1568
1569  dualityGap(gap, z, g, s, Av, lambda, nn);
1570
1571  /*
1572  printf("\n iterStep=%d, numS=%d, gap=%e, diff=%e",iterStep, numS, *gap, *gap -gapp);
1573
1574 */
1575
1576  if (*gap <=tol){
1577  //tFlag=1;
1578  break;
1579  }
1580
1581  if ( (*gap <1) && (numS==numSp) && fabs(gapp == *gap) ){
1582  //tFlag=1;
1583  break;
1584  /*we terminate the program is *gap <1
1585  numS==numSP
1586  and gapp==*gap
1587  */
1588  }
1589
1590  }/*end of if tau*/
1591
1592  }/*end for */
1593
1594  free (S);
1595
1596  * activeS=numS;
1597  return(iterStep);
1598
1599 }
1600
1601
1602 /*
1603
1604  We do one gradient descent, and then restart the program
1605  */
1606
1607
1608 int sfa_one(double *x, double *gap, int * activeS,
1609  double *z, double * v, double * Av,
1610  double lambda, int nn, int maxStep,
1611  double *s, double *g,
1612  double tol, int tau){
1613
1614  int i, iterStep, m;
1615  int tFlag=0;
1616  //int n=nn+1;
1617  double temp;
1618  int* S=(int *) malloc(sizeof(int)*nn);
1619  double gapp=-1, gappp=-2; /*gapp denotes the previous gap*/
1620  int numS=-100, numSp=-200, numSpp=-300;
1621  /*
1622  numS denotes the number of elements in the Support Set S
1623  numSp denotes the number of elements in the previous Support Set S
1624  */
1625
1626  *gap=-1; /*initialize *gap a value*/
1627
1628  /*
1629  The main algorithm by Nesterov's method
1630
1631  B is an nn x nn tridiagonal matrix.
1632
1633  The nn eigenvalues of B are 2- 2 cos (i * PI/ n), i=1, 2, ..., nn
1634  */
1635
1636
1637  /*
1638  we first do a gradient step based on z
1639  */
1640
1641
1642  /*
1643  ---------------------------------------------------
1644  A gradient step begins
1645  */
1646  g[0]=z[0] + z[0] - z[1] - Av[0];
1647  for (i=1;i<nn-1;i++){
1648  g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1649  }
1650  g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1651
1652
1653  /*
1654  do a gradient step based on z to get the new z
1655  */
1656  m=nn%5;
1657  if (m!=0){
1658  for(i=0;i<m; i++)
1659  z[i]=z[i] - g[i]/4;
1660  }
1661  for (i=m;i<nn; i+=5){
1662  z[i] = z[i] - g[i] /4;
1663  z[i+1] = z[i+1] - g[i+1]/4;
1664  z[i+2] = z[i+2] - g[i+2]/4;
1665  z[i+3] = z[i+3] - g[i+3]/4;
1666  z[i+4] = z[i+4] - g[i+4]/4;
1667  }
1668
1669  /*
1670  project z onto the L_{infty} ball with radius lambda
1671
1672  z is the new approximate solution
1673  */
1674  for (i=0;i<nn; i++){
1675  if (z[i]>lambda)
1676  z[i]=lambda;
1677  else
1678  if (z[i]<-lambda)
1679  z[i]=-lambda;
1680  }
1681
1682  /*
1683  ---------------------------------------------------
1684  A gradient descent step ends
1685  */
1686
1687
1688  /*compute the gradient at z*/
1689
1690  g[0]=z[0] + z[0] - z[1] - Av[0];
1691  for (i=1;i<nn-1;i++){
1692  g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1693  }
1694  g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1695
1696  for (iterStep=1; iterStep<=maxStep; iterStep++){
1697
1698
1699  /*
1700  ---------------------------------------------------
1701  restart the algorithm with x=omega(z)
1702  */
1703
1704  numSpp=numSp;
1705  numSp=numS; /*record the previous numS*/
1706  numS = supportSet(x, v, z, g, S, lambda, nn);
1707
1708
1709  /*With x, we compute z via
1710  AA^T z = Av - Ax
1711  */
1712
1713  /*
1714  compute s= Av -Ax
1715  */
1716
1717  for (i=0;i<nn; i++)
1718  s[i]=Av[i] - x[i+1] + x[i];
1719
1720
1721  /*
1722  Apply Rose Algorithm for solving z
1723  */
1724
1725  Thomas(&temp, z, s, nn);
1726
1727  /*
1728  Rose(&temp, z, s, nn);
1729  */
1730
1731  /*
1732  project z to [-lambda, lambda]
1733  */
1734
1735  for(i=0;i<nn;i++){
1736  if (z[i]>lambda)
1737  z[i]=lambda;
1738  else
1739  if (z[i]<-lambda)
1740  z[i]=-lambda;
1741  }
1742
1743  /*
1744  ---------------------------------------------------
1745  restart the algorithm with x=omega(z)
1746
1747  we have computed a new z, based on the above relationship
1748  */
1749
1750
1751  /*
1752  ---------------------------------------------------
1753  A gradient step begins
1754  */
1755  g[0]=z[0] + z[0] - z[1] - Av[0];
1756  for (i=1;i<nn-1;i++){
1757  g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1758  }
1759  g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1760
1761
1762  /*
1763  do a gradient step based on z to get the new z
1764  */
1765  m=nn%5;
1766  if (m!=0){
1767  for(i=0;i<m; i++)
1768  z[i]=z[i] - g[i]/4;
1769  }
1770  for (i=m;i<nn; i+=5){
1771  z[i] = z[i] - g[i] /4;
1772  z[i+1] = z[i+1] - g[i+1]/4;
1773  z[i+2] = z[i+2] - g[i+2]/4;
1774  z[i+3] = z[i+3] - g[i+3]/4;
1775  z[i+4] = z[i+4] - g[i+4]/4;
1776  }
1777
1778  /*
1779  project z onto the L_{infty} ball with radius lambda
1780
1781  z is the new approximate solution
1782  */
1783  for (i=0;i<nn; i++){
1784  if (z[i]>lambda)
1785  z[i]=lambda;
1786  else
1787  if (z[i]<-lambda)
1788  z[i]=-lambda;
1789  }
1790
1791  /*
1792  ---------------------------------------------------
1793  A gradient descent step ends
1794  */
1795
1796  /*compute the gradient at z*/
1797
1798  g[0]=z[0] + z[0] - z[1] - Av[0];
1799  for (i=1;i<nn-1;i++){
1800  g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i];
1801  }
1802  g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1];
1803
1804
1805  if (iterStep % tau==0){
1806  gappp=gapp;
1807  gapp=*gap; /*record the previous gap*/
1808
1809  dualityGap2(gap, z, g, s, Av, lambda, nn);
1810  /*g, the gradient of z should be computed before calling this function*/
1811
1812
1813  /*
1814  printf("\n iterStep=%d, numS=%d, gap=%e",iterStep, numS, *gap);
1815  */
1816
1817
1818  /*
1819  printf("\n %d & %d & %2.0e \\\\ \n \\hline ",iterStep, numS, *gap);
1820  */
1821
1822
1823  /*
1824  printf("\n %e",*gap);
1825  */
1826
1827  /*
1828
1829  printf("\n %d",numS);
1830
1831 */
1832
1833  if (*gap <=tol){
1834  //tFlag=1;
1835  break;
1836  }
1837
1838  m=1;
1839  if (nn > 1000000)
1840  m=5;
1841  else
1842  if (nn > 100000)
1843  m=3;
1844
1845  if ( abs( numS-numSp) <m ){
1846
1847  /*
1848  printf("\n numS=%d, numSp=%d",numS,numSp);
1849  */
1850
1851  m=generateSolution(x, z, gap, v, Av,
1852  g, s, S, lambda, nn);
1853  /*g, the gradient of z should be computed before calling this function*/
1854
1855  if (*gap < tol){
1856
1857  numS=m;
1858  tFlag=2;
1859  break;
1860  }
1861
1862
1863  if ( (*gap ==gappp) && (numS==numSpp) ){
1864
1865  tFlag=2;
1866  break;
1867
1868  }
1869
1870  } /*end of if*/
1871
1872  }/*end of if tau*/
1873
1874
1875  } /*end of for*/
1876
1877
1878
1879  if (tFlag!=2){
1880  numS=generateSolution(x, z, gap, v, Av, g, s, S, lambda, nn);
1881  /*g, the gradient of z should be computed before calling this function*/
1882  }
1883
1884  free(S);
1885
1886  *activeS=numS;
1887  return(iterStep);
1888 }

SHOGUN Machine Learning Toolbox - Documentation