nikcleju@10: function s=SL0_approx(A, x, eps, sigma_min, sigma_decrease_factor, mu_0, L, A_pinv, true_s) nikcleju@10: % nikcleju@10: % SL0(A, x, sigma_min, sigma_decrease_factor, mu_0, L, A_pinv, true_s) nikcleju@10: % nikcleju@10: % Returns the sparsest vector s which satisfies underdetermined system of nikcleju@10: % linear equations A*s=x, using Smoothed L0 (SL0) algorithm. Note that nikcleju@10: % the matrix A should be a 'wide' matrix (more columns than rows). The nikcleju@10: % number of the rows of matrix A should be equal to the length of the nikcleju@10: % column vector x. nikcleju@10: % nikcleju@10: % The first 3 arguments should necessarily be provided by the user. The nikcleju@10: % other parameters have defult values calculated within the function, or nikcleju@10: % may be provided by the user. nikcleju@10: % nikcleju@10: % Sequence of Sigma (sigma_min and sigma_decrease_factor): nikcleju@10: % This is a decreasing geometric sequence of positive numbers: nikcleju@10: % - The first element of the sequence of sigma is calculated nikcleju@10: % automatically. The last element is given by 'sigma_min', and the nikcleju@10: % change factor for decreasing sigma is given by 'sigma_decrease_factor'. nikcleju@10: % - The default value of 'sigma_decrease_factor' is 0.5. Larger value nikcleju@10: % gives better results for less sparse sources, but it uses more steps nikcleju@10: % on sigma to reach sigma_min, and hence it requires higher nikcleju@10: % computational cost. nikcleju@10: % - There is no default value for 'sigma_min', and it should be nikcleju@10: % provided by the user (depending on his/her estimated source noise nikcleju@10: % level, or his/her desired accuracy). By `noise' we mean here the nikcleju@10: % noise in the sources, that is, the energy of the inactive elements of nikcleju@10: % 's'. For example, by the noiseless case, we mean the inactive nikcleju@10: % elements of 's' are exactly equal to zero. As a rule of tumb, for the nikcleju@10: % noisy case, sigma_min should be about 2 to 4 times of the standard nikcleju@10: % deviation of this noise. For the noiseless case, smaller 'sigma_min' nikcleju@10: % results in better estimation of the sparsest solution, and hence its nikcleju@10: % value is determined by the desired accuracy. nikcleju@10: % nikcleju@10: % mu_0: nikcleju@10: % The value of mu_0 scales the sequence of mu. For each vlue of nikcleju@10: % sigma, the value of mu is chosen via mu=mu_0*sigma^2. Note that this nikcleju@10: % value effects Convergence. nikcleju@10: % The default value is mu_0=2 (see the paper). nikcleju@10: % nikcleju@10: % L: nikcleju@10: % number of iterations of the internal (steepest ascent) loop. The nikcleju@10: % default value is L=3. nikcleju@10: % nikcleju@10: % A_pinv: nikcleju@10: % is the pseudo-inverse of matrix A defined by A_pinv=A'*inv(A*A'). nikcleju@10: % If it is not provided, it will be calculated within the function. If nikcleju@10: % you use this function for solving x(t)=A s(t) for different values of nikcleju@10: % 't', it would be a good idea to calculate A_pinv outside the function nikcleju@10: % to prevent its re-calculation for each 't'. nikcleju@10: % nikcleju@10: % true_s: nikcleju@10: % is the true value of the sparse solution. This argument is for nikcleju@10: % simulation purposes. If it is provided by the user, then the function nikcleju@10: % will calculate the SNR of the estimation for each value of sigma and nikcleju@10: % it provides a progress report. nikcleju@10: % nikcleju@10: % Authors: Massoud Babaie-Zadeh and Hossein Mohimani nikcleju@10: % Version: 1.4 nikcleju@10: % Last modified: 4 April 2010. nikcleju@10: % nikcleju@10: % nikcleju@10: % Web-page: nikcleju@10: % ------------------ nikcleju@10: % http://ee.sharif.ir/~SLzero nikcleju@10: % nikcleju@10: % Code History: nikcleju@10: %-------------- nikcleju@10: % Version 2.0: 4 April 2010 nikcleju@10: % Doing a few small modifications that enable the code to work also nikcleju@10: % for complex numbers (not only for real numbers). nikcleju@10: % nikcleju@10: % Version 1.3: 13 Sep 2008 nikcleju@10: % Just a few modifications in the comments nikcleju@10: % nikcleju@10: % Version 1.2: Adding some more comments in the help section nikcleju@10: % nikcleju@10: % Version 1.1: 4 August 2008 nikcleju@10: % - Using MATLAB's pseudo inverse function to generalize for the case nikcleju@10: % the matrix A is not full-rank. nikcleju@10: % nikcleju@10: % Version 1.0 (first official version): 4 July 2008. nikcleju@10: % nikcleju@10: % First non-official version and algorithm development: Summer 2006 nikcleju@10: nikcleju@10: if nargin < 5 nikcleju@10: sigma_decrease_factor = 0.5; nikcleju@10: A_pinv = pinv(A); nikcleju@10: mu_0 = 2; nikcleju@10: L = 3; nikcleju@10: ShowProgress = logical(0); nikcleju@10: elseif nargin == 5 nikcleju@10: A_pinv = pinv(A); nikcleju@10: mu_0 = 2; nikcleju@10: L = 3; nikcleju@10: ShowProgress = logical(0); nikcleju@10: elseif nargin == 6 nikcleju@10: A_pinv = pinv(A); nikcleju@10: L = 3; nikcleju@10: ShowProgress = logical(0); nikcleju@10: elseif nargin == 7 nikcleju@10: A_pinv = pinv(A); nikcleju@10: ShowProgress = logical(0); nikcleju@10: elseif nargin == 8 nikcleju@10: ShowProgress = logical(0); nikcleju@10: elseif nargin == 9 nikcleju@10: ShowProgress = logical(1); nikcleju@10: else nikcleju@10: error('Error in calling SL0_approx function'); nikcleju@10: end nikcleju@10: nikcleju@10: nikcleju@10: % Initialization nikcleju@10: %s = A\x; nikcleju@10: s = A_pinv*x; nikcleju@15: sigma = 2*max(abs(s)); nikcleju@10: nikcleju@10: % Main Loop nikcleju@10: while sigma>sigma_min nikcleju@10: for i=1:L nikcleju@10: delta = OurDelta(s,sigma); nikcleju@10: % Update s in the direction of steepest ascent nikcleju@10: s = s - mu_0*delta; nikcleju@10: % At this point, s no longer exactly satisfies x = A*s nikcleju@10: % The original SL0 algorithm projects s onto {s | x = As} with nikcleju@10: %s = s - A_pinv*(A*s-x); % Projection nikcleju@10: % We want to project s onto {s | |x-As| < eps} nikcleju@10: % We move onto the direction -A_pinv*(A*s-x), but only with a nikcleju@10: % smaller step: nikcleju@10: dir = A_pinv*(A*s-x); nikcleju@10: if norm(A*dir) >= eps nikcleju@10: s = s - (1-eps/norm(A*dir)) * dir; nikcleju@10: end nikcleju@12: assert(norm(x - A*s) < eps + 1e-6) nikcleju@10: end nikcleju@10: nikcleju@10: if ShowProgress nikcleju@10: fprintf(' sigma=%f, SNR=%f\n',sigma,estimate_SNR(s,true_s)) nikcleju@10: end nikcleju@10: nikcleju@10: sigma = sigma * sigma_decrease_factor; nikcleju@10: end nikcleju@10: nikcleju@10: nikcleju@10: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nikcleju@10: function delta=OurDelta(s,sigma) nikcleju@10: nikcleju@10: delta = s.*exp(-abs(s).^2/sigma^2); nikcleju@10: nikcleju@10: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nikcleju@10: function SNR=estimate_SNR(estim_s,true_s) nikcleju@10: nikcleju@10: err = true_s - estim_s; nikcleju@10: SNR = 10*log10(sum(abs(true_s).^2)/sum(abs(err).^2));