idamnjanovic@40: % Copyright ©2010. The Regents of the University of California (Regents). idamnjanovic@40: % All Rights Reserved. Contact The Office of Technology Licensing, idamnjanovic@40: % UC Berkeley, 2150 Shattuck Avenue, Suite 510, Berkeley, CA 94720-1620, idamnjanovic@40: % (510) 643-7201, for commercial licensing opportunities. idamnjanovic@40: idamnjanovic@40: % Authors: Arvind Ganesh, Allen Y. Yang, Zihan Zhou. idamnjanovic@40: % Contact: Allen Y. Yang, Department of EECS, University of California, idamnjanovic@40: % Berkeley. idamnjanovic@40: idamnjanovic@40: % IN NO EVENT SHALL REGENTS BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, idamnjanovic@40: % SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, idamnjanovic@40: % ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF idamnjanovic@40: % REGENTS HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. idamnjanovic@40: idamnjanovic@40: % REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED idamnjanovic@40: % TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A idamnjanovic@40: % PARTICULAR PURPOSE. THE SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF ANY, idamnjanovic@40: % PROVIDED HEREUNDER IS PROVIDED "AS IS". REGENTS HAS NO OBLIGATION TO idamnjanovic@40: % PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. idamnjanovic@40: idamnjanovic@40: %% This function is modified from Matlab code proximal_gradient_bp idamnjanovic@40: idamnjanovic@40: function [x_hat,nIter] = SolveFISTA(A,b, varargin) idamnjanovic@40: idamnjanovic@40: % b - m x 1 vector of observations/data (required input) idamnjanovic@40: % A - m x n measurement matrix (required input) idamnjanovic@40: % idamnjanovic@40: % tol - tolerance for stopping criterion. idamnjanovic@40: % - DEFAULT 1e-7 if omitted or -1. idamnjanovic@40: % maxIter - maxilambdam number of iterations idamnjanovic@40: % - DEFAULT 10000, if omitted or -1. idamnjanovic@40: % lineSearchFlag - 1 if line search is to be done every iteration idamnjanovic@40: % - DEFAULT 0, if omitted or -1. idamnjanovic@40: % continuationFlag - 1 if a continuation is to be done on the parameter lambda idamnjanovic@40: % - DEFAULT 1, if omitted or -1. idamnjanovic@40: % eta - line search parameter, should be in (0,1) idamnjanovic@40: % - ignored if lineSearchFlag is 0. idamnjanovic@40: % - DEFAULT 0.9, if omitted or -1. idamnjanovic@40: % lambda - relaxation parameter idamnjanovic@40: % - ignored if continuationFlag is 1. idamnjanovic@40: % - DEFAULT 1e-3, if omitted or -1. idamnjanovic@40: % outputFileName - Details of each iteration are dumped here, if provided. idamnjanovic@40: % idamnjanovic@40: % x_hat - estimate of coeeficient vector idamnjanovic@40: % numIter - number of iterations until convergence idamnjanovic@40: % idamnjanovic@40: % idamnjanovic@40: % References idamnjanovic@40: % "Robust PCA: Exact Recovery of Corrupted Low-Rank Matrices via Convex Optimization", J. Wright et al., preprint 2009. idamnjanovic@40: % "An Accelerated Proximal Gradient Algorithm for Nuclear Norm Regularized Least Squares problems", K.-C. Toh and S. Yun, preprint 2009. idamnjanovic@40: % idamnjanovic@40: % Arvind Ganesh, Summer 2009. Questions? abalasu2@illinois.edu idamnjanovic@40: idamnjanovic@40: DEBUG = 0 ; idamnjanovic@40: idamnjanovic@40: STOPPING_GROUND_TRUTH = -1; idamnjanovic@40: STOPPING_DUALITY_GAP = 1; idamnjanovic@40: STOPPING_SPARSE_SUPPORT = 2; idamnjanovic@40: STOPPING_OBJECTIVE_VALUE = 3; idamnjanovic@40: STOPPING_SUBGRADIENT = 4; idamnjanovic@40: STOPPING_DEFAULT = STOPPING_SUBGRADIENT; idamnjanovic@40: idamnjanovic@40: stoppingCriterion = STOPPING_DEFAULT; idamnjanovic@40: maxIter = 1000 ; idamnjanovic@40: tolerance = 1e-3; idamnjanovic@40: [m,n] = size(A) ; idamnjanovic@40: x0 = zeros(n,1) ; idamnjanovic@40: xG = []; idamnjanovic@40: idamnjanovic@40: %% Initializing optimization variables idamnjanovic@40: t_k = 1 ; idamnjanovic@40: t_km1 = 1 ; idamnjanovic@40: L0 = 1 ; idamnjanovic@40: G = A'*A ; idamnjanovic@40: nIter = 0 ; idamnjanovic@40: c = A'*b ; idamnjanovic@40: lambda0 = 0.99*L0*norm(c,inf) ; idamnjanovic@40: eta = 0.6 ; idamnjanovic@40: lambda_bar = 1e-4*lambda0 ; idamnjanovic@40: xk = zeros(n,1) ; idamnjanovic@40: lambda = lambda0 ; idamnjanovic@40: L = L0 ; idamnjanovic@40: beta = 1.5 ; idamnjanovic@40: idamnjanovic@40: % Parse the optional inputs. idamnjanovic@40: if (mod(length(varargin), 2) ~= 0 ), idamnjanovic@40: error(['Extra Parameters passed to the function ''' mfilename ''' lambdast be passed in pairs.']); idamnjanovic@40: end idamnjanovic@40: parameterCount = length(varargin)/2; idamnjanovic@40: idamnjanovic@40: for parameterIndex = 1:parameterCount, idamnjanovic@40: parameterName = varargin{parameterIndex*2 - 1}; idamnjanovic@40: parameterValue = varargin{parameterIndex*2}; idamnjanovic@40: switch lower(parameterName) idamnjanovic@40: case 'stoppingcriterion' idamnjanovic@40: stoppingCriterion = parameterValue; idamnjanovic@40: case 'groundtruth' idamnjanovic@40: xG = parameterValue; idamnjanovic@40: case 'tolerance' idamnjanovic@40: tolerance = parameterValue; idamnjanovic@40: case 'linesearchflag' idamnjanovic@40: lineSearchFlag = parameterValue; idamnjanovic@40: case 'lambda' idamnjanovic@40: lambda_bar = parameterValue; idamnjanovic@40: case 'maxiteration' idamnjanovic@40: maxIter = parameterValue; idamnjanovic@40: case 'isnonnegative' idamnjanovic@40: isNonnegative = parameterValue; idamnjanovic@40: case 'continuationflag' idamnjanovic@40: continuationFlag = parameterValue; idamnjanovic@40: case 'initialization' idamnjanovic@40: xk = parameterValue; idamnjanovic@40: if ~all(size(xk)==[n,1]) idamnjanovic@40: error('The dimension of the initial xk does not match.'); idamnjanovic@40: end idamnjanovic@40: case 'eta' idamnjanovic@40: eta = parameterValue; idamnjanovic@40: if ( eta <= 0 || eta >= 1 ) idamnjanovic@40: disp('Line search parameter out of bounds, switching to default 0.9') ; idamnjanovic@40: eta = 0.9 ; idamnjanovic@40: end idamnjanovic@40: otherwise idamnjanovic@40: error(['The parameter ''' parameterName ''' is not recognized by the function ''' mfilename '''.']); idamnjanovic@40: end idamnjanovic@40: end idamnjanovic@40: clear varargin idamnjanovic@40: idamnjanovic@40: if stoppingCriterion==STOPPING_GROUND_TRUTH && isempty(xG) idamnjanovic@40: error('The stopping criterion must provide the ground truth value of x.'); idamnjanovic@40: end idamnjanovic@40: idamnjanovic@40: keep_going = 1 ; idamnjanovic@40: nz_x = (abs(xk)> eps*10); idamnjanovic@40: f = 0.5*norm(b-A*xk)^2 + lambda_bar * norm(xk,1); idamnjanovic@40: xkm1 = xk; idamnjanovic@40: while keep_going && (nIter < maxIter) idamnjanovic@40: nIter = nIter + 1 ; idamnjanovic@40: idamnjanovic@40: yk = xk + ((t_km1-1)/t_k)*(xk-xkm1) ; idamnjanovic@40: idamnjanovic@40: stop_backtrack = 0 ; idamnjanovic@40: idamnjanovic@40: temp = G*yk - c ; % gradient of f at yk idamnjanovic@40: idamnjanovic@40: while ~stop_backtrack idamnjanovic@40: idamnjanovic@40: gk = yk - (1/L)*temp ; idamnjanovic@40: idamnjanovic@40: xkp1 = soft(gk,lambda/L) ; idamnjanovic@40: idamnjanovic@40: temp1 = 0.5*norm(b-A*xkp1)^2 ; idamnjanovic@40: temp2 = 0.5*norm(b-A*yk)^2 + (xkp1-yk)'*temp + (L/2)*norm(xkp1-yk)^2 ; idamnjanovic@40: idamnjanovic@40: if temp1 <= temp2 idamnjanovic@40: stop_backtrack = 1 ; idamnjanovic@40: else idamnjanovic@40: L = L*beta ; idamnjanovic@40: end idamnjanovic@40: idamnjanovic@40: end idamnjanovic@40: idamnjanovic@40: switch stoppingCriterion idamnjanovic@40: case STOPPING_GROUND_TRUTH idamnjanovic@40: keep_going = norm(xG-xkp1)>tolerance; idamnjanovic@40: case STOPPING_SUBGRADIENT idamnjanovic@40: sk = L*(yk-xkp1) + G*(xkp1-yk) ; idamnjanovic@40: keep_going = norm(sk) > tolerance*L*max(1,norm(xkp1)); idamnjanovic@40: case STOPPING_SPARSE_SUPPORT idamnjanovic@40: % compute the stopping criterion based on the change idamnjanovic@40: % of the number of non-zero components of the estimate idamnjanovic@40: nz_x_prev = nz_x; idamnjanovic@40: nz_x = (abs(xkp1)>eps*10); idamnjanovic@40: num_nz_x = sum(nz_x(:)); idamnjanovic@40: num_changes_active = (sum(nz_x(:)~=nz_x_prev(:))); idamnjanovic@40: if num_nz_x >= 1 idamnjanovic@40: criterionActiveSet = num_changes_active / num_nz_x; idamnjanovic@40: keep_going = (criterionActiveSet > tolerance); idamnjanovic@40: end idamnjanovic@40: case STOPPING_OBJECTIVE_VALUE idamnjanovic@40: % compute the stopping criterion based on the relative idamnjanovic@40: % variation of the objective function. idamnjanovic@40: prev_f = f; idamnjanovic@40: f = 0.5*norm(b-A*xkp1)^2 + lambda_bar * norm(xk,1); idamnjanovic@40: criterionObjective = abs(f-prev_f)/(prev_f); idamnjanovic@40: keep_going = (criterionObjective > tolerance); idamnjanovic@40: case STOPPING_DUALITY_GAP idamnjanovic@40: error('Duality gap is not a valid stopping criterion for PGBP.'); idamnjanovic@40: otherwise idamnjanovic@40: error('Undefined stopping criterion.'); idamnjanovic@40: end idamnjanovic@40: idamnjanovic@40: lambda = max(eta*lambda,lambda_bar) ; idamnjanovic@40: idamnjanovic@40: idamnjanovic@40: t_kp1 = 0.5*(1+sqrt(1+4*t_k*t_k)) ; idamnjanovic@40: idamnjanovic@40: t_km1 = t_k ; idamnjanovic@40: t_k = t_kp1 ; idamnjanovic@40: xkm1 = xk ; idamnjanovic@40: xk = xkp1 ; idamnjanovic@40: end idamnjanovic@40: idamnjanovic@40: x_hat = xk ; idamnjanovic@40: idamnjanovic@40: function y = soft(x,T) idamnjanovic@40: if sum(abs(T(:)))==0 idamnjanovic@40: y = x; idamnjanovic@40: else idamnjanovic@40: y = max(abs(x) - T, 0); idamnjanovic@40: y = sign(x).*y; idamnjanovic@40: end