SHOGUN
v2.0.0
|
00001 /* This program is free software: you can redistribute it and/or modify 00002 * it under the terms of the GNU General Public License as published by 00003 * the Free Software Foundation, either version 3 of the License, or 00004 * (at your option) any later version. 00005 * 00006 * This program is distributed in the hope that it will be useful, 00007 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00008 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00009 * GNU General Public License for more details. 00010 * 00011 * You should have received a copy of the GNU General Public License 00012 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00013 * 00014 * Copyright (C) 2009 - 2012 Jun Liu and Jieping Ye 00015 */ 00016 00017 #include <shogun/lib/slep/flsa/sfa.h> 00018 #include <stdlib.h> 00019 #include <stdio.h> 00020 #include <time.h> 00021 #include <math.h> 00022 00023 #define delta 1e-10 00024 00025 /* 00026 Revision History 00027 00028 First Version available on October 10, 2009 00029 00030 A runnable version on October 15, 2009 00031 00032 Major revision on October 29, 2009 00033 (Some functions appearing in a previous version have deleted, please refer to the previous version for the old functions. 00034 Some new functions have been added as well) 00035 00036 */ 00037 00038 /* 00039 00040 Files contained in this header file sfa.h: 00041 00042 1. Algorithms for solving the linear system A A^T z0 = Av (see the description of A from the following context) 00043 00044 void Thomas(double *zMax, double *z0, 00045 double * Av, int nn) 00046 00047 void Rose(double *zMax, double *z0, 00048 double * Av, int nn) 00049 00050 int supportSet(double *x, double *v, double *z, 00051 double *g, int * S, double lambda, int nn) 00052 00053 void dualityGap(double *gap, double *z, 00054 double *g, double *s, double *Av, 00055 double lambda, int nn) 00056 00057 void dualityGap2(double *gap, double *z, 00058 double *g, double *s, double *Av, 00059 double lambda, int nn) 00060 00061 00062 2. The Subgraident Finding Algorithm (SFA) for solving problem (4) (refer to the description of the problem for detail) 00063 00064 int sfa(double *x, double *gap, 00065 double *z, double *z0, double * v, double * Av, 00066 double lambda, int nn, int maxStep, 00067 double *s, double *g, 00068 double tol, int tau, int flag) 00069 00070 int sfa_special(double *x, double *gap, 00071 double *z, double * v, double * Av, 00072 double lambda, int nn, int maxStep, 00073 double *s, double *g, 00074 double tol, int tau) 00075 00076 int sfa_one(double *x, double *gap, 00077 double *z, double * v, double * Av, 00078 double lambda, int nn, int maxStep, 00079 double *s, double *g, 00080 double tol, int tau) 00081 00082 00083 */ 00084 00085 00086 /* 00087 00088 Some mathematical background. 00089 00090 In this file, we discuss how to solve the following subproblem, 00091 00092 min_x 1/2 \|x-v\|^2 + lambda \|A x\|_1, (1) 00093 00094 which is a key problem used in the Fused Lasso Signal Approximator (FLSA). 00095 00096 Also, note that, FLSA is a building block for solving the optimation problmes with fused Lasso penalty. 00097 00098 In (1), x and v are n-dimensional vectors, 00099 and A is a matrix with size (n-1) x n, and is defined as follows (e.g., n=4): 00100 A= [ -1 1 0 0; 00101 0 -1 1 0; 00102 0 0 -1 1] 00103 00104 The above problem can be reformulated as the following equivalent min-max optimization problem 00105 00106 min_x max_z 1/2 \|x-v\|^2 + <A x, z> 00107 subject to \|z\|_{infty} \leq lambda (2) 00108 00109 00110 It is easy to get that, at the optimal point 00111 00112 x = v - AT z, (3) 00113 00114 where z is the optimal solution to the following optimization problem 00115 00116 min_z 1/2 z^T A AT z - < z, A v>, 00117 subject to \|z\|_{infty} \leq lambda (4) 00118 00119 00120 00121 Let B=A A^T. It is easy to get that B is a (n-1) x (n-1) tridiagonal matrix. 00122 When n=5, B is defined as: 00123 B= [ 2 -1 0 0; 00124 -1 2 -1 0; 00125 0 -1 2 -1; 00126 0 0 -1 2] 00127 00128 Let z0 be the solution to the linear system: 00129 00130 A A^T * z0 = A * v (5) 00131 00132 The problem (5) can be solve by the Thomas Algorithm, in about 5n multiplications and 4n additions. 00133 00134 It can also be solved by the Rose's Algorithm, in about 2n multiplications and 2n additions. 00135 00136 Moreover, considering the special structure of the matrix A (and B), 00137 it can be solved in about n multiplications and 3n additions 00138 00139 If lambda \geq \|z0\|_{infty}, x_i= mean(v), for all i, 00140 the problem (1) admits near analytical solution 00141 00142 00143 We have also added the restart technique, please refer to our paper for detail! 00144 00145 */ 00146 00147 00148 /* 00150 */ 00151 00152 void Thomas(double *zMax, double *z0, double * Av, int nn){ 00153 00154 /* 00155 00156 We apply the Tomas algorithm for solving the following linear system 00157 B * z0 = Av 00158 Thomas algorithm is also called the tridiagonal matrix algorithm 00159 00160 B=[ 2 -1 0 0; 00161 -1 2 -1 0; 00162 0 -1 2 -1; 00163 0 0 -1 2] 00164 00165 z0 is the result, Av is unchanged after the computation 00166 00167 00168 c is a precomputed nn dimensional vector 00169 c=[-1/2, -2/3, -3/4, -4/5, ..., -nn/(nn+1)] 00170 00171 c[i]=- (i+1) / (i+2) 00172 c[i-1]=- i / (i+1) 00173 00174 z0 is an nn dimensional vector 00175 00176 */ 00177 00178 int i; 00179 double tt, z_max; 00180 00181 /* 00182 Modify the coefficients in Av (copy to z0) 00183 */ 00184 z0[0]=Av[0]/2; 00185 for (i=1;i < nn; i++){ 00186 tt=Av[i] + z0[i-1]; 00187 z0[i]=tt - tt / (i+2); 00188 } 00189 00190 /*z0[i]=(Av[i] + z0[i-1]) * (i+1) / (i+2);*/ 00191 00192 /*z0[i]=(Av[i] + z0[i-1])/ ( 2 - i / (i+1));*/ 00193 00194 00195 /* 00196 Back substitute (obtain the result in z0) 00197 */ 00198 z_max= fabs(z0[nn-1]); 00199 00200 for (i=nn-2; i>=0; i--){ 00201 00202 z0[i]+= z0[i+1] - z0[i+1]/ (i+2); 00203 00204 /*z0[i]+= z0[i+1] * (i+1) / (i+2);*/ 00205 00206 tt=fabs(z0[i]); 00207 00208 if (tt > z_max) 00209 z_max=tt; 00210 00211 } 00212 *zMax=z_max; 00213 00214 } 00215 00216 00217 00218 00219 /* 00221 */ 00222 00223 void Rose(double *zMax, double *z0, double * Av, int nn){ 00224 00225 /* 00226 We use the Rose algorithm for solving the following linear system 00227 B * z0 = Av 00228 00229 00230 B=[ 2 -1 0 0; 00231 -1 2 -1 0; 00232 0 -1 2 -1; 00233 0 0 -1 2] 00234 00235 z0 is the result, Av is unchanged after the computation 00236 00237 z0 is an nn dimensional vector 00238 00239 */ 00240 00241 int i, m; 00242 double s=0, z_max; 00243 00244 00245 /* 00246 We follow the style in CLAPACK 00247 */ 00248 m= nn % 5; 00249 if (m!=0){ 00250 for (i=0;i<m; i++) 00251 s+=Av[i] * (i+1); 00252 } 00253 for(i=m;i<nn;i+=5) 00254 s+= Av[i] * (i+1) 00255 + Av[i+1] * (i+2) 00256 + Av[i+2] * (i+3) 00257 + Av[i+3] * (i+4) 00258 + Av[i+4] * (i+5); 00259 s/=(nn+1); 00260 00261 00262 /* 00263 from nn-1 to 0 00264 */ 00265 z0[nn-1]=Av[nn-1]- s; 00266 for (i=nn-2;i >=0; i--){ 00267 z0[i]=Av[i] + z0[i+1]; 00268 } 00269 00270 /* 00271 from 0 to nn-1 00272 */ 00273 z_max= fabs(z0[0]); 00274 for (i=0; i<nn; i++){ 00275 00276 z0[i]+= z0[i-1]; 00277 00278 s=fabs(z0[i]); 00279 00280 if (s > z_max) 00281 z_max=s; 00282 00283 } 00284 *zMax=z_max; 00285 00286 } 00287 00288 00289 00290 /* 00292 00293 x=omega(z) 00294 00295 v: the vector to be projected 00296 z: the approximate solution 00297 g: the gradient at z (g should be computed before calling this function 00298 00299 nn: the length of z, g, and S (maximal length for S) 00300 00301 n: the length of x and v 00302 00303 S: records the indices of the elements in the support set 00304 */ 00305 00306 int supportSet(double *x, double *v, double *z, double *g, int * S, double lambda, int nn){ 00307 00308 int i, j, n=nn+1, numS=0; 00309 double temp; 00310 00311 00312 /* 00313 we first scan z and g to obtain the support set S 00314 */ 00315 00316 /*numS: number of the elements in the support set S*/ 00317 for(i=0;i<nn; i++){ 00318 if ( ( (z[i]==lambda) && (g[i] < delta) ) || ( (z[i]==-lambda) && (g[i] >delta) )){ 00319 S[numS]=i; 00320 numS++; 00321 } 00322 } 00323 00324 /* 00325 printf("\n %d",numS); 00326 */ 00327 00328 if (numS==0){ /*this shows that S is empty*/ 00329 temp=0; 00330 for (i=0;i<n;i++) 00331 temp+=v[i]; 00332 00333 temp=temp/n; 00334 for(i=0;i<n;i++) 00335 x[i]=temp; 00336 00337 return numS; 00338 } 00339 00340 00341 /* 00342 Next, we deal with numS >=1 00343 */ 00344 00345 /*process the first block 00346 00347 j=0 00348 */ 00349 temp=0; 00350 for (i=0;i<=S[0]; i++) 00351 temp+=v[i]; 00352 /*temp =sum (v [0: s[0] ]*/ 00353 temp=( temp + z[ S[0] ] ) / (S[0] +1); 00354 for (i=0;i<=S[0]; i++) 00355 x[i]=temp; 00356 00357 00358 /*process the middle blocks 00359 00360 If numS=1, it belongs the last block 00361 */ 00362 for (j=1; j < numS; j++){ 00363 temp=0; 00364 for (i= S[j-1] +1; i<= S[j]; i++){ 00365 temp+=v[i]; 00366 } 00367 00368 /*temp =sum (v [ S[j-1] +1: s[j] ]*/ 00369 00370 temp=(temp - z[ S[j-1] ] + z[ S[j] ])/ (S[j]- S[j-1]); 00371 00372 for (i= S[j-1] +1; i<= S[j]; i++){ 00373 x[i]=temp; 00374 } 00375 } 00376 00377 /*process the last block 00378 j=numS-1; 00379 */ 00380 temp=0; 00381 for (i=S[numS-1] +1 ;i< n; i++) 00382 temp+=v[i]; 00383 /*temp =sum (v [ (S[numS-1] +1): (n-1) ]*/ 00384 00385 temp=( temp - z[ S[numS-1] ] ) / (nn - S[numS-1]); /*S[numS-1] <= nn-1*/ 00386 00387 for (i=S[numS-1] +1 ;i< n; i++) 00388 x[i]=temp; 00389 00390 return numS; 00391 00392 } 00393 00394 00395 00396 /* 00397 00399 00400 we compute the duality corresponding the solution z 00401 00402 z: the approximate solution 00403 g: the gradient at z (we recompute the gradient) 00404 s: an auxiliary variable 00405 Av: A*v 00406 00407 nn: the lenght for z, g, s, and Av 00408 00409 The variables g and s shall be revised. 00410 00411 The variables z and Av remain unchanged. 00412 */ 00413 00414 void dualityGap(double *gap, double *z, double *g, double *s, double *Av, double lambda, int nn){ 00415 00416 int i, m; 00417 double temp; 00418 00419 00420 g[0]=z[0] + z[0] - z[1] - Av[0]; 00421 for (i=1;i<nn-1;i++){ 00422 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i]; 00423 } 00424 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1]; 00425 00426 00427 for (i=0;i<nn;i++) 00428 if (g[i]>0) 00429 s[i]=lambda + z[i]; 00430 else 00431 s[i]=-lambda + z[i]; 00432 00433 00434 temp=0; 00435 m=nn%5; 00436 00437 if (m!=0){ 00438 for(i=0;i<m;i++) 00439 temp+=s[i]*g[i]; 00440 } 00441 00442 for(i=m;i<nn;i+=5) 00443 temp=temp + s[i] *g[i] 00444 + s[i+1]*g[i+1] 00445 + s[i+2]*g[i+2] 00446 + s[i+3]*g[i+3] 00447 + s[i+4]*g[i+4]; 00448 *gap=temp; 00449 } 00450 00451 00452 /* 00453 Similar to dualityGap, 00454 00455 The difference is that, we assume that g has been computed. 00456 */ 00457 00458 void dualityGap2(double *gap, double *z, double *g, double *s, double *Av, double lambda, int nn){ 00459 00460 int i, m; 00461 double temp; 00462 00463 00464 /* 00465 g[0]=z[0] + z[0] - z[1] - Av[0]; 00466 for (i=1;i<nn-1;i++){ 00467 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i]; 00468 } 00469 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1]; 00470 00471 */ 00472 00473 for (i=0;i<nn;i++) 00474 if (g[i]>0) 00475 s[i]=lambda + z[i]; 00476 else 00477 s[i]=-lambda + z[i]; 00478 00479 00480 temp=0; 00481 m=nn%5; 00482 00483 if (m!=0){ 00484 for(i=0;i<m;i++) 00485 temp+=s[i]*g[i]; 00486 } 00487 00488 for(i=m;i<nn;i+=5) 00489 temp=temp + s[i] *g[i] 00490 + s[i+1]*g[i+1] 00491 + s[i+2]*g[i+2] 00492 + s[i+3]*g[i+3] 00493 + s[i+4]*g[i+4]; 00494 *gap=temp; 00495 } 00496 00497 00498 /* 00499 generateSolution: 00500 00501 generate the solution x based on the information of z and g 00502 (!!!!we assume that g has been computed as the gradient of z!!!!) 00503 00504 */ 00505 00506 int generateSolution(double *x, double *z, double *gap, 00507 double *v, double *Av, 00508 double *g, double *s, int *S, 00509 double lambda, int nn){ 00510 00511 int i, m, numS, n=nn+1; 00512 double temp, funVal1, funVal2; 00513 00514 /* 00515 z is the appropriate solution, 00516 and g contains its gradient 00517 */ 00518 00519 00520 /* 00521 We assume that n>=3, and thus nn>=2 00522 00523 We have two ways for recovering x. 00524 The first way is x = v - A^T z 00525 The second way is x =omega(z) 00526 */ 00527 00528 temp=0; 00529 m=nn%5; 00530 if (m!=0){ 00531 for (i=0;i<m;i++) 00532 temp+=z[i]*(g[i] + Av[i]); 00533 } 00534 for (i=m;i<nn;i+=5) 00535 temp=temp + z[i] *(g[i] + Av[i]) 00536 + z[i+1]*(g[i+1] + Av[i+1]) 00537 + z[i+2]*(g[i+2] + Av[i+2]) 00538 + z[i+3]*(g[i+3] + Av[i+3]) 00539 + z[i+4]*(g[i+4] + Av[i+4]); 00540 funVal1=temp /2; 00541 00542 temp=0; 00543 m=nn%5; 00544 if (m!=0){ 00545 for (i=0;i<m;i++) 00546 temp+=fabs(g[i]); 00547 } 00548 for (i=m;i<nn;i+=5) 00549 temp=temp + fabs(g[i]) 00550 + fabs(g[i+1]) 00551 + fabs(g[i+2]) 00552 + fabs(g[i+3]) 00553 + fabs(g[i+4]); 00554 funVal1=funVal1+ temp*lambda; 00555 00556 00557 /* 00558 we compute the solution by the second way 00559 */ 00560 00561 numS= supportSet(x, v, z, g, S, lambda, nn); 00562 00563 /* 00564 we compute the objective function of x computed in the second way 00565 */ 00566 00567 temp=0; 00568 m=n%5; 00569 if (m!=0){ 00570 for (i=0;i<m;i++) 00571 temp+=(x[i]-v[i]) * (x[i]-v[i]); 00572 } 00573 for (i=m;i<n;i+=5) 00574 temp=temp + (x[i]- v[i]) * ( x[i]- v[i]) 00575 + (x[i+1]-v[i+1]) * (x[i+1]-v[i+1]) 00576 + (x[i+2]-v[i+2]) * (x[i+2]-v[i+2]) 00577 + (x[i+3]-v[i+3]) * (x[i+3]-v[i+3]) 00578 + (x[i+4]-v[i+4]) * (x[i+4]-v[i+4]); 00579 funVal2=temp/2; 00580 00581 temp=0; 00582 m=nn%5; 00583 if (m!=0){ 00584 for (i=0;i<m;i++) 00585 temp+=fabs( x[i+1]-x[i] ); 00586 } 00587 for (i=m;i<nn;i+=5) 00588 temp=temp + fabs( x[i+1]-x[i] ) 00589 + fabs( x[i+2]-x[i+1] ) 00590 + fabs( x[i+3]-x[i+2] ) 00591 + fabs( x[i+4]-x[i+3] ) 00592 + fabs( x[i+5]-x[i+4] ); 00593 funVal2=funVal2 + lambda * temp; 00594 00595 00596 /* 00597 printf("\n funVal1=%e, funVal2=%e, diff=%e\n", funVal1, funVal2, funVal1-funVal2); 00598 */ 00599 00600 00601 00602 00603 if (funVal2 > funVal1){ /* 00604 we compute the solution by the first way 00605 */ 00606 x[0]=v[0] + z[0]; 00607 for(i=1;i<n-1;i++) 00608 x[i]= v[i] - z[i-1] + z[i]; 00609 x[n-1]=v[n-1] - z[n-2]; 00610 } 00611 else{ 00612 00613 /* 00614 the solution x is computed in the second way 00615 the gap can be further reduced 00616 (note that, there might be numerical error) 00617 */ 00618 00619 *gap=*gap - (funVal1- funVal2); 00620 if (*gap <0) 00621 *gap=0; 00622 } 00623 00624 return (numS); 00625 } 00626 00627 00628 void restartMapping(double *g, double *z, double * v, 00629 double lambda, int nn) 00630 { 00631 00632 int i, n=nn+1; 00633 double temp; 00634 int* S=(int *) malloc(sizeof(int)*nn); 00635 double *x=(double *)malloc(sizeof(double)*n); 00636 double *s=(double *)malloc(sizeof(double)*nn); 00637 double *Av=(double *)malloc(sizeof(double)*nn); 00638 //int numS=-1; 00639 00640 /* 00641 for a given input z, 00642 we compute the z0 after restarting 00643 00644 The elements in z lie in [-lambda, lambda] 00645 00646 The returned value is g 00647 */ 00648 00649 00650 for (i=0;i<nn; i++) 00651 Av[i]=v[i+1]-v[i]; 00652 00653 00654 00655 g[0]=z[0] + z[0] - z[1] - Av[0]; 00656 for (i=1;i<nn-1;i++){ 00657 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i]; 00658 } 00659 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1]; 00660 00661 00662 //numS = supportSet(x, v, z, g, S, lambda, nn); 00663 00664 00665 /*With x, we compute z via 00666 AA^T z = Av - Ax 00667 */ 00668 00669 /* 00670 compute s= Av -Ax 00671 */ 00672 00673 for (i=0;i<nn; i++) 00674 s[i]=Av[i] - x[i+1] + x[i]; 00675 00676 00677 /* 00678 Apply Rose Algorithm for solving z 00679 */ 00680 00681 Thomas(&temp, g, s, nn); 00682 00683 /* 00684 Rose(&temp, g, s, nn); 00685 */ 00686 00687 /* 00688 project g to [-lambda, lambda] 00689 */ 00690 00691 for(i=0;i<nn;i++){ 00692 if (g[i]>lambda) 00693 g[i]=lambda; 00694 else 00695 if (g[i]<-lambda) 00696 g[i]=-lambda; 00697 } 00698 00699 00700 free (S); 00701 free (x); 00702 free (s); 00703 free (Av); 00704 00705 } 00706 00707 00708 00709 /* 00710 00712 00713 Our objective is to solve the fused Lasso signal approximator (flsa) problem: 00714 00715 min_x g(x) 1/2 \|x-v\|^2 + lambda \|A x\|_1, (1) 00716 00717 Let x* be the solution (which is unique), it satisfies 00718 00719 0 in x* - v + A^T * lambda *SGN(Ax*) (2) 00720 00721 To solve x*, it suffices to find 00722 00723 y* in A^T * lambda *SGN(Ax*) (3) 00724 that satisfies 00725 00726 x* - v + y* =0 (4) 00727 which leads to 00728 x*= v - y* (5) 00729 00730 Due to the uniqueness of x*, we conclude that y* is unique. 00731 00732 As y* is a subgradient of lambda \|A x*\|_1, 00733 we name our method as Subgradient Finding Algorithm (sfa). 00734 00735 y* in (3) can be further written as 00736 00737 y*= A^T * z* (6) 00738 where 00739 00740 z* in lambda* SGN (Ax*) (7) 00741 00742 From (6), we have 00743 z* = (A A^T)^{-1} A * y* (8) 00744 00745 Therefore, from the uqniueness of y*, we conclude that z* is also unique. 00746 Next, we discuss how to solve this unique z*. 00747 00748 The problem (1) can reformulated as the following equivalent problem: 00749 00750 min_x max_z f(x, z)= 1/2 \|x-v\|^2 + <A x, z> 00751 subject to \|z\|_{infty} \leq lambda (9) 00752 00753 At the saddle point, we have 00754 00755 x = v - AT z, (10) 00756 00757 which somehow concides with (5) and (6) 00758 00759 Plugging (10) into (9), we obtain the problem 00760 00761 min_z 1/2 z^T A AT z - < z, A v>, 00762 subject to \|z\|_{infty} \leq lambda, (11) 00763 00764 In this program, we apply the Nesterov's method for solving (11). 00765 00766 00767 Duality gap: 00768 00769 At a given point z0, we compute x0= v - A^T z0. 00770 It is easy to show that 00771 min_x f(x, z0) = f(x0, z0) <= max_z f(x0, z) (12) 00772 00773 Moreover, we have 00774 max_z f(x0, z) - min_x f(x, z0) 00775 <= lambda * \|A x0\|_1 - < z0, Av - A A^T z0> (13) 00776 00777 It is also to get that 00778 00779 f(x0, z0) <= f(x*, z*) <= max_z f(x0, z) (14) 00780 00781 g(x*)=f(x*, z*) (15) 00782 00783 g(x0)=max_z f(x0, z) (17) 00784 00785 Therefore, we have 00786 00787 g(x0)-g(x*) <= lambda * \|A x0\|_1 - < z0, Av - A A^T z0> (18) 00788 00789 00790 We have applied a restarting technique, which is quite involved; and thus, we do not explain here. 00791 00793 */ 00794 00795 00796 /* 00798 00799 For sfa, the stepsize of the Nesterov's method is fixed to 1/4, so that no line search is needed. 00800 00801 00802 00803 Explanation of the parameters: 00804 00805 Output parameters 00806 x: the solution to the primal problem 00807 gap: the duality gap (pointer) 00808 00809 Input parameters 00810 z: the solution to the dual problem (before calling this function, z contains a starting point) 00811 !!!!we assume that the starting point has been successfully initialized in z !!!! 00812 z0: a variable used for multiple purposes: 00813 1) the previous solution z0 00814 2) the difference between z and z0, i.e., z0=z- z0 00815 00816 lambda: the regularization parameter (and the radius of the infity ball, see (11)). 00817 nn: the length of z, z0, Av, g, and s 00818 maxStep: the maximal number of iterations 00819 00820 v: the point to be projected (not changed after the program) 00821 Av: A*v (not changed after the program) 00822 00823 s: the search point (used for multiple purposes) 00824 g: the gradient at g (and it is also used for multiple purposes) 00825 00826 tol: the tolerance of the gap 00827 tau: the duality gap or the restarting technique is done every tau steps 00828 flag: if flag=1, we apply the resart technique 00829 flag=2, just run the SFA algorithm, terminate it when the absolution change is less than tol 00830 flag=3, just run the SFA algorithm, terminate it when the duality gap is less than tol 00831 flag=4, just run the SFA algorithm, terminate it when the relative duality gap is less than tol 00832 00833 00834 We would like to emphasis that the following assumptions 00835 have been checked in the functions that call this function: 00836 1) 0< lambda < z_max 00837 2) nn >=2 00838 3) z has been initialized with a starting point 00839 4) z0 has been initialized with all zeros 00840 00841 The termination condition is checked every tau iterations. 00842 00843 For the duality gap, please refer to (12-18) 00844 */ 00845 00846 int sfa(double *x, double *gap, int * activeS, 00847 double *z, double *z0, double * v, double * Av, 00848 double lambda, int nn, int maxStep, 00849 double *s, double *g, 00850 double tol, int tau, int flag){ 00851 00852 int i, iterStep, m, tFlag=0, n=nn+1; 00853 double alphap=0, alpha=1, beta=0, temp; 00854 int* S=(int *) malloc(sizeof(int)*nn); 00855 double gapp=-1, gappp=-1; /*gapp denotes the previous gap*/ 00856 int numS=-1, numSp=-2, numSpp=-3;; 00857 /* 00858 numS denotes the number of elements in the Support Set S 00859 numSp denotes the number of elements in the previous Support Set S 00860 */ 00861 00862 *gap=-1; /*initial a value -1*/ 00863 00864 /* 00865 The main algorithm by Nesterov's method 00866 00867 B is an nn x nn tridiagonal matrix. 00868 00869 The nn eigenvalues of B are 2- 2 cos (i * PI/ n), i=1, 2, ..., nn 00870 */ 00871 00872 for (iterStep=1; iterStep<=maxStep; iterStep++){ 00873 00874 00875 /*------------- Step 1 ---------------------*/ 00876 00877 beta=(alphap -1 ) / alpha; 00878 /* 00879 compute search point 00880 00881 s= z + beta * z0 00882 00883 We follow the style of CLAPACK 00884 */ 00885 m=nn % 5; 00886 if (m!=0){ 00887 for (i=0;i<m; i++) 00888 s[i]=z[i]+ beta* z0[i]; 00889 } 00890 for (i=m;i<nn;i+=5){ 00891 s[i] =z[i] + beta* z0[i]; 00892 s[i+1] =z[i+1] + beta* z0[i+1]; 00893 s[i+2] =z[i+2] + beta* z0[i+2]; 00894 s[i+3] =z[i+3] + beta* z0[i+3]; 00895 s[i+4] =z[i+4] + beta* z0[i+4]; 00896 } 00897 00898 /* 00899 s and g are of size nn x 1 00900 00901 compute the gradient at s 00902 00903 g= B * s - Av, 00904 00905 where B is an nn x nn tridiagonal matrix. and is defined as 00906 00907 B= [ 2 -1 0 0; 00908 -1 2 -1 0; 00909 0 -1 2 -1; 00910 0 0 -1 2] 00911 00912 We assume n>=3, which leads to nn>=2 00913 */ 00914 g[0]=s[0] + s[0] - s[1] - Av[0]; 00915 for (i=1;i<nn-1;i++){ 00916 g[i]= - s[i-1] + s[i] + s[i] - s[i+1] - Av[i]; 00917 } 00918 g[nn-1]= -s[nn-2] + s[nn-1] + s[nn-1] - Av[nn-1]; 00919 00920 00921 /* 00922 z0 stores the previous -z 00923 */ 00924 m=nn%7; 00925 if (m!=0){ 00926 for (i=0;i<m;i++) 00927 z0[i]=-z[i]; 00928 } 00929 for (i=m; i <nn; i+=7){ 00930 z0[i] = - z[i]; 00931 z0[i+1] = - z[i+1]; 00932 z0[i+2] = - z[i+2]; 00933 z0[i+3] = - z[i+3]; 00934 z0[i+4] = - z[i+4]; 00935 z0[i+5] = - z[i+5]; 00936 z0[i+6] = - z[i+6]; 00937 } 00938 00939 00940 /* 00941 do a gradient step based on s to get z 00942 */ 00943 m=nn%5; 00944 if (m!=0){ 00945 for(i=0;i<m; i++) 00946 z[i]=s[i] - g[i]/4; 00947 } 00948 for (i=m;i<nn; i+=5){ 00949 z[i] = s[i] - g[i] /4; 00950 z[i+1] = s[i+1] - g[i+1]/4; 00951 z[i+2] = s[i+2] - g[i+2]/4; 00952 z[i+3] = s[i+3] - g[i+3]/4; 00953 z[i+4] = s[i+4] - g[i+4]/4; 00954 } 00955 00956 /* 00957 project z onto the L_{infty} ball with radius lambda 00958 00959 z is the new approximate solution 00960 */ 00961 for (i=0;i<nn; i++){ 00962 if (z[i]>lambda) 00963 z[i]=lambda; 00964 else 00965 if (z[i]<-lambda) 00966 z[i]=-lambda; 00967 } 00968 00969 /* 00970 compute the difference between the new solution 00971 and the previous solution (stored in z0=-z_p) 00972 00973 the difference is written to z0 00974 */ 00975 00976 m=nn%5; 00977 if (m!=0){ 00978 for (i=0;i<m;i++) 00979 z0[i]+=z[i]; 00980 } 00981 for(i=m;i<nn; i+=5){ 00982 z0[i] +=z[i]; 00983 z0[i+1]+=z[i+1]; 00984 z0[i+2]+=z[i+2]; 00985 z0[i+3]+=z[i+3]; 00986 z0[i+4]+=z[i+4]; 00987 } 00988 00989 00990 alphap=alpha; 00991 alpha=(1+sqrt(4*alpha*alpha+1))/2; 00992 00993 /* 00994 check the termination condition 00995 */ 00996 if (iterStep%tau==0){ 00997 00998 00999 /* 01000 The variables g and s can be modified 01001 01002 The variables x, z0 and z can be revised for case 0, but not for the rest 01003 */ 01004 switch (flag){ 01005 case 1: 01006 01007 /* 01008 01009 terminate the program once the "duality gap" is smaller than tol 01010 01011 compute the duality gap: 01012 01013 x= v - A^T z 01014 Ax = Av - A A^T z = -g, 01015 where 01016 g = A A^T z - A v 01017 01018 01019 the duality gap= lambda * \|Ax\|-1 - <z, Ax> 01020 = lambda * \|g\|_1 + <z, g> 01021 01022 In fact, gap=0 indicates that, 01023 if g_i >0, then z_i=-lambda 01024 if g_i <0, then z_i=lambda 01025 */ 01026 01027 gappp=gapp; 01028 gapp=*gap; /*record the previous gap*/ 01029 numSpp=numSp; 01030 numSp=numS; /*record the previous numS*/ 01031 01032 dualityGap(gap, z, g, s, Av, lambda, nn); 01033 /*g is computed as the gradient of z in this function*/ 01034 01035 01036 /* 01037 printf("\n Iteration: %d, gap=%e, numS=%d", iterStep, *gap, numS); 01038 */ 01039 01040 /* 01041 If *gap <=tol, we terminate the iteration 01042 Otherwise, we restart the algorithm 01043 */ 01044 01045 if (*gap <=tol){ 01046 tFlag=1; 01047 break; 01048 01049 } /* end of *gap <=tol */ 01050 else{ 01051 01052 /* we apply the restarting technique*/ 01053 01054 /* 01055 we compute the solution by the second way 01056 */ 01057 numS = supportSet(x, v, z, g, S, lambda, nn); 01058 /*g, the gradient of z should be computed before calling this function*/ 01059 01060 /*With x, we compute z via 01061 AA^T z = Av - Ax 01062 */ 01063 01064 /* 01065 printf("\n iterStep=%d, numS=%d, gap=%e",iterStep, numS, *gap); 01066 */ 01067 01068 01069 m=1; 01070 if (nn > 1000000) 01071 m=10; 01072 else 01073 if (nn > 100000) 01074 m=5; 01075 01076 if ( abs(numS-numSp) < m){ 01077 01078 numS=generateSolution(x, z, gap, v, Av, 01079 g, s, S, lambda, nn); 01080 /*g, the gradient of z should be computed before calling this function*/ 01081 01082 01083 if (*gap <tol){ 01084 tFlag=2; /*tFlag =2 shows that the result is already optimal 01085 There is no need to call generateSolution for recomputing the best solution 01086 */ 01087 break; 01088 } 01089 01090 if ( (*gap ==gappp) && (numS==numSpp) ){ 01091 01092 tFlag=2; 01093 break; 01094 01095 } 01096 01097 /*we terminate the program is *gap <1 01098 numS==numSP 01099 and gapp==*gap 01100 */ 01101 } 01102 01103 /* 01104 compute s= Av -Ax 01105 */ 01106 for (i=0;i<nn; i++) 01107 s[i]=Av[i] - x[i+1] + x[i]; 01108 01109 /* 01110 apply Rose Algorithm for solving z 01111 */ 01112 01113 Thomas(&temp, z, s, nn); 01114 01115 /* 01116 Rose(&temp, z, s, nn); 01117 */ 01118 01119 /* 01120 printf("\n Iteration: %d, %e", iterStep, temp); 01121 */ 01122 01123 /* 01124 project z to [-lambda2, lambda2] 01125 */ 01126 for(i=0;i<nn;i++){ 01127 if (z[i]>lambda) 01128 z[i]=lambda; 01129 else 01130 if (z[i]<-lambda) 01131 z[i]=-lambda; 01132 } 01133 01134 01135 01136 m=nn%7; 01137 if (m!=0){ 01138 for (i=0;i<m;i++) 01139 z0[i]=0; 01140 } 01141 for (i=m; i<nn; i+=7){ 01142 z0[i] = z0[i+1] 01143 = z0[i+2] 01144 = z0[i+3] 01145 = z0[i+4] 01146 = z0[i+5] 01147 = z0[i+6] 01148 =0; 01149 } 01150 01151 01152 alphap=0; alpha=1; 01153 01154 /* 01155 we restart the algorithm 01156 */ 01157 01158 } 01159 01160 break; /*break case 1*/ 01161 01162 case 2: 01163 01164 /* 01165 The program is terminated either the summation of the absolution change (denoted by z0) 01166 of z (from the previous zp) is less than tol * nn, 01167 or the maximal number of iteration (maxStep) is achieved 01168 Note: tol indeed measures the averaged per element change. 01169 */ 01170 temp=0; 01171 m=nn%5; 01172 if (m!=0){ 01173 for(i=0;i<m;i++) 01174 temp+=fabs(z0[i]); 01175 } 01176 for(i=m;i<nn;i+=5) 01177 temp=temp + fabs(z0[i]) 01178 + fabs(z0[i+1]) 01179 + fabs(z0[i+2]) 01180 + fabs(z0[i+3]) 01181 + fabs(z0[i+4]); 01182 *gap=temp / nn; 01183 01184 if (*gap <=tol){ 01185 01186 tFlag=1; 01187 } 01188 01189 break; 01190 01191 case 3: 01192 01193 /* 01194 01195 terminate the program once the "duality gap" is smaller than tol 01196 01197 compute the duality gap: 01198 01199 x= v - A^T z 01200 Ax = Av - A A^T z = -g, 01201 where 01202 g = A A^T z - A v 01203 01204 01205 the duality gap= lambda * \|Ax\|-1 - <z, Ax> 01206 = lambda * \|g\|_1 + <z, g> 01207 01208 In fact, gap=0 indicates that, 01209 if g_i >0, then z_i=-lambda 01210 if g_i <0, then z_i=lambda 01211 */ 01212 01213 01214 g[0]=z[0] + z[0] - z[1] - Av[0]; 01215 for (i=1;i<nn-1;i++){ 01216 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i]; 01217 } 01218 01219 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1]; 01220 01221 for (i=0;i<nn;i++) 01222 if (g[i]>0) 01223 s[i]=lambda + z[i]; 01224 else 01225 s[i]=-lambda + z[i]; 01226 01227 temp=0; 01228 m=nn%5; 01229 if (m!=0){ 01230 for(i=0;i<m;i++) 01231 temp+=s[i]*g[i]; 01232 } 01233 for(i=m;i<nn;i+=5) 01234 temp=temp + s[i] *g[i] 01235 + s[i+1]*g[i+1] 01236 + s[i+2]*g[i+2] 01237 + s[i+3]*g[i+3] 01238 + s[i+4]*g[i+4]; 01239 *gap=temp; 01240 01241 /* 01242 printf("\n %e", *gap); 01243 */ 01244 01245 01246 if (*gap <=tol) 01247 tFlag=1; 01248 01249 break; 01250 01251 case 4: 01252 01253 /* 01254 01255 terminate the program once the "relative duality gap" is smaller than tol 01256 01257 01258 compute the duality gap: 01259 01260 x= v - A^T z 01261 Ax = Av - A A^T z = -g, 01262 where 01263 g = A A^T z - A v 01264 01265 01266 the duality gap= lambda * \|Ax\|-1 - <z, Ax> 01267 = lambda * \|g\|_1 + <z, g> 01268 01269 In fact, gap=0 indicates that, 01270 if g_i >0, then z_i=-lambda 01271 if g_i <0, then z_i=lambda 01272 01273 01274 Here, the "relative duality gap" is defined as: 01275 duality gap / - 1/2 \|A^T z\|^2 + < z, Av> 01276 01277 We efficiently compute - 1/2 \|A^T z\|^2 + < z, Av> using the following relationship 01278 01279 - 1/2 \|A^T z\|^2 + < z, Av> 01280 = -1/2 <z, A A^T z - Av -Av> 01281 = -1/2 <z, g - Av> 01282 */ 01283 01284 01285 g[0]=z[0] + z[0] - z[1] - Av[0]; 01286 for (i=1;i<nn-1;i++){ 01287 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i]; 01288 } 01289 01290 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1]; 01291 01292 for (i=0;i<nn;i++) 01293 if (g[i]>0) 01294 s[i]=lambda + z[i]; 01295 else 01296 s[i]=-lambda + z[i]; 01297 01298 temp=0; 01299 m=nn%5; 01300 if (m!=0){ 01301 for(i=0;i<m;i++) 01302 temp+=s[i]*g[i]; 01303 } 01304 for(i=m;i<nn;i+=5) 01305 temp=temp + s[i] *g[i] 01306 + s[i+1]*g[i+1] 01307 + s[i+2]*g[i+2] 01308 + s[i+3]*g[i+3] 01309 + s[i+4]*g[i+4]; 01310 *gap=temp; 01311 /* 01312 Now, *gap contains the duality gap 01313 Next, we compute 01314 - 1/2 \|A^T z\|^2 + < z, Av> 01315 =-1/2 <z, g - Av> 01316 */ 01317 01318 temp=0; 01319 m=nn%5; 01320 if (m!=0){ 01321 for(i=0;i<m;i++) 01322 temp+=z[i] * (g[i] - Av[i]); 01323 } 01324 for(i=m;i<nn;i+=5) 01325 temp=temp + z[i] * (g[i] - Av[i]) 01326 + z[i+1]* (g[i+1]- Av[i+1]) 01327 + z[i+2]* (g[i+2]- Av[i+2]) 01328 + z[i+3]* (g[i+3]- Av[i+3]) 01329 + z[i+4]* (g[i+4]- Av[i+4]); 01330 temp=fabs(temp) /2; 01331 01332 if (temp <1) 01333 temp=1; 01334 01335 *gap/=temp; 01336 /* 01337 *gap now contains the relative gap 01338 */ 01339 01340 01341 if (*gap <=tol){ 01342 tFlag=1; 01343 } 01344 01345 break; 01346 01347 default: 01348 01349 /* 01350 The program is terminated either the summation of the absolution change (denoted by z0) 01351 of z (from the previous zp) is less than tol * nn, 01352 or the maximal number of iteration (maxStep) is achieved 01353 Note: tol indeed measures the averaged per element change. 01354 */ 01355 temp=0; 01356 m=nn%5; 01357 if (m!=0){ 01358 for(i=0;i<m;i++) 01359 temp+=fabs(z0[i]); 01360 } 01361 for(i=m;i<nn;i+=5) 01362 temp=temp + fabs(z0[i]) 01363 + fabs(z0[i+1]) 01364 + fabs(z0[i+2]) 01365 + fabs(z0[i+3]) 01366 + fabs(z0[i+4]); 01367 *gap=temp / nn; 01368 01369 if (*gap <=tol){ 01370 01371 tFlag=1; 01372 } 01373 01374 break; 01375 01376 }/*end of switch*/ 01377 01378 if (tFlag) 01379 break; 01380 01381 }/* end of the if for checking the termination condition */ 01382 01383 /*-------------- Step 3 --------------------*/ 01384 01385 } 01386 01387 /* 01388 for the other cases, except flag=1, compute the solution x according the first way (the primal-dual way) 01389 */ 01390 01391 if ( (flag !=1) || (tFlag==0) ){ 01392 x[0]=v[0] + z[0]; 01393 for(i=1;i<n-1;i++) 01394 x[i]= v[i] - z[i-1] + z[i]; 01395 x[n-1]=v[n-1] - z[n-2]; 01396 } 01397 01398 if ( (flag==1) && (tFlag==1)){ 01399 01400 /* 01401 We assume that n>=3, and thus nn>=2 01402 01403 We have two ways for recovering x. 01404 The first way is x = v - A^T z 01405 The second way is x =omega(z) 01406 */ 01407 01408 /* 01409 We first compute the objective function value of the first choice in terms f(x), see our paper 01410 */ 01411 01412 /* 01413 for numerical reason, we do a gradient descent step 01414 */ 01415 01416 /* 01417 --------------------------------------------------- 01418 A gradient step begins 01419 */ 01420 g[0]=z[0] + z[0] - z[1] - Av[0]; 01421 for (i=1;i<nn-1;i++){ 01422 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i]; 01423 } 01424 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1]; 01425 01426 01427 /* 01428 do a gradient step based on z to get the new z 01429 */ 01430 m=nn%5; 01431 if (m!=0){ 01432 for(i=0;i<m; i++) 01433 z[i]=z[i] - g[i]/4; 01434 } 01435 for (i=m;i<nn; i+=5){ 01436 z[i] = z[i] - g[i] /4; 01437 z[i+1] = z[i+1] - g[i+1]/4; 01438 z[i+2] = z[i+2] - g[i+2]/4; 01439 z[i+3] = z[i+3] - g[i+3]/4; 01440 z[i+4] = z[i+4] - g[i+4]/4; 01441 } 01442 01443 /* 01444 project z onto the L_{infty} ball with radius lambda 01445 01446 z is the new approximate solution 01447 */ 01448 for (i=0;i<nn; i++){ 01449 if (z[i]>lambda) 01450 z[i]=lambda; 01451 else 01452 if (z[i]<-lambda) 01453 z[i]=-lambda; 01454 } 01455 01456 /* 01457 --------------------------------------------------- 01458 A gradient descent step ends 01459 */ 01460 01461 /*compute the gradient at z*/ 01462 01463 g[0]=z[0] + z[0] - z[1] - Av[0]; 01464 for (i=1;i<nn-1;i++){ 01465 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i]; 01466 } 01467 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1]; 01468 01469 01470 numS=generateSolution(x, z, gap, v, Av, 01471 g, s, S, lambda, nn); 01472 /*g, the gradient of z should be computed before calling this function*/ 01473 01474 } 01475 01476 free (S); 01477 /* 01478 free the variables S 01479 */ 01480 01481 *activeS=numS; 01482 return (iterStep); 01483 01484 } 01485 01486 01487 /* 01488 01489 Refer to sfa for the defintions of the variables 01490 01491 In this file, we restart the program every step, and neglect the gradient step. 01492 01493 It seems that, this program does not converge. 01494 01495 This function shows that the gradient step is necessary. 01496 */ 01497 01498 int sfa_special(double *x, double *gap, int * activeS, 01499 double *z, double * v, double * Av, 01500 double lambda, int nn, int maxStep, 01501 double *s, double *g, 01502 double tol, int tau){ 01503 01504 int i, iterStep; 01505 //int tFlag=0; 01506 //int n=nn+1; 01507 double temp; 01508 int* S=(int *) malloc(sizeof(int)*nn); 01509 double gapp=-1; /*gapp denotes the previous gap*/ 01510 int numS=-1, numSp=-1; 01511 /* 01512 numS denotes the number of elements in the Support Set S 01513 numSp denotes the number of elements in the previous Support Set S 01514 */ 01515 01516 *gap=-1; /*initialize *gap a value*/ 01517 01518 for (iterStep=1; iterStep<=maxStep; iterStep++){ 01519 01520 01521 g[0]=z[0] + z[0] - z[1] - Av[0]; 01522 for (i=1;i<nn-1;i++){ 01523 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i]; 01524 } 01525 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1]; 01526 01527 numSp=numS; /*record the previous numS*/ 01528 numS = supportSet(x, v, z, g, S, lambda, nn); 01529 01530 01531 /*With x, we compute z via 01532 AA^T z = Av - Ax 01533 */ 01534 01535 /* 01536 compute s= Av -Ax 01537 */ 01538 01539 for (i=0;i<nn; i++) 01540 s[i]=Av[i] - x[i+1] + x[i]; 01541 01542 01543 /* 01544 Apply Rose Algorithm for solving z 01545 */ 01546 01547 Thomas(&temp, z, s, nn); 01548 01549 /* 01550 Rose(&temp, z, s, nn); 01551 */ 01552 01553 /* 01554 project z to [-lambda, lambda] 01555 */ 01556 01557 for(i=0;i<nn;i++){ 01558 if (z[i]>lambda) 01559 z[i]=lambda; 01560 else 01561 if (z[i]<-lambda) 01562 z[i]=-lambda; 01563 } 01564 01565 01566 if (iterStep%tau==0){ 01567 gapp=*gap; /*record the previous gap*/ 01568 01569 dualityGap(gap, z, g, s, Av, lambda, nn); 01570 01571 /* 01572 printf("\n iterStep=%d, numS=%d, gap=%e, diff=%e",iterStep, numS, *gap, *gap -gapp); 01573 01574 */ 01575 01576 if (*gap <=tol){ 01577 //tFlag=1; 01578 break; 01579 } 01580 01581 if ( (*gap <1) && (numS==numSp) && fabs(gapp == *gap) ){ 01582 //tFlag=1; 01583 break; 01584 /*we terminate the program is *gap <1 01585 numS==numSP 01586 and gapp==*gap 01587 */ 01588 } 01589 01590 }/*end of if tau*/ 01591 01592 }/*end for */ 01593 01594 free (S); 01595 01596 * activeS=numS; 01597 return(iterStep); 01598 01599 } 01600 01601 01602 /* 01603 01604 We do one gradient descent, and then restart the program 01605 */ 01606 01607 01608 int sfa_one(double *x, double *gap, int * activeS, 01609 double *z, double * v, double * Av, 01610 double lambda, int nn, int maxStep, 01611 double *s, double *g, 01612 double tol, int tau){ 01613 01614 int i, iterStep, m; 01615 int tFlag=0; 01616 //int n=nn+1; 01617 double temp; 01618 int* S=(int *) malloc(sizeof(int)*nn); 01619 double gapp=-1, gappp=-2; /*gapp denotes the previous gap*/ 01620 int numS=-100, numSp=-200, numSpp=-300; 01621 /* 01622 numS denotes the number of elements in the Support Set S 01623 numSp denotes the number of elements in the previous Support Set S 01624 */ 01625 01626 *gap=-1; /*initialize *gap a value*/ 01627 01628 /* 01629 The main algorithm by Nesterov's method 01630 01631 B is an nn x nn tridiagonal matrix. 01632 01633 The nn eigenvalues of B are 2- 2 cos (i * PI/ n), i=1, 2, ..., nn 01634 */ 01635 01636 01637 /* 01638 we first do a gradient step based on z 01639 */ 01640 01641 01642 /* 01643 --------------------------------------------------- 01644 A gradient step begins 01645 */ 01646 g[0]=z[0] + z[0] - z[1] - Av[0]; 01647 for (i=1;i<nn-1;i++){ 01648 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i]; 01649 } 01650 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1]; 01651 01652 01653 /* 01654 do a gradient step based on z to get the new z 01655 */ 01656 m=nn%5; 01657 if (m!=0){ 01658 for(i=0;i<m; i++) 01659 z[i]=z[i] - g[i]/4; 01660 } 01661 for (i=m;i<nn; i+=5){ 01662 z[i] = z[i] - g[i] /4; 01663 z[i+1] = z[i+1] - g[i+1]/4; 01664 z[i+2] = z[i+2] - g[i+2]/4; 01665 z[i+3] = z[i+3] - g[i+3]/4; 01666 z[i+4] = z[i+4] - g[i+4]/4; 01667 } 01668 01669 /* 01670 project z onto the L_{infty} ball with radius lambda 01671 01672 z is the new approximate solution 01673 */ 01674 for (i=0;i<nn; i++){ 01675 if (z[i]>lambda) 01676 z[i]=lambda; 01677 else 01678 if (z[i]<-lambda) 01679 z[i]=-lambda; 01680 } 01681 01682 /* 01683 --------------------------------------------------- 01684 A gradient descent step ends 01685 */ 01686 01687 01688 /*compute the gradient at z*/ 01689 01690 g[0]=z[0] + z[0] - z[1] - Av[0]; 01691 for (i=1;i<nn-1;i++){ 01692 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i]; 01693 } 01694 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1]; 01695 01696 for (iterStep=1; iterStep<=maxStep; iterStep++){ 01697 01698 01699 /* 01700 --------------------------------------------------- 01701 restart the algorithm with x=omega(z) 01702 */ 01703 01704 numSpp=numSp; 01705 numSp=numS; /*record the previous numS*/ 01706 numS = supportSet(x, v, z, g, S, lambda, nn); 01707 01708 01709 /*With x, we compute z via 01710 AA^T z = Av - Ax 01711 */ 01712 01713 /* 01714 compute s= Av -Ax 01715 */ 01716 01717 for (i=0;i<nn; i++) 01718 s[i]=Av[i] - x[i+1] + x[i]; 01719 01720 01721 /* 01722 Apply Rose Algorithm for solving z 01723 */ 01724 01725 Thomas(&temp, z, s, nn); 01726 01727 /* 01728 Rose(&temp, z, s, nn); 01729 */ 01730 01731 /* 01732 project z to [-lambda, lambda] 01733 */ 01734 01735 for(i=0;i<nn;i++){ 01736 if (z[i]>lambda) 01737 z[i]=lambda; 01738 else 01739 if (z[i]<-lambda) 01740 z[i]=-lambda; 01741 } 01742 01743 /* 01744 --------------------------------------------------- 01745 restart the algorithm with x=omega(z) 01746 01747 we have computed a new z, based on the above relationship 01748 */ 01749 01750 01751 /* 01752 --------------------------------------------------- 01753 A gradient step begins 01754 */ 01755 g[0]=z[0] + z[0] - z[1] - Av[0]; 01756 for (i=1;i<nn-1;i++){ 01757 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i]; 01758 } 01759 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1]; 01760 01761 01762 /* 01763 do a gradient step based on z to get the new z 01764 */ 01765 m=nn%5; 01766 if (m!=0){ 01767 for(i=0;i<m; i++) 01768 z[i]=z[i] - g[i]/4; 01769 } 01770 for (i=m;i<nn; i+=5){ 01771 z[i] = z[i] - g[i] /4; 01772 z[i+1] = z[i+1] - g[i+1]/4; 01773 z[i+2] = z[i+2] - g[i+2]/4; 01774 z[i+3] = z[i+3] - g[i+3]/4; 01775 z[i+4] = z[i+4] - g[i+4]/4; 01776 } 01777 01778 /* 01779 project z onto the L_{infty} ball with radius lambda 01780 01781 z is the new approximate solution 01782 */ 01783 for (i=0;i<nn; i++){ 01784 if (z[i]>lambda) 01785 z[i]=lambda; 01786 else 01787 if (z[i]<-lambda) 01788 z[i]=-lambda; 01789 } 01790 01791 /* 01792 --------------------------------------------------- 01793 A gradient descent step ends 01794 */ 01795 01796 /*compute the gradient at z*/ 01797 01798 g[0]=z[0] + z[0] - z[1] - Av[0]; 01799 for (i=1;i<nn-1;i++){ 01800 g[i]= - z[i-1] + z[i] + z[i] - z[i+1] - Av[i]; 01801 } 01802 g[nn-1]= -z[nn-2] + z[nn-1] + z[nn-1] - Av[nn-1]; 01803 01804 01805 if (iterStep % tau==0){ 01806 gappp=gapp; 01807 gapp=*gap; /*record the previous gap*/ 01808 01809 dualityGap2(gap, z, g, s, Av, lambda, nn); 01810 /*g, the gradient of z should be computed before calling this function*/ 01811 01812 01813 /* 01814 printf("\n iterStep=%d, numS=%d, gap=%e",iterStep, numS, *gap); 01815 */ 01816 01817 01818 /* 01819 printf("\n %d & %d & %2.0e \\\\ \n \\hline ",iterStep, numS, *gap); 01820 */ 01821 01822 01823 /* 01824 printf("\n %e",*gap); 01825 */ 01826 01827 /* 01828 01829 printf("\n %d",numS); 01830 01831 */ 01832 01833 if (*gap <=tol){ 01834 //tFlag=1; 01835 break; 01836 } 01837 01838 m=1; 01839 if (nn > 1000000) 01840 m=5; 01841 else 01842 if (nn > 100000) 01843 m=3; 01844 01845 if ( abs( numS-numSp) <m ){ 01846 01847 /* 01848 printf("\n numS=%d, numSp=%d",numS,numSp); 01849 */ 01850 01851 m=generateSolution(x, z, gap, v, Av, 01852 g, s, S, lambda, nn); 01853 /*g, the gradient of z should be computed before calling this function*/ 01854 01855 if (*gap < tol){ 01856 01857 numS=m; 01858 tFlag=2; 01859 break; 01860 } 01861 01862 01863 if ( (*gap ==gappp) && (numS==numSpp) ){ 01864 01865 tFlag=2; 01866 break; 01867 01868 } 01869 01870 } /*end of if*/ 01871 01872 }/*end of if tau*/ 01873 01874 01875 } /*end of for*/ 01876 01877 01878 01879 if (tFlag!=2){ 01880 numS=generateSolution(x, z, gap, v, Av, g, s, S, lambda, nn); 01881 /*g, the gradient of z should be computed before calling this function*/ 01882 } 01883 01884 free(S); 01885 01886 *activeS=numS; 01887 return(iterStep); 01888 }