nikcleju@45: function [xk,niter,residuals,outputData,opts] =NESTA(A,At,b,muf,delta,opts) nikcleju@45: % [xk,niter,residuals,outputData] =NESTA(A,At,b,muf,delta,opts) nikcleju@45: % nikcleju@45: % Solves a L1 minimization problem under a quadratic constraint using the nikcleju@45: % Nesterov algorithm, with continuation: nikcleju@45: % nikcleju@45: % min_x || U x ||_1 s.t. ||y - Ax||_2 <= delta nikcleju@45: % nikcleju@45: % Continuation is performed by sequentially applying Nesterov's algorithm nikcleju@45: % with a decreasing sequence of values of mu0 >= mu >= muf nikcleju@45: % nikcleju@45: % The primal prox-function is also adapted by accounting for a first guess nikcleju@45: % xplug that also tends towards x_muf nikcleju@45: % nikcleju@45: % The observation matrix A is a projector nikcleju@45: % nikcleju@45: % Inputs: A and At - measurement matrix and adjoint (either a matrix, in which nikcleju@45: % case At is unused, or function handles). m x n dimensions. nikcleju@45: % b - Observed data, a m x 1 array nikcleju@45: % muf - The desired value of mu at the last continuation step. nikcleju@45: % A smaller mu leads to higher accuracy. nikcleju@45: % delta - l2 error bound. This enforces how close the variable nikcleju@45: % must fit the observations b, i.e. || y - Ax ||_2 <= delta nikcleju@45: % If delta = 0, enforces y = Ax nikcleju@45: % Common heuristic: delta = sqrt(m + 2*sqrt(2*m))*sigma; nikcleju@45: % where sigma=std(noise). nikcleju@45: % opts - nikcleju@45: % This is a structure that contains additional options, nikcleju@45: % some of which are optional. nikcleju@45: % The fieldnames are case insensitive. Below nikcleju@45: % are the possible fieldnames: nikcleju@45: % nikcleju@45: % opts.xplug - the first guess for the primal prox-function, and nikcleju@45: % also the initial point for xk. By default, xplug = At(b) nikcleju@45: % opts.U and opts.Ut - Analysis/Synthesis operators nikcleju@45: % (either matrices of function handles). nikcleju@45: % opts.normU - if opts.U is provided, this should be norm(U) nikcleju@45: % otherwise it will have to be calculated (potentially nikcleju@45: % expensive) nikcleju@45: % opts.MaxIntIter - number of continuation steps. nikcleju@45: % default is 5 nikcleju@45: % opts.maxiter - max number of iterations in an inner loop. nikcleju@45: % default is 10,000 nikcleju@45: % opts.TolVar - tolerance for the stopping criteria nikcleju@45: % opts.stopTest - which stopping criteria to apply nikcleju@45: % opts.stopTest == 1 : stop when the relative nikcleju@45: % change in the objective function is less than nikcleju@45: % TolVar nikcleju@45: % opts.stopTest == 2 : stop with the l_infinity norm nikcleju@45: % of difference in the xk variable is less nikcleju@45: % than TolVar nikcleju@45: % opts.TypeMin - if this is 'L1' (default), then nikcleju@45: % minimizes a smoothed version of the l_1 norm. nikcleju@45: % If this is 'tv', then minimizes a smoothed nikcleju@45: % version of the total-variation norm. nikcleju@45: % The string is case insensitive. nikcleju@45: % opts.Verbose - if this is 0 or false, then very nikcleju@45: % little output is displayed. If this is 1 or true, nikcleju@45: % then output every iteration is displayed. nikcleju@45: % If this is a number p greater than 1, then nikcleju@45: % output is displayed every pth iteration. nikcleju@45: % opts.fid - if this is 1 (default), the display is nikcleju@45: % the usual Matlab screen. If this is the file-id nikcleju@45: % of a file opened with fopen, then the display nikcleju@45: % will be redirected to this file. nikcleju@45: % opts.errFcn - if this is a function handle, nikcleju@45: % then the program will evaluate opts.errFcn(xk) nikcleju@45: % at every iteration and display the result. nikcleju@45: % ex. opts.errFcn = @(x) norm( x - x_true ) nikcleju@45: % opts.outFcn - if this is a function handle, nikcleju@45: % then then program will evaluate opts.outFcn(xk) nikcleju@45: % at every iteration and save the results in outputData. nikcleju@45: % If the result is a vector (as opposed to a scalar), nikcleju@45: % it should be a row vector and not a column vector. nikcleju@45: % ex. opts.outFcn = @(x) [norm( x - xtrue, 'inf' ),... nikcleju@45: % norm( x - xtrue) / norm(xtrue)] nikcleju@45: % opts.AAtinv - this is an experimental new option. AAtinv nikcleju@45: % is the inverse of AA^*. This allows the use of a nikcleju@45: % matrix A which is not a projection, but only nikcleju@45: % for the noiseless (i.e. delta = 0) case. nikcleju@45: % opts.USV - another experimental option. This supercedes nikcleju@45: % the AAtinv option, so it is recommended that you nikcleju@45: % do not define AAtinv. This allows the use of a matrix nikcleju@45: % A which is not a projection, and works for the nikcleju@45: % noisy ( i.e. delta > 0 ) case. nikcleju@45: % opts.USV should contain three fields: nikcleju@45: % opts.USV.U is the U from [U,S,V] = svd(A) nikcleju@45: % likewise, opts.USV.S and opts.USV.V are S and V nikcleju@45: % from svd(A). S may be a matrix or a vector. nikcleju@45: % nikcleju@45: % Outputs: nikcleju@45: % xk - estimate of the solution x nikcleju@45: % niter - number of iterations nikcleju@45: % residuals - first column is the residual at every step, nikcleju@45: % second column is the value of f_mu at every step nikcleju@45: % outputData - a matrix, where each row r is the output nikcleju@45: % from opts.outFcn, if supplied. nikcleju@45: % opts - the structure containing the options that were used nikcleju@45: % nikcleju@45: % Written by: Jerome Bobin, Caltech nikcleju@45: % Email: bobin@acm.caltech.edu nikcleju@45: % Created: February 2009 nikcleju@45: % Modified (version 1.0): May 2009, Jerome Bobin and Stephen Becker, Caltech nikcleju@45: % Modified (version 1.1): Nov 2009, Stephen Becker, Caltech nikcleju@45: % nikcleju@45: % NESTA Version 1.1 nikcleju@45: % See also Core_Nesterov nikcleju@45: nikcleju@45: nikcleju@45: if nargin < 6, opts = []; end nikcleju@45: if isempty(opts) && isnumeric(opts), opts = struct; end nikcleju@45: nikcleju@45: %---- Set defaults nikcleju@45: fid = setOpts('fid',1); nikcleju@45: Verbose = setOpts('Verbose',true); nikcleju@45: function printf(varargin), fprintf(fid,varargin{:}); end nikcleju@45: MaxIntIter = setOpts('MaxIntIter',5,1); nikcleju@45: TypeMin = setOpts('TypeMin','L1'); nikcleju@45: TolVar = setOpts('tolvar',1e-5); nikcleju@45: [U,U_userSet] = setOpts('U', @(x) x ); nikcleju@45: if ~isa(U,'function_handle') nikcleju@45: Ut = setOpts('Ut',[]); nikcleju@45: else nikcleju@45: Ut = setOpts('Ut', @(x) x ); nikcleju@45: end nikcleju@45: xplug = setOpts('xplug',[]); nikcleju@45: normU = setOpts('normU',[]); % so we can tell if it's been set nikcleju@45: nikcleju@45: residuals = []; outputData = []; nikcleju@45: AAtinv = setOpts('AAtinv',[]); nikcleju@45: USV = setOpts('USV',[]); nikcleju@45: if ~isempty(USV) nikcleju@45: if isstruct(USV) nikcleju@45: Q = USV.U; % we can't use "U" as the variable name nikcleju@45: % since "U" already refers to the analysis operator nikcleju@45: S = USV.S; nikcleju@45: if isvector(S), s = S; %S = diag(s); nikcleju@45: else s = diag(S); end nikcleju@45: %V = USV.V; nikcleju@45: else nikcleju@45: error('opts.USV must be a structure'); nikcleju@45: end nikcleju@45: end nikcleju@45: nikcleju@45: % -- We can handle non-projections IF a (fast) routine for computing nikcleju@45: % the psuedo-inverse is available. nikcleju@45: % We can handle a nonzero delta, but we need the full SVD nikcleju@45: if isempty(AAtinv) && isempty(USV) nikcleju@45: % Check if A is a partial isometry, i.e. if AA' = I nikcleju@45: z = randn(size(b)); nikcleju@45: if isa(A,'function_handle'), AAtz = A(At(z)); nikcleju@45: else AAtz = A*(A'*z); end nikcleju@45: if norm( AAtz - z )/norm(z) > 1e-8 nikcleju@45: error('Measurement matrix A must be a partial isometry: AA''=I'); nikcleju@45: end nikcleju@45: end nikcleju@45: nikcleju@45: % -- Find a initial guess if not already provided. nikcleju@45: % Use least-squares solution: x_ref = A'*inv(A*A')*b nikcleju@45: % If A is a projection, the least squares solution is trivial nikcleju@45: if isempty(xplug) || norm(xplug) < 1e-12 nikcleju@45: if ~isempty(USV) && isempty(AAtinv) nikcleju@45: AAtinv = Q*diag( s.^(-2) )*Q'; nikcleju@45: end nikcleju@45: if ~isempty(AAtinv) nikcleju@45: if delta > 0 && isempty(USV) nikcleju@45: error('delta must be zero for non-projections'); nikcleju@45: end nikcleju@45: if isa(AAtinv,'function_handle') nikcleju@45: x_ref = AAtinv(b); nikcleju@45: else nikcleju@45: x_ref = AAtinv * b; nikcleju@45: end nikcleju@45: else nikcleju@45: x_ref = b; nikcleju@45: end nikcleju@45: nikcleju@45: if isa(A,'function_handle') nikcleju@45: x_ref=At(x_ref); nikcleju@45: else nikcleju@45: x_ref = A'*x_ref; nikcleju@45: end nikcleju@45: nikcleju@45: if isempty(xplug) nikcleju@45: xplug = x_ref; nikcleju@45: end nikcleju@45: % x_ref itself is used to calculate mu_0 nikcleju@45: % in the case that xplug has very small norm nikcleju@45: else nikcleju@45: x_ref = xplug; nikcleju@45: end nikcleju@45: nikcleju@45: % use x_ref, not xplug, to find mu_0 nikcleju@45: if isa(U,'function_handle') nikcleju@45: Ux_ref = U(x_ref); nikcleju@45: else nikcleju@45: Ux_ref = U*x_ref; nikcleju@45: end nikcleju@45: switch lower(TypeMin) nikcleju@45: case 'l1' nikcleju@45: mu0 = 0.9*max(abs(Ux_ref)); nikcleju@45: case 'tv' nikcleju@45: mu0 = ValMUTv(Ux_ref); nikcleju@45: end nikcleju@45: nikcleju@45: % -- If U was set by the user and normU not supplied, then calcuate norm(U) nikcleju@45: if U_userSet && isempty(normU) nikcleju@45: % simple case: U*U' = I or U'*U = I, in which case norm(U) = 1 nikcleju@45: z = randn(size(xplug)); nikcleju@45: if isa(U,'function_handle'), UtUz = Ut(U(z)); else UtUz = U'*(U*z); end nikcleju@45: if norm( UtUz - z )/norm(z) < 1e-8 nikcleju@45: normU = 1; nikcleju@45: else nikcleju@45: z = randn(size(Ux_ref)); nikcleju@45: if isa(U,'function_handle') nikcleju@45: UUtz = U(Ut(z)); nikcleju@45: else nikcleju@45: UUtz = U*(U'*z); nikcleju@45: end nikcleju@45: if norm( UUtz - z )/norm(z) < 1e-8 nikcleju@45: normU = 1; nikcleju@45: end nikcleju@45: end nikcleju@45: nikcleju@45: if isempty(normU) nikcleju@45: % have to actually calculate the norm nikcleju@45: if isa(U,'function_handle') nikcleju@45: [normU,cnt] = my_normest(U,Ut,length(xplug),1e-3,30); nikcleju@45: if cnt == 30, printf('Warning: norm(U) may be inaccurate\n'); end nikcleju@45: else nikcleju@45: [mU,nU] = size(U); nikcleju@45: if mU < nU, UU = U*U'; else UU = U'*U; end nikcleju@45: % last resort is to call MATLAB's "norm", which is slow nikcleju@45: if norm( UU - diag(diag(UU)),'fro') < 100*eps nikcleju@45: % this means the matrix is diagonal, so norm is easy: nikcleju@45: normU = sqrt( max(abs(diag(UU))) ); nikcleju@45: elseif issparse(UU) nikcleju@45: normU = sqrt( normest(UU) ); nikcleju@45: else nikcleju@45: if min(size(U)) > 2000 nikcleju@45: % norm(randn(2000)) takes about 5 seconds on my PC nikcleju@45: printf('Warning: calculation of norm(U) may be slow\n'); nikcleju@45: end nikcleju@45: normU = sqrt( norm(UU) ); nikcleju@45: end nikcleju@45: end nikcleju@45: end nikcleju@45: opts.normU = normU; nikcleju@45: end nikcleju@45: nikcleju@45: nikcleju@45: niter = 0; nikcleju@45: Gamma = (muf/mu0)^(1/MaxIntIter); nikcleju@45: mu = mu0; nikcleju@45: Gammat= (TolVar/0.1)^(1/MaxIntIter); nikcleju@45: TolVar = 0.1; nikcleju@45: nikcleju@45: for nl=1:MaxIntIter nikcleju@45: nikcleju@45: mu = mu*Gamma; nikcleju@45: TolVar=TolVar*Gammat; opts.TolVar = TolVar; nikcleju@45: opts.xplug = xplug; nikcleju@45: if Verbose, printf('\tBeginning %s Minimization; mu = %g\n',opts.TypeMin,mu); end nikcleju@45: [xk,niter_int,res,out,optsOut] = Core_Nesterov(... nikcleju@45: A,At,b,mu,delta,opts); nikcleju@45: nikcleju@45: xplug = xk; nikcleju@45: niter = niter_int + niter; nikcleju@45: nikcleju@45: residuals = [residuals; res]; nikcleju@45: outputData = [outputData; out]; nikcleju@45: nikcleju@45: end nikcleju@45: opts = optsOut; nikcleju@45: nikcleju@45: nikcleju@45: %---- internal routine for setting defaults nikcleju@45: function [var,userSet] = setOpts(field,default,mn,mx) nikcleju@45: var = default; nikcleju@45: % has the option already been set? nikcleju@45: if ~isfield(opts,field) nikcleju@45: % see if there is a capitalization problem: nikcleju@45: names = fieldnames(opts); nikcleju@45: for i = 1:length(names) nikcleju@45: if strcmpi(names{i},field) nikcleju@45: opts.(field) = opts.(names{i}); nikcleju@45: opts = rmfield(opts,names{i}); nikcleju@45: break; nikcleju@45: end nikcleju@45: end nikcleju@45: end nikcleju@45: if isfield(opts,field) && ~isempty(opts.(field)) nikcleju@45: var = opts.(field); % override the default nikcleju@45: userSet = true; nikcleju@45: else nikcleju@45: userSet = false; nikcleju@45: end nikcleju@45: % perform error checking, if desired nikcleju@45: if nargin >= 3 && ~isempty(mn) nikcleju@45: if var < mn nikcleju@45: printf('Variable %s is %f, should be at least %f\n',... nikcleju@45: field,var,mn); error('variable out-of-bounds'); nikcleju@45: end nikcleju@45: end nikcleju@45: if nargin >= 4 && ~isempty(mx) nikcleju@45: if var > mx nikcleju@45: printf('Variable %s is %f, should be at least %f\n',... nikcleju@45: field,var,mn); error('variable out-of-bounds'); nikcleju@45: end nikcleju@45: end nikcleju@45: opts.(field) = var; nikcleju@45: end nikcleju@45: nikcleju@45: nikcleju@45: nikcleju@45: nikcleju@45: %---- internal routine for setting mu0 in the tv minimization case nikcleju@45: function th=ValMUTv(x) nikcleju@45: nikcleju@45: N = length(x);n = floor(sqrt(N)); nikcleju@45: Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ... nikcleju@45: reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N); nikcleju@45: Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ... nikcleju@45: reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N); nikcleju@45: D = sparse([Dh;Dv]); nikcleju@45: nikcleju@45: nikcleju@45: Dhx = Dh*x; nikcleju@45: Dvx = Dv*x; nikcleju@45: nikcleju@45: sk = sqrt(abs(Dhx).^2 + abs(Dvx).^2); nikcleju@45: th = max(sk); nikcleju@45: nikcleju@45: end nikcleju@45: nikcleju@45: end %-- end of NESTA function nikcleju@45: nikcleju@45: %%%%%%%%%%%% POWER METHOD TO ESTIMATE NORM %%%%%%%%%%%%%%% nikcleju@45: % Copied from MATLAB's "normest" function, but allows function handles, not just sparse matrices nikcleju@45: function [e,cnt] = my_normest(S,St,n,tol, maxiter) nikcleju@45: %MY_NORMEST Estimate the matrix 2-norm via power method. nikcleju@45: if nargin < 4, tol = 1.e-6; end nikcleju@45: if nargin < 5, maxiter = 20; end nikcleju@45: if isempty(St) nikcleju@45: St = S; % we assume the matrix is symmetric; nikcleju@45: end nikcleju@45: x = ones(n,1); nikcleju@45: cnt = 0; nikcleju@45: e = norm(x); nikcleju@45: if e == 0, return, end nikcleju@45: x = x/e; nikcleju@45: e0 = 0; nikcleju@45: while abs(e-e0) > tol*e && cnt < maxiter nikcleju@45: e0 = e; nikcleju@45: Sx = S(x); nikcleju@45: if nnz(Sx) == 0 nikcleju@45: Sx = rand(size(Sx)); nikcleju@45: end nikcleju@45: e = norm(Sx); nikcleju@45: x = St(Sx); nikcleju@45: x = x/norm(x); nikcleju@45: cnt = cnt+1; nikcleju@45: end nikcleju@45: end