annotate DL/RLS-DLA/SolveFISTA.m @ 40:6416fc12f2b8

(none)
author idamnjanovic
date Mon, 14 Mar 2011 15:35:24 +0000
parents
children
rev   line source
idamnjanovic@40 1 % Copyright ©2010. The Regents of the University of California (Regents).
idamnjanovic@40 2 % All Rights Reserved. Contact The Office of Technology Licensing,
idamnjanovic@40 3 % UC Berkeley, 2150 Shattuck Avenue, Suite 510, Berkeley, CA 94720-1620,
idamnjanovic@40 4 % (510) 643-7201, for commercial licensing opportunities.
idamnjanovic@40 5
idamnjanovic@40 6 % Authors: Arvind Ganesh, Allen Y. Yang, Zihan Zhou.
idamnjanovic@40 7 % Contact: Allen Y. Yang, Department of EECS, University of California,
idamnjanovic@40 8 % Berkeley. <yang@eecs.berkeley.edu>
idamnjanovic@40 9
idamnjanovic@40 10 % IN NO EVENT SHALL REGENTS BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT,
idamnjanovic@40 11 % SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS,
idamnjanovic@40 12 % ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
idamnjanovic@40 13 % REGENTS HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
idamnjanovic@40 14
idamnjanovic@40 15 % REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED
idamnjanovic@40 16 % TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
idamnjanovic@40 17 % PARTICULAR PURPOSE. THE SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF ANY,
idamnjanovic@40 18 % PROVIDED HEREUNDER IS PROVIDED "AS IS". REGENTS HAS NO OBLIGATION TO
idamnjanovic@40 19 % PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
idamnjanovic@40 20
idamnjanovic@40 21 %% This function is modified from Matlab code proximal_gradient_bp
idamnjanovic@40 22
idamnjanovic@40 23 function [x_hat,nIter] = SolveFISTA(A,b, varargin)
idamnjanovic@40 24
idamnjanovic@40 25 % b - m x 1 vector of observations/data (required input)
idamnjanovic@40 26 % A - m x n measurement matrix (required input)
idamnjanovic@40 27 %
idamnjanovic@40 28 % tol - tolerance for stopping criterion.
idamnjanovic@40 29 % - DEFAULT 1e-7 if omitted or -1.
idamnjanovic@40 30 % maxIter - maxilambdam number of iterations
idamnjanovic@40 31 % - DEFAULT 10000, if omitted or -1.
idamnjanovic@40 32 % lineSearchFlag - 1 if line search is to be done every iteration
idamnjanovic@40 33 % - DEFAULT 0, if omitted or -1.
idamnjanovic@40 34 % continuationFlag - 1 if a continuation is to be done on the parameter lambda
idamnjanovic@40 35 % - DEFAULT 1, if omitted or -1.
idamnjanovic@40 36 % eta - line search parameter, should be in (0,1)
idamnjanovic@40 37 % - ignored if lineSearchFlag is 0.
idamnjanovic@40 38 % - DEFAULT 0.9, if omitted or -1.
idamnjanovic@40 39 % lambda - relaxation parameter
idamnjanovic@40 40 % - ignored if continuationFlag is 1.
idamnjanovic@40 41 % - DEFAULT 1e-3, if omitted or -1.
idamnjanovic@40 42 % outputFileName - Details of each iteration are dumped here, if provided.
idamnjanovic@40 43 %
idamnjanovic@40 44 % x_hat - estimate of coeeficient vector
idamnjanovic@40 45 % numIter - number of iterations until convergence
idamnjanovic@40 46 %
idamnjanovic@40 47 %
idamnjanovic@40 48 % References
idamnjanovic@40 49 % "Robust PCA: Exact Recovery of Corrupted Low-Rank Matrices via Convex Optimization", J. Wright et al., preprint 2009.
idamnjanovic@40 50 % "An Accelerated Proximal Gradient Algorithm for Nuclear Norm Regularized Least Squares problems", K.-C. Toh and S. Yun, preprint 2009.
idamnjanovic@40 51 %
idamnjanovic@40 52 % Arvind Ganesh, Summer 2009. Questions? abalasu2@illinois.edu
idamnjanovic@40 53
idamnjanovic@40 54 DEBUG = 0 ;
idamnjanovic@40 55
idamnjanovic@40 56 STOPPING_GROUND_TRUTH = -1;
idamnjanovic@40 57 STOPPING_DUALITY_GAP = 1;
idamnjanovic@40 58 STOPPING_SPARSE_SUPPORT = 2;
idamnjanovic@40 59 STOPPING_OBJECTIVE_VALUE = 3;
idamnjanovic@40 60 STOPPING_SUBGRADIENT = 4;
idamnjanovic@40 61 STOPPING_DEFAULT = STOPPING_SUBGRADIENT;
idamnjanovic@40 62
idamnjanovic@40 63 stoppingCriterion = STOPPING_DEFAULT;
idamnjanovic@40 64 maxIter = 1000 ;
idamnjanovic@40 65 tolerance = 1e-3;
idamnjanovic@40 66 [m,n] = size(A) ;
idamnjanovic@40 67 x0 = zeros(n,1) ;
idamnjanovic@40 68 xG = [];
idamnjanovic@40 69
idamnjanovic@40 70 %% Initializing optimization variables
idamnjanovic@40 71 t_k = 1 ;
idamnjanovic@40 72 t_km1 = 1 ;
idamnjanovic@40 73 L0 = 1 ;
idamnjanovic@40 74 G = A'*A ;
idamnjanovic@40 75 nIter = 0 ;
idamnjanovic@40 76 c = A'*b ;
idamnjanovic@40 77 lambda0 = 0.99*L0*norm(c,inf) ;
idamnjanovic@40 78 eta = 0.6 ;
idamnjanovic@40 79 lambda_bar = 1e-4*lambda0 ;
idamnjanovic@40 80 xk = zeros(n,1) ;
idamnjanovic@40 81 lambda = lambda0 ;
idamnjanovic@40 82 L = L0 ;
idamnjanovic@40 83 beta = 1.5 ;
idamnjanovic@40 84
idamnjanovic@40 85 % Parse the optional inputs.
idamnjanovic@40 86 if (mod(length(varargin), 2) ~= 0 ),
idamnjanovic@40 87 error(['Extra Parameters passed to the function ''' mfilename ''' lambdast be passed in pairs.']);
idamnjanovic@40 88 end
idamnjanovic@40 89 parameterCount = length(varargin)/2;
idamnjanovic@40 90
idamnjanovic@40 91 for parameterIndex = 1:parameterCount,
idamnjanovic@40 92 parameterName = varargin{parameterIndex*2 - 1};
idamnjanovic@40 93 parameterValue = varargin{parameterIndex*2};
idamnjanovic@40 94 switch lower(parameterName)
idamnjanovic@40 95 case 'stoppingcriterion'
idamnjanovic@40 96 stoppingCriterion = parameterValue;
idamnjanovic@40 97 case 'groundtruth'
idamnjanovic@40 98 xG = parameterValue;
idamnjanovic@40 99 case 'tolerance'
idamnjanovic@40 100 tolerance = parameterValue;
idamnjanovic@40 101 case 'linesearchflag'
idamnjanovic@40 102 lineSearchFlag = parameterValue;
idamnjanovic@40 103 case 'lambda'
idamnjanovic@40 104 lambda_bar = parameterValue;
idamnjanovic@40 105 case 'maxiteration'
idamnjanovic@40 106 maxIter = parameterValue;
idamnjanovic@40 107 case 'isnonnegative'
idamnjanovic@40 108 isNonnegative = parameterValue;
idamnjanovic@40 109 case 'continuationflag'
idamnjanovic@40 110 continuationFlag = parameterValue;
idamnjanovic@40 111 case 'initialization'
idamnjanovic@40 112 xk = parameterValue;
idamnjanovic@40 113 if ~all(size(xk)==[n,1])
idamnjanovic@40 114 error('The dimension of the initial xk does not match.');
idamnjanovic@40 115 end
idamnjanovic@40 116 case 'eta'
idamnjanovic@40 117 eta = parameterValue;
idamnjanovic@40 118 if ( eta <= 0 || eta >= 1 )
idamnjanovic@40 119 disp('Line search parameter out of bounds, switching to default 0.9') ;
idamnjanovic@40 120 eta = 0.9 ;
idamnjanovic@40 121 end
idamnjanovic@40 122 otherwise
idamnjanovic@40 123 error(['The parameter ''' parameterName ''' is not recognized by the function ''' mfilename '''.']);
idamnjanovic@40 124 end
idamnjanovic@40 125 end
idamnjanovic@40 126 clear varargin
idamnjanovic@40 127
idamnjanovic@40 128 if stoppingCriterion==STOPPING_GROUND_TRUTH && isempty(xG)
idamnjanovic@40 129 error('The stopping criterion must provide the ground truth value of x.');
idamnjanovic@40 130 end
idamnjanovic@40 131
idamnjanovic@40 132 keep_going = 1 ;
idamnjanovic@40 133 nz_x = (abs(xk)> eps*10);
idamnjanovic@40 134 f = 0.5*norm(b-A*xk)^2 + lambda_bar * norm(xk,1);
idamnjanovic@40 135 xkm1 = xk;
idamnjanovic@40 136 while keep_going && (nIter < maxIter)
idamnjanovic@40 137 nIter = nIter + 1 ;
idamnjanovic@40 138
idamnjanovic@40 139 yk = xk + ((t_km1-1)/t_k)*(xk-xkm1) ;
idamnjanovic@40 140
idamnjanovic@40 141 stop_backtrack = 0 ;
idamnjanovic@40 142
idamnjanovic@40 143 temp = G*yk - c ; % gradient of f at yk
idamnjanovic@40 144
idamnjanovic@40 145 while ~stop_backtrack
idamnjanovic@40 146
idamnjanovic@40 147 gk = yk - (1/L)*temp ;
idamnjanovic@40 148
idamnjanovic@40 149 xkp1 = soft(gk,lambda/L) ;
idamnjanovic@40 150
idamnjanovic@40 151 temp1 = 0.5*norm(b-A*xkp1)^2 ;
idamnjanovic@40 152 temp2 = 0.5*norm(b-A*yk)^2 + (xkp1-yk)'*temp + (L/2)*norm(xkp1-yk)^2 ;
idamnjanovic@40 153
idamnjanovic@40 154 if temp1 <= temp2
idamnjanovic@40 155 stop_backtrack = 1 ;
idamnjanovic@40 156 else
idamnjanovic@40 157 L = L*beta ;
idamnjanovic@40 158 end
idamnjanovic@40 159
idamnjanovic@40 160 end
idamnjanovic@40 161
idamnjanovic@40 162 switch stoppingCriterion
idamnjanovic@40 163 case STOPPING_GROUND_TRUTH
idamnjanovic@40 164 keep_going = norm(xG-xkp1)>tolerance;
idamnjanovic@40 165 case STOPPING_SUBGRADIENT
idamnjanovic@40 166 sk = L*(yk-xkp1) + G*(xkp1-yk) ;
idamnjanovic@40 167 keep_going = norm(sk) > tolerance*L*max(1,norm(xkp1));
idamnjanovic@40 168 case STOPPING_SPARSE_SUPPORT
idamnjanovic@40 169 % compute the stopping criterion based on the change
idamnjanovic@40 170 % of the number of non-zero components of the estimate
idamnjanovic@40 171 nz_x_prev = nz_x;
idamnjanovic@40 172 nz_x = (abs(xkp1)>eps*10);
idamnjanovic@40 173 num_nz_x = sum(nz_x(:));
idamnjanovic@40 174 num_changes_active = (sum(nz_x(:)~=nz_x_prev(:)));
idamnjanovic@40 175 if num_nz_x >= 1
idamnjanovic@40 176 criterionActiveSet = num_changes_active / num_nz_x;
idamnjanovic@40 177 keep_going = (criterionActiveSet > tolerance);
idamnjanovic@40 178 end
idamnjanovic@40 179 case STOPPING_OBJECTIVE_VALUE
idamnjanovic@40 180 % compute the stopping criterion based on the relative
idamnjanovic@40 181 % variation of the objective function.
idamnjanovic@40 182 prev_f = f;
idamnjanovic@40 183 f = 0.5*norm(b-A*xkp1)^2 + lambda_bar * norm(xk,1);
idamnjanovic@40 184 criterionObjective = abs(f-prev_f)/(prev_f);
idamnjanovic@40 185 keep_going = (criterionObjective > tolerance);
idamnjanovic@40 186 case STOPPING_DUALITY_GAP
idamnjanovic@40 187 error('Duality gap is not a valid stopping criterion for PGBP.');
idamnjanovic@40 188 otherwise
idamnjanovic@40 189 error('Undefined stopping criterion.');
idamnjanovic@40 190 end
idamnjanovic@40 191
idamnjanovic@40 192 lambda = max(eta*lambda,lambda_bar) ;
idamnjanovic@40 193
idamnjanovic@40 194
idamnjanovic@40 195 t_kp1 = 0.5*(1+sqrt(1+4*t_k*t_k)) ;
idamnjanovic@40 196
idamnjanovic@40 197 t_km1 = t_k ;
idamnjanovic@40 198 t_k = t_kp1 ;
idamnjanovic@40 199 xkm1 = xk ;
idamnjanovic@40 200 xk = xkp1 ;
idamnjanovic@40 201 end
idamnjanovic@40 202
idamnjanovic@40 203 x_hat = xk ;
idamnjanovic@40 204
idamnjanovic@40 205 function y = soft(x,T)
idamnjanovic@40 206 if sum(abs(T(:)))==0
idamnjanovic@40 207 y = x;
idamnjanovic@40 208 else
idamnjanovic@40 209 y = max(abs(x) - T, 0);
idamnjanovic@40 210 y = sign(x).*y;
idamnjanovic@40 211 end