[xk,niter,residuals,outputData,opts] =Core_Nesterov(A,At,b,mu,delta,opts)
Solves a L1 minimization problem under a quadratic constraint using the
Nesterov algorithm, without continuation:
min_x || U x ||_1 s.t. ||y - Ax||_2 <= delta
If continuation is desired, see the function NESTA.m
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)
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.
If the SVD of A is U*S*V', then AAtinv = U*(S^{-2})*U'.
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: May 2009, Jerome Bobin and Stephen Becker, Caltech
Modified: Nov 2009, Stephen Becker
NESTA Version 1.1
See also NESTA0001 function [xk,niter,residuals,outputData,opts] = Core_Nesterov(... 0002 A,At,b,mu,delta,opts) 0003 % [xk,niter,residuals,outputData,opts] =Core_Nesterov(A,At,b,mu,delta,opts) 0004 % 0005 % Solves a L1 minimization problem under a quadratic constraint using the 0006 % Nesterov algorithm, without continuation: 0007 % 0008 % min_x || U x ||_1 s.t. ||y - Ax||_2 <= delta 0009 % 0010 % If continuation is desired, see the function NESTA.m 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 % opts.maxiter - max number of iterations in an inner loop. 0039 % default is 10,000 0040 % opts.TolVar - tolerance for the stopping criteria 0041 % opts.stopTest - which stopping criteria to apply 0042 % opts.stopTest == 1 : stop when the relative 0043 % change in the objective function is less than 0044 % TolVar 0045 % opts.stopTest == 2 : stop with the l_infinity norm 0046 % of difference in the xk variable is less 0047 % than TolVar 0048 % opts.TypeMin - if this is 'L1' (default), then 0049 % minimizes a smoothed version of the l_1 norm. 0050 % If this is 'tv', then minimizes a smoothed 0051 % version of the total-variation norm. 0052 % The string is case insensitive. 0053 % opts.Verbose - if this is 0 or false, then very 0054 % little output is displayed. If this is 1 or true, 0055 % then output every iteration is displayed. 0056 % If this is a number p greater than 1, then 0057 % output is displayed every pth iteration. 0058 % opts.fid - if this is 1 (default), the display is 0059 % the usual Matlab screen. If this is the file-id 0060 % of a file opened with fopen, then the display 0061 % will be redirected to this file. 0062 % opts.errFcn - if this is a function handle, 0063 % then the program will evaluate opts.errFcn(xk) 0064 % at every iteration and display the result. 0065 % ex. opts.errFcn = @(x) norm( x - x_true ) 0066 % opts.outFcn - if this is a function handle, 0067 % then then program will evaluate opts.outFcn(xk) 0068 % at every iteration and save the results in outputData. 0069 % If the result is a vector (as opposed to a scalar), 0070 % it should be a row vector and not a column vector. 0071 % ex. opts.outFcn = @(x) [norm( x - xtrue, 'inf' ),... 0072 % norm( x - xtrue) / norm(xtrue)] 0073 % opts.AAtinv - this is an experimental new option. AAtinv 0074 % is the inverse of AA^*. This allows the use of a 0075 % matrix A which is not a projection, but only 0076 % for the noiseless (i.e. delta = 0) case. 0077 % If the SVD of A is U*S*V', then AAtinv = U*(S^{-2})*U'. 0078 % opts.USV - another experimental option. This supercedes 0079 % the AAtinv option, so it is recommended that you 0080 % do not define AAtinv. This allows the use of a matrix 0081 % A which is not a projection, and works for the 0082 % noisy ( i.e. delta > 0 ) case. 0083 % opts.USV should contain three fields: 0084 % opts.USV.U is the U from [U,S,V] = svd(A) 0085 % likewise, opts.USV.S and opts.USV.V are S and V 0086 % from svd(A). S may be a matrix or a vector. 0087 % Outputs: 0088 % xk - estimate of the solution x 0089 % niter - number of iterations 0090 % residuals - first column is the residual at every step, 0091 % second column is the value of f_mu at every step 0092 % outputData - a matrix, where each row r is the output 0093 % from opts.outFcn, if supplied. 0094 % opts - the structure containing the options that were used 0095 % 0096 % Written by: Jerome Bobin, Caltech 0097 % Email: bobin@acm.caltech.edu 0098 % Created: February 2009 0099 % Modified: May 2009, Jerome Bobin and Stephen Becker, Caltech 0100 % Modified: Nov 2009, Stephen Becker 0101 % 0102 % NESTA Version 1.1 0103 % See also NESTA 0104 0105 %---- Set defaults 0106 % opts = []; 0107 fid = setOpts('fid',1); 0108 function printf(varargin), fprintf(fid,varargin{:}); end 0109 maxiter = setOpts('maxiter',10000,0); 0110 TolVar = setOpts('TolVar',1e-5); 0111 TypeMin = setOpts('TypeMin','L1'); 0112 Verbose = setOpts('Verbose',true); 0113 errFcn = setOpts('errFcn',[]); 0114 outFcn = setOpts('outFcn',[]); 0115 stopTest = setOpts('stopTest',1,1,2); 0116 U = setOpts('U', @(x) x ); 0117 if ~isa(U,'function_handle') 0118 Ut = setOpts('Ut',[]); 0119 else 0120 Ut = setOpts('Ut', @(x) x ); 0121 end 0122 xplug = setOpts('xplug',[]); 0123 normU = setOpts('normU',1); 0124 0125 if delta < 0, error('delta must be greater or equal to zero'); end 0126 0127 if isa(A,'function_handle') 0128 Atfun = At; 0129 Afun = A; 0130 else 0131 Atfun = @(x) A'*x; 0132 Afun = @(x) A*x; 0133 end 0134 Atb = Atfun(b); 0135 0136 AAtinv = setOpts('AAtinv',[]); 0137 USV = setOpts('USV',[]); 0138 if ~isempty(USV) 0139 if isstruct(USV) 0140 Q = USV.U; % we can't use "U" as the variable name 0141 % since "U" already refers to the analysis operator 0142 S = USV.S; 0143 if isvector(S), s = S; S = diag(s); 0144 else s = diag(S); end 0145 V = USV.V; 0146 else 0147 error('opts.USV must be a structure'); 0148 end 0149 if isempty(AAtinv) 0150 AAtinv = Q*diag( s.^(-2) )*Q'; 0151 end 0152 end 0153 % --- for A not a projection (experimental) 0154 if ~isempty(AAtinv) 0155 if isa(AAtinv,'function_handle') 0156 AAtinv_fun = AAtinv; 0157 else 0158 AAtinv_fun = @(x) AAtinv * x; 0159 end 0160 0161 AtAAtb = Atfun( AAtinv_fun(b) ); 0162 0163 else 0164 % We assume it's a projection 0165 AtAAtb = Atb; 0166 AAtinv_fun = @(x) x; 0167 end 0168 0169 if isempty(xplug) 0170 xplug = AtAAtb; 0171 end 0172 0173 %---- Initialization 0174 N = length(xplug); 0175 wk = zeros(N,1); 0176 xk = xplug; 0177 0178 0179 %---- Init Variables 0180 Ak= 0; 0181 Lmu = normU/mu; 0182 yk = xk; 0183 zk = xk; 0184 fmean = realmin/10; 0185 OK = 0; 0186 0187 %---- Computing Atb 0188 Atb = Atfun(b); 0189 Axk = Afun(xk);% only needed if you want to see the residuals 0190 % Axplug = Axk; 0191 0192 0193 %---- TV Minimization 0194 if strcmpi(TypeMin,'TV') 0195 Lmu = 8*Lmu; 0196 global T_frames ROW COL 0197 if isempty(T_frames) 0198 T_frames = 1; 0199 N = length(xplug); n = floor(sqrt(N)); 0200 ROW = n; COL = n; 0201 end 0202 0203 % Dv = spdiags(repmat([reshape([-ones(n-1,n); zeros(1,n)],N,1) ... 0204 % reshape([zeros(1,n); ones(n-1,n)],N,1)], T_frames,1), [0 1], T_frames*N, T_frames*N); 0205 % Dh = spdiags(repmat([reshape([-ones(n,n-1) zeros(n,1)],N,1) ... 0206 % reshape([zeros(n,1) ones(n,n-1)],N,1)],T_frames,1), [0 n], T_frames*N, T_frames*N); 0207 % D = sparse([Dh;Dv]); 0208 0209 N = ROW*COL; 0210 Dv = spdiags(repmat([reshape([-ones(ROW-1,COL); zeros(1,COL)],N,1) ... 0211 reshape([zeros(1,COL); ones(ROW-1,COL)],N,1)], T_frames,1), [0 1], T_frames*N, T_frames*N); 0212 Dh = spdiags(repmat([reshape([-ones(ROW,COL-1) zeros(ROW,1)],N,1) ... 0213 reshape([zeros(ROW,1) ones(ROW,COL-1)],N,1)],T_frames,1), [0 ROW], T_frames*N, T_frames*N); 0214 D = sparse([Dh;Dv]); 0215 % Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ... 0216 % reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N); 0217 % Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ... 0218 % reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N); 0219 % D = sparse([Dh;Dv]); 0220 end 0221 0222 0223 Lmu1 = 1/Lmu; 0224 % SLmu = sqrt(Lmu); 0225 % SLmu1 = 1/sqrt(Lmu); 0226 lambdaY = 0; 0227 lambdaZ = 0; 0228 0229 %---- setup data storage variables 0230 [DISPLAY_ERROR, RECORD_DATA] = deal(false); 0231 outputData = deal([]); 0232 residuals = zeros(maxiter,2); 0233 if ~isempty(errFcn), DISPLAY_ERROR = true; end 0234 if ~isempty(outFcn) && nargout >= 4 0235 RECORD_DATA = true; 0236 outputData = zeros(maxiter, size(outFcn(xplug),2) ); 0237 end 0238 0239 for k = 0:maxiter-1, 0240 0241 %---- Dual problem 0242 clear Acp AtAcp df fx uk; 0243 if strcmpi(TypeMin,'L1') [df,fx,val,uk,U,Ut,normU] = Perform_L1_Constraint(xk,mu,U,Ut,length(b));end 0244 0245 if strcmpi(TypeMin,'TV') [df,fx] = Perform_TV_Constraint(xk,mu,Dv,Dh,D,U,Ut);end 0246 0247 %---- Primal Problem 0248 0249 %---- Updating yk 0250 0251 % 0252 % yk = Argmin_x Lmu/2 ||x - xk||_l2^2 + <df,x-xk> s.t. ||b-Ax||_l2 < delta 0253 % Let xp be sqrt(Lmu) (x-xk), dfp be df/sqrt(Lmu), bp be sqrt(Lmu)(b- Axk) and deltap be sqrt(Lmu)delta 0254 % yk = xk + 1/sqrt(Lmu) Argmin_xp 1/2 || xp ||_2^2 + <dfp,xp> s.t. || bp - Axp ||_2 < deltap 0255 % 0256 0257 0258 cp = xk - 1/Lmu*df; % this is "q" in eq. (3.7) in the paper 0259 0260 Acp = Afun( cp ); 0261 if ~isempty(AAtinv) && isempty(USV) 0262 AtAcp = Atfun( AAtinv_fun( Acp ) ); 0263 else 0264 AtAcp = Atfun( Acp ); 0265 end 0266 0267 residuals(k+1,1) = norm( b-Axk); % the residual 0268 residuals(k+1,2) = fx; % the value of the objective 0269 %--- if user has supplied a function, apply it to the iterate 0270 if RECORD_DATA 0271 outputData(k+1,:) = outFcn(xk); 0272 end 0273 0274 if delta > 0 0275 if ~isempty(USV) 0276 % there are more efficient methods, but we're assuming 0277 % that A is negligible compared to U and Ut. 0278 % Here we make the change of variables x <-- x - xk 0279 % and df <-- df/L 0280 dfp = -Lmu1*df; Adfp = -(Axk - Acp); 0281 bp = b - Axk; 0282 deltap = delta; 0283 % Check if we even need to project: 0284 if norm( Adfp - bp ) < deltap 0285 lambdaY = 0; projIter = 0; 0286 % i.e. projection = dfp; 0287 yk = xk + dfp; 0288 Ayk = Axk + Adfp; 0289 else 0290 lambdaY_old = lambdaY; 0291 [projection,projIter,lambdaY] = fastProjection(Q,S,V,dfp,bp,... 0292 deltap, .999*lambdaY_old ); 0293 if lambdaY > 0, disp('lambda is positive!'); keyboard; end 0294 yk = xk + projection; 0295 Ayk = Afun(yk); 0296 % DEBUGGING 0297 % if projIter == 50 0298 % fprintf('\n Maxed out iterations at y\n'); 0299 % keyboard 0300 % end 0301 end 0302 else 0303 lambda = max(0,Lmu*(norm(b-Acp)/delta - 1));gamma = lambda/(lambda + Lmu); 0304 yk = lambda/Lmu*(1-gamma)*Atb + cp - gamma*AtAcp; 0305 % for calculating the residual, we'll avoid calling A() 0306 % by storing A(yk) here (using A'*A = I): 0307 Ayk = lambda/Lmu*(1-gamma)*b + Acp - gamma*Acp; 0308 end 0309 else 0310 % if delta is 0, the projection is simplified: 0311 yk = AtAAtb + cp - AtAcp; 0312 Ayk = b; 0313 end 0314 0315 % DEBUGGING 0316 % if norm( Ayk - b ) > (1.05)*delta 0317 % fprintf('\nAyk failed projection test\n'); 0318 % keyboard; 0319 % end 0320 0321 %--- Stopping criterion 0322 qp = abs(fx - mean(fmean))/mean(fmean); 0323 0324 switch stopTest 0325 case 1 0326 % look at the relative change in function value 0327 if qp <= TolVar && OK; break;end 0328 if qp <= TolVar && ~OK; OK=1; end 0329 case 2 0330 % look at the l_inf change from previous iterate 0331 if k >= 1 && norm( xk - xold, 'inf' ) <= TolVar 0332 break 0333 end 0334 end 0335 fmean = [fx,fmean]; 0336 if (length(fmean) > 10) fmean = fmean(1:10);end 0337 0338 0339 0340 %--- Updating zk 0341 0342 apk =0.5*(k+1); 0343 Ak = Ak + apk; 0344 tauk = 2/(k+3); 0345 0346 wk = apk*df + wk; 0347 0348 % 0349 % zk = Argmin_x Lmu/2 ||b - Ax||_l2^2 + Lmu/2||x - xplug ||_2^2+ <wk,x-xk> 0350 % s.t. ||b-Ax||_l2 < delta 0351 % 0352 0353 cp = xplug - 1/Lmu*wk; 0354 0355 Acp = Afun( cp ); 0356 if ~isempty( AAtinv ) && isempty(USV) 0357 AtAcp = Atfun( AAtinv_fun( Acp ) ); 0358 else 0359 AtAcp = Atfun( Acp ); 0360 end 0361 0362 if delta > 0 0363 if ~isempty(USV) 0364 % Make the substitution wk <-- wk/K 0365 0366 % dfp = (xplug - Lmu1*wk); % = cp 0367 % Adfp= (Axplug - Acp); 0368 dfp = cp; Adfp = Acp; 0369 bp = b; 0370 deltap = delta; 0371 % dfp = SLmu*xplug - SLmu1*wk; 0372 % bp = SLmu*b; 0373 % deltap = SLmu*delta; 0374 0375 % See if we even need to project: 0376 if norm( Adfp - bp ) < deltap 0377 zk = dfp; 0378 Azk = Adfp; 0379 else 0380 [projection,projIter,lambdaZ] = fastProjection(Q,S,V,dfp,bp,... 0381 deltap, .999*lambdaZ ); 0382 if lambdaZ > 0, disp('lambda is positive!'); keyboard; end 0383 zk = projection; 0384 % zk = SLmu1*projection; 0385 Azk = Afun(zk); 0386 0387 % DEBUGGING: 0388 % if projIter == 50 0389 % fprintf('\n Maxed out iterations at z\n'); 0390 % keyboard 0391 % end 0392 end 0393 else 0394 lambda = max(0,Lmu*(norm(b-Acp)/delta - 1));gamma = lambda/(lambda + Lmu); 0395 zk = lambda/Lmu*(1-gamma)*Atb + cp - gamma*AtAcp; 0396 % for calculating the residual, we'll avoid calling A() 0397 % by storing A(zk) here (using A'*A = I): 0398 Azk = lambda/Lmu*(1-gamma)*b + Acp - gamma*Acp; 0399 end 0400 else 0401 % if delta is 0, this is simplified: 0402 zk = AtAAtb + cp - AtAcp; 0403 Azk = b; 0404 end 0405 0406 % DEBUGGING 0407 % if norm( Ayk - b ) > (1.05)*delta 0408 % fprintf('\nAzk failed projection test\n'); 0409 % keyboard; 0410 % end 0411 0412 %--- Updating xk 0413 0414 xkp = tauk*zk + (1-tauk)*yk; 0415 xold = xk; 0416 xk=xkp; 0417 Axk = tauk*Azk + (1-tauk)*Ayk; 0418 0419 if ~mod(k,10), Axk = Afun(xk); end % otherwise slowly lose precision 0420 % DEBUG 0421 % if norm(Axk - Afun(xk) ) > 1e-6, disp('error with Axk'); keyboard; end 0422 0423 %--- display progress if desired 0424 if ~mod(k+1,Verbose ) 0425 printf('Iter: %3d ~ fmu: %.3e ~ Rel. Variation of fmu: %.2e ~ Residual: %.2e',... 0426 k+1,fx,qp,residuals(k+1,1) ); 0427 %--- if user has supplied a function to calculate the error, 0428 % apply it to the current iterate and dislay the output: 0429 if DISPLAY_ERROR, printf(' ~ Error: %.2e',errFcn(xk)); end 0430 printf('\n'); 0431 0432 0433 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0434 end 0435 if abs(fx)>1e20 || abs(residuals(k+1,1)) >1e20 || isnan(fx) 0436 error('Nesta: possible divergence or NaN. Bad estimate of ||A''A||?'); 0437 end 0438 0439 end 0440 0441 niter = k+1; 0442 0443 %-- truncate output vectors 0444 residuals = residuals(1:niter,:); 0445 if RECORD_DATA, outputData = outputData(1:niter,:); end 0446 0447 0448 0449 %---- internal routine for setting defaults 0450 function [var,userSet] = setOpts(field,default,mn,mx) 0451 var = default; 0452 % has the option already been set? 0453 if ~isfield(opts,field) 0454 % see if there is a capitalization problem: 0455 names = fieldnames(opts); 0456 for i = 1:length(names) 0457 if strcmpi(names{i},field) 0458 opts.(field) = opts.(names{i}); 0459 opts = rmfield(opts,names{i}); 0460 break; 0461 end 0462 end 0463 end 0464 0465 if isfield(opts,field) && ~isempty(opts.(field)) 0466 var = opts.(field); % override the default 0467 userSet = true; 0468 else 0469 userSet = false; 0470 end 0471 0472 % perform error checking, if desired 0473 if nargin >= 3 && ~isempty(mn) 0474 if var < mn 0475 printf('Variable %s is %f, should be at least %f\n',... 0476 field,var,mn); error('variable out-of-bounds'); 0477 end 0478 end 0479 if nargin >= 4 && ~isempty(mx) 0480 if var > mx 0481 printf('Variable %s is %f, should be at least %f\n',... 0482 field,var,mn); error('variable out-of-bounds'); 0483 end 0484 end 0485 opts.(field) = var; 0486 end 0487 0488 0489 end %% end of main Core_Nesterov routine 0490 0491 0492 %%%%%%%%%%%% PERFORM THE L1 CONSTRAINT %%%%%%%%%%%%%%%%%% 0493 0494 function [df,fx,val,uk, U, Ut, normU] = Perform_L1_Constraint(xk,mu,U,Ut,M) 0495 0496 if isa(U,'function_handle') 0497 uk = U(xk); 0498 else 0499 uk = U*xk; 0500 end 0501 fx = uk; 0502 0503 uk = uk./max(mu,abs(uk)); 0504 val = real(uk'*fx); 0505 fx = real(uk'*fx - mu/2*norm(uk)^2); 0506 0507 if isa(Ut,'function_handle') 0508 df = Ut(uk); 0509 else 0510 df = U'*uk; 0511 end 0512 0513 gamma = find(xk); 0514 ewt_b = 2*M*(norm(xk,2)/norm(xk,1))^2; 0515 W_new = ones(length(xk),1); 0516 W_new(gamma) = min([W_new(gamma) 1./abs(xk(gamma))/ewt_b],[],2); 0517 0518 U = U; %spdiags(W_new,0,length(xk),length(xk)); 0519 Ut = U; 0520 normU = max(W_new); 0521 end 0522 0523 %%%%%%%%%%%% PERFORM THE TV CONSTRAINT %%%%%%%%%%%%%%%%%% 0524 0525 function [df,fx] = Perform_TV_Constraint(xk,mu,Dv,Dh,D,U,Ut) 0526 if isa(U,'function_handle') 0527 x = U(xk); 0528 else 0529 x = U*xk; 0530 end 0531 df = zeros(size(x)); 0532 0533 Dhx = Dh*x; 0534 Dvx = Dv*x; 0535 0536 tvx = sum(sqrt(abs(Dhx).^2+abs(Dvx).^2)); 0537 w = max(mu,sqrt(abs(Dhx).^2 + abs(Dvx).^2)); 0538 uh = Dhx ./ w; 0539 uv = Dvx ./ w; 0540 u = [uh;uv]; 0541 fx = real(u'*D*x - mu/2 * 1/numel(u)*sum(u'*u)); 0542 if isa(Ut,'function_handle') 0543 df = Ut(D'*u); 0544 else 0545 df = U'*(D'*u); 0546 end 0547 end 0548