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