采用凸松弛算法 BP、LASSO、和贪婪算法 OMP、stOMP总共四种重建算法实现压缩感知重建

function [x,y,z,inform,PDitns,CGitns,time] = ... pdco( Fname,Aname,b,bl,bu,d1,d2,options,x0,y0,z0,xsize,zsize ) %----------------------------------------------------------------------- % pdco.m: Primal-Dual Barrier Method for Convex Objectives (23 Sep 2003) %----------------------------------------------------------------------- % [x,y,z,inform,PDitns,CGitns,time] = ... % pdco(Fname,Aname,b,bl,bu,d1,d2,options,x0,y0,z0,xsize,zsize); % % solves optimization problems of the form % % minimize phi(x) + 1/2 norm(D1*x)^2 + 1/2 norm(r)^2 % x,r % subject to A*x + D2*r = b, bl <= x <= bu, r unconstrained, % % where % phi(x) is a smooth convex function defined by function Fname; % A is an m x n matrix defined by matrix or function Aname; % b is a given m-vector; % D1, D2 are positive-definite diagonal matrices defined from d1, d2. % In particular, d2 indicates the accuracy required for % satisfying each row of Ax = b. % % D1 and D2 (via d1 and d2) provide primal and dual regularization % respectively. They ensure that the primal and dual solutions % (x,r) and (y,z) are unique and bounded. % % A scalar d1 is equivalent to d1 = ones(n,1), D1 = diag(d1). % A scalar d2 is equivalent to d2 = ones(m,1), D2 = diag(d2). % Typically, d1 = d2 = 1e-4. % These values perturb phi(x) only slightly (by about 1e-8) and request % that A*x = b be satisfied quite accurately (to about 1e-8). % Set d1 = 1e-4, d2 = 1 for least-squares problems with bound constraints. % The problem is then equivalent to % % minimize phi(x) + 1/2 norm(d1*x)^2 + 1/2 norm(A*x - b)^2 % subject to bl <= x <= bu. % % More generally, d1 and d2 may be n and m vectors containing any positive % values (preferably not too small, and typically no larger than 1). % Bigger elements of d1 and d2 improve the stability of the solver. % % At an optimal solution, if x(j) is on its lower or upper bound, % the corresponding z(j) is positive or negative respectively. % If x(j) is between its bounds, z(j) = 0. % If bl(j) = bu(j), x(j) is fixed at that value and z(j) may have % either sign. % % Also, r and y satisfy r = D2 y, so that Ax + D2^2 y = b. % Thus if d2(i) = 1e-4, the i-th row of Ax = b will be satisfied to % approximately 1e-8. This determines how large d2(i) can safely be. % % % EXTERNAL FUNCTIONS: % options = pdcoSet; provided with pdco.m % [obj,grad,hess] = Fname( x ); provided by user % y = Aname( name,mode,m,n,x ); provided by user if pdMat % is a string, not a matrix % % INPUT ARGUMENTS: % Fname may be an explicit n x 1 column vector c, % or a string containing the name of a function Fname.m %%%!!! Revised 12/16/04 !!! % (Fname cannot be a function handle) % such that [obj,grad,hess] = Fname(x) defines % obj = phi(x) : a scalar, % grad = gradient of phi(x) : an n-vector, % hess = diag(Hessian of phi): an n-vector. % Examples: % If phi(x) is the linear function c'*x, Fname could be % be the vector c, or the name or handle of a function % that returns % [obj,grad,hess] = [c'*x, c, zeros(n,1)]. % If phi(x) is the entropy function E(x) = sum x(j) log x(j), % Fname should return % [obj,grad,hess] = [E(x), log(x)+1, 1./x]. % Aname may be an explicit m x n matrix A (preferably sparse!), % or a string containing the name of a function Aname.m %%%!!! Revised 12/16/04 !!! % (Aname cannot be a function handle) % such that y = aname( name,mode,m,n,x ) % returns y = A*x (mode=1) or y = A'*x (mode=2). % The input parameter "name" will be the string 'Aname' % or whatever the name of the actual function is. % b is an m-vector. % bl is an n-vector of lower bounds. Non-existent bounds % may be represented by bl(j) = -Inf or bl(j) <= -1e+20. % bu is an n-vector of upper bounds. Non-existent bounds % may be represented by bu(j) = Inf or bu(j) >= 1e+20. % d1, d2 may be positive scalars or positive vectors (see above). % options is a structure that may be set and altered by pdcoSet % (type help pdcoSet). % x0, y0, z0 provide an initial solution. % xsize, zsize are estimates of the biggest x and z at the solution. % They are used to scale (x,y,z). Good estimates % should improve the performance of the barrier method. % % % OUTPUT ARGUMENTS: % x is the primal solution. % y is the dual solution associated with Ax + D2 r = b. % z is the dual solution associated with bl <= x <= bu. % inform = 0 if a solution is found; % = 1 if too many iterations were required; % = 2 if the linesearch failed too often. % = 3 if the step lengths became too small. % PDitns is the number of Primal-Dual Barrier iterations required. % CGitns is the number of Conjugate-Gradient iterations required % if an iterative solver is used (LSQR). % time is the cpu time used. %---------------------------------------------------------------------- % PRIVATE FUNCTIONS: % pdxxxbounds % pdxxxdistrib % pdxxxlsqr % pdxxxlsqrmat % pdxxxmat % pdxxxmerit % pdxxxresid1 % pdxxxresid2 % pdxxxstep % % GLOBAL VARIABLES: % global pdDDD1 pdDDD2 pdDDD3 % % % NOTES: % The matrix A should be reasonably well scaled: norm(A,inf) =~ 1. % The vector b and objective phi(x) may be of any size, but ensure that % xsize and zsize are reasonably close to norm(x,inf) and norm(z,inf) % at the solution. % % The files defining Fname and Aname % must not be called Fname.m or Aname.m!! % % % AUTHOR: % Michael Saunders, Systems Optimization Laboratory (SOL), % Stanford University, Stanford, California, USA. % saunders@stanford.edu % % CONTRIBUTORS: % Byunggyoo Kim, Samsung, Seoul, Korea. % hightree@samsung.co.kr % % DEVELOPMENT: % 20 Jun 1997: Original version of pdsco.m derived from pdlp0.m. % 29 Sep 2002: Original version of pdco.m derived from pdsco.m. % Introduced D1, D2 in place of gamma*I, delta*I % and allowed for general bounds bl <= x <= bu. % 06 Oct 2002: Allowed for fixed variabes: bl(j) = bu(j) for any j. % 15 Oct 2002: Eliminated some work vectors (since m, n might be LARGE). % Modularized residuals, linesearch % 16 Oct 2002: pdxxx..., pdDDD... names rationalized. % pdAAA eliminated (global copy of A). % Aname is now used directly as an explicit A or a function. % NOTE: If Aname is a function, it now has an extra parameter. % 23 Oct 2002: Fname and Aname can now be function handles. % 01 Nov 2002: Bug fixed in feval in pdxxxmat. % 19 Apr 2003: Bug fixed in pdxxxbounds. % 07 Aug 2003: Let d1, d2 be scalars if input that way. % 10 Aug 2003: z isn't needed except at the end for output. % 10 Aug 2003: mu0 is now an absolute value -- the initial mu. % 13 Aug 2003: Access only z1(low) and z2(upp) everywhere. % stepxL, stepxU introduced to keep x within bounds. % (With poor starting points, dx may take x outside, % where phi(x) may not be defined. % Entropy once gave complex values for the gradient!) % 16 Sep 2003: Fname can now be a vector c, implying a linear obj c'*x. % 19 Sep 2003: Large system K4 dv = rhs implemented. % 23 Sep 2003: Options LSproblem and LSmethod replaced by Method. % 18 Nov 2003: stepxL, stepxU gave trouble on lptest (see 13 Aug 2003). % Disabled them for now. Nonlinear problems need good x0. % 19 Nov 2003: Bugs with x(fix) and z(fix). % In particular, x(fix) = bl(fix) throughout, so Objective % in iteration log is correct for LPs with explicit c vector. %----------------