SHOGUN  4.1.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 #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) && (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 }
sfa_one
int sfa_one(double *x, double *gap, int *activeS, double *z, double *v, double *Av, double lambda, int nn, int maxStep, double *s, double *g, double tol, int tau)
Definition: sfa.cpp:1608
Rose
void Rose(double *zMax, double *z0, double *Av, int nn)
Definition: sfa.cpp:223
generateSolution
int generateSolution(double *x, double *z, double *gap, double *v, double *Av, double *g, double *s, int *S, double lambda, int nn)
Definition: sfa.cpp:506
restartMapping
void restartMapping(double *g, double *z, double *v, double lambda, int nn)
Definition: sfa.cpp:628
Thomas
void Thomas(double *zMax, double *z0, double *Av, int nn)
Definition: sfa.cpp:152
sfa_special
int sfa_special(double *x, double *gap, int *activeS, double *z, double *v, double *Av, double lambda, int nn, int maxStep, double *s, double *g, double tol, int tau)
Definition: sfa.cpp:1498
dualityGap2
void dualityGap2(double *gap, double *z, double *g, double *s, double *Av, double lambda, int nn)
Definition: sfa.cpp:458
dualityGap
void dualityGap(double *gap, double *z, double *g, double *s, double *Av, double lambda, int nn)
Definition: sfa.cpp:414
sfa
int sfa(double *x, double *gap, int *activeS, double *z, double *z0, double *v, double *Av, double lambda, int nn, int maxStep, double *s, double *g, double tol, int tau, int flag)
Definition: sfa.cpp:846
supportSet
int supportSet(double *x, double *v, double *z, double *g, int *S, double lambda, int nn)
Definition: sfa.cpp:306
delta
#define delta
Definition: sfa.cpp:23
SHOGUN Machine Learning Toolbox - Documentation