压缩感知和稀疏恢复中的L1同伦算法

% Copyright �2010. The Regents of the University of California (Regents). % All Rights Reserved. Contact The Office of Technology Licensing, % UC Berkeley, 2150 Shattuck Avenue, Suite 510, Berkeley, CA 94720-1620, % (510) 643-7201, for commercial licensing opportunities. % Authors: Arvind Ganesh, Allen Y. Yang, Zihan Zhou. % Contact: Allen Y. Yang, Department of EECS, University of California, % Berkeley. <yang@eecs.berkeley.edu> % IN NO EVENT SHALL REGENTS BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, % SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, % ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF % REGENTS HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. % REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED % TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A % PARTICULAR PURPOSE. THE SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF ANY, % PROVIDED HEREUNDER IS PROVIDED "AS IS". REGENTS HAS NO OBLIGATION TO % PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. %% This function is modified from Matlab Package SpaRSA function [x,iter]= SolveSpaRSA(A,y,varargin) % SpaRSA version 2.0, December 31, 2007 % % This function solves the convex problem % % arg min_x = 0.5*|| y - A x ||_2^2 + lambda 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. % % lambda: 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,lambda) = arg min_x (1/2)*(x - y)^2 + lambda phi(x). % Default: in the absence of any Phi given by the user, % it is assumed that phi(x) = ||x||_1 thus % Psi(y,lambda) = soft(y,lambda) % 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 'tolerance' % 1 = stop when the relative % change in the objective function % falls below 'tolerance' % 2 = stop when relative duality gap % falls below 'tolerance' % 3 = stop when relative noise magnitude % falls below 'tolerance' % 4 = stop when the objective function % becomes equal or less than toleranceA. % Default = 3 % % 'Tolerance' = 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