[xk,niter,residuals,outputData] =NESTA(A,At,b,muf,delta,opts)
Solves a L1 minimization problem under a quadratic constraint using the
Nesterov algorithm, with continuation:
min_x || U x ||_1 s.t. ||y - Ax||_2 <= delta
Continuation is performed by sequentially applying Nesterov's algorithm
with a decreasing sequence of values of mu0 >= mu >= muf
The primal prox-function is also adapted by accounting for a first guess
xplug that also tends towards x_muf
The observation matrix A is a projector
Inputs: A and At - measurement matrix and adjoint (either a matrix, in which
case At is unused, or function handles). m x n dimensions.
b - Observed data, a m x 1 array
muf - The desired value of mu at the last continuation step.
A smaller mu leads to higher accuracy.
delta - l2 error bound. This enforces how close the variable
must fit the observations b, i.e. || y - Ax ||_2 <= delta
If delta = 0, enforces y = Ax
Common heuristic: delta = sqrt(m + 2*sqrt(2*m))*sigma;
where sigma=std(noise).
opts -
This is a structure that contains additional options,
some of which are optional.
The fieldnames are case insensitive. Below
are the possible fieldnames:
opts.xplug - the first guess for the primal prox-function, and
also the initial point for xk. By default, xplug = At(b)
opts.U and opts.Ut - Analysis/Synthesis operators
(either matrices of function handles).
opts.normU - if opts.U is provided, this should be norm(U)
otherwise it will have to be calculated (potentially
expensive)
opts.MaxIntIter - number of continuation steps.
default is 5
opts.maxiter - max number of iterations in an inner loop.
default is 10,000
opts.TolVar - tolerance for the stopping criteria
opts.stopTest - which stopping criteria to apply
opts.stopTest == 1 : stop when the relative
change in the objective function is less than
TolVar
opts.stopTest == 2 : stop with the l_infinity norm
of difference in the xk variable is less
than TolVar
opts.TypeMin - if this is 'L1' (default), then
minimizes a smoothed version of the l_1 norm.
If this is 'tv', then minimizes a smoothed
version of the total-variation norm.
The string is case insensitive.
opts.Verbose - if this is 0 or false, then very
little output is displayed. If this is 1 or true,
then output every iteration is displayed.
If this is a number p greater than 1, then
output is displayed every pth iteration.
opts.fid - if this is 1 (default), the display is
the usual Matlab screen. If this is the file-id
of a file opened with fopen, then the display
will be redirected to this file.
opts.errFcn - if this is a function handle,
then the program will evaluate opts.errFcn(xk)
at every iteration and display the result.
ex. opts.errFcn = @(x) norm( x - x_true )
opts.outFcn - if this is a function handle,
then then program will evaluate opts.outFcn(xk)
at every iteration and save the results in outputData.
If the result is a vector (as opposed to a scalar),
it should be a row vector and not a column vector.
ex. opts.outFcn = @(x) [norm( x - xtrue, 'inf' ),...
norm( x - xtrue) / norm(xtrue)]
opts.AAtinv - this is an experimental new option. AAtinv
is the inverse of AA^*. This allows the use of a
matrix A which is not a projection, but only
for the noiseless (i.e. delta = 0) case.
opts.USV - another experimental option. This supercedes
the AAtinv option, so it is recommended that you
do not define AAtinv. This allows the use of a matrix
A which is not a projection, and works for the
noisy ( i.e. delta > 0 ) case.
opts.USV should contain three fields:
opts.USV.U is the U from [U,S,V] = svd(A)
likewise, opts.USV.S and opts.USV.V are S and V
from svd(A). S may be a matrix or a vector.
Outputs:
xk - estimate of the solution x
niter - number of iterations
residuals - first column is the residual at every step,
second column is the value of f_mu at every step
outputData - a matrix, where each row r is the output
from opts.outFcn, if supplied.
opts - the structure containing the options that were used
Written by: Jerome Bobin, Caltech
Email: bobin@acm.caltech.edu
Created: February 2009
Modified (version 1.0): May 2009, Jerome Bobin and Stephen Becker, Caltech
Modified (version 1.1): Nov 2009, Stephen Becker, Caltech
NESTA Version 1.1
See also Core_Nesterov0001 function [xk,niter,residuals,outputData,opts] =NESTA_adpW(A,At,b,muf,delta,opts) 0002 % [xk,niter,residuals,outputData] =NESTA(A,At,b,muf,delta,opts) 0003 % 0004 % Solves a L1 minimization problem under a quadratic constraint using the 0005 % Nesterov algorithm, with continuation: 0006 % 0007 % min_x || U x ||_1 s.t. ||y - Ax||_2 <= delta 0008 % 0009 % Continuation is performed by sequentially applying Nesterov's algorithm 0010 % with a decreasing sequence of values of mu0 >= mu >= muf 0011 % 0012 % The primal prox-function is also adapted by accounting for a first guess 0013 % xplug that also tends towards x_muf 0014 % 0015 % The observation matrix A is a projector 0016 % 0017 % Inputs: A and At - measurement matrix and adjoint (either a matrix, in which 0018 % case At is unused, or function handles). m x n dimensions. 0019 % b - Observed data, a m x 1 array 0020 % muf - The desired value of mu at the last continuation step. 0021 % A smaller mu leads to higher accuracy. 0022 % delta - l2 error bound. This enforces how close the variable 0023 % must fit the observations b, i.e. || y - Ax ||_2 <= delta 0024 % If delta = 0, enforces y = Ax 0025 % Common heuristic: delta = sqrt(m + 2*sqrt(2*m))*sigma; 0026 % where sigma=std(noise). 0027 % opts - 0028 % This is a structure that contains additional options, 0029 % some of which are optional. 0030 % The fieldnames are case insensitive. Below 0031 % are the possible fieldnames: 0032 % 0033 % opts.xplug - the first guess for the primal prox-function, and 0034 % also the initial point for xk. By default, xplug = At(b) 0035 % opts.U and opts.Ut - Analysis/Synthesis operators 0036 % (either matrices of function handles). 0037 % opts.normU - if opts.U is provided, this should be norm(U) 0038 % otherwise it will have to be calculated (potentially 0039 % expensive) 0040 % opts.MaxIntIter - number of continuation steps. 0041 % default is 5 0042 % opts.maxiter - max number of iterations in an inner loop. 0043 % default is 10,000 0044 % opts.TolVar - tolerance for the stopping criteria 0045 % opts.stopTest - which stopping criteria to apply 0046 % opts.stopTest == 1 : stop when the relative 0047 % change in the objective function is less than 0048 % TolVar 0049 % opts.stopTest == 2 : stop with the l_infinity norm 0050 % of difference in the xk variable is less 0051 % than TolVar 0052 % opts.TypeMin - if this is 'L1' (default), then 0053 % minimizes a smoothed version of the l_1 norm. 0054 % If this is 'tv', then minimizes a smoothed 0055 % version of the total-variation norm. 0056 % The string is case insensitive. 0057 % opts.Verbose - if this is 0 or false, then very 0058 % little output is displayed. If this is 1 or true, 0059 % then output every iteration is displayed. 0060 % If this is a number p greater than 1, then 0061 % output is displayed every pth iteration. 0062 % opts.fid - if this is 1 (default), the display is 0063 % the usual Matlab screen. If this is the file-id 0064 % of a file opened with fopen, then the display 0065 % will be redirected to this file. 0066 % opts.errFcn - if this is a function handle, 0067 % then the program will evaluate opts.errFcn(xk) 0068 % at every iteration and display the result. 0069 % ex. opts.errFcn = @(x) norm( x - x_true ) 0070 % opts.outFcn - if this is a function handle, 0071 % then then program will evaluate opts.outFcn(xk) 0072 % at every iteration and save the results in outputData. 0073 % If the result is a vector (as opposed to a scalar), 0074 % it should be a row vector and not a column vector. 0075 % ex. opts.outFcn = @(x) [norm( x - xtrue, 'inf' ),... 0076 % norm( x - xtrue) / norm(xtrue)] 0077 % opts.AAtinv - this is an experimental new option. AAtinv 0078 % is the inverse of AA^*. This allows the use of a 0079 % matrix A which is not a projection, but only 0080 % for the noiseless (i.e. delta = 0) case. 0081 % opts.USV - another experimental option. This supercedes 0082 % the AAtinv option, so it is recommended that you 0083 % do not define AAtinv. This allows the use of a matrix 0084 % A which is not a projection, and works for the 0085 % noisy ( i.e. delta > 0 ) case. 0086 % opts.USV should contain three fields: 0087 % opts.USV.U is the U from [U,S,V] = svd(A) 0088 % likewise, opts.USV.S and opts.USV.V are S and V 0089 % from svd(A). S may be a matrix or a vector. 0090 % 0091 % Outputs: 0092 % xk - estimate of the solution x 0093 % niter - number of iterations 0094 % residuals - first column is the residual at every step, 0095 % second column is the value of f_mu at every step 0096 % outputData - a matrix, where each row r is the output 0097 % from opts.outFcn, if supplied. 0098 % opts - the structure containing the options that were used 0099 % 0100 % Written by: Jerome Bobin, Caltech 0101 % Email: bobin@acm.caltech.edu 0102 % Created: February 2009 0103 % Modified (version 1.0): May 2009, Jerome Bobin and Stephen Becker, Caltech 0104 % Modified (version 1.1): Nov 2009, Stephen Becker, Caltech 0105 % 0106 % NESTA Version 1.1 0107 % See also Core_Nesterov 0108 0109 0110 if nargin < 6, opts = []; end 0111 if isempty(opts) && isnumeric(opts), opts = struct; end 0112 0113 %---- Set defaults 0114 fid = setOpts('fid',1); 0115 Verbose = setOpts('Verbose',true); 0116 function printf(varargin), fprintf(fid,varargin{:}); end 0117 MaxIntIter = setOpts('MaxIntIter',5,1); 0118 TypeMin = setOpts('TypeMin','L1'); 0119 TolVar = setOpts('tolvar',1e-5); 0120 [U,U_userSet] = setOpts('U', @(x) x ); 0121 if ~isa(U,'function_handle') 0122 Ut = setOpts('Ut',[]); 0123 else 0124 Ut = setOpts('Ut', @(x) x ); 0125 end 0126 xplug = setOpts('xplug',[]); 0127 normU = setOpts('normU',[]); % so we can tell if it's been set 0128 0129 residuals = []; outputData = []; 0130 AAtinv = setOpts('AAtinv',[]); 0131 USV = setOpts('USV',[]); 0132 if ~isempty(USV) 0133 if isstruct(USV) 0134 Q = USV.U; % we can't use "U" as the variable name 0135 % since "U" already refers to the analysis operator 0136 S = USV.S; 0137 if isvector(S), s = S; %S = diag(s); 0138 else s = diag(S); end 0139 %V = USV.V; 0140 else 0141 error('opts.USV must be a structure'); 0142 end 0143 end 0144 0145 % -- We can handle non-projections IF a (fast) routine for computing 0146 % the psuedo-inverse is available. 0147 % We can handle a nonzero delta, but we need the full SVD 0148 if isempty(AAtinv) && isempty(USV) 0149 % Check if A is a partial isometry, i.e. if AA' = I 0150 z = randn(size(b)); 0151 if isa(A,'function_handle'), AAtz = A(At(z)); 0152 else AAtz = A*(A'*z); end 0153 if norm( AAtz - z )/norm(z) > 1e-8 0154 % error('Measurement matrix A must be a partial isometry: AA''=I'); 0155 end 0156 end 0157 0158 % -- Find a initial guess if not already provided. 0159 % Use least-squares solution: x_ref = A'*inv(A*A')*b 0160 % If A is a projection, the least squares solution is trivial 0161 if isempty(xplug) || norm(xplug) < 1e-12 0162 if ~isempty(USV) && isempty(AAtinv) 0163 AAtinv = Q*diag( s.^(-2) )*Q'; 0164 end 0165 if ~isempty(AAtinv) 0166 if delta > 0 && isempty(USV) 0167 error('delta must be zero for non-projections'); 0168 end 0169 if isa(AAtinv,'function_handle') 0170 x_ref = AAtinv(b); 0171 else 0172 x_ref = AAtinv * b; 0173 end 0174 else 0175 x_ref = b; 0176 end 0177 0178 if isa(A,'function_handle') 0179 x_ref=At(x_ref); 0180 else 0181 x_ref = A'*x_ref; 0182 end 0183 0184 if isempty(xplug) 0185 xplug = x_ref; 0186 end 0187 % x_ref itself is used to calculate mu_0 0188 % in the case that xplug has very small norm 0189 else 0190 x_ref = xplug; 0191 end 0192 0193 % use x_ref, not xplug, to find mu_0 0194 if isa(U,'function_handle') 0195 Ux_ref = U(x_ref); 0196 else 0197 Ux_ref = U*x_ref; 0198 end 0199 switch lower(TypeMin) 0200 case 'l1' 0201 mu0 = 0.9*max(abs(Ux_ref)); 0202 case 'tv' 0203 mu0 = ValMUTv(Ux_ref); 0204 end 0205 0206 % -- If U was set by the user and normU not supplied, then calcuate norm(U) 0207 if U_userSet && isempty(normU) 0208 % simple case: U*U' = I or U'*U = I, in which case norm(U) = 1 0209 z = randn(size(xplug)); 0210 if isa(U,'function_handle'), UtUz = Ut(U(z)); else UtUz = U'*(U*z); end 0211 if norm( UtUz - z )/norm(z) < 1e-8 0212 normU = 1; 0213 else 0214 z = randn(size(Ux_ref)); 0215 if isa(U,'function_handle') 0216 UUtz = U(Ut(z)); 0217 else 0218 UUtz = U*(U'*z); 0219 end 0220 if norm( UUtz - z )/norm(z) < 1e-8 0221 normU = 1; 0222 end 0223 end 0224 0225 if isempty(normU) 0226 % have to actually calculate the norm 0227 if isa(U,'function_handle') 0228 [normU,cnt] = my_normest(U,Ut,length(xplug),1e-3,30); 0229 if cnt == 30, printf('Warning: norm(U) may be inaccurate\n'); end 0230 else 0231 [mU,nU] = size(U); 0232 if mU < nU, UU = U*U'; else UU = U'*U; end 0233 % last resort is to call MATLAB's "norm", which is slow 0234 if norm( UU - diag(diag(UU)),'fro') < 100*eps 0235 % this means the matrix is diagonal, so norm is easy: 0236 normU = sqrt( max(abs(diag(UU))) ); 0237 elseif issparse(UU) 0238 normU = sqrt( normest(UU) ); 0239 else 0240 if min(size(U)) > 2000 0241 % norm(randn(2000)) takes about 5 seconds on my PC 0242 printf('Warning: calculation of norm(U) may be slow\n'); 0243 end 0244 normU = sqrt( norm(UU) ); 0245 end 0246 end 0247 end 0248 opts.normU = normU; 0249 end 0250 0251 0252 niter = 0; 0253 Gamma = (muf/mu0)^(1/MaxIntIter); 0254 mu = mu0; 0255 Gammat= (TolVar/0.1)^(1/MaxIntIter); 0256 TolVar = 0.1; 0257 0258 for nl=1:MaxIntIter 0259 0260 mu = mu*Gamma; 0261 TolVar=TolVar*Gammat; opts.TolVar = TolVar; 0262 opts.xplug = xplug; 0263 if Verbose, printf('\tBeginning %s Minimization; mu = %g\n',opts.TypeMin,mu); end 0264 [xk,niter_int,res,out,optsOut] = Core_Nesterov_adpW(... 0265 A,At,b,mu,delta,opts); 0266 0267 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0268 % Adaptive weights selection [sasif] % 0269 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0270 gamma = find(xk); 0271 ewt_b = 2*length(b)*(norm(xk,2)/norm(xk,1))^2; 0272 W_new = ones(length(xk),1); 0273 W_new(gamma) = min([W_new(gamma) 1./abs(xk(gamma))/ewt_b],[],2); 0274 0275 opts.U = spdiags(W_new,0,length(xk),length(xk)); 0276 opts.Ut = opts.U; 0277 opts.normU = max(W_new); 0278 0279 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0280 0281 xplug = xk; 0282 niter = niter_int + niter; 0283 0284 residuals = [residuals; res]; 0285 outputData = [outputData; out]; 0286 0287 end 0288 opts = optsOut; 0289 0290 0291 %---- internal routine for setting defaults 0292 function [var,userSet] = setOpts(field,default,mn,mx) 0293 var = default; 0294 % has the option already been set? 0295 if ~isfield(opts,field) 0296 % see if there is a capitalization problem: 0297 names = fieldnames(opts); 0298 for i = 1:length(names) 0299 if strcmpi(names{i},field) 0300 opts.(field) = opts.(names{i}); 0301 opts = rmfield(opts,names{i}); 0302 break; 0303 end 0304 end 0305 end 0306 if isfield(opts,field) && ~isempty(opts.(field)) 0307 var = opts.(field); % override the default 0308 userSet = true; 0309 else 0310 userSet = false; 0311 end 0312 % perform error checking, if desired 0313 if nargin >= 3 && ~isempty(mn) 0314 if var < mn 0315 printf('Variable %s is %f, should be at least %f\n',... 0316 field,var,mn); error('variable out-of-bounds'); 0317 end 0318 end 0319 if nargin >= 4 && ~isempty(mx) 0320 if var > mx 0321 printf('Variable %s is %f, should be at least %f\n',... 0322 field,var,mn); error('variable out-of-bounds'); 0323 end 0324 end 0325 opts.(field) = var; 0326 end 0327 0328 0329 0330 0331 %---- internal routine for setting mu0 in the tv minimization case 0332 function th=ValMUTv(x) 0333 0334 global T_frames ROW COL 0335 if isempty(T_frames) 0336 T_frames = 1; 0337 N = length(x)/T_frames; n = floor(sqrt(N)); 0338 ROW = n; COL = n; 0339 end 0340 N = ROW*COL; 0341 Dv = spdiags(repmat([reshape([-ones(ROW-1,COL); zeros(1,COL)],N,1) ... 0342 reshape([zeros(1,COL); ones(ROW-1,COL)],N,1)], T_frames,1), [0 1], T_frames*N, T_frames*N); 0343 Dh = spdiags(repmat([reshape([-ones(ROW,COL-1) zeros(ROW,1)],N,1) ... 0344 reshape([zeros(ROW,1) ones(ROW,COL-1)],N,1)],T_frames,1), [0 ROW], T_frames*N, T_frames*N); 0345 D = sparse([Dh;Dv]); 0346 0347 % Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ... 0348 % reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N); 0349 % Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ... 0350 % reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N); 0351 % D = sparse([Dh;Dv]); 0352 0353 Dhx = Dh*x; 0354 Dvx = Dv*x; 0355 0356 sk = sqrt(abs(Dhx).^2 + abs(Dvx).^2); 0357 th = max(sk); 0358 0359 end 0360 0361 end %-- end of NESTA function 0362 0363 %%%%%%%%%%%% POWER METHOD TO ESTIMATE NORM %%%%%%%%%%%%%%% 0364 % Copied from MATLAB's "normest" function, but allows function handles, not just sparse matrices 0365 function [e,cnt] = my_normest(S,St,n,tol, maxiter) 0366 %MY_NORMEST Estimate the matrix 2-norm via power method. 0367 if nargin < 4, tol = 1.e-6; end 0368 if nargin < 5, maxiter = 20; end 0369 if isempty(St) 0370 St = S; % we assume the matrix is symmetric; 0371 end 0372 x = ones(n,1); 0373 cnt = 0; 0374 e = norm(x); 0375 if e == 0, return, end 0376 x = x/e; 0377 e0 = 0; 0378 while abs(e-e0) > tol*e && cnt < maxiter 0379 e0 = e; 0380 Sx = S(x); 0381 if nnz(Sx) == 0 0382 Sx = rand(size(Sx)); 0383 end 0384 e = norm(Sx); 0385 x = St(Sx); 0386 x = x/norm(x); 0387 cnt = cnt+1; 0388 end 0389 end