SpaRSA version 2.0, December 31, 2007
This function solves the convex problem
arg min_x = 0.5*|| y - A x ||_2^2 + tau phi(x)
using the SpaRSA algorithm, which is described in "Sparse Reconstruction
by Separable Approximation" by S. Wright, R. Nowak, M. Figueiredo,
IEEE Transactions on Signal Processing, 2009 (to appear).
The algorithm is related GPSR (Figueiredo, Nowak, Wright) but does not
rely on the conversion to QP form of the l1 norm, because it is not
limited to being used with l1 regularization. Instead it forms a separable
approximation to the first term of the objective, which has the form
d'*A'*(A x - y) + 0.5*alpha*d'*d
where alpha is obtained from a BB formula. In a monotone variant, alpha is
increased until we see a decreasein the original objective function over
this step.
-----------------------------------------------------------------------
Copyright (2007): Mario Figueiredo, Robert Nowak, Stephen Wright
GPSR is distributed under the terms
of the GNU General Public License 2.0.
Permission to use, copy, modify, and distribute this software for
any purpose without fee is hereby granted, provided that this entire
notice is included in all copies of any software which is or includes
a copy or modification of this software and in all copies of the
supporting documentation for such software.
This software is being provided "as is", without any express or
implied warranty. In particular, the authors do not make any
representation or warranty of any kind concerning the merchantability
of this software or its fitness for any particular purpose."
----------------------------------------------------------------------
Please check for the latest version of the code and paper at
www.lx.it.pt/~mtf/SpaRSA
===== Required inputs =============
y: 1D vector or 2D array (image) of observations
A: if y and x are both 1D vectors, A can be a
k*n (where k is the size of y and n the size of x)
matrix or a handle to a function that computes
products of the form A*v, for some vector v.
In any other case (if y and/or x are 2D arrays),
A has to be passed as a handle to a function which computes
products of the form A*x; another handle to a function
AT which computes products of the form A'*x is also required
in this case. The size of x is determined as the size
of the result of applying AT.
tau: regularization parameter (scalar)
===== Optional inputs =============
'AT' = function handle for the function that implements
the multiplication by the conjugate of A, when A
is a function handle. If A is an array, AT is ignored.
'Psi' = handle to the denoising function, that is, to a function
that computes the solution of the densoing probelm
corresponding to the desired regularizer. That is,
Psi(y,tau) = arg min_x (1/2)*(x - y)^2 + tau phi(x).
Default: in the absence of any Phi given by the user,
it is assumed that phi(x) = ||x||_1 thus
Psi(y,tau) = soft(y,tau)
Important: if Psi is given, phi must also be given,
so that the algorithm may also compute
the objective function.
'StopCriterion' = type of stopping criterion to use
0 = algorithm stops when the relative
change in the number of non-zero
components of the estimate falls
below 'ToleranceA'
1 = stop when the relative
change in the objective function
falls below 'ToleranceA'
2 = stop when relative duality gap
falls below 'ToleranceA'
3 = stop when LCP estimate of relative
distance to solution falls below ToleranceA
4 = stop when the objective function
becomes equal or less than toleranceA.
5 = stop when the norm of the difference between
two consecutive estimates, divided by the norm
of one of them falls below toleranceA
Default = 2
'ToleranceA' = stopping threshold; Default = 0.01
'Debias' = debiasing option: 1 = yes, 0 = no.
Default = 0.
'ToleranceD' = stopping threshold for the debiasing phase:
Default = 0.0001.
If no debiasing takes place, this parameter,
if present, is ignored.
'MaxiterA' = maximum number of iterations allowed in the
main phase of the algorithm.
Default = 1000
'MiniterA' = minimum number of iterations performed in the
main phase of the algorithm.
Default = 5
'MaxiterD' = maximum number of iterations allowed in the
debising phase of the algorithm.
Default = 200
'MiniterD' = minimum number of iterations to perform in the
debiasing phase of the algorithm.
Default = 5
'Initialization' must be one of {0,1,2,array}
0 -> Initialization at zero.
1 -> Random initialization.
2 -> initialization with A'*y.
array -> initialization provided by the user.
Default = 0;
'BB_variant' specifies which variant of Barzila-Borwein to use, or not.
0 -> don't use a BB rule - instead pick the starting alpha
based on the successful value at the previous iteration
1 -> standard BB choice s'r/s's
2 -> inverse BB variant r'r/r's
Default = 1
'BB_cycle' specifies the cycle length - the number of iterations between
recalculation of alpha. Requires integer value at least
1. Relevant only if a **nonmonotone BB rule** is used
(BB_variant = 1 or 2 and Monotone=0).
Default = 1
'Monotone' = enforce monotonic decrease in f, or not?
any nonzero -> enforce monotonicity (overrides 'Safeguard')
0 -> don't enforce monotonicity.
Default = 0;
'Safeguard' = enforce a "sufficient decrease" over the largest
objective value of the past M iterations.
any nonzero -> safeguard
0 -> don't safeguard
Default = 0.
'M' = number of steps to look back in the safeguarding process.
Ignored if Safeguard=0 or if Monotone is nonzero.
(positive integer. Default = 5)
'sigma' = sigma value used in Safeguarding test for sufficient
decrease. Ignored unless 'Safeguard' is nonzero. Must be
in (0,1). Drfault: .01.
'Eta' = factor by which alpha is multiplied within an iteration,
until a decrease in the objective function is
obtained.
Default = 2;
'Alpha_factor' = factor by which to reduce the successful value of
alpha at iteration k, to give the first value of alpha
to be tried at iteration k+1.
If a Barzilai-Borwein rule is specified (BB_variant > 0),
this parameter is ignored.
Default = 0.8;
'Continuation' = Continuation or not (1 or 0)
Specifies the choice for a continuation scheme,
in which we start with a large value of tau, and
then decrease tau until the desired value is
reached. At each value, the solution obtained
with the previous values is used as initialization.
Default = 0
'ContinuationSteps' = Number of steps in the continuation procedure;
ignored if 'Continuation' equals zero.
If -1, an adaptive continuation procedure is used.
Default = -1.
'FirstTauFactor' = Initial tau value, if using continuation, is
obtained by multiplying the given tau by
this factor. This parameter is ignored if
'Continuation' equals zero or
'ContinuationSteps' equals -1.
Default = 10.
'True_x' = if the true underlying x is passed in
this argument, MSE plots are generated.
'AlphaMin' = the alphamin parameter of the BB method.
Default = 1e-30;
'AlphaMax' = the alphamax parameter of the BB method.
Default = 1e30;
'Verbose' = work silently (0) or verbosely (1)
===================================================
============ Outputs ==============================
x = solution of the main algorithm
x_debias = solution after the debiasing phase;
if no debiasing phase took place, this
variable is empty, x_debias = [].
objective = sequence of values of the objective function
times = CPU time after each iteration
debias_start = iteration number at which the debiasing
phase started. If no debiasing took place,
this variable is returned as zero.
mses = sequence of MSE values, with respect to True_x,
if it was given; if it was not given, mses is empty,
mses = [].
========================================================0001 function [x,x_debias,objective,times,debias_start,mses,taus, numA, numAt]= ... 0002 SpaRSA_adpW(y,A,tau,varargin) 0003 0004 % SpaRSA version 2.0, December 31, 2007 0005 % 0006 % This function solves the convex problem 0007 % 0008 % arg min_x = 0.5*|| y - A x ||_2^2 + tau phi(x) 0009 % 0010 % using the SpaRSA algorithm, which is described in "Sparse Reconstruction 0011 % by Separable Approximation" by S. Wright, R. Nowak, M. Figueiredo, 0012 % IEEE Transactions on Signal Processing, 2009 (to appear). 0013 % 0014 % The algorithm is related GPSR (Figueiredo, Nowak, Wright) but does not 0015 % rely on the conversion to QP form of the l1 norm, because it is not 0016 % limited to being used with l1 regularization. Instead it forms a separable 0017 % approximation to the first term of the objective, which has the form 0018 % 0019 % d'*A'*(A x - y) + 0.5*alpha*d'*d 0020 % 0021 % where alpha is obtained from a BB formula. In a monotone variant, alpha is 0022 % increased until we see a decreasein the original objective function over 0023 % this step. 0024 % 0025 % ----------------------------------------------------------------------- 0026 % Copyright (2007): Mario Figueiredo, Robert Nowak, Stephen Wright 0027 % 0028 % GPSR is distributed under the terms 0029 % of the GNU General Public License 2.0. 0030 % 0031 % Permission to use, copy, modify, and distribute this software for 0032 % any purpose without fee is hereby granted, provided that this entire 0033 % notice is included in all copies of any software which is or includes 0034 % a copy or modification of this software and in all copies of the 0035 % supporting documentation for such software. 0036 % This software is being provided "as is", without any express or 0037 % implied warranty. In particular, the authors do not make any 0038 % representation or warranty of any kind concerning the merchantability 0039 % of this software or its fitness for any particular purpose." 0040 % ---------------------------------------------------------------------- 0041 % 0042 % Please check for the latest version of the code and paper at 0043 % www.lx.it.pt/~mtf/SpaRSA 0044 % 0045 % ===== Required inputs ============= 0046 % 0047 % y: 1D vector or 2D array (image) of observations 0048 % 0049 % A: if y and x are both 1D vectors, A can be a 0050 % k*n (where k is the size of y and n the size of x) 0051 % matrix or a handle to a function that computes 0052 % products of the form A*v, for some vector v. 0053 % In any other case (if y and/or x are 2D arrays), 0054 % A has to be passed as a handle to a function which computes 0055 % products of the form A*x; another handle to a function 0056 % AT which computes products of the form A'*x is also required 0057 % in this case. The size of x is determined as the size 0058 % of the result of applying AT. 0059 % 0060 % tau: regularization parameter (scalar) 0061 % 0062 % ===== Optional inputs ============= 0063 % 0064 % 0065 % 'AT' = function handle for the function that implements 0066 % the multiplication by the conjugate of A, when A 0067 % is a function handle. If A is an array, AT is ignored. 0068 % 0069 % 'Psi' = handle to the denoising function, that is, to a function 0070 % that computes the solution of the densoing probelm 0071 % corresponding to the desired regularizer. That is, 0072 % Psi(y,tau) = arg min_x (1/2)*(x - y)^2 + tau phi(x). 0073 % Default: in the absence of any Phi given by the user, 0074 % it is assumed that phi(x) = ||x||_1 thus 0075 % Psi(y,tau) = soft(y,tau) 0076 % Important: if Psi is given, phi must also be given, 0077 % so that the algorithm may also compute 0078 % the objective function. 0079 % 0080 % 'StopCriterion' = type of stopping criterion to use 0081 % 0 = algorithm stops when the relative 0082 % change in the number of non-zero 0083 % components of the estimate falls 0084 % below 'ToleranceA' 0085 % 1 = stop when the relative 0086 % change in the objective function 0087 % falls below 'ToleranceA' 0088 % 2 = stop when relative duality gap 0089 % falls below 'ToleranceA' 0090 % 3 = stop when LCP estimate of relative 0091 % distance to solution falls below ToleranceA 0092 % 4 = stop when the objective function 0093 % becomes equal or less than toleranceA. 0094 % 5 = stop when the norm of the difference between 0095 % two consecutive estimates, divided by the norm 0096 % of one of them falls below toleranceA 0097 % Default = 2 0098 % 0099 % 'ToleranceA' = stopping threshold; Default = 0.01 0100 % 0101 % 'Debias' = debiasing option: 1 = yes, 0 = no. 0102 % Default = 0. 0103 % 0104 % 'ToleranceD' = stopping threshold for the debiasing phase: 0105 % Default = 0.0001. 0106 % If no debiasing takes place, this parameter, 0107 % if present, is ignored. 0108 % 0109 % 'MaxiterA' = maximum number of iterations allowed in the 0110 % main phase of the algorithm. 0111 % Default = 1000 0112 % 0113 % 'MiniterA' = minimum number of iterations performed in the 0114 % main phase of the algorithm. 0115 % Default = 5 0116 % 0117 % 'MaxiterD' = maximum number of iterations allowed in the 0118 % debising phase of the algorithm. 0119 % Default = 200 0120 % 0121 % 'MiniterD' = minimum number of iterations to perform in the 0122 % debiasing phase of the algorithm. 0123 % Default = 5 0124 % 0125 % 'Initialization' must be one of {0,1,2,array} 0126 % 0 -> Initialization at zero. 0127 % 1 -> Random initialization. 0128 % 2 -> initialization with A'*y. 0129 % array -> initialization provided by the user. 0130 % Default = 0; 0131 % 0132 % 'BB_variant' specifies which variant of Barzila-Borwein to use, or not. 0133 % 0 -> don't use a BB rule - instead pick the starting alpha 0134 % based on the successful value at the previous iteration 0135 % 1 -> standard BB choice s'r/s's 0136 % 2 -> inverse BB variant r'r/r's 0137 % Default = 1 0138 % 0139 % 'BB_cycle' specifies the cycle length - the number of iterations between 0140 % recalculation of alpha. Requires integer value at least 0141 % 1. Relevant only if a **nonmonotone BB rule** is used 0142 % (BB_variant = 1 or 2 and Monotone=0). 0143 % Default = 1 0144 % 0145 % 'Monotone' = enforce monotonic decrease in f, or not? 0146 % any nonzero -> enforce monotonicity (overrides 'Safeguard') 0147 % 0 -> don't enforce monotonicity. 0148 % Default = 0; 0149 % 0150 % 'Safeguard' = enforce a "sufficient decrease" over the largest 0151 % objective value of the past M iterations. 0152 % any nonzero -> safeguard 0153 % 0 -> don't safeguard 0154 % Default = 0. 0155 % 0156 % 'M' = number of steps to look back in the safeguarding process. 0157 % Ignored if Safeguard=0 or if Monotone is nonzero. 0158 % (positive integer. Default = 5) 0159 % 0160 % 'sigma' = sigma value used in Safeguarding test for sufficient 0161 % decrease. Ignored unless 'Safeguard' is nonzero. Must be 0162 % in (0,1). Drfault: .01. 0163 % 0164 % 'Eta' = factor by which alpha is multiplied within an iteration, 0165 % until a decrease in the objective function is 0166 % obtained. 0167 % Default = 2; 0168 % 0169 % 'Alpha_factor' = factor by which to reduce the successful value of 0170 % alpha at iteration k, to give the first value of alpha 0171 % to be tried at iteration k+1. 0172 % If a Barzilai-Borwein rule is specified (BB_variant > 0), 0173 % this parameter is ignored. 0174 % Default = 0.8; 0175 % 0176 % 'Continuation' = Continuation or not (1 or 0) 0177 % Specifies the choice for a continuation scheme, 0178 % in which we start with a large value of tau, and 0179 % then decrease tau until the desired value is 0180 % reached. At each value, the solution obtained 0181 % with the previous values is used as initialization. 0182 % Default = 0 0183 % 0184 % 'ContinuationSteps' = Number of steps in the continuation procedure; 0185 % ignored if 'Continuation' equals zero. 0186 % If -1, an adaptive continuation procedure is used. 0187 % Default = -1. 0188 % 0189 % 'FirstTauFactor' = Initial tau value, if using continuation, is 0190 % obtained by multiplying the given tau by 0191 % this factor. This parameter is ignored if 0192 % 'Continuation' equals zero or 0193 % 'ContinuationSteps' equals -1. 0194 % Default = 10. 0195 % 0196 % 'True_x' = if the true underlying x is passed in 0197 % this argument, MSE plots are generated. 0198 % 0199 % 'AlphaMin' = the alphamin parameter of the BB method. 0200 % Default = 1e-30; 0201 % 0202 % 'AlphaMax' = the alphamax parameter of the BB method. 0203 % Default = 1e30; 0204 % 0205 % 'Verbose' = work silently (0) or verbosely (1) 0206 % 0207 % =================================================== 0208 % ============ Outputs ============================== 0209 % x = solution of the main algorithm 0210 % 0211 % x_debias = solution after the debiasing phase; 0212 % if no debiasing phase took place, this 0213 % variable is empty, x_debias = []. 0214 % 0215 % objective = sequence of values of the objective function 0216 % 0217 % times = CPU time after each iteration 0218 % 0219 % debias_start = iteration number at which the debiasing 0220 % phase started. If no debiasing took place, 0221 % this variable is returned as zero. 0222 % 0223 % mses = sequence of MSE values, with respect to True_x, 0224 % if it was given; if it was not given, mses is empty, 0225 % mses = []. 0226 % ======================================================== 0227 0228 %---------------------------------------------------------------- 0229 % Some modification for solving weighted L1 problem 0230 % Added an option to perform adaptive reweighting using flag 0231 % 0232 % adp_wt = adaptive reweighting 0233 % 1 - adaptive reweighting, 0234 % 0 - standard SpaRSA (default) 0235 % 0236 % W_new assigns the values to the weights 0237 % Termination criteria 2 and 3 modified to accomodate weighted L1 norm 0238 % 0239 % Modified: June 2012 by Salman Asif @ Georgia Tech. 0240 %---------------------------------------------------------------- 0241 0242 % start the clock 0243 t0 = cputime; 0244 times(1) = cputime - t0; 0245 0246 % test for number of required parametres 0247 if (nargin-length(varargin)) ~= 3 0248 error('Wrong number of required parameters'); 0249 end 0250 0251 % Set the defaults for the optional parameters 0252 stopCriterion = 2; 0253 tolA = 0.01; 0254 tolD = 0.0001; 0255 debias = 0; 0256 maxiter = 10000; 0257 maxiter_debias = 200; 0258 miniter = 5; 0259 miniter_debias = 0; 0260 init = 0; 0261 bbVariant = 1; 0262 bbCycle = 1; 0263 enforceMonotone = 0; 0264 enforceSafeguard = 0; 0265 M = 5; 0266 sigma = .01; 0267 alphamin = 1e-30; 0268 alphamax = 1e30; 0269 compute_mse = 0; 0270 AT = 0; 0271 verbose = 1; 0272 adp_wt = 0; % adaptive reweighting 0273 W_new = 1; 0274 continuation = 0; 0275 cont_steps = -1; 0276 psi_ok = 0; 0277 phi_ok = 0; 0278 % amount by which to increase alpha after an unsuccessful step 0279 eta = 2.0; 0280 % amount by which to decrease alpha between iterations, if a 0281 % Barzilai-Borwein rule is not used to make the initial guess at each 0282 % iteration. 0283 alphaFactor = 0.8; 0284 phi_l1 = 0; 0285 0286 % Set the defaults for outputs that may not be computed 0287 debias_start = 0; 0288 x_debias = []; 0289 mses = []; 0290 0291 % number of matrix multiplies 0292 numA = 0; numAt = 0; 0293 0294 % Read the optional parameters 0295 if (rem(length(varargin),2)==1) 0296 error('Optional parameters should always go by pairs'); 0297 else 0298 for i=1:2:(length(varargin)-1) 0299 switch upper(varargin{i}) 0300 case 'PSI' 0301 psi_function = varargin{i+1}; 0302 case 'PHI' 0303 phi_function = varargin{i+1}; 0304 case 'STOPCRITERION' 0305 stopCriterion = varargin{i+1}; 0306 case 'TOLERANCEA' 0307 tolA = varargin{i+1}; 0308 case 'TOLERANCED' 0309 tolD = varargin{i+1}; 0310 case 'DEBIAS' 0311 debias = varargin{i+1}; 0312 case 'MAXITERA' 0313 maxiter = varargin{i+1}; 0314 case 'MAXITERD' 0315 maxiter_debias = varargin{i+1}; 0316 case 'MINITERA' 0317 miniter = varargin{i+1}; 0318 case 'MINITERD' 0319 miniter_debias = varargin{i+1}; 0320 case 'INITIALIZATION' 0321 if prod(size(varargin{i+1})) > 1 % we have an initial x 0322 init = 33333; % some flag to be used below 0323 x = varargin{i+1}; 0324 else 0325 init = varargin{i+1}; 0326 end 0327 case 'BB_VARIANT' 0328 bbVariant = varargin{i+1}; 0329 case 'BB_CYCLE' 0330 bbCycle = varargin{i+1}; 0331 case 'MONOTONE' 0332 enforceMonotone = varargin{i+1}; 0333 case 'SAFEGUARD' 0334 enforceSafeguard = varargin{i+1}; 0335 case 'M' 0336 M = varargin{i+1}; 0337 case 'SIGMA' 0338 sigma = varargin{i+1}; 0339 case 'ETA' 0340 eta = varargin{i+1}; 0341 case 'ALPHA_FACTOR' 0342 alphaFactor = varargin{i+1}; 0343 case 'CONTINUATION' 0344 continuation = varargin{i+1}; 0345 case 'CONTINUATIONSTEPS' 0346 cont_steps = varargin{i+1}; 0347 case 'FIRSTTAUFACTOR' 0348 firstTauFactor = varargin{i+1}; 0349 case 'TRUE_X' 0350 compute_mse = 1; 0351 true = varargin{i+1}; 0352 case 'ALPHAMIN' 0353 alphamin = varargin{i+1}; 0354 case 'ALPHAMAX' 0355 alphamax = varargin{i+1}; 0356 case 'AT' 0357 AT = varargin{i+1}; 0358 case 'VERBOSE' 0359 verbose = varargin{i+1}; 0360 case 'ADP_WT' 0361 adp_wt = varargin{i+1}; 0362 case 'W_NEW' 0363 W_new = varargin{i+1}; 0364 otherwise 0365 % Hmmm, something wrong with the parameter string 0366 error(['Unrecognized option: ''' varargin{i} '''']); 0367 end; 0368 end; 0369 end 0370 %%%%%%%%%%%%%% 0371 0372 % it makes no sense to ask for a nonmonotone variant of a non-BB method 0373 if ~enforceMonotone && bbVariant==0 0374 error(['non-monotone, non-BBmethod requested']); 0375 end 0376 0377 if (sum(stopCriterion == [0 1 2 3 4 5])==0) 0378 error(['Unknown stopping criterion']); 0379 end 0380 0381 % if A is a function handle, we have to check presence of AT, 0382 if isa(A, 'function_handle') && ~isa(AT,'function_handle') 0383 error(['The function handle for transpose of A is missing']); 0384 end 0385 0386 % if A is a matrix, we find out dimensions of y and x, 0387 % and create function handles for multiplication by A and A', 0388 % so that the code below doesn't have to distinguish between 0389 % the handle/not-handle cases 0390 if ~isa(A, 'function_handle') 0391 AT = @(x) A'*x; 0392 A = @(x) A*x; 0393 end 0394 % from this point down, A and AT are always function handles. 0395 0396 % Precompute A'*y since it'll be used a lot 0397 Aty = AT(y); numAt = numAt+1; 0398 0399 % if phi was given, check to see if it is a handle and that it 0400 % accepts two arguments 0401 if exist('psi_function','var') 0402 if isa(psi_function,'function_handle') 0403 try % check if phi can be used, using Aty, which we know has 0404 % same size as x 0405 dummy = psi_function(Aty,tau); 0406 psi_ok = 1; 0407 catch 0408 error(['Something is wrong with function handle for psi']) 0409 end 0410 else 0411 error(['Psi does not seem to be a valid function handle']); 0412 end 0413 else %if nothing was given, use soft thresholding 0414 psi_function = @(x,tau) soft(x,tau); 0415 end 0416 0417 % if psi exists, phi must also exist 0418 if (psi_ok == 1) 0419 if exist('phi_function','var') 0420 if isa(phi_function,'function_handle') 0421 try % check if phi can be used, using Aty, which we know has 0422 % same size as x 0423 dummy = phi_function(Aty); 0424 catch 0425 error(['Something is wrong with function handle for phi']) 0426 end 0427 else 0428 error(['Phi does not seem to be a valid function handle']); 0429 end 0430 else 0431 error(['If you give Psi you must also give Phi']); 0432 end 0433 else % if no psi and phi were given, simply use the l1 norm. 0434 phi_function = @(x) sum(abs(x(:))); 0435 phi_l1 = 1; 0436 end 0437 phi_l1 = 1; 0438 0439 % Initialization 0440 switch init 0441 case 0 % initialize at zero, using AT to find the size of x 0442 x = AT(zeros(size(y))); 0443 case 1 % initialize randomly, using AT to find the size of x 0444 x = randn(size(AT(zeros(size(y))))); 0445 case 2 % initialize x0 = A'*y 0446 x = Aty; 0447 case 33333 0448 % initial x was given as a function argument; just check size 0449 if size(A(x)) ~= size(y) 0450 error(['Size of initial x is not compatible with A']); 0451 end 0452 otherwise 0453 error(['Unknown ''Initialization'' option']); 0454 end 0455 0456 % if the true x was given, check its size 0457 if compute_mse & (size(true) ~= size(x)) 0458 error(['Initial x has incompatible size']); 0459 end 0460 0461 % if tau is large enough, in the case of phi = l1, thus psi = soft, 0462 % the optimal solution is the zero vector 0463 if phi_l1 0464 aux = AT(y); numAt = numAt+1; 0465 max_tau = max(abs(aux(:))); 0466 firstTauFactor = 0.8*max_tau / tau; 0467 if (tau >= max_tau) && (psi_ok==0) 0468 x = zeros(size(aux)); 0469 if debias 0470 x_debias = x; 0471 end 0472 objective(1) = 0.5*(y(:)'*y(:)); 0473 times(1) = 0; 0474 if compute_mse 0475 mses(1) = sum(true(:).^2); 0476 end 0477 return 0478 end 0479 end 0480 0481 % define the indicator vector or matrix of nonzeros in x 0482 nz_x = (x ~= 0.0); 0483 num_nz_x = sum(nz_x(:)); 0484 0485 % store given tau, because we're going to change it in the 0486 % continuation procedure 0487 final_tau = tau; 0488 % if we choose to use adaptive continuation, need to reset tau to realmax to 0489 % make things work (don't ask...) 0490 if cont_steps == -1 0491 tau = realmax; 0492 end 0493 0494 % store given stopping criterion and threshold, because we're going 0495 % to change them in the continuation procedure 0496 final_stopCriterion = stopCriterion; 0497 final_tolA = tolA; 0498 0499 % set continuation factors 0500 0501 if (continuation && phi_l1 && (cont_steps > 1)) 0502 % If tau is scalar, first check top see if the first factor is 0503 % too large (i.e., large enough to make the first 0504 % solution all zeros). If so, make it a little smaller than that. 0505 if prod(size(tau)) == 1 0506 if firstTauFactor*tau >= max_tau 0507 firstTauFactor = 0.5 * max_tau / tau; 0508 if verbose 0509 fprintf(1,'\n setting parameter FirstTauFactor\n') 0510 end 0511 end 0512 cont_factors = 10.^[log10(firstTauFactor):... 0513 log10(1/firstTauFactor)/(cont_steps-1):0]; 0514 end 0515 else 0516 if ( ~continuation ) 0517 cont_factors = 1; 0518 cont_steps = 1; 0519 end 0520 end 0521 0522 keep_continuation = 1; 0523 cont_loop = 1; 0524 iter = 1; 0525 taus = []; 0526 0527 % loop for continuation 0528 while keep_continuation 0529 0530 % initialize the count of steps since last update of alpha 0531 % (for use in cyclic BB) 0532 iterThisCycle = 0; 0533 0534 % Compute the initial residual and gradient 0535 resid = A(x) - y; numA = numA+1; 0536 gradq = AT(resid); numAt = numAt+1; 0537 0538 if cont_steps == -1 0539 0540 temp_tau = max(final_tau,0.2*max(abs(gradq(:)))); 0541 0542 if temp_tau > tau 0543 tau = final_tau; 0544 else 0545 tau = temp_tau; 0546 end 0547 0548 if tau == final_tau 0549 stopCriterion = final_stopCriterion; 0550 tolA = final_tolA; 0551 keep_continuation = 0; 0552 else 0553 stopCriterion = 1; 0554 tolA = 1e-5; 0555 end 0556 else 0557 tau = final_tau * cont_factors(cont_loop); 0558 if cont_loop == cont_steps 0559 stopCriterion = final_stopCriterion; 0560 tolA = final_tolA; 0561 keep_continuation = 0; 0562 else 0563 stopCriterion = 1; 0564 tolA = 1e-5; 0565 end 0566 end 0567 0568 taus = [taus tau]; 0569 0570 if verbose 0571 fprintf('\n Regularization parameter tau = %10.6e\n',tau) 0572 end 0573 0574 % compute and store initial value of the objective function 0575 % for this tau 0576 alpha = 1; %1/eps; 0577 0578 f = 0.5*(resid(:)'*resid(:)) + tau * phi_function(x); 0579 if enforceSafeguard 0580 f_lastM = f; 0581 end 0582 0583 % if we are at the very start of the process, store the initial mses and 0584 % objective in the plotting arrays 0585 if cont_loop==1 0586 objective(1) = f; 0587 if compute_mse 0588 mses(1) = (x(:)-true(:))'*(x(:)-true(:)); 0589 end 0590 if verbose 0591 fprintf(1,'Initial obj=%10.6e, alpha=%6.2e, nonzeros=%7d\n',... 0592 f,alpha,num_nz_x); 0593 end 0594 end 0595 0596 %---------------------------------------------------------------------- 0597 % adaptive weight selection 0598 % [modified by Salman Asif @ Georgia Tech June 2012] 0599 %---------------------------------------------------------------------- 0600 if adp_wt 0601 if exist('W_new','var') 0602 gamma = find(abs(x)>tau); 0603 ewt_b = length(y)*(norm(x,2)/norm(x,1))^2; 0604 W_old = W_new; 0605 W_new = ones(length(x),1); 0606 if tau < 5*final_tau; 0607 % [alpha beta epsilon] = weight_param(5,1,x,length(y)); 0608 % W_new = 1/alpha./(beta*(abs(x))+epsilon); 0609 0610 W_new(gamma) = min([W_old(gamma) 1./abs(x(gamma))/ewt_b],[],2); 0611 else 0612 W_new(gamma) = min([W_old(gamma) 1./abs(x(gamma))/ewt_b],[],2); 0613 % W_new(gamma) = W_old(gamma)/2; % 1/2 0614 end 0615 psi_function = @(x,tau) soft(x,tau*W_new); 0616 phi_function = @(x) sum(abs(W_new.*x)); 0617 % figure(115); 0618 % subplot(211); plot([true x]); 0619 % subplot(212); semilogy([W_new]); pause(1); 0620 else 0621 W_new = 1; 0622 psi_function = @(x,tau) soft(x,tau*W_new); 0623 phi_function = @(x) sum(abs(W_new.*x)); 0624 end 0625 end 0626 0627 0628 % initialization of alpha 0629 % alpha = 1/max(max(abs(du(:))),max(abs(dv(:)))); 0630 % or just do a dumb initialization 0631 %alphas(iter) = alpha; 0632 0633 % control variable for the outer loop and iteration counter 0634 keep_going = 1; 0635 0636 while keep_going 0637 0638 % compute gradient 0639 gradq = AT(resid); numAt = numAt+1; 0640 0641 % save current values 0642 prev_x = x; 0643 prev_f = f; 0644 prev_resid = resid; 0645 0646 0647 % computation of step 0648 0649 cont_inner = 1; 0650 while cont_inner 0651 x = psi_function(prev_x - gradq*(1/alpha),tau/alpha); 0652 dx = x - prev_x; 0653 Adx = A(dx); numA = numA+1; 0654 resid = prev_resid + Adx; 0655 f = 0.5*(resid(:)'*resid(:)) + tau * phi_function(x); 0656 if enforceMonotone 0657 f_threshold = prev_f; 0658 elseif enforceSafeguard 0659 f_threshold = max(f_lastM) - 0.5*sigma*alpha*(dx(:)'*dx(:)); 0660 else 0661 f_threshold = inf; 0662 end 0663 % f_threshold 0664 0665 if f <= f_threshold 0666 cont_inner=0; 0667 else 0668 % not good enough, increase alpha and try again 0669 alpha = eta*alpha; 0670 if verbose 0671 fprintf(1,' f=%10.6e, increasing alpha to %6.2e\n', f, alpha); 0672 end 0673 end 0674 end % of while cont_inner 0675 0676 if enforceSafeguard 0677 if length(f_lastM)<M+1 0678 f_lastM = [f_lastM f]; 0679 else 0680 f_lastM = [f_lastM(2:M+1) f]; 0681 end 0682 end 0683 0684 % print stuff 0685 if verbose 0686 fprintf(1,'t=%4d, obj=%10.6e, alpha=%e ', iter, f, alpha ); 0687 end 0688 0689 if bbVariant==1 0690 % standard BB choice of initial alpha for next step 0691 if iterThisCycle==0 | enforceMonotone==1 0692 dd = dx(:)'*dx(:); 0693 dGd = Adx(:)'*Adx(:); 0694 alpha = min(alphamax,max(alphamin,dGd/(realmin+dd))); 0695 end 0696 elseif bbVariant==2 0697 % alternative BB choice of initial alpha for next step 0698 if iterThisCycle==0 | enforceMonotone==1 0699 dd = dx(:)'*dx(:); 0700 dGd = Adx(:)'*Adx(:); 0701 ATAdx=AT(Adx); numAt = numAt+1; 0702 dGGd = ATAdx(:)'*ATAdx(:); 0703 alpha = min(alphamax,max(alphamin,dGGd/(realmin+dGd))); 0704 end 0705 else 0706 % reduce current alpha to get initial alpha for next step 0707 alpha = alpha * alphaFactor; 0708 end 0709 0710 % update iteration counts, store results and times 0711 iter=iter+1; 0712 iterThisCycle=mod(iterThisCycle+1,bbCycle); 0713 objective(iter) = f; 0714 times(iter) = cputime-t0; 0715 % alphas(iter) = alpha; 0716 if compute_mse 0717 err = true - x; 0718 mses(iter) = (err(:)'*err(:)); 0719 end 0720 0721 % compute stopping criteria and test for termination 0722 switch stopCriterion 0723 case 0, 0724 % compute the stopping criterion based on the change 0725 % of the number of non-zero components of the estimate 0726 nz_x_prev = nz_x; 0727 nz_x = (abs(x)~=0.0); 0728 num_nz_x = sum(nz_x(:)); 0729 num_changes_active = (sum(nz_x(:)~=nz_x_prev(:))); 0730 if num_nz_x >= 1 0731 criterionActiveSet = num_changes_active / num_nz_x; 0732 keep_going = (criterionActiveSet > tolA); 0733 end 0734 if verbose 0735 fprintf(1,'Delta nz = %d (target = %e)\n',... 0736 criterionActiveSet , tolA) 0737 end 0738 case 1, 0739 % compute the stopping criterion based on the relative 0740 % variation of the objective function. 0741 criterionObjective = abs(f-prev_f)/(prev_f); 0742 keep_going = (criterionObjective > tolA); 0743 if verbose 0744 fprintf(1,'Delta obj. = %e (target = %e)\n',... 0745 criterionObjective , tolA) 0746 end 0747 case 2, 0748 % compute the "duality" stopping criterion - actually based on the 0749 % iterate PRIOR to the step just taken. Make it relative to the primal 0750 % function value. 0751 scaleFactor = norm(gradq(:)./W_new,inf); 0752 w = tau.*prev_resid(:) / scaleFactor; 0753 criterionDuality = 0.5* (prev_resid(:)'*prev_resid(:)) + ... 0754 tau * phi_function(prev_x) + 0.5*w(:)'*w(:) + y(:)'*w(:); 0755 criterionDuality = criterionDuality / prev_f; 0756 keep_going = (criterionDuality > tolA); 0757 if verbose 0758 fprintf(1,'Duality = %e (target = %e)\n',... 0759 criterionDuality , tolA) 0760 end 0761 case 3, 0762 % compute the "LCP" stopping criterion - again based on the previous 0763 % iterate. Make it "relative" to the norm of x. 0764 w = [ min(tau.*W_new + gradq(:), max(prev_x(:),0.0)); ... 0765 min(tau.*W_new - gradq(:), max(-prev_x(:),0.0))]; 0766 criterionLCP = norm(w(:), inf); 0767 criterionLCP = criterionLCP / max(1.0e-6, norm(prev_x(:),inf)); 0768 keep_going = (criterionLCP > tolA); 0769 if verbose 0770 fprintf(1,'LCP = %e (target = %e)\n',criterionLCP,tolA) 0771 end 0772 case 4, 0773 % continue if not yeat reached target value tolA 0774 keep_going = (f > tolA); 0775 if verbose 0776 fprintf(1,'Objective = %e (target = %e)\n',f,tolA) 0777 end 0778 case 5, 0779 % stopping criterion based on relative norm of step taken 0780 delta_x_criterion = sqrt(dx(:)'*dx(:))/(x(:)'*x(:)); 0781 keep_going = (delta_x_criterion > tolA); 0782 if verbose 0783 fprintf(1,'Norm(delta x)/norm(x) = %e (target = %e)\n',... 0784 delta_x_criterion,tolA) 0785 end 0786 otherwise, 0787 error(['Unknown stopping criterion']); 0788 end % end of the stopping criteria switch 0789 0790 % overrule the stopping decision to ensure we take between miniter and 0791 % maxiter iterations 0792 if iter<=miniter 0793 % take no fewer than miniter... 0794 keep_going = 1; 0795 elseif iter > maxiter 0796 % and no more than maxiter iterations 0797 keep_going = 0; 0798 end 0799 0800 end % end of the main loop of the GPBB algorithm (while keep_going) 0801 0802 cont_loop = cont_loop + 1; 0803 0804 end % end of the continuation loop (while keep_continuation) 0805 0806 % Print results 0807 if verbose 0808 fprintf(1,'\nFinished the main algorithm! Results:\n') 0809 fprintf(1,'Number of iterations = %d\n',iter) 0810 fprintf(1,'0.5*||A x - y ||_2^2 = %10.3e\n',0.5*resid(:)'*resid(:)) 0811 fprintf(1,'tau * Penalty = %10.3e\n',tau * phi_function(x)) 0812 fprintf(1,'Objective function = %10.3e\n',f); 0813 fprintf(1,'Number of non-zero components = %d\n',sum(x(:)~=0)); 0814 fprintf(1,'CPU time so far = %10.3e\n', times(iter)); 0815 fprintf(1,'\n'); 0816 end 0817 0818 % If the 'Debias' option is set to 1, we try to 0819 % remove the bias from the l1 penalty, by applying CG to the 0820 % least-squares problem obtained by omitting the l1 term 0821 % and fixing the zero coefficients at zero. 0822 0823 if (debias & (sum(x(:)~=0)~=0)) 0824 if verbose 0825 fprintf(1,'\nStarting the debiasing phase...\n\n') 0826 end 0827 0828 x_debias = x; 0829 zeroind = (x_debias~=0); 0830 cont_debias_cg = 1; 0831 debias_start = iter; 0832 0833 % calculate initial residual 0834 resid = A(x_debias); numA = numA+1; 0835 resid = resid-y; 0836 prev_resid = eps*ones(size(resid)); 0837 0838 rvec = AT(resid); numAt = numAt+1; 0839 0840 % mask out the zeros 0841 rvec = rvec .* zeroind; 0842 rTr_cg = rvec(:)'*rvec(:); 0843 0844 % set convergence threshold for the residual || RW x_debias - y ||_2 0845 tol_debias = tolD * (rvec(:)'*rvec(:)); 0846 0847 % initialize pvec 0848 pvec = -rvec; 0849 0850 % main loop 0851 while cont_debias_cg 0852 0853 % calculate A*p = Wt * Rt * R * W * pvec 0854 RWpvec = A(pvec); numA = numA+1; 0855 Apvec = AT(RWpvec); numAt = numAt+1; 0856 0857 % mask out the zero terms 0858 Apvec = Apvec .* zeroind; 0859 0860 % calculate alpha for CG 0861 alpha_cg = rTr_cg / (pvec(:)'* Apvec(:)); 0862 0863 % take the step 0864 x_debias = x_debias + alpha_cg * pvec; 0865 resid = resid + alpha_cg * RWpvec; 0866 rvec = rvec + alpha_cg * Apvec; 0867 0868 rTr_cg_plus = rvec(:)'*rvec(:); 0869 beta_cg = rTr_cg_plus / rTr_cg; 0870 pvec = -rvec + beta_cg * pvec; 0871 0872 rTr_cg = rTr_cg_plus; 0873 0874 iter = iter+1; 0875 0876 objective(iter) = 0.5*(resid(:)'*resid(:)) + ... 0877 tau * phi_function(x_debias(:)); 0878 times(iter) = cputime - t0; 0879 0880 if compute_mse 0881 err = true - x_debias; 0882 mses(iter) = (err(:)'*err(:)); 0883 end 0884 0885 % in the debiasing CG phase, always use convergence criterion 0886 % based on the residual (this is standard for CG) 0887 if verbose 0888 fprintf(1,'t = %5d, debias resid = %13.8e, convergence = %8.3e\n', ... 0889 iter, resid(:)'*resid(:), rTr_cg / tol_debias); 0890 end 0891 cont_debias_cg = ... 0892 (iter-debias_start <= miniter_debias )| ... 0893 ((rTr_cg > tol_debias) & ... 0894 (iter-debias_start <= maxiter_debias)); 0895 0896 end 0897 if verbose 0898 fprintf(1,'\nFinished the debiasing phase! Results:\n') 0899 fprintf(1,'Final number of iterations = %d\n',iter); 0900 fprintf(1,'0.5*||A x - y ||_2 = %10.3e\n',0.5*resid(:)'*resid(:)) 0901 fprintf(1,'tau * penalty = %10.3e\n',tau * phi_function(x)) 0902 fprintf(1,'Objective function = %10.3e\n',f); 0903 fprintf(1,'Number of non-zero components = %d\n',... 0904 sum((x_debias(:)~=0.0))); 0905 fprintf(1,'CPU time so far = %10.3e\n', times(iter)); 0906 fprintf(1,'\n'); 0907 end 0908 end 0909 0910 if compute_mse 0911 mses = mses/length(true(:)); 0912 end 0913